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

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

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.

Best regard

Department of Chemistry,

Tehran University,

Islamic Republic of Iran,

Phone : +98 – 21 – 61113358

Fax :  +98 – 21 – 66409348

On Thu, Aug 3, 2017 at 4:45 PM, Axel Kohlmeyer wrote:

Dear Axel Kohlmeyer

My purpose is also researching and I have known that the mailing list goal is validating and developing models to make confidence about them.

​where is that written? i am answering on mailing lists, and particularly this one for many years, and other people and me have consistently ​told people, that this mailing list is not to validate and approve other people's inputs. if a command does not behave as documented or if someone has a specific question about how something is implemented, that is a topic for the list.

whether your simulation setup is valid for your research is *your* responsibility. it would become my responsibility only in the case of you submitting a paper to a journal and me being asked to review it or if i was part of your thesis comittee or similar.

I have been investigating on molecular dynamics simulation (exactly on the reactive force field ReaxFF) since four years ago. In my MSc thesis I spend the whole time to investigate the ReaxFF potential functions, it applications. Finally I did some reactive simulations on the Flame Synthesis processes (with and without imposing the constant external electric field on my simulation box). As imposing the external electric field, I have recognized the total energy of the system is not conserved in the ReaxFF Fortran source code. I have extensive research about this phenomenon (e.g. elicit the equation of the electric field from the source code “efih=vfieldx*23.02*c1c*ch(i1)*c(i1,1)”, investigation of this equation and its constant conversion factors, how the external electric field effect on the reactive simulations, contact to experts of this field and also to the developer of the ReaxFF). And in my PhD thesis I have been deciding to use the Lammps/reax/c for my research since about one year ago. So, all of us have similar tasks.

​what has that to do with the topic at hand? i can list much more experience than you, did research and collaborated with researchers on much more complex issues involving electric fields and dielectric properties, and also most likely have written/modified more lines in the USER-REAXC​ code in LAMMPS than anybody except for the original author. these facts, do not relieve you to do your due diligence, also do not oblige anybody else to validate your work, and most certainly do not entitle you to what you have been asking.

I had learned to install and run simulations on Lammps in my own limit and I have read the whole document about electric field command, fix addforce command, compute command, variable command, region command and other related commands which I thing they are relating to the external electric field. In my last email I told you “I know that the external electric field in my simulations is a classical approximation and it just impose electric forces on my charged elements” so I knew that we are not exactly applying an electric field!!

​here you are wrong. the issue i am raising is not about whether this is a classical approximation or not. electrostatic interactions for any empirical model in LAMMPS (and that includes ReaxFF) are *always* classical. however, what fix efield does is a crude approximation even within that model, compared to, for instance, modeling the same potential difference through placing explicit charged particles  arranged in planes behind a repulsive wall.​

1- Is my electric field command in the input file correct or not?

​whether it is syntactically correct or not, you can easily gather from the LAMMPS documentation. the documentation also describes the physics that fix efield implements. if there is a discrepancy between the documentation and the actu​al behavior, this would be a topic for this list (or the issue tracker on github), however, whether this is adequate for *your* research is *your* job to determine and not anybody elses.

2- Is the total energy conserve in the lammps/reax/c simulations if we apply an external electric field or not? And what is a solution for this probably problem?

​again, whether the reax/c pair style sufficiently conserves energy or not, is something that you need to test and verify. since ReaxFF has more adjustable parameters, it may take more effort than for a simpler model. outside of the usual suspects that affect energy conservation (system setup, time step, neighbor list update settings, system manipulations through fixes), you also have to tweak multiple reaxff specific parameters, most importantly the convergence of the charge equilibration.

that you ask about energy​ conservation with an external field (or force) applied significantly discredits your claim of experience in MD simulations, as that question can be simply answered from basic considerations. if you are as seasoned and well trained a person as you claim, you should know the answer to that question.

3- Is the result of an external electric field in lammps/reax/c runs from this equation “efih=vfieldx*23.02*c1c*ch(i1)*c(i1,1)” in subroutine efield reax code or “F = E.Q” in common lammps command?

​what fix efield does is documented in the LAMMPS manual, it does exactly that and nothing else. i already stated in my previous criticism, that one of the problems in using fix efield with reax/c is, that the field contribution is not considered beyond adding the force term in the fix.

I think the second one is the answer since when I compare my results from the serial ReaxFF Fortarn code, the electric field in the lammps run is far weaker than the electric field in the ReaxFF Fortarn code in the absolutely similar conditions.

4- If the answer of the last question is the first one, why the electric field in the lammps run is much weaker?

​you are being redundant here.​
if you want to be certain, set up a simple test system, where you can compute the outcome by hand or by other external means and compare.

5- It is correct that the external electric field cannot include in the charge equilibration scheme in the ReaxFF force field potential but I think the resulted forces on the partial charged particles can cleavage or deform bonds and molecular structures and finally can affect resulted fragments or products of the chemical reactions. So, as the products of a chemical reaction are changed, the charge equilibration of the ReaxFF produce other partial charges for the atoms in molecules due to the change of atomic chemical environments. On the other words, the external electric field can indirectly affect the atomic partial charges in the molecules and fragments!

​yes, it can. but ​having an effect at all. does not mean, that the magnitude is correct.

If you can, please answer to this questions and likely other questions that I will ask in lammps-mailing-list and if you cannot, you don’t have to answer my questions.

the mailing list is a volunteer effort and not a service that you are entitled to. ​if you cannot handle (constructive) criticism, you should make an ​effort to avoid questions that raise such criticism. please keep in mind, that criticism actually is a sign that somebody cares about you doing the right thing. it would be easier to just let your questions unanswered. as with all volunteer efforts, the quality of responses you get significantly depends on what and how you ask. having your ego dented a little bit, when you are crossing the lines of what is considered appropriate, is a small price to pay for getting help at all.

