【LAMMPS如何系列】计算体积模量
(在阅读本文前,建议你先阅读:【LAMMPS如何系列】计算平衡晶格常数)
体积模量(Bulk Modulus)是材料很常用的一个属性,下面是维基百科里的解释。
The bulk modulus (
) 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)
即使一种材料体积减小一定程度时需要的增加的压强。
推导体积模量另一种形式的表达式
对于立方原胞,压强可以有如下定义:
上式中,E 是晶胞总能量,M是晶胞中的原子数,a是晶格常数。
将上式代入体积模量的定义式,可以得到体积模量的另外一种表达:
其中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 % Poweed by Xianbao Duan @ Lab. of Advanced Material Design % Email: Xianbao.d at gmail.com % ********************************************************************** 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。
已有 7 条评论
博主好 謝謝您的分享
最近在做硬度分析
如果只需要計算壓強
最後幾行的程式碼該如何更改呢?
admin您好,我想请教一个问题,我拟合了一个二元作用势,如何用lammps计算c11、c12、c44弹性常数 来验证 这个作用势的准确性呢?
通过计算“ 平衡态晶格常数、结合能及本文所说的体积模量”,来与实验值进行对比,可否作为验证作用势的方法呢?
谢谢回复!
还有这个单位的转换我觉得有点问题
你是指哪个地方啊?
cvt_factor 是将eV/A^3 转换成GPa,就是这个值,你可以再算一下。
我算了一下,感觉对拟合阶数的依赖比较大。不知道这个问题能否更好的解决。
这个是正常的,多项式拟合就是存在这个特点。
如果能够将拟合的范围缩小到0.1Angstrom以内,然后采用3~6次拟合,效果应该会好一些。当然还有一些更加精确的方法,因为比较麻烦,这里就不介绍了。