Preliminary implementation of colinear lonepairs 22/4522/4
authorradakb <brian.radak@gmail.com>
Fri, 24 Aug 2018 18:20:09 +0000 (14:20 -0400)
committerDavid Hardy <dhardy@ks.uiuc.edu>
Fri, 31 Aug 2018 17:41:43 +0000 (12:41 -0500)
The CHARMM Drude and additive force fields now support colinear as
well as relative and bisector lonepairs. The lonepair facility is
now extended to enable PSFs with these specifications. Other lone
pair types exist in CHARMM but only as part of the RTF standard and
these are still NOT supported, since no extant force field makes use
of those lonepair types.

Summary of changes:

o The lphosts struct now contains a numhosts attribute (i.e. the
  number of particles needed to specify the position of the lonepair).
  This was previously assumed to be three, but can now also be two.

o The read_lphosts() routine now stores the number of hosts for
  each lonepair and can read variable numbers of indices in the
  index section (this previously assumed four indices per lonepair).

o The reposition and redistribute routines inside HomePatch now
  have multiple code paths depending on the number of lphosts.
  These could probably be better identified and the if/else tree
  likely needs refinement in order to permit other types of
  lonepairs.

o Preliminary doxygen comments for the various routines mentioned
  above.

NB: All lonepair specifications in a PSF are expected to have three
    associated values. However, the meaning of these values depends
    on the lonepair type. Nonetheless, for simplicity with the old
    code, the struct values are still called "distance", "angle",
    and "dihedral", even when this is not what the value signifies.

NB: The "weight" specification in the PSF remains hardcoded as "F"
    (for false) inside NAMD. Again, no extant force field uses
    lonepairs with weighted placement, so this can really only be
    specified inside an RTF, which NAMD never reads anyway.

Change-Id: I4751d597ba4d43b078c1da052e56f8a5252481a4

src/HomePatch.C
src/HomePatch.h
src/Molecule.C
src/Molecule.h
src/structures.h

index 16c46a8..152b4d3 100644 (file)
@@ -1274,7 +1274,37 @@ void HomePatch::saveForce(const int ftag)
 #undef DEBUG_REDISTRIB_FORCE 
 #undef DEBUG_REDISTRIB_FORCE_VERBOSE
 //#define DEBUG_REDISTRIB_FORCE
-/*
+
+/**
+ * Redistribute forces from lone pair site onto its colinear host atoms.
+ *
+ * /param fi the force vector on lonepair i
+ * /param fj the force vector on the first host j
+ * /param fk the force vector on the second host k
+ * /param ri the position vector of the lone pair i
+ * /param rj the position vector of the first host j
+ * /param rk the position vector of the second host k
+ * /param distance the distance to be placed from the reference point along the
+ *     vector from k to j
+ * /param scale the reference point is placed using a scaled value of the
+ *     vector from k to j.
+ */
+void HomePatch::redistrib_colinear_lp_force(
+    Vector& fi, Vector& fj, Vector& fk,
+    const Vector& ri, const Vector& rj, const Vector& rk,
+    Real distance_f, Real scale_f) {
+  BigReal distance = distance_f;
+  BigReal scale = scale_f;
+  Vector r_jk = rj - rk;
+  // TODO: Add a check for small distances?
+  BigReal r_jk_rlength = r_jk.rlength();
+  distance *= r_jk_rlength;
+  BigReal fdot = distance*(fi*r_jk)*r_jk_rlength*r_jk_rlength;
+  fj += (1. + scale + distance)*fi - r_jk*fdot;
+  fk -= (scale + distance)*fi + r_jk*fdot;
+}
+
+/**
  * Redistribute forces from lone pair site onto its host atoms.
  *
  * Atoms are "labeled" i, j, k, l, where atom i is the lone pair.
@@ -1300,7 +1330,7 @@ void HomePatch::saveForce(const int ftag)
  * the cross(r,s) zero.
  */
 #define FIX_FOR_WATER
