分子动力学模拟的Lammps实现

1、第五讲 分子动力学模拟的Lammps实现A.常用分子动力学模拟程序 NAMD：免费，/Research/namd/ 主要针对与生物和化学软材料体系，程序设计水平高，计算效率高。有很好的分析辅助软件VMD。 AMBER 主要针对生物体系，也适当兼容一般化学分子。有很好的内置势能模型，自定义新模型和新分子很方便，有很完善的维护网站。计算效率不高运算速度慢。 CHARMM / 主要针对生物体系，也包含部分化学体系。势能模型更新很快自定义新模型比较方便。计算效率低。 GRO

2、MACS 免费 / 主要针对生物体系，也适当照顾一般化学体系。算法好，计算效率高。界面友好，维护服务好。 TINKER 免费 /tinker/一般性分子动力学软件，对生物体系略有偏重。优点支持多种模型。仍在开发中，某些方面还不完善。 DL-POLY http:/www.cse.clrc.ac.uk/msi/software/DL_POLY/ 一般性分子模拟软件，界面友好，计算效率高。维护服务很好.A.常用分子动力学模拟程序 LAMMPS 免费 / 一般性分子模拟软

3、件。兼容当前大多数的势能模型，编程水平高，计算效率高。可以模拟软材料和固体物理系统。 Materials Explorer 立足于Windows平台的多功能分子动力学软件。拥有强大的分子动力学计算及Monte Carlo软件包，是结合应用领域来研究材料工程的有力工具。Materials Explorer可以用来研究有机物、高聚物、生物大分子、金属、陶瓷材料、半导体等晶体、非晶体、溶液，流体，液体和气体相变、膨胀、压缩系数、抗张强度、缺陷等。Materials Explorer软件中包含2Body，3Body，EAM，AMBER等63个力场可供用户选择。Materials Explorer软件拥

4、有完美的图形界面，方便使用者操作。A.常用分子动力学模拟程序B、使用L-J势模拟裂纹的扩展 L-J 无量纲化计算中的单位B、使用L-J势模拟裂纹的扩展裂纹lowerupperleftupperleftlowerxyLammps计算输入文件 # 2d LJ crack simulation(问题的基本初始化) dimension 2 #2维分子动力学模拟 boundary s s p #x,y方向不存在周期边界条件，z方向周期边界条件。 atom_style atomic #原子类型（金属） Neighbor 0.3 bin #建近邻表参数（rc+0.3）bin表示为近邻表类型。 neigh_m

5、odify delay 5 #间隔多少载荷步重新形成近邻表Lammps计算输入文件# create geometry创建初始几何构形Lattice hex 0.93#指定晶格类型（二维hex）和晶格常数Region box block 0 100 0 40 -0.25 0.25#定义一个区域create_box5 box#在指定区域建立一个simulation box,5表示原子类型的种类数create_atoms 1 box#在simulation box中创建类型为1的原子（原子位置初始化）Mass 1 1.0Mass 2 1.0Mass 3 1.0Mass 4 1.0Mass 5 1.0

6、#指定单个原子的质量Lammps计算输入文件# LJ potentials(指定原子作用势)pair_stylelj/cut 2.5#指定lj势，截断半径为2.5pair_coeff * * 1.0 1.0 2.5#指定lj势参数# define groups（便于加载）Region 1 block INF INF INF 1.25 INF INFGroup lower region 1#定义lower组（便于施加外加速度）Region 2 block INF INF 38.75 INF INF INFGroup upper region 2#定义upper组（便于施加外加速度）Group b

7、oundary union lower upper#定义总边界组Group mobile subtract all boundary#定义可动原子组（便于统计温度）Lammps计算输入文件 regionleftupper block INF 20 20 INF INF INF regionleftlower block INF 20 INF 20 INF INF groupleftupper region leftupper groupleftlower region leftlower #定义左上、左下原子组（便于指定裂纹的存在） setgroup leftupper type 2 setg

8、roup leftlower type 3 setgroup lower type 4 setgroup upper type 5 #指定原子类型（便于指定裂纹的存在）Lammps计算输入文件# initial velocities初始化速度computenew mobile temp#定义温度的计算（可动区域内统计平均）compute new2 mobile stress/atom #定义原子应力的计算（整个区域）Velocity mobile create 0.01 887723 temp new#按指定的温度（0.01）计算方法，初始化原子的速度Velocity upper set 0.

9、0 0.3 0.0#upper原子组y方向的速度为0.3Velocity mobile ramp vy 0.0 0.3 y 1.25 38.75 sum yes#mobile原子的速初始度从0到0.3线性变化# fixes施加约束fix1 all nve#nve系综的积分算法fix2 boundary setforce NULL 0.0 0.0 #边界boundary上力条件，钢化原子，便于加载！Lammps计算输入文件 # run运行计算 timestep0.003 #时间间隔步 Thermo 200 #每200步输出热动力学统计量 thermo_modify temp new #计算温度通

10、过new指示的方法计算 neigh_modify exclude type 2 3 #原子2，3之间作用取消（也就是通过不使他们在近邻表中出现实现） Dump 1 all atom 500 dump.crack #每隔500步将原子信息写入文件dump.crack Dump 2 mobile custom 500 dump2.crack tag x y z c_new22 run5000 #进行5000步的模拟C.使用EAM势模拟Ni的循环行为uuuzxy常温30K条件下金属材料模拟中Lammps的单位 # 3d metal 拉伸模拟 #模拟条件的初始化 Units metal #指定模拟中的

