I am running a phase equilibrium simulation where after compression and equilibration in a cubic box, the box is elongated 4 times the original length in the z-direction and temperature is quenched.
Instantly, simulation was terminated with error message “ERROR on proc 1: Bond atom missing in image check”. I did some reading in the forum and found that I might need to include the following line:
comm_modify mode single cutoff some number
to increase the ghost cutoff distance (because my box dimension is significantly larger than my pairwise cutoff: in real units, initial 90x90x90, compressed to 60x60x60, elongate in the z-direction to 60x60x240, with lj/cut 14 and neighbor skin = 2).
I have been reading a lot that the default settings in lammps should be able to handle most well set-up simulations. I am doubting if I am doing the right thing and hence would like to have some advice from the experts and community.
The box dimensions as such do not cause the problem with “bond atom missing”, but rather if you are stretching a bond too far.
So the issue is probably more likely due to what happens when you elongate the box. By default (group argument is “all”) the entire system is scaled to the new box size, but it sounds that perhaps you only want to increase the box without expanding the - previously compressed - system accordingly. This can be changed by creating an empty group group none empty and then use that “none” group for the change_box command.
If my guess is not correct, then you have to reconsider how you expand the system because stretching bonds larger than the cutoff distance will be causing high potential energy and very unphysical behavior when the atoms are going to “snap back” assuming you work around the cutoff issue by expanding the communication cutoff.
I used the deform command with remap none (I have attached my run command lines below). I read the change_box command a while ago, but I was unsure as to how it should be used in between fix and unfix.
fix 1 all nvt temp 800 800 100.0
run 1000000
unfix 1
fix 1 all nvt temp 800.0 800.0 100.0
fix 2 all deform 1 x final 0 65 y final 0 65 z final 0 65 units box
run 5000000
unfix 1
unfix 2
fix 2 all deform 1 z final -65.0 130.0 remap none units box
run 1
unfix 2
fix 1 all nvt temp 400 400 100.0
run 15000000
unfix 1
what you are doing with fix deform followed by a run 1 is “abuse” of the fix deform command. The change_box command was meant to be used for exactly those purposes (and doesn’t require a run because it is immediate and not gradual).
the “remap” keyword is not preventing the scaling of coordinates. That should apply to the fix group ID atoms same as for the change box group ID.
in either case, it would be best practice to visualize the system before and after the change to see whether what you want to do as a system manipulation is actually implemented as expected. In fact visualization is often the fastest and simplest approach to confirm that a command works as expected (and confirming those steps should become second nature. working under false expectations can result in a huge waste of effort and resources. like the saying goes in other crafts “measure twice, cut one”, you have here “visualize twice, run once”).
Thank you for the advice. The reason why I brought this to the forum for help is because it actually worked fine for 1-bead and 2-bead molecules in LJ and real units. I have been getting the abovementioned error after changing into real units and tested it on larger (> 20 atoms) linear or cyclic molecules.
I also just tested using change_box as below, but the same error arises. (I think I am using the command right)
fix 1 all nvt temp 800.0 800.0 100.0
fix 2 all deform 1 x final 0 65 y final 0 65 z final 0 65 remap none units box
run 100000
unfix 1
unfix 2
change_box none z final -65.0 130.0
fix 1 all nvt temp 400 400 100.0
run 5
unfix 1
Total # of neighbors = 10894921
Ave neighs/atom = 581
Ave special neighs/atom = 8.5
[data.vle|attachment](upload://z84Ww3SqgnqKNGZFYMIooldskpF.vle) (2.8 MB)
[in.lammps|attachment](upload://yrfIZ3lqStWRgStmYwC0KDd7171.lammps) (2.6 KB)
Neighbor list builds = 11626
Dangerous builds = 0
Changing box ...
orthogonal box = (0 0 -65) to (65 65 130)
Setting up Verlet run ...
Unit style : real
Current step : 100000
Time step : 1
ERROR on proc 33: Bond atom missing in image check (src/domain.cpp:765)
Last command: run 5
ERROR on proc 37: Bond atom missing in image check (src/domain.cpp:765)
Last command: run 5
ERROR on proc 1: Bond atom missing in image check (src/domain.cpp:765)
Last command: run 5
ERROR on proc 41: Bond atom missing in image check (src/domain.cpp:765)
Last command: run 5
...
There is one additional twist that you have to keep in mind, that may cause issues here.
When you equilibrate your system and have periodic boundary conditions, you can have bonds crossing periodic boundaries. That will result in atoms participating in a bond having different values for the image flags. Those atoms are still be displaced by the change_box or fix deform command even though their coordinates are unchanged. This is because the actual unwrapped coordinate includes adding box*imageflag to the position that that is changed by changing the box.
This is not so easy to correct. Perhaps you can use write dump and write out unwrapped coordinates on the compressed system, then change the box dimensions with change_box and use read_dump to read the coordinates back.
Failing that, it may be possible with some creative use of atom style variables, or post-processing of a trajectory or data file.
However the core of the problem is that your workflow should be implemented differently. Rather than compressing the box itself, you should start with a box in the desired dimensions and then add walls and move the walls to get your compressed/condensed configuration.
Thanks for all the suggestions. I will try out the wall method.
In cases where the “change_box” method work (for example, simple 1 or 2 or 3-bead molecules), in order to restrain the bulk liquid phase from drifting away from the center of the box, is fix recenter the correct way to do this? (I don’t think I should be using fix spring after reading the doc.)
...
(change_box)
...
fix 1 all nvt temp 400 400 100.0
fix 2 all recenter INIT INIT INIT units box
run 5000000
unfix 1
unfix 2
It will still create high potential energy configurations, which invalidates the whole point of your system preparation up to that point.
No. Fix recenter translates the point of reference. That is not what you want.
Fix spring, but only in directions with vacuum and with a very weak force constant should be the correct choice.
If you set up and equilibrate your system correctly, you should only be required to remove a center of mass drift once with the velocity command at the end of equilibration and it should retain the position.
I agree with Axel – the main issue here is almost certainly that change_box will blow up any molecules that happen to straddle boundaries in the dimension you change.
If I were you / your supervisor, I might worry about whether the simulations are fundamentally physical. A system simulated in a periodic box of width 1 (whichever unit) is fundamentally different from a system simulated between walls 1 unit apart, since inter-particle RDFs are different for a particle near a wall as compared to a particle in “the bulk”.
As another issue, there is no “natural” center of mass in a PBC system, without choices about what happens when a particle crosses the boundaries. It is even less “natural” to define a center of mass in a molecular PBC system (consider a molecule straddling a boundary). So a fix spring or fix recenter may not do what you think it does.
Thank you for giving me a deeper understanding of this - as I was just trying to implement the methods/workflow I mentioned above which I got from literature and figuring out the right commands to execute them using LAMMPS (because they either used GROMACS, D-POLY or HOOMD).
If I may ask you as an expert in MD and coarse-graining, where researchers in this area compute phase equilibria using the compress-then-expand-and-quench or temperature-quench-only method, is there one that is preferred over the other?
There is also another issue on “bond atoms missing on proc”, which I have been reading up and down the forum. In a temperature-quench-only simulation (with pairwise cutoff = 14 A, timestep = 1, as in the input file attached previously), things work fine. However, the error would sometimes arise if the pairwise cutoff is increased to 18 or 20 (because literature mentions that pairwise cutoff should be at least 6\sigma_{max}. What I learned from the forum so far is that the main cause of such error is due to bad dynamics (either the forcefield is incompatible, or neighbor checks do not happen frequently enough). I couldn’t understand why increasing the cutoff would make the dynamics bad that bonds could potentially be stretched beyond this cutoff…
Errors like bond atoms missing or lost atoms can happen for multiple reasons and it would be a mistake to collate them all into one. What is being explained in most discussions is that the most common cause is a bad choice of potential parameters and/or simulation settings (e.g. too large a timestep), that will cause atoms to come very close and then be attracted or repelled with very large forces and then - occasionally - so fast that they move through the “ghost atom region” between two neighbor list updates and then cannot be communicated between subdomains or periodic images.
Since it is a requirement that the ghost atom region (= communicaton cutoff) is at least the size of the largest pairwise cutoff plus neighbor list skin, it is usually large enough to contain all constituent atoms of bonds, angles, dihedrals and impropers, even if the atom “owning” that interaction is local and the others are ghosts.
Problems can arise when the box dimensions change a lot, e.g. when you have shrinkwrap boundary conditions but leave a lot of empty space, so that the simulation cell will shrink a lot in the first MD step, or when you expand a box by a large margin without enforcing relative atom distances within bonded interactions (e.g. when you scale all coordinates or have different image flags on different atoms).
Your description is too unspecific to make an assessment as to why you see the unexpected loss of bond atoms. The fact that the most common reasons for missing atoms do not apply does not mean that your input is correct.
I don’t regularly calculate phase diagrams in my work so I can’t give you detailed advice on how to do so – and even if I did, what someone says on an internet forum is not a substitute for a good literature review. When someone reviewing your paper asks why you did what you did, “someone told me on Matsci.org” will not be an acceptable answer!
You would be best off contacting the authors of the papers from which you got those workflows. Every MD software has its own quirks and features and there is a lot of implicit knowledge that is not communicated in a Methods section (which is bad for science – that’s a different story).
For example, in GROMACS, you cannot change box sizes so dramatically in the course of a single simulation. So anybody trying to do what you seem to be doing in this workflow would have run a simulation at one box size, stopped the run, post processed the configuration in a separate tool to change the box, and then run a second simulation. Such a workflow would ensure that the configuration is “good to go” for the quench simulation – and by trying to do both simulations in one LAMMPS script you are robbing yourself of the opportunity to do the same checks.
Before discussing specifics any further – have you visualised your system? Do you know what it looks like before you perform change_box and do you know what it looks like just after? A picture can be worth a thousand words.
Many thanks to both Axel and Shern for the kind and detailed explanation. This has been very very helpful and insightful especially when I have to do some MD with a systems engineering background.