Error in Molecular Dynamic Simulation: Particle leaving simulation territory!

[sukhito teh]

Dear ATK staffs and users,
I am trying to implement self-defined CLAY force field in molecular dynamic simulation (ATK2022), but keep getting the error "Particle leaving simulation territory!". If I replace the Si and O force field with Pedone_2006Fe2, then the simulation can run smoothly without problem. A similar problem had been posted earlier (https://forum.quantumatk.com/index.php?topic=6641.msg27814#msg27814), does the errors occur because bond potential is not suitable for molecular dynamic simulation?

[Anders Blom]

Maybe there are some parameter in the forcefield that causes very large forces on the atoms? The "leaving the domain" only is an issue when running an MD calculation in parallel over MPI, and using domain decomposition. This can provide a nice performance boost, if each domain contains 1000+ atoms, as in your case. But, I would first start to check if the potential is generally stable by just running in serial. With 27000 atoms it will still be relatively fast, esp. if you run threaded on 16 cores or so. I don't see your submit script or output, it would be useful, but general advise is only use MPI if you go beyond a single node, since the MPI and OpenMP performance intranode is about equivalent, and OpenMP is easier and potentially more stable. Also, you can save some code by using

Code

from AddOns.VASPPlugins.IO import readPOSCAR
bulk_configuration = readPOSCAR("./POSCAR")

[Anders Blom]

I checked the other post you mentioned, there was no resolution there. Are you using the same forcefield parameters? Where do these parameters come from? Is it well established that these values are correct and good? 

[sukhito teh]

Maybe there are some parameter in the forcefield that causes very large forces on the atoms? The "leaving the domain" only is an issue when running an MD calculation in parallel over MPI, and using domain decomposition. This can provide a nice performance boost, if each domain contains 1000+ atoms, as in your case. But, I would first start to check if the potential is generally stable by just running in serial. With 27000 atoms it will still be relatively fast, esp. if you run threaded on 16 cores or so. I don't see your submit script or output, it would be useful, but general advise is only use MPI if you go beyond a single node, since the MPI and OpenMP performance intranode is about equivalent, and OpenMP is easier and potentially more stable. Also, you can save some code by using

Code

from AddOns.VASPPlugins.IO import readPOSCAR
bulk_configuration = readPOSCAR("./POSCAR")

Thank you for your reply. I tried to use single domain but surprisingly the problem persists(input and output files are included). The parameters of force  fields are retrieved from this paper(https://pubs.acs.org/doi/10.1021/la504178g), although the paper use these parameters for cristobalite phase (while I used it for quartz phase),. https://pubs.acs.org/doi/10.1021/la504178g

[sukhito teh]

Maybe there are some parameter in the forcefield that causes very large forces on the atoms? The "leaving the domain" only is an issue when running an MD calculation in parallel over MPI, and using domain decomposition. This can provide a nice performance boost, if each domain contains 1000+ atoms, as in your case. But, I would first start to check if the potential is generally stable by just running in serial. With 27000 atoms it will still be relatively fast, esp. if you run threaded on 16 cores or so. I don't see your submit script or output, it would be useful, but general advise is only use MPI if you go beyond a single node, since the MPI and OpenMP performance intranode is about equivalent, and OpenMP is easier and potentially more stable. Also, you can save some code by using

Code

from AddOns.VASPPlugins.IO import readPOSCAR
bulk_configuration = readPOSCAR("./POSCAR")

I think you are right, I tried to relax the structure with the force field and the force is unreasonably large. I would look into the parameter used.

[Anders Blom]

I note that you didn't add the torsion parts of the potential in the papers, but maybe that is just a small term.

However, upon opening the forcefield from the HDF5 file in the GUI (in the Editor) I think it might be possible that your forcefield generation script is not working as intended, with respect to the tags etc. There are terms generates for "Mo" with tag "sio2", and "Si" in "mos2" and other mismatched cross-terms it seems.

[sukhito teh]

I note that you didn't add the torsion parts of the potential in the papers, but maybe that is just a small term.

However, upon opening the forcefield from the HDF5 file in the GUI (in the Editor) I think it might be possible that your forcefield generation script is not working as intended, with respect to the tags etc. There are terms generates for "Mo" with tag "sio2", and "Si" in "mos2" and other mismatched cross-terms it seems.

Thank you so much for helping me looking into the details. The torsion term (for sio2) was ignored in the referenced paper as well, and I am aware of those mismatched cross-terms, they are there because I use TremoloXPotentialSet.setTags() (similar to example given in https://docs.quantumatk.com/tutorials/combining_potentials/combining_potentials.html), I think the wrong cross-term would simply be ignored as long as I tag my structure correctly.

And if I optimize the structure before running the NVT MD, then the calculation can run smoothly (so I considered the problem solved). However,  when I later use the results of NVT MD (configuration + velocities) to continue as NVE-MD simulation, the energy is not conserved. Again, such problem doesn't exist if if use Pedone_2006Fe2 force field.

I would try using bonded potential on other system to see if the problem persist; before that, do you have any idea what can be the cause for non-conserving energy in NVE MD run?

Thank you for your attention.

[Anders Blom]

Yes, I think you are right about the tags.

I don't personally have so much experience with NVE, but have you checked with smaller time steps?
Maybe also use a tighter accuracy for the Coulomb term or try some of the other charge equilibration methods like CoulombDebye which is actually the one used in Pedone

[sukhito teh]

Yes, I think you are right about the tags.

I don't personally have so much experience with NVE, but have you checked with smaller time steps?
Maybe also use a tighter accuracy for the Coulomb term or try some of the other charge equilibration methods like CoulombDebye which is actually the one used in Pedone

Thank you for your suggestion. I had tried the above methods and the problem (energy decreasing in NVE simulation) persists. I recently used the same force field for NVE simulation using LAMMPS, though I haven't evaluate the accuracy of the results yet,  the energy of the NVE simulation is conserved, so I think the forcefield isn't the cause of problem above. It would be great if QuantumATK staffs can check if similar problems happen for MD using other bonded potentials.

[Anders Blom]

Our MD experts have looked at the details of this case. All cases where we observed that energy was not conserved in NVE, it was caused by a too large time step for the chosen potential, when it e.g. uses a spring that is too stiff, or it can also happen with ReaxFF. In the case here, we didn’t see any obvious reason, but also noticed that it is a slow drift which happens over 100s of ps. On this timescale it is not necessarily something unusual given that the numerical integration with a finite time step is an approximation and does not guarantee to conserve the energy 100% exactly. And your potential uses harmonic springs, which have fast vibrations, which would normally indicate a relatively fast time step to accurately conserve the energy.

So, the main suggestion is to lower the time step.

[sukhito teh]

Our MD experts have looked at the details of this case. All cases where we observed that energy was not conserved in NVE, it was caused by a too large time step for the chosen potential, when it e.g. uses a spring that is too stiff, or it can also happen with ReaxFF. In the case here, we didn’t see any obvious reason, but also noticed that it is a slow drift which happens over 100s of ps. On this timescale it is not necessarily something unusual given that the numerical integration with a finite time step is an approximation and does not guarantee to conserve the energy 100% exactly. And your potential uses harmonic springs, which have fast vibrations, which would normally indicate a relatively fast time step to accurately conserve the energy.

So, the main suggestion is to lower the time step.

Thank you for your advice, I would play with the time step parameters to see if the results improve.