LAMMPS教程-学习指南.pdf

LAMMPS教程-学习指南.pdf

《LAMMPS教程-学习指南.pdf》由会员分享，可在线阅读，更多相关《LAMMPS教程-学习指南.pdf（110页珍藏版）》请在冰豆网上搜索。

LAMMPSLAMMPS教教程程王王延延颋颋中国科学院理论物理研究所中科院超算中心培训北京2012年9月17-18日121.1.分分子子动动力力学学模模拟拟基基本本概概念念3计算科学理论科学实验科学模型框架求解模型参数指导1.1.1.1.计算科学计算科学理论与实验科学的桥梁理论与实验科学的桥梁求解解析理论无法求解的问题模拟多体问题，得到更贴近实际体系的结果模拟实验做起来困难或不可能做的条件计算用到大量近似和模拟理想化实验条件，所以还是更贴近理论计算物理的作用计算物理的作用41.2.1.2.分子模拟分子模拟研究原子分子层面的物性研究原子分子层面的物性计算物理包括数值求解和计算机模拟两大类对物质的模拟大致分为宏观、介观、微观三个层次分子模拟方法分为MonteCarlo（MC）和MolecularDynamics（MD）两大类分子模拟对Boltzmann分布的重要性采样：

统计物理是理论基础.5解多体薛定谔方程：

H（r）=E（r）方程中包含所有原子核和电子的哈密顿量和波函数原子间相互作用由主要由电子决定计算系统最低能量构型，原子光谱等计算化学反应势垒，速率等Gaussian软件包1.3.1.3.第一性计算第一性计算电子结构方法电子结构方法6解多体薛定谔方程：

H（r）=E（r）用电子密度代替电子坐标，从而大大减少计算自由度适于进行较大体系的计算VASP软件包仍然难以做有限温度的模拟1.4.1.4.第一性计算第一性计算密度泛函理论密度泛函理论7http:

/cmm.cit.nih.gov/intro_simulation/node25.html只计算势能，不用计算力基于Boltzmann分布，只能模拟平衡态体系模拟步长可以很大1.5.1.5.蒙特卡罗（蒙特卡罗（MCMC）模拟）模拟8iijijjttFFriiittmFa1iiittttvva11iiittttrrv通过拟合实验数据或第一性计算获得的数据确定经验力场经验力场的精度直接影响模拟结果的质量并行计算环境下可模拟几百万个原子，时间尺度约几十纳秒ijF数值方法求解多体牛顿方程.1.6.1.6.全原子分子动力学模拟全原子分子动力学模拟热耦调节速度以实现恒温模拟9分子动力学模拟软件分子动力学模拟软件AMBERCHARMM10目的目的：

把整个原子作为一个质点进行模拟（全原子模拟，AtomisticSimulation或All-atomSimulation，也叫做ForceFieldMethod或MolecularMechanics），以减少计算自由度，加大可计算的体系的空间和时间尺度，简化数据处理和分析。

方法方法：

一般的做法是根据原子间相互作用的物理特性，预先设定一个有待定参数的二体或多体的相互作用的经验势的函数形式，然后根据第一性计算的数据或实验结果拟合经验势的参数。

误差误差：

因为描述体系的自由度被大大减少，全原子模型不可能重建系统的所有性质。

拟合参数时，往往选择一组最关心的物理性质进行拟合，以求误差尽量小，而放松对其它性质的要求。

所以要根据待研究的物理问题适当选取全原子模型。

1.7.1.7.分分子子建建模模11约化单位的换算约化单位的换算1.7.1.约化单位（约化单位（reducedunit）数值模拟时使用的内部单位，需要乘上常数才能对应于实际体系的真实物理单位（国际单位制，SI-units）。

1）给定四个基本物理量的单位:

长度L，质量M，时间t，电荷电量Q2）计算其它物理量：

能量温度压力质量密度数量密度介电常数B/TEk31/nL3/ML2ANQLE22*/EMLt3/PEL其中kB是Boltzmann常数，NA是Avogadro常数。

12最常用的描述原子间范德华力的经验势。

最广泛使用的是12-6LJ：

126126（）4Vrrr126126（）242rVrrrrFr惰性气体的原子间相互作用仅用LJ就基本可以完全描述。

氩原子之间的相互作用，wikipedia截断距离（截断距离（cutoffdistance）：

对于短程作用，大于cutoff的贡献是常数。

24crVrrdr三维空间中，以上积分收敛的为短程作用，发散的为长程作用。

1.7.2.Lennard-Jones势13因为金属中的价电子可以自由运动，所以一般要用多体作用描述金属体系的力场。

