分子模擬-水分子擴散系數(shù)計算教程LAMMPS

1.基礎理論擴散系數(shù)是衡量物質擴散程度的物理量，是表征物質傳質的重要參數(shù)。自擴散系數(shù)是一個張量，在均相系統(tǒng)中通常取張量的對角線元素1D=_(D+D+D)3xxyyzzxyz表示方向，且DDD三者相等。通常非對角線元素為0。xxyyzz分子模擬中，計算自擴散系數(shù)的方法比較常用的有Einsteinequation和Green-Kuboequation。此教程使用Einsteinequation，即D二__1_limd｛丄Xr(t)-r(0)2xx 2Dimt*dt\N 丨， 1\ i=1其中Dim是維度，r(t)即為i粒子t時刻的位置，N為粒子總數(shù)，尖括號表示系綜平均。其中MSD其中MSD=穆習叩)-ri(o)i2表示粒子的均方位移(MeanSquareDisplacement)。i=1模型構建模型構建使用的工具為:moltemplate一般在lammps/tools/目錄下已有安裝(參考了官網上的范例)需要的文件為spce.lt，包含水分子模型的參數(shù)°System.lt,指定水分子的排列方式。InmoltemplateeachmoleculetypeisstoredinafilewhoseformatmimicstheformatofaLAMMPS"data"file,with"Atoms","Bonds",and"Angles"sections.Oncedefined,thesemoleculescanbecopiedandmovedtobuildlargersystems.#"SPCE"("extendedsimplepointcharge")isoneofseveralpopularwatermodels(SeeBerendsenH.J.Cetal.,J.Phys.Chem.Vol91pp.6269-6271,1987fordetails)#file"spce.lt"##H1H2#\/#OSPCE{# Definitionofthe"SPCE"(watermoleculetype)

#The"DataAtoms"sectionisalistofatomattributes(type,charge,x,y,z)write("DataAtoms"){$atom:o$molwrite("DataAtoms"){$atom:o$mol$atom:h1$mol$atom:h2$mol@atom:O-0.8476@atom:O-0.8476@atom:H0.42380.81649040.000000.5773590@atom:H0.4238-0.81649040.000000.5773590}#Note:LAMMPSexpectsanintegerinthe2ndcolumn(theMolecule-IDnumber).# Ifweput"$mol"there,moltemplatewillgeneratethisintegerforyou#Alistofthebondsinthemolecule:#BondIDBondTypeAtomID1AtomID2write("DataBonds"){$bond:OH1@bond:OH$atom:o$atom:h1$bond:OH2@bond:OH$atom:o$atom:h2}#Alistoftheanglesinthemolecule:#AngleIDAngleTypeAtomID1AtomID2AtomID3write("DataAngles"){$angle:HOH@angle:HOH$atom:h1$atom:o$atom:h2}#Themassofeachatomtype:write_once("DataMasses"){@atom:O15.9994@atom:H1.008}}#endofdefinitionof"SPCE"watermoleculetype以及system.lt,創(chuàng)建模擬所需要的水分子盒子。import"spce.lt"#<-ThisdefinestheSPCEwatermolecule.#Periodicboundaryconditions:write_once("DataBoundary"){0.034.50xloxhi

0.034.50yloyhi0.034.50zlozhi}#Thenextcommandcreates1000(10x10x10)copiesoftheSPCEwatermolecule#andarrangesthemina(ratherdense)cubiclatticewith#spacing3.45Angstroms.(Thepressuremustbeequilibratedlater.)wat=newSPCE[10].move(0.00,0.00,3.45)[10].move(0.00,3.45,0.01)[10].move(3.45,0.01,0.01)運行:$moltemplate.shsystem.lt即可得到水分子的data文件：system.data。參數(shù)設置以下為in文件:unitsreal#指定計算用的單位atom_stylefull#data文件的數(shù)據(jù)類型bond_styleharmonic#計算鍵能的方式angle_styleharmonic#計算鍵角相互作用的方式pair_stylelj/charmm/coul/long9.010.010.0#非鍵相互作用的類型pair_modifymixarithmetic#非鍵相互作用的混合方法此范例所使用的LAMMPS版本為11Aug2017。運行LAMMPS所需的最基本文件為in文件，在in文件中也可以創(chuàng)建簡單的結構模型。此范例所需要的文件為in以及包含水分子模型的data文件。#計算長程庫侖力的方式read_datasystem.data#讀入結構文件bond_coeffread_datasystem.data#讀入結構文件bond_coeff11000.01.0#鍵參數(shù)angle_coeff11000.0109.47#鍵角參數(shù)pair_coeff110.15533.166#O-O非鍵相互作用pair_coeff220.02.058#H-H非鍵相互作用groupgSPCEtype12#將1，2原子設定為組fixfShakeSPCEgSPCEshake0.000110100b1a1#SHAKE算法令水分子為剛性kspace_stylepppm0.0001#時間步長#輸出格式dump1allcustom1000traj.lammpstrjidmoltypeqxyzdump_modify1sortid #按原子id排序fixfxnptallnpttemp300.0300.0100.0iso1.01.01000.0drag1.0 #npt系綜run50000unfixfxnpt #解除nptfixfxnvtallnvttemp300.0300.0100.0 #nvt系綜nose-hoover熱浴thermo100 #log文件中輸出信息的頻率computemsdgSPCEmsd #計算msdfixfxmsdgSPCEave/time100011000c_msd[1]c_msd[2]c_msd[3]c_msd[4]filemsd.txt #輸出msdrun600000 #模擬時長可自行調整，此處為0.6納秒運行模擬將data文件，in文件放置在同一文件夾下，運行模擬。此處為1000個水分子，耗時較長。推薦使用并行計算。當然你也可以自行調整水分子的數(shù)目以及模擬時長。計算結果模擬得到水分子的MSD數(shù)據(jù)（msd.txt）,將其作圖。XXyyZZ01000-1S00-200-Slopexx=4.98217e-4Slopeyy=5.01173e-4Slopezz=5.35622e-4模擬得到水分子的MSD數(shù)據(jù)（msd.txt）,將其作圖。XXyyZZ01000-1S00-200-Slopexx=4.98217e-4Slopeyy=5.01173e-4Slopezz=5.35622e-41D00002000003000004000005DOOOO600000700000Time(fs)600-400-Slope=0.00154▼totalFig.水分子的MSD隨時間變化因為xx,yy,zz三個方向擴散系數(shù)相同，因此我們取其斜率的平均值5?1167e-4。此時維度為1,由E