分子动力学基本分析-数密度二维云图及一维分布计算

Mr. Material

已于 2024-07-16 13:25:55 修改

阅读量1k

点赞数

分类专栏：分子动力学—LAMMPS后处理及编程技巧文章标签： Lammps 数密度计算

于 2024-07-14 22:14:54 首次发布

本文链接：https://blog.csdn.net/qq_43689832/article/details/140424299

版权

分子动力学—LAMMPS后处理及编程技巧专栏收录该内容

57 篇文章 ¥99.90 ¥299.90

订阅专栏

超级会员免费看

关注 $\color{Violet} \rm Mr.material\ ,$ $\color{red}{更}$ $\color{blue}{多}$ $\color{orange}{精}$ $\color{green}{彩}$ ！

主要专栏内容包括：
$\textbf{ \underline{\dag《LAMMPS小技巧》：}}$ 主要介绍采用分子动力学（ $L amm p s$ ）模拟相关安装教程、原理以及模拟小技巧（难度： $\bigstar$ ）
$\textbf{ \underline{\dag\dag《LAMMPS实例教程—In文件详解》：}}$ 主要介绍采用分子动力学（ $L amm p s$ ）模拟相关物理过程模拟。（包含：热导率计算、定压比热容计算，难度： $\bigstar$ $\bigstar$ $\bigstar$ ）
$\textbf{ \underline{\dag\dag\dag《Lammps编程技巧及后处理程序技巧》：}}$ 主要介绍针对分子模拟的动力学过程（轨迹文件）进行后相关的处理分析（需要一定编程能力。难度： $\bigstar$ $\bigstar$ $\bigstar$ $\bigstar$ $\bigstar$ ）。
$\textbf{ \underline{\dag\dag\dag\dag《分子动力学后处理集成函数—Matlab》：}}$ 主要介绍针对后处理过程中指定函数，进行包装，方便使用者直接调用（需要一定编程能力，难度： $\bigstar$ $\bigstar$ $\bigstar$ $\bigstar$ ）。
$\textbf{ \underline{\dag\dag\dag\dag\dag《SCI论文绘图—Python绘图常用模板及技巧》：}}$ 主要介绍针对处理后的数据可视化，并提供对应的绘图模板（需要一定编程能力，难度： $\bigstar$ $\bigstar$ $\bigstar$ $\bigstar$ ）。
$\textbf{ \underline{\dag\dag\dag\dag\dag\dag《分子模拟—Ovito渲染案例教程》：}}$ 主要采用 $\rm Ovito$ 软件，对 $\rm Lammps$ 生成的轨迹文件进行渲染（难度： $\bigstar$ $\bigstar$ ）。

$\color{red}{\textbf{ \underline{专栏说明（订阅后可浏览对应专栏全部博文）：}}}$
$\color{red} 注意：$ 如需只订阅某个单独博文，请联系博主邮箱咨询。 $\rm lammps\_materials@163.com$

$\spadesuit$ $\dag$ 开源后处理集成程序：请关注专栏《LAMMPS后处理——MATLAB子函数合集整理》
$\spadesuit$ $\dag$ $\dag$ 需要付费定制后处理程序请邮件联系： $\rm lammps\_materials@163.com$

请添加图片描述

一、数密度计算思路

数密度可以定义为：原子的数量/所在空间的体积，单位为 $\rm \AA^{-3}$ 。

根据上面的分析可以看出，计算流程为：

1）对体系进行划分；

针对一维数密度分布的计算，需要沿着该方向进行划分，统计原子数量，冰并除以对应的空间体积。

针对二维数密度分布的计算，需要先将三维空间划分为二维空间，统计每一个格子内原子的数量，并除以体积。

2）计算统计划分后的空间内原子数量；

二、读取dump轨迹文件获得坐标

$\rm \color{red} 注意：请按照：ID\ TYPE\ X\ Y\ Z\ 的顺序在lammps中输出dump$

clc;clear;
% 这里写出需要读取的dump文件及路径
file ="test.dump";
%========================================%
%  这里开始读dump文件，这里不用修改程序！ %
%========================================%