GlueModelEAM（EmbeddedAtomModel）12iijijijijVFrriijijjijVrUr其中rij是两个原子间的距离，是类型为和的原子之间的二体势，是类型为的原子j产生的电子电量密度在i处的值，F是一个嵌入函数，代表把类型为的原子i嵌入电子云中需要的能量。

可以用于单一金属或者合金体系。

被广泛应用于多种金属及其合金。

只适用于单一金属。

较好地平衡了表面和内部的结构和能量。

1.7.3.金属体系的力场14成键作用（成键作用（BondedInteractions）：

）：

Bonds，ValenceAngles，DihedralAngles（TorsionalAngles）,ImproperDihedralAngles非成键作用（非成键作用（NonbondedInteractions）：

）：

范德华力和静电力PandeGroupinStanfordAMBERCHARMMOPLSGROMOS1.7.4.化学和生物体系的力场15BondDihedralAngleValenceAngleImproperDihedralAngleKaptonunitEMIm+16类似于从第一性原理层面到全原子层面，粗粒化方法力图从全原子层面进一步简化到粗粒化（coarse-graining，CG）层面，以期大大提高计算的时间和空间尺度。

困难在于全原子层面上，原子间相互作用并不集中在局部。

而在第一性层面上，电子及其相互作用基本局限在相应的原子核周围。

不同的粗粒化方法着重于重建不同的物性，如结构或扩散特性等。

一些粗粒化方法假定作用势的函数形式，然后用全原子模拟的结果定参数。

另一类从结构函数（如径向分布函数）出发，反推出作用势。

我们的方法从全原子作用势出发，通过数学变换较严格地得到粗粒化力场。

1.7.5.粗粒化方法17事先做全原子模拟，得到全原子力场。

假设中心二体力的粗粒化力场形式。

最小化全原子力场和粗粒化力场之间的差值。

*W.Noid,P.Liu,Y.Wangetal.J.Chem.Phys.128,244115（2008）.从全原子力场出发严格建立粗粒化力场，不预先设定力场的函数形式。

大大减少计算自由度。

系统扩散过程加快。

可以在粗粒化层面去除一些原子自由度（如水分子）保证较好地重建结构性质理论中没有考虑可移植性！

粗粒化方法I：

MultiscaleCoarse-Graining（MS-CG）*2,113CGINNIAAICGICGIFFNN18粗粒化方法II：

EffectiveForceCoarse-Graining（EF-CG）*Y.Wang,W.Noid,P.Liu,G.A.VothPhys.Chem.Chem.Phys.11,2002（2009）.显式计算原子间作用力具有比MS-CG更好的可移植性11MNijijijrr1111MNijDijRDDDMNijijijDijDrRRRRFRRRRRFrRRRRRrR19势能面（势能面（potentialenergysurface）：

由不同构型形成的势能的集合。

系综（系综（ensemble）：

系统在给定宏观条件下所有状态的集合。

两个基本假设：

等几率原理与各态历经。

等几率原理等几率原理（principleofequalweights）：

一个热力学体系有相同的几率访问每一个微观态（注意：

不是能量的等几率！

一个能量一般会对应很多微观态）。

由等几率原理推导得出Boltzmann分布：

exp（）/jjPEQexp（）iBiQEkT其中配分函数（partitionfunction）1.8.1.8.平衡态统计物理基本概念平衡态统计物理基本概念各态历经各态历经（egodicity）：

只要系统演化无穷长时间，总有几率历经势能面上的所有点。

即在极限情况下，系综平均和时间平均是等价的。

20exp/jjjAAEQ时间平均时间平均：

分子动力学模拟（MolecularDynamics,MD）系综平均系综平均：

蒙特卡罗模拟（MonteCarlo,MC）0111limMtitiAAtdtAttM微正则系综微正则系综（MicrocanonicalEnsemble）:

NVE皆为常数。

正则系综正则系综（CanonicalEnsemble）:

NVT皆为常数。

巨正则系综巨正则系综（GrandcanonicalEnsemble）:

VT皆为常数，粒子数不固定。

等压等压-等温系综等温系综（Isobaric-IsothermalEnsemble）:

NPT皆为常数。

等张力等张力-等温系综等温系综（Isotension-IsothermalEnsemble）:

模拟盒子的形状可变。

常用系综常用系综21pp1ENiiE=H=E+pV2k112NiiiEmv动能动能B1ijijijkTNp=VdVfr常用热力学量常用热力学量温度温度势能势能21B1NiiiTmdNkv其中d是空间维数焓焓可以理解为NPT下的有效总内能压强压强lnBS=kN,V,E熵熵其中是系统的总微观状态数22BlnF=ETS=kTQHelmholtz自由能自由能G=F+pVE-TS+pVGibbs自由能自由能NPT下的自由能NVT下的自由能T,pT,VGF=NN化学势化学势23空间的连续性：

