Fix TIP4P water models when using AMBER inputs 24/4424/2
authorradakb <brian.radak@gmail.com>
Tue, 31 Jul 2018 19:44:41 +0000 (15:44 -0400)
committerDavid Hardy <dhardy@ks.uiuc.edu>
Tue, 31 Jul 2018 22:34:48 +0000 (17:34 -0500)
The effective deprecation of the "lonepairs" keyword caused some
unexpected behavior when using a parm7 in place of a psf. The main
difference is that lonepairs are now _assumed_ until they are not
found. Since parm7s do not flag lonepairs, the flag is now always
set to false.

Note that TIP4P is a special case of lonepairs, so this is not a
problem. However, the fix does seem a bit counterintuitive, so later
a more involved fix is in order.

The TIP4P lonepairs are treated as a special case.  The rigid bond
constraints routine HomePatch::rattle1old() calls tip4_omrepos()
routine to reposition the lonepair site after the atom positions
are updated on each water molecule.

Change-Id: I91f430f79c810ad0ff13999de4cdcd453f6171ff

src/Molecule.C
src/SimParameters.C

index 7a2cb91..96069c3 100644 (file)
@@ -10144,6 +10144,16 @@ Molecule::Molecule(SimParameters *simParams, Parameters *param, Ambertoppar *amb
 {
   initialize(simParams,param);
 
+  // This is AMBER file, so it is not a lonepairs PSF.
+  // Needs to be set FALSE to avoid crash.
+  is_lonepairs_psf = 0;
+
+  // General lonepairs flag must also be set FALSE to avoid crash.
+  // TIP4P lonepairs are still supported by setting config watmodel=tip4.
+  // The TIP4P lonepair repositioning is called from rigid bond constraints
+  // routine rattle1old().
+  simParams->lonepairs = 0;
+
   read_parm(amber_data);
 
 #ifndef MEM_OPT_VERSION
index 532cbe2..f3757b1 100644 (file)
@@ -2978,8 +2978,10 @@ void SimParameters::check_config(ParseOptions &opts, ConfigList *config, char *&
        watmodel = WAT_TIP4;
      } else if (!strncasecmp(s, "tip3", 4)) {
        iout << iINFO << "Using TIP3P water model.\n" << endi;
+       watmodel = WAT_TIP3;
      } else if (!strncasecmp(s, "swm4", 4)) {
        iout << iINFO << "Using SWM4-DP water model.\n" << endi;
+       watmodel = WAT_SWM4;
      } else {
        char err_msg[128];
        sprintf(err_msg,
@@ -3126,7 +3128,17 @@ void SimParameters::check_config(ParseOptions &opts, ConfigList *config, char *&
         NAMD_die(err_msg);
       }
    }
-   
+
+   // TIP4P and SWM4-DP water models require rigid water
+   if ((watmodel == WAT_TIP4 || watmodel == WAT_SWM4)
+       && rigidBonds == RIGID_NONE) {
+     char err_msg[256];
+     sprintf(err_msg,
+         "Water model %s requires rigidBonds set to \"all\" or \"water\"",
+         (watmodel == WAT_TIP4 ? "TIP4P" : "SWM4-DP"));
+     NAMD_die(err_msg);
+   }
+
    //  Take care of switching stuff
    if (switchingActive)
    {