(1)lammps计算(2)VMD计算(3)OVITO计算(4)ISAACS软件计算(5)自编程计算总结
在统计力学中,均方位移(MSD,均方位移或均方波动)是粒子随时间移动后的位置相对于参考位置的偏差的量度。它是随机运动中空间范围的最常见度量,可以被认为是对随机行者“探索”的系统部分进行度量。在生物物理学和环境工程领域,均方根位移是随时间测量的,以确定粒子是否仅由于扩散而扩散,或者对流力是否也在起作用。
统
计
力
学
中
,
M
S
D
定
义
为
t
时
刻
的
系
综
平
均
:
color{red}统计力学中,MSD定义为t时刻的系综平均:
统计力学中,MSD定义为t时刻的系综平均:
下面分享几种常用的MSD计算方法,便于大家选择使用。
(1)lammps计算
lammps程序包通过compute MSD命令计算MSD,计算结果有四个,保存在C_msd[1~4]数组中,前三个向量分别为x,y,z方向的MSD,第四个是总均方位移。
如果 com 选项设置为 yes,则在计算每个原子的位移之前,将减去原子组质心中任何漂移的影响。
fix 1 all nvt temp 300 300 1 compute msd all msd com yes fix 2 all vector 10 c_msd[4] fix 3 all ave/time 10 1000 10000 c_msd[*] file MSD.out mode vector thermo_style custom step temp c_msd[4](2)VMD计算
VMD中可进行RMSD的计算,RMSD的计算公式如下:
VMD的使用:文件导入→Extensions→Analysis→RMSD Trajectory Tool
网页教学网址 link.
视频教学网址 link.
OVITO这个功能贼傻*,一直报错,xyz文件好像必须用标准格式,最好跟案例一样用dump文件,懒得搞了。
而且如果计算某类原子的MSD需要在次代码上修改,对于不懂OVITO python算法的上手有些困难,总之放弃这个选择!
哪位仁兄用这个比较上手求指教。
from ovito.io import import_file, export_file
from ovito.modifiers import CalculateDisplacementsModifier
import numpy
# Load input data and create a data pipeline.
pipeline = import_file(r"C:UsersLenlovoDesktopMSD_final.xyz")
# Calculate per-particle displacements with respect to initial simulation frame:
pipeline.modifiers.append(CalculateDisplacementsModifier())
# Define the custom modifier function:
def calculate_msd(frame, data):
# Access the per-particle displacement magnitudes computed by the
# CalculateDisplacementsModifier that precedes this user-defined modifier in the
# data pipeline:
displacement_magnitudes = data.particles['Displacement Magnitude']
# Compute MSD:
msd = numpy.sum(displacement_magnitudes ** 2) / len(displacement_magnitudes)
# Output MSD value as a global attribute:
data.attributes["MSD"] = msd
# Insert user-defined modifier function into the data pipeline.
pipeline.modifiers.append(calculate_msd)
# Export calculated MSD value to a text file and let OVITO's data pipeline do the rest:
export_file(pipeline, r"C:UsersLenlovoDesktopmsd_data.txt",
format = "txt/attr",
columns = ["Timestep", "MSD"],
multiple_frames = True)
OVITO官方MSD计算教学链接 link.
(4)ISAACS软件计算ISAACS软阿介绍及安装可以看我上一篇博客,需将lammps导出的XYZ文件转换为标准格式才可以使用,写个python脚本处理一下即可。
以NaCl为例,计算了Na以及Cl的MSD,如下图
与lammps计算结果如下图,???好像有点不太对,lammps和ISAACS软件计算的Na的MSD对不上,搞错了再来。
之后又用了Ar进行测试,也是不行,与lammps对不上差距较大,不能理解。
lammps导出201帧的dump文件,使用MATLAB读取位置信息,MSD的平均使用的是all pairs方法,因为这种平均使用了所有的数据,而且它对数据点的权重几乎相同。计算过程通过两次平均,第一次是对该帧下所有原子位移量进行平均(原子平均),第二次是对不同lag下的帧平均(帧平均)。
此外,在计算过程中需要利用前一帧对后一帧进行位置修正。
最终,计算结果与LAMMPS对比如下,可见计算结果与lammps计算结果基本一致。
在分子动力学学习中,虽然现在有很多软件以及分析工具包可以帮助我们快速的计算某些参数,但后处理工具越来越多也带来使用复杂(例如isaacs就必须用xyz标准格式,还要python处理一下),计算可靠度难衡量的问题。因此,学会自编程开展计算或结构判断是形成独特且有深度研究成果的必要技能,虽然MSD计算公式看似很简单,但在自编程计算中需要进行多个细节处理(也可能我太菜了。。。),也学习编程处理MD原子信息的方法,收获良多。
这严格意义来说编出了我第一个MD计算程序,特此记录一下,继续踏实学习。



