[lammps-users] Bond missing in system containing Drude oscillators

LAMMPS LAMMPS Mailing List Mirror
1

Hello,

I am trying to simulate a system containing SPC/E water and polarizable ions using the DRUDE package. First, I did a non-polarizable simulation, which works well. When I add drudes to the atoms of the ions (I do that using the polarizer script “https://github.com/agiliopadua/clandpol” and then I modify the input files), however, the total energy of the system and delta of the shake statistics for angle constraint “H-O-H” increase and I receive the “Bond missing” error for a bond between a core atom and its drude particle at the beginning of the simulation. I know that it should be due to bad dynamics, but I couldn’t fix the problem neither by decreasing the time step nor by checking the neighbor list every step. Do you have any suggestions?

I use:

  • hybrid/overlay pair style (lj/cut/coul/long , coul/long/cs , thole)

  • pppm kspace style

  • cg minimization algorithm

  • tgnvt/drude thermostat (I also tried the combination of fix langevin/drude and fix nve)

  • the group of the fix shake does not include the drude particles

  • the parameters for the harmonic bonds between core atoms and their DPs are k=500kcal/A^2 and r=0

Best regards, Sakor

2

It is difficult to comment on this and make suggestions from just vague descriptions.
What may be a problem is that you have non-zero image flags from the simulation without the drude particles and that polarizer script may get confused by those.
Please carefully check for warnings when you start after adding the drude particles.

3

Dear Axel,

Thanks for the reply.
Indeed, I don’t receive any warnings when I start my polarizable simulation.
I visualized my xyz output using VMD. The strange thing is that I see some undesirable bonds that are not defined in my data file (for example, between oxygens of each water molecule).

Best, Sakor

4

A trajectory file does not contain any topology information, so tools like VMD will use heuristics to guess the bonds and those are sometimes wrong.
For VMD it is therefore recommended to create a suitable .psf format file from the data file (or first load the data file) to have the correct topology data.

A way to spot problematic interactions would be to add compute pe/atom to your input and create a custom dump with positions atom ids and the pe/atom value and then look for atoms with unusually high potential energy from a run 0.

If the adding of the drude particles has added some unwanted close contacts, you would have to try a minimization. Since minimization does not support fix shake, you would have to change the bond and angle style force constants temporarily to very high values to maintain the molecular structures.