attorney第⼀性原理分⼦动⼒学（AIMD）结果分析教程
与经典分⼦动⼒学不同，第⼀性原理分⼦动⼒学不需要提供⼒场参数，只需要提供原⼦初始结构，就能根据电⼦波函数正交化产⽣的虚拟⼒，求解⽜顿运动⽅程。在运⾏优化任务时，VASP⽣成的XDATCAR记录的是优化步骤的离⼦构型；在运⾏AIMD任务时，记录的就是运动轨迹。⽽现阶段读取XDATCAR轨迹分析性质的后处理软件并不多，能读取的兼容性也并不好。VASPKIT0.72版本之后⽀持了将XDATCAR转换成通⽤的多帧PDB⽂件的功能（504）以便可视化并进⾏后处理分析。但是并没有提供后处理分析接⼝，因此我们开发了⼀个Python脚本XDATCAR_toolkit.py，除了实现了选择⼀定范围内的帧数转换成PDB⽂件的功能，还可以提取分⼦动⼒学模拟过程中的能量，温度并做出变化趋势图。这对判断动⼒学是否平衡很有帮助。另外本脚本预留了接⼝，可以调⽤读取每⼀帧的晶格信息和原⼦坐标，以便进⾏后续扩展编程。此脚本需要安装了numpy包的python环境，以及matplotlib包以便于画图。
在得到通⽤轨迹PDB⽂件后，就可以利⽤现⽤的分⼦动⼒学后处理软件进⾏处理分析，⽐如VMD，MDtraj，MD Analysis， Pymol等。本教程将演⽰通过VMD和MD Analysis软件包分析RDF（径向分布函数）和RMSD（均⽅根偏差），前者可以⽤来分析结构性质，后者对判断结构是否稳定以及模拟是否平衡很有帮助。
将XDATCAR转换成PDB⽂件
以VASP官⽹中单个⽔分⼦的AIMD模拟为例。模拟的输⼊⽂件如下，模拟的步长是0.5 fs，模拟步数1000步，模拟时间500 fs。脚本和测试例⼦可以在我的Github仓库下载：
图1. 模拟的盒⼦
INCAR
PREC = Normal ! standard precision
ENMAX = 400 ! cutoff should be t manually
ISMEAR = 0 ; SIGMA = 0.1
ISYM = 0 ! strongly recommened for MD
IBRION = 0 ! molecular dynamics
NSW = 1000 ! 1000 steps
POTIM = 0.5 ! timestep 0.5 fs
SMASS = -3 ! No Hoover thermostat
notatall什么意思TEBEG = 2000 ; TEEND = 2000 ! temperature
exhibitionistNBANDS = 8
POSCAR
H2O _2
0.52918 ! scaling parameter
12 0 0
0 12 0
0 0 12
1 2
lect
cart
0.00 0.00 0.00 T T F
1.10 -1.43 0.00 T T F
1.10 1.43 0.00 T T F
KPOINTS
Gamma-point only
1 ! one k-point
rec ! in units of the reciprocal lattice vector
0 0 0 1 ! 3 coordinates and weight
模拟完成后将XDATCAR_toolkit.py上传到⽂件夹中（或者置于环境变量的路径⽂件夹中并赋予可执⾏权限即可直接调⽤命令
XDATCAR_toolkit.py运⾏脚本），在shell环境中运⾏以下命令:
python XDATCAR_toolkit.py -p -t 0.5 --pbc
即可将XDATCAR的全部帧也就是0~499.5 fs的轨迹转化成PDB格式。其中-p⽤于开启PDB转换功能，-t 0.5⽤于指定时间步为0.5
fs，–pbc⽤于获取基于第⼀帧演变的连续轨迹。
Now reading vasp MD energies and temperature.
Now reading vasp XDATCAR.
Total frames 1000, NpT is Fal
Finish reading XDATCAR.
Selected time-range:0.0~499.5fs
[debug] Now entering function plotfigure…
锻炼口才
运⾏完成后，将会在⽂件夹内⽣成Temperature.dat，Energy.dat，ENERGY.png和XDATCAR.pdb四个⽂件，前⾯两个分别为温度和能量随着模拟时间的的变化数据，第三个是使⽤matplotlib绘制的趋势图（如下图），最后⼀个是转换得到的轨迹PDB⽂件，可以⽤于可视化轨迹，亦可⽤于后处理分析。-
b参数⽤于指定转换从哪⼀帧开始，-e参数⽤于指定转换到哪⼀帧结束。经刘锦程博⼠建议，增加⼀个–pbc 的选项，⽤于处理周期性获取连续的轨迹。当分⼦穿过盒⼦边界时，记录真实的位置坐标（尽管它出了边界）⽽不是从盒⼦另⼀边穿⼊的ghost原⼦的坐标。这对于分析与时间相关性的量(⽐如RMSD)很有帮助。所谓连续指的是后⾯的轨迹都是从第⼀帧演变得到的真实坐标，但是并不能保证第⼀帧的分⼦是完整的，由于周期性的缘故，第⼀帧内摆放的分⼦可能分处于盒⼦两侧。李继存⽼师有篇博⽂(见下⾯链接)讲的很明⽩，可以参考。如果发现第⼀帧内分⼦不完整，可以通过添加-i 1参数将分⼦向第⼀个原⼦靠近平移以获得完整的分⼦。如果发现不理想，可以通过调整-i的参数获得完整的分⼦。
图2. 温度和系统能量的变化趋势图
RDF径向分布函数分析
得到PDB⽂件后，可以使⽤VMD，MD Analysis等分⼦动⼒学后处理软件进⾏分析。
1、使⽤VMD分析⼯具分析
打开VMD，将PDB⽂件拖⼊显⽰窗⼝，在主菜单VMD Main中选择Extensions-Analysis-Radial Pair Distribution Function g®，选择分析H(type H)在O(type O)周围的概率分布。值得注意的是分析RDF时，横坐标也就是max r不能超过盒⼦最⼩边长的⼀半，也就是得满⾜最⼩映像约定。如图4所⽰，在计算RDF时，如果max r的取值⼤于盒⼦最⼩边长的⼀半，就有可能重复算到⼀个粒⼦和它的映像粒⼦，这使得程序的周期性判断失准。将⽣成的dat⽂件的第⼀列和第⼆列作图即可得到RDF图。
图3. VMD中计算RDF
图4. 最⼩映像约定⽰意图
2、使⽤ MD Analysis分析 RDF
法语学习MD Analysis是⼀个成熟的分⼦动⼒学后处理软件，使⽤Python编写，开源。其教程不仅步骤详细还会给出背景理论知识。可以通过conda或者pip⼯具在线安装。
conda config --add channels conda-forge
conda install mdanalysis
#or
pip install --upgrade MDAnalysis
RDF分析的介绍和使⽤⽅法在⽹页(见下⾯链接)上查看。使⽤以下的脚本得到在O原⼦周围找到H原⼦的概率，并调⽤matplotlib绘制RDF 图。在1.0
Å
处出现⼀个尖峰，也就是对应了O-H键的平衡键长(0.96
Å
you are not alone 歌词)。
import MDAnalysis
import MDAnalysis.analysis.rdf
import matplotlib.pyplot as plt
u = MDAnalysis.Univer(‘XDATCAR.pdb’, permissive=True)
g1= u.lect_atoms(‘type O’)
倒装句的用法
g2= u.lect_atoms(‘type H’)
法硕辅导班
rdf = MDAnalysis.analysis.rdf.InterRDF(g1,g2,nbins=75, range=(0.0, min(u.dimensions[:3])/2.0))
rdf.run()
fig = plt.figure(figsize=(5,4))
ax = fig.add_subplot(111)
谷歌翻译插件
ax.plot(rdf.bins, rdf.rdf, ‘k-’, label=“rdf”)
ax.legend(loc=“best”)
ax.t_xlabel(r"Distance (KaTeX par error: Can't u function '\r' in math mode at position 1: \r A)")
ax.t_ylabel(r"RDF")
fig.savefig(“RDF_all.png”)
#plt.show()
图5. RDF_O_H
RMSD均⽅根偏差分析
1、VMD分析RMSD
确保使⽤了-pbc参数以获取连续的轨迹，将⽣成的XDATCAR.pdb⽂件拖⼊显⽰窗⼝。如图6右所⽰，第⼀帧内⽔的三个原⼦不在同⼀个镜像内，分⼦不完整。在进⾏RMSD分析时，尽管轨迹是连续的，但是在对齐分⼦时就会出现问题。因此在本例中需要选择第⼀个原⼦作为中⼼将分⼦平移完整，在图6左中，分⼦已经在同⼀个镜像中了。
python XDATCAR_toolkit.py -p -t 0.5 --pbc -i 1
图6. 完整和不完整的⽔分⼦
秋天服装将重新⽣成的PDB⽂件拖⼊显⽰窗⼝，在主菜单VMD Main中选择Extensions-Analysis-Analysis-RMSD Trajectory Tool，在计算RMSD前必须先做Align（对齐），这会使得每⼀帧结构进⾏平移、旋转来与参考帧的结构尽可能贴近，从⽽使得RMSD最⼩化。刘锦程提到研究⽣物法分⼦的RMSD时需要对齐操作，⽽研究⼩分⼦时不需要对齐分⼦。
