brucefan1983的个人博客分享 http://blog.sciencenet.cn/u/brucefan1983

博文

A C code for calculating diffusion coefficient in MD

已有 3955 次阅读 2016-7-20 22:54 |个人分类:计算物理基础|系统分类:科研笔记|关键词:分子动力学模拟

GK_vs_Einstein.pdf

该pdf附件给出了一个爱因斯坦公式与格林-久保公式等价的较为详细的证明。

md_diffusion.c

plot_results.m

vac.txt

msd.txt

以上是C语言源代码、matlab后处理脚本、以及我的计算数据。C语言编译命令可以为:

gcc -std=c99 -O3 -Wall md_diffusion.c -lm

对应的运行命令则为

./a.out

C语言程序在我的笔记本电脑中运行时间为25秒(扩散系数的计算在七八十年代研究得比较多。当时需要在大型计算机花几天甚至几个月完成的工作现在可能只需要在一个笔记本电脑上花几分钟就能搞定)。读者可先运行C语言程序,再运行matlab画图脚本,但也可以先不运行C语言程序,而直接运行matlab脚本,因为我也上传了我计算的数据。由于用了随机数,读者得到的数据不一定与我得到的一致,但应该与我的相当。实际上,多次模拟求平均值是很多模拟中必须做的事情,我就不多说了。



                          不以结婚为目的的谈恋爱都是耍流氓?

                          不给源代码的计算物理书都是耍流氓!

                用分子动力学模拟计算扩散系数

一、引言 

  还记得物理系大一《热学》课中学过的三种典型的输运现象吗?它们是扩散、热传导、和粘滞现象。我在上一篇博文介绍了用平衡态分子动力学模拟计算热传导系数(即热导率)的方法。本文则介绍用平衡态分子动力学模拟计算扩散系数的方法。我将比较两种方法:一种是基于爱因斯坦公式的;一种是基于格林-久保(Green-Kubo)公式的。在爱因斯坦公式中,跑动扩散系数(RDC  =  running diffusion coefficient) 可以写为均方位移  (MSD  = mean  square  displacement)的时间导数;在格林-久保公式中,跑动扩散系数可以写为速度自关联(VAC = velocity auto-correlation)函数的时间积分。 我将用数学推导(见开头的pdf附件)和数值实验证明:这两种方法是等价的。

二、什么是扩散和扩散系数?

  我们知道,不管何种系统,在不受外力(准确地说是外场)作用时,若其内部有热力学性质的不均匀性 (inhomogeneity),则它一定处于非平衡的状态,并有向平衡态靠近的趋势。这是热力学第二定律的结果。这种由热力学性质的不均匀性导致的热力学过程叫做输运过程 (transport process),相应的现象叫做输运现象。前一篇博文讨论的热输运是由温度的不均匀性导致的;本文要讨论的扩散过程是由(粒子数)密度n(x,y,z)的不均匀性导致的。

  大家一定还记得有一个描述扩散现象的经验定律,那就是Fick定律。这个定律是说粒子流J_N——即单位时间穿过单位面积的粒子数——正比于数密度梯度dn/dx(为简单起见,我们假定输运沿着x方向):

         J_N = -D dn/dx         

这里,D是扩散的输运系数,叫做扩散系数 (diffusion coefficient),其中负号表示粒子向数密度小的地方扩散。

三、计算扩散系数的爱因斯坦公式  
  在分子动力学模拟中,我们可以用爱因斯坦公式计算自扩散系数(self-diffusion coefficient)。为简单起见,以下简称自扩散系数为扩散系数。

  爱因斯坦的公式那么多,究竟是哪个爱因斯坦公式呢?就是那个关于布朗运动的爱因斯坦公式啦。对于x方向的扩散,扩散系数可以表达为均方位移的时间导数:

    D_x (t) = (1/2)  d MSD(t) / dt.
其中,MSD(t) 是均方位移,它是关联时间 t 的函数。因为这里定义的扩散系数是关联时间的函数,所以我们叫它跑动扩散系数(running diffusion coefficient),意思是这个扩散系数的值是随关联时间的增大而“跑动”的。上面公式中的1/2不能丢掉。要特别注意的是,这里的t是关联时间(correlation time),不是分子动力学模拟中的模拟时间。

  均方位移的定义如下:

  MSD(t) = < (1/N) * sum_i^N  [x_i(t+t_0) - x_i(0+t_0)]^2 >

