轨迹分析—Matlab计算均方位移

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

---

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

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

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

---

前言

$热力学性质的不均匀性导致的热力学过程叫做输运过程，相应的现象叫做输运现象。$
$由于粒子数密度的不均匀性导致粒子的输运（扩散现象）.$$均方根位移是随时间测量的，以确定粒子是否仅由于扩散而扩散。$
$\mathrm{那么如何根据轨迹文件计算}\mathrm{MSD}呢？$

$Dr.Xian\mathrm{总结了计算}\mathrm{MSD}的常用方法：$
扩展阅读：《5种方法计算均方位移MSD》

一、请先试用旧代码

clc;clear;
t_start_all = cputime;
%----------Step 1-------------
% Step 1: read dump
file ="NaCl_DuiBi.xyz";
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');                     % cancel按钮设置
set(hBtn_3, 'String', 'cancle', 'FontSize', 10);
i=1;                                                                   
while feof(dump) == 0                                                       % read to end of file
    compT = i/202;                                                          % 这里需要手动设置一下分母，为帧数+1
    waitbar(compT, ReadDump_Waitbar, ['Dump info reading...', num2str(round(compT, 2) * 100), '%']);
    pause(0.01); 
    id = fgetl(dump);                                                       % back the next line content
     if (strncmpi(id,'ITEM: TIMESTEP',numel('ITEM: TIMESTEP')))
            timestep(i) = str2num(fgetl(dump));                             % current timestep
    else
     if (strncmpi(id,'ITEM: NUMBER OF ATOMS',numel('ITEM: NUMBER OF ATOMS')))
            Natoms(i) = str2num(fgetl(dump));                               % number of atoms
     else
      if (strncmpi(id,'ITEM: BOX BOUNDS',numel('ITEM: BOX BOUNDS')))        
            x_bound(i,:) = str2num(fgetl(dump));                            % n*2
            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)
frame_num = i-1;
num_atoms = j;
xl = x_bound(1,2)-x_bound(1,1);                                             % x length of box
yl = y_bound(1,2)-y_bound(1,1);                                             % y length of box      
zl = z_bound(1,2)-z_bound(1,1);                                             % z length of box
% ----------Step 2-------------
% 位移修正，记录为real_data
% 设置DispCor_waitbar
DispCor_Waitbar = waitbar(0, 'Displacement correction...', 'CreateCancelBtn', 'delete(gcbf)');% 进度条初始设置
set(DispCor_Waitbar, 'Color', [0.9, 0.9, 0.9]);
hBtn_1 = findall(DispCor_Waitbar, 'type', 'Uicontrol');                     % cancel按钮设置
set(hBtn_1, 'String', 'cancle', 'FontSize', 10);

% 位移修正
real_data = atom_data;
for i1 = 1:frame_num
    real_data(:,:,i1) = sortrows(real_data(:,:,i1),1);                      % 先调整一下real_data的排序
end

for k = 1: frame_num-1
    compT_1 = k/(frame_num-1);
    waitbar(compT_1, DispCor_Waitbar, ['Displacement correction...', num2str(round(compT_1, 2) * 100), '%']);
    pause(0.01);
    current_data    = atom_data(:,1:5,k+1);
    current_data    = sortrows(current_data,1);  
    current_x = current_data(:,3);
    current_y = current_data(:,4);
    current_z = current_data(:,5);
    last_data = atom_data(:,1:5,k);
    last_data = sortrows(last_data,1);  
    last_x = last_data(:,3);
    last_y = last_data(:,4);
    last_z = last_data(:,5);
    for n = 1:num_atoms
        d_x(n) = current_x(n) - last_x(n);
        d_y(n) = current_y(n) - last_y(n);
        d_z(n) = current_z(n) - last_z(n);
        [d_x(n),d_y(n),d_z(n)] = PBC(d_x(n),d_y(n),d_z(n),xl,yl,zl);

        real_data(n,3,k+1) = real_data(n,3,k)+ d_x(n);                      % 修正x
        real_data(n,4,k+1) = real_data(n,4,k)+ d_y(n);                      % 修正y
        real_data(n,5,k+1) = real_data(n,5,k)+ d_z(n);                      % 修正z        
    end
end
close(DispCor_Waitbar)
% ----------Step 3-------------
% MSD计算
mean_squard_dx_2 = zeros(frame_num-1,1);
mean_squard_dy_2 = zeros(frame_num-1,1);
mean_squard_dz_2 = zeros(frame_num-1,1);
mean_squard_all_2 = zeros(frame_num,1);
% 设置MSD_waitbar
MSD_Waitbar = waitbar(0, 'MSD computing...', 'CreateCancelBtn', 'delete(gcbf)');% 进度条初始设置
set(MSD_Waitbar, 'Color', [0.9, 0.9, 0.9]);
hBtn_2 = findall(MSD_Waitbar, 'type', 'Uicontrol');                             % cancel按钮设置
set(hBtn_2, 'String', 'cancle', 'FontSize', 10);

