Update to lonepair and Drude particle bond behavior 55/5155/6
authorradakb <brian.radak@gmail.com>
Wed, 1 May 2019 20:09:06 +0000 (16:09 -0400)
committerDavid Hardy <dhardy@ks.uiuc.edu>
Thu, 30 May 2019 21:09:29 +0000 (16:09 -0500)
The CHARMM RTF/PSF standard for lonepairs no longer includes explicit bonds in
the BOND section. This would break NAMD, which needs explicit bond connectivity
to build inherited exclusions (and migration groups?). The development version
of Psfgen that builds equivalent PSFs also only writes Drude bonds as an option.

This update yields the following behavior:

1) When reading a PSF, Ignore all explicit BOND declarations that include
  lonepairs or Drude particles (flagged by a mass heuristic). This is a little
  dicey, but not actually that bad because NAMD then checks lonepair
  assignments against the lphost entries and Drude particles against non-zero
  polarizability -- a mismatch causes the program to stop.

2) If no wildecard DRUD-X bond is defined in the parameter files then, if
  CHARMM format Drude is on, a default bond is allocated (force constant of 500
  -- a global constant in Psfgen and all extant Drude RTFs).

3) Lonepair bonds are not "real" and now get assigned dummy parameters. This
  happens regardless of any definitions in the prm file (which usually don't
  exist and thus we never check for them).

4) After reading the lphost section, allocate new bonds between lonepairs and
  their primary host with dummy parameters. This fixes the missing connectivity
  problem used to build exclusions.

5) After reading a PSF, allocate new bonds between Drude particles and their
  preceding atom (which is already required to be their parent).

This has been extensively tested against Drude/lonepair containing systems
built with explicit bonds -- the old code and new code (w/o explicit bonds)
yield identical energies.

Possible issues/TODO:

1) Verify that MEMOPT builds are not effected.

Change-Id: I83467807a8cb71cb87fcd32c6a0ba3d09876f7e6

src/Molecule.C
src/Molecule.h
src/Parameters.C
src/Parameters.h

index 88248a3..e04f441 100644 (file)
@@ -320,8 +320,6 @@ void Molecule::initialize(SimParameters *simParams, Parameters *param)
   numAngles=0;
   numDihedrals=0;
   numImpropers=0;
-  numTholes=0;
-  numAnisos=0;
   numCrossterms=0;
   // JLai
   numLJPair=0;
@@ -333,8 +331,10 @@ void Molecule::initialize(SimParameters *simParams, Parameters *param)
   // DRUDE
   numLonepairs=0;
   numDrudeAtoms=0;
-  numLphosts=0;
+  numTholes=0;
   numAnisos=0;
+  numLphosts=0;
+  numZeroMassAtoms=0;
   // DRUDE
 
   numConstraints=0;
@@ -1273,6 +1273,27 @@ void Molecule::read_psf_file(char *fname, Parameters *params)
   /*  Close the .psf file.  */
   Fclose(psf_file);
 
+  if (is_drude_psf && numDrudeAtoms) {
+    // Automatically build Drude bonds that were previously ignored.
+    Bond *newbonds = new Bond[numBonds+numDrudeAtoms];
+    memcpy(newbonds, bonds, numBonds*sizeof(Bond));
+    delete [] bonds;
+    bonds = newbonds;
+    int origNumBonds = numBonds;
+    for (i=0; i < numAtoms; i++) {
+      if (!is_drude(i)) continue;
+      Bond *b = &(bonds[numBonds++]);
+      b->atom1 = i-1;
+      b->atom2 = i;
+      params->assign_bond_index(
+          atomNames[i-1].atomtype, atomNames[i].atomtype, b
+      );
+    }
+    if (numBonds-origNumBonds != numDrudeAtoms) {
+      NAMD_die("must have same number of Drude particles and parents");
+    }
+  }
+
   //  analyze the data and find the status of each atom
   numRealBonds = numBonds;
   build_atom_status();
@@ -1351,6 +1372,7 @@ void Molecule::read_atoms(FILE *fd, Parameters *params)
     read_count=sscanf(buffer, "%d %s %s %s %s %s %f %f",
        &atom_number, segment_name, residue_number,
        residue_name, atom_name, atom_type, &charge, &mass);