空间的连续性：

离散模型，如伊辛（Ising）模型，连续模型边界条件：

边界条件：

自由、刚性、周期周期性边界条件周期性边界条件（PeriodicBoundaryCondition,PBC）：

模拟的盒子中的粒子与无穷多的镜像中的粒子有相互作用，从而可以用103-106个粒子模拟1023个粒子的体系。

totij,1EL2i,jEnrn1.9.1.9.模拟与采样模拟与采样1.9.1.基本概念24特征长度特征长度（characteristiclength）：

某一特定物理量在空间的相关性的长度。

原则上，模拟盒子的边长应该大于所关心的物理量的特征长度。

具体操作上，可以通过变化模拟尺寸来了解有限尺度效应（finitesizeeffect）的影响。

作用势的截断距离作用势的截断距离（cutoffdistance）：

小于模拟盒子的边长的一半以避免与同一粒子的两个镜像同时作用。

有简单截断、截断平移、最小镜像法三种处理方法。

采样采样（sampling）：

本质在于在有限时间内进行重要性采样（importancesampling），即采样对系综平均贡献最大的瞬时量的子集。

一般采用均匀时间间隔的采样。

初始构型初始构型（InitialConfiguration）：

尽量接近平衡态。

一般需要一段初始的模拟过程以让初始构型达到平衡。

在这段初始的模拟过程中不采样。

需要某些参数来量化观察系统是否平衡（如液体的体积很容易平衡，势能其次，而扩散系数则较难）。

样本的相关度样本的相关度（Correlation）：

离得越近的采样样本相关度越大。

相关的样本不影响平均值，但是影响误差范围。

25Isokineticsthermostat:

每一步对速度进行修正，得到需要的温度。

Berendsenthermostat：

Andersenthermostat:

每一步从温度为T的Maxwell-Boltzmann分布中随机产生一个速度赋予被选中的粒子。

Nos-Hooverthermostat:

增加额外的自由度，在扩展的微正则系综内实现实际体系的正则系综。

1.9.2.热耦（ThermostatorHeatBath）2scaleB031,22iiiiiTNkTmvvvT12011TtTT26周期边界条件下长程力的计算方法。

把点电荷用屏蔽电荷分为实空间的短程力和长程力，短程力在实空间做截断，长程力傅里叶变换到K空间里做截断。

1.9.3.EwaldSummationByDr.JohnA.BoardatDukeUniversityParticlemeshEwald（PME）和Particle-ParticleParticle-Mesh（PPPM）算法：

把空间离散化，从而可以使用快速傅里叶变换（FFT）把计算复杂度从O（N2）降为O（Nlog（N）。

271.9.4.邻居列表（neighborlist）每两个粒子都计算距离的遍历计算复杂度是O（N2），而实际上因为短程力可以做截断，大量远距离的力为零。

为了加快计算速度，把空间划分成立方块，每个立方块的边长为截断距离。

计算粒子间距离时，对于给定粒子，只遍历最临近的27个立方块（包括自身所在的立方块）内的粒子。

计算复杂度降为O（N）。

281.9.4.势能列表为了加快计算速度，在程序开始时，把解析的力和势能的表达式在离散的坐标位置上的数值计算出来，并存入数组。

模拟过程中采用样条插值法求得精确的力和势能值。

291.10.1.10.简单的数据分析简单的数据分析1.10.1.能量数据动能Ek,势能Ep，总能量E比热容22ppB2B32NVTVEECNkkT瞬时温度kB2TEdk其中d是维数瞬时压强1ijijijPkBTdVfrr301.10.2.相关系数22（,）,0,1AABBcABcAABB时间相关函数0000220000（,）AttAttBtBtcABtAttAttBtBt时间自相关函数0000200（,）AttAttAtAtcAtAtAt311.10.3.结构特征Radialdistributionfunction（RDF）:

在距离r处找到粒子的几率。

2（）ijijiVgrrrN理想气体（）1grStructurefactor:

RDF在K空间的傅里叶变换，实验可直接观测。

20sin14krSkrgrdrkr321.10.4.扩散系数Meansquaredisplacement（MSD）:

时间t间隔内粒子运动的平均距离的平方。

22211NiirtrtrtNDiffusioncoefficient:

22rtdD其中d是维数Velocityautocorrelationfunction（VACF）:

