Reactive Molecular Dynamics
ReMoDy
A Molecular dynamics model with chemical reactions and electrochemistry.
The Concept
A common drawback of molecular dynamics (MD) simulations is their inability to describe bond formation,
which is essential in modeling of chemical reactions
This drawback can be
remedied if an MD simulation could incorporate not only the intermolecular force potentials, but also
allow for the possibility of different states of molecular excitation and bonding. These effects can in
principle be modeled in a probabilistic sense, if the MD simulation takes advantage of the relevant
probability data from ab-initio quantum mechanical simulations or experimental results.
In this study a molecular dynamics code was developed for computing reactive gas mixtures, including
catalytic reactions at active surfaces (electrochemistry) based on the
Collision Theory.
In simulations of a
gaseous phase and gas-solid interactions, each molecule is advanced following the classical laws of
conservation of its linear and angular momenta. At the same time a reaction possibility is considered,
which determines the transition event for the two species to enter into a reaction with the
production of new species or forming a molecular bond. The reaction event occurs when the impact energy
of the molecules plus their internal energy exceeds the known activation energy for the reaction. Each
molecule of a newly formed specie is again treated as a rigid body subjected to Newton's laws of motion.
Energy distribution among the molecules which result from the collision are calculated in
accordance with the molecular degrees of freedom determined by molecule's specific heat constants.
To be able to trace the motion and interactions of millions of molecules efficiently, a domain segmentation algorithm was implemented, which enabled to reduce the time of processing interactions from N2 to a near-linear dependence. The code was implemented in C++ programming language. The figure below shows a snapshot of a simulation of hydrogen interaction with the reactive surface with the production of H2O. The simulation inside one micron pore under atmospheric conditions performed on a 2GHz CPU, 4GB RAM workstation took about one week. The program enables one to simulate up to 10 million molecules per each Gigabyte of RAM.
Features
The main features of the code include:
- Specification of species and reactions data. These data are supplied in an XML file
and contain the necessary infomation for molecular dynamics simulations and chemical reactions, similar
to CHEMKIN input.
- Simulations of a gaseous phase and gas-solid
surface reactions based on
Collision Theory approximation.
- Specification of physical and chemical boundary conditions, such as: elastic, sticky, adiabatic, isothermal, as well as open boundary with the possibility to assign gas compositions, temperatures and pressures across the open boundaries.
- Graphical output based on OpenGL, as well as bach job execution with data dumps and restart options.
- Parallel implementation (OpenMP).
Source Code and Documentation
Utilities
Pre-Processing
The init.cc utility can be used to generate the initial uniform distribution of molecules. For
example, the command:
./init init
will generate file init.dat.gz which can be used as input to remody, like:
./job init.dat.gz
To se the parameters one should edit the init.cc file and compile the program as:
g++ -lz init.cc -o init
Post-Processing
Script readlog.py is used to parse the output file
generated by the remody job, when run as:
./job -f remody.xml init.dat.gz > job.log &
The above command will start the job process, which will dump screen output into the job.log file.
Note that the 'xterm' tag in the config file remody.xml should be set to 1 as
<xterm>1</xterm> to generate the log file in the appropriate format.
To extract the sequence of molecular concentrations of, say, hydrogen (H2) for the whole domain at different times, one should run the following script:
python readlog.py H2 < job.log > H2.dat
File H2.dat will contain time distribution of concentration of H2.
Program hist.cc is used to parse the dump file generated by
remody when run as
./job -f remody.xml job-02.dat.gz &
where job-02.dat.gz is the dump file from the previous run used to restart the run.
Usually the run will generate the output dump files at regular time intervals as specified in the
remody.xml file inside the <time><output>...<time><output> tag.
The output files will have names as: job-XX.dat.gz where XX is the next consecutive number of
the output dump. To retrieve a histogram of molecular distributions inside the domain, one can run the
hist utility as:
gunzip -c job-XX.dat.gz | grep ^H2 | grep -v ^H2O | cut -f 3 | hist > H2.dat
This will produce file H2.dat with the column of numbers corresponding to the number of molecules for
different slices of the domain. Note, that the repeated grep functions were used to capture only
the lines beginning with H2 and exclude those beginning with H2O. Other more efficient line retrieval
can be arranged using awk, perl, or python.
To change the domain limits, slice width, the number of slices, and the direction
of slicing one should edit the hist.cc file accordingly and recompile it
as:
g++ -lm hist.cc -o hist
Licencing and Distribution
ReMoDy code was developed at NIFT under the DOE EPSCOR program.