Problem with using external electric field in reactive MD simulations

[Mohammad_Izadi]

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. 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.

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.

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.

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.

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

[akohlmey]

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).

​axel​