关注 $Mr.material,$ $更$$多$$精$$彩$！

---

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

$\text{ 专栏说明（订阅后可浏览对应专栏全部博文）：}$
$注意：$如需只订阅某个单独博文，请联系博主邮箱咨询。$[email\; protected]$

\spadesuit $†$ 开源后处理集成程序：请关注专栏《LAMMPS后处理——MATLAB子函数合集整理》
\spadesuit $†$$†$ 需要付费定制后处理程序请邮件联系：$[email\; protected]$

---

\spadesuit \spadesuit $Lammps$ 后处理及编程技巧第一次线上分享论坛（$\mathrm{2021.11.20}$）：
$视频链接$（请点击）
$资料下载$ (请点击)
$提取码：oaxk$

$\mathrm{欢迎关注公众号}:\mathrm{分子模拟}\mathrm{MD}\mathrm{Club}(Mr\_Materials)$

---

关注 $Mr.material,$ $更$$多$$精$$彩$！

博客介绍

这一部分的内容是独立于博客专栏（《LAMMPS后处理以及编程技巧》）之外的分支。在《LAMMPS后处理以及编程技巧》专栏中博文内容是针对于单个特殊的问题的，并且代码只是提供编程思路，并不是一个集成的子函数。
本博客将会介绍一些集成后子函数，可直接调用，用于Lammps的后处理，子函数内容和功能会及时更新，请关注博客更新时间。

更新于 2021-05-24
更新dump文件的读写函数：readdump_all(“your dump file”)，write_dump(fid,martix_name,all_dump,out_matrix,frame_flag);
更新于 2021-06-13
更新热力学量分析函数：thermo_output(all_data,target_name,dim,bin_number)
更新于 2021-06-14
更新热力学量分析函数：thermo_output(all_data,target_name,dim,bin_number) 目前包括，温度、动能，势能和数密度的一维和二维计算
更新于 2021-08-01
更新热力学量分析函数：thermo_output(all_data,target_name,dim,bin_number) ，修改了部分输出，增加了时间统计
更新于 2021-08-02
更新热力学量分析函数：thermo_output(all_data,target_name,dim,bin_number) ，修改了温度NAN bug以及输出问题

需要PDF资料请下载：

链接：请点击
提取码：b1nd

一、Dump轨迹的读入、写出

（1）读取Dump轨迹

NOTE!!!
Example:
Dump command must be (see below)

# 计算KE
compute          KE all    ke/atom 
# Dump出KE
dump        1    all custom 10000 dump.all                   &
                 id type x y z vx vy vz fx fy fz c_KE 

 all_dump = readdump_all("your dump file")

