분자 운동 모사를 이용한 상호 확산 계수의 계산

화학공학의 이론과 응용 제8권 제1호 2002년

Theories and Applications of Chem. Eng., 2002, Vol. 8, No. 1

전동혁, 이태용

한국과학기술원 생명화학공학과

Calculation of Mutual Diffusion Coefficients Using Molecular Dynamics Simulation

Dong-Hyuk Chun and Tai-yong Lee

Department of Chemical and Biomolecular Engineering, KAIST

Introduction

Molecular dynamics simulation can be used to compute physico-chemical quantities of practical interest. Kinetic energy, potential energy and total energy can be obtained by averaging the expressions for these quantities over the trajectory. Temperature and pressure can be obtained from the equipartition theorem. Fluctuations of various quantities give the expansion coefficients, specific heats, and the compressibilities.[1]

In the present study, diffusion coefficient will be calculated among many physical properties, because diffusion phenomena strongly depends on the molecular motion. There are several kinds of diffusion - self diffusion and mutual diffusion. Self diffusion is diffusion under single-particle system and mutual diffusion is diffusion under binary system.

Self diffusion coefficient is computed by molecular dynamics simulation from Einstein relation. Calculation of mutual diffusion coefficient is more difficult than that of self diffusion coefficient since mutual diffusion is involved with molecular interaction between different molecules. In this study, mutual diffusion coefficients are computed using molecular dynamics simulation.

Theory

1. Molecular Dynamics Simulation for Constant Temperature

A simple method of fixing the kinetic temperature of a system in molecular dynamics simulation is to rescale the velocities at each time step by a factor of [eqn: chi =LEFT ( T / IMAG RIGHT ) ^ 1/2] where [eqn: IMAG] is the current kinetic temperature, T is the desired thermodynamic temperature[1]. Constant kinetic temperature dynamics is generated by the equations of motion

[eqn: rm dot boldr =bold p /itm~~~~~~~~~]

