【LAMMPS如何系列】计算热膨胀系数

热膨胀系数(thermal expansion coefficent)是物质因温度改变时,体积发生变化的趋势。

维基百科的定义:Thermal expansion is the tendency of matter to change in volume in response to a change in temperature.  (Link:thermal expansion)

热膨胀系数有体膨胀系数 β=ΔV/(V*ΔT) 和 线膨胀系数 α=ΔL/(L*ΔT) 。

从这个定义大致就可以获得计算热膨胀系数的思路:计算不同温度下物质的体积或长度,那么体积或长度随着温度的变化曲线的斜率就是ΔL/ΔT了。所以,在模拟的过程中,体积是变化的,那么自然就会想到使用NPT系综。在使用NPT系综之前和之后,建议使用NVT系综进行平衡一定的步数。

举例计算热膨胀系数

下面仍然以铜为例,具体介绍如何使用LAMMPS计算线膨胀系数α。

输入脚本:in.thermal.expansion

# This input script is used to calculate
# the thermal expansion of copper.
# Powered by Xianbao Duan
# Email: xianbao.d@gmail.com
# Website: http://www.52souji.net/

units           metal
boundary        p p p
atom_style      atomic

variable	i loop 10
variable        x equal 200+100*$i

lattice		fcc 3.62
region		box block 0 8 0 8 0 8
create_box	1 box
create_atoms 	1 box

pair_style      eam
pair_coeff      1 1 Cu_u3.eam

variable        N equal step
variable        pote equal pe
variable        Etotal equal etotal
variable        T equal temp
variable        Press equal press
variable        V equal vol

velocity        all create 2.5 825577 dist gaussian

timestep	0.002
thermo		1000

fix		1 all nvt temp 2.5 2.5 0.2
run		2000

unfix		1

fix             2 all npt temp 2.5 $x 0.2 iso 0 0 10
run             120000

unfix		2

fix             extra all print 100 "${N} ${T} ${V} ${pote} ${Etotal} ${Press}" append data

fix		3 all nvt temp $x $x 0.2
run		2000

unfix		3
unfix		extra

clear
next		i
jump		in.thermal.expansion

对以上脚本进行简单解释:

1. 第12行定义的变量x为温度,即计算这些温度下体系的晶格常数。

2. 第34-35行将体系在2.5K下进行NVT平衡,为NPT升温过程做准备。

3. 第39-40行将体系在NPT系综下进行升温。

4. 第44行定义输出有关变量。注意这里使用的是append关键字,而不是file,是为了让脚本在循环执行时不重写这个文件。

5. 第46-47行将升温后的体系在NVT下进行适当平衡,在这个过程中输出了有关变量到文件中。

得到的data文件即为结果数据。将其中的温度和体积复制出来列表,绘图,如下图。

thermal expansion coefficient - temperature vs lattice constant

图中曲线的斜率是6.89444e-5,即ΔL/ΔT= 6.89444e-5 A/K。又初始的晶格常数大致为3.6346A,所以可以算出线膨胀系数α=ΔL/(L*ΔT) =18.97 10-6 K-1

实验结果为17.5 10-6 K-1,两者基本吻合。

材料的热膨胀系数可以参考:热膨胀系数查询表

标签: lammps, input file

相关文章推荐

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

已有 31 条评论

  1. 王勇俊

    博主你好,我问一下运行多少步会达到多少温度是怎么确定的?

  2. hkt

    我想问一下,为什么不需要提供原子质量信息,谢谢。

    1. eam势里已经包括,不需要提供质量。

  3. 11

    data文件里面没有数据

  4. 可否加一下联系方式,诚心向大神请教lammps相关问题

  5. 我用lammps跑了一下,为什么结果跟博主的不一致,向博主请教

  6. 没有晶格常数怎么办

    1. 非晶体物质怎么算呢

      1. 我没有做过非晶物质的建模,可能需要使用一些软件或者编程来产生。

  7. guang

    博主你好,问下为什么到最后出现错误:cannot open input script in.thermal.expansion

    1. 是不是你把脚本的文件名改了? 找不到文件了,就会出现这个提示。

  8. 我套用你的模版去做其他元素的,可是结果与实际值差好多。是什么原因

  9. mjj

    请问博主,体积和晶格常数如何换算?多谢了

  10. zhongyanan

    请问可以由温度和体积计算体膨胀系数吗

    1. zhongyanan

      博主您好,我用你的程序算了一下金元素的,得到了十个温度和十个体积,初始温度是273K,得到的体积和温度的斜率是44.145,初始的体积是35319,那斜率和初始体积相乘得到的数据也太大了吧,这个体积单位是不是要换算一下呢?请博主指导一下

  11. minlei

    请问如何用得到的温度和体积,来得到下图的温度和晶格常数的关系曲线,

    1. 我爱搜集网博主

      将体积换算成晶格常数就可以啦~

  12. 冯闯

    用最新版的LAMMPS跑了下你的程序,跑一个循环后显示错误“unbalanced quotes in input file",这个什么原因?

    1. 冯闯

      少输入了一个引号,见笑了。

  13. charles

    代码非常简洁明了。想请教博主几个问题:

    图中的纵坐标晶格常数是输出data文件中的那个参数?是体积吗? 这个例子当中的热胀系数跟温度没有关系,如果有的材料热胀系数跟温度有关系,该怎么处理呢?
    期待您的回复。多谢了。
    1. 我爱搜集网博主
      可以通过体积换算,也可以通过修改脚本直接输出晶格常数 可以考虑将温度控制在一个较小的范围内
      1. charles

        多谢,我试试。

  14. chy

    您好,我现在在学习计算弹性模量,不知道用lammps计算,请教一下请您。

    1. 我爱搜集网博主

      请问你想要计算什么模量?
      如果是弹性常数,请参考本文。
      如果是体积模量请参考
      http://www.52souji.net/how-to-calculate-bulk-modulus-using-lammps/

      1. chy

        谢谢,我是想计算杨氏模量,和剪切模量,但是一直出问题。

  15. 付佩

    脚本写的简洁明了,但是我想问下,我运行之后,改变npt的步数,发现温度下对应的体积会发生变化,然后热膨胀系数也会变化,请问那么如何确定正确的npt步数呢?

  16. 高士才

    为什么没写质量呢?

    1. 我爱搜集网博主

      EAM势函数不需要设置质量的。
      the EAM potential files list atomic masses; thus you do not need to use the mass command to specify them.
      link: http://lammps.sandia.gov/doc/pair_eam.html

  17. 高士才

    脚本写的简洁好看,但是为什么没写质量呢?

  18. donghaipark

    为什么NVT和NPT的温度是2.5K?

    1. 我爱搜集网博主

      这个温度差不多是随便给的,是为了将体系从一个比较低的温度开始升温。