分子动力学后处理自编程系列（3）---纳米颗粒外基液数密度分布

数密度分布同密度分布一样作为一种重要的结构描述参数，在微观研究方法中成为判断局部结构变化的重要手段。本程序利用LAMMPS输出的dump文件，目标在于计算基液原子在距离纳米颗粒表面d埃位置的数密度分布情况。

本程序的优势主要表现在：
💖输入明确：dump文件（包含Type、mol-ID以及x，y，z）、纳米颗粒的粒径R
💖计算球形颗粒外数密度：LAMMPS中可以轻易进行分区划分并统计原子数目，进而获得数密度。但是该分区为矩形，本程序针对球形颗粒外区域。
💖简单易修改：程序思路清晰，可根据自己需求在本程序基础上进行修改实现特性需求计算。

1、编程思路

（1）读取lammps输出dump文件，将原子信息存储为变量atom_data；
（2）在帧循环条件下，计算每一帧下纳米颗粒的质心位置；
（3）计算基液粒子到纳米颗粒质心的距离，结合纳米颗粒的粒径判断该基液粒子处于表面外什么位置；
（4）统计出每一帧下距离纳米颗粒表面不同位置的基液粒子个数，并对帧平均获得平均数密度
（5）完成绘图。

2、计算步骤

(1) 读取dump文件

dump的读取和该系列之前的读取代码一致，请移步之前的教程，复制该部分代码，为了避免教程冗长，在后续的教程中该部分一并省略。

(2) 输入计算相关参数

% 读取计算相关参数
all_frame = size(atom_data,3);                                               
Nlocal = Natoms(1,1);                                                        
R = 20.5;  							% 纳米颗粒半径                                                         
Square_R = R^2;
N_layers = zeros(10,all_frame);   	% 用于统计每层的原子数                                           

(3) 基液粒子位置判断并完成粒子数统计

该部分计算了纳米颗粒质心位置，之后计算了基液粒子与质心的距离并依据距离d将粒子统计到不同的空间层中。

for frame_2 = 1:all_frame
    now_frame = atom_data(:,:,frame_2);                                               
    ID     = now_frame(:,1);
    TYPE   = now_frame(:,2);
    XYZ    = now_frame(:,3:5);
    [xnew ynew znew]=mass_center(now_frame);
    n=1;
    for i = 1:Nlocal
        if(TYPE(i)==1 || TYPE(i)==2 || TYPE(i)==3)
            Square_dist(n,frame_2) = (XYZ(i,1)-xnew)^2+(XYZ(i,2)-ynew)^2+(XYZ(i,3)-znew)^2; 
            if(Square_dist(n,frame_2)>Square_R && Square_dist(n,frame_2)<= (R+1)^2)
                N_layers(1,frame_2) = N_layers(1,frame_2)+1;
            elseif(Square_dist(n,frame_2)>(R+1)^2 && Square_dist(n,frame_2)<= (R+2)^2)
                N_layers(2,frame_2) = N_layers(2,frame_2)+1;
            elseif(Square_dist(n,frame_2)>(R+2)^2 && Square_dist(n,frame_2)<= (R+3)^2)
                N_layers(3,frame_2) = N_layers(3,frame_2)+1;
            elseif(Square_dist(n,frame_2)>(R+3)^2 && Square_dist(n,frame_2)<= (R+4)^2)
                N_layers(4,frame_2) = N_layers(4,frame_2)+1;
            elseif(Square_dist(n,frame_2)>(R+4)^2 && Square_dist(n,frame_2)<= (R+5)^2)
                N_layers(5,frame_2) = N_layers(5,frame_2)+1;
            elseif(Square_dist(n,frame_2)>(R+5)^2 && Square_dist(n,frame_2)<= (R+6)^2)
                N_layers(6,frame_2) = N_layers(6,frame_2)+1;
            elseif(Square_dist(n,frame_2)>(R+6)^2 && Square_dist(n,frame_2)<= (R+7)^2)
                N_layers(7,frame_2) = N_layers(7,frame_2)+1;
            elseif(Square_dist(n,frame_2)>(R+7)^2 && Square_dist(n,frame_2)<= (R+8)^2)
                N_layers(8,frame_2) = N_layers(8,frame_2)+1;            
            elseif(Square_dist(n,frame_2)>(R+8)^2 && Square_dist(n,frame_2)<= (R+9)^2)
                N_layers(9,frame_2) = N_layers(9,frame_2)+1;            
            elseif(Square_dist(n,frame_2)>(R+9)^2 && Square_dist(n,frame_2)<= (R+10)^2)
                N_layers(10,frame_2) = N_layers(10,frame_2)+1;            
            end
            n=n+1;
        end
    end
end

mass_center函数

function [x_center y_center z_center] = mass_center(matrix)
 % NOTE: [x_center y_center z_center] = mass_center(matrix)
 % This code is to calculate center of mass for silica
 xnew = 0;ynew = 0;znew=0;silica=0;
 for i = 1:length(matrix)
     if(matrix(i,2)==4)                                 % Si
         xnew = xnew+matrix(i,3)*28;
         ynew = ynew+matrix(i,4)*28;
         znew = znew+matrix(i,5)*28;
         silica = silica + 1;
     elseif(matrix(i,2)==5)
         xnew = xnew+matrix(i,3)*16;                    % O
         ynew = ynew+matrix(i,4)*16;
         znew = znew+matrix(i,5)*16;
         silica = silica + 1;
     end
 end
 masstotal = silica*(28+16*2)/3;
 x_center = xnew/masstotal;
 y_center = ynew/masstotal;
 z_center = znew/masstotal;