其中,N是粒子数,sum_i^N代表对粒子求和,x_i(t)是粒子i在t时刻的x坐标,尖括号<>代表系综平均,但在分子动力学模拟中由对“时间原点” t_0 的平均替代。

四、计算扩散系数的格林-久保公式

  从以上的爱因斯坦公式,可以推导出如下格林-久保公式(具体推导我下次传个pdf上来):

  D_x (t) = int_0^t dt  VAC(t).

其中,int_0^t dt代表对关联时间从0到t积分,而VAC(t)就是关联时间为t时的速度自关联函数。

  速度自关联函数的定义如下:

    VAC(t) = < (1/N) * sum_i^N  vx_i(t+t_0) * vx_i(0+t_0) >

这里,vx_i(t)是粒子i在t时刻的x方向的速度。其它符号的意义类似于MSD的情形。

  这里的格林不是“格林函数”中的那个格林,而久保(Kubo)就是那个发明线性响应理论的日本科学家(还是蛮有名的;固体物理和量子多体理论的书都少不了要介绍他的理论的)。

五、分子动力学模拟代码

  根据以上对均方位移和速度自关联函数的定义,再结合我上一篇博文给出的计算热导率的分子动力学模拟程序,很容易写出一个计算均方位移和速度自关联函数的分子动力学模拟程序。有一点要特别注意的是,在模拟流体(气体和液体)时,要定义两套坐标,对其中一套坐标施加一个周期边界条件(即每当某个粒子跑出模拟盒子,就把它拉回来;对于固体,这是没有必要的),对另一套坐标不用施加周期边界条件(即不用将跑出盒子的粒子拉回盒子,而是任其自由跑动,即允许体系“扩散”)。求力的时候,一定要用那套施加了周期边界条件的坐标(不然求的力是不对的);而计算MSD时,则用那套没有施加周期边界条件的坐标(不然计算的MSD是错的)。好了,废话少说,直接上代码吧。下面是一个完整的约500行的C语言程序。该程序可计算由Lennard-Jones势描述的流体系统的均方位移和速度自关联函数并将结果写入文件。稍后,我将给出一个matlab脚本,用来计算跑动扩散系数并画图。


/*****************************************************************************80
   Copyright 2016 Zheyong Fan <brucenju@gmail.com>

   This is the source code of MD_DIFFUSION, which can
   calculate the velocity auto-correlation (VAC) using the velocity data and
   calculate the mean square displacement (MSD) using the position data.

   MD_DIFFUSION is a free software: you can redistribute it and/or modify
   it under the terms of the GNU General Public License as published by
   the Free Software Foundation, either version 3 of the License, or
   (at your option) any later version.

   MD_DIFFUSION is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   GNU General Public License for more details.

   You can find a copy of the GNU General Public License at
   <http://www.gnu.org/licenses/>.
*******************************************************************************/


/*----------------------------------------------------------------------------80
   Author: Zheyong Fan <brucenju@gmail.com>
------------------------------------------------------------------------------*/


/*----------------------------------------------------------------------------80
The units used in this program are as follows:
   (1) The basic units are chosen to be
       (A) Mechanics
           -- energy: eV (electron volt)
           -- length: A (Angstrom)
           -- mass: amu (atomic mass unit)
       (B) Thermodynamics
           -- temperature: K (Kelvin)
   (2) Some derived units are
       (A) Mechanics
           -- time: A amu^(1/2) eV^(-1/2) ~ 10 fs
           -- velocity: eV^(1/2) amu^(-1/2) ~ 0.1 A fs^(-1)
           -- acceleration: eV A^(-1) amu^(-1) ~ 0.01 A fs^(-2)
           -- force: eV A^(-1)
       (B) Thermodynamics
           Boltzmann constant: K_B ~ 0.8625 * 10^(-4) eV K^(-1)
------------------------------------------------------------------------------*/


#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <time.h>

#define K_B                   8.625e-5 // Boltzmann's constant in my unit
#define TIME_UNIT_CONVERSION  1.018e+1 // fs <-> my unit