% 循环不同的时间间隔，最大为最后一帧减第一帧为frame_num-1
for tau = 1:frame_num-1 
    compT_2 = tau/(frame_num-1);
    waitbar(compT_2, MSD_Waitbar, ['Different tau value Computing...', num2str(round(compT_2, 2) * 100), '%']);
    pause(0.01);
    % 循环计算公式中的不同的分子项数，与时间间隔tau有关，最大等于frame_num-tau
    for frame = 1: frame_num-tau
        now_frame    = real_data(:,1:5,frame + tau);
        % now_frame    = sortrows(now_frame,1);  
        now_x = now_frame(:,3);
        now_y = now_frame(:,4);
        now_z = now_frame(:,5);
        before_frame = real_data(:,1:5,frame);
        % before_frame = sortrows(before_frame,1);  
        before_x = before_frame(:,3);
        before_y = before_frame(:,4);
        before_z = before_frame(:,5);
        for n = 1:num_atoms
            dx(n) = now_x(n) - before_x(n);
            dy(n) = now_y(n) - before_y(n);
            dz(n) = now_z(n) - before_z(n);
            squard_dx(n,frame,tau) =  dx(n).^2;
            squard_dy(n,frame,tau) = dy(n).^2;
            squard_dz(n,frame,tau) = dz(n).^2;                              
        end
        
        % 第一次平均是对原子平均
        mean_squard_dx_1(frame,tau) = mean(squard_dx(:,frame,tau)); 
        mean_squard_dy_1(frame,tau) = mean(squard_dy(:,frame,tau)); 
        mean_squard_dz_1(frame,tau) = mean(squard_dz(:,frame,tau)); 
    end
   
    % 第二次是对分子中加项数平均
    mean_squard_dx_2(tau) = sum(mean_squard_dx_1(:,tau))./(frame_num - tau);
    mean_squard_dy_2(tau) = sum(mean_squard_dy_1(:,tau))./(frame_num - tau);
    mean_squard_dz_2(tau) = sum(mean_squard_dz_1(:,tau))./(frame_num - tau);
    mean_squard_all_2(tau+1) = mean_squard_dx_2(tau) + mean_squard_dy_2(tau) + mean_squard_dz_2(tau);
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("-------------------");

$以下代码先不要看需要检查一下$

一、读取Dump轨迹

这里不再赘述，其他读取方法$（\mathrm{Python}\mathrm{和}\mathrm{C}++）$可参考博文《LAMMPS后处理以及编程技巧》的编程技巧篇。

clc;clear;
file ="test.xyz";

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

i=1;now_frame=0;
while feof(dump) == 0
    id = fgetl(dump);
     if (strncmpi(id,'ITEM: TIMESTEP',numel('ITEM: TIMESTEP')))
            now_frame = now_frame+1;
                
            if(rem(now_frame,100)==0)
                disp('%------------------------------------------------%');
                fprintf('Now the reading frame is: %d.\n',now_frame);
            end
            
            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("----ALL DONE!!!----");
disp("-------------------");

二、坐标矩阵处理子函数

$1.函数调用$

$\mathrm{其中选择原子}Type10,11,12,13$;
$\mathrm{xl};\mathrm{yl};\mathrm{zl}\mathrm{分别为模拟盒子大小};$

[real_data] = re_define_matrix(atom_data,1,[10,11,12,13],xl,yl,zl);

function [real_data] = re_define_matrix(atom_data,nc,type,xl,yl,zl)
%
%
    real_data=zeros(size(atom_data));
    TYPE=atom_data(:,2,1);
    N_local=0;
    for i = 1:length(atom_data)
       if((TYPE(i)==type)~=0)
        N_local = N_local+1;
       end
    end

    for m = 1: size(atom_data,3)-nc
        
       if((TYPE(m)==type)~=0)             
           continue;
       end
   
        current_data    = atom_data(:,1:5,m+1);
        current_data    = sortrows(current_data,1);  
        current_x = current_data(:,3);
        current_y = current_data(:,4);
        current_z = current_data(:,5);
        last_data = atom_data(:,1:5,m);
        last_data = sortrows(last_data,1);  
        last_x = last_data(:,3);
        last_y = last_data(:,4);
        last_z = last_data(:,5);
        
        for n = 1:length(atom_data)
           if((TYPE(n)==type)~=0)             
               continue;
           end              
            d_x(n) = current_x(n) - last_x(n);
            d_y(n) = current_y(n) - last_y(n);
            d_z(n) = current_z(n) - last_z(n);
            [d_x(n),d_y(n),d_z(n)] = PBC(d_x(n),d_y(n),d_z(n),xl,yl,zl);

            real_data(n,3,m+1) = real_data(n,3,m)+ d_x(n);                     
            real_data(n,4,m+1) = real_data(n,4,m)+ d_y(n);                     
            real_data(n,5,m+1) = real_data(n,5,m)+ d_z(n);                       
        end
    end 
