返回首页

LAMMPS如何计算体积模量

时间:2012-11-08 13:15来源:知行网www.zhixing123.cn 编辑:麦田守望者

体积模量(Bulk Modulus)是材料很常用的一个属性,下面是维基百科里的解释。

The bulk modulus (K) of a substance measures the substance’s resistance to uniform compression. It is defined as the pressure increase needed to decrease the volume by a factor of 1/e.    (From wikipedia: Bulk Modulus)

definition of bulk modulus

即使一种材料体积减小一定程度时需要的增加的压强。

推导体积模量另一种形式的表达式

对于立方原胞,压强可以有如下定义:

definition of pressure for cubic cell

上式中,E 是晶胞总能量,M是晶胞中的原子数,a是晶格常数。

将上式代入体积模量的定义式,可以得到体积模量的另外一种表达:

another definition of bulk modulus

其中a0是平衡晶格常数。

分析以上表达式,晶胞中的原子数M很容易得到,平衡晶格常数在前面的文章已经介绍如何计算,而d2E/da2|a0只是在计算平衡晶格常数时进一步求二次导数就可以获得。因此,在计算了平衡晶格常数后,并不需要进行新的计算就可以获得体积模量。

计算体积模量

下面接着计算平衡晶格常数的计算,进一步计算体积模量,仍然以金刚石为例。

金刚石为diamond类型,每个原胞中有8个原子,即M=8。仍然使用计算平衡晶格常数中的晶格常数-结合能数据。

下面的MATLAB计算程序可以直接给出体积模量。

cal_bulk_modulus.m

% calcuate the bulk modulus according to "lat_const cohesive energy"
%
% Paramters:
% N: input. int. order of the polynomial fitting.
% M: number of atoms in the unit cell
%       sc: M=1;    bcc: M=2;   fcc: M=4;   dc: M=8
% inFileNmae: input. string. name of the file storing "lat_const cohesive energy"
%
% Example:
% cal_bulk_modulus(6,8,'Au')
%
% Create: 2012-1-9     Complete: 2012-1-9

% *****************************

 function cal_bulk_modulus(N,M,inFileName)

cvt_factor = 160.22;        % convert the modulus from eV/A^2 to GPa

data = load(inFileName,'-ascii');       % read in data from the file
x = data(:,1);  y = data(:,2);
% [x y] =  textread(inFileName,'%f %f');
bindEnergy = polyfit(x,y,N);    % polynomial fitting
dbindEnergy = polyder(bindEnergy);  % derivation of the polynomial equation
zero_points = roots(dbindEnergy);    % solve the zero points
for i = 1: length(zero_points)
    if isreal(zero_points(i))
        if zero_points(i) > x(1)
            if zero_points(i) < x(end)
                lat_const = zero_points(i)
                coh_energy = spline(x,y,lat_const)
            end
        end
    end
end
d2 = polyder(dbindEnergy);
d2_da = polyval(d2,lat_const);
bulk_modulus = M*d2_da*cvt_factor/(9*lat_const)

已经注释的比较清楚了,不做更多解释了。

如果使用4阶拟合(即N=4),得到金刚石的体积模量为:425.7265GPa,实验值为442 GPa。

------分隔线----------------------------
标签(Tag):MATLAB Matlab技巧 Matlab实例教程 matlab源代码 matlab基础教程
------分隔线----------------------------
推荐内容
猜你感兴趣