+    if (mass <= 0.05) ++numZeroMassAtoms;
 
     /*  Check to make sure we found what we were expecting  */
     if (read_count != 8)
@@ -1385,6 +1407,7 @@ void Molecule::read_atoms(FILE *fd, Parameters *params)
       }
       drudeConsts[atom_number-1].alpha = alpha;
       drudeConsts[atom_number-1].thole = thole;
+      if (fabs(alpha) >= 1e-6) ++numDrudeAtoms;
     }
     // DRUDE
 
@@ -1483,6 +1506,9 @@ void Molecule::read_bonds(FILE *fd, Parameters *params)
   register int j;      // Loop counter
   int num_read=0;    // Number of bonds read so far
   int origNumBonds = numBonds;   // number of bonds in file header
+  int numZeroFrcBonds = 0;
+  int numLPBonds = 0;
+  int numDrudeBonds = 0;
 
   /*  Allocate the array to hold the bonds      */
   bonds=new Bond[numBonds];
@@ -1529,23 +1555,34 @@ void Molecule::read_bonds(FILE *fd, Parameters *params)
     /*  Make sure this isn't a fake bond meant for shake in x-plor.  */
     Real k, x0;
     params->get_bond_params(&k,&x0,b->bond_type);
-    if (is_lonepairs_psf) {
-      // need to retain Lonepair bonds for Drude
-      if ( k == 0. && !is_lp(b->atom1) && !is_lp(b->atom2)) --numBonds;
-      else ++num_read;
-    }
-    else {
-      if ( k == 0. ) --numBonds;  // fake bond
-      else ++num_read;  // real bond
-    }
+
+    Bool is_lp_bond = (is_lp(b->atom1) || is_lp(b->atom2));
+    Bool is_drude_bond = (is_drude(b->atom1) || is_drude(b->atom2));
+    numZeroFrcBonds += (k == 0.);
+    numLPBonds += is_lp_bond;
+    numDrudeBonds += is_drude_bond;
+    if (k == 0. || is_lp_bond || is_drude_bond) --numBonds;  // fake bond
+    else ++num_read;  // real bond
   }
 
   /*  Tell user about our subterfuge  */
   if ( numBonds != origNumBonds ) {
-    iout << iWARN << "Ignored " << origNumBonds - numBonds <<
-            " bonds with zero force constants.\n" << endi;
-    iout << iWARN <<
-       "Will get H-H distance in rigid H2O from H-O-H angle.\n" << endi;
+    if (numZeroFrcBonds) {
+      iout << iWARN << "Ignored " << numZeroFrcBonds <<
+          " bonds with zero force constants.\n" <<
+          iWARN << "Will get H-H distance in rigid H2O from H-O-H angle.\n" <<
+          endi;
+    }
+    if (numLPBonds) {
+      iout << iWARN << "Ignored " << numLPBonds <<
+          " bonds with lone pairs.\n" <<
+          iWARN << "Will infer lonepair bonds from LPhost entries.\n" << endi;
+    }
+    if (numDrudeBonds) {
+      iout << iWARN << "Ignored " << numDrudeBonds <<
+          " bonds with Drude particles.\n" <<
+          iWARN << "Will use polarizability to assign Drude bonds.\n" << endi;
+    }
   }
 
   return;
@@ -2363,6 +2400,11 @@ void Molecule::read_lphosts(FILE *fd)
       lphosts[i].dihedral = value3 * (M_PI/180);  // convert to radians
     }
   }
+  // Resize bonds to accommodate the lonepairs.
+  Bond *newbonds = new Bond[numBonds+numLphosts];
+  memcpy(newbonds, bonds, numBonds*sizeof(Bond));
+  delete [] bonds;
+  bonds = newbonds;
   for (i = 0;  i < numLphosts;  i++) {
     // Subtract 1 to get 0-based atom index
     lphosts[i].atom1 = NAMD_read_int(fd, "LPHOSTS")-1;
@@ -2371,6 +2413,11 @@ void Molecule::read_lphosts(FILE *fd)
     // For numhosts==2, set unused atom4 to atom1
     lphosts[i].atom4 = ( lphosts[i].numhosts == 3 ?
         NAMD_read_int(fd, "LPHOSTS")-1 : lphosts[i].atom1 );
+    // Add dummy bond entry for connectivity.
+    Bond *b = &(bonds[numBonds++]);
+    b->atom1 = lphosts[i].atom1;
+    b->atom2 = lphosts[i].atom2;
+    b->bond_type = -1; // dummy index, never used
   }
 }
 /*      END OF FUNCTION read_lphosts    */
