先说一个我见过的真实场景。一个做金属材料的师弟从网上下了一个Cu的EAM势函数文件对照论文配了一份in文件运行后界面上只有一行血红色的ERROR: Unrecognized fix style。他跑到群里问得到的回复几乎千篇一律“检查in文件。”问题是他根本不知道in文件里那几十行到底哪些能改、哪些不能改、改完会怎样。这事我遇到过太多次因为每个刚接触LAMMPS的人都会卡在同一道坎上你以为自己在填一份模拟参数表其实你是在写一份驱动分子动力学引擎运行的脚本。而这份脚本就是in文件。这篇文章不做安装教程也不聊后处理只干一件事把金属分子动力学模拟的in文件从头到尾拆开从零开始一行一行手把手教你配置。适合刚装好LAMMPS、面对空文件无从下手的新手也适合已经能跑通例子、但不清楚每个参数为什么要这么设的人。我尽量把命令背后的原理、常见的坑、我实际踩过的教训都写进去读完你能独立写出一份可运行的金属in文件并且知道怎么改材料、改温度、改拉伸条件。1. LAMMPS的in文件为什么是模拟的真正核心1.1 in文件的本质它给LAMMPS下的是“指令”而非“数据”很多新手第一次接触LAMMPS时会下意识把in文件类比成配置文件觉得里面是一堆参数跟Word里的设置项差不多。这个理解会让你在排查问题时走很大弯路。LAMMPS的本质是一个原子级别的积分推进器。你把原子坐标给它它根据力场计算出原子受力然后推进一个时间步更新坐标和速度如此循环。in文件的本质不是“参数表”而是一段按顺序执行的指令序列。每一行都是一条命令LAMMPS从上到下逐行读入、逐行执行。命令的执行顺序会影响后续的结果同样的参数换一个摆放顺序得到的体系可能就是错的。可以打个比方data文件是演员名单势函数是演员的表演规则而in文件是导演的剧本。剧本决定这场戏怎么演演员什么时候上场、用什么情绪、演多久、录制成什么格式。所以当你模拟结果不对时大部分人第一反应是怀疑势函数、怀疑data文件但真正的源头往往是剧本写错了。1.2 in文件的基本框架建模、设置、执行三阶段一份完整的金属模拟in文件不管体系多复杂骨架都逃不开下面这个顺序初始化告诉LAMMPS用哪套单位制、盒子边界条件、原子模型风格。构建模型用命令生成原子或者从data文件读入已有原子坐标。设置物理场包括势函数、邻居列表、速度初始化、控温控压方式、输出内容。执行模拟minimize做能量最小化或者run跑指定的步数。这个顺序不是随便定的。前两步决定“模拟的对象是谁”第三步决定“模拟的环境是什么”第四步才真正开始推进时间。新手常见的错误是上来就把fix nvt写在建模之前或者还没设定pair_style就急着run这些都会直接报错或者得到毫无物理意义的结果。1.3 新手容易走偏的两个认知误区第一个误区以为只要会读data文件in文件可以随便复制粘贴。跑通官方例子很容易但实际课题里你要自己建模型比如带一个空位的缺陷结构、某个特定取向的单晶、一个纳米线这时候你不会写建模段就只能一直等别人给data文件非常被动。第二个误区以为in文件里命令越少越好。现实中确实存在极简的in文件但那往往是因为很多设置用了默认值。默认值不等于最优值比如neighbor默认的skin距离、thermo默认输出频率在小测试里看不出来一放大体系就出问题。所以真正有用的in文件每一条命令都得知道自己为什么写。2. 金属体系模拟前的三件大事单位制、原子类型和势函数选择2.1 units metal不只是一行命令它决定了整个文件里所有量的“文法”金属分子动力学模拟里第一行几乎永远是units metal这一行决定了后面所有命令里物理量的单位。LAMMPS内置了lammps、real、metal、si、cgs等多套单位系统切换之后时间、长度、能量、力、压力的量纲全部跟着变。metal单位制下物理量单位长度埃Å时间皮秒ps能量电子伏特eV质量克/摩尔g/mol力eV/Å速度Å/ps温度开尔文K压力巴bar为什么金属模拟通常用metal单位最直接的原因是EAM势函数文件里的参数几乎都是按metal单位拟合的。你如果手滑把单位写成real能量单位变成kcal/mol时间单位变成飞秒fs那EAM势函数文件里的嵌入函数、对势参数全都对不上轻则结果差几个数量级重则直接报错。我在第5部分会专门讲这个坑。新手最容易忽略的一点单位不只是一个数值标签它直接影响mass命令里你填的数字。在metal单位制下mass 1 63.546里的63.546对应的是铜原子的摩尔质量g/mol。如果换成real单位虽然real的质量也是g/mol但EAM势文件不兼容如果换成si单位质量单位就变成了kg/mol你填63.546就错得离谱。2.2 atom_style与boundary你的模拟“房间”长什么样初始化阶段还有两条命令容易糊弄过去但它们是物理图像的基础atom_style atomic boundary p p patom_style atomic的意思是每个原子只有坐标、速度、原子类型这些最基本属性不带电荷。金属体系没有净电荷用atomic完全够。如果你模拟的是带部分电荷的体系比如氧化物、离子液体才需要改成charge。把atomic记成“这个体系里原子是什么物种”就行。boundary p p p定义的是模拟盒子三个方向的外墙属性。p代表周期性边界意思是这个方向上是无限重复的原子从一个边界飞出会从对面边界进来。块体金属材料用三个方向的周期性边界是最合理的。但如果你模拟的是纳米线线方向用p横截面的两个方向通常改成f固定边界或者s收缩边界否则原子会在真空层里飘走。这个选择会影响后面你建盒子时尺寸怎么取。需要特别提醒boundary必须在创建盒子之前设置盒子一旦建好边界类型就不能改。2.3 EAM势为什么金属模拟几乎离不开它接下来是势函数。金属体系最常用的组合是pair_style eam/alloy pair_coeff * * Cu_u3.eam Cu先说pair_coeff里的Cu这是在告诉LAMMPS“势函数文件里那个元素名叫做Cu的对应我的原子类型1”。如果你的体系是纯Cu这里写一个Cu就行如果是CuAl合金就要对应原子类型顺序写Cu Al。为什么金属模拟多用EAM而不是简单的Lennard-Jones你可以把金属原子想象成泡在一锅电子汤里的球。每个原子感受到的能量不仅来自它与邻居之间“一对一”的吸引排斥还取决于周围电子密度的整体厚度。EAM把总能量拆成两部分对势项负责原子核之间的短程排斥嵌入能项负责原子“泡”在局域电子密度里需要付出的能量。这个模型能同时给出合理的弹性常数、空位形成能和表面能而纯两体势很难做到。所以在做纯金属、合金的熔化、拉伸、辐照、位错模拟时EAM几乎是默认选择。EAM势函数文件一般有funcfl和setfl两种格式Cu_u3.eam这类纯金属势文件就属于funcfl格式。pair_style eam/alloy可以同时兼容这两种格式所以用起来很方便。但要注意势函数文件里的元素符号必须和pair_coeff里写的名字一致。有些势文件里面把元素写成数字编号或者在文件头部声明了元素列表这时pair_coeff的写法会不一样。新手最容易在这里栽跟头。3. 建模段lattice、region和create_atoms是怎么搭起一颗完美单晶的3.1 lattice晶格类型和晶格常数必须一次写对命令行建模是LAMMPS最方便的特性之一尤其适合从头构建单晶。核心命令是这个lattice fcc 3.615这行命令的意思是告诉LAMMPS我要一个面心立方fcc晶格晶格常数是3.615埃。这个数字就是室温下铜的实验晶格常数不是随便写的。你模拟Cu这里就必须是3.615或者你自己从实验/文献里查到的值模拟Al就要写4.05模拟Ni写3.52。很多人不清楚lattice的真正作用它不只是“声明”晶体结构它还建立了一套坐标系后续用region划分盒子、用create_atoms填充原子时都会基于这套晶格来算位置。如果你在lattice里声明fcc 3.615然后region box block 0 10 0 10 0 10那么盒子的边长10不是指10埃而是指10个晶格常数也就是36.15埃。这是LAMMPS里很特殊的设计新手刚接触时容易懵。如果你需要构建一个旋转取向的单晶比如把[100]方向转到z轴可以在lattice命令里加orient选项lattice fcc 3.615 orient x 1 1 0 orient y -1 1 0 orient z 0 0 1这样做纳米线或者倾斜晶界的时候经常会用到后面第6部分再展开。3.2 region、create_atoms和mass填充盒子前先想清楚的事有了晶格模板接下来画盒子、填原子region box block 0 20 0 10 0 10 create_box 1 box create_atoms 1 box mass 1 63.546region定义了一个名为box的长方体区域create_box根据region创建模拟盒子并声明这个体系里只有1种原子类型create_atoms按照前面lattice规定的晶格把原子填进box区域里mass给原子类型1赋予质量。这里有个容易翻车的地方盒子边长必须是晶格常数的整数倍。因为周期性边界条件下盒子里的原子排布必须满足周期性否则会在盒子界面处产生人为的应力。20×10×10个晶胞的fcc Cu盒子每个方向边长是72.3 Å、36.15 Å、36.15 Å一共20×10×10×48000个原子。这个量级对新手练习很合适跑起来也快。如果你在建模前就知道自己需要指定边长的埃数也可以写成region box block 0 72.3 0 36.15 0 36.15 units box这种情况下units box表示用真实的长度单位来定义盒子但需要自己确保盒子和晶格匹配。新手阶段我建议直接用lattice单位的写法让LAMMPS替你算尺寸省心得多。3.3 命令行建模与read_data的取舍纯LAMMPS命令行建模适合生成完美单晶、简单超胞、带空位或间隙原子的缺陷结构。但如果你要模拟多晶、复杂形状或者从实验表征里重建的构型通常需要用第三方工具生成data文件再通过read_data box.data读入。两种方式的取舍我列个表建模方式优点缺点适用场景lattice create_atoms命令简单、可重复性强、方便参数扫描只适合晶格规则体系单晶拉伸、单晶熔化、缺陷结构read_data读取data文件可以处理任意复杂构型需要额外工具生成data文件多晶、纳米线、真实微观组织我个人的建议是哪怕你以后主要用data文件建模也先把lattice、create_atoms这套命令练熟。因为在调试一个异常构型时最快的验证方法往往就是用命令行生成一个完美晶格当作对照确认力场和in文件没问题再回去查data文件的问题。4. 跑起来之前的最后冲刺速度初始化、控温控压和输出设置4.1 速度初始化给原子一个符合物理的“起跑速度”原子坐标建好了但一开始所有原子都是静止的这时候体系处于一个完全非物理的状态。做分子动力学模拟必须先给原子赋速度让体系温度对应目标温度。命令是velocity all create 300 12345 dist gaussian velocity all zero linear第一行给所有原子赋予一个符合300K麦克斯韦速度分布的初速度后面的12345是随机数种子你换成任意整数都可以但不同种子会得到不同的微观初态这会影响模拟结果所以做对比研究时要记录下来。dist gaussian的意思是速度分布采用高斯分布金属MD里这是最常用的设置。第二行zero linear必须紧跟其后作用是消除体系的整体平动动量。如果不加这一行所有原子可能会带着一个很小的整体漂移速度体系质心会持续移动温度统计也会失真。这个小细节对照实验里相当于样品在拉伸机上没夹紧一拉就开始平移数据全是废的。4.2 fix nvt与fix npt控温控压不是越复杂越好原子有了初速度接下来就要引入热浴。金属模拟里最常见的两把“枪”fix 1 all nvt temp 300 300 0.1 fix 1 all npt temp 300 300 0.1 x 0 0 0.1 y 0 0 0.1 z 0 0 0.1fix nvt把体系控制在恒温恒容状态适合做热力学性质统计或者用于等体积弛豫。fix npt在控温的基础上还允许盒子伸缩从而控制压力。命令里temp 300 300 0.1的意思分别是起始温度300K、终止温度300K、温度阻尼系数0.1 ps。阻尼系数越小控温越“硬”也就是原子会被更快拽回目标温度但也会更严重地干扰体系真实的动力学演化。一个非常关键的物理问题模拟的应该是NVT还是NPT取决于你的实验条件。块体金属在常压下的加热、冷却、拉伸通常更接近恒压环境所以用NPT更合理。而固定盒子体积、只研究某个温度下的结构演化或者模拟周期性边界下的单轴拉伸且不希望侧向收缩时才用NVT。拉伸模拟里我会用NPT在y、z方向控压具体写法在第6部分。还有个新手上路常犯的错想控制温度就同时开两个fix一个nvt一个npt结果LAMMPS直接报错。同一个原子组上同一时刻只能有一个控温fix控温和控压可以搭配但控温只能一次。4.3 dump、thermo、run让体系真正“走起来”的三板斧in文件最后通常长这样thermo 100 thermo_style custom step temp pe ke etotal press vol dump 1 all custom 1000 dump.cu id type x y z fx fy fz dump_modify 1 sort id timestep 0.001 run 5000thermo 100表示每100步在屏幕上输出一次热力学量thermo_style指定输出哪些量温度、势能、动能、总能量、压力、体积。dump负责把原子轨迹写到文件后面用OVITO可视化或者做位错分析DXA时全靠这个文件。timestep 0.001把时间步设为0.001 ps也就是1 fs这是金属MD里最常用的时间步长因为EAM势函数下原子的振动频率决定了步长不能太大否则积分不稳定。这些设置组合起来LAMMPS的执行逻辑就是按步推进每100步汇总一次热力学量每1000步输出一次原子坐标直到跑完run指定的步数。到这里一份基本的金属模拟in文件已经完整了。5. 我踩过的坑新手配置金属in文件时的5个高频翻车点5.1 单位混用EAM报错与结果荒谬都源于此这是所有LAMMPS新手遇到的第一个“玄学”问题。明明照着官方例子写的结果运行时报ERROR: Pair style eam/alloy requires metal units原因就一句话pair_style eam/alloy只接受metal单位制你在in文件里要么写成了units real要么压根没写units metalLAMMPS默认是lj单位。EAM势函数文件里的参数是在metal单位下拟合的强行在real单位下用等于拿英寸的尺子量厘米的图纸。即使不报错单位混用也会让结果变得不可理喻。比如温度一下变成几万K或者能量数量级完全对不上。排查方法很简单检查units那行是不是metal检查mass数值和单位对不对。metal单位下Cu的质量是63.546 g/mol不是0.063546 kg/mol。5.2 周期性盒子尺寸与晶格不匹配导致的假应力这个问题藏得很深。假设你写lattice fcc 3.615 region box block 0 13.83 0 6.91 0 6.91 units box这个盒子的边长13.83埃除以3.615大约是3.825个晶胞不是整数。LAMMPS不会报错它只是忠实地把原子填进region但在周期性边界下盒子边界处的原子排布会产生人为的应力导致体系一弛豫压力就异常高或者出现畸变。判断标准其实很简单fcc体系的周期盒在三个方向上必须是晶格常数的整数倍。如果你从外面读入了一个实验测得的尺寸一定要检查是否能整除。不能整除就微调盒子尺寸或者改用unit cell来构建千万别硬跑。5.3 pair_coeff路径、元素顺序和势函数文件格式我见过很多次这种报错ERROR: No pair_coeff for atom type 1原因可能是pair_coeff完全没写对也可能是势函数文件路径不对。LAMMPS查找势函数文件是相对于当前运行目录来找的。如果你的in文件和势文件不在同一目录就必须写相对路径或绝对路径pair_coeff * * /home/user/potentials/Cu_u3.eam Cu另一个问题是元素顺序。pair_coeff命令里元素名的顺序必须和create_box时定义的原子类型顺序一一对应。如果你定义了两个原子类型type 1是Cutype 2是Al那pair_coeff * * CuAl.eam Cu Al里的“Cu Al”顺序不能反反过来等于把Cu原子的参数赋给了Al。还有一点Windows和Linux系统之间拷贝in文件或势文件时容易引入额外的回车符\r导致LAMMPS解析出乱码。遇到莫名其妙的“unknown command”或“invalid pair_coeff”报错先执行一下sed -i s/\r$// in.cu把行尾清理干净再跑能省很多时间。5.4 fix的堆积一个in文件里多次run的隐藏炸弹写长一点的模拟流程时你可能要在同一个in文件里先弛豫、再拉伸、再退火。如果只是简单地在后面追加fix 1 all nvt temp 300 300 0.1 run 5000 fix 1 all nvt temp 1000 300 0.1 run 5000第二次运行会报告ERROR: Fix ID 1 is already used或者更隐蔽的情况是你第二次用了新的fix ID但忘了删掉旧的fix。旧fix会一直在那儿继续控制体系温度结果就是两套热浴同时作用在原子组上。正确做法是每次阶段结束后用unfix删除不再需要的fixfix 1 all nvt temp 300 300 0.1 run 5000 unfix 1建议把in文件按阶段分段建模、弛豫、拉伸、退火分开放到不同in文件里或者在一个in文件里明确标注fix的“生命区间”避免累积。5.5 邻居列表参数静态模拟和变形模拟不能共用一套写法金属模拟里还有一行命令看起来不起眼但直接影响计算速度和精度neighbor 0.3 bin neigh_modify every 1 delay 0 check yesneighbor 0.3 bin里的0.3是skin距离单位埃意思是构建邻居列表时多留出0.3埃的余量原子稍微移动也不会立刻破坏列表。对EAM金属体系0.3是常用默认值不建议乱改。真正的坑在第二行。neigh_modify every 1 delay 0 check yes表示每个时间步都重新检查邻居列表是否需要更新。如果是静态平衡或熔化模拟用默认的every 10 delay 5 check yes也能跑而且在性能上更友好但如果你在做拉伸变形原子位移很快邻居列表更新不及时会导致某些原子对没被计入受力结果就是应力曲线毛刺特别多甚至出现原子“穿模”。拉伸变形模拟里我通常固定用every 1 delay 0 check yes牺牲一点性能换精度。6. 完整实例一份可直接跑通的Cu单晶拉伸in文件逐一讲解6.1 完整in文件下面是我实际调试过、能直接运行的Cu单晶单轴拉伸in文件。盒子是20×10×10个晶胞一共8000个原子初始温度300K应变率1e-4 ps⁻¹拉伸方向为x。# ---------- 初始化 ---------- units metal boundary p p p atom_style atomic # ---------- 建模 ---------- lattice fcc 3.615 region box block 0 20 0 10 0 10 create_box 1 box create_atoms 1 box mass 1 63.546 # ---------- 势函数 ---------- pair_style eam/alloy pair_coeff * * Cu_u3.eam Cu # ---------- 邻居列表 ---------- neighbor 0.3 bin neigh_modify every 1 delay 0 check yes # ---------- 弛豫阶段 ---------- velocity all create 300 12345 dist gaussian velocity all zero linear fix 1 all nvt temp 300 300 0.1 thermo 100 thermo_style custom step temp pe ke etotal press vol run 5000 unfix 1 # ---------- 拉伸阶段 ---------- reset_timestep 0 variable strain equal (lx-72.3)/72.3 fix 2 all npt temp 300 300 0.05 y 0 0 0.1 z 0 0 0.1 fix 3 all deform 1 x erate 1e-4 remap x dump 1 all custom 1000 dump.cu id type x y z fx fy fz dump_modify 1 sort id thermo 100 thermo_style custom step temp v_strain press pxx pyy pzz pe ke run 200000把这个文件保存成in.cu_tensile然后在装有Cu势函数的目录下运行lmp -in in.cu_tensile就能看到LAMMPS开始跑步输出热力学量。跑完之后得到dump.cu轨迹文件用OVITO打开就能看到原子构型的演化过程。6.2 逐行拆解与作用初始化部分前几行和第2部分讲的一样不再重复。重点说几个和拉伸模拟强相关的选择。fix 2 all npt temp 300 300 0.05 y 0 0 0.1 z 0 0 0.1这个命令很关键。它的意思是温度控制在300Ky和z方向分别用恒压器把压力控制在0 bar阻尼因子0.1 ps但是x方向不施加任何压力控制。为什么x方向不控压因为x方向的变形由fix deform直接控制如果x方向再挂一个压力控制两个fix就会在x方向上抢控制权造成盒长振荡应力曲线没法看。fix 3 all deform 1 x erate 1e-4 remap x是拉伸的核心。erate 1e-4表示应变率为1e-4 ps⁻¹这意味着在1 ps内盒子x方向长度按应变率增加remap x表示每步变形时原子的x坐标也要按盒子变形的比例重新标定这样才能保证原子跟着盒子一起拉伸否则盒子变形但原子不移动等于白拉。200000步乘以0.001 ps等于200 ps对应约2%的工程应变。这个应变幅值对观察弹塑性转变起始阶段是够的但要看颈缩和断裂需要把run步数再加大或者把erate调大。注意应变率越大非平衡效应越明显应力响应会偏高所以不要为了速度快盲目调大。6.3 输出和数据处理怎么从lammps日志里得到应力应变曲线thermo_style custom step temp v_strain press pxx pyy pzz pe ke里v_strain就是前面定义的变量strain对应当前工程的应变。pxx pyy pzz是三个方向的压力分量。注意LAMMPS输出的压力是“正压为压”拉伸时x方向受拉应力pxx会是负值。大多数人做应力应变曲线时习惯把拉应力画成正的所以后处理时要取-pxx单位是bar如果换算到GPa就除以10000。更精确的做法是直接输出每个原子的应力张量然后按原子体积平均。但对新手阶段直接用thermo_style里的pxx做初步的应力应变曲线完全可以等需要发表级别的结果再去用compute stress/atom。数据处理时把终端里输出的step和pxx两列抓出来用任意脚本画图就行。我习惯在跑完后加一行variable stress equal -pxx/10000然后在thermo_style里输出v_stress这样直接得到单位是GPa的拉应力省去了一步后处理。6.4 从Cu换成其他金属或添加缺陷这份in文件换材料异常简单。比如换成Al把lattice改成fcc 4.05mass改成26.98pair_coeff里的元素名和势函数文件换成Al的比如Al03.eam或Al_u3.eam。换成Ni也一样晶格常数3.52质量58.69。这里唯一要注意的是势函数文件里的元素名必须和pair_coeff里写的一致。要加空位也很简单。比如想在盒子中心去掉一个原子region hole sphere 36.15 18.08 18.08 2.0 delete_atoms region hole这样就在中心半径2埃的球形区域内删除了原子形成空位团。如果想建单个空位就把球半径设小一点例如1.0甚至0.5保证只删掉中心的一个原子即可。6.5 我个人实际操作中的经验如果让我给这份经验做个总结我会说in文件不是写一次就完的它是一张可以反复修改的实验计划表。我自己习惯把模拟拆成三个文件保存build.in只负责建模并把原子写进data文件equil.in负责弛豫prod.in负责拉伸、加热等生产模拟。这样每改一个阶段不用把前面全重跑一遍特别是建了一个复杂的缺陷模型之后重跑建模段的代价很高。这个习惯帮我少走了很多弯路也建议你从第一次正式做金属模拟开始就养成。