lammps 計算實例——Cu 熱膨脹係數和熔化溫度的計算

熱膨脹係數（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) 。

舉例計算熱膨脹係數

下面仍然以銅為例，具體介紹如何使用LAMMPS計算線膨脹係數α。

輸入腳本：in.thermal.expansion

# This input script is used to calculaten# the thermal expansion of copper.n# Powered by Xianbao Duann# Email: xianbao.d@gmail.comn# Website: http://www.52souji.net/nnunits ttttttmetalnboundary tttttp p pnatom_style tttttatomicnnvariabletttttttti loop 18nvariable ttttttx equal 200+100*$innlatticettttttttfcc 3.62nregionttttttttbox block 0 8 0 8 0 8ncreate_boxttttttt1 boxncreate_atoms tttttt1 boxnnpair_style tttttt eamnpair_coeff tttttt1 1 Cu_u3.eamnnvariable ttttttN equal stepnvariable ttttttpote equal penvariable ttttttEtotal equal etotalnvariable ttttttT equal tempnvariable ttttttPress equal pressnvariable ttttttV equal volnnvelocity ttttttall create 2.5 825577 dist gaussiannntimestepttttttt0.002nthermotttttttt1000nnfixtttttttttt1 all nvt temp 2.5 2.5 0.2nrunttttttttt2000nnunfixttttttttt1nnfix ttttttt2 all npt temp 2.5 $x 0.2 iso 0 0 10nrun tttttt200000nnunfixttttttttt2nnfix tttttttextra all print 100 "${N} ${T} ${V} ${pote} ${Etotal} ${Press}" append datannfixttttttttt3 all nvt temp $x $x 0.2nruntttttttt2000nnunfixtttttttt3nunfixttttttttextrannclearnnextttttttttinjumpttttttttin.thermal.expansionn

得到的log日誌文件包含結果數據。將其中的溫度和體積複製出來列表，繪圖，如下圖。

Cu 熱膨脹係數

擬合參數

Equation y = a + b*x nWeight No Weighting nResidual Sum of Squares 8.27051E-4 nPearsons r 0.99397 nAdj. R-Square 0.98787 n Value Standard ErrornB Intercept 3.61253 4.77796E-4nB Slope 6.82149E-5 6.89885E-7n

圖中曲線的斜率是6.82149e-5，即ΔL/ΔT= 6.82149e-5 A/K。又初始的晶格常數大致為3.62A。

所以可以算出線膨脹係數α=ΔL/(L*ΔT) = 18.84*10-6 K-1。實驗結果為17.5 10-6 K-1，兩者基本吻合。

熔化溫度

將腳本的計算溫度設定到2000K，可以計算Cu的熔化溫度。

計算得到Cu的熔化溫度為1614K，Cu的實驗測定溫度為1357.77K，相差18.87%。

下面我們再來計算一下在不同溫度下體系中Cu原子的方均根位移

腳本如下

# LAMMPS Melt_Cu_tempnunits ttttttttmetalnboundary ttttttp p pnatom_style tttttatomicnvariable ttttttx index 10 500 1000 1800 2000nprint tttttttt"-------------------------------Temperature=$x K---------"nlattice tttttttfcc 3.62nregion tttttttbox block 0 8 0 8 0 8ncreate_box ttttt1 boxncreate_atoms ttttt1 boxntimestep tttttt0.01nthermo ttttttt1000npair_style tttttteamnpair_coeff tttttt1 1 Cu_u3.eamnneighbor tttttt0.6 binnneigh_modify ttttevery 5 delay 0 check yesnvelocity ttttttall create $x 825577 dist gaussian #初始化速度，按高斯分布nfix tttttttt1 all npt temp 2.0 $x 0.2 iso 0 0 10 tncompute tttttt3 all pe/atomncompute tttttt4 all ke/atomncompute tttttt5 all coord/atom cutoff 3.0nrun tttttttt50000nunfix ttttttt1nfix tttttttt2 all nvt temp $x $x 2.0 drag 0.2ndump ttttttt1 all custom 50 dump_$x.atom id xs ys zs c_3 c_4 c_5nthermo ttttttt100ncompute tttttt6 all msd com yes average yes nprinttttttttt"msd_Cu_$x.dat" #輸出msd 文件nrun tttttttt1000nclearnnext ttttttttxnjump tttttttin.melt_Cu_tempn

將數據提取，得到可視化文件。

對於鋁，LAMMPS 建立8×8×5 的FCC 晶格體系;充分弛豫後利用Nose-Hover 方法，保

持壓強為零，使體系從T=2.5K 開始加熱，直至發生熔化轉變。

以下是輸入腳本

# LAMMPS Melt_Cu or Alnunits ttttttmetal # 單位，指定為lammps 里的金屬類的單位，長度為?，能量為eV。nboundary ttttp p p # 周期性邊界條件natom_style tttatomic # 原子模式nn#variable ttti loop 150n#variable tttx equal 2.5+10*$invariable tttx equal 2.5n#variable ttttx equal 25 # 定義變數 x 為初始溫度n#lattice tttttfcc 3.61 # Cu 的晶格常數3.61nlattice tttttfcc 4.045nregion tttttbox block 0 8 0 8 0 5 # x,y,z 各方向上的晶胞重複單元數，也即區域大小ncreate_box ttt1 box # 將上述區域指定為模擬的盒子ncreate_atoms ttt1 box # 將原子按晶格填滿盒子ntimestep tttt0.01 # 步長 0.005fsn#pair_style tttteam/alloy # 選取 Cu 的EAM 勢作為模型n#pair_coeff tttt1 1 Cu_u3.eam # EAM 勢文件名稱nnpair_style ttttteam/fsnpair_coeff ttttt* * Al_mm.eam.fs Alnthermo ttttt1000 #每隔1000 步輸出熱力學結果nneighbor ttttt0.5 binnneigh_modify tttevery 5 delay 0 check yesn#velocity tttttall create $x 825577 dist gaussiannfix ttttttt1 all nvt temp $x $x 1.0 drag 0.2 #保持初始溫度，在NVT 下弛豫n#compute ttttt3 all pe/atomn#compute ttttt4 all ke/atomn#compute ttttt5 all coord/atom cutoff 3.0n#dump tttttt1 all custom 1 dump.atom id xs ys zs c_3 c_4 c_5nrun ttttttt10000 # 運行10000 步nunfix tttttt1n#fix ttttttt1 all npt $x 2000 4.00 xyz 0.0 0.0 6.0 drag 0.2 #在NPT 下加熱至2000Knfix ttttttt1 all npt temp $x 1500 4.0 couple xyz iso 0.0 0.0 6.0 drag 0.2nrun ttttttt1200000 #運行1200000 步n

運行得到log日誌文件，將數據提取。

得到單個原子平均體積與溫度的關係。在T=1110K時，體積發生躍變。根據熱力學，我們知道在該溫度下，體系發生了一級相變。