I need to do a probability analysis for measuring the importance of the reaction pathways in a reactive system. Therefore, I have to perform a certain number of unique MD simulations. May I ask you if anyone can help me how I can do this in LAMMPS? I am using ReaxFF as forcefield. Should I start with unique conformations (for example, by changing the dihedral angles) or can I start do each of the unique simulations with identical conformations but different initial velocities? Any help is really appreciated.
Your question is somewhat vague, but you can do multiple
simulations easily in LAMMPS. It's up to you to prepare
your multiple input data files or input scripts. You can also
use variables in input scripts (see variable doc page for
world and universe variables), to have one script that runs
a multitude of simulations on different processors. Also
see the -partition command line switch.
Thanks a lot for your prompt response. With that question, I intended to ask for advice from expert users about tow aspects of statistical analysis of reaction pathways using "REAX" as the forefield, though my intention might not be conveyed properly in the last post. The first aspect is an appropriate method for doing this probability analysis and the second one is a proper implementation of this method into LAMMPS input script. I got a quite helpful guide for the latter aspect from you response.
it's the reaction pathways part of your Q that is vague to me. That can
mean a lot of things. Just running MD on a system doesn't produce
any direct info about a reaction pathway, so you'd have to describe
more directly what you want to calculate and how.
one thing that is not clear from your description is whether
you simply want to know which is the more likely reaction,
or whether you want to follow the reaction mechanism.
for the former, the free energies differences of the various
possible reaction scenarios is what you are looking for and
there is a wealth of literature on that subject describing a
variety of (mostly complementary) methods that can produce
the answer you are looking for.
brute force statistical analysis is typically the worst way, due
to the restricted statistical sampling in the simulations. some
kind of accelerated method is usually needed.
i suggest to search the literature for "transition path sampling"
if you are more interested into the mechanism itself.
the biggest challenge for this kind of calculation is to
find out whether your results are meaningful or not.
most of the time, it is very difficult to tell and requires
a lot of experience. i don't think that there is any approach
that comes with a guarantee to be successful.
Thanks a lot for your useful remarks. To be clear, my goal is to study the time evolution of a molecular system containing certain molecules (e.g. Oxygen and Methane molecules) which are allowed to form new bonds and to have bond cleavages due to the usage of REAXFF as a reactive potential. Since Reax is a bond-order potential, I would obtain the time history of reactions from the MD trajectory of the system after a sufficiently long time simulation. I need to obtain the reaction species and the kinetic rate constants of the reactions as a function of temperature.
But the critical issue which is frequently referenced in the literature is that these simulations involve rare reactive events for which multiple simulations starting from independent configurations are required to check the statistical reliability of the simulation results. As Chenoweth et. al have mentioned "One could perform multiple simulations and analyze the trajectories to get a statistical distribution of the observed reactions to provide a measure of the importance of each pathway." in their paper "ReaxFF Reactive Force Field for Molecular Dynamics Simulation of Hydrocarbon Oxidation". Other authors also make similar notes, but it is not clear that how many simulations are considered sufficient for collecting statistics to calculate the rate constants of reactions. Axel's suggestion for "transition path sampling" is good and I'm doing search for it.
So based on the description above, since I am interested in detecting most likely reactions and their respective rate constants, Axel and Steve do you believe that it is not so much fruitful to go through the procedure mentioned above to reach my goal and as Axel mentioned in the previous email I have to calculate the free energy differences of various reaction scenarios?
Another critical issue concerned with these kind of simulations is the time-scale problem associated with MD simulations. Because of the relatively short simulation time frame, high temperatures (in range of 2000-3000 K) are often used to accelerate the chemical reactions which would take a prohibitively long simulation time to happen at lower temperatures (below 1000 K). But how can we extrapolate the reaction rates calculated at high temperatures to obtain the rates for lower temperatures while those reaction mechanisms can just be activated at high temperatures? In analogous situation, we cannot extrapolate the highly random fluid characteristics in turbulent regime to laminar regime.
These are all good questions, but they're not really LAMMPS questions,
and I'm not a reaction modeling expert.
So all I will say is that you can track ReaxFF bond formation/destruction
with the fix reax/bonds command and that you easily write input scripts
that will do multiple simulations while varying initial configs or temperatures
or other parameters.
If you are not already doing it, you should be looking at the "prd
command in LAMMPS" (http://lammps.sandia.gov/doc/prd.html) and the two
papers at the bottom of the link by Voter.
I think the reaction rates calculated at high temperatures can be
extrapolated to lower temperatures IF AND ONLY IF the reactions follow
the Arrhenius relation. When the reaction follows the Arrhenius
relation, the probability of the reaction is proportional to
exp(-E_{act}/K_B T), where E_{act} is the activation energy of the
reaction, K_B is the Boltzmann Constant and T is the temperature.
I hope I am right in emphasizing the "IF AND IF ONLY".
good tip - I had forgotten about PRD in LAMMPS.
We also plan to add a NEB capability soon, which
allows you to search for barriers and compute barrier
heights. But you have to know something about the
reaction end points.