// minimum image convention
void apply_mic
(
   double lx, double ly, double lz, double lxh, double lyh,  
   double lzh, double *x12, double *y12, double *z12
)
{
   if (*x12 < - lxh) {*x12 += lx;} else if (*x12 > + lxh) {*x12 -= lx;}
   if (*y12 < - lyh) {*y12 += ly;} else if (*y12 > + lyh) {*y12 -= ly;}
   if (*z12 < - lzh) {*z12 += lz;} else if (*z12 > + lzh) {*z12 -= lz;}
}


// pull the atoms back to the simulation box (needed for simulating fluids)
void apply_pbc
(
   int N, double lx, double ly, double lz, double *x, double *y, double *z
)
{
   for (int n = 0; n < N; n++)
   {
       if (x[n] < 0) {x[n] += lx;} else if (x[n] > + lx) {x[n] -= lx;}
       if (y[n] < 0) {y[n] += ly;} else if (y[n] > + ly) {y[n] -= ly;}
       if (z[n] < 0) {z[n] += lz;} else if (z[n] > + lz) {z[n] -= lz;}
   }
}


// neighbor list construction (slow for large system, but ok for teaching)
void find_neighbor
(
   int N, int *NN, int *NL, double *x, double *y, double *z,  
   double lx, double ly, double lz, int MN, double cutoff
)              
{
   double lxh = lx * 0.5;
   double lyh = ly * 0.5;
   double lzh = lz * 0.5;  
   double cutoff_square = cutoff * cutoff;
   for (int n = 0; n < N; n++) {NN[n] = 0;}
   for (int n1 = 0; n1 < N - 1; n1++)
   {  
       for (int n2 = n1 + 1; n2 < N; n2++)
       {    
           double x12 = x[n2] - x[n1];
           double y12 = y[n2] - y[n1];
           double z12 = z[n2] - z[n1];
           apply_mic(lx, ly, lz, lxh, lyh, lzh, &x12, &y12, &z12);
           double  distance_square = x12 * x12 + y12 * y12 + z12 * z12;
           if (distance_square < cutoff_square)
           {        
               NL[n1 * MN + NN[n1]] = n2;
               NN[n1]++;
               NL[n2 * MN + NN[n2]] = n1;
               NN[n2]++;
           }
           if (NN[n1] > MN)
           {
               printf("Error: cutoff for neighbor list is too large.n");
               exit(1);
           }
       }
   }  
}


// position initialization (FCC crystal)
void initialize_position  
(
   int n0, int nx, int ny, int nz, double ax, double ay, double az,  
   double *x, double *y, double *z
)
{
   double x0[4] = {0.0, 0.0, 0.5, 0.5}; // only for simple FCC lattice
   double y0[4] = {0.0, 0.5, 0.0, 0.5};  
   double z0[4] = {0.0, 0.5, 0.5, 0.0};
   int n = 0;
   for (int ix = 0; ix < nx; ++ix)
   {
       for (int iy = 0; iy < ny; ++iy)
       {
           for (int iz = 0; iz < nz; ++iz)
           {
               for (int i = 0; i < n0; ++i)
               {
                   x[n] = (ix + x0[i]) * ax;
                   y[n] = (iy + y0[i]) * ay;
                   z[n] = (iz + z0[i]) * az;
                   n++;
               }
           }
       }
   }
}  


// velocity scaling (the simplest way of controlling temperature)
void scale_velocity
(int N, double T_0, double *m, double *vx, double *vy, double *vz)
{  
   double temperature = 0.0;
   for (int n = 0; n < N; ++n)  
   {
       double v2 = vx[n] * vx[n] + vy[n] * vy[n] + vz[n] * vz[n];      
       temperature += m[n] * v2;  
   }
   temperature /= 3.0 * K_B * N;
   double scale_factor = sqrt(T_0 / temperature);
   for (int n = 0; n < N; ++n)
   {  
       vx[n] *= scale_factor;
       vy[n] *= scale_factor;
       vz[n] *= scale_factor;
   }
}  
   

// velocity initialization (with zero total momentum)  
void initialize_velocity
(int N, double T_0, double *m, double *vx, double *vy, double *vz)
{
   double momentum_average[3] = {0.0, 0.0, 0.0};
   for (int n = 0; n < N; ++n)
   {  
       vx[n] = -1.0 + (rand() * 2.0) / RAND_MAX;  
       vy[n] = -1.0 + (rand() * 2.0) / RAND_MAX;  
       vz[n] = -1.0 + (rand() * 2.0) / RAND_MAX;    
       
       momentum_average[0] += m[n] * vx[n] / N;
       momentum_average[1] += m[n] * vy[n] / N;
       momentum_average[2] += m[n] * vz[n] / N;
   }  
   for (int n = 0; n < N; ++n)  
   {  
       vx[n] -= momentum_average[0] / m[n];
       vy[n] -= momentum_average[1] / m[n];
       vz[n] -= momentum_average[2] / m[n];  
   }
   scale_velocity(N, T_0, m, vx, vy, vz);
}