0tvv00Dtdtvv332.2.LinuxLinux操操作作及及编编程程342.1.2.1.LinuxLinux基基本本命命令令man命令帮助信息。

ls看目录里文件的信息。

cd换目录mkdir创建目录rmdir删除目录rm删除文件mv移动或重命名cp拷贝文件cat、more、less、head、tail查看文件内容wc统计文件字符数、行数等grep在文件中查找特定字符cut、paste、sed、awk文本操作35chmod改权限；多个命令顺序执行&命令放到后台操作，等效于Z，然后bg。

管道操作，把前一个命令的输出作为后一个命令的输入。

前面命令的输出存入后面文件36#!

/bin/bash#LinuxShellScriptingTutorial1.05r3,Summer-2002#WrittenbyVivekG.Gite#Latestversioncanbefoundathttp:

/Starsfor（i=1;i=5;i+）dofor（j=1;j=i;j+）doecho-n*doneechodone脚本:

Linuxshell可以执行的批处理文件。

例子：

372.2.2.2.文本编辑器文本编辑器vivi和和emacsemacsvi缺省是命令模式，按“i”或者insert键进入插入模式，按Esc回到命令模式。

命令模式按“r”替换字符，“x”删除字符，“dd”删除行，“:

w”存盘，“:

q”退出。

emacs有菜单的编辑界面，会根据文件扩展名自动识别文件类型并给出格式模板。

有很多快捷键组合，X+S存盘，X+C退出。

382.3.2.3.分子构型绘图软件分子构型绘图软件xmakemolxmakemol和和VMDVMDxmakemol非常简单，能方便地画矢量。

只能用XYZ输入文件格式。

VMDUIUC的KlausSchulten组开发，与NAMD配套。

功能强大，有自己的脚本语言，能完成一些简单的MD模拟辅助功能。

尤其适合于针对生物体系作图。

392.4.2.4.C/C+C/C+编编译译gcc-otesttest.c-lmg+-otesttest.C更复杂的编译连接可以写Makefile完成。

402.4.2.4.并并行行作作业业管管理理命命令令#BSUB-W15#BSUBaintelmpi#BSUBqx64_3950dbg#BSUBn24#BSUBR“spanptile=12”#BSUBoout#BSUBeerrmpirun.lsf./mpi_openmp_hello.impi脚本文件bsubmpiopenmpexportOMP_NUM_THREADS=12bsub107V/mDiffusionbecomesanisotropicCagestructureplaysanessentialroleindeterminingthedynamicsMobilityandEffectiveDiffusion22eff（）（）lim2xxD/driftuvEMobility:

EffectiveDiffusion:

*RuiShiandYantingWang,tobesumittted.100EinsteinRelationBreaksDownCationdiffusesfasterAniondriftsfasterPossibleMechanisms:

FarfromequilibriumComplexfluidBDukTze101细胞膜K+，Na+Betaine+，K+，Na+细胞外细胞内BetP5.4.5.4.调节渗透压膜蛋白的工作机制调节渗透压膜蛋白的工作机制102S.Ressl,A.vanScheltinga,C.Vonrhein,V.Ott&C.Ziegler.Nature458,4752（2009）594594个氨基酸残基，个氨基酸残基，1212个跨膜片段，两个亲水端，三个完全对称的单体个跨膜片段，两个亲水端，三个完全对称的单体BetP的结构103生生物物物物理理工工作作机机制制BetaineGlycines在溶液中随机行走溶液中的K+或Na+使得BetP的C-terminal拉紧，引起BetP上部形变，从而对Betaine打开Na+结合到BetP上的结合位点上，引起BetP进一步形变，摄入Betaine104分子动力学模拟离子溶液引起的形变分子动力学模拟离子溶液引起的形变溶液中的Na+和K+离子会使C-terminal贴在另一单体上105PotentialofMeanForcePotentialofMeanForce计算计算C-terminal在离子溶液中比在纯水中受的吸引力大很多106弹弹簧簧网网络络方方法法计计算算变变形形计算C-terminal拉紧所引起的BetP整体的变形107变变形形易易于于BetaineBetaine结结合合W194与Y197会产生0.1nm的位移，使Betaine更容易结合108量化计算量化计算Na+Na+结合位点结合位点Na+的大小正好稳定嵌入结合位点，而K+太大Na+结合能:

423.8219KJ/molK+结合能:

334.1701KJ/mol结合能差值:

89.6518KJ/molkBT:

2.5KJ/molBSSE算法MP2/HF6-311G*109结结论论分子动力学模拟对于研究化学和生物分子体系的结构和动力学的微观机制是必不可少的理论工具！

110