与经典分子动力学不同，第一性原理分子动力学不需要提供力场参数，只需要提供原子初始结构，就能根据电子波函数正交化产生的虚拟力，求解牛顿运动方程。在运行优化任务时，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仓库下载：
https://github.com/tamaswells/VASP_script/tree/master/XDATCAR_tookit
图1. 模拟的盒子
INCAR
PREC = Normal ! standard precision
ENMAX = 400 ! cutoff should be set 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 ! Nose Hoover thermostat
TEBEG = 2000 ; TEEND = 2000 ! temperature
NBANDS = 8
POSCAR
H2O _2
0.52918 ! scaling parameter
12 0 0
0 12 0
0 0 12
1 2
select
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 False
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的参数获得完整的分子。
http://jerkwin.github.io/2016/05/31/GROMACS%E8%BD%A8%E8%BF%B9%E5%91%A8%E6%9C%9F%E6%80%A7%E8%BE%B9%E7%95%8C%E6%9D%A1%E4%BB%B6%E7%9A%84%E5%A4%84%E7%90%86/
图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
Å
)。
https://www.mdanalysis.org/docs/documentation_pages/analysis/rdf.html#radial-distribution-functions-mdanalysis-analysis-rdf
import MDAnalysis
import MDAnalysis.analysis.rdf
import matplotlib.pyplot as plt
u = MDAnalysis.Universe(‘XDATCAR.pdb’, permissive=True)
g1= u.select_atoms(‘type O’)
g2= u.select_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.set_xlabel(r"Distance (
KaTeX parse error: Can't use function '\r' in math mode at position 1: \̲r̲ ̲A
)")
ax.set_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时需要对齐操作，而研究小分子时不需要对齐分子。
图7. VMD中计算RMSD
把左上角文本框里的默认的Protein改成all（代表所有原子都纳入考虑），然后把noh复选框的勾去掉（否则将忽略氢原子）。然后点右上角的ALIGN按钮，此时所有帧的结构就已经对齐了。本例中演示以模拟的第一帧为参考，分析氧原子位置的均方根偏差。因此在Reference mol那里选top作为参考结构，左上角文本框由all改为type O（代表计算O原子的RMSD），然后勾上Plot复选框，最后点击RMSD按钮即可得到O原子的RMSD图。在File菜单栏可以选择导出dat数据。
图8. VMD中未对齐轨迹计算的RMSD
图9. VMD中对齐了轨迹后计算的RMSD
2、使用 MD Analysis分析 RMSD
RMSD分析的介绍和使用方法在网页(见下面链接)上查看。使用以下的脚本可以分别得到所有原子，氢原子，氧原子的RMSD，并调用matplotlib绘制RMSD图。网页中有一段话(Note If you use trajectory data from simulations performed under periodic boundary conditions then you must make your molecules whole before performing RMSD calculations so that the centers of mass of the selected and reference structure are properly superimposed)也就是在计算RMSD的时候选择的分子必须是完整的，不能分处于盒子的两边。这与我们之前的描述是一致的。MD Analysis默认对齐了分子。
https://www.mdanalysis.org/docs/documentation_pages/analysis/rms.html?highlight=average
使用以下脚本可以绘制对齐了轨迹后所有原子，氧原子和氢原子的RMSD。
import MDAnalysis
import MDAnalysis.analysis.rms
import matplotlib.pyplot as plt
u = MDAnalysis.Universe(‘XDATCAR.pdb’, permissive=True)
ref = MDAnalysis.Universe(‘XDATCAR.pdb’, permissive=True) # reference (with the default ref_frame=0)
ref.trajectory[0] #use first frame as reference
R = MDAnalysis.analysis.rms.RMSD(u, ref,
select=“all”, # superimpose on whole backbone of all atoms # align based on all atoms
groupselections=[“type H”,“type O”],
filename=“rmsd_all.dat”,center=True)#, # CORE
timestep=0.0005 #0.5fs from fs to ps as Reader has no dt information, set to 1.0 ps
R.run()
rmsd = R.rmsd.T # transpose makes it easier for plotting
time = rmsd[1]*timestep
fig = plt.figure(figsize=(5,4))
ax = fig.add_subplot(111)
ax.plot(time, rmsd[2], ‘k-’, label=“all”)
ax.plot(time, rmsd[3], ‘r–’, label=“type H”)
ax.plot(time, rmsd[4], ‘b–’, label=“type O”)
ax.legend(loc=“best”)
ax.set_xlabel(“time (ps)”)
ax.set_ylabel(r"RMSD (
KaTeX parse error: Can't use function '\r' in math mode at position 1: \̲r̲ ̲A
)")
fig.savefig(“rmsd_md_analysis.png”)
图10. MD Analysis计算的对齐了轨迹后计算的所有原子，H原子和O原子的RMSD
原文链接：
https://mp.weixin.qq.com/s?src=11×tamp=1589009633&ver=2327&signature=1Yk-QhQFJ6hyc6uoGEig1CS*X3ybv24G15dcs*26lxP27eCuA-JrcuqgTJTmM5TrDCuDzLAFFm65FKf8l1O2f9NJ0H0LLA3oM4hdkp4MNa6GeotfRV5B40Dl7xh-z1nj&new=1
VASP_XDATCAR
用于从
分子
动力学
解析VASP XDATCAR的python类。
XDATCAR文件包含
分子
动力学
运行期间的轨迹，即在每个时间步上所有原子的位置。 根据这些信息，我们可以计算出以下物理量
系统随时间变化的温度
速度自相关函数（VAF）和声子密度
对相关函数（PCF）
Set NBLOCK = 1 in the INCAR so that all the configuration in the MD run is
wrtten to XDATCAR.
The element mass (POMASS) and the MD time step (POTIM) is read from OUTCAR.
随着
分子
动力学
模拟在物理、化学和材料科学中的日益广泛应用，处理和
分析
模拟产生的大量数据成为了一个关键步骤。为此，应运而生，为科研工作者提供了一个强大而便捷的MATLAB工具包，能够从VASP输出文件中提取轨迹，并以VESTA可读的格式展现核密度。本文将详细介绍如何使用
AIMD
_Toolkit进行相关操作。
AIMD
_Toolkit无疑是
分子
动力学
研究中的一大利器。本文为您提供了从安装、使用到深入理解的全面
教程
，希望能帮助您在科研之路上更加顺利。
随着计算机网络的迅速发展，高效的网络拥塞控制算法对于保证网络传输质量至关重要。
AIMD
（Additive Increase Multiplicative Decrease）算法作为一种经典的网络拥塞控制算法，在当前网络技术中得到广泛应用。本文将深入探讨
AIMD
算法的原理、工作过程以及优势，并通过实例和案例
分析
，揭示其在现代计算机网络中的重要地位。
VASP 6.3 正式发布了机器学习力场功能，在IBRION=0 进行从头算
分子
动力学
模拟（
AIMD
）时只需在INCAR中控制两个标签ML_LMLFF ；ML_ISTART便能在跑
分子
动力学
时选择是否引入即时（on-the-fly）机器学习力场。官方宣称在从头算
分子
动力学
模拟中引入on-the-fly机器学习力场将会显著降低模拟时间
之前的
教程
里我们介绍了一个XDATCAR_toolkit.py工具，将VASP的
AIMD
的轨迹转成PDB的格式，再借用VMD和MD Analysis软件包
分析
RDF(径向分布函数)和RMSD(均方根偏差)。可惜的是，PDB文件不能承载速度信息，因此我们用c++重新写了一个VASP2GRO程序，输出可以记录原子位置和原子速度的GRO轨迹，它是 GROMACS的私有文件格式，但是也被目前的...
本系列CP2K
教程
是《CP2K菜根谭》的升级版，在旧版基础上添加了如何结合Pymatflow工具简化计算流程的内容。话不多说，本文将为您带来CP2K系列
教程
中的
第一性原理
分子
动力学
(
AIMD
)部分。简介CP2K支持
第一性原理
分子
动力学
(
AIMD
)以及混合量子经典
分子
动力学
(QM/MM)。除了波恩-奥本海默
分子
动力学
以及艾伦菲斯特
分子
动力学
以外，CP2K还实现了一种类似于传统的Car-Parrine...
前情提要无论是锂电池电极材料，还是固态电解质，锂离子的迁移速率都直接决定了其倍率性能，而倍率性能又决定了手机/汽车是1小时充满电还是1分钟冲满电。通过理论计算可以直接算出Li+的扩散系数，从而直接衡量倍率性能。昨天！锂离子电池终于拿到了诺贝尔化学奖。小编也是早就预测并编辑好了这三个人的资料，在直播公布的第一分钟就发出了新闻。2019年诺贝尔化学奖公布--锂离子电池可以预测，未来一段时间锂离子电池研...