// force evaluation for Lennard-Jones potential
void find_force
(
   int N, int *NN, int *NL, int MN, double lx, double ly, double lz,
   double *x, double *y, double *z, double *fx, double *fy, double *fz,
   double cutoff
)
{
   // precompute something (a trick to make the code faster)
   const double epsilon = 1.032e-2;
   const double sigma = 3.405;
   const double sigma_3 = sigma * sigma * sigma;
   const double sigma_6 = sigma_3 * sigma_3;
   const double sigma_12 = sigma_6 * sigma_6;
   const double factor_1 = 24.0 * epsilon * sigma_6;  
   const double factor_2 = 48.0 * epsilon * sigma_12;

   // initialize force before accumulation
   for (int n = 0; n < N; ++n) { fx[n]=fy[n]=fz[n]=0.0; }

   // do somthing before the loop (a trick to make the code faster)
   double cutoff_square = cutoff * cutoff;
   double lxh = lx * 0.5;
   double lyh = ly * 0.5;
   double lzh = lz * 0.5;

   for (int i = 0; i < N; ++i)
   {
       for (int k = 0; k < NN[i]; k++)
       {
           int j = NL[i * MN + k];
           if (j < i) { continue; } // will use Newton's 3rd law

           double x_ij = x[j] - x[i];
           double y_ij = y[j] - y[i];
           double z_ij = z[j] - z[i];
           apply_mic(lx, ly, lz, lxh, lyh, lzh, &x_ij, &y_ij, &z_ij);

           double r_2 = x_ij * x_ij + y_ij * y_ij + z_ij * z_ij;
           if (r_2 > cutoff_square) { continue; }

           double r_4  = r_2 * r_2;
           double r_8  = r_4 * r_4;
           double r_14 = r_2 * r_4 * r_8;
           double f_ij = factor_1 / r_8 - factor_2 / r_14;
           fx[i] += f_ij * x_ij; fx[j] -= f_ij * x_ij; // use Newton's 3rd law
           fy[i] += f_ij * y_ij; fy[j] -= f_ij * y_ij;
           fz[i] += f_ij * z_ij; fz[j] -= f_ij * z_ij;  
       }
   }
}


// velocity-Verlet method of integration
void integrate
(
   int N, double time_step, double *m, double *fx, double *fy, double *fz,  
   double *vx, double *vy, double *vz, double *x, double *y, double *z,  
   double *x_msd, double *y_msd, double *z_msd, int flag
)
{
   double time_step_half = time_step * 0.5;
   for (int n = 0; n < N; ++n)
   {
       double mass_inv = 1.0 / m[n];
       double ax = fx[n] * mass_inv;
       double ay = fy[n] * mass_inv;
       double az = fz[n] * mass_inv;
       vx[n] += ax * time_step_half;
       vy[n] += ay * time_step_half;
       vz[n] += az * time_step_half;
       if (flag == 1)  
       {  
           x[n] += vx[n] * time_step;  
           y[n] += vy[n] * time_step;  
           z[n] += vz[n] * time_step;  
           x_msd[n] += vx[n] * time_step;  
           y_msd[n] += vy[n] * time_step;  
           z_msd[n] += vz[n] * time_step;  
       }
   }
}


