Gaussian动力学时间步长与轨迹处理
按相同物理时长比较BOMD与ADMP,读实际时间表,用ASE转换并查看轨迹。
阅读本篇目录
下载21份完整输入、实际轨迹、CSV与6张图的绘图脚本。在新目录只放本次下载的一个ZIP,然后:
unzip gaussian-dynamics-pbc-practical*.zip -d gaussian-dynamics
cd gaussian-dynamics
bash run.sh
文件名末尾可带下载版本和哈希,星号同时兼容这类文件名。默认入口使用随附实际数值生成 analysis/ 与 figures/,不启动量子化学计算。Python环境需要NumPy、Matplotlib、ASE;集群默认使用ASE 3.29.0的环境,其他电脑可设置 PYTHON=python3 bash run.sh。
重新计算前,先完成个人有权使用的Gaussian环境初始化:
yskit gaussian-init
按Gaussian输入与CPU提交申请CPU资源,加载个人环境后检查入口。下面的计算命令均应在已分配的CPU资源内执行:
command -v g16
command -v formchk
本例使用个人Gaussian16 C.01 AVX厂商预编译CPU程序,每例1核心、1GB,没有公共Gaussian module或自行编译的MPI版本。bash run.sh compute 会顺序重算全部21例;单例可用 bash compute.sh 案例名。它不会替你申请队列资源,已有同名work目录时会拒绝覆盖。
1. 先把总时长对齐#
本组水分子初始坐标和明确初始速度完全相同,分别采用0.1、0.05和0.025fs,所有轨迹都覆盖1fs。
| 步长 / fs | ADMP MaxPoints | 梯度BOMD MaxPoints | 实际帧数 |
|---|---|---|---|
| 0.1 | 10 | 11 | 11 |
| 0.05 | 20 | 21 | 21 |
| 0.025 | 40 | 41 | 41 |
在本版本实际输出中,ADMP把MaxPoints解释为传播步数;本组BOMD计数包含初始点。以原生Trajectory summary的时间表为准,不靠日志中的Step Number或相同关键词值猜测总时长。
下面为0.025fs的完整BOMD输入:
%nprocshared=1
%mem=1GB
%chk=result.chk
#p HF/STO-3G BOMD=(GradOnly,MaxPoints=41,StepSize=250,ReadVelocity) SCF=Tight NoSymm
Small native dynamics or periodic teaching probe
0 1
O 0.000000000000 0.000000000000 0.000000000000
H 0.000000000000 0.757000000000 0.586000000000
H 0.000000000000 -0.757000000000 0.586000000000
0
0.000000000000e+00 0.000000000000e+00 -2.520363643572e+12
0.000000000000e+00 1.000000000000e+13 2.000000000000e+13
0.000000000000e+00 -1.000000000000e+13 2.000000000000e+13
bash compute.sh water-bomd-dt25
grep -A 44 '^ Trajectory summary$' work/water-bomd-dt25/output.log
2. 实际短时能量变化#
最大总能量变化(约数,microHartree;1microHartree=10⁻⁶Hartree):
| 步长 / fs | 密度传播ADMP,扩展能量 | FullSCF ADMP | 梯度BOMD |
|---|---|---|---|
| 0.1 | 2.7 | 3.1 | 3.1 |
| 0.05 | 0.7 | 0.7 | 0.7 |
| 0.025 | 0.1 | 0.2 | 0.2 |

减小步长总体改善本组短时守恒,但打印精度也会影响最后几位;不要据两个很小的变化值强行拟合积分器阶数。ADMP还要区分扩展能量和核动能加势能,见ADMP能量实操。
3. 导出、转换和查看轨迹#
默认 bash run.sh 已从实际坐标生成XYZ,并用ASE写出TRAJ,保留每帧时间和适用的周期晶胞。只想手动转换一条普通水轨迹,可用:
module load ase/3.29.0
ase convert analysis/water-bomd-dt50.xyz water.traj
这个命令实际转换了21帧,坐标逐帧相同。下载到自己的图形工作站后,可以用ASE图形界面或支持XYZ/TRAJ的可视化软件打开。例如在本地图形环境执行:
ase gui water.traj
周期轨迹优先使用 analysis/h2-periodic-admp-gamma.traj。普通XYZ本身不够表达完整晶胞和周期方向;配套转换脚本从原生TV向量恢复这些信息,再检查位置和周期标记。
head -n 6 analysis/trajectories.csv
ls analysis/*traj
cat analysis/ASE-ROUNDTRIP.json
4. IRC与带时间的动力学不同#
IRC用于从过渡态沿反应路径追踪两端结构,路径点不是以上1fs轨迹中的物理时间帧。先找过渡态、检查相应虚频,再核对两端极小值;具体操作见过渡态与IRC实操。不要把IRC路径图直接改名为反应的实时分子动力学轨迹。
命令格式参见ASE转换命令;时间步长定义参见Gaussian BOMD和ADMP手册。