-void HomePatch::redistrib_lp_force(
+void HomePatch::redistrib_relative_lp_force(
     Vector& fi, Vector& fj, Vector& fk, Vector& fl,
     const Vector& ri, const Vector& rj, const Vector& rk, const Vector& rl,
     Tensor *virial, int midpt) {
@@ -1535,7 +1565,54 @@ void HomePatch::redistrib_lp_water_force(
 #endif /* DEBUG */
 }
 
-void HomePatch::reposition_lonepair(
+/**
+ * Reposition lonepair i in a colinear fashion relative to its hosts j and k.
+ *
+ * The lonepair is placed at a fixed (possibly negative) distance along the
+ * vector from k to j relative to a reference point. The reference point is
+ * computed by multiplying the vector from k to j by a (possibly negative)
+ * scale factor. For example, a scale value of 0 (the default) sets the
+ * reference point as rj, while a value of -1 sets the reference point as rk.
+ * A scale value of -0.5 would use the center of the two hosts.
+ *
+ * /param ri the position vector of the lonepair, i
+ * /param rj the position vector of the first host, j
+ * /param rk the position vector of the second host, k
+ * /param distance the distance to be placed from the reference point along the
+ *     vector from k to j
+ * /param scale the reference point is placed using a scaled value of the
+ *     vector from k to j.
+ */
+void HomePatch::reposition_colinear_lonepair(
+    Vector& ri, const Vector& rj, const Vector& rk,
+    Real distance_f, Real scale_f)
+{
+  BigReal distance = distance_f;
+  BigReal scale = scale_f;
+  Vector r_jk = rj - rk;
+  BigReal r2 = r_jk.length2();
+  if (r2 < 1e-10 || 100. < r2) { // same low tolerance as used in CHARMM
+    iout << iWARN << "Large/small distance between lonepair reference atoms: ("
+         << rj << ") (" << rk << ")\n" << endi;
+  }
+  ri = rj + (scale + distance*r_jk.rlength())*r_jk;
+}
+
+/**
+ * Reposition lonepair i relative to its hosts j, k, and l.
+ *
+ * The lonepair is placed in a bond-angle-torsion coordinate frame. An angle
+ * and torsion of 0 corresponds to a bisector placement, as in TIP4P water.
+ *
+ * /param ri the position vector of the lonepair, i
+ * /param rj the position vector of the first host, j
+ * /param rk the position vector of the second host, k
+ * /param rl the position vector of the third host, j
+ * /param distance the negative distance from j
+ * /param angle the angle made by i-j-k
+ * /param dihedral the dihedral made by i-j-k-l
+ */
+void HomePatch::reposition_relative_lonepair(
     Vector& ri, const Vector& rj, const Vector& rk, const Vector& rl,
     Real distance, Real angle, Real dihedral)
 {
@@ -1600,9 +1677,15 @@ void HomePatch::reposition_all_lonepairs(void) {
         NAMD_die(errmsg);
       }
       // reposition this lone pair
-      reposition_lonepair(atom[i].position, atom[j.index].position,
-          atom[k.index].position, atom[l.index].position,
-          lph->distance, lph->angle, lph->dihedral);
+      if (lph->numhosts == 2) {
+        reposition_colinear_lonepair(atom[i].position, atom[j.index].position,
+            atom[k.index].position, lph->distance, lph->angle);
+      }
+      else if (lph->numhosts == 3) {
+        reposition_relative_lonepair(atom[i].position, atom[j.index].position,
+            atom[k.index].position, atom[l.index].position,
+            lph->distance, lph->angle, lph->dihedral);
+      }
     }
   }
 }
@@ -1670,11 +1753,18 @@ void HomePatch::redistrib_lonepair_forces(const int ftag, Tensor *virial) {
         NAMD_die(errmsg);
       }
       // redistribute forces from this lone pair
-      int midpt = (lph->distance < 0);
-      redistrib_lp_force(f_mod[ftag][i], f_mod[ftag][j.index],
-          f_mod[ftag][k.index], f_mod[ftag][l.index],
-          atom[i].position, atom[j.index].position,
-          atom[k.index].position, atom[l.index].position, virial, midpt);
+      if (lph->numhosts == 2) {
+        redistrib_colinear_lp_force(f_mod[ftag][i], f_mod[ftag][j.index],
+            f_mod[ftag][k.index], atom[i].position, atom[j.index].position,
+            atom[k.index].position, lph->distance, lph->angle);
+      }
+      else if (lph->numhosts == 3) {
+        int midpt = (lph->distance < 0);
+        redistrib_relative_lp_force(f_mod[ftag][i], f_mod[ftag][j.index],
+            f_mod[ftag][k.index], f_mod[ftag][l.index],
+            atom[i].position, atom[j.index].position,
+            atom[k.index].position, atom[l.index].position, virial, midpt);
+      }
     }
   }
 }
index 53865d0..373a8b5 100644 (file)
@@ -396,7 +396,10 @@ private:
   BigReal settle_mOrmT; BigReal settle_mHrmT; BigReal settle_ra;
   BigReal settle_rb; BigReal settle_rc; BigReal settle_rra;
 
-  // Drude lone pairs
+  /**
+   * Redistribute all lonepair forces (of any kind). This may include a direct
+   * correction to the virial.
+   */
   void redistrib_lonepair_forces(const int, Tensor *);
 
   // PLF -- for TIP4P