11、单位类型 boundary s s s #指定模拟的边界条件 atom_style atomic #指定原子类型，原子的属性 Lattice fcc 3.52 #指定材料的晶格类型和晶格常数 region box block -5 5 -5 5 -15 15 #xlo,xhi,ylo,yhi,zlo,zhi =区域box的上下限 create_box 1 box #建立只有1中原子类型的simulation box create_atoms 1 box #在simulation box中创建类型为1的原子（原子位置初始化）Lammps计算输入文件 pair_styleeam #定义原子之间相互

12、作用势 pair_coeff* * Ni_u3.eam #定义作用势的参数，通过文件输入 Neighbor 2.0 bin #定义近邻表参数(r+rc),rc=2.0埃 neigh_modify delay 5 #定义更新近邻表的加载步间隔Lammps计算输入文件 regionlower block INF INF INF INF INF -13 #生成lower区域，便于加载 regionupper block INF INF INF INF 13 INF #生成upper区域，便于加载 grouplower region lower #定义lower组 groupupper region

13、upper #定义upper组 groupboundary union lower upper #定义boundary组 groupmobile subtract all boundary #定义mobile组，便于计算统计温度（调温）Lammps计算输入文件 #定义计算和边界钢化 compute new2 all centro/atom #定义计算原子的适配度（针对所有的原子） compute new3 mobile temp 定义计算统计温度（针对可动原子） Fix 1 all nve #时间积分算法（微正则系综） fix 2 lower setforce NULL NULL 0.0 #使

14、得lower group之间没有z方向的力，z方向钢化 fix 3 upper setforce NULL NULL 0.0 #使得upper group之间没有z方向的力，z方向钢化Lammps计算输入文件 #定义控温和原子速度初始化 fix 4 mobile temp/rescale 10 30.0 30.0 5.0 1.0 #使用直接调温法调整mobile原子的速度，每10步一调 #使得目标温度为30+-1K，温度在30+-5rescale之外调 fix_modify 4 temp new3 #调温时，计算温度采用news计算 Velocity mobile create 30 8877

15、23 temp new3 #初始化mobile group中原子的初速度 #通过温度30K初始化，温度通过new3计算，887723为随机数Lammps计算输入文件 dump 2 all custom 200 d_t.txt tag x y z c_new2 #定义输出，每两百步输出一次 #tag为原子标号 #x，y，z为原子坐标 #c_new2为new2所计算的原子适配度 #这一参数可以判断原子是否是理想晶格 #或者在位错和堆垛层错上 timestep 0.001 #时间步长间隔（ps=10e-12秒） #时间步长在飞秒量级10e-15秒 #小于原子自振周期的十分之一（0.1ps/10）La

16、mmps计算输入文件Lammps计算输入文件 thermo_style custom step temp pe etotal press vol #屏幕输出梁的控制，step加载步 #temp温度，pe总势能，etotal总能量 #press为压强,vol为simulation box的体积 thermo100 # 每100步在屏幕上打印统计热动力学量 thermo_modify temp new3 #输出的温度为new3的计算值 #弛豫和循环加载 run 2000 #不加任何载荷计算2000步 velocityupper set 0 0 0.4 velocitylower set 0 0 -

17、0.4 run 10000 velocityupper set 0 0 -0.4 velocitylower set 0 0 0.4 run 20000 velocityupper set 0 0 0.4 velocitylower set 0 0 -0.4 run 10000 #循环载荷模拟Lammps计算输入文件D.使用EAM势模拟Ni的剪切行为vupperlowermobilezxyNi FCC Ni FCC 晶格晶格D.使用EAM势模拟Ni的剪切行为ZxYyz # 3d metal shear simulation 3维金属剪切模拟 Units metal #采用金属材料单位 boun

18、dary s s p #施加z方向的周期边界条件 atom_style atomic #定义原子之间的相互作用（没有键键作用） Lattice fcc 3.52 #指定晶格为fcc，晶格常数3.52埃 Region box block 0 16.0 0 10.0 0 2.828427 #定义区域box （x上下限，y上下限，z上下限） create_box 3 box #在box区域内创建包含3种原子类型的simulation boxLammps计算输入文件Lammps计算输入文件 Lattice fcc 3.52 orient x 1 0 0 orient y 0 1 1 orient z

19、& 0 -1 1 origin 0.5 0 0 #指定晶格类型、晶格常数、x,y,z坐标的晶向 # &续行符、晶胞的起始点（0.5表示半个晶格常数） create_atoms 1 box #在box区域内建立类型为1的原子 pair_styleeam #定义原子之间作用势为嵌入原子势 pair_coeff * * Ni_u3.eam #指定嵌入势参数文件Lammps计算输入文件 neighbor0.3 bin #指定原子近邻列表参数rc+r1中的r1 neigh_modify delay 5 #指定近邻表更新频率 regionlower block INF INF INF 0.

20、9 INF INF regionupper block INF INF 6.1 INF INF INF grouplower region lower groupupper region upper groupboundary union lower upper #指定边界组，用于施加边界条件 groupmobile subtract all boundary #指定可动部分组，用于计算温度Lammps计算输入文件 setgroup lower type 2 setgroup upper type 3 #设置原子类型lower为2，upper为3 Compute new3d mobile te

21、mp #设定温度计算new3d Compute new2d mobile temp/partial 0 1 1 #设定温度计算new2d，忽略x方向速度 compute new1d all stress/atom #设定原子应力计算 Velocity mobile create 300.0 5812775 temp new3d #初始化mobile区域原子的速度 #使new3d方法计算的温度为300k？ Fix 1 all nve #使系综为nve系综，同时确定了方程的积分方法 Fix 2 boundary setforce 0.0 0.0 0.0 #使边界上原子x，y，z方向的受力均为0 Fix 3 mobile temp/rescale 10 300.0 300.0 10.0 1.0 #使用直接调温法，每10步准备调一次 #如果在300+-10k内不调 #如调的