Shake atoms lost on proc X (bad dynamics)

[clf_x]

Dear all,

I am working with a system that contains X number of polymer chains and Y number of water molecules.

I wanted to insert an additional 1000 water molecules into the simulation box. To do this, I first changed the box size to the closest (larger) integer value for each dimension (without remapping). I then used ‘create_atoms’ to insert all new water molecules randomly in the box - using a molecule template. I deleted all water molecules within a 1 Angstrom cut-off of any other atom, repeatedly using ‘create_atoms’ until all waters were inserted with no further deletions.

This seemed at first successful. I minimised the system to a low potential energy and force tolerance and a further simulation at 300K with fix NPT ran without error (or warning). The problems started when heating to 400K. I get the error ‘lost shake atoms A B C on proc X’ semi-consistently somewhere in the middle of the simulation (i.e., after 75000 1fs timesteps).

I have included two snapshots of the simulation box (each from a trajectory file) below. The first is of the system before additional waters are inserted and the second is after. Both have been wrapped by residue/molecule and I am just showing water molecules for clarity here.

It’s clear that there is some issue with incorrect image flags for atoms water molecules at the periodic boundary leading to the appearance of bonding between different molecules when the trajectory is visualised. It’s worth noting that VMD reads the bonding for the trajectory files incorrectly (assigning a bond between water hydrogens for example) and that the topology file itself shows the correct bonding.

Would greatly appreciate any help/advice with this.

Thanks

[simongravelle]

Hello,

This is a risky move on a system with periodic boundary condition (molecules located at the border will be completely deformed):

@clf_x

I first changed the box size to the closest (larger) integer value for each dimension (without remapping)

Does your procedure works if you skip this step, and simply add the water molecules?

Simon

[akohlmey]

@clf_x

Would greatly appreciate any help/advice with this.

When reporting problems, please always report which LAMMPS version you are using and what platform you are on. It also often helps to spot issues, if you would provide some of the log file output, e.g. from the startup and close the the point where the run errors out.

@clf_x

I get the error ‘lost shake atoms A B C on proc X’ semi-consistently somewhere in the middle of the simulation (i.e., after 75000 1fs timesteps).

First you need to observe the trajectory and see, if there are any unusual developments with respect to geometry or the involved molecules. If not, you need to check what happens briefly before the error. Since this happens after some time, they would suggest that there is some issue with your force field parameters, i.e. the LJ part is not sufficiently repulsive for the Coulomb and then atoms can get too close (and closer due to heating) and then forces get too large and atoms “explode”. This often happens when mixing matching different force fields and water models.

@clf_x

It’s worth noting that VMD reads the bonding for the trajectory files incorrectly (assigning a bond between water hydrogens for example) and that the topology file itself shows the correct bonding.

The “incorrect” bonding would happen only when you load the trajectory without loading the data file first. The trajectory dump does not contain any bonding information and VMD’s heuristics are limited. When loading the data file first and then the trajectory dump on top of that, you should retain the correct bonding (unless you have incorrectly wrapped positions or inconsistent image flags).