1. Can I use lammps to chemical reaction systems?
Please note that you can only get as good an answer as you question is. nobody can tell without knowing what kind of chemical reaction. most likely not.
In most cases that i know of you need to do ab initio MD, or something similar.
2. I know lammps can be used to chemical reaction systems(e.g. intramolecular cross-linking of BCB/styrene copolymers , Macromelecules,2009,42,8534-8542).I can't just add some bonds to the system directly and say that the reaction is completed. I need to let itself form new bonds. How can I do it?
There is a fix bond/create command which will cross-link, but it is for idealized (e.g. bead-spring) models. Adding cross-links to all-atom polymer chains would require a new, more sophisticated fix, e.g. to
turn on angle and dihedral interactions when the bond is formed.
----------------------------------------------------------------------
Thanks, steve. I am also curious about simulating all-atom cross-linking system with the ReaxFF in LAMMPS. I am wondering if the bond breaking and forming could go on automatically with ReaxFF.
As it is developed, ReaxFF is meant to take angle and dihedral interactions as a function of bond orders. I perceive it as an spontaneous process of reaction. Is that ture? If not, what we could actually do with the ReaxFF in LAMMPS since no commands have come up to assist the reaction process? Looking forward to your reply. Thanks in advance.
--------------------------------------------------
Bond breaking and creation in automatic in ReaxFF and in all bond-order potentials. There is a fix reax/bonds diagnostic command to monitor what happens, but you don't control it.
3. 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.
------------------------------------------------------------------
A few comments in addition to what steve wrote.
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.
---------------------------------------------
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.
4. I have a ionic crystal (LiMn2O4) as a chatode material in a battery.i define this structure in lammps and equilibrate it. In the second stage I have to simulate the electrochemical processes that occure In a battery.for example in battery charging with applying a electrical potential difference, the Li ions migrate from LiMn2O4 structure. My problem is that: How can i apply a virtual electrical potential difference on this crystal as a cathode? Actully I want to see Li ions migration from this structure .Are the Fix efield or fix addforce or any other commands effective for my simulations?
You should study the literature to see, if there are existing studies that have done similar things and follow that examples.
It has been a while since i have talked to people studying electrochemical systems like yours, but my impression was that the process as a whole is happening on time and size scales far beyond the reach of classical MD (discounting the issue of how you would model the electrochemical reactions which could blow up the computational cost by additional order(s) of magnitude). typically MD would be used to gain atomic level insight into individual steps and provide parameters for models like KMC that would then be used to study the chemical reaction and transport related parts of the process.
Again, the best place to look for insight are text books and the published literature. you also need to discuss with your adviser. this is not so much a question of how to technically realize this, but rather of planning a research project and following a viable strategy.
参考链接：
[1] [lammps-users] lammps chemical reaction; [lammps-users] use lammps to chemical reaction systems;lammps-users Digest
[2] [lammps-users] Statistical analysis of reaction pathways; [lammps-users] MD simulation of Electrochemical Processes