【LAMMPS如何系列】计算体积模量

(在阅读本文前,建议你先阅读:【LAMMPS如何系列】计算平衡晶格常数

体积模量(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
% 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。

标签: matlab, lammps, input file

相关文章推荐

添加新评论 (无需注册,可直接评论)

已有 7 条评论

  1. eason1021

    博主好 謝謝您的分享
    最近在做硬度分析
    如果只需要計算壓強
    最後幾行的程式碼該如何更改呢?

  2. lammps_vasp

    admin您好,我想请教一个问题,我拟合了一个二元作用势,如何用lammps计算c11、c12、c44弹性常数 来验证 这个作用势的准确性呢?

    通过计算“ 平衡态晶格常数、结合能及本文所说的体积模量”,来与实验值进行对比,可否作为验证作用势的方法呢?
    谢谢回复!

    1. 我爱搜集网博主
      如何计算c11,c12,c44可以参考lammps提供的例子。 计算这些物理性质可以作为验证势准确与否的方法。但对于特定的应用,可能还会需要更多的性质。
  3. 周艳光

    还有这个单位的转换我觉得有点问题

    1. 我爱搜集网博主

      你是指哪个地方啊?
      cvt_factor 是将eV/A^3 转换成GPa,就是这个值,你可以再算一下。

  4. 周艳光

    我算了一下,感觉对拟合阶数的依赖比较大。不知道这个问题能否更好的解决。

    1. 我爱搜集网博主

      这个是正常的,多项式拟合就是存在这个特点。
      如果能够将拟合的范围缩小到0.1Angstrom以内,然后采用3~6次拟合,效果应该会好一些。当然还有一些更加精确的方法,因为比较麻烦,这里就不介绍了。