Transcript Document

Molecular dynamics study of gas phase and gas-surface
reaction using “MD Trajectory” software complex
Michael Pogosbekian
Valery Kovalev
Institute of Mechanics, Moscow State University
Dept. Mechanics and Mathematics, Moscow State University
MD Trajectory
Three levels of simulation

N  X


X

N

i
i
i 1 i
i 1 i
AVOGADRO structure
N
Reaction rate : w  k   [ X i ]
i 1
Equilibrium chemical kinetics
One-temperature approach
Experimental DataBase
k  T   A  T n exp  Ea / kT 
Non-equilibrium chemical kinetics
Two-temperature approach
Level chemical kinetics
where
 i
Theoretical models
k  T ,Tv 
kv  T 
kvw  T 
Computarized handbook
Models’ catalogue
“MD Trajectory”
T,Tv – translational and vibrational temperatures
v, w – vibrational states of reagents and products
MD Trajectory
Two typical cases of the nonequilibrium conditions
Boeing winged version of the Orbital
Space Plane during reentry (T >> Tv)
Gas discharge (T << Tv)
MD Trajectory
Software complex “MD Trajectory”
Triatomic collisions
A  BC( v )  A  BC( w )
A  BC( v )  XY ( w )  P
A  BC( v )  A  B  C
vibrational relaxation
exchange reaction
dissociation
where XY  {AB, AC} and P - residual atom
Tetratomic collisions
AB( v )  CD( v'
AB( v )  CD( v'
AB( v )  CD( v'
AB( v )  CD( v'
)  AB( w )  CD( w' )
)  XY ( w )  WZ ( w' )
)  MN( w )  P  Q
)  A  B C  D
vibrational relaxation
exchange reaction
partial dissociation
full dissociation
where XY and WZ - diatomic molecule  {AB, CD}, MN - any
diatomic molecule and P,Q - residual atoms
classical trajectory method
Collision of A atom with BC molecule
Jacobi coordinate
qi
 xi 6  xi 3
i  1,2 ,3
relative B-C motion
mB xi 3  mC xi 6
mB  mC
i  1,2 ,3
relative A-BC motion
mA xi  mB xi 3  mC xi 6
M
i  1,2 ,3
ABC motion as whole
qi 3  xi 
qi 6 
where
M  mA  mB  mC
3 

mC qi
R1   
 qi  3 

i 1  mB  mC
2
R1  RAB
R2 
3

i 1
qi2
R2  RBC
3 
R3  RAC

mB qi
R3   
 qi  3 

i 1  mB  mC
2
classical trajectory method
Collision of A atom with BC molecule
Motion equations
d qi  H

dt
 pi
i  1, 2, ...6
3 U  R
d pi
H
U
k




dt
 qi
 qi
k 1  Rk  qi
i  1, 2, ...6
d q i 6 1
 pi 6
dt
M
i  1,2 ,3
d pi 6
0
dt
i  1,2 ,3
3 

1
1
1
 pi2 
 pi2 3 
 pi26   U R1 , R2 , R3 
where H   
2  A BC
2M

i 1  2 BC
BC 
mB mC
mB  mC
 A BC 
mA mB  mC 
M
classical trajectory method
Definition of initial conditions
Collision’s scheme
Modified parameters
b 
 
 
 
 
v
j
Et
0, bmax 
0,  
0,2  
0,2 
0,2 
MD Trajectory
Numerical integration scheme of motion equations
Kutta-Merson method
y  f t , y 
yn
y n 1
tn
tn  
k1n  f tn , y n 

 

kn2  f  tn  , y n  k1n 
3
3 



 

kn3  f  tn  , y n  k1n  kn2 
3
6
6 



3 

kn4  f  tn  , y n  k1n  kn3 
2
8
8



3


kn5  f  t n  , y n  k1n  kn3  2 kn4 
2
2


forth order approximation
automatic selection of integration step
(reduce calculation time for 10-30%)
2
1 
1
y n 1     k1n  kn4  kn5   y n
3
6 
6
1 1 3 3
4
y nerror



k

k

2
k

n
n
n   yn
1
2
2


error  y n 1  y nerror
1
error > max :  =  / 2
error < min :  = 1.5 
max  error  min :  = 
MD Trajectory
Random number generator uniformly distributed in (0,1)
Main generator
X ( N  1) 
 X ( N )  517  mod 240
Secondary generator


Y ( N  1)   Y ( N )  23 mod 108  1
M.D. MacLaren – G.Marsaglia method

V [ 0 ], V [ 1], , V [ k  1]  X / 240


j  kY / 108  1


out  V [ j ]; V [ j ]  X / 240


classical trajectory method
Results of trajectory calculations
Reaction cross-section
2
  Et ,v , j    bmax
P  Et ,v , j 
P  Et ,v , j  
1
2
 2 3 bmax
bmax

2
2
2
0
0
0
0
0
 b d b  sin  d   d  d  Pr Et ,v , j , b, , ,,  d
Monte-Carlo method
P  Et ,v , j  


lim  Pr Et ,v , j ; bi  , i  , i  , i  , i   lim
N
N  i 1
   Et ,v , j     Et ,v , j 
N  Nr
NNr
N 
Nr  Et ,v , j 
N  Et ,v , j 
classical trajectory method
Results of trajectory calculations
Level rate constants
1/ 2
 8 kT 

kv , j  T   
   ABC 
2
 1 
  Et 



E
,
v
,
j
E
exp

 

 dEt
t
t
 kT  0
 kT 
Two-temperature rate constant


exp Ev kTv  g j  2 j  1 exp  E j kT
k  T ,Tv   
kv , j T 

Q
Q
v
r
v
j
PES
Potential energy surface
Semiempirical methods
Generalized LEPS model
Method of diatomic complexes in molecules
Bond energy- bond order method
Ab-initio calculations
GAUSSIAN
MOLCAS
GAMESS (special version for INTEL platform - PC GAMESS)
MD Trajectory
Generalized LEPS model

3
1 3  J i
Jj
Qi

U 
 


1

c
2
1

c
i  i 1 
i 1cj
i 1
 j  i
Qi  J i
 Di
1  ci
Qi  Ji
D
 i
1  ci
2
1/ 2
2
 
 

 

 exp   2i  Ri  Ri0    2 exp   i  Ri  Ri0   
 exp   2 i  Ri  Ri0    2 exp   i  Ri  Ri0   
Di - dissociation energy
 i - Morse parameter
Ri0 - equilibrium distance
ci - adjusted Sato parameter
MD Trajectory
Analytical representation of PES
Many-body expansion method
 VA1  VB1  VC1 
VABC R1, R2 , R3 
2 R   V 2 R   V 2  R   V 3  R , R , R 
VAB
1
BC 2
AC 3
ABC 1 2 3
One-particle term
VA1,VB1,VC1
– define the energy of
electronically excited atom
Two-particles term
2 ,V 2 ,V 2  – describe the potential curve of
VAB
BC AC
diatomic molecule
Three-particles term
3 
VABC
– define the interaction
at close internuclear distances
MD Trajectory
Sorbie-Murrell function
Two-particles term
M



2  
V
R  D  1   ak  k   e a1 

k 1
where   R  Re

– extended Rydberg function
– displacement from equilibrium distance
Three-particles term
V 3 R1, R2 , R3   P  T
3
P  C0   Ci i   Cij i  j   Cijk i  j k 
i 1
where
i j
i , j , k, l ,... 1,2,3
    
T   1  tanh i i  
 2 
i 1 
3 
i  j k

i  j k l
– bond number
– switching function
Cijkl i  j k l  ...
MD Trajectory
Garcia-Lagana (bond-order) function
Two-particles term
N

2
V R   D   a j n j
j 1

b R Re
where n  e

– polynomial of N-th order
– bond order coordinate
Three-particles term
M

3
V R1, R2 , R3    c jkl n1j n2k n3l
j ,k ,l
where j  k  l  j  k  l
j k l M
– polynomial of M-th order
MD Trajectory
Aguado-Paniagua function
Two-particles term
c e
V 2 R   0
 R
R
where   R  e
N
  ck  k
k 1
  l R
l = 2 or 3 for two- or three-particles term
Three-particles term
M

3
V R1, R2 , R3    dijk 1i 2j 3k
i , j ,k
where i  j  k  i  j  k
i  j k M
– polynomial of M-th order
MD Trajectory
Features of software complex
Wide set of PES analytical functions
State-to-state rate constants
kv w T 
v, w - vibrational levels of reagent and product
Angle distribution of reaction products
Distribution products by vibrational and rotational numbers
Save coordinate and pulses along of trajectory for subsequent
demonstration purposes
Optimization of trajectory code for usage of cluster technologies
based on MPI (Message Passing Interface)
Object oriented C++ code
XML – style for input & output data (LibXML library)
MD Trajectory
High performance supercomputer facilities
Moscow State University, cluster SCI - 36 CPUs total
Node configuration: Dual Pentium III/500MHz, 1Gb RAM, 3.2 Gb HDD
Network environment : SCI + Fast Ethernet
Russian Academy of Sciences, cluster MVS-1000M - 768 CPUs total
Node configuration: Dual Alpha 21264A/667MHz, 1Gb RAM, 15 Gb HDD
Network environment : Myrinet (2 Gb/s) + Fast Ethernet
MD Trajectory
Parallel version of “MD Trajectory”
Service layer - 1 muster process
Module of
task definition
Module of
data allocation
Module of
results visualization
Module of
results processing
Calculation layer - N slave processes
Trajectories calculation
module #1
Trajectories calculation
module #N
Data access layer
XML – files (LibXML2 library)
SQL Server (MySQL is planned)
PES, Tasks, Results
MD Trajectory
Data Access Layer
Three kinds of input & output files
Input data – fixed size – rather simple structure – less than 1 Mb
Control Group describes control function;
Molecule Group contains spectroscopic characteristic of diatomic
molecules
PES Group
describes PES of the investigated system
Log file – dynamic size – rather simple structure – less than 10 Mb
Auxiliary information which can be required for calculation control
and calculation continuation at the next time
Result file – dynamic size – very complex structure – up to 100 Mb
Two ways for realization of Data Access Layer
XML files. It is realized due to LibXML2 library
XML parser and toolkit of Gnome
WWW . XMLSOFT . ORG
DataBase connectivity module for MySQL Server is planned
MD Trajectory
MSU cluster
MD Trajectory
MSU cluster
CO + N  CN + O
Leading reaction in CN formation behind the shock wave front
in Martian atmosphere
A single experimental work for "short" temperature range
Absence of data at high temperatures
Experiment, L.B. Ibragimova
CO + N  CN + O
Potential energy surface
2

  
   
A" PES for CO 1   N 2D  NO  2  C 3P
Modified LEPS model [1]
Sorbie-Murrell function [3]
based on ab-initio data [2,4]
Aguado-Paniagua function [4]
4

  

  
A" PES for CO 1  N 4S  CN 2  O 3P
Generalized LEPS model [1]
References
1. K.J.Schmatjko and J.Wolfrum, Ber. Bunsen Phys. Chem., 1975, 79, pp.696-707
2. P.Halvick, J.C.Rayez, E.M.Evleth, J. Chem. Phys., 1984, 81, pp.728-737
3. SM.Simonson, N.Markovic, S.Nordholm and B.J.Persson, Chem. Phys., 1995,
200, pp.141-160
4. Andersson, N.Markovic and G.Nyman, Phys. Chem. Chem. Phys., 2000, 2,
pp.613-620
CO + N  CN + O
Generalized LEPS model
Equipotential contour map
CO + N  CN + O
Generalized LEPS model
3D View
CO + N  CN + O
Reaction cross-sections for CO(v,j)
CO + N  CN + O
Comparison QCT calculations with experimental data
One-temperature rate constant
CO + N  CN + O
Comparison QCT calculations with experimental data
One-temperature rate constant
CO + N  CN + O
Two-temperature rate constants
CO + N  CN + O
Nonequilibrium factor Z T ,Tv  
k T ,Tv 
k T ,Tv  T 
CO + N  CN + O
Level rate constants
MD Trajectory
Theoretical models for exchange reactions
 - model

 Ea   Ev 
 for  Ev  Ea
k0 T   exp  
kT
kv T   


k T 
for  Ev  Ea
 0
Generalized Marrone-Treanor model (CVCV)

 D Ev


k
T

exp

 0
kU

kv T   
k T   exp   D Ev
 0
kU


 Ea Ev 
exp
for Ev   Ea



kT 



 1   Ea 
 exp  
 for Ev   Ea
kT



Theoretically informational model
CO + N  CN + O
Comparison QCT calculations with theoretical models
Level factor
CO + N  CN + O
Comparison QCT calculations with theoretical models
Level factor
CO + N  CN + O
Comparison QCT calculations with theoretical models
Level factor
CO + N  CN + O
Comparison QCT calculations with theoretical models
Level factor
MD Trajectory
Further development of “MD Trajectory”
Investigation of Gas-Surface processes
Design of thermal protection systems in space vehicles
Microelectronics applications
Heterogenous combustion
Recombination processes
Eley-Rideal
Langmuir-Hinshelwood
Main objectives
Recombination coefficient

Accomodation coefficient of chemical energy

Cite-specific effects and influence of top-layer surface structure
MD Trajectory
Classical molecular dynamics
Atoms are divided in two groups:
1. i = 1, … n (gas-phase atoms)
2. k = 1, … N (lattice atoms)
Total hamiltonian is:
1 2 N
1 2
H  
Pi   
Pk  V11 Rij  V22 Rkl   V12 Rik 
i 1   x , y , z M i
k 1   x , y , z M k
i j
k l
ik
n
The hamiltonian equations of motion are:
i  Pi M i ;
Rij

H

V
V12 Rik
11

Pi  
 

 i
j i Rij  i
k Rik  i
k  Pk M k ;
H
V R
V R
Pk  
  22 kl   12 ik
 k
l  k Rkl  k
k Rik  k
(  x, y, z )
MD Trajectory
Definition of initial conditions
Collision scheme
Assumptions
flat surface
instead rough one
monocrystal
instead polycrystal
clear surface
without adsorbed
layer
Detailed description of classical molecular dynamics is represented in:
Gert D. Billing Dynamics of Molecular Surface Interactions. New York,
John Wiley&Sons, 2000, chapter 6, pp.93-102
MD Trajectory
Definition of initial conditions
For incident gas atom B:
X B  R sin cos   X Bs
PX  PR sin  cos 
YB  R sin  sin   YBs
PY  PR sin  sin 
ZB  R cos
PZ  PR cos 
where
PR  2MBEcoll
X Bs ,YBs - randomly distributed on the surface
For adsorbed gas atom A:
X As ,YAs - randomly distributed on the surface
PX , PY , PZ
- the same as for atom B, where
For lattice atoms:
 
Rk  Rk0  2k BTs / Fk cos  k0
0
where Rk - equilibrium position
PR  2MBkBTs
 
Pk  2M k k BTs / Fk sin  k0
Ts - surface temperature
Fk
- force constant for atom k
 k0
- phase angle, randomly distributed in 0,2 
Oad + Ogas  O2
PES for -cristobalite
B.P.Feuston, S.H.Garofalini Empirical three-body potential for vitreous silica.
Journal of Chemical Physics, Vol. 89, No. 9, 1988. pp. 5818-5824


V r1, r2 , r3 ,..., rN    v 1ri    v 2 ri , r j 
i j
i
 v 3 ri , r j , rk   ...  v N r1, r2 ,..., rN 
i  j k
Modified form of the Born-Mayer-Huggins (BMH) potential:


 



 

v 2 ri , r j  v 2 rij  Aij exp  rij /   Z i Z j e 2 / rij erfc rij /  ij
Aij  1  Z i / ni  Z j / n j  b  exp  i   j / 
where rij - separation distance, Z i - formal ionic charge, b – constant,



 
 ,  ij - adjustable parameters, ni - number of valence shell electrons,

 
 
 
v 3 ri , r j , rk  h rij , rik , jik  h r jk , r ji , kji  h rki , rkj , ikj


 




where h rij , rik ,  jik  i exp  i / rij  r c   i / rik  r c cos  jik  cos  c
i
i
jik
c
c
for ( rij  ri and rik  ri ); h rij , rik , jik  0 in other case.


2
c
c
where i ,  i , ri , cos  jik - constants,  jik - angle subtended by rij and rik
Oad + Ogas  O2
Unit cell of  – cristobalite lattice
Ralph W.G. Wyckoff The crystal structure of the high temperature form of
cristobalite (SiO2), American Journal of Science, Ser.5, Vol.9, 1925, pp.448-459
Oad + Ogas  O2
Top layer structure of the  – cristobalite surface
M.Cacciatore, M.Rutigliano, G.D.Billing Eley-Rideal and Lengmuir-Hinshelwod
Recombination Coefficients for Oxygen on Silica Surfaces, Journal of
Thermophysics and Heat Transfers, Vol. 13, No. 2, 1999, pp.195-203.
Oad + Ogas  O2
Comparison of MD calculations results for SiO2
Eley-Rideal recombination probability
Oad + Ogas  O2
Comparison of MD calculations results for SiO2
Chemical energy accomodation coefficient in the Eley-Rideal recombination
Oad + Ogas  O2
MD calculations results for SiO2 surface
Vibrational distribution of the formed O2 molecules in the Eley-Rideal reaction
Oad + Ogas  O2
MD calculations results for SiO2 surface
Vibrational distribution of the formed O2 molecules in the Eley-Rideal reaction
Oad + Ogas  O2
Fragment of crystal lattice of 3C-SiC and top layer
structure of the surface
Oad + Ogas  O2
Comparison of MD calculations results
Eley-Rideal recombination probability
Oad + Ogas  O2
Comparison of MD calculations results
Chemical energy accomodation coefficient in the Eley-Rideal recombination
Oad + Ogas  O2
MD calculations results for SiC surface
Vibrational distribution of the formed O2 molecules in the Eley-Rideal reaction
Oad + Ogas  O2
MD calculations results for SiC surface
Vibrational distribution of the formed O2 molecules in the Eley-Rideal reaction
Thank you for your attention