分子动力学后处理自编程系列（2）---角分布函数
本程序的优势主要表现在：
💖输入明确：dump文件（包含Type、mol-ID以及x，y，z）
💖简单易修改：程序思路清晰，可根据自己需求在本程序基础上进行修改实现特性需求计算，特别是针对聚合物以及蛋白质分子等复杂结构。
回转半径作为一种重要的结构表征方法，在微观研究方法中成为确定聚合物或蛋白质形态及大小的重要手段。目前对回转半径计算的后处理软件比较少，本程序可助力快速实现Rg计算。
1、编程思路
（1）读取lammps输出dump文件，将原子信息存储为变量atom_data；
（2）在帧循环条件下，依据分子ID选出属于同一分子（链单元）的原子；
（3）在链循环下，依据分子ID，计算每条链的质心以及组成该链的链单元的质心，本示例一共25条链，每个链13个链单元，共325个链单元；
（4）在链循环和链单元循环下，计算组成该链每一个链单元的质心与链质心的距离，并平均计算得出每一个链的Rg；
（5）平均计算出每一帧25条链的平均Rg，并获得偏差；
（6）完成绘图。
2、计算步骤
lammps输出命令，注意这里要输出分子id
dump 1 all custom 2000 nvt.dump id mol type x y z
(1) 读取dump文件
clc;clear;
% ------------Read Dump Process------------ %%
file = "dump.340";
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
now_frame = atom_data(:,:,frame_1);
ID = now_frame(:,1);
MOL = now_frame(:,2);
TYPE = now_frame(:,3);
XYZ = now_frame(:,4:6);
% 计算当前帧链单元和链的质心
[x_center_unit,y_center_unit,z_center_unit] = mass_center_unit(now_frame);
[x_center_chain,y_center_chain,z_center_chain] = mass_center_chain(now_frame);
for i = 1:25 % 链数
k = 0;
for j = 1:325 % 总链单元数
if ((i-1)*13<j&&j<=13*i)
k = k + 1;
square_dist_x(k,i,frame_1) = (x_center_unit(j)-x_center_chain(i))^2;
square_dist_y(k,i,frame_1) = (y_center_unit(j)-y_center_chain(i))^2;
square_dist_z(k,i,frame_1) = (z_center_unit(j)-z_center_chain(i))^2;
distance(k,i,frame_1) = sqrt(square_dist_x(k,i,frame_1)+square_dist_y(k,i,frame_1)+square_dist_z(k,i,frame_1));
end
end
end
end
计算组成链的unit的质心的函数如下：
function [x_center_unit, y_center_unit, z_center_unit] = mass_center_unit(matrix)
% NOTE: [x_center y_center z_center] = mass_center(matrix)
% This code is to calculate center of mass for unit
xnew = zeros(325,1);ynew =zeros(325,1);znew=zeros(325,1);
for j = 1:325
n_c = 0;n_f = 0;n_o = 0;n_s = 0;
for i = 1:length(matrix)
if (matrix(i,2)==j)
if (matrix(i,3)==1||matrix(i,3)==3||matrix(i,3)==6)
xnew(j,1) = xnew(j,1) + matrix(i,4)*12;
ynew(j,1) = ynew(j,1) + matrix(i,5)*12;
znew(j,1) = znew(j,1) + matrix(i,6)*12;
n_c = n_c + 1;
elseif (matrix(i,3)==2||matrix(i,3)==5||matrix(i,3)==9)
xnew(j,1) = xnew(j,1) + matrix(i,4)*19;
ynew(j,1) = ynew(j,1) + matrix(i,5)*19;
znew(j,1) = znew(j,1) + matrix(i,6)*19;
n_f = n_f + 1;
elseif (matrix(i,3)==4||matrix(i,3)==8)
xnew(j,1) = xnew(j,1) + matrix(i,4)*16;
ynew(j,1) = ynew(j,1) + matrix(i,5)*16;
znew(j,1) = znew(j,1) + matrix(i,6)*16;
n_o = n_o + 1;
elseif (matrix(i,3)==7)
xnew(j,1) = xnew(j,1) + matrix(i,4)*32;
ynew(j,1) = ynew(j,1) + matrix(i,5)*32;
znew(j,1) = znew(j,1) + matrix(i,6)*32;
n_s = n_s + 1;
end
end
end
masstotal(j) = n_c*12 + n_f*19 + n_o*16 + n_s*32;
x_center_unit(j) = xnew(j,1)/masstotal(j);
y_center_unit(j) = ynew(j,1)/masstotal(j);
z_center_unit(j) = znew(j,1)/masstotal(j);
end
计算组成链的质心的函数如下：
function [x_center_chain, y_center_chain, z_center_chain] = mass_center_chain(matrix)
% NOTE: [x_center y_center z_center] = mass_center(matrix)
% This code is to calculate center of mass for polymer chain
xnew = zeros(25,1);ynew =zeros(25,1);znew=zeros(25,1);
for j = 1:25
n_c = 0;n_f = 0;n_o = 0;n_s = 0;
for i = 1:length(matrix)
if ((j-1)*13<matrix(i,2) && matrix(i,2)<=j*13)
if (matrix(i,3)==1||matrix(i,3)==3||matrix(i,3)==6)
xnew(j,1) = xnew(j,1) + matrix(i,4)*12;
ynew(j,1) = ynew(j,1) + matrix(i,5)*12;
znew(j,1) = znew(j,1) + matrix(i,6)*12;
n_c = n_c + 1;
elseif (matrix(i,3)==2||matrix(i,3)==5||matrix(i,3)==9)
xnew(j,1) = xnew(j,1) + matrix(i,4)*19;
ynew(j,1) = ynew(j,1) + matrix(i,5)*19;
znew(j,1) = znew(j,1) + matrix(i,6)*19;
n_f = n_f + 1;
elseif (matrix(i,3)==4||matrix(i,3)==8)
xnew(j,1) = xnew(j,1) + matrix(i,4)*16;
ynew(j,1) = ynew(j,1) + matrix(i,5)*16;
znew(j,1) = znew(j,1) + matrix(i,6)*16;
n_o = n_o + 1;
elseif (matrix(i,3)==7)
xnew(j,1) = xnew(j,1) + matrix(i,4)*32;
ynew(j,1) = ynew(j,1) + matrix(i,5)*32;
znew(j,1) = znew(j,1) + matrix(i,6)*32;
n_s = n_s + 1;
end
end
end
masstotal_chain(j) = n_c*12 + n_f*19 + n_o*16 + n_s*32;
x_center_chain(j) = xnew(j,1)/masstotal_chain(j);
y_center_chain(j) = ynew(j,1)/masstotal_chain(j);
z_center_chain(j) = znew(j,1)/masstotal_chain(j);
end
(3) 统计处理
% --------------------------------------------------- %
% --------------------统计处理------------------------ %
for frame_2 = 1:all_frame
singleframe_dist_all = 0;
for i = 1:25
singlechain_dist_all = 0;
for k = 1:13
singlechain_dist_all = singlechain_dist_all + distance(k,i,frame_2);
singleframe_dist_all = singleframe_dist_all + distance(k,i,frame_2);
end
ave_singlechain = singlechain_dist_all/13;
fianl_chain_dist(i,frame_2) = ave_singlechain;
end
ave_singleframe = singleframe_dist_all/325;
fianl_frame_dist(frame_2) = ave_singleframe;
end
(4) 绘图
示例dump仅11帧，绘制出每一帧下25条链平均Rg（平均回转半径）
plot(fianl_frame_dist);
xlim([1,11]);
xlabel('Frame');
ylabel('Rg');
set(gca, 'fontsize', 16);
set(gcs,'LineWidth',20);
set(gca,'fontname','Times');
3、全部代码
%-------------------------------------------------------------------------%
clc;clear;
% ------------Read Dump Process------------ %%
file = "dump.340";
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("---------------------------------");
% --------------------------------------------------- %
% ------------计算每个链单元的质心和链的质心------------ %
all_frame = size(atom_data,3);
for frame_1 = 1:all_frame
now_frame = atom_data(:,:,frame_1);
ID = now_frame(:,1);
MOL = now_frame(:,2);
TYPE = now_frame(:,3);
XYZ = now_frame(:,4:6);
% 计算当前帧链单元和链的质心
[x_center_unit,y_center_unit,z_center_unit] = mass_center_unit(now_frame);
[x_center_chain,y_center_chain,z_center_chain] = mass_center_chain(now_frame);
for i = 1:25 % 链数
k = 0;
for j = 1:325 % 总链单元数
if ((i-1)*13<j&&j<=13*i)
k = k + 1;
square_dist_x(k,i,frame_1) = (x_center_unit(j)-x_center_chain(i))^2;
square_dist_y(k,i,frame_1) = (y_center_unit(j)-y_center_chain(i))^2;
square_dist_z(k,i,frame_1) = (z_center_unit(j)-z_center_chain(i))^2;
distance(k,i,frame_1) = sqrt(square_dist_x(k,i,frame_1)+square_dist_y(k,i,frame_1)+square_dist_z(k,i,frame_1));
end
end
end
end
% --------------------------------------------------- %
% --------------------统计处理------------------------ %
for frame_2 = 1:all_frame
singleframe_dist_all = 0;
for i = 1:25
singlechain_dist_all = 0;
for k = 1:13
singlechain_dist_all = singlechain_dist_all + distance(k,i,frame_2);
singleframe_dist_all = singleframe_dist_all + distance(k,i,frame_2);
end
ave_singlechain = singlechain_dist_all/13;
fianl_chain_dist(i,frame_2) = ave_singlechain;
end
ave_singleframe = singleframe_dist_all/325;
fianl_frame_dist(frame_2) = ave_singleframe;
end
plot(fianl_frame_dist);
xlim([1,11]);
xlabel('Frame');
ylabel('Rg');
set(gca, 'fontsize', 16);
set(gcs,'LineWidth',20);
set(gca,'fontname','Times');
3981
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