t_start_all = cputime;
try
    dump = fopen(file,'r');
catch
    error('Dumpfile not found!');
end

i=1;
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
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("-------------------");
%========================================%

三、数密度一维分布统计

将原子按三个方向划分，依次统计每个空间内的原子数量。

%------------------------------------------------%          
            if((XYZ(i,1)>(block-1)*(xl/bin_num)) && (XYZ(i,1)<=block*(xl/bin_num)))
                  
                DEN_bins_x(block,frame)=((DEN_bins_x(block,frame)+1));  

            end     
% XYZ(i,2) y方向           
%------------------------------------------------%            
            if((XYZ(i,2)>(block-1)*(yl/bin_num)) && (XYZ(i,2)<=block*(yl/bin_num)))
                  
                DEN_bins_y(block,frame)=((DEN_bins_y(block,frame)+1));  

            end  
% XYZ(i,3) z方向            
%------------------------------------------------%            
            if((XYZ(i,3)>(block-1)*(zl/bin_num)) && (XYZ(i,3)<=block*(zl/bin_num)))
                  
                DEN_bins_z(block,frame)=((DEN_bins_z(block,frame)+1));  

            end

四、数密度二维分布统计

将体系划分为二维空间，根据原子距离模拟体系边界的位置统计个数

                if(1)   
                    
                    x_n = floor((XYZ(i,1)-x_bound(1,1))/bin_ALL_mesh_xsize)+1;
                    y_n = floor((XYZ(i,2)-x_bound(1,1))/bin_ALL_mesh_ysize)+1;                    
                    z_n = floor((XYZ(i,3)-z_bound(1,1))/bin_ALL_mesh_zsize)+1;          
                    
                    if(x_n<=0)
                       x_n=1;
                    elseif(x_n>=bin_ALL_mesh)
                        x_n=bin_ALL_mesh;
                    end
                    
                    if(y_n<=0)
                       y_n=1;
                    elseif(y_n>=bin_ALL_mesh)
                        y_n=bin_ALL_mesh;
                    end          
                    
                    if(z_n<=0)
                       z_n=1;
                    elseif(z_n>=bin_ALL_mesh)
                       z_n=bin_ALL_mesh;   
                    end
                    
                     DEN_mesh_xy(x_n,y_n)=(DEN_mesh_xy(x_n,y_n)+1);
                     DEN_mesh_yz(y_n,z_n)=(DEN_mesh_yz(y_n,z_n)+1);
                     DEN_mesh_xz(x_n,z_n)=(DEN_mesh_xz(x_n,z_n)+1);
                end

