本文介绍了一种使用MATLAB自编程计算分子动力学的键角分布函数的方法，详细阐述了编程思路和计算步骤，适用于复杂结构如聚合物和蛋白质分子。程序优势在于输入明确、计算准确且易于修改。
分子动力学后处理自编程系列（1）---角分布函数
本程序的优势主要表现在：
💖输入明确：dump文件、中心原子i的type ID、构成键角的原子j 和k的type ID以及原子i-k与原子i-j键长（若存在键为键长，不存在键则为近邻距离）
💖计算准确：在程序中加入了配位原子的判断，确保所计算键角为目标键角
💖简单易修改：程序思路清晰，可根据自己需求在本程序基础上进行修改实现特性需求计算，特别是针对聚合物以及蛋白质分子等复杂结构。
键角分布函数作为一种重要的结构表征方法，在微观研究方法中成为确定物质结构演化规律以及拓扑结构状态的重要手段。虽然现有分子动力学后处理软件及程序包大都有键角分布函数计算的功能，但ADF的计算依旧存在一些问题，主要表现在：（1）操作简单的软件功能有限，不能满足需求；（2）部分软件ADF计算非常不准确；（3）能实现特点结构计算的程序包使用起来复杂，大都需要python一顿改，对新手不友好。
因此，基于以上存在的问题，本博文主要介绍自编程ADF的MATLAB程序，该程序的优势如前所述，希望通过本程序能够解决计算问题。
1、编程思路
（1）读取dump文件，将每一帧的原子信息存储为变量atom_data；
（2）在帧循环条件下，选出每一帧下的中心原子的坐标信息；
（3）构建中心原子的近邻列表，通过选择非中心原子的类型并计算目标原子与中心原子的距离，依据判断条件(若存在键，则处于距离在键长正负0.02Å距离左右；若没有键则为小于第一溶剂层厚度)的原子归为近邻原子；
（4）由于部分中心原子只有一个近邻原子，不构成键角，因此将这些中心原子去除
（5）分别计算每一帧下每一个中心原子与近邻原子形成的键角，注意是每一个中心原子与紧邻原子的组合；
（6）记录所有帧所有中心原子所有键角的大小以及总键角数，判断键角大小分布区间(每2°一区间)并计算分布概率；
（7）完成绘图。
2、计算步骤
以SPCE水为例，计算H-O-H键角分布情况。
dump 1 all custom 2000 nvt.dump id type x y z
(1) 读取dump文件
clc;clear;
% ------------Read Dump Process------------ %%
file = "nvt.dump";
tic
try
dump = fopen(file,'r');
catch
error('Dumpfile not found!');
end
while feof(dump) == 0
id = fgetl(dump);
if (strncmpi(id,'ITEM: TIMESTEP',numel('ITEM: TIMESTEP')))
timestep(i) = str2num(fgetl(dump));
else
if (strncmpi(id,'ITEM: NUMBER OF ATOMS',numel('ITEM: NUMBER OF ATOMS')))
Natoms(i) = str2num(fgetl(dump));
else
if (strncmpi(id,'ITEM: BOX BOUNDS',numel('ITEM: BOX BOUNDS')))
x_bound(i,:) = str2num(fgetl(dump));
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
disp("---------------------------------");
disp("----Read Dump Process Done!!!----");
disp("---------------------------------");
(2) 选出中心原子并构建中心原子的配位原子近邻列表
all_frame = size(atom_data,3);
for frame_1 = 1:all_frame
xl(frame_1) = x_bound(frame_1,2)-x_bound(frame_1,1);
yl(frame_1) = y_bound(frame_1,2)-y_bound(frame_1,1);
zl(frame_1) = z_bound(frame_1,2)-z_bound(frame_1,1);
end
lengthx = mean(xl);
lengthy = mean(yl);
lengthz = mean(zl);
Nlocal = Natoms(1,1);
% -------------------- H-O-H ---------------------- %
% ------------Step 1. 选出中心原子的近邻原子----------- %
% 设置Cal_waitbar
Distance_waitbar = waitbar(0, 'Distance calculating...', 'CreateCancelBtn', 'delete(gcbf)');
set(Distance_waitbar, 'Color', [0.9, 0.9, 0.9]);
hBtn_3 = findall(Distance_waitbar, 'type', 'Uicontrol');
set(hBtn_3, 'String', 'cancle', 'FontSize', 10);
for frame_1 = 1:all_frame
compC = frame_1/all_frame;
waitbar(compC, Distance_waitbar, ['Distance Calculating...', num2str(round(compC, 2) * 100), '%']);
pause(0.01);
now_frame = atom_data(:,:,frame_1);
ID = now_frame(:,1);
TYPE = now_frame(:,2);
XYZ = now_frame(:,3:5);
n_1 = 0; % 中心原子编号
for i = 1:Nlocal
if TYPE(i) == 1 % 中心原子类型O
n_1 = n_1+1;
XYZ_center(n_1,:) = XYZ(i,:);
end
end
for j = 1:n_1 % 选择非中心原子的类型并确定每个中心原子的近邻列表
centeratomX(:,j,frame_1) = XYZ_center(j,1);
centeratomY(:,j,frame_1) = XYZ_center(j,2);
centeratomZ(:,j,frame_1) = XYZ_center(j,3);
n_2 = 0; % 中心原子周围的近邻原子编号
for k = 1:Nlocal
if (TYPE(k) == 2) % 判断近邻的原子类型2-H
distance = (centeratomX(:,j,frame_1)-XYZ(k,1))^2+(centeratomY(:,j,frame_1)-XYZ(k,2))^2+(centeratomZ(:,j,frame_1)-XYZ(k,3))^2;
if (distance>0 && distance<=1.69)
n_2 = n_2+1;
neighbor_type(n_2,j,frame_1) = TYPE(k); % 近邻原子的type
neighbor_distance(n_2,j,frame_1) = distance; % 某帧某中心原子外某近邻原子与它的距离
neighbor_X(n_2,j,frame_1) = XYZ(k,1);
neighbor_Y(n_2,j,frame_1) = XYZ(k,2);
neighbor_Z(n_2,j,frame_1) = XYZ(k,3);
end
end
end
end
end
close(Distance_waitbar)
(3) 去除近邻原子数小于2的中心原子
all_centeratom = size(neighbor_type,2);
for frame_2 = 1:all_frame
neighbor_atom_num(:,:,frame_2) = sum(neighbor_type(:,:,frame_2)~=0); %计算每列非0的个数
N_select_centeratom(frame_2) = 0;
for i = 1:all_centeratom
if (neighbor_atom_num(:,i,frame_2)>=2)
N_select_centeratom(frame_2) = N_select_centeratom(frame_2) +1;
neighbor_finaltype(:,N_select_centeratom(frame_2),frame_2) = neighbor_type(:,i,frame_2);
neighbor_finaldistance(:,N_select_centeratom(frame_2),frame_2) = neighbor_distance(:,i,frame_2);
neighbor_finalX(:,N_select_centeratom(frame_2),frame_2) = neighbor_X(:,i,frame_2);
neighbor_finalY(:,N_select_centeratom(frame_2),frame_2) = neighbor_Y(:,i,frame_2);
neighbor_finalZ(:,N_select_centeratom(frame_2),frame_2) = neighbor_Z(:,i,frame_2);
center_fianlX(:,N_select_centeratom(frame_2),frame_2) = centeratomX(:,i,frame_2);
center_fianlY(:,N_select_centeratom(frame_2),frame_2) = centeratomY(:,i,frame_2);
center_fianlZ(:,N_select_centeratom(frame_2),frame_2) = centeratomZ(:,i,frame_2);
end
end
end
(4) 计算近邻中任意两个原子与中心原子的夹角
number_angle = 0;
for frame_3 = 1:all_frame
neighbor_atom_fianlnum(:,:,frame_3) = sum(neighbor_finaltype(:,:,frame_3)~=0);
for i = 1:N_select_centeratom(frame_3)
N = neighbor_atom_fianlnum(:,i,frame_3);
fid = 1:1:N;
fid_choose = nchoosek(fid,2);
X_center = center_fianlX(:,i,frame_3);
Y_center = center_fianlY(:,i,frame_3);
Z_center = center_fianlZ(:,i,frame_3);
for j = 1:size(fid_choose,1)
number_angle = number_angle + 1;
X_1 = neighbor_finalX(fid_choose(j,1),i,frame_3);
Y_1 = neighbor_finalY(fid_choose(j,1),i,frame_3);
Z_1 = neighbor_finalZ(fid_choose(j,1),i,frame_3);
X_2 = neighbor_finalX(fid_choose(j,2),i,frame_3);
Y_2 = neighbor_finalY(fid_choose(j,2),i,frame_3);
Z_2 = neighbor_finalZ(fid_choose(j,2),i,frame_3);
dis_c1 = sqrtm((X_center-X_1)^2+(Y_center-Y_1)^2+(Z_center-Z_1)^2);
dis_c2 = sqrtm((X_center-X_2)^2+(Y_center-Y_2)^2+(Z_center-Z_2)^2);
dis_12 = sqrtm((X_2-X_1)^2+(Y_2-Y_1)^2+(Z_2-Z_1)^2);
costheta(number_angle) = (dis_c1^2 + dis_c2^2 - dis_12^2)/(2*dis_c1*dis_c2);
theta(number_angle) = acos(costheta(number_angle))*(180/pi); % 计算键角的值
end
end
end
(5) 角分布统计及绘图
% 分区，每2°一个区间，共36个区间（0~180）
angle_interval = 2;
N_zones = 180/angle_interval;
angle_distribution = zeros(N_zones+1,3);
for i = 1: length(theta)
for j = 1:N_zones
angle(j) = j*angle_interval;
angle_distribution(j+1,1) = angle(j);
if j==1
if theta(i)<angle(j)
angle_distribution(j+1,2) = angle_distribution(j+1,2)+1;
end
elseif (theta(i)>angle(j-1)&&theta(i)<angle(j))
angle_distribution(j+1,2) = angle_distribution(j+1,2)+1;
end
end
end
for i = 1:N_zones
angle_distribution(i,3) = angle_distribution(i,2)/length(theta);
end
plot(angle_distribution(:,1),angle_distribution(:,3));
xlim([0,180]);
xlabel('angle(°)');
ylabel('Frequency');
set(gca, 'fontsize', 16);
set(gcs,'LineWidth',20);
set(gca,'fontname','Times');
3、全部代码
clc;clear;
% ------------Read Dump Process------------ %%
file = "nvt.dump";
tic
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');
set(hBtn_3, 'String', 'cancle', 'FontSize', 10);
i=1;
while feof(dump) == 0
compT = i/6; % 这里需要手动设置一下分母，为帧数+1
waitbar(compT, ReadDump_Waitbar, ['Dump info reading...', num2str(round(compT, 2) * 100), '%']);
pause(0.01);
id = fgetl(dump);
if (strncmpi(id,'ITEM: TIMESTEP',numel('ITEM: TIMESTEP')))
timestep(i) = str2num(fgetl(dump));
else
if (strncmpi(id,'ITEM: NUMBER OF ATOMS',numel('ITEM: NUMBER OF ATOMS')))
Natoms(i) = str2num(fgetl(dump));
else
if (strncmpi(id,'ITEM: BOX BOUNDS',numel('ITEM: BOX BOUNDS')))
x_bound(i,:) = str2num(fgetl(dump));
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)
disp("---------------------------------");
disp("----Read Dump Process Done!!!----");
disp("---------------------------------");
% ------------Calculation Process------------ %%
all_frame = size(atom_data,3);
%xl,yl,zl
for frame_1 = 1:all_frame
xl(frame_1) = x_bound(frame_1,2)-x_bound(frame_1,1);
yl(frame_1) = y_bound(frame_1,2)-y_bound(frame_1,1);
zl(frame_1) = z_bound(frame_1,2)-z_bound(frame_1,1);
end
lengthx = mean(xl);
lengthy = mean(yl);
lengthz = mean(zl);
Nlocal = Natoms(1,1);
% ------------------ H-O-H -------------------- %
% ------------Step 1. 选出中心原子的近邻原子------------ %
% 设置Cal_waitbar
Distance_waitbar = waitbar(0, 'Distance calculating...', 'CreateCancelBtn', 'delete(gcbf)');
set(Distance_waitbar, 'Color', [0.9, 0.9, 0.9]);
hBtn_3 = findall(Distance_waitbar, 'type', 'Uicontrol');
set(hBtn_3, 'String', 'cancle', 'FontSize', 10);
for frame_1 = 1:all_frame
compC = frame_1/all_frame;
waitbar(compC, Distance_waitbar, ['Distance Calculating...', num2str(round(compC, 2) * 100), '%']);
pause(0.01);
now_frame = atom_data(:,:,frame_1);
ID = now_frame(:,1);
TYPE = now_frame(:,2);
XYZ = now_frame(:,3:5);
n_1 = 0;
for i = 1:Nlocal
if TYPE(i) == 1
n_1 = n_1+1;
XYZ_center(n_1,:) = XYZ(i,:);
end
end
for j = 1:n_1
centeratomX(:,j,frame_1) = XYZ_center(j,1);
centeratomY(:,j,frame_1) = XYZ_center(j,2);
centeratomZ(:,j,frame_1) = XYZ_center(j,3);
n_2 = 0;
for k = 1:Nlocal
if (TYPE(k) == 2)
distance = (centeratomX(:,j,frame_1)-XYZ(k,1))^2+(centeratomY(:,j,frame_1)-XYZ(k,2))^2+(centeratomZ(:,j,frame_1)-XYZ(k,3))^2;
if (distance>0 && distance<=1.69)
n_2 = n_2+1;
neighbor_type(n_2,j,frame_1) = TYPE(k);
neighbor_distance(n_2,j,frame_1) = distance;
neighbor_X(n_2,j,frame_1) = XYZ(k,1);
neighbor_Y(n_2,j,frame_1) = XYZ(k,2);
neighbor_Z(n_2,j,frame_1) = XYZ(k,3);
end
end
end
end
end
close(Distance_waitbar)
% 选出近邻原子大于等于2个的中心原子
all_centeratom = size(neighbor_type,2);
for frame_2 = 1:all_frame
neighbor_atom_num(:,:,frame_2) = sum(neighbor_type(:,:,frame_2)~=0);
N_select_centeratom(frame_2) = 0;
for i = 1:all_centeratom
if (neighbor_atom_num(:,i,frame_2)>=2)
N_select_centeratom(frame_2) = N_select_centeratom(frame_2) +1;
neighbor_finaltype(:,N_select_centeratom(frame_2),frame_2) = neighbor_type(:,i,frame_2);
neighbor_finaldistance(:,N_select_centeratom(frame_2),frame_2) = neighbor_distance(:,i,frame_2);
neighbor_finalX(:,N_select_centeratom(frame_2),frame_2) = neighbor_X(:,i,frame_2);
neighbor_finalY(:,N_select_centeratom(frame_2),frame_2) = neighbor_Y(:,i,frame_2);
neighbor_finalZ(:,N_select_centeratom(frame_2),frame_2) = neighbor_Z(:,i,frame_2);
center_fianlX(:,N_select_centeratom(frame_2),frame_2) = centeratomX(:,i,frame_2);
center_fianlY(:,N_select_centeratom(frame_2),frame_2) = centeratomY(:,i,frame_2);
center_fianlZ(:,N_select_centeratom(frame_2),frame_2) = centeratomZ(:,i,frame_2);
end
end
end
% ------------Step 2. 计算任意两个原子与中心原子的夹角------------ %
number_angle = 0;
for frame_3 = 1:all_frame
neighbor_atom_fianlnum(:,:,frame_3) = sum(neighbor_finaltype(:,:,frame_3)~=0);
for i = 1:N_select_centeratom(frame_3)
N = neighbor_atom_fianlnum(:,i,frame_3);
fid = 1:1:N;
fid_choose = nchoosek(fid,2);
X_center = center_fianlX(:,i,frame_3);
Y_center = center_fianlY(:,i,frame_3);
Z_center = center_fianlZ(:,i,frame_3);
for j = 1:size(fid_choose,1)
number_angle = number_angle + 1;
X_1 = neighbor_finalX(fid_choose(j,1),i,frame_3);
Y_1 = neighbor_finalY(fid_choose(j,1),i,frame_3);
Z_1 = neighbor_finalZ(fid_choose(j,1),i,frame_3);
X_2 = neighbor_finalX(fid_choose(j,2),i,frame_3);
Y_2 = neighbor_finalY(fid_choose(j,2),i,frame_3);
Z_2 = neighbor_finalZ(fid_choose(j,2),i,frame_3);
dis_c1 = sqrtm((X_center-X_1)^2+(Y_center-Y_1)^2+(Z_center-Z_1)^2);
dis_c2 = sqrtm((X_center-X_2)^2+(Y_center-Y_2)^2+(Z_center-Z_2)^2);
dis_12 = sqrtm((X_2-X_1)^2+(Y_2-Y_1)^2+(Z_2-Z_1)^2);
costheta(number_angle) = (dis_c1^2 + dis_c2^2 - dis_12^2)/(2*dis_c1*dis_c2);
theta(number_angle) = acos(costheta(number_angle))*(180/pi);
end
end
end
% 分区，每2°一个区间，共36个区间（0~180）
angle_interval = 2;
N_zones = 180/angle_interval;
angle_distribution = zeros(N_zones+1,3);
for i = 1: length(theta)
for j = 1:N_zones
angle(j) = j*angle_interval;
angle_distribution(j+1,1) = angle(j);
if j==1
if theta(i)<angle(j)
angle_distribution(j+1,2) = angle_distribution(j+1,2)+1;
end
elseif (theta(i)>angle(j-1)&&theta(i)<angle(j))
angle_distribution(j+1,2) = angle_distribution(j+1,2)+1;
end
end
end
for i = 1:N_zones
angle_distribution(i,3) = angle_distribution(i,2)/length(theta);
end
plot(angle_distribution(:,1),angle_distribution(:,3));
xlim([0,180]);
xlabel('angle(°)');
ylabel('Frequency');
set(gca, 'fontsize', 16);
set(gcs,'LineWidth',20);
set(gca,'fontname','Times');
toc
disp(['运行时间: ',num2str(toc)]);
1973
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
