Calculation freezing when moving from 'boundary p p p' to 'boundary p p s'

Dear community, I would really appreciate your help in solving the problem with calculation freezing when moving from ‘boundary p p p’ to ‘boundary p p s’.

I study the impact of N2 molecule on diamond using many-body (Tersoff type) and ReaxFF pair potentials.

The diamond target is divided into middle and edge groups, where the edge group (shown as blue atoms in the figure) is used to prevent vibrational energy from being reflected from the edges.

First, periodic (p p p) boundary conditions is used. At that conditions minimization of the energy, then NPT, and NVT equilibration are performed. Then boundary is changed to p p s, which is necessary for bombarding the surface in the z direction.

If I use the Tersoff type many-body potential, everything works fine. As soon as I switch to ReaxFF, the computation freezes when going from the p p p boundary to the p p s boundary.

This means that the task keeps running forever without any computation results. Please see the attached input and output files for more details. Note, I tried to change the fix group from ‘middle’ to ‘all’ but it did not help.

Following Axcel’s (to Fernanda_S_Teixeira advice from May 2017) I changed from “boundary p p s” to “boundary p p m”. Thus, everything works fine. I just want to remind that there are no any problems when using the Tersoff type potential.

There are several questions that I would very please for your assistance. What probably I do wrong when use “boundary p p s” and what could be a principial difference when switch to “boundary p p m” so the tasks run with no freezing.

I don’t see this in the LAMMPS manual. What am I missing? In the LAMMPS manual I could see that not only “f” or “m” can be used, but also their combination “fm”. Unfortunately, I don’t understand the physical meaning of the combination “fm”. Thus, could I probably use “sm” or “sf”? But I’d rather understand the combinations before using them.

Thank you and Kind regards,

Victor


N2dia.in (6.8 KB)
N2dia.out (17.1 KB)
N2.dat (198 Bytes)
lmp_control (1.1 KB)
log.N2dia (11.1 KB)
ffield.reax.chon2019 (8.6 KB)

Why?

The main difference is that with p p s your box will instantly shrink where you have extra vacuum. If that vacuum region is large and you are running in parallel, this can lead to lost atoms. With p p m this will not happen, since the m boundary is like s with a minimum, so the box will not shrink below the initial value. This is all documented in the boundary command documentation.

There is a significant difference between the Tersoff and ReaxFF implementation. ReaxFF uses a different memory management method, since it needs to store lots of internal data about the dynamics bond topology including bonds, angles, dihedrals, hydrogen bonds and more. Tersoff doesn’t require this, only a “shortened” neighbor lists for close pairs. Thus ReaxFF can have corrupted memory if the environment changes by a large margin and that is exactly what happens when you are switching from periodic to shrinkwrap and run in parallel.

The different boundary conditions are explained in detail in the documentation for the boundary command.

The sensitivity of LAMMPS’ ReaxFF implementation to significant changes of the environment has been discussed many times, often on a weekly basis. Just check out older posts here. The remedy for that is to use the KOKKOS implementation (can also be compiled without GPU support and with or without multi-threading support.

Again, study the documentation for the boundary command carefully. All the information you need is there and in great detail.

1 Like

Dear Axel,
First, thank you very much for your assistance and the detailed explanations!

Could I ask you to explain me some details that are not entirely clear to me yet?

You say “Why?” to my note on the need to change ‘p p p’ to ‘p p s’.
The reason of that was that I suggested that a periodic boundary on Z bombardment direction would not be a correct approach. It probably shouldn’t be ‘s’, but it shouldn’t be ‘p’ either. Am I correct?

Does it follow from your explanation that since ‘p p m’ in the case of my simulation task is fine for ReaxFF it might also be preferable when using Tersoff?

“Again, study the documentation for the boundary command carefully. All the information you need is there and in great detail.”
Yes, you are absolutely right! I missed the first paragraph, so unwise.

Would you say that if the lower face (bound) is supposed to be periodic (as in my simulations) and the upper one is not, it would be reasonable to use “ps” or (as in my situation with ReaxFF) “pm”?

Thank you and kind regards,
Victor

This doesn’t answer my question. What does “correct” mean in this context? Any choice you make in your input has to be done for a reason. Choices in MD simulations are rarely of the “correct or incorrect” kind, and if they are you have to be able to tell what makes the choice correct.

Overall, why choose periodic boundaries at all, when you create an isolated system? …and is an isolated system the correct model for your simulation. … and running fix npt on such a system makes no sense at all and doing isotropic box deformations with it even less so.

Bottom line, you are making lots of choices in your input, but many seem arbitrary and conflicting with best practices or generally the physics of your model. It seems you have not discussed this with somebody that has sufficient experience and have not received sufficient training to do this on your own. This is not something that can be addressed from remote and by discussing here in the forum. You need proper in-person tutoring.

Most of the remaining questions are also along the same lines and are off-topic for this forum, since they are not about LAMMPS but about how to do your research.

A bunch of very useful advices, thank you!