@@ -410,20 +413,49 @@ private:
   void swm4_omrepos(Vector*, Vector*, Vector*, BigReal);
   void init_swm4();
 
-  // reposition a lone pair using its host atoms and additional parameters
-  void reposition_lonepair(
+  /**
+   * Reposition lonepair i in a colinear fashion relative to its hosts j and k
+   * and according to a fixed distance and scaled vector magnitude between the
+   * two hosts.
+   */
+  void reposition_colinear_lonepair(
+      Vector& ri, const Vector& rj, const Vector& rk, Real distance,
+      Real scale);
+
+  /**
+   * Reposition a lonepair i relative to its hosts j, k, and l according to a
+   * given distance, angle, and dihedral formed with the three hosts.
+   */
+  void reposition_relative_lonepair(
       Vector& ri, const Vector& rj, const Vector& rk, const Vector& rl,
       Real distance, Real angle, Real dihedral);
 
+  /**
+   * Reposition all lonepairs (of any kind).
+   */
   void reposition_all_lonepairs(void);
 
-  // general redistribution of lone pair forces to host atoms
-  void redistrib_lp_force(
+  /**
+   * Redistribute the force on a colinear lonepair onto its hosts.
+   */
+  void redistrib_colinear_lp_force(
+       Vector& fi, Vector& fj, Vector& fk,
+       const Vector& ri, const Vector& rj, const Vector& rk,
+       Real distance, Real scale);
+
+  /**
+   * Redistribute the force on a relative lonepair onto its hosts.
+   */
+  void redistrib_relative_lp_force(
       Vector& fi, Vector& fj, Vector& fk, Vector& fl,
       const Vector& ri, const Vector& rj, const Vector& rk, const Vector& rl,
       Tensor *virial, int midpt);
 
-  // use for both TIP4P and SWM4 water
+  /**
+   * Redistribute the force on a water (TIP4P, SWM4) lonepair onto its hosts.
+   * This is similar to redistrib_relative_lp_force but specialized for the
+   * bisector case.
+   */
   void redistrib_lp_water_force(
       Vector& f_ox, Vector& f_h1, Vector& f_h2, Vector& f_lp,
       const Vector& p_ox, const Vector& p_h1, const Vector& p_h2,
index 96069c3..d8d13ea 100644 (file)
@@ -2354,50 +2354,79 @@ void Molecule::read_exclusions(int* atom_i, int* atom_j, int num_exclusion)
 }
 /*      END OF FUNCTION read_exclusions      */
 
-/************************************************************************/
-/*                  */
-/*        FUNCTION read_lphosts    */
-/*                  */
-/*   INPUTS:                */
-/*  fd - file pointer to the .psf file        */
-/*                  */
-/*  this function reads in the lone pair host section of the .psf file. */
-/*                  */
+/************************************************************************
+ *
+ *        FUNCTION read_lphosts
+ *
+ *   INPUTS:
+ *  fd - file pointer to the .psf file
+ *
+ *  This function reads in the lone pair host section of the .psf file.
+ *  Each lone pair is supported by 2 or 3 host atoms.
+ *
+ *  All lonepair specifications in a PSF are expected to have three
+ *  associated values.  However, the meaning of these values depends
+ *  on the lonepair type.  Nonetheless, for simplicity with the old
+ *  code, the struct values are still called "distance", "angle",
+ *  and "dihedral", even when this is not what the value signifies.
+ *
+ ************************************************************************/
 void Molecule::read_lphosts(FILE *fd)
 {
   char buffer[512];  // Buffer for reading from file
-  char lptype[8];
-  int numhosts, index, i, read_count;
-  Real distance, angle, dihedral;
-
+  char weight[8];    // Weighting type identified by string --
+                     // unused/unsupported outside of CHARMM RTF
+                     // so we read it, but ignore it.
+  int numhosts;  // Refers to the number of atoms that support the
+                 // given lone pair, either 2 or 3.
+  int index;  // 1-based index into the stream of numbers (8 per line)
+              // that indicates the start of each LP host record.
+  int next_index = 1;  // Forecast next index value as an error check.
+  int i, read_count;
+  Real value1, value2, value3; 
+
+  // called only if numLphosts > 0
   lphosts = new Lphost[numLphosts];
-  if (lphosts == NULL)
-  {
+  if (lphosts == NULL) {
     NAMD_die("memory allocation failed in Molecule::read_lphosts");
   }
-  for (i = 0;  i < numLphosts;  i++)
-  {
+  for (i = 0;  i < numLphosts;  i++) {
     NAMD_read_line(fd, buffer);
     if ( (NAMD_blank_string(buffer)) || (buffer[0] == '!') ) continue;
     read_count=sscanf(buffer, "%d %d %6s %f %f %f",
-        &numhosts, &index, lptype, &distance, &angle, &dihedral);
-    if (read_count != 6 || numhosts != 3 || index != 4*i + 1
-        || strcmp(lptype,"F") != 0)
-    {
+        &numhosts, &index, weight, &value1, &value2, &value3);
+    // The "weight" specification in PSF remains hardcoded as "F"
+    // (for false) inside NAMD.  Again, no extant force field uses
+    // lonepairs with weighted placement, so this can really only be
+    // specified inside an RTF, which NAMD never reads anyway.
+    if (read_count != 6 || index != next_index ||
+        strcmp(weight,"F") != 0 || numhosts < 2 || 3 < numhosts) {
       char err_msg[128];
       sprintf(err_msg, "BAD FORMAT FOR LPHOST LINE %d IN PSF FILE LINE\n"
           "LINE=%s\n", i+1, buffer);
       NAMD_die(err_msg);
     }
-    lphosts[i].distance = distance;
-    lphosts[i].angle = angle * (M_PI/180);        // convert to radians
-    lphosts[i].dihedral = dihedral * (M_PI/180);  // convert to radians
+    lphosts[i].numhosts = numhosts; // currently must be 2 or 3
+    next_index += numhosts + 1;  // add 1 to account for LP index
+    if (numhosts == 2) {
+      lphosts[i].distance = value1;
+      lphosts[i].angle = value2;
+      lphosts[i].dihedral = 0.0;  // ignore value3
+    }
+    else {  // numhosts == 3
+      lphosts[i].distance = value1;
+      lphosts[i].angle = value2 * (M_PI/180);     // convert to radians
+      lphosts[i].dihedral = value3 * (M_PI/180);  // convert to radians
+    }
   }
   for (i = 0;  i < numLphosts;  i++) {
+    // Subtract 1 to get 0-based atom index
     lphosts[i].atom1 = NAMD_read_int(fd, "LPHOSTS")-1;
     lphosts[i].atom2 = NAMD_read_int(fd, "LPHOSTS")-1;
     lphosts[i].atom3 = NAMD_read_int(fd, "LPHOSTS")-1;
-    lphosts[i].atom4 = NAMD_read_int(fd, "LPHOSTS")-1;
+    // For numhosts==2, set unused atom4 to atom1
+    lphosts[i].atom4 = ( lphosts[i].numhosts == 3 ?
+        NAMD_read_int(fd, "LPHOSTS")-1 : lphosts[i].atom1 );
   }
 }
 /*      END OF FUNCTION read_lphosts    */
@@ -9409,6 +9438,8 @@ void Molecule::build_atom_status(void) {
   }
 
   // determine migration groups based on lone pair hosts
+  // Find the lowest atom index in each migration group for all
+  // lone pair support atoms.  This value marks the migration group.
   for (i=0; i<numLphosts; ++i) {
     int a1 = lphosts[i].atom1;
     int a2 = lphosts[i].atom2;
index dfed468..8a8fb47 100644 (file)
@@ -573,7 +573,7 @@ public:
   int numDrudeAtoms;  // Number of Drude particles
   int numTholes;  // Number of Thole terms
   int numAnisos;  // Number of anisotropic terms
-  int numLphosts;  // Number of lone pair hosts
+  int numLphosts;  ///< Number of lone pair host records in PSF
   // DRUDE
   
   int numConstraints; //  Number of atoms constrained
index 1bdccbe..7f636cd 100644 (file)
@@ -106,12 +106,22 @@ typedef struct drude_constants  // supplement Atom data
   Real thole;
 } DrudeConst;
 
-typedef struct lphost    // lone pair host
+/** Lonepair host record
+ *
+ * Maintains record of LP index and supporting atom indices.
+ * Field numhosts is either 2 or 3, depending on LP type.
+ * For the case of 2 supporting atoms, index atom4 is set to repeat atom1.
+ * The lonepair parameters make geometric sense only in 3-numhosts case.
+ * For the 2-numhosts case, "distance" and "angle" both represent
+ * distances and set "dihedral" to zero.
+ */
+typedef struct lphost
 {
-  int32 atom1;
+  int32 atom1;     ///< First index is to the LP
   int32 atom2;
   int32 atom3;
   int32 atom4;
+  int32 numhosts;  ///< Either 2 or 3 host atoms, depending on LP type
   Real distance;
   Real angle;
   Real dihedral;