分子模拟三件套二次开发教程04第一个程序——三个 Hello WorldLAMMPS 输入脚本 / Gaussian 单点 / MaterialsScript版本声明块工具/软件LAMMPS 2Sep2026特性版/ 22Jul2025稳定版Gaussian 16 Rev. C.01 / Rev. C.02BIOVIA Materials Studio 2026Materials Scripting Editor语言/环境LAMMPS 输入脚本语法、Gaussian.gjf语法、PerlMS 2026 内置 5.40.2与 PythonMS 2026 起原生支持、bash本文目标读完能独立写出并跑通三个最小算例且能用客观判据判断跑对了而不是跑完了一句话结论三个 Hello World 的共同骨架是定义系统 → 定义交互 → 推进 → 用客观判据验收LAMMPS 侧是units ljlattice/region/create_box/create_atomspair_style lj/cut/pair_coeffvelocity/fix nverununits lj下时间步默认 0.005、邻区 skin 默认 0.3thermo的广延量默认按原子数归一化Gaussian 侧是五段输入 g16 job读job.gjf写job.log或g16 in out验收判据是日志中的Normal terminationMS 侧没有命令行 Hello World正确入口是 Materials Scripting Editor 配合 CopyScript 生成官方可用的脚本骨架。〇、本篇要解决的认知问题Q1三个 Hello World 有没有共同的骨架如果有这个骨架是哪四步Q2LAMMPS 的最小输入脚本里哪些命令是必须的顺序上有什么硬约束尤其是unitsQ3units lj下有哪几个反直觉的默认值不写就会踩Q4Gaussian 单点作业跑对了的判据是什么g16的两种命令行形态有什么区别Q5为什么 Materials Studio 的最小脚本不能是一个print hello world而要在 Scripting Editor 里起步一、机制解析1.1 三个 Hello World 的共同骨架把它们并排看会发现都是同一个四步结构① 定义系统 ② 定义交互 ③ 推进 ④ 客观判据 LAMMPS units/atom_style/lattice pair_stylepair_coeff velocityfix nve 能量漂移 / dump 帧数 region/create_box/create_atoms run Gaussian Link 0 与 Route Route 里的方法/基组 g16 进程 Normal termination Molecule specification .chk 生成 MS 打开/新建结构文档 模块里设参数 模块运行 结果表里的数值为什么这个骨架值得单独记住二次开发时你的脚本就是在自动化这四步中的某一环。谁在第 ②步定义交互上偷懒例如用错截断、漏设质量第 ④ 步的判据就会以结果看起来合理但物理错了的方式惩罚你——这是计算材料工作里最难查的一类错误。1.2 LAMMPS 最小闭环命令顺序与两条硬约束一个完整的 LJ 熔体输入脚本命令按序如下每一条都在官方手册里有自己的页面顺序命令作用硬约束1units设定单位制决定所有量的单位与两个默认值不能在盒子定义之后使用官方 Restrictions 明确2atom_style声明原子携带哪些属性atomic只有类型与坐标必须与 data/dump 的列语义一致3latticeregioncreate_boxcreate_atoms生成规则晶格、定义盒子、真正创建原子create_box一旦执行盒子就锁定4masspair_stylepair_coeff给类型赋质量与相互作用参数质量必须在动理学之前设定5velocity按温度分布赋初速度用固定的随机种子才能复现6fix nve时间积分微正则系综没有 fix 就不会发生任何运动7thermorun输出频率与推进步数步数 × 时间步长 物理时长反直觉默认值units lj相关全部取自官方units命令页时间步与 skin 是跟着单位制走的units lj下默认dt 0.005 τ、邻区skin 0.3 σ换成metal就变成dt 0.001 ps、skin 2.0 Å。所以我的温度和能量为什么震荡得这么厉害的第一诊断就是确认单位制与时间步的搭配。thermo的广延量默认按原子数归一化官方说明在units lj下thermo_style的默认行为是把所有广延量按原子数归一化输出势能按每原子给出可用thermo_modify norm改变。这意味着你看到的pe数值可能不是体系总势能而是人均值。做总量统计时这是必须处理的一步。units只能出现在盒子定义之前官方 Restrictions 明确本命令不能在read_data或create_box定义模拟盒子之后使用。把units写在脚本中段是新手高频错误而且报错信息往往不直观。经验法则任何 LAMMPS 脚本的第 1、2 行永远固定是units与atom_style。让这两行成为肌肉记忆能省掉后面 80% 的单位/属性错配类问题。1.3 Gaussian 单点五段输入 两种运行形态水分子单点是最经典的第一个 Gaussian 作业。它的 Route 段只写方法与基组如HF/6-31G(d)不加任务关键字因此只做一次能量计算。运行形态官方登记了两种形态 A作业名 g16 water → 读 water.gjf写 water.log 形态 B重定向 g16 water.gjf water.log两者等价但工程含义不同形态 A 依赖文件命名约定输入必须叫water.gjf适合批处理流水线第 08 篇会用它形态 B 允许输入来自管道适合在脚本里用 heredoc 内联输入第 02 篇的冒烟作业就是这么写的。跑对了的判据只有一条是可靠的日志里出现Normal termination。配套的三个副产品则决定后续能不能继续做多步作业.chk是否在 Link 0 里写了%Chk、日志首部的版本行用于溯源、以及能量值本身是否合理例如水的 HF/6-31G(d) 单点能应在负几十 Hartree 量级——具体数值必须与你方的实际计算比对本文不给出臆造的参考值铁律 9。1.4 MaterialsScript为什么最小脚本要在 Scripting Editor 里起步MS 的Hello World和另外两件有本质差别MS 没有公开的命令行接口也没有公开发布的脚本 API 文档。可核验的事实只有两类官方发行说明与数据表例如 MS 2026 起 MaterialsScript 原生支持 Perl 与 Python、内置 Perl 5.40.2以及客户端内的功能Materials Scripting Editor、CopyScript。因此学习路径只有一条是正确的打开 Materials Scripting Editor它自带脚本编写与调试能力在任意计算对话框里用 CopyScript 把当前设置导出成脚本——这是官方给出的从 GUI 设置到脚本的正规捷径导出的脚本里就包含该版本可用的写法以导出的写法为模板改造成批处理脚本遇到不确定的对象名与方法签名查 MS 安装后的内置在线帮助MaterialsScript 文档而不是靠搜索或猜。最佳实践把CopyScript 导出的原始脚本当成项目的第一个版本提交进仓库例如scripts/_reference/后续所有自写脚本都从它派生。这样当 MS 升级到下一个年份发行版时你能立刻用 diff 看出 API 是否变化。1.5 三者运行方式对照维度LAMMPSGaussian 16Materials Studio 2026交互式运行命令行lmp -in in.lj可执行名随构建而异不推荐交互式用g16 job不适用作业从客户端发起批处理运行shell/队列脚本包一层多份 input 批量提交NQSsubg16或用 PBS/Slurm 包装g16经 Gateway 派发到计算服务器输入形态input 脚本 data 文件.gjf可include内联外部文件文档对象 模块参数脚本化设置脚本层入口lammpsPython 模块、python命令shell 输入生成器 formchkMaterialsScriptPerl / Python成功判据正常结束且观测量守恒Normal termination结果表出现目标数值且无误标二、完整代码与逐行剖析代码 2-1in.lj——LAMMPS 的 LJ 熔体最小闭环# in.lj —— 3d Lennard-Jones 熔体教学最小闭环 # 适用版本LAMMPS 2Sep2026 / 22Jul2025 # 说明晶格常数与温度是本文为“能跑通且物理自洽”而设定的教学值不是官方推荐参数铁律 9 units lj # 单位制必须最先设定且不得晚于盒子定义官方 Restrictions atom_style atomic # 原子只带类型与坐标与代码 03-2-1 生成的 data 文件语义一致 # ---- ① 定义系统 ---- lattice fcc 1.71 # fcc 晶格常数 a1.71σ → 数密度约 4/a^3 ≈ 0.80液态量级 region box block 0 6 0 6 0 6 # 盒子 6×6×6 个晶格单元 create_box 1 box # 1 种原子类型执行后盒子尺寸被锁定 create_atoms 1 box # 在盒子内按 lattice 填满原子本例 4×6^3 864 个 # ---- ② 定义交互 ---- mass 1 1.0 # 质量必须早于任何动理学量计算 pair_style lj/cut 2.5 # LJ 截断 2.5σ势文件/表势亦可此处用解析式 pair_coeff 1 1 1.0 1.0 2.5 # ε1.0、σ1.0、截断 2.5lj 单位下均为无量纲 # ---- ③ 推进 ---- velocity all create 1.0 87287 loop geom # T*1.0 的初速度固定种子 87287 保证可复现 neighbor 0.3 bin # 与 units lj 的默认 skin 0.3 一致官方 units 页所列默认值 neigh_modify delay 0 every 20 check no # 每 20 步重建邻区表关闭每次重建时的检查以提速 fix 1 all nve # 微正则NVE积分没有 fix 原子不会运动 timestep 0.005 # 显式写出等于 units lj 的默认值把隐含默认变成可见输入 thermo 100 # 每 100 步输出一行热力学数据 thermo_style custom step temp pe ke etotal press # 注意units lj 下 thermo 的广延量默认按原子数归一化官方 units 页 # 因此上面的 pe/ke/etotal 是“每原子”值需要总量时用 thermo_modify norm no dump 1 all custom 200 dump.lj id type x y z # 每 200 步写一帧 dump_modify 1 sort id # 显式按 id 排序否则多帧之间原子顺序不确定官方 dump 页警告 # ---- ④ 客观判据 ---- run 2000 # 2000 × 0.005 τ 10 τ 的物理时长逐段剖析lattice fcc 1.71与region box block 0 6 0 6 0 6的组合给出 864 个原子、盒长 6×1.71 10.26σ、数密度约 0.80。这个数字是算出来的不是抄来的——你应当能自己复述这条换算fcc 每晶格单元 4 个原子。create_box 1 box之后再改盒子尺寸是做不到的这就是官方read_data页提到的若干系统参数在第一次定义后被锁定的同一机制。理解这一点就理解了为什么先决定盒子再填原子是唯一合理的顺序。neighbor 0.3 bin与units lj的默认 skin 0.3 一致写出来是为了让默认值可见。当你把单位制换成metal时这个数字必须跟着改skin 2.0 Å显式书写能提醒你这一点。dump_modify 1 sort id是硬要求而非优化官方dump页明确警告若 dump 不含原子 ID 又不排序可视化与分析工具无法在不同帧之间一致地分配数据。这也是第 03 篇的纪律在代码里的落地。timestep 0.005与thermo的归一化注释是本文特意留的两处反直觉默认值提示。代码 2-2water.gjf——Gaussian 水分子单点%Chkwater.chk %Mem1GB %NProcShared2 # HF/6-31G(d) water single point - first hello world 0 1 O 0.000000 0.000000 0.117300 H 0.000000 0.757000 -0.469200 H 0.000000 -0.757000 -0.469200逐段剖析第 1~4 行是 Link 0 Commands以%开头%Chk保留检查点、%Mem申请 1 GB官方默认 800 MB、%NProcShared指定共享内存并行核数。注意 Link 0 段不以空行结尾紧接着就是 Route 行。第 5 行是以#开头的 Route sectionHF/6-31G(d)只声明方法与基组因此这是一个单点能量计算。方法/基组的选择要结合体系与许可本文不做推荐。第 7 行是 Title自由文本唯一要求是与后面的段之间留空行。第 9 行0 1是电荷与多重度中性单重态水。这一行写错例如电荷不是 0会直接导致Error termination是新手最常见的一处。第 10~12 行是笛卡尔坐标本文用的是水的实验平衡几何O–H 0.9572 Å、H–O–H 104.52°换算得到的坐标你应当在纸上复核出至少一个 O–H 距离以确认自己会读这段。文件末尾的空行是官方示例的通行写法不要省。代码 2-3run_three_hello.sh——跑前两个闭环并做客观验收#!/usr/bin/env bash# 适用版本LAMMPS 2Sep2026/22Jul2025Gaussian 16 Rev. C.01/C.02# 用途跑通 LAMMPS 与 Gaussian 两个 Hello World并用客观判据验收MS 侧只能人工在客户端完成set-euopipefailROOT$PWD/hello_$(date%Y%m%d-%H%M%S)mkdir-p$ROOT/lammps$ROOT/g16echo 第一部分LAMMPS cd$ROOT/lammpscp${IN_LJ:?请用 IN_LJ 指向你的 in.lj}./in.lj# 可执行文件名随构建而异常见 lmp / lmp_mpi / lmp_serial以本机安装名为准LMP_BIN${LMP_BIN:-lmp}command-v$LMP_BIN/dev/null||{echo[FAIL] 未找到$LMP_BIN请用 LMP_BIN可执行文件路径 指定;exit1;}$LMP_BIN-inin.lj-loglog.lj-screennone# 屏幕输出关闭只看日志test-flog.lj||{echo[FAIL] 未生成 log.lj;exit1;}echo[PASS] LAMMPS 已完成运行log.lj$(wc-llog.lj)行# 判据 1dump 帧数每帧由“ITEM:”开头的若干段组成见官方 dump 页FRAMES$(grep-c^ITEM: TIMESTEPdump.lj)echo[INFO] dump 帧数 $FRAMES期望$((2000/2001)) 11# 判据 2能量漂移。取 thermo 表中的 etotal 列units lj 下为每原子值python3 -PY import re rows [] for ln in open(log.lj, encodingutf-8, errorsignore): p ln.split() # thermo 表的数据行特征首列是整数步数且共 7 列与 thermo_style custom 的字段数对应 if len(p) 7 and p[0].isdigit(): rows.append([float(v) for v in p]) if len(rows) 3: raise SystemExit([FAIL] log.lj 中的 thermo 数据行不足检查 run/thermo 设置) et [r[6] for r in rows] # step temp pe ke etotal press → 下标 6 为 etotal drift max(et) - min(et) print(f[INFO] etotal 采样 {len(rows)} 行首值 {et[0]:.6f} 末值 {et[-1]:.6f} 漂移 {drift:.6g}) print([PASS] NVE 能量漂移 1e-3每原子 if drift 1e-3 else [WARN] 漂移偏大请缩小 timestep 复核) PYechoecho 第二部分Gaussian cd$ROOT/g16cp${WATER_GJF:?请用 WATER_GJF 指向你的 water.gjf}./water.gjfecho-- 形态 Ag16 job读 job.gjf 写 job.log--g16 water||true# 成败由日志判定不依赖退出码grep-qNormal terminationwater.logecho[PASS] 形态 ANormal termination\||{echo[FAIL] 形态 A 未正常终止;tail-n20water.log;}test-fwater.chkecho[PASS] 生成 water.chk%Chk 生效||echo[WARN] 未生成 water.chkecho-- 形态 Bg16 in out --g16water.gjfwater_stdin.log||truegrep-qNormal terminationwater_stdin.logecho[PASS] 形态 BNormal termination\||echo[FAIL] 形态 B 未正常终止echoecho 第三部分Materials Studio人工 catTXT MS 侧无命令行入口请在客户端完成 1) 打开 Materials Scripting Editor 2) 在任一计算对话框用 CopyScript 导出当前设置 → 得到该版本可用的脚本写法 3) 以导出脚本为模板补全“打开结构文档 → 设置模块参数 → 运行 → 读取结果”四步。 脚本对象名与方法签名以 MS 内置的 MaterialsScript 文档为准铁律 3。 TXT逐行剖析-screen none把屏幕输出关掉、只留-log log.lj批处理里必须这样做否则日志会分散在两个地方事后排查时容易只看到一半。判据 1 用dump 帧数grep -c ^ITEM: TIMESTEP这是时间真的推进了的硬证据。期望值2000/200 1 11是可以事先算出来的这就是可验证的客观标准。判据 2 用thermo 表的 etotal 漂移NVE 系综下总能量应当守恒漂移越小越好。脚本里阈值取1e-3每原子只是教学量级真实体系应按你的物理要求与采样长度设定阈值不要照抄。p[0].isdigit()这个过滤条件很关键log.lj里既有 thermo 表也有大量提示文字。用首列是整数且列数匹配来锁定数据行比按行号取更稳跨度不同的运行会输出不同的提示行数。Gaussian 两段都写|| true再查日志成败判据是Normal termination不是退出码。同时跑形态 A 与形态 B是为了让你亲眼确认两者等价。第三部分用cat TXT直接打印人工步骤把无法自动化的部分显式写进脚本比留一句注释更不容易被后来的自己忽略。代码 2-4MaterialsScript 骨架Perl 与 Python 两版#!/usr/bin/env perl# 适用版本BIOVIA Materials Studio 2026内置 Perl 5.40.2# 重要声明铁律 3本段是“骨架 注释”。MS 的脚本 API 未公开发布# 所有对象族与方法的实际名字、参数、返回结构必须以 MS 内置的# MaterialsScript 文档Materials Scripting Editor 的帮助为准。# 请勿把下列注释里的功能描述当作可用 API 直接调用。usestrict;usewarnings;# 第 1 步导入 MaterialsScript 相关的模块。# —— 具体模块名与导入形式请照抄 Materials Scripting Editor 中# CopyScript 导出的脚本首行不同发行版可能不同。# use MaterialsScript qw(:all); # ← 本行仅为位置示意须以官方文档核对后启用# 第 2 步取到当前结构文档对象功能描述结构/集合文档族。# —— “如何取当前文档”“结构对象有哪些属性”以官方 MaterialsScript 文档为准。# 第 3 步把计算模块封装成对象并设置参数功能描述模块封装族如 Forcite。# —— 参数名与取值集合请从 CopyScript 导出的脚本中对照。# 第 4 步运行模块并读取结果功能描述结果表族。# —— 结果表的字段名与读取方式以官方文档为准。print[skeleton] 请在 Materials Scripting Editor 中按 CopyScript 导出的写法补全本脚本。\n;exit0;#!/usr/bin/env python3# -*- coding: utf-8 -*-# 适用版本BIOVIA Materials Studio 20262026 起 MaterialsScript 原生支持 Python# 重要声明铁律 3以下为骨架。MS 脚本 API 未公开发布对象名/方法签名一律# 以 MS 内置的 MaterialsScript 文档为准在 2026 之前的发行版上本骨架的# Python 前提并不成立那时只有 Perl。defmain()-int:# 第 1 步导入 MS 脚本模块模块名以官方文档与 CopyScript 导出内容为准# import 官方模块名 as ms# 第 2 步获取结构文档对象结构/集合文档族# doc 取得当前结构文档的官方写法# 第 3 步配置并运行计算模块模块封装族# module 取得模块对象的官方写法# 按官方文档设置模块参数# 第 4 步读取结果表结果表族# for row in 遍历结果表的官方写法:# print(按官方字段名取值)print([skeleton] 请以 Materials Scripting Editor 的 CopyScript 导出内容为准补全。)return0if__name____main__:raiseSystemExit(main())逐段剖析两版骨架都把官方文档边界写在文件最前面且把导入行注释掉并标注须以官方文档核对后启用这不是偷懒而是把铁律 3 变成代码里的可见约束——任何人拿到这份骨架都不会误以为里面的写法可以直接运行。骨架的四步导入 → 取文档 → 配模块 → 读结果与 §1.1 的四步骨架严格对应也与第 06 篇将要展开的 MaterialsScript 对象模型Documents / Modules / StudyTable 三个对象族一致。为什么还是坚持给骨架因为它定义了你必须向官方文档提什么问题。有骨架的人会去查四个具体问题没有骨架的人会在编辑器里瞎试。三、常见报错与排查报错 3-1LAMMPS 报Unknown command或不认识的命令。现象脚本第 N 行命令无法识别。根因三类——排版问题同一行写了多条命令而没有用续行、拼写与官方命令名不一致、或该命令来自你未启用的包如用了只有特定包才有的pair_style。解法一行只写一条命令命令名以官方手册页为准铁律 2涉及可选包时先确认构建时启用了它第 02 篇的包选择。报错 3-2LAMMPS 报units相关错误命令位置不合法。现象把units写在脚本中段后报错。根因官方 Restrictions 明确units不能在read_data或create_box定义盒子之后使用。解法把units与atom_style提到脚本最前两行如果脚本是模板生成的就在生成器里强制这两行置顶第 08 篇的模板会这么做。报错 3-3LAMMPS 报Lost atoms。现象运行若干步后提示丢原子。根因体系能量/力异常导致原子跑飞——常见于初速度与时间步不匹配units lj下dt默认 0.005若你按其他单位制的直觉给温度就会严重超速、原子初始位置重叠、或盒子过小导致周期镜像自相互作用。解法先确认单位制与timestep的搭配再检查盒子尺寸是否给非键截断留了足够空间第 03 篇代码 2-1 的margin就是为此必要时用velocity重新赋一个更温和的初速度。报错 3-4Gaussian 报Error termination日志指向电荷/多重度或 Link 0。现象作业未正常终止。根因0 1这一行与体系不符电子数为奇数却写单重态等、Link 0 行写进了 Route 段、或%Mem申请超限。解法核对电荷与多重度的奇偶性一致性确认%开头的行全部在 Route 行之前%Mem从官方默认 800 MB 起逐步上调并观察节点内存。判据始终是Normal termination。报错 3-5MS 脚本在 Scripting Editor 中报对象或方法不存在。现象脚本无法运行。根因三类——对象名/方法名与当前发行版不一致API 未公开发布只能以内置帮助为准、脚本语言栈不对把 Python 脚本用在只有 Perl 支持、即 2026 之前的发行版上、或作业需要 Gateway 而当前未连接。解法先跑 CopyScript 导出的原始脚本确认环境可用再从该模板派生API 细节一律查内置 MaterialsScript 文档铁律 3。四、动手练习练习 1LAMMPS 最小闭环跑通代码 2-1。判定dump.lj中ITEM: TIMESTEP出现 11 次2000/200 1thermo 表中temp在 1.0 附近波动etotal每原子漂移小于你在脚本里设定的阈值用上面两条判据自动给出 PASS。练习 2自己算一遍参数把lattice fcc 1.71改成两个不同值各跑run 0从输出反推数密度。判定能写出数密度 4/a³并验证改后的体系在run 0下势能有限能解释盒子大小为什么受截断 2.5σ 约束。练习 3Gaussian 两种形态用代码 2-3 分别以形态 A、形态 B 运行水分子单点。判定两份日志都出现Normal terminationwater.chk存在能写出形态 A 依赖job.gjf命名约定形态 B 允许输入来自管道因此可用于 heredoc 内联这一区别。五、小结与下一篇预告三个 Hello World 的意义不在跑通而在四步骨架 客观判据定义系统、定义交互、推进、验收。LAMMPS 侧请记住三条反直觉默认值units lj下dt0.005、skin0.3、thermo广延量按原子数归一化与units的位置约束Gaussian 侧请记住唯一可靠的判据是Normal termination、以及 Link 0 段不以空行结尾MS 侧请记住从 Scripting Editor 的 CopyScript 起步是唯一正路API 细节一律查内置帮助。第 05 篇《LAMMPS Python 接口》把第 04 篇的输入脚本升级成程序用lammps.lammps构造实例、用command/file/commands_list/commands_string/expand驱动计算、用extract_global/extract_atom/extract_compute边跑边取数、用numpy包装对象直接拿到数组并把整条链路写成一份可落盘的 CSV 结果。本篇认知问题回显FAQQ1三个 Hello World 的共同骨架是哪四步A四步为定义系统、定义交互、推进、客观判据LAMMPS 用 units 与 create_atoms、pair_coeff、fix nve 与 run、能量漂移与 dump 帧数Gaussian 用 Link 0 与 Route、分子规范、g16、Normal termination。Q2LAMMPS 最小输入脚本中哪些命令必须有有什么顺序硬约束A主线命令是 units、atom_style、lattice/region/create_box/create_atoms或 read_data、mass、pair_style 与 pair_coeff、velocity、fix nve、thermo 与 run硬约束是 units 不得晚于盒子定义create_box 后盒子与类型上限即锁定。Q3units lj 有哪几个反直觉默认值A一是时间步与 skin 随单位制变化units lj 下默认 dt0.005 τ、skin0.3 σ二是 thermo 的广延量默认按原子数归一化pe 是每原子值三是 units 不能在盒子定义之后使用。Q4Gaussian 单点作业的成败判据是什么g16 的两种命令行形态有什么区别A唯一可靠判据是日志中的 Normal termination 而非退出码g16 job 读 job.gjf 写 job.log适合批处理g16 in out 从标准输入读适合脚本内 heredoc 内联输入两者等价。Q5为什么 Materials Studio 的最小脚本不能是简单的打印语句要在 Scripting Editor 里起步A因为 MS 无公开命令行入口、脚本 API 未在发行说明中定义完整类名与方法签名只在安装后的内置帮助中正确路径是在 Materials Scripting Editor 中用 CopyScript 导出可用脚本作模板再派生。