五、MATLAB绘图

  
    drawnow;
    
    x_real_length = ((1:length(DEN_bins_x)).*xl./bin_num)';
    y_real_length = ((1:length(DEN_bins_y)).*yl./bin_num)';
    z_real_length = ((1:length(DEN_bins_z)).*zl./bin_num)';
    %------------------------------------------------%  
    set(gcf, 'Position', [300 200 1500 800]); % gcf 获取当前figure的句柄

    subplot(2,3,1);  box on;hold on;
    title(strcat("frame: ",num2str(frame)));
    plot(x_real_length,(DEN_bins_x(:,frame)./bin_size./yl./zl),"ko-");
    xlim([0,max(x_real_length)]);
    xlabel('Length /Å','fontsize',20);
    ylabel('Density /Å^-3','fontsize',20);
    set(gca,'fontsize',20,'linewidth',1.5);
    legend ('Along X-axis','Location','southEast');        
    %------------------------------------------------%  
    set(gcf, 'Position', [300 200 1500 800]); % gcf 获取当前figure的句柄    
    subplot(2,3,2);  box on;hold on;
    plot(y_real_length,DEN_bins_y(:,frame)./bin_size./xl./zl,"r^-");
    xlim([0,max(x_real_length)]);
    xlabel('Length /Å','fontsize',20);
    ylabel('Density /Å^-3','fontsize',20);
    set(gca,'fontsize',20,'linewidth',1.5);
    legend ('Along Y-axis','Location','southEast');        
    %------------------------------------------------%  
    subplot(2,3,3);  box on;hold on;
    plot(z_real_length,DEN_bins_z(:,frame)./bin_size./xl./yl,"b*-");
    xlim([0,max(x_real_length)]);
    xlabel('Length /Å','fontsize',20);
    ylabel('Density /Å^-3','fontsize',20);
    set(gca,'fontsize',20,'linewidth',1.5);
    legend ('Along Z-axis','Location','southEast');
    %------------------------------------------------%  
    
    subplot(2,3,4);   box on;hold on;
    surf(DEN_mesh_xz./bin_ALL_mesh_xsize./bin_ALL_mesh_zsize./yl);
    shading interp;
    xlabel('Bin number','fontsize',20);
    ylabel('Bin number','fontsize',20);
    zlabel('Bin number','fontsize',20);
    
    
    set(gca,'fontsize',20,'linewidth',1.5);
    legend ('X-Z','Location','southEast');
    %------------------------------------------------%      
    subplot(2,3,5);  box on;hold on;
    surf(DEN_mesh_xy./bin_ALL_mesh_xsize./bin_ALL_mesh_ysize./zl);
    shading interp;
    
    xlabel('Bin number','fontsize',20);
    ylabel('Bin number','fontsize',20);
    zlabel('Bin number','fontsize',20);
    
    set(gca,'fontsize',20,'linewidth',1.5);
    legend ('X-Y','Location','southEast');
    %------------------------------------------------%  
    
    subplot(2,3,6);  box on;hold on;
    surf(DEN_mesh_yz./bin_ALL_mesh_ysize./bin_ALL_mesh_zsize./xl);
    shading interp;    
    
    
    xlabel('Bin number','fontsize',20);
    ylabel('Bin number','fontsize',20);
    zlabel('Bin number','fontsize',20);
    
    
    set(gca,'fontsize',20,'linewidth',1.5);
    legend ('Y-Z','Location','southEast');
    %------------------------------------------------%

六、全部代码

clc;clear;
% 这里写出需要读取的dump文件及路径
file ="test.dump";
%========================================%
%  这里开始读dump文件，这里不用修改程序！ %
%========================================%

t_start_all = cputime;
try
    dump = fopen(file,'r');
catch
    error('Dumpfile not found!');
end

i=1;
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
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("-------------------");
%========================================%

% 请修改这里！
bin_num   = 100;        % 将体系划分block,计算沿着某一方向的密度分布（一维）

%%----------------------------------------------%%

bin_ALL_mesh = 100;   % 将体系划分block,计算沿着某一方向的密度分布（二维）
%%----------------------------------------------%%

%（不用修改）
all_frame = size(atom_data,3);  % 全部需要计算的帧
%------------------------------------------------% 
DEN_mesh_xz = zeros(bin_ALL_mesh,bin_ALL_mesh);%
DEN_mesh_xy= zeros(bin_ALL_mesh,bin_ALL_mesh);%
DEN_mesh_yz = zeros(bin_ALL_mesh,bin_ALL_mesh);%


DEN_bins_x = zeros(bin_num,all_frame);
DEN_bins_y = zeros(bin_num,all_frame);
DEN_bins_z = zeros(bin_num,all_frame);
Block_Size_bins = zeros(all_frame,1);

%%----------------------------------------------%%
number_check = zeros(bin_num,all_frame);
%%----------------------------------------------%%

Nlocal = Natoms(1,1);
t_all_start=cputime;

for frame = 1:all_frame

    
%计算盒子大小（不用修改）
%------------------------------------------------% 
xl = x_bound(1,2)-x_bound(1,1);
yl = y_bound(1,2)-y_bound(1,1);
zl = z_bound(1,2)-z_bound(1,1);
%------------------------------------------------% 

%------------------------------------------------% 
    t_start = cputime;
    if(rem(frame,50)==0)
        disp('%------------------------------------------------%');
        fprintf('Now the frame %d is completed!\n',frame);
    end