​axel.​

Dear Axel Kohlmeyer

My purpose is also researching and I have known that the mailing list goal is validating and developing models to make confidence about them. I have been investigating on molecular dynamics simulation (exactly on the reactive force field ReaxFF) since four years ago. In my MSc thesis I spend the whole time to investigate the ReaxFF potential functions, it applications. Finally I did some reactive simulations on the Flame Synthesis processes (with and without imposing the constant external electric field on my simulation box). As imposing the external electric field, I have recognized the total energy of the system is not conserved in the ReaxFF Fortran source code. I have extensive research about this phenomenon (e.g. elicit the equation of the electric field from the source code “efih=vfieldx*23.02*c1c*ch(i1)*c(i1,1)”, investigation of this equation and its constant conversion factors, how the external electric field effect on the reactive simulations, contact to experts of this field and also to the developer of the ReaxFF). And in my PhD thesis I have been deciding to use the Lammps/reax/c for my research since about one year ago. So, all of us have similar tasks.

I had learned to install and run simulations on Lammps in my own limit and I have read the whole document about electric field command, fix addforce command, compute command, variable command, region command and other related commands which I thing they are relating to the external electric field. In my last email I told you “I know that the external electric field in my simulations is a classical approximation and it just impose electric forces on my charged elements” so I knew that we are not exactly applying an electric field!!

1- Is my electric field command in the input file correct or not?

2- Is the total energy conserve in the lammps/reax/c simulations if we apply an external electric field or not? And what is a solution for this probably problem?

3- Is the result of an external electric field in lammps/reax/c runs from this equation “efih=vfieldx*23.02*c1c*ch(i1)*c(i1,1)” in subroutine efield reax code or “F = E.Q” in common lammps command?

I think the second one is the answer since when I compare my results from the serial ReaxFF Fortarn code, the electric field in the lammps run is far weaker than the electric field in the ReaxFF Fortarn code in the absolutely similar conditions.

4- If the answer of the last question is the first one, why the electric field in the lammps run is much weaker?

5- It is correct that the external electric field cannot include in the charge equilibration scheme in the ReaxFF force field potential but I think the resulted forces on the partial charged particles can cleavage or deform bonds and molecular structures and finally can affect resulted fragments or products of the chemical reactions. So, as the products of a chemical reaction are changed, the charge equilibration of the ReaxFF produce other partial charges for the atoms in molecules due to the change of atomic chemical environments. On the other words, the external electric field can indirectly affect the atomic partial charges in the molecules and fragments!

If you can, please answer to this questions and likely other questions that I will ask in lammps-mailing-list and if you cannot, you don’t have to answer my questions.

On Thu, Aug 3, 2017 at 4:02 AM, Axel Kohlmeyer wrote:

Dear lammps users,

I’m using lammps/reax/c for simulating decomposition of an organic compound under a constant electric field. From outputs (I mean dump.xyz and specie) I recognized that the electric field is worked but I’m not sure that it’s exactly correct or not. Please look into my input file and make recommends if it needs. This is my input file:

​you seem to be confusing the purpose of a mailing list. it is most certainly not an input file approval or debugging/correction service​.
in other words, it is part of the job of a researcher to validate that a chosen model is giving results at the expected level of accuracy and usefulness.

please let me repeat what i have pointed out in an earlier e-mail response. when you are using fix efield, you are not exactly applying an electric field, you are adding a force to all charged atoms, that depends on the product of the (partial) charge of each atom and a given force constant with conversion factor.
this does not include the electric field in the charge equilibration (only indirectly over response to the force generated by fix efield) and there is no screening or other effects taken under consideration.

​axel.​

echo            both

units                real

dimension 3

boundary  p p p

atom_style       charge

restart 500      restart11 restart22

pair_style         reax/c NULL

pair_coeff        * * ffield.reax.rdx C H O N

compute reax all pair reax/c

variable eb        equal c_reax[1]

variable ea        equal c_reax[2]

variable et         equal c_reax[9]

variable ew       equal c_reax[11]

variable efi       equal c_reax[13]

variable eqeq    equal c_reax[14]

#define a region for efield

region    elecf block 0.0 100.0 0.0 100.0 0.0 100.0 units box

group     elf region elecf

neighbor          2 bin

neigh_modify  every 10 delay 0 check no

velocity        all create 300 235485 mom yes rot yes

fix                    1 all nvt temp 300 300 25    ###########timestep*100=25

fix             2 all qeq/reax 1 0.0 10.0 1e-6 reax/c

fix             3 all reax/c/species 1 50 250 species.txt

fix             4 elf efield 2.0 0.0 0.0 region elecf

fix_modify      4 energy yes

thermo_style    custom   step    atoms   temp   press   vol   etotal   ke   pe   v_eb   v_ea  v_et   v_ew v_efi   v_eqeq

thermo          250

timestep           0.25

dump               1 all atom 30 dump.reax.rdx

dump               2 all xyz  30 dumpnvt.xyz

run    30000

------------------------------------------------------------------------------
Check out the vibrant tech community on one of the world's most
_______________________________________________
lammps-users mailing list
lammps-users@...655....net
https://lists.sourceforge.net/lists/listinfo/lammps-users

--
Dr. Axel Kohlmeyer  akohlmey@...24...  http://goo.gl/1wk0
College of Science & Technology, Temple University, Philadelphia PA, USA
International Centre for Theoretical Physics, Trieste. Italy.

--
Dr. Axel Kohlmeyer  akohlmey@...24...  http://goo.gl/1wk0
College of Science & Technology, Temple University, Philadelphia PA, USA
International Centre for Theoretical Physics, Trieste. Italy.