@@ -9170,31 +9217,22 @@ void Molecule::build_atom_status(void) {
   int a1, a2, a3;
   int numDrudeWaters = 0;
 
-  // if any atoms have a mass of zero set to 0.001 and warn user
-  int numZeroMassAtoms = 0;
-  for (i=0; i < numAtoms; i++) {
-    if ( atoms[i].mass <= 0. ) {
-      if (simParams->watmodel == WAT_TIP4 || is_lonepairs_psf) {
-        ++numLonepairs;
-      } else {
-        atoms[i].mass = 0.001;
-        ++numZeroMassAtoms;
-      }
-    }
-    else if (atoms[i].mass < 1.) {
-      ++numDrudeAtoms;
-    }
-  }
-  // DRUDE: verify number of LPs
-  if (is_lonepairs_psf && numLonepairs != numLphosts) {
-    NAMD_die("must have same number of LP hosts as lone pairs");
-  }
-  // DRUDE
-  if ( ! CkMyPe() ) {
-    if (simParams->watmodel == WAT_TIP4 || is_lonepairs_psf) {
+  if (simParams->watmodel == WAT_TIP4 || is_lonepairs_psf) {
+    numLonepairs = numZeroMassAtoms; // These MUST be lonepairs.
+    if ( ! CkMyPe() ) {
       iout << iWARN << "CORRECTION OF ZERO MASS ATOMS TURNED OFF "
         "BECAUSE LONE PAIRS ARE USED\n" << endi;
-    } else if ( numZeroMassAtoms ) {
+    }
+    // Compare the number of massless particles against the number of lonepair
+    // entries in the PSF -- these must match.
+    if (numLonepairs != numLphosts) {
+      NAMD_die("must have same number of LP hosts as lone pairs");
+    }
+  } else if (numZeroMassAtoms) {
+    for (i=0; i < numAtoms; i++) {
+      if ( atoms[i].mass < 0.001 ) atoms[i].mass = 0.001;
+    }
+    if ( ! CkMyPe() ) {
       iout << iWARN << "FOUND " << numZeroMassAtoms <<
         " ATOMS WITH ZERO OR NEGATIVE MASSES!  CHANGED TO 0.001\n" << endi;
     }
index 02d636e..e27af34 100644 (file)
@@ -569,11 +569,12 @@ public:
   int64 numTotalExclusions; //  Real Total Number of Exclusions // hack
 
   // DRUDE
-  int numLonepairs; // Number of lone pairs
-  int numDrudeAtoms;  // Number of Drude particles
-  int numTholes;  // Number of Thole terms
-  int numAnisos;  // Number of anisotropic terms
-  int numLphosts;  ///< Number of lone pair host records in PSF
+  int numLonepairs;     ///< Number of lone pairs
+  int numDrudeAtoms;    ///< Number of Drude particles
+  int numTholes;        ///< Number of Thole terms
+  int numAnisos;        ///< Number of anisotropic terms
+  int numLphosts;       ///< Number of lone pair host records in PSF
+  int numZeroMassAtoms; ///< Number of atoms with zero mass
   // DRUDE
   
   int numConstraints; //  Number of atoms constrained
index 1c29b75..8525387 100644 (file)
@@ -295,7 +295,7 @@ Parameters::Parameters(SimParameters *simParams, StringList *f)
       f = f->next;
     } while ( f != NULL );
 
-    done_reading_files();
+    done_reading_files(simParams->drudeOn && paramType == paraCharmm);
   }
 
 }
@@ -457,7 +457,7 @@ void Parameters::read_parameter_file(char *fname)
       /*  on the type of parameter we have    */
       if (strncasecmp(first_word, "bond", 4)==0)
       {
-        add_bond_param(buffer);
+        add_bond_param(buffer, TRUE);
         NumBondParams++;
       }
       else if (strncasecmp(first_word, "angl", 4)==0)
@@ -715,7 +715,7 @@ void Parameters::read_charmm_parameter_file(char *fname)
       /*  I know, this should really be a switch ...  */
       if (par_type == 1)
       {
-        add_bond_param(buffer);
+        add_bond_param(buffer, TRUE);
         NumBondParams++;
       }
       else if (par_type == 2)
@@ -839,7 +839,7 @@ void Parameters::skip_stream_read(char *buf, FILE *fd) {
 /*                  */
 /************************************************************************/
 
-void Parameters::add_bond_param(char *buf)
+void Parameters::add_bond_param(const char *buf, Bool overwrite)
 
 {
   char atom1name[11];    //  Atom type for atom 1
@@ -914,7 +914,7 @@ void Parameters::add_bond_param(char *buf)
 
   /*  Make call to recursive call to actually add the node to the */
   /*  tree              */
-  bondp=add_to_bond_tree(new_node, bondp);
+  bondp=add_to_bond_tree(new_node, bondp, overwrite);
 
   return;
 }
@@ -939,7 +939,7 @@ void Parameters::add_bond_param(char *buf)
 /************************************************************************/
 
 struct bond_params *Parameters::add_to_bond_tree(struct bond_params *new_node,
-             struct bond_params *tree)
+             struct bond_params *tree, Bool overwrite)
 
 {
   int compare_code;  //  Results from strcasecmp
@@ -962,13 +962,14 @@ struct bond_params *Parameters::add_to_bond_tree(struct bond_params *new_node,
     /*  If atom 1 AND atom 2 are the same, we have a duplicate */
     if (compare_code == 0)
     {
+
       /*  We have a duplicate.  So print out a warning*/
       /*  message.  Then assign the new values to the */
       /*  tree and free the new_node      */
       //****** BEGIN CHARMM/XPLOR type changes
       /* we do not care about identical replacement */
       if ((tree->forceconstant != new_node->forceconstant) || 
-          (tree->distance != new_node->distance))
+          (tree->distance != new_node->distance) && overwrite)
       {
         iout << "\n" << iWARN << "DUPLICATE BOND ENTRY FOR "
           << new_node->atom1name << "-"
@@ -995,11 +996,11 @@ struct bond_params *Parameters::add_to_bond_tree(struct bond_params *new_node,
   /*  otherwise add it to the right child        */
   if (compare_code < 0)
   {
-    tree->left = add_to_bond_tree(new_node, tree->left);
+    tree->left = add_to_bond_tree(new_node, tree->left, overwrite);
   }
   else
   {
-    tree->right = add_to_bond_tree(new_node, tree->right);
+    tree->right = add_to_bond_tree(new_node, tree->right, overwrite);
   }
 
   return(tree);
@@ -2955,11 +2956,16 @@ void Parameters::add_to_nbthole_pair_list(struct nbthole_pair_params *new_node)
 /*                  */
 /************************************************************************/
 
-void Parameters::done_reading_files()
+void Parameters::done_reading_files(Bool addDrudeBond)
 
 {
   AllFilesRead = TRUE;
 
+  if (addDrudeBond) {
+    // default definition for Drude bonds if none given
+    NumBondParams++;
+    add_bond_param("X DRUD 500.0 0.0\n", FALSE);
+  }
   //  Allocate space for all of the arrays
   if (NumBondParams)
   {
index 15f9408..fcf1ee2 100644 (file)
@@ -279,9 +279,9 @@ private:
 
        void skip_stream_read(char *, FILE *);  // skip part of stream file
 
-       void add_bond_param(char *);            //  Add a bond parameter
+       void add_bond_param(const char *, Bool);//  Add a bond parameter
        struct bond_params *add_to_bond_tree(struct bond_params * , 
-                                    struct bond_params *);
+                                    struct bond_params *, Bool overwrite);
 
        void add_angle_param(char *);           //  Add an angle parameter
        struct angle_params *add_to_angle_tree(struct angle_params * , 
@@ -400,7 +400,7 @@ public:
 
        //  Signal the parameter object that all of
        //  the parameter files have been read in
-       void done_reading_files();
+       void done_reading_files(Bool);
 
        //  Signal the parameter object that the
        //  structure file has been read in