Reactive Molecular Dynamics


A Molecular dynamics model with chemical reactions and electrochemistry.

Documentation (PDF) (HTML)
Presentation (ODP) (PPT)

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.

Bulk Reactions
Bulk reaction model
Surface Reactions
Surface reaction model

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.

1 cubic micron
H2 + O(s) = H2O reaction of 14 million molecules in 13µ


The main features of the code include:

  1. 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.
  2. Simulations of a gaseous phase and gas-solid surface reactions based on Collision Theory approximation.
  3. 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.
  4. Graphical output based on OpenGL, as well as bach job execution with data dumps and restart options.
  5. Parallel implementation (OpenMP).

Source Code and Documentation



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


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.

Valid XHTML 1.0 Strict