[eqn: rm dot bold p =bold f - xi (bold r , bold p ) bold p ```.]

The quantity [eqn: rm xi ( bold r, bold p)] is a kind of friction coefficient which varies so as to constrain [eqn: IMAG] to a constant value. The constraint is chosen so as to perturb as little as possible the classical equations of motion, by making [eqn: xi] a Lagrange multiplier which minimizes the difference between the constrained and Newtonian trajectories. Finally a modified velocity equation in leap-frog scheme can be written as followed:

[eqn: rm dot boldr (it t + 1over2 delta t ) = (2 chi -1) rm dot boldr (it t- 1over2 delta t ) + chi rm bold f it (t) delta t / m] .

2. Simulation Technique

Generally, there are a huge number of molecule in a real system. It is extremely difficult to consider all molecules in the molecular dynamics simulation. Some simulation techniques are used in molecular dynamics simulation in order to reduce calculation time. In this study, periodic boundary condition and linked list method are applied and they are frequently used in ordinary molecular dynamics simulation.

1) Periodic Boundary Condition : In the course of the simulation, as a molecule moves in the original box, its periodic image in each of the neighboring boxes moves in exactly the same way. A two-dimensional version of such a periodic system is shown in Fig 1. When we simulate the system, we just store the coordinates of the molecules in the central box.[1]

Fig . A two dimensional periodic system

2) Cutoff Method : In this method, potential of two molecules is not considered whose distance is greater than a certain distance, so called cutoff distance. It can reduce calculation time. In this study, linked list is used for using cutoff method.[1]

3. Diffusion Coefficient

The self diffusion coefficient [eqn: D] is described (for three dimension)

[eqn: D= 1over3 int from 0 to inf < rm boldv (it t) cdot rm boldv (0) > it dt]

where [eqn: rm bold v ( it t)] is the center-of-mass velocity of a single molecule. It is matched with corresponding Einstein relation

[eqn: 2tD = 1 over 3 lim_t->inf <|rm boldr (it t) - rm boldr (0) |^2 >]

where [eqn: rm boldr (it t)] is the molecule position.

The mutual diffusion coefficient [eqn: D_12][6] is written as a product of the thermodynamic factor [eqn: Q] and the kinematic factor [eqn: L] :

[eqn: D_12 = QL] .

The thermodynamic factor [eqn: Q] is defined by

[eqn: Q=[1+x_1 x_2 (Gamma _11 + Gamma_22 -2Gamma_12 )]^-1]

with

[eqn: Gamma _ij = 4 pi int from 0 to inf r^2 [g_ij (r) -1 ] dr]

[eqn: g_ij] is the radial distribution function for the pair of species ij. The kinematic factor [eqn: L] is represented by Adamson relation[5],

[eqn: L=phi_2 D_1 + phi_1 D_2]

where [eqn: phi_i] is the volume fraction of component [eqn: i].

Model and method of molecular dynamics simulation

In present study we conduct simulation in the canonical ensemble at 300K. The molecules interact via Lennard-Jones potential[3],

[eqn: u(r) = 4epsilon _ij left[ left( sigma_ij over r right) ^12 - left( sigma_ij over r right) ^6 right ]] .

The Lennard-Jones potential parameters for the model cyclohexane and carbon tetrachloride molecules are given in Table 1[6].

Table . Lennard-Jones potential parameter

Molecule	[eqn: epsilon / k_B](K)	[eqn: sigma ( ANGSTROM )]
[eqn: rm C_6 H_12]	324	6.093
[eqn: rm CCl_4]	327	5.881

Molecules are confined in a cubic box and randomly placed in initial time. The particle velocity follows the Maxwell-Boltzmann distribution. Periodic boundary condition and linked list are employed. 365 molecules are traced with leap-frog algorithm. The integration time step is 0.0001ps.

Result and Discussion

The mean square displacements (MSD) of cyclohexane and carbon tetrachloride molecules are plotted in Fig 2 and Fig 3. The self diffusion coefficients were calculated from the slopes of the curves.  The calculation data of self diffusion coefficient and relative deviations with respect to experimental data are listed in Table 2.

Fig . MSD for cyclohexane and carbon tetrachloride

Fig 4 shows results for an ideal diffusion system cyclohexane-carbon tetrachloride. Ideal mutual diffusion coefficients are calculated by Darken's model[2],

[eqn: D_12 = x_2 D_1 + x_1 D_2]

where [eqn: x_i] is the mole fraction of component. Excess value of diffusion coefficients can be corrected by thermodynamic factor [eqn: Q], because, in the normal solution, [eqn: L] is not so different from [eqn: x_2 D_1 + x_1 D_2][4].

Fig 3. Mutual diffusion coefficients of cyclohexane-carbon tetrachloride. Solid line : this work; solid triangle : experimental data.

Table . Results of calculation and comparison

Molecule	Calculation Data	Experimental Data	Error(%)
[eqn: rm C_6 H_12]	1.5245	1.48	3%
[eqn: rm CCl_4]	1.1379	1.20	5%

Conclusion

The molecular dynamics simulation is a powerful method of physical property prediction and results from molecular dynamics simulation are very similar to experimental data. Therefore, we can predict many kinds of useful quantities using molecular dynamics simulation without experiment.

Acknowledgment

This work was also partially supported by the Brain Korea 21 project.

화학공학의 이론과 응용 제8권 제1호 2002년

Reference

[1] Allen M. P. and D. J. Tildesley, Computer Simulation of Liquids, Clarendon Press, Oxford 1988.

화학공학의 이론과 응용 제8권 제1호 2002년

[2] Dysthe D. K. and B. Hafskjold, Int. J. of Thermophysics 16(5) p. 1213 1995.

[3] Hail J. M., Molecular Dynamics Simulation, Wiley & Sons, New York 1992.

[4] Kataoka Yosuke, Fluid Phase Equilibria 144 p. 257 1998.

[5] Takaba H. and K. Mizukami, J. Membrane Science 121(2) p. 21 1996.

[6] Welty J. R., Charles E. Wicks and Robert E. Wilson, Foundamentals of Momentum, Heat, and Mass Transfer, Wiley & Sons, New York 1983.