function [varargout] = readdump_all(varargin)
try
    dump = fopen(varargin{
    1},'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
%----------Outputs-------------
%OUTPUTS IN SAME VARIABLE STRUCTURE
varargout{
    1}.timestep = timestep;
varargout{
    1}.Natoms = Natoms;
varargout{
    1}.x_bound = x_bound;
varargout{
    1}.y_bound = y_bound;
varargout{
    1}.z_bound = z_bound;
varargout{
    1}.atom_data = atom_data;
%------------------------------
fclose(dump);

（2）写出Dump轨迹

Example:

martix_name = ["water_free_tag"];
frame_flag = frame;
all_dump.out_1 = water_curve_tag;
out_matrix = all_dump.atom_data(:,1:5);
out_matrix(:,6)=water_free_tag;
write_dump(fid,martix_name,all_dump,out_matrix,frame_flag);

function dump = write_dump(file_fid,martix_name,all_dump,out_matrix,frame_flag)
%{
    
%------------------------------------------------%
%      Author_Email: lammps_materials@163.com    %  
%------------------------------------------------%
file_fid   : the file fid for output by MATLAB;
martix_name : new column for dump data, like put new two columns into
dump,namely "watertag1 watertag2". Here the martix_name is
[watertag1;watertag2].Note that the name for each column is divided by ";".The
matrix faultly contains "id type x y z",therefore, the string name is
following: "ITEM: ATOMS id type x y z..."
matrix_data : all dump file includes the info: atom-id, type x y z... 
timestep    : you need to run "readlammps.m" fristly.Then you get timestep
matrix for all frame";
frame_flag  : Here is flag for each frame.please put the loop parameter
into there.
x_bound,y_bound,z_bound: Here is doundary info from readlammps.m. 
%------------------------------------------------%
%}
    title_first    = "ITEM: TIMESTEP";
    title_number   = "ITEM: NUMBER OF ATOMS";
    title_boundary = "ITEM: BOX BOUNDS pp pp pp";
    title_final = [];
    for ii = 1:max(size(martix_name))
        martix_name(ii)= martix_name(ii)+' ';
        title_final = strcat(title_final,martix_name(ii)) ; 
    end
    title_final      = strcat("ITEM: ATOMS id type x y z ",title_final);
       
    fprintf(file_fid,'%s\n',title_first);
    fprintf(file_fid,'%d\n',all_dump.timestep(frame_flag));
    
    fprintf(file_fid,'%s\n',title_number);
    fprintf(file_fid,'%d\n',length(all_dump.atom_data));
    
    fprintf(file_fid,'%s\n',title_boundary);
    fprintf(file_fid,'%d %d\n',[all_dump.x_bound(1,1),all_dump.x_bound(1,2)]);
    fprintf(file_fid,'%d %d\n',[all_dump.y_bound(1,1),all_dump.y_bound(1,2)]);
    fprintf(file_fid,'%d %d\n',[all_dump.z_bound(1,1),all_dump.z_bound(1,2)]);     
        
    fprintf(file_fid,'%s\n',title_final);
    [r,c]=size(out_matrix(:,:));            
        for i=1:r
            for j=1:c
                fprintf(file_fid,'%f\t',out_matrix(i,j));
            end
                fprintf(file_fid,'\r\n');
        end 
   
end

二、热力学量（温度、数密度、动能和势能）计算（一维线图、二维云图）

测试 case 和 code链接（MATLAB 2020a）：https://pan.baidu.com/s/1wXmMF9ZsCyCO2jMd8jgJNg
提取码：clw2

基本概念：（具体请参考）

Example
Keywords of “terget_name” can be “TEMP”, “KE”, “DEN” for temperature, kinetic energy，potential energy and number density for both one-dimension plot and two-dimension plot.

here is three atom type of SI, O and N, the keyowrd type should be slelected first. For example, the “KE” for x-z axis can be calcluated as " thermo_output(all_data,“KE”,[2,1,0,1],[50,50],[1,2,3]) ".
It means that all_data is your dump matrix, you need to run " all_dump = readdump_all(“your dump file”) " first. “KE” is keywod for kinetic energy.
[2,1,0,1] is dim, first 2 is dimension, [1,0,1] cooresponds the three axis, here x and z are selected.
[50,50] is mesh grid number for x and z axis.
[1,2,3] is all type form dump, which is defaultly form 0 to the max of atom type.

For simplicity, the results form matlab are compared with same paremeters of KE which is calculated by Ovito.

EXAMPLE:
all_data = readdump_all("your dump file")
out = thermo_output(all_data,"KE",[2,1,0,1],[50,50],[0,0,3])

function [varargout] = thermo_output(all_data,target_name,dim,bin_number,type)
%{
    
%------------------------------------------------%
%      Author_Email: lammps_materials@163.com    %  
%------------------------------------------------%
The dump file must be "ITEM: ATOMS id type x y z vx vy vz fx fy fz c_ke"
you should run "readdump_all.m" first!
The keywords for target_name can be: TEMP,KE,DEN;
dim: dimension,it should contain 4 parameter,e.i. [2,1,1,0]. 
     the first parameter is dimension,the next three parameters are calculating
     on which axis like [2,1,1,0] is dimension of 2 and x-y axis.If dim is
     1,the next three parameters should only one parameter is equal to 1,like
     [1,1,0,0], which means the simulation box will be sliced along x-axis
bin_number: bin_number is 1*2 matrix, If the dim is 2,the bin_number(2) is no-zero
    number; If the dim is 1,the bin_number(2) is zero number;
type: target atom type, If you have 5 types which are "Si,O,H,C,N" for
    example, the only Si and O is contributed to thermal value. the keywords 
    "type" should be [1,1,0,0,0]
out = thermo_output(all_data,"KE",[2,1,0,1],[50,50],[0,0,3])
%}
t_start_all = cputime;   % record time

if nargin < 5
   error('Error:Input parameter is smaller than 5!');      
elseif(size(bin_number)<2)
   error('Error:All dimension should be given!'); 
elseif(dim(1)==1 && ((bin_number(2))~=0)  )
   error('Error:the dim2 parameter of binnumber should be zero for the one dimension!');
elseif(dim(1)==1 && ((sum(dim(2:4)))>1)  )   
   error('Error:the only one axis should be given in one-dimension!'); 
 
   
elseif(dim(1)==2 && ((bin_number(2))==0)  ) 
   error('Error:the dim2 parameter of binnumber should be given for the two dimension!');    
elseif(dim(1)==2 && ((sum(dim(2:4)))>2)  )   
   error('Error:the two axis should be given in two-dimension!');   
end

if(target_name~="KE"&&target_name~="DEN"&&target_name~="TEMP")
   error('Error:target_name must be "KE", "DEN" or "TEMP"!');
end

disp("-----------------------------------------------------------------------------");

fprintf("Note:the dump file must be ITEM: ATOMS id type x y z vx vy vz fx fy fz c_ke!\n");
fprintf("Note:All unit must be in metal!\n");

disp("-----------------------------------------------------------------------------");
% KB 
KB        = 8.625*1e-5; % ev/k
%all frame number 
all_frame = size(all_data.atom_data,3);
% simulation box region
xl = all_data.x_bound(1,2)-all_data.x_bound(1,1);
yl = all_data.y_bound(1,2)-all_data.y_bound(1,1);
zl = all_data.z_bound(1,2)-all_data.z_bound(1,1);
%bin length for one dimension
unit_dim1_dx = floor(xl/bin_number(1));
unit_dim1_dy = floor(yl/bin_number(1));
unit_dim1_dz = floor(zl/bin_number(1));

unit_dim1    = floor(xl/bin_number(1));

%bin length for two dimension
unit_dim2_1  = floor(bin_number(1));
unit_dim2_2  = floor(bin_number(2));

%define plot matrix 
dim1_plot_matrix_x  = zeros(unit_dim1_dx,1);
dim1_plot_matrix_y  = zeros(unit_dim1_dy,1);
dim1_plot_matrix_z  = zeros(unit_dim1_dz,1);

%define plot temp matrix for one dimension
dim1_plot_matrix      = zeros(bin_number(1),all_frame);
dim1_plot_matrix_temp = zeros(bin_number(1),all_frame);

%define plot temp matrix for two dimension
dim2_mesh_matrix  = zeros(unit_dim2_1,unit_dim2_2);

dim2_mesh_matrix_temp     = zeros(unit_dim2_1,unit_dim2_2,all_frame);
dim2_mesh_matrix_temp_num = zeros(unit_dim2_1,unit_dim2_2,all_frame);

% mesh size
bin_mesh_xsize    = xl/bin_number(1);
bin_mesh_ysize    = yl/bin_number(2);
bin_mesh_zsize    = zl/bin_number(2);

flag_dim1=0;flag_dim2=0;

type_flag = find(type==1);

now_test       = all_data.atom_data(:,:,1);
TYPE_test      = now_test(:,2);

if(max(TYPE_test)<max(type))
   error('Error:Type is out of range!');
end

for frame = 1:all_frame
    
    now_frame = all_data.atom_data(:,:,frame);
    ID        = now_frame(:,1);
    TYPE      = now_frame(:,2);
    XYZ       = now_frame(:,3:5);
    KE        = now_frame(:,12);
    
  for i = 1:length(XYZ) 
      if((TYPE(i)==type)==0)
         continue; 
      end
      if (dim(1)==1)
          %------------------------------------------------------------%
          block = bin_number(1);
          if(dim(2)==1 && dim(3)==0 && dim(4)==0)    %x-axis
             bin_size = xl/block;
             flag_dim1=1;
             for block =1:bin_number(1)
                 if((XYZ(i,1)>(block-1)*bin_size) && (XYZ(i,1)<=block*bin_size))
                     if (strcmp(target_name,'TEMP'))
                         TEMP_now = KE(i)/1.5/KB;
                         dim1_plot_matrix(block,frame)=dim1_plot_matrix(block,frame)+TEMP_now;
                         dim1_plot_matrix_temp(block,frame)=dim1_plot_matrix_temp(block,frame)+1;
                     elseif(strcmp(target_name,'KE'))
                         dim1_plot_matrix(block,frame)=dim1_plot_matrix(block,frame)+KE(i);
                     elseif(strcmp(target_name,'DEN'))
                         dim1_plot_matrix(block,frame)=(dim1_plot_matrix(block,frame)+1)./(zl*yl*bin_size);
                     end  
                 end 
             end
          elseif(dim(2)==0 && dim(3)==1 && dim(4)==0)%y-axis
            flag_dim1=2;
            bin_size = yl/block;
            for block =1:bin_number(1)  
                 if((XYZ(i,2)>(block-1)*bin_size) && (XYZ(i,2)<=block*bin_size))
                     if (strcmp(target_name,'TEMP'))
                         TEMP_now = KE(i)/1.5/KB;
                         dim1_plot_matrix(block,frame)=dim1_plot_matrix(block,frame)+TEMP_now;
                         dim1_plot_matrix_temp(block,frame)=dim1_plot_matrix_temp(block,frame)+1;
                     elseif(strcmp(target_name,'KE'))
                         dim1_plot_matrix(block,frame)=dim1_plot_matrix(block,frame)+KE(i);
                     elseif(strcmp(target_name,'DEN'))
                         dim1_plot_matrix(block,frame)=(dim1_plot_matrix(block,frame)+1)./(xl*zl*bin_size);
                     end  
                 end
            end
          elseif(dim(2)==0 && dim(3)==0 && dim(4)==1)%z-axis
            flag_dim1=3;bin_size = zl/block;
            for  block =1:bin_number(1) 
                 if((XYZ(i,3)>(block-1)*bin_size) && (XYZ(i,3)<=block*bin_size))
                     if (strcmp(target_name,'TEMP'))
                         TEMP_now = KE(i)/1.5/KB;
                         dim1_plot_matrix(block,frame)=dim1_plot_matrix(block,frame)+TEMP_now;
                         dim1_plot_matrix_temp(block,frame)=dim1_plot_matrix_temp(block,frame)+1;
                     elseif(strcmp(target_name,'KE'))
                         dim1_plot_matrix(block,frame)=dim1_plot_matrix(block,frame)+KE(i);
                     elseif(strcmp(target_name,'DEN'))
                         dim1_plot_matrix(block,frame)=(dim1_plot_matrix(block,frame)+1)./(xl*yl*bin_size);
                     end  
                 end
            end
          end
          %------------------------------------------------------------%
      elseif(dim(1)==2)
          x_n = floor((XYZ(i,1)-all_data.x_bound(1,1))/bin_mesh_xsize);
          y_n = floor((XYZ(i,2)-all_data.y_bound(1,1))/bin_mesh_ysize);
          z_n = floor((XYZ(i,3)-all_data.z_bound(1,1))/bin_mesh_zsize);
          if(x_n<=0)
             x_n=1;
          end
          if(x_n>=bin_number(1))
             x_n =bin_number(1);
          end
          if(y_n<=0)
             y_n=1;
          end
          if(y_n>=bin_number(2))
             y_n =bin_number(2);
          end
          if(z_n<=0)
             z_n=1;
          end
          if(z_n>=bin_number(2))
             z_n=bin_number(2); 
          end
          if(dim(2)==1 && dim(3)==1 && dim(4)==0)    %x-y axis
             flag_dim2=1;
             if (strcmp(target_name,'TEMP'))
                 TEMP_now = KE(i)/1.5/KB;
                 dim2_mesh_matrix_temp(x_n,y_n,frame)=(dim2_mesh_matrix(x_n,y_n,frame)+TEMP_now);
                 dim2_mesh_matrix_temp_num(x_n,y_n,frame)=dim2_mesh_matrix_temp(x_n,y_n,frame)+1;
             elseif(strcmp(target_name,'KE'))
                 dim2_mesh_matrix(x_n,y_n)=(dim2_mesh_matrix(x_n,y_n)+KE(i));
             elseif(strcmp(target_name,'DEN'))
                 dim2_mesh_matrix(x_n,y_n)=((dim2_mesh_matrix(x_n,y_n)+1))./zl;
             end
          elseif(dim(2)==0 && dim(3)==1 && dim(4)==1)%y-z axis
             flag_dim2=2;
             if (strcmp(target_name,'TEMP'))
                 TEMP_now = KE(i)/1.5/KB;
                 dim2_mesh_matrix_temp(y_n,z_n,frame)=(dim2_mesh_matrix(y_n,z_n,frame)+TEMP_now);
                 dim2_mesh_matrix_temp_num(y_n,z_n,frame)=dim2_mesh_matrix_temp(y_n,z_n,frame)+1;
             elseif(strcmp(target_name,'KE'))
                 dim2_mesh_matrix(y_n,z_n)=(dim2_mesh_matrix(y_n,z_n)+KE(i));
             elseif(strcmp(target_name,'DEN'))
                 dim2_mesh_matrix(y_n,z_n)=(dim2_mesh_matrix(y_n,z_n)+1)./xl;
             end
          elseif(dim(2)==1 && dim(3)==0 && dim(4)==1)%x-z axis
               flag_dim2=3;
             if (strcmp(target_name,'TEMP'))
                 TEMP_now = KE(i)/1.5/KB;
                 
                 
                 
                 
                 dim2_mesh_matrix_temp(x_n,z_n,frame)=(dim2_mesh_matrix_temp(x_n,z_n,frame)+TEMP_now);
                 dim2_mesh_matrix_temp_num(x_n,z_n,frame)=dim2_mesh_matrix_temp_num(x_n,z_n,frame)+1;
             
             elseif(strcmp(target_name,'KE'))
                 dim2_mesh_matrix(x_n,z_n)=(dim2_mesh_matrix(x_n,z_n)+KE(i));
                 
             elseif(strcmp(target_name,'DEN'))
                 dim2_mesh_matrix(x_n,z_n)=(dim2_mesh_matrix(x_n,z_n)+1)./ly;
             end
          end
      end
  end
  
end
% output dim1_plot_matrix_temp(block,frame)=dim1_plot_matrix_temp(block,frame)+1;
    if(dim(1)==1&& strcmp(target_name,'TEMP')~=1)
        dim1_plot_matrix_new =sum(dim1_plot_matrix,2)./size(dim1_plot_matrix,2);
        x = [0 length(dim1_plot_matrix_new)];
        y = [mean(dim1_plot_matrix_new) mean(dim1_plot_matrix_new)];
        plot(dim1_plot_matrix_new,'bo-');
        line(x,y,'Color','red','LineStyle','--','Linewidth',2);
        
        legend(target_name,'Averaged value');
        xlabel('Axis /Å');
        ylabel(target_name);
    elseif(dim(1)==1 && strcmp(target_name,'TEMP')==1)    
        dim1_plot_matrix_new =sum(dim1_plot_matrix./size(dim1_plot_matrix,2)./dim1_plot_matrix_temp,2);
        dim1_plot_matrix_new(find(isnan(dim1_plot_matrix_new)==1)) = 0;
        plot(dim1_plot_matrix_new,'bo-');hold on;
        
        x = [0 length(dim1_plot_matrix_new)];
        y = [mean(mean(dim1_plot_matrix_new)) mean(mean(dim1_plot_matrix_new))];
        line(x,y,'Color','red','LineStyle','--','Linewidth',2);
        if (flag_dim1==1)
            xlabel('X-Axis /bin\_number');
        elseif(flag_dim1==2)
            xlabel('Y-Axis /bin\_number');
        elseif(flag_dim1==3)
            xlabel('Z-Axis /bin\_number');
        end
        ylabel(target_name);    
        legend('Temperature','Averaged value');
        set(gca, 'fontsize', 16);
       
        set(gca,'fontname','Times');
    elseif(dim(1)==2 && strcmp(target_name,'TEMP')==1)
        
        dim2_mesh_matrix = sum(dim2_mesh_matrix_temp,3)./all_frame;
        surf(dim2_mesh_matrix);view(2);
        if (flag_dim2==1)
            ylabel('X-Axis /bin\_number');
            xlabel('Y-Axis /bin\_number');
        elseif(flag_dim2==2)
            ylabel('Y-Axis /bin\_number');
            xlabel('Z-Axis /bin\_number');
        elseif(flag_dim2==3)
            ylabel('X-Axis /bin\_number');
            xlabel('Z-Axis /bin\_number');
        end
    elseif(dim(1)==2 && strcmp(target_name,'TEMP')~=1)   
        dim2_mesh_matrix = dim2_mesh_matrix./all_frame;
        surf(dim2_mesh_matrix);
        if (flag_dim2==1)
            ylabel('X-Axis /bin\_number');
            xlabel('Y-Axis /bin\_number');
        elseif(flag_dim2==2)
            ylabel('Y-Axis /bin\_number');
            xlabel('Z-Axis /bin\_number');
        elseif(flag_dim2==3)
            ylabel('X-Axis /bin\_number');
            xlabel('Z-Axis /bin\_number');
        end        
        shading interp; view(2);
    end
    
varargout{
    1}.all_data=all_data;     
if(bin_number(2)==0) 
varargout{
    1}.one_dim_mesh_plot=dim1_plot_matrix;
else
varargout{
    1}.two_dim_mesh_plot=dim2_mesh_matrix;    
end
% time record
%--------------------------------------------------%
t_start_stop = (cputime-t_start_all);
if(t_start_stop<1)
    fprintf('Now the used time is: %.1f s.\n',t_start_stop);
else
    fprintf('Now the used time is: %.1f mins.\n',t_start_stop/60);
end
%--------------------------------------------------%
disp("----------------");
fprintf("----ALL DONE----\n");
disp("----------------");
end

（Continue…）

三、结构分析

径向分布

回旋半径

配位数

(Continue…)

四、动力学分析

速度自关联分析

(continue…)