end

三、计算MSD子函数

[mean_squard_all_2] = find_msd(redefine_data,nc);

function [out_mean_squard] = find_msd(redefine_data,nc)

    frame_num=size(redefine_data,3);
 for tau = 1:frame_num-nc
     
    frame = 1: frame_num-tau;      
    dx = redefine_data(:,3,frame+tau) - redefine_data(:,3,frame);
    dy = redefine_data(:,4,frame+tau) - redefine_data(:,4,frame);
    dz = redefine_data(:,5,frame+tau) - redefine_data(:,5,frame);
    squard_dx(:,frame,tau) =  dx.^2;
    squard_dy(:,frame,tau) =  dy.^2;
    squard_dz(:,frame,tau) =  dz.^2;  

    mean_squard_dx_1(frame,tau) = mean(squard_dx(:,frame,tau)); 
    mean_squard_dy_1(frame,tau) = mean(squard_dy(:,frame,tau)); 
    mean_squard_dz_1(frame,tau) = mean(squard_dz(:,frame,tau)); 

    mean_squard_dx_2(tau) = sum(mean_squard_dx_1(:,tau))./(frame_num - tau);
    mean_squard_dy_2(tau) = sum(mean_squard_dy_1(:,tau))./(frame_num - tau);
    mean_squard_dz_2(tau) = sum(mean_squard_dz_1(:,tau))./(frame_num - tau);
    mean_squard_all_2(tau+nc) = mean_squard_dx_2(tau) + mean_squard_dy_2(tau) + mean_squard_dz_2(tau);
    out_mean_squard(tau+nc) = mean_squard_all_2(tau+1);
 end
 out_mean_squard=out_mean_squard';
end

四、周期性边界处理

%pbc 
function [dx,dy,dz] = PBC(dx,dy,dz,xl,yl,zl)
%-------------------------------------------------%
            if (dx>xl/2)
                dx = dx - xl;
            elseif  (dx<(-xl/2))
                dx = dx + xl;
            end
            if (dy>yl/2)
                dy = dy - yl;
            elseif  (dy<(-yl/2))
                dy = dy + yl;
            end
            if (dz>zl/2)
                dz = dz - zl;
            elseif  (dz<(-zl/2))
                dz = dz + zl;
            end
%-------------------------------------------------%
end

五、主函数

%{
MSD = 1/N*sum_1_to_N(r(t)-r(0))^2
%}

clc;clear;
file ="test.xyz";

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

i=1;now_frame=0;
while feof(dump) == 0
    id = fgetl(dump);
     if (strncmpi(id,'ITEM: TIMESTEP',numel('ITEM: TIMESTEP')))
            now_frame = now_frame+1;
                
            if(rem(now_frame,100)==0)
                disp('%------------------------------------------------%');
                fprintf('Now the reading frame is: %d.\n',now_frame);
            end
            
            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("----ALL DONE!!!----");
disp("-------------------");

% input_matrix 
%====================================================%
frame_num   = size(atom_data,3);

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);
msd_all_save = zeros(frame_num,1);
num_atoms = length(atom_data);

for nc = 1

    [real_data] = re_define_matrix(atom_data,1,[10,11,12,13],xl,yl,zl);
    %====================================================%
    % here is to calculate msd 
    redefine_data=real_data;

    [mean_squard_all_2] = find_msd(redefine_data,nc);
    msd_all_save(:,nc)=mean_squard_all_2;
    %====================================================%
end

六、LAMMPS VS MATLAB

七、Python绘图代码及数据

import os
def make_save_file(filepath,filename,savename):
    # filepath 存储路径； filename：创建文件夹的名字 savename：存储图片的名字
    save_path = filepath+os.sep+filename+os.sep+savename
    all_path  = filepath+os.sep+filename
    if not os.path.exists(all_path):
        os.mkdir(all_path)
    plt.savefig(save_path,dpi=1000)

filepath = r"MSD"

font1 = {'family': 'Times New Roman',
             'weight': 'bold',
             'size': 22,
             }
