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

[Justin Lemkul]

On 2/26/17 11:29 AM, Dayhoff, Guy wrote:
> 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.
>

Your approach to "fix" the Drude coordinates is inappropriate.  You're negating 
(1) the purpose of equilibration by resetting the coordinates and (2) the 
polarizable model entirely.

The coordinate file provided on our website is already equilibrated and runs 
indefinitely for me.  I don't know why your system is failing.  You are using 
the latest checkout of the Drude branch, and you are running using only OpenMP, 
not domain decomposition?

If that is the case, please send me a tarball of all your files, a shell script 
to show your exact sequence of commands, and we can debug off-list.

-Justin

-- 
==================================================

Justin A. Lemkul, Ph.D.
Ruth L. Kirschstein NRSA Postdoctoral Fellow

Department of Pharmaceutical Sciences
School of Pharmacy
Health Sciences Facility II, Room 629
University of Maryland, Baltimore
20 Penn St.
Baltimore, MD 21201

jalemkul at outerbanks.umaryland.edu | (410) 706-7441
http://mackerell.umaryland.edu/~jalemkul

==================================================