end

(4) 数密度计算

for frame_3 = 1:all_frame
    for j = 1:10
        den_layers(j,frame_3)= N_layers(j,frame_3)./((4/3)*pi*((R+j)^3-(R+j-1)^3));
    end
end
den_layers_final = mean(den_layers,2); 

(5) 绘图

plot(den_layers_final);
xlim([1,10]);
ylim([0,0.032]);
xlabel('distance d');
ylabel('number density');
set(gca, 'fontsize', 16);                                              
set(gcs,'LineWidth',20);                                                
set(gca,'fontname','Times');   

3、全部代码

本程序计算的是所有TYPE的基液粒子的数密度分布情况，若有选择某类型原子的数密度分布需求的可在计算部分TYPE(i)选择处进行修改。程序仅提供思路，可在此基础上开展更复杂问题的研究，如有疑问欢迎交流。

clc;clear;
% ------------Read Dump Process------------ %%
file = "nvt.dump";
t_start_all = cputime;
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/51;      			% 这里需要手动设置一下分母，为帧数+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)
t_start_stop = (cputime-t_start_all)/60;
fprintf('Now the used time is: %.1f mins.\n',t_start_stop);
disp("---------------------------------");
disp("----Read Dump Process Done!!!----");
disp("---------------------------------");

% ------------Calculation Process------------ %%
all_frame = size(atom_data,3);                                               
Nlocal = Natoms(1,1);                                                        
R = 20.5;                                                                    Square_R = R^2;
N_layers = zeros(10,all_frame);                                              

%-------------------计算基液与纳米颗粒质心距离的平方------------------------%
for frame_2 = 1:all_frame
    now_frame = atom_data(:,:,frame_2);                                               
    ID     = now_frame(:,1);
    TYPE   = now_frame(:,2);
    XYZ    = now_frame(:,3:5);
    [xnew ynew znew]=mass_center(now_frame);
    n=1;
    for i = 1:Nlocal
        if(TYPE(i)==1 || TYPE(i)==2 || TYPE(i)==3)
            Square_dist(n,frame_2) = (XYZ(i,1)-xnew)^2+(XYZ(i,2)-ynew)^2+(XYZ(i,3)-znew)^2;
            if(Square_dist(n,frame_2)>Square_R && Square_dist(n,frame_2)<= (R+1)^2)
                N_layers(1,frame_2) = N_layers(1,frame_2)+1;
            elseif(Square_dist(n,frame_2)>(R+1)^2 && Square_dist(n,frame_2)<= (R+2)^2)
                N_layers(2,frame_2) = N_layers(2,frame_2)+1;
            elseif(Square_dist(n,frame_2)>(R+2)^2 && Square_dist(n,frame_2)<= (R+3)^2)
                N_layers(3,frame_2) = N_layers(3,frame_2)+1;
            elseif(Square_dist(n,frame_2)>(R+3)^2 && Square_dist(n,frame_2)<= (R+4)^2)
                N_layers(4,frame_2) = N_layers(4,frame_2)+1;
            elseif(Square_dist(n,frame_2)>(R+4)^2 && Square_dist(n,frame_2)<= (R+5)^2)
                N_layers(5,frame_2) = N_layers(5,frame_2)+1;
            elseif(Square_dist(n,frame_2)>(R+5)^2 && Square_dist(n,frame_2)<= (R+6)^2)
                N_layers(6,frame_2) = N_layers(6,frame_2)+1;
            elseif(Square_dist(n,frame_2)>(R+6)^2 && Square_dist(n,frame_2)<= (R+7)^2)
                N_layers(7,frame_2) = N_layers(7,frame_2)+1;
            elseif(Square_dist(n,frame_2)>(R+7)^2 && Square_dist(n,frame_2)<= (R+8)^2)
                N_layers(8,frame_2) = N_layers(8,frame_2)+1;            
            elseif(Square_dist(n,frame_2)>(R+8)^2 && Square_dist(n,frame_2)<= (R+9)^2)
                N_layers(9,frame_2) = N_layers(9,frame_2)+1;            
            elseif(Square_dist(n,frame_2)>(R+9)^2 && Square_dist(n,frame_2)<= (R+10)^2)
                N_layers(10,frame_2) = N_layers(10,frame_2)+1;            
            end
            n=n+1;
        end
    end
end

% --------------------数密度计算--------------------%
for frame_3 = 1:all_frame
    for j = 1:10
        den_layers(j,frame_3)= N_layers(j,frame_3)./((4/3)*pi*((R+j)^3-(R+j-1)^3));
    end
end
den_layers_final = mean(den_layers,2); 

% --------------------绘图--------------------%
plot(den_layers_final);
xlim([1,10]);
ylim([0,0.032]);
xlabel('distance d');
ylabel('number density');
set(gca, 'fontsize', 16);                                              
set(gcs,'LineWidth',20);                                                
set(gca,'fontname','Times');