// calculate MSD from position data
void find_msd
(int N, int Nd, int Nc, double delta_tau, double *x, double *y, double *z)
{  
   FILE *fid_out = fopen("msd.txt", "w");
   int M = Nd - Nc;
   for (int nc = 0; nc < Nc; nc++)
   {
       double msd_x = 0.0; double msd_y = 0.0; double msd_z = 0.0;
       for (int m = 0; m < M; m++)
       {
           double sum_x = 0.0; double sum_y = 0.0; double sum_z = 0.0;
           for (int i = 0; i < N; i++)
           {
               double dx = (x[m * N + i] - x[(nc + 1 + m) * N + i]);
               double dy = (y[m * N + i] - y[(nc + 1 + m) * N + i]);
               double dz = (z[m * N + i] - z[(nc + 1 + m) * N + i]);
               sum_x += dx * dx; sum_y += dy * dy; sum_z += dz * dz;
           }
           msd_x += sum_x; msd_y += sum_y; msd_z += sum_z;
       }
       msd_x /= N * M; msd_y /= N * M; msd_z /= N * M;
       fprintf
       (
           fid_out, "%25.15e%25.15e%25.15e%25.15en",  
           nc * delta_tau, msd_x, msd_y, msd_z
       );
   }
   fclose(fid_out);
}


// calculate VAC from velocity data
void find_vac
(int N, int Nd, int Nc, double delta_tau, double *vx, double *vy, double *vz)
{    
   FILE *fid_out = fopen("vac.txt", "w");
   int M = Nd - Nc;
   for (int nc = 0; nc < Nc; nc++)
   {
       double vac_x = 0.0; double vac_y = 0.0; double vac_z = 0.0;
       for (int m = 0; m < M; m++)
       {
           for (int i = 0; i < N; i++)
           {
               vac_x += vx[m * N + i] * vx[(m + nc) * N + i];  
               vac_y += vy[m * N + i] * vy[(m + nc) * N + i];
               vac_z += vz[m * N + i] * vz[(m + nc) * N + i];
           }
       }
       vac_x /= N * M; vac_y /= N * M; vac_z /= N * M;
       fprintf
       (
           fid_out, "%25.15e%25.15e%25.15e%25.15en",  
           nc * delta_tau, vac_x, vac_y, vac_z
       );
   }
   fclose(fid_out);
}


