LAMMPS WWW Site - LAMMPS Documentation - LAMMPS Mailing List Archives
Re: [lammps-users] Problem with using external electric field in reactive MD simulations
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: [lammps-users] Problem with using external electric field in reactive MD simulations

From: Axel Kohlmeyer <akohlmey@...24...>
Date: Fri, 4 Aug 2017 07:50:33 -0400

On Fri, Aug 4, 2017 at 1:31 AM, Mohammad Izadi <izadi0511@...1125.....> wrote:

Dear Axel Kohlmeyer

Thank you for your fruitful reply to my questions and also about your advices. I am indebted to you. I think our discussion about the external electric field in the reactive MD simulations can help to clear the subject and also find good solutions to the probably problems.  

About the energy conservation along with the external electric field (or force) I think when we are imposing an external electric field (or force) continually to the simulation box (if we have charged particles in the system, this is for the efield), actually we are injecting some energy continually to the system and it is clear that the total energy of the system is not conserved during the simulation. In the other words, the system is not isolated and it is closed.

there should be no doubt at all, that adding an external force will gradually increase the kinetic energy of the system. this easily follows from newton's laws of motion: consider a single charged particle at rest with no other interactions. by adding a constant force, it will be continuously accelerated and thus pick up kinetic energy. 

I questioned this subject (in my last email) because of this reasons:

1. We had written an article (entitled: a simulation with ReaxFF by applying an external efield), but the referee made a major revise due to two reasons that one of them was about the problem with conservation energy along with the ReaxFF simulation when we use external efield. He/she told “The most important thing in reactive molecular dynamics simulation is that the total energy of the system is conserved. If the ReaxFF program uses the Verlet method for the integrator, it is impossible to preserve the total energy of the system, and the calculation results shown in this paper cannot be verified.The author need to use the sixth symplectic integrator with a time step of 0.1 fs. the symplectic integrator is the numerical integraton scheme for Hamiltonian system, which conserve the symplectic two-form exactly, so that (q(0), p(0)) → (q(τ), p(τ)) is canonical transformation. This algorithm is accurate and has no accumulation of numerical errors for total energy in constrast to the other common algorithm to solve the Hamiltonial equation of motion.   We know that ReaxFF uses the velocity Verlet, not the ‘normal’ Verlet integrator.

​i find it confusing to argue about energy conservation, when you are discussing a simulation setup, where energy conservation is not possible.
whether you are using ​ReaxFF​ or some other force field, that is independent from the choice of integrator. LAMMPS indeed uses a velocity verlet algorithm, which is a second order symplectic algorithm, but so is the regular Verlet algorithm.

i cannot help you argue your review since i neither know your research, your paper, any references that support the point your reviewer is making nor am i interested in learning about it. i do have my opinions, but without spending the time to learn the entire context and support them with suitable references, they have no place in a scientific discussion.

2. In reactive MD simulations by applying continually constant electric field, when charged particles are crossing the PBC I seen a jump in kinetic energy. I think this is because the opposite sign of the efield at the in front wall in direction to the efield. I did this with ReaxFF Fortran source code.

​this observation and explanation does not make sense to me, however - as noted above - i cannot argue your research based on only the minimal and vague information your provide. there are lots of ways to have simulations behave in unexpected ways and lots of ways to make mistakes.

but since we were discussing ​using fix efield and pair style reax/c this is also not relevant.

3. In a simple system I checked the differential total energy in alternatively similar delta t during the simulation with continually constant efield but I didn’t gain the equal total energy differences. I think this is due to continually charge equilibration (and changing the atomic partial charges) in each time step and also change the position of the particles in the box. So, we have different forces and energies due to the external electric field and also different delta total energy in each delta t. I did this with ReaxFF Fortran source code.

​i don't know the force reaxff implementation and it is deprecated anyway. the kind of test you were making only makes sense, when you turn off charge equilibration. i don't understand what you are trying to prove here.​

It is not clear for me that exactly/truly what happens for the energy when we use external electric field in the Reactive MD simulation! If you have any suggestion or criticism to the causes and effects raised above, I gladly welcomed.

​sorry, but i don't have much additional insight outside of criticizing your confusing ways of presenting your arguments. on one hand, you want to discuss fix efield with reax/c, but then your always bring up the fortran version of ReaxFF. also, it is impossible to provide specific answers on vague descriptions. it seems to me, you need to go back and do some serious basic testing. these tests should be done with simple systems (just a few handfuls of atoms) and then you should first establish what kind of settings are required to maintain a suitable energy conservation with your system *without* the external field. then you need to verify, that adding the electric field is actually reproducing the physics that you want to study. in case you encounter inconsistencies or unphysical behavior, you should come back to this list and provide the re​levant files and information to reproduce it so it can be investigated, whether this is a bug, a mistake in your input, or the correct behavior (i.e. a mistake in your expectations).



Thanks in advance for your help and feedback.

Best regard


Mohammad Ebrahim izadi,

Department of Chemistry,

Tehran University,

Islamic Republic of Iran,

Phone : +98 – 21 – 61113358

Fax :  +98 – 21 – 66409348

Dr. Axel Kohlmeyer  akohlmey@...24...
College of Science & Technology, Temple University, Philadelphia PA, USA
International Centre for Theoretical Physics, Trieste. Italy.