%------------------------------------------------%    
% ITEM: id type x y z 
    now_frame = atom_data(:,:,frame);
    ID     = now_frame(:,1);
    TYPE   = now_frame(:,2);
    XYZ    = now_frame(:,3:5);
    %xl,yl,zl
%%----------------------------------------------%%
%%----------------------------------------------%%

bin_size  = xl/bin_num; % 

Block_Size_bins(frame) =bin_size;

bin_ALL_mesh_xsize = xl/bin_ALL_mesh;      % 能量/动能/势能-在x_bin格子大小
bin_ALL_mesh_ysize = yl/bin_ALL_mesh;      % 能量/动能/势能-在x_bin格子大小
bin_ALL_mesh_zsize = zl/bin_ALL_mesh;      % 能量/动能/势能-在z_bin格子大小

%%----------------------------------------------%% 
  for block = 1:bin_num   
      number_bin=0;
      for i = 1:Nlocal
%------------------------------------------------%
%------------------------------------------------%

          if((block==1))
%mesh temp of each atom in bin
                if(1)   
                    
                    x_n = floor((XYZ(i,1)-x_bound(1,1))/bin_ALL_mesh_xsize)+1;
                    y_n = floor((XYZ(i,2)-x_bound(1,1))/bin_ALL_mesh_ysize)+1;                    
                    z_n = floor((XYZ(i,3)-z_bound(1,1))/bin_ALL_mesh_zsize)+1;          
                    
                    if(x_n<=0)
                       x_n=1;
                    elseif(x_n>=bin_ALL_mesh)
                        x_n=bin_ALL_mesh;
                    end
                    
                    if(y_n<=0)
                       y_n=1;
                    elseif(y_n>=bin_ALL_mesh)
                        y_n=bin_ALL_mesh;
                    end          
                    
                    if(z_n<=0)
                       z_n=1;
                    elseif(z_n>=bin_ALL_mesh)
                       z_n=bin_ALL_mesh;   
                    end
                    
                     DEN_mesh_xy(x_n,y_n)=(DEN_mesh_xy(x_n,y_n)+1);
                     DEN_mesh_yz(y_n,z_n)=(DEN_mesh_yz(y_n,z_n)+1);
                     DEN_mesh_xz(x_n,z_n)=(DEN_mesh_xz(x_n,z_n)+1);
                end
          end %select atom
          
% here is to claculate density of all blocks

% XYZ(i,1) x方向
%------------------------------------------------%          
            if((XYZ(i,1)>(block-1)*(xl/bin_num)) && (XYZ(i,1)<=block*(xl/bin_num)))
                  
                DEN_bins_x(block,frame)=((DEN_bins_x(block,frame)+1));  

            end     
% XYZ(i,2) y方向           
%------------------------------------------------%            
            if((XYZ(i,2)>(block-1)*(yl/bin_num)) && (XYZ(i,2)<=block*(yl/bin_num)))
                  
                DEN_bins_y(block,frame)=((DEN_bins_y(block,frame)+1));  

            end  
% XYZ(i,3) z方向            
%------------------------------------------------%            
            if((XYZ(i,3)>(block-1)*(zl/bin_num)) && (XYZ(i,3)<=block*(zl/bin_num)))
                  
                DEN_bins_z(block,frame)=((DEN_bins_z(block,frame)+1));  

            end     

