Bond atoms missing in SWM4-NDP water with alkane system

LAMMPS LAMMPS General Discussion
1

Hi All,

I have a system of polarizable water, SWM4-NDP with methane which I am trying to equilibrate. I am using thole damping since polarization catastrophe is common with this water model, but I still keep getting “Bond Atoms Missing” error. The offending atoms are the oxygen and its drude. I have reduced the timestep to even a ridiculous amount (0.0001 fs) and that doesn’t fix it. I have even increased the cut-off from 12 to15 Angstrom and that just delayed the error by a few more steps. I don’t know what else I can try and need some suggestions. I have attached my input and data files. Thanks a lot.

in.equil-run5 (1.9 KB)
data.water+methane-from-restart (9.8 MB)

2

The data file cannot be downloaded. Also, for debugging purposes it would be much better, if you created a much smaller test system, that exhibits the same behavior.

3

Thank you very much, Axel. I took your advice and created a much smaller system as a test: 33 methane and 3375 water molecules. It took some time, but I was able to equilibrate the system and it’s even stable at 1.0 fs. Meanwhile, the system of interest: 117 methane and 11638 water still remains unstable.

I have attached the data files for the test system and re-attached the main system for perusal.
data.water+methane-from-restart.bak (9.8 MB)

Edit: I am having trouble with the test data file. I attach it and then get file unavailable when I want to verify that it’s correctly attached.

data.water+methane-small-test (2.8 MB)

5

Try to rename it. e.g. from data.xxxx to xxxx.data.

6

I think I have found the proverbial “smoking gun™” that could be the cause of your problems.

Your input uses the “delete_bonds” command and you can see the following output in your log:

Deleting bonds ...
  47020 total bonds, 12106 turned on, 34914 turned off
  12340 total angles, 702 turned on, 11638 turned off
  0 total dihedrals, 0 turned on, 0 turned off
  0 total impropers, 0 turned on, 0 turned off

As you can see from the information, the bonds and angles are turned off but not really removed. This is done internally by flipping the sign on the bond and angle type. Those will not be contribute to the particle forces, but the definition remains in place. Since the bond forces are not computed, atoms may move further apart and eventually they can separate by more than a sub-domain’s length and then you get the error you see. With a larger system, this situation is more likely, hence your success of simulating the smaller system.

Now, if I add the “remove” keyword to that command the output changes to:

Deleting bonds ...
  12106 total bonds, 12106 turned on, 0 turned off
  702 total angles, 702 turned on, 0 turned off
  0 total dihedrals, 0 turned on, 0 turned off
  0 total impropers, 0 turned on, 0 turned off

and the bonds and angles are actually gone and not present and turned off.

But if you study the entire documentation of the delete_bonds documentation, you should learn that you also need to add the special keyword to trigger recomputing the excluded non-bonded
interactions. And by looking at the thermo output you can see that this has a significant impact, considering how many interactions you have removed.

I will change the documentation to make the paragraphs discussing those issues more visible and also mention the possibility of turning off without removing has to cause errors.

7

Thanks a lot, Axel. I appreciate your looking into this. I have updated the “delete_bonds” command with the “remove” and 'special" keywords. I definitely can see a change in the thermo output. But the issue still persists. I found that the temperature of the drude atoms increased exponentially, getting up to 50,000 K when its supposed to be cold at 1 K. This of course, causes the kinetic energy to go up and the whole thing crashes after a few timesteps.