// the main function
int main(int argc, char *argv[])
{
   srand(time(NULL));

   int nx = 4; // number of unit cells in the x-direction
   int ny = 4; // number of unit cells in the y-direction
   int nz = 4; // number of unit cells in the z-direction
   int n0 = 4; // number of particles in the unit cell
   int N = n0 * nx * ny * nz; // total number of particles

   int Ne = 100000;  // number of steps in the equilibration stage
   int Np = 100000;  // number of steps in the production stage
   int Ns = 10;      // sampling interval
   int Nd = Np / Ns; // number of position/velocity data
   int Nc = Nd / 10; // number of correlation data
   double time_step = 10.0 / TIME_UNIT_CONVERSION; // time step
   
   double T_0 = 1.475 * 119.8; // temperature prescribed
   double ax = 11.6449;        // lattice constant in the x direction
   double ay = ax;             // lattice constant in the y direction
   double az = ax;             // lattice constant in the z direction
   double lx = ax * nx;        // box length in the x direction
   double ly = ay * ny;        // box length in the y direction
   double lz = az * nz;        // box length in the z direction

   double rcn = 15.0;     // cutoff distance for neighbor list
   double rcf = 10.0;     // cutoff distance for force
   int neighbor_update_interval = 10; // interval of neighbor list update
   int MN = 100;                      // amximal number of neighbors  
   
   // memory for neighbor list
   int *NN = (int*) malloc(N * sizeof(int));
   int *NL = (int*) malloc(N * MN * sizeof(int));

   // major data for the particles
   double *m  = (double*) malloc(N * sizeof(double)); // mass
   double *x  = (double*) malloc(N * sizeof(double)); // position
   double *y  = (double*) malloc(N * sizeof(double));
   double *z  = (double*) malloc(N * sizeof(double));
   double *vx = (double*) malloc(N * sizeof(double)); // velocity
   double *vy = (double*) malloc(N * sizeof(double));
   double *vz = (double*) malloc(N * sizeof(double));
   double *fx = (double*) malloc(N * sizeof(double)); // force
   double *fy = (double*) malloc(N * sizeof(double));
   double *fz = (double*) malloc(N * sizeof(double));

   // copy of position data for msd calculations (do not apply PBC to them)
   double *x_msd = (double*) malloc(N * sizeof(double));  
   double *y_msd = (double*) malloc(N * sizeof(double));  
   double *z_msd = (double*) malloc(N * sizeof(double));  

   // data used for correlation function calculations
   double *x_all = (double*) malloc(Nd * N * sizeof(double));  
   double *y_all = (double*) malloc(Nd * N * sizeof(double));  
   double *z_all = (double*) malloc(Nd * N * sizeof(double));  

   double *vx_all = (double*) malloc(Nd * N * sizeof(double));  
   double *vy_all = (double*) malloc(Nd * N * sizeof(double));  
   double *vz_all = (double*) malloc(Nd * N * sizeof(double));  
   
   // initialize mass, position, and velocity
   for (int n = 0; n < N; ++n) { m[n] = 40.0; } // mass for argon atom
   initialize_position(n0, nx, ny, nz, ax, ay, az, x, y, z);
   for (int n = 0; n < N; n++) // make a copy
   {  
       x_msd[n] = x[n];  
       y_msd[n] = y[n];
       z_msd[n] = z[n];  
   }
   initialize_velocity(N, T_0, m, vx, vy, vz);

   // initialize neighbor list and force
   find_neighbor(N, NN, NL, x, y, z, lx, ly, lz, MN, rcn);
   find_force(N, NN, NL, MN, lx, ly, lz, x, y, z, fx, fy, fz, rcf);
 
   // equilibration
   clock_t time_begin = clock();
   for (int step = 0; step < Ne; ++step)
   {  
       if (0 == step % neighbor_update_interval)
       {
           find_neighbor(N, NN, NL, x, y, z, lx, ly, lz, MN, rcn);
       }

       integrate
       (
           N, time_step, m, fx, fy, fz, vx, vy, vz, x, y, z,  
           x_msd, y_msd, z_msd, 1
       );

       find_force(N, NN, NL, MN, lx, ly, lz, x, y, z, fx, fy, fz, rcf);

       integrate
       (
           N, time_step, m, fx, fy, fz, vx, vy, vz, x, y, z,  
           x_msd, y_msd, z_msd, 2
       );

       scale_velocity(N, T_0, m, vx, vy, vz); // control temperature

       apply_pbc(N, lx, ly, lz, x, y, z); // needed for simulating fluids
   }  
   clock_t time_finish = clock();
   double time_used = (time_finish - time_begin) / (double) CLOCKS_PER_SEC;
   fprintf(stderr, "time use for equilibration = %f sn", time_used);  

   // production
   time_begin = clock();
   for (int step = 0; step < Np; ++step)
   {  
       if (0 == step % neighbor_update_interval)
       {
           find_neighbor(N, NN, NL, x, y, z, lx, ly, lz, MN, rcn);
       }

       integrate
       (
           N, time_step, m, fx, fy, fz, vx, vy, vz, x, y, z,  
           x_msd, y_msd, z_msd, 1
       );

       find_force(N, NN, NL, MN, lx, ly, lz, x, y, z, fx, fy, fz, rcf);

       integrate
       (
           N, time_step, m, fx, fy, fz, vx, vy, vz, x, y, z,  
           x_msd, y_msd, z_msd, 2
       );

       apply_pbc(N, lx, ly, lz, x, y, z); // needed for simulating fluids

       // record data
       if (0 == step % Ns)  
       {
           int offset = (step / Ns) * N;
           for (int n = 0; n < N; n++)
           {  
               x_all[n + offset]  = x_msd[n]; // record x_msd, not x
               y_all[n + offset]  = y_msd[n];
               z_all[n + offset]  = z_msd[n];
               vx_all[n + offset] = vx[n];
               vy_all[n + offset] = vy[n];
               vz_all[n + offset] = vz[n];
           }
       }
   }  
   time_finish = clock();
   time_used = (time_finish - time_begin) / (double) CLOCKS_PER_SEC;
   fprintf(stderr, "time use for production = %f sn", time_used);  

   // free some memory
   free(NN); free(NL); free(m);  free(x);  free(y);  free(z);
   free(vx); free(vy); free(vz); free(fx); free(fy); free(fz);
   free(x_msd);  free(y_msd);  free(z_msd);

   // calculate MSD and VAC
   find_msd(N, Nd, Nc, Ns * time_step, x_all, y_all, z_all);
   find_vac(N, Nd, Nc, Ns * time_step, vx_all, vy_all, vz_all);

   // free some other memory
   free(x_all);  free(y_all);  free(z_all);
   free(vx_all); free(vy_all); free(vz_all);

   return 0;
}

六、matlab后处理脚本