%------------------------------------------------%            
      end % loop all atom
  end 
  
    drawnow;
    
    x_real_length = ((1:length(DEN_bins_x)).*xl./bin_num)';
    y_real_length = ((1:length(DEN_bins_y)).*yl./bin_num)';
    z_real_length = ((1:length(DEN_bins_z)).*zl./bin_num)';
    %------------------------------------------------%  
    set(gcf, 'Position', [300 200 1500 800]); % gcf 获取当前figure的句柄

    subplot(2,3,1);  box on;hold on;
    title(strcat("frame: ",num2str(frame)));
    plot(x_real_length,(DEN_bins_x(:,frame)./bin_size./yl./zl),"ko-");
    xlim([0,max(x_real_length)]);
    xlabel('Length /Å','fontsize',20);
    ylabel('Density /Å^-3','fontsize',20);
    set(gca,'fontsize',20,'linewidth',1.5);
    legend ('Along X-axis','Location','southEast');        
    %------------------------------------------------%  
    set(gcf, 'Position', [300 200 1500 800]); % gcf 获取当前figure的句柄    
    subplot(2,3,2);  box on;hold on;
    plot(y_real_length,DEN_bins_y(:,frame)./bin_size./xl./zl,"r^-");
    xlim([0,max(x_real_length)]);
    xlabel('Length /Å','fontsize',20);
    ylabel('Density /Å^-3','fontsize',20);
    set(gca,'fontsize',20,'linewidth',1.5);
    legend ('Along Y-axis','Location','southEast');        
    %------------------------------------------------%  
    subplot(2,3,3);  box on;hold on;
    plot(z_real_length,DEN_bins_z(:,frame)./bin_size./xl./yl,"b*-");
    xlim([0,max(x_real_length)]);
    xlabel('Length /Å','fontsize',20);
    ylabel('Density /Å^-3','fontsize',20);
    set(gca,'fontsize',20,'linewidth',1.5);
    legend ('Along Z-axis','Location','southEast');
    %------------------------------------------------%  
    
    subplot(2,3,4);   box on;hold on;
    surf(DEN_mesh_xz./bin_ALL_mesh_xsize./bin_ALL_mesh_zsize./yl);
    shading interp;
    xlabel('Bin number','fontsize',20);
    ylabel('Bin number','fontsize',20);
    zlabel('Bin number','fontsize',20);
    
    
    set(gca,'fontsize',20,'linewidth',1.5);
    legend ('X-Z','Location','southEast');
    %------------------------------------------------%      
    subplot(2,3,5);  box on;hold on;
    surf(DEN_mesh_xy./bin_ALL_mesh_xsize./bin_ALL_mesh_ysize./zl);
    shading interp;
    
    xlabel('Bin number','fontsize',20);
    ylabel('Bin number','fontsize',20);
    zlabel('Bin number','fontsize',20);
    
    set(gca,'fontsize',20,'linewidth',1.5);
    legend ('X-Y','Location','southEast');
    %------------------------------------------------%  
    
    subplot(2,3,6);  box on;hold on;
    surf(DEN_mesh_yz./bin_ALL_mesh_ysize./bin_ALL_mesh_zsize./xl);
    shading interp;    
    
    
    xlabel('Bin number','fontsize',20);
    ylabel('Bin number','fontsize',20);
    zlabel('Bin number','fontsize',20);
    
    
    set(gca,'fontsize',20,'linewidth',1.5);
    legend ('Y-Z','Location','southEast');
    %------------------------------------------------%  

    
end

%----------------------------------------------%
%                  计算一维分布
%----------------------------------------------%
DEN_bins_x = DEN_bins_x./bin_size./yl./zl;
DEN_bins_y = DEN_bins_y./bin_size./xl./zl;
DEN_bins_z = DEN_bins_z./bin_size./xl./yl;

DEN_bins_x_ave = sum(DEN_bins_x,2)./all_frame;
DEN_bins_y_ave = sum(DEN_bins_y,2)./all_frame;
DEN_bins_z_ave = sum(DEN_bins_z,2)./all_frame;

%----------------------------------------------%
%                  计算二维分布
%----------------------------------------------%
% 注意这里计算的是X-Z平面上的投影，因此，计算数密度的时候需要除以yl
DEN_mesh_xz_ave =DEN_mesh_xz./bin_ALL_mesh_xsize./bin_ALL_mesh_zsize./yl./all_frame;

%计算数密度的时候需要除以zl
DEN_mesh_xy_ave =DEN_mesh_xy./bin_ALL_mesh_xsize./bin_ALL_mesh_ysize./zl./all_frame;

%计算数密度的时候需要除以xl
DEN_mesh_yz_ave =DEN_mesh_yz./bin_ALL_mesh_ysize./bin_ALL_mesh_zsize./xl./all_frame;
%----------------------------------------------%

disp("-------------------");
disp("----Finished!!!----");
disp("-------------------");