font2 = {'family': 'Times New Roman',
             'weight': 'normal',
             'size': 11,
             }
# Global setting
# ----------------------------------------------------------------------#
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
from scipy import optimize 
# coding:utf-8
import pylab
import matplotlib.ticker as plticker
from matplotlib.pyplot import MultipleLocator
import xlrd
import matplotlib.pyplot as plt
from matplotlib.lines import Line2D
plt.style.use('science')
# ----------------------------------------------------------------------#
def setup(ax,x_label,y_label):
    # 边框线的线宽设置
    width = 1.5
    ax.spines['top'].set_linewidth(width)
    ax.spines['bottom'].set_linewidth(width)
    ax.spines['left'].set_linewidth(width)
    ax.spines['right'].set_linewidth(width)
    # 边框上的ticks的出现
    ax.tick_params(top='on', bottom='on', left='on', right='on', direction='in')
    # 边框上的ticks对应的lable,即数字的尺寸
    ax.tick_params(labelsize=20,pad = 6)
    #ax.yaxis.set_ticks_position('right')
    ax.tick_params(which='major', width=1.00, length=5)
    ax.tick_params(which='minor', width=0.75, length=2)
    # 边框上的ticks的出现的间隔
    locx = plticker.MultipleLocator(base=100000) # this locator puts ticks at regular intervals
    ax.xaxis.set_major_locator(locx)
    locy = plticker.MultipleLocator(base=2000) # this locator puts ticks at regular intervals
    ax.yaxis.set_major_locator(locy)
    # 边框上的注释labels
    labels = ax.get_xticklabels() + ax.get_yticklabels()
    [label.set_fontname('Times New Roman') for label in labels]
    # 边框上的注释labels距离坐标轴的距离
    xAxisLable = x_label 
    yAxisLable = y_label
    ax.set_xlabel(xAxisLable, font1, labelpad=5)
    ax.set_ylabel(yAxisLable, font1, labelpad=5)
# ----------------------------------------------------------------------#
# 读取excel
excel_path = r"MSD.xlsx"
work_sheet = xlrd.open_workbook(excel_path);

# sheet的数量
sheet_num = work_sheet.nsheets;
# 获取sheet name
sheet_name = []
for sheet in work_sheet.sheets():
    sheet_name.append(sheet.name)
# ----------------------------------------------------------------------#
# ----------------------------------------------------------------------#
for sheet_num in range(1):
    # 选取sheet
    fig, axe = plt.subplots(1, 1, figsize=(5, 5))
    now_sheet = work_sheet.sheets()[0]

    # sheet的行
    row = now_sheet.nrows
    # sheet的列
    ncol = now_sheet.ncols
    #----------------------------------------------------------------------#
    #----------------------------------------------------------------------#
    font1 = {'family': 'Times New Roman',
             'weight': 'bold',
             'size': 18,
             }
    font2 = {'family': 'Times New Roman',
             'weight': 'normal',
             'size': 11,
             }

    time          = now_sheet.col_values(0)
    time_matrix   = time[1:] 
    
    lammps_msd          = now_sheet.col_values(1)
    lammps_msd_matrix   = lammps_msd[1:]    
    
    matlab_msd          = now_sheet.col_values(2)
    matlab_msd_matrix   = matlab_msd[1:]               

       
    #----------------------------------------------------------------------#

    
    markersize = 10                                          
    linewidth = 2
    markevery =10
    alpha =0.5 
    
    
    p_1 = dict(marker='o',color = 'b',linestyle = 'none',markevery=markevery,markerfacecolor='none',markersize=12,linewidth=linewidth,alpha=1)
    p_2 = dict(marker='o',color = 'r',linestyle = 'none',markevery=markevery,markerfacecolor='r',markersize=markersize,linewidth=linewidth,alpha=alpha)


    #----------------------------------------------------------------------#    
    markersize = 10
    
    axe.plot(time_matrix, lammps_msd_matrix, **p_1,label='Lammps')   
    axe.plot(time_matrix, matlab_msd_matrix, **p_2,label='Matlab')   

    
    
    
    axe.legend(loc='best', frameon=False, labelspacing=0.3)
    axe.legend(bbox_to_anchor=(0.5, 1), loc='upper left', borderaxespad=0.,fontsize='x-large')
    axe.set_xlim(0,200000)
    axe.set_ylim(0, 4000)
    setup(axe,r'$\rm Timestep$ ',r'$\rm MSD$') 
    # ============================================================#
    filename = "python"
    savename = "msd.jpg"
    
    make_save_file(filepath,filename,savename)

plt.show()    

八、代码下载

链接：请点击
提取码：2u8y