[gmx-users] Drude-2013 FF, hardwall constraint errors with swm4-ndp wate

[Dayhoff, Guy]

Hi,

   I'm experiencing an issue with the hardwall constraint used with the Drude FF
during NVT runs employing Lagrangian dynamics which I have yet to rectify on my
own.

In short....

I'm using the GROMACS dev release obtained via:
 git clone git://git.gromacs.org/gromacs.git
 cd gromacs
 git checkout drude

I'm using the drude-2013b_16-sep2016 ff available from MacKerell's site here.<http://mackerell.umaryland.edu/download.php?filename=CHARMM_ff_params_files/drude-2013b_16sep2016.tgz>

I have downloaded the equilibrated swm4-ndp waterbox available from MacKerell's
site here<http://mackerell.umaryland.edu/download.php?filename=CHARMM_ff_params_files/swm4-ndp_1024.gro>.

I have also created 3 .mdp files, dmin.mdp, dnpt.mdp and dnvt.mdp. These were all
simple copies of the recommended settings from the supporting material (here<http://onlinelibrary.wiley.com/store/10.1002/jcc.23937/asset/supinfo/jcc23937-sup-0001-suppinfo01.pdf?v=1&s=f470e09c8821d34bbf1ab29c8ab2f1aa991f374a>) of the
paper on "Implementation of Extended Lagrangian Dynamics in GROMACS
for Polarizable Simulations Using the Classical Drude Oscillator Model"

Next I generat a waterbox with:
gmx solvate -box 3.0 3.0 3.0 -cs swm4-ndp_1024.gro -o swm4box.gro

Minimization, and NPT proceed without issue:
gmx grompp -v -f dmin.mdp -c swm4box.gro -o swm4box-EM.tpr -p swm4box.top
gmx mdrun -v -deffnm swm4box-EM
gmx grompp -v -f dnpt.mdp -c swm4box-EM.gro -o swm4box-NPT.tpr -p swm4box.top
gmx mdrun -v -deffnm swm4box-NPT

Next I generate, without issue, an index file to id the two new groups, Drudes and Atoms:
gmx make_ndx -f swm4box-NPT.gro -o swm4box.ndx
a D* | a LP* (makes new group #3)
!3
name 3 Drudes
name 4 Atoms
q

Finally I generate the .tpr file and attempt the NVT run:

gmx grompp -v -f dnvt.mdp -c swm4box-NPT.gro -o swm4box-NVT.tpr -p swm4box.top -n swm4box.ndx
gmx mdrun -v -deffnm swm4box-NVT

At this point I receive errors derived from the hardwall constraint.

HARDWALL: step 0 - imposing constraint on atom 4275 rab = 1.245179
HARDWALL: step 0 - imposing constraint on atom 1235 rab = 0.591511
HARDWALL: step 0 - imposing constraint on atom 865 rab = 0.798873
HARDWALL: step 0 - imposing constraint on atom 3260 rab = 2.016762
HARDWALL: step 0 - imposing constraint on atom 3375 rab = 1.581781

Fatal error:
Drude atom 4275 is too far (r = 1.245179) from its heavy atom 4271.
Cannot apply hardwall.

Fatal error:
Drude atom 1235 is too far (r = 0.591511) from its heavy atom 1231.
Cannot apply hardwall.

Taking into consideration the error I looked at my working drude model of glucose which does not
have this error with the hardwall constraints and wrote a program to fix the positions of the drude
particles causing the errors (specifically the DOH2s). Here I copy the position of the heavy atom (OH2)
to the respective drude particle and alter the y value by + or - 0.001, generating a "fixed" version
of the aforementioned swm4-ndp_1024.gro, named swm4-ndp-fixed.gro.

Repeating the above steps exactly, using -cs swm4-ndp-fixed.gro in place of -cs swm4-ndp_1024.gro
when generating the waterbox produces similar, yet different errors.

HARDWALL: step 0 - imposing constraint on atom 1295 rab = 1.924643

Fatal error:
Drude atom 1295 is too far (r = 1.924643) from its heavy atom 1291.
Cannot apply hardwall.

Interestingly, when I look at swm4box-NPT.gro and read the position of the drude attom and it's respective
heavy atom (1295 and 1291) I see:

  259SOL   DOH2 1295   1.059   1.475   0.637  0.0000  0.0000  0.0000
  259SOL    OH2 1291   1.066   1.480   0.632  0.4306  0.2402 -0.8126

I’m not surprised here as I have ran the system through minimization and NPT at this point so
again here I employed my program to fix the drude positions and obtain a "fixed" file showing:

  259SOL   DOH2 1295   1.066   1.481   0.632  0.0000  0.0000  0.0000
  259SOL    OH2 1291   1.066   1.480   0.632  0.4306  0.2402 -0.8126

Generating a new .tpr from this "fixed" NPT.gro and running still fails and produces similar errors.

HARDWALL: step 0 - imposing constraint on atom 3330 rab = 1.448386

Fatal error:
Drude atom 3330 is too far (r = 1.448386) from its heavy atom 3326.
Cannot apply hardwall.

Once again I check on the position of the implicated drude & heavy atom and see:
  666SOL   DOH2 3330   1.539   2.952   0.319  0.0000  0.0000  0.0000
  666SOL    OH2 3326   1.539   2.951   0.319 -0.0639 -0.0928  0.8176

My fix is intact here, and from what I can ascertain they should not have an r = 1.45...

At this point i'm not sure where to go next, and recognize that may be due to my
own limits with respect to knowledge of how GROMACS works. If this is simply a case
of a bad starting waterbox and there is a working one out there i'd appreciate a link to
download it.

I've uploaded the files used to generated these errors in case they prove useful. They can be
found: here<https://www.dropbox.com/sh/974n6p3d2mrrrgc/AADVChq4h4FP77tePM0K8n9Ka?dl=0>, where swm4box0-NPT.gro was made using the swm4-ndp_1024.gro (as it was
when downloaded), swm4box-NPT.gro was made using the swm4-ndp-fixed.gro I created, and
swm4box-NPT-FIXED.gro was made from swm4box-NPT.gro using my program to adjust drude
positions.

Thank you for your time.

   My Best,
     Guy Dayhoff