编译并运行(在我的笔记本上运行时间不到半分钟)上面的C语言程序,将会在当前目录产生两个文件,msd.txt和vac.txt,分别记录了均方位移数据和速度自关联函数数据。运行我给的matlab脚本,即可通过爱因斯坦公式和格林-久保公式分别计算一个跑动扩散系数,然后画图显示结果(注意脚本里面的单位换算)。matlab脚本如下:


clear all; close all;
load msd.txt;
load vac.txt;

N = size(msd,1);
t = msd(:,1) * 10.18 / 1000; %ps
time_step = t(2) - t(1);
msd = [0; mean(msd(:,2:4),2) / 100]; % nm^2
vac = mean(vac(:,2:4),2) * (1.6/1.66*100); % nm^2/ps^2

D_Einstein = zeros(N, 1);
D_Green_Kubo = zeros(N, 1);

for n = 1 : N
   D_Einstein(n) = (msd(n + 1) - msd(n)) / time_step / 2;
   D_Green_Kubo(n) = sum(vac(1:n)) * time_step;
end
D_Green_Kubo = D_Green_Kubo - 0.5 * time_step * vac(1);

figure(1);
plot(t+time_step, msd(2:end), 'b', 'linewidth', 3);
xlabel('correlation time (ps)', 'fontsize', 20);
ylabel('MSD (nm^2)', 'fontsize', 20);
set(gca, 'fontsize', 20)

figure(2);
plot(t, vac, 'r', 'linewidth', 3);
hold off;
xlabel('correlation time (ps)', 'fontsize', 20);
ylabel('VAC (nm^2/ps^2)', 'fontsize', 20);
set(gca, 'fontsize', 20)

figure(3);
plot(t + time_step, D_Einstein, 'b-', 'linewidth', 4);
hold on;
plot(t + time_step, D_Green_Kubo, 'r--', 'linewidth', 3);
hold off;
ylim([0,0.15]);
xlabel('correlation time (ps)', 'fontsize', 20);
ylabel('D (nm^2/ps)', 'fontsize', 20);
legend('Einstein', 'Green-Kubo');
set(gca, 'fontsize', 20)

运行该脚本将得到如下三幅图:




上图显示了均方位移随关联时间的变化趋势。可以看出,关联时间比较大时,均方位移曲线变得很直,这代表正常扩散。关联时间较小时(几个皮秒内),均方位移从关联时间的二次函数过渡到一次函数,代表粒子的扩散从无障碍的弹道(ballistic)行为过渡到稳定的扩散行为。




上图显示了速度自关联函数随关联时间的变化趋势。可以看出,关联时间比较大时,速度自关联函数衰减为零,这也代表正常扩散。



上图显示了两条跑动扩散系数:一条是对均方位移求导(用爱因斯坦公式)得到的(蓝色的实线);一条是对速度自关联函数求积分(用格林-久保公式)得到的(红色的虚线)。可以看出,它们几乎是重合的。理论上,它们应该完全重合,但由于数值积分和数值求导都有计算误差,它们还是有一点点差别。但是,这个差别还是比统计误差要小。这里说的统计误差指的是用不同的初始速度得到的不同结果之间的差别。


七、与前人结果的比较

      我忘了我模拟的状态(温度和密度)是参考哪个文献的了(好像是Haile的书里面的;这本书的信息请参考我前一篇博文)。下次找到再说。不过我记得我与文献中的结果比较过,是一致的。


八、小结

     本文介绍了在分子动力学模拟中用爱因斯坦公式和格林-久保公式求(自)扩散系数的方法并给出了完整的C语言程序以及相应的matlab后处理脚本。本文讨论的是最简单的扩散现象,属于教科书水平的问题。文献中有不少研究“反常”扩散现象的,属于比较高级与活跃的话题。我对此暂无研究,故不再进一步讨论。感兴趣的读者可深入研究之。






http://blog.sciencenet.cn/blog-3102863-991891.html

上一篇:A C code for thermal conductivity calculations in EMD
下一篇:CUDA编程与应用之(一)——在Ubuntu下安装CUDA工具箱
收藏 分享 举报

1 王福昌

该博文允许注册用户评论 请点击登录 评论 (12 个评论)

数据加载中...

Archiver|手机版|小黑屋|科学网 ( 京ICP备14006957 )

GMT+8, 2017-9-20 00:58

Powered by ScienceNet.cn

Copyright © 2007-2017 中国科学报社

返回顶部