- （1）lammps计算
- （2）VMD计算
- （3）OVITO计算
- （4）ISAACS软件计算
- （5）自编程计算-附MATLAB代码
- 总结
在统计力学中，均方位移（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]
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:\Users\Lenlovo\Desktop\MSD_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:\Users\Lenlovo\Desktop\msd_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下的帧平均（帧平均）。
此外，由于计算过程的周期性边界的作用，任何一帧下原子坐标范围均在box范围内，为了反映原子真实的扩散位置情况，在计算过程中需要利用前一帧对后一帧进行位置修正，并用修正后的位置信息去计算MSD。
最终，计算结果与LAMMPS对比如下，可见计算结果与lammps计算结果基本一致，且相较于lammps计算结果线性度更好。
MSD计算matlab代码，代码不够简练，但思路较为清晰哈哈，计算速度可以接受。
%--------------------------------------------%
%----------Author:磊子的MD记录簿---------------%
%-------Email:Victor_Xian285@163.com---------%
%--------------------------------------------%
clc;clear;
t_start_all = cputime;
%----------Step 1-------------
% Step 1: read dump
file ="NaCl_DuiBi.xyz";
try
dump = fopen(file,'r');
catch
error('Dumpfile not found!');
end
% 设置ReadDump_waitbar
ReadDump_Waitbar = waitbar(0, 'Dump info reading...', 'CreateCancelBtn', 'delete(gcbf)');% 进度条初始设置
set(ReadDump_Waitbar, 'Color', [0.9, 0.9, 0.9]);
hBtn_3 = findall(ReadDump_Waitbar, 'type', 'Uicontrol'); % cancel按钮设置
set(hBtn_3, 'String', 'cancle', 'FontSize', 10);
i=1;
while feof(dump) == 0 % read to end of file
compT = i/202; % 这里需要手动设置一下分母，为帧数+1
waitbar(compT, ReadDump_Waitbar, ['Dump info reading...', num2str(round(compT, 2) * 100), '%']);
pause(0.01);
id = fgetl(dump); % back the next line content
if (strncmpi(id,'ITEM: TIMESTEP',numel('ITEM: TIMESTEP')))
timestep(i) = str2num(fgetl(dump)); % current timestep
else
if (strncmpi(id,'ITEM: NUMBER OF ATOMS',numel('ITEM: NUMBER OF ATOMS')))
Natoms(i) = str2num(fgetl(dump)); % number of atoms
else
if (strncmpi(id,'ITEM: BOX BOUNDS',numel('ITEM: BOX BOUNDS')))
x_bound(i,:) = str2num(fgetl(dump)); % n*2
y_bound(i,:) = str2num(fgetl(dump));
z_bound(i,:) = str2num(fgetl(dump));
else
if (strcmpi(id(1:11),'ITEM: ATOMS'))
for j = 1 : 1: Natoms
atom_data(j,:,i) = str2num(fgetl(dump));
end
i=i+1;
end
end
end
end
end
close(ReadDump_Waitbar)
frame_num = i-1;
num_atoms = j;
xl = x_bound(1,2)-x_bound(1,1); % x length of box
yl = y_bound(1,2)-y_bound(1,1); % y length of box
zl = z_bound(1,2)-z_bound(1,1); % z length of box
% ----------Step 2-------------
% 位移修正，记录为real_data
% 设置DispCor_waitbar
DispCor_Waitbar = waitbar(0, 'Displacement correction...', 'CreateCancelBtn', 'delete(gcbf)');% 进度条初始设置
set(DispCor_Waitbar, 'Color', [0.9, 0.9, 0.9]);
hBtn_1 = findall(DispCor_Waitbar, 'type', 'Uicontrol'); % cancel按钮设置
set(hBtn_1, 'String', 'cancle', 'FontSize', 10);
% 位移修正
real_data = atom_data;
for i1 = 1:frame_num
real_data(:,:,i1) = sortrows(real_data(:,:,i1),1); % 先调整一下real_data的排序
end
for k = 1: frame_num-1
compT_1 = k/(frame_num-1);
waitbar(compT_1, DispCor_Waitbar, ['Displacement correction...', num2str(round(compT_1, 2) * 100), '%']);
pause(0.01);
current_data = atom_data(:,1:5,k+1);
current_data = sortrows(current_data,1);
current_x = current_data(:,3);
current_y = current_data(:,4);
current_z = current_data(:,5);
last_data = atom_data(:,1:5,k);
last_data = sortrows(last_data,1);
last_x = last_data(:,3);
last_y = last_data(:,4);
last_z = last_data(:,5);
for n = 1:num_atoms
d_x(n) = current_x(n) - last_x(n);
d_y(n) = current_y(n) - last_y(n);
d_z(n) = current_z(n) - last_z(n);
[d_x(n),d_y(n),d_z(n)] = PBC(d_x(n),d_y(n),d_z(n),xl,yl,zl);
real_data(n,3,k+1) = real_data(n,3,k)+ d_x(n); % 修正x
real_data(n,4,k+1) = real_data(n,4,k)+ d_y(n); % 修正y
real_data(n,5,k+1) = real_data(n,5,k)+ d_z(n); % 修正z
end
end
close(DispCor_Waitbar)
% ----------Step 3-------------
% MSD计算
mean_squard_dx_2 = zeros(frame_num-1,1);
mean_squard_dy_2 = zeros(frame_num-1,1);
mean_squard_dz_2 = zeros(frame_num-1,1);
mean_squard_all_2 = zeros(frame_num,1);
% 设置MSD_waitbar
MSD_Waitbar = waitbar(0, 'MSD computing...', 'CreateCancelBtn', 'delete(gcbf)');% 进度条初始设置
set(MSD_Waitbar, 'Color', [0.9, 0.9, 0.9]);
hBtn_2 = findall(MSD_Waitbar, 'type', 'Uicontrol'); % cancel按钮设置
set(hBtn_2, 'String', 'cancle', 'FontSize', 10);
% 循环不同的时间间隔，最大为最后一帧减第一帧为frame_num-1
for tau = 1:frame_num-1
compT_2 = tau/(frame_num-1);
waitbar(compT_2, MSD_Waitbar, ['Different tau value Computing...', num2str(round(compT_2, 2) * 100), '%']);
pause(0.01);
% 循环计算公式中的不同的分子项数，与时间间隔tau有关，最大等于frame_num-tau
for frame = 1: frame_num-tau
now_frame = real_data(:,1:5,frame + tau);
% now_frame = sortrows(now_frame,1);
now_x = now_frame(:,3);
now_y = now_frame(:,4);
now_z = now_frame(:,5);
before_frame = real_data(:,1:5,frame);
% before_frame = sortrows(before_frame,1);
before_x = before_frame(:,3);
before_y = before_frame(:,4);
before_z = before_frame(:,5);
for n = 1:num_atoms
dx(n) = now_x(n) - before_x(n);
dy(n) = now_y(n) - before_y(n);
dz(n) = now_z(n) - before_z(n);
squard_dx(n,frame,tau) = dx(n).^2;
squard_dy(n,frame,tau) = dy(n).^2;
squard_dz(n,frame,tau) = dz(n).^2;
end
% 第一次平均是对原子平均
mean_squard_dx_1(frame,tau) = mean(squard_dx(:,frame,tau));
mean_squard_dy_1(frame,tau) = mean(squard_dy(:,frame,tau));
mean_squard_dz_1(frame,tau) = mean(squard_dz(:,frame,tau));
end
% 第二次是对分子中加项数平均
mean_squard_dx_2(tau) = sum(mean_squard_dx_1(:,tau))./(frame_num - tau);
mean_squard_dy_2(tau) = sum(mean_squard_dy_1(:,tau))./(frame_num - tau);
mean_squard_dz_2(tau) = sum(mean_squard_dz_1(:,tau))./(frame_num - tau);
mean_squard_all_2(tau+1) = mean_squard_dx_2(tau) + mean_squard_dy_2(tau) + mean_squard_dz_2(tau);
end
t_start_stop = (cputime-t_start_all)/60;
fprintf('Now the used time is: %.1f mins.\n',t_start_stop);
disp("-------------------");
disp("----ALL DONE!!!----");
disp("-------------------");
在分子动力学学习中，虽然现在有很多软件以及分析工具包可以帮助我们快速的计算某些参数，但后处理工具越来越多也带来使用复杂（例如isaacs就必须用xyz标准格式，还要python处理一下），计算可靠度难衡量的问题。因此，学会自编程开展计算或结构判断是形成独特且有深度研究成果的必要技能，虽然MSD计算公式看似很简单，但在自编程计算中需要进行多个细节处理（也可能我太菜了。。。），也学习编程处理MD原子信息的方法，收获良多。
这严格意义来说编出了我第一个MD计算程序，特此记录一下，继续踏实学习。
欢迎分享，转载请注明来源：内存溢出
微信扫一扫
支付宝扫一扫
评论列表（0条）