lammps实例

Project #5 熔化与凝固：氩，铜，铝 铜和铝的熔化转变: 对于铜和铝，LAMMPS 建立 8×8×5 的 FCC 晶格体系;充分弛豫后利用 Nose-Hover 方法，保 持压强为零，使体系从 T=2.5K 开始加热，直至发生熔化转变。 下面是铜熔化的输入文件： # LAMMPS Melt_Cu or Al units metal # 单位，指定为 lammps 里的金属类的单位，长度为 Å，能量为 eV。 boundary p p p # 周期性边界条件 atom_style atomic # 原子模式 variable x equal 2.5 # 定义变量 x 为初始温度 lattice fcc 3.61 # Cu 的晶格常数 3.61 #lattice fcc 4.05 region box block 0 8 0 8 0 5 # x,y,z 各方向上的晶胞重复单元数，也即区域大小 create_box 1 box # 将上述区域指定为模拟的盒子 create_atoms 1 box # 将原子按晶格填满盒子 timestep 0.01 # 步长 0.005fs thermo 1000 #每隔 1000 步输出热力学结果 pair_style eam/alloy # 选取 Cu 的 EAM 势作为模型 pair_coeff * * jin_copper_lammps.setfl Cu # EAM 势文件名称 #pair_style eam/fs #pair_coeff * * Al_FM.eam.fs Al neighbor 0.5 bin neigh_modify every 5 delay 0 check yes #velocity all create $x 825577 dist gaussian fix 1 all nvt $x $x 1.0 drag 0.2 #保持初始温度，在 NVT下弛豫 #compute 3 all pe/atom #compute 4 all ke/atom #compute 5 all coord/atom 3.0 #dump 1 all custom 1 dump.atom id xs ys zs c_3 c_4 c_5 run 10000 # 运行 10000 步 unfix 1 fix 1 all npt $x 2000 4.00 xyz 0.0 0.0 6.0 drag 0.2 #在 NPT 下加热至 2000K #fix 1 all npt $x 1500 4.00 xyz 0.0 0.0 6.0 drag 0.2 run 1200000 #运行 1200000 步 对于铜，其熔点为 1357.77K，但在我们的模拟中其在 1609K 附近发生一级相变，比其平 衡时熔点增大了 18.7%。 对于铝，其熔点为 934.477K, 在我们的模拟中，其在 1110.0K 附近发生一级相变，单位体 积发生突变。其相对于平衡熔点增大了了 18.9%。 从铜和铝的熔化过热，以及后面所涉及的氩的熔化过热与凝固过冷，我们可以看出，在利 用分子动力学方法模拟熔化与凝固时往往会发生过热与过冷，其值基本在 10%-30%之间。 产生过热与过冷的因素基本上可以从热力学与动力学的方面阐述，例如均匀形核而导致。 铜在特定温度下的性质: 利用 LAMMPS 建立 8×8×8 的 FCC 格子，分别在 10K，500K，1000K， 1800K 和 2000K 下 保持零外压弛豫，得到在不同温度下原子运动的情况，以及不同情况下的均方根位移。 下面是 LAMMPS 的输入文件 in.melt_Cu_temp # LAMMPS Melt_Cu_temp units metal boundary p p p atom_style atomic variable x index 10 500 1000 1800 2000 print "-------------------------------Temperature=$x K---------" lattice fcc 3.62 region box block 0 8 0 8 0 8 create_box 1 box create_atoms 1 box timestep 0.01 thermo 1000 pair_style eam/alloy pair_coeff * * jin_copper_lammps.setfl Cu neighbor 0.6 bin neigh_modify every 5 delay 0 check yes velocity all create $x 825577 dist gaussian #初始化速度，按高斯分布 fix 1 all npt $x $x 2.0 xyz 0.0 0.0 6.0 drag 0.2 compute 3 all pe/atom compute 4 all ke/atom compute 5 all coord/atom 3.0 run 50000 unfix 1 fix 1 all nvt $x $x 2.0 drag 0.2 dump 1 all custom 50 dump_$x.atom id xs ys zs c_3 c_4 c_5 thermo 100 fix 2 all msd 1 msd_Cu_$x.dat #输出 msd 文件 run 1000 clear next x jump in.melt_Cu_temp 原子在不同温度下的运动 T=10K T=500K T=1000K T=2000K 均方根位移: 模拟体系中的均方根位移可以通过如下公式求得： 是指相应量的统计平均值 均方根位移的量与原子的扩散系数存在对应的关系。当体系是固态时，即体系温度处于熔点 之下时，均方根位移存在上限值；而当体系处于液态时，均方根位移呈线性关系，而且其斜 率与原子的扩散系数存在如下关系： 在 2 维体系中上式的 6 应该被 4 所取代。 下图为 T=2000K，即体系处于液态时的均方根位移图。 氩的熔化与凝固： LAMMPS 中对氩取 LJ 约化单位，其与国际单位制的转换如下: • mass = mass or m • distance = sigma, where x* = x / sigma • time = tau, where tau = t* = t (Kb T / m / sigma^2)^1/2 • energy = epsilon, where E* = E / epsilon • velocity = sigma/tau, where v* = v tau / sigma • force = epsilon/sigma, where f* = f sigma / epsilon • temperature = reduced LJ temperature, where T* = T Kb / epsilon • pressure = reduced LJ pressure, where P* = P sigma^3 / epsilon for argon: mass= 6.6362126e-26(kg) sigma=3.405e-10(m); epsilon=1.6545e-21(J); Kb=8.314/(6.02e23); tau=2.156472211e-12(s) reduced LJ temperature=119.799(K) reduced LJ velocity=157.8967679(m/s) reduced LJ pressure=41909784.02(Pa) 氩的熔化转变: 对于氩，LAMMPS 建立 8×8×5 的 FCC 晶格体系;充分弛豫后利用 Nose-Hover 方法，保持压 强为零，使体系从 T=0.01 开始加热，直至发生熔化转变。 如下是输入文件 in.melt_Ar # 2d Lennard-Jones melt units lj atom_style atomic boundary p p p #processors 2 2 1 lattice fcc 1.073 region box block 0 8 0 8 0 5 create_box 1 box create_atoms 1 box mass 1 1.0 velocity all create 0.01 872877 timestep 0.01 dump 1 all xyz 1000 melt.xyz pair_style lj/cut 2.5 pair_coeff 1 1 1.0 1.0 2.5 neighbor 0.3 bin neigh_modify every 10 delay 0 check yes thermo 1000 fix 1 all npt 0.01 0.01 1.00 xyz 0.0 0.0 1.0 drag 0.2 run 50000 unfix 1 thermo 1000 fix 1 all npt 0.01 0.85 2.0 xyz 0.0 0.0 1.0 drag 0.2 run 1000000 在 Tm=0.775 左右发生一级相变，原子单位体积发生跃变；而氩的平衡熔化温度约为 Tm=0.66， 故其在分子动力学方法模拟下过热约 20%。 氩的凝固转变： 对于氩，LAMMPS 建立 8×8×5 的 FCC 晶格体系;充分弛豫后利用 Nose-Hover 方法，保持压 强为零，使体系从 T=0.85 开始降温，发生凝固转变，直至温度降至 0 附近。 其下为输入文件 in.quench # 2d Lennard-Jones quench units lj atom_style atomic boundary p p p lattice fcc 0.851 region box block 0 8 0 8 0 5 create_box 1 box create_atoms 1 box mass 1 1.0 velocity all create 0.85 872877 timestep 0.01 pair_style lj/cut 2.5 pair_coeff 1 1 1.0 1.0 2.5 neighbor 0.3 bin neigh_modify every 10 delay 0 check yes thermo 1000 fix 1 all npt 0.85 0.85 2.0 xyz 0.0 0.0 1.00 run 50000 unfix 1 dump 1 all xyz 1000 quench.xyz fix 1 all npt 0.85 0.01 2.0 xyz 0.0 0.0 1.00 run 1000000 如下图所示，在分子动力学模拟下，氩在 T=0.44 附近发生一级相变，单位原子体积突然下 降：相对与其平衡凝固温度 Tm=0.44，约有 30%左右的过冷。 特定温度下的性质: 利用 LAMMPS 建立 8×8×8 的 FCC 格子，分别在 T=0.1，0.4，0.6，0.8 下保持零外压弛豫， 得到在不同温度下原子运动的情况, 以及不同径向分布函数和速度自相关函数。 下面是 LAMMPS 的输入文件 in.melt_Ar_temp # 2d Lennard-Jones melt_temp units lj atom_style atomic boundary p p p variable x index 0.80 0.60 0.40 0.10 lattice fcc 0.888 region box block 0 8 0 8 0 8 create_box 1 box create_atoms 1 box mass 1 1.0 print "----------------------------Temperature=$x--------------------" velocity all create $x 872877 timestep 0.01 pair_style lj/cut 2.5 pair_coeff 1 1 1.0 1.0 2.5 neighbor 0.3 bin neigh_modify every 20 delay 0 check no thermo 1000 fix 1 all npt $x $x 2.0 xyz 0.0 0.0 1.0 run 50000 unfix 1 fix 1 all nvt $x $x 1.0 compute 3 all pe/atom compute 4 all ke/atom compute 5 all coord/atom 3.0 dump 1 all custom 50 dump_$x.atom id xs ys zs c_3 c_4 c_5 run 1000 clear next x jump in.melt_temp 径向分布函数： 径向分布函数即是原子径向上的原子密度与体系总密度的比值函数（具体参见 Project_5.1）。 下图为四个温度下的氩体系的径向分布函数。可以发现在温度较低的情况下，径向分布函数 的峰比较尖锐，各个峰所对应的径向值，分别对应最近邻、次近邻等配位的位置；随着温度 的升高，径向分布函数的峰变宽，一些位置上的峰消失；直至到达液相，此时分布函数的值 不再 表 关于同志近三年现实表现材料材料类招标技术评分表图表与交易pdf视力表打印pdf用图表说话 pdf 示配位情况，而是反映了此时其它原子相对于中心原子的位置的概率分布。    速度自相关函数： 速度自相关函数的物理意义为原子某一时候的速度与其在初始时速度的函数关系，如下（具 体参见 Project_5.1） Cv(t) = < vi(0) . vi(t) >. 下图为 4 个不同温度下的速度自相关函数。可以发现如下规律：随着时间的变大，原子的速 度基本上与初始值没有相关关系（即 Cv(t)‐>0）。另外，随着温度的提高，原子的振动幅度与 频率上升，自相关函数的第一个波谷深度变小，而且到达无关速度所花费的时间减少。        1    Project #1 硅的晶格常数和体弹模量的计算   一、平衡晶格常数和内聚能 自然条件下硅为金刚石结构(dc)。计算模拟时，我们可以假定它为各种结构，fcc, bcc, sc, dc. 可以预测，模拟的dc结构的硅的体系能量最低，也即最稳定。下面我们将运用LAMMPS 来对硅的各种结构进行模拟。 定义晶格能量为Φ， 数密度为 ρ： potE N   N V   其中 Epot为势能，N 为体系总原子数，V 为体系的体积。选取 Stillinger-Weber (SW)， 以下面命令执行 lammps 运算： 其中，lmp_serial 为 lammps 命令；”<” 符号为读取符；in.Silicon 为输入文件，里 面包含运算所需要的各种数据和命令；-log 指定输出文件的名称。 可以看到屏幕上显示出 lammps 运行的信息。这个计算量很小，所以很快就结束。接下 来以如下命令来查看计算得到的数据： grep 是 linux 中一个很重要的命令，用来搜索文本，读取匹配的行并打印出来。这里是 搜索 dc.log 文件，将 @ 开头的行打印出来。如下： 晶格参数为 5.4305 埃,数密度为 0.0499540303，每个原子的能量为-4.336599609eV. @@@@ (lattice parameter, rho, energy per atom): 5.4305 0.0499540303 -4.336599609 @@@@ (lattice parameter, rho, energy per atom): 5.4306 0.04995127076 -4.336599763 @@@@ (lattice parameter, rho, energy per atom): 5.4307 0.04994851143 -4.336599879 $ grep ^@ dc.log $ lmp_serial < in.Silicon –log dc.log 2    下面具体来看刚才给的输入文件，in.Silicon. dc.log 文件中有原子总数的信息， 每个金刚石晶胞中有 8 个原子， 38 3 216  ，所以是 216 个原子。如下给出各种结构 下的体系的原子数： ^^^^^^^^^^^^^^^^^ This is a new loop ^ ^^^^^^^^^^^^^ # 注释行，标志新循环开始 Lattice spacing in x,y,z = 5.4315 5.4315 5.4315 # 晶胞大小，也即为指定的晶格参数 x Created orthogonal box = (0 0 0) to (16.2945 16.2945 16.2945) # 构造的正交盒子的坐标 1 by 1 by 1 processor grid # 一个核在运算，串行 Created 216 atoms # 构造了 216 个原子 # bulk Silicon lattice # 注释行，随便给 units metal # 单位，指定为 lammps 里的金属类的单位，长度为 Å，能量为 eV。 boundary p p p # 周期性边界条件 atom_style atomic # 原子模式 variable x index 5.4305 5.4306 5.4307 5.4308 5.4309 5.4310 5.4311 5.4312 5.4313 5.4314 5.4315 # 定义变量 x，在运行中 x 逐一取这些值。本例中为各个晶格常数。 lattice diamond $x # 晶格，指定金刚石结构的晶格，晶格常数为 x 的值 #lattice diamond 5.431 #lattice fcc 3.615 # 如果计算 fcc 结构的晶格，则将晶格常数取在 3.615 附近 #lattice bcc 3.28 # 同上 #lattice sc 2.60 # 同上 region box block 0 3 0 3 0 3 # 划定区域，x∈[0,3], y∈[0,3] z∈[0,3]，单位为晶胞 create_box 1 box # 在上面这个区域里创建一个模拟的盒子 create_atoms 1 box # 将这个盒子按晶格填满一种原子 pair_style sw # 选取 sw 势 pair_coeff * * Si.sw Si # 势文件名为 Si.sw mass 1 28 # 给定硅的质量，此处与势对应 neighbor 1.0 bin neigh_modify every 1 delay 5 check yes variable P equal pe/216 # 定义 P 为每个原子的势能。216 为原子数，pe 为体系总势能 variable r equal 216/($x*3)^3 # 定义 r 为数密度，单位体积里的原子数 timestep 0.005 # 步长为 0.005 飞秒 thermo 10 # 每 10 步在屏幕上打印一次热力学状态信息 min_style sd # 能量最小化模式，lammps 提供 sd 和 cg 可选 minimize 1.0e-12 1.0e-12 1000 1000 # 能量最小化参数设置，第一项和第二项为能量和力的 判据，指数越高，最小化程度越高，效果越好，但计算量也多。后两项为限制的最多迭代次数。 compute 3 all pe/atom # 计算每个原子的势能 compute 4 all ke/atom # 计算每个原子的动能 compute 5 all coord/atom 3.0 # 计算每个原子的近邻原子数 dump 1 all custom 1 dump.atom id xs ys zs c_3 c_4 c_5 # 输出到 dump 文件中 print "@@@@ (lattice parameter, rho, energy per atom): $x $r $P" # 打印 clear # 清除该次循环里的数据信息 next x # 跳转到下一个 x jump in.Silicon # 跳到 in.Silicon 文件，再从头读起。也即循环。 3    晶体结构类型 晶胞中的原子数 总原子数 简单立方 SC 1 27 体心立方 BCC 2 54 面心立方 FCC 4 108 金刚石 DC 8 216 表 1. 不同晶体结构中的原子数 下图是计算模拟得出的各种结构下的数密度与每个原子能量的关系图。横坐标为数密度， 以金刚石为例，ρ= 8/5.4315^3=0.049926，也即我们直接通过 grep 命令得到的第二 项值；纵坐标为每个原子的能量，为第三项值。 金刚石之外，还需计算其他结构。只需对 in.Silicon 做稍微改动： 首先，将 in.Silicon 复制成 in.fcc: 然后编辑 in. fcc 改动如下几项： 然后如下命令执行： 相应的，如下命令查看 log 文件中的数据： $ grep ^@ fcc.log $ lmp_serial < in.fcc –log fcc.log variable x index 3.610 3.611 3.612 3.613 3.614 3.615 3.616 3.617 3.618 # 设置 x 值，为平衡晶格常数 3.615 附近的值。 #lattice diamond $x # 将金刚石这项注释掉 #lattice diamond 5.431 lattice fcc 3.615 # 启用此项，fcc 结构 #lattice bcc 3.28 #lattice sc 2.60 variable P equal pe/108 # 此时原子总数为 108 variable r equal 108/($x*3)^3 # 相应修改为 108 $ gedit in.fcc $ cp in.Silicon in.fcc 4    以同样方法编辑 in.bcc, in.sc，计算不同晶格参数时的体系能量值，并绘制下图： 图 1. 不同结构下的硅的晶格能。可以看出金刚石结构对应最低能量，最为稳定 下图更为细致地画出金刚石结构中，不同晶格参数所对应的内聚能。内聚能(cohesive energy Ecoh)的定义是，最小的晶格能。由图可以得到，平衡晶格常数为 a0 = 5.431 Å， 内聚能为 Ecoh = −4.3365 eV. 图 2. 金刚石结构中的晶格能 VS 晶格参数。五阶拟合得到平衡晶格常数 5.43095(Å). 5    2 3 4 5E 104799.57911 27112.82866a 15604.99601a 7592.99351a 1134.25198a 57.82091a      二、体弹模量 我们同时可以从晶格能曲线在最低处得到体弹模量的信息。体弹模量定义为： /   dPB dV V V 和 P 分别为晶胞的体积和体系的压强。 我们已经得到了内聚能与晶格参数的函数关系，对于立方晶胞而言， 23    d M dEP dV a da 由此， 0 2 2 09  a M d EB a da 其中 a0 为平衡晶格常数，M 为体积为 3V a 的晶胞中的原子数目。 从多项式拟合，可以得到 0 2 2 a d E da , 此例中为 3.87081 eV/ Å2 (a0 =5.431 Å). 由上面 公式，计算得到硅的体弹模量为 B = 101.366 GPa. 文献中的实验数据为 99 GPa. 可视化 每一次 lammps 运行后，会生成一个 dump.atom 文件。可以通过如下命令转换为 Atomeye 可读取的 .cfg 文件： 这时当前目录中就会生成 001.cfg 之类的 .cfg 文件，然后通过 Atomeye 查看： Tab 键可以改变观察方位 上下左右箭头键可以转动原子 PgUp 和 PgDn 改变原子大小 Alt + 1 和 Alt + 2 以及 Alt + 3 改变原子颜色 先按一次9这个键，转动时就是每次以90度转动 $ A.i686 001.cfg $ lammps2cfg 1    Project #2 金属中的点缺陷：空位和间隙原子   一、空位 从晶体中移去一个原子，即可形成空位。本例将运用 LAMMPS 计算空位形成能, Ev. LAMMPS 输入文件为 in.vacancy 1) 在 fcc 结构的完整 Cu 晶体中引入一个空位 沿<100>方向构造一个 4 ×N×N×N 的晶体。N 为 input 文件中 lattice 命令指定的个 方向上的晶胞重复单元数。 2) 弛豫 当一个原子从晶体中移走之后，周围的原子将相应地调整位置以降低体系势能。为得到 稳定的构型，需要对体系进行弛豫，relaxation. LAMMPS提供两种能量最小化方式， cg 和 sd。本例中选用 sd 方式进行能量最小化。 如下是输入文件， in.vacancy: units metal # 单位为 lammps中的 metel类型 boundary p p p # 周期性边界条件 atom_style atomic # 原子模式 lattice fcc 3.61 # Cu 的晶格常数 3.61 region box block 0 6 0 6 0 6 # x,y,z 各方向上的晶胞重复单元数，也即区域大小 create_box 1 box # 将上述区域指定为模拟的盒子 create_atoms 1 box # 将原子按晶格填满盒子 pair_style eam/alloy # 选取 Cu 的 EAM 势作为模型 pair_coeff * * jin_copper_lammps.setfl Cu # EAM 势文件名称 run 0 # 运行 0 步，仅为启动 lammps 的热力学数据计算 variable E equal pe # 定义变量 E 为系统总势能 variable N equal atoms # 定义变量 N 为系统总原子数 print "the number of atoms & system energy now are $N $E" # 打印信息 region centerpoint block 3 3.05 3 3.05 3 3.05 # 指定一个原子大小的区域 delete_atoms region centerpoint # 删除这个区域的原子 min_style sd # 能量最小化模式，sd minimize 1.0e-12 1.0e-12 1000 1000 # 能量最小化参数，指数越大最小化程度越深 print "the number of atoms & system energy now are $N $E" # 打印信息 compute 3 all pe/atom # 计算每个原子的势能 compute 4 all ke/atom # 计算每个原子的动能 compute 5 all coord/atom 3.0 # 计算每个原子的近邻原子数 dump 1 all custom 1 dump.atom id xs ys zs c_3 c_4 c_5 # 将指定的各原子信息写入 dump.atom. timestep 0.005 # 步长 0.005fs run 1 # 运行 1 步 2    3) 运行 lammps 4) 计算空位形成能 空位浓度由下式给出： [n] = exp( − Fv / kBT). 其中 Fv = Ev − TSv 为形成一个空位所需要的 Helmholtz 自由能. 忽略熵 Sv, 空位浓度公式简化为 [n] = exp( −Ev / kBT). 设 E1 为完整晶体能量，含 N 个原子；E2 为弛豫后的晶体能量，含 N – 1 个原子。空 位形成能 Ev 为： 2 1 1 v NE E E N   或  2 1v cohE E N E   , 其中 Ecoh = E1 / N, 为完整晶体的内聚能。 本例中以 EAM 模型计算 4×(20×20×20)=32000 个原子的体系，得到空位形成能 Ev ~ 1.26 eV，文献中的实验值为~1.28 eV，符合较好。 另由上式计算得到，300K 温度下的空位浓度为~ 7.59×10-22 ，1350 K (Tm) 时的空位浓 度~ 2.2×10-5（文献中的实验值为~2×10-4 ）。换算时注意(1 eV/kB = 1.1604×10−4 K) 图 1. 空位处于 4×(6×6×6) 的 FCC 晶体中心， 1 06c a , 2 06c a , 3 06c a . 颜色依据原子势能标注。 3    二、间隙原子 向完整晶体中插入一个原子，即形成间隙原子。如果新插入的原子和晶体原子相同，则 为自间隙原子，self-interstitial。 与空位计算类似，我们用如下式子计算金属 Cu 中的自间隙原子形成能 Ei 2 1 1 i NE E E N   Ei 可能取决于间隙原子引入的初始位置。但对应最小的 Ei值的，应该是唯一的位置。 在 Cu, Ni 和 Pt 等 FCC 金属中，最稳定的自间隙原子构型均为[100]方向的哑铃型。 如下图： 图 2. FCC 金属中的哑铃状自间隙原子构型。右图为 LAMMPS 得到的 Cu 的构型。 晶体中引入间隙原子后，周围原子将做相应的位置调整以期达到最低能量状态。为了得 到弛豫后的构型，我们采用 LAMMPS 里的 cg 和 sd 的能量最小化方法。本例中采用 的是 sd 方法。 相对于空位，间隙原子的引入需要更大程度的弛豫。结合能量最小化方法，我们采用 NVT 或 NVE 系统的热力学平衡方法。给体系升温，让原子充分动起来，找到最稳定 的位置，得到最稳定的构型。然后淬火 quench 到 0K。最后再运用能量最小化。 The input file in.interstitial looks like 4    计算 32000 个原子的体系，得到 Cu 的自间隙原子形成能为 Ei ~ 3.1 eV. 类似上述计算，Cu 在 T=300 K 和 1350 K (Tm) 时的间隙原子浓度分别为 ~ 8.4×10-53 和 ~ 2.7×10-12. 讨论 如果体系未得到充分弛豫，可以得到各种不同的间隙原子构型，如图 3 和图 4. 1. 观察 Ei 与模拟体系大小的关系。改变盒子大小。 2. 改变间隙原子的引入位置，计算可能的间隙原子构型，并指出最稳定的间隙原子构 型和形成能。 units metal # 单位为 lammps中的 metel类型 boundary p p p # 周期性边界条件 atom_style atomic # 原子模式 lattice fcc 3.61 # Cu 的晶格常数 3.61 region box block 0 4 0 4 0 4 # x,y,z 各方向上的晶胞重复单元数，也即区域大小 create_box 1 box # 将上述区域指定为模拟的盒子 create_atoms 1 box # 将原子按晶格填满盒子 pair_style eam/alloy # 选取 Cu 的 EAM 势作为模型 pair_coeff * * jin_copper_lammps.setfl Cu # EAM 势文件名称 run 0 # 运行 0 步，仅为启动 lammps 的热力学数据计算 variable E equal pe # 定义变量 E 为系统总势能 variable N equal atoms # 定义变量 N 为系统总原子数 print "the number of atoms & system energy now are $N $E" # 打印信息 create_atoms 1 single 2.45 2.05 2.05 # 在该位置插入一个原子 min_style sd # 能量最小化模式，sd minimize 1.0e-12 1.0e-12 1000 1000 # 能量最小化参数，指数越大最小化程度越深 print "interstitial introduced, minimized: $N atoms, energy is $E" fix 1 all nvt 100 100 100 drag 0.2 # nvt 系综，原子数、体积和温度保持不变;T=100K timestep 0.005 # 步长 0.005fs run 1000 # 运行 1000 步 print "nvt performed, temperature up: $N atoms, total energy is $E" fix 1 all nvt 100 0.0001 100 drag 0.2 # nvt 系综，温度由 100K 到 0.0001K run 1000 # 运行 1000 步 print "nvt performed, temperature down: $N atoms, total energy is $E" compute 3 all pe/atom # 计算每个原子的势能 compute 4 all ke/atom # 计算每个原子的动能 compute 5 all coord/atom 3.0 # 计算每个原子的近邻原子数 dump 1 all custom 1 dump.atom id xs ys zs c_3 c_4 c_5 # 将信息写入 dump.atom min_style sd minimize 1.0e-12 1.0e-12 10000 10000 # 再次能量最小化 print "the final state: $N atoms, total energy is $E" # 打印信息 5    (a) (b) (c) 图 3. (a) FCC 晶体中的八面体（红色）和四面体（蓝色）间隙位置；(b) LAMMPS 计算所得最稳定的间隙原 子位置，恰为八面体中心；(c) LAMMPS 计算所得的另一构型，为四面体中心，体系能量比较高。 6    图 4. LAMMPS 计算得到的五种可能的间隙原子构型。 Project #3 on Molecular Dynamics Simulations FCC 金属中的面缺陷 FCC 晶体中，密排面为 {111}，它既是滑移面也是共格孪晶面。孪晶关于这个面成镜面对称。 {111} 的另一种面缺陷是层错。层错有两种，本征 (intrinsic) 和非本征 (extrinsic) 。抽出一层 原子形成本征层错，插入一层原子形成非本征层错。需要注意的是，低能量层错都可以由该面上 的剪切 (shearing) 操作得到，比如本征层错就是将某层原子上方所有的原子整体移动 1 112 . 6 < > 层错是密排面上的原子错排，层错能 (SFE) 是 材料 关于××同志的政审材料调查表环保先进个人材料国家普通话测试材料农民专业合作社注销四查四问剖析材料 的本征属性，可度量晶面滑移发生的难易。 接下来我们将运用 LAMMPS 计算 FCC 金属 Cu 和 Al 的层错能和孪晶形成能。层错和 孪晶的构型由其他代码生成，LAMMPS 通过 read_data 命令读取构型，进行计算。 例如，计算 Cu 中层错的输入文件 in.isfCu 如下： units metal boundary p p p atom_style atomic read_data isf-Cu region middle block INF INF INF INF INF 87.659 units box group middle region middle pair_style eam/alloy pair_coeff * * jin_copper_lammps.setfl Cu min_style sd minimize 1.0e-8 1.0e-8 10000 10000 compute 3 all pe/atom compute 4 all ke/atom compute 5 all coord/atom 3.0 dump 1 middle custom 1 dump.atom id xs ys zs c_3 c_4 c_5 dump_modify 1 format "%d %16.9g %16.9g %16.9g %16.9g %16.9g %g" run 0 shell mkdir isfCu_dump shell mv dump.atom isfCu_dump/dump.atom variable E equal pe print "--------------------- ISF in Cu, E = $E ------------------" (a) (b) (c) (d) 图1 (a) Cu 中的 ISF 层错; (b) Al 中的 ISF 层错；(c) Cu 中的孪晶； (d) Al 中的孪晶 层错能和孪晶形成能可由下式计算得到： 0( ) /E E Aγ = − 其中，E 为引入层错或孪晶后的体系的能量，E0 为完整晶体体系的能量，A 为层错面或孪晶面 的面积。如下是 EAM 模型计算所得的 Cu 和 Al 的层错能和孪晶形成能： Table 1 Fault energies calculated with LAMMPS Project #4 表面与界面能 铜的表面能 当物体形成表面时，表面上的原子键发生断裂，接近表面的几层原子不再如之前处于平衡状 态，从而导致能量的升高，升高的温度便是物体的表面能。 利用 LAMMPS 做出 20*20*40 fcc 的盒子，删去边缘的原子制造出一段真空层；算出此时 体系的总能量 0E ，然后从中间把盒子切成两半并移至足够远的距离，此时的体系总能量为 E final，  从而表面能：                                            0 2 final surface E E A γ = −   A为表面的面积 (100) 面与 (111) 面 如下是输入文件 in.surface_Cu_100 # LAMMPS Cu _Surface_100 units metal boundary p p p atom_style atomic lattice fcc 3.61 region box block 0 20 0 20 0 40 create_box 1 box create_atoms 1 box timestep 0.005 thermo 5 pair_style eam/alloy pair_coeff * * jin_copper_lammps.setfl Cu region boundary1 block INF INF INF INF 29.9 INF region boundary2 block INF INF INF INF INF 9.9 group boundary1 region boundary1 group boundary2 region boundary2 group boundary union boundary1 boundary2 delete_atoms group boundary neighbor 0.6 bin neigh_modify every 5 delay 0 check yes compute 3 all pe/atom compute 4 all ke/atom compute 5 all coord/atom 3.0 dump 1 all custom 100 dump.atom id xs ys zs c_3 c_4 c_5 dump_modify 1 format "%d %16.9g %16.9g %16.9g %16.9g %16.9g %g" min_style sd minimize 1.0e-30 1.0e-15 1000 10000 variable E equal pe print "---------------------- E=$E---------------------------" run 0 region down block INF INF INF INF INF 19.94 region up block INF INF INF INF 19.95 INF group up region up group down region down displace_box all z delta 0 40 units lattice remap none displac