Update default behavior and documentation for Drude 43/5043/2
authorradakb <brian.radak@gmail.com>
Thu, 28 Mar 2019 19:21:09 +0000 (15:21 -0400)
committerDavid <dhardy@ks.uiuc.edu>
Fri, 29 Mar 2019 22:08:26 +0000 (17:08 -0500)
Changes to default behavior:

- drudeHardWall is now documented in the user guide and is the
  default option.
- drudeBondLen, drudeBondConst, and drudeNbTholeCut now have default
  options equal to the recommended values from the user guide. These
  are still "optional" values but are essentially required for any
  sensible simulation.
- drudeBondConst is ignored unless drudeHardWall off (not default).

Cosmetic updates:

- INFO message is now upated to avoid incorrect duplicate message for
  quartic potential (mutually exclusive with drudeHardWall).
- Cleaned up the user guide to improve description of hyperpolarization
  events and options for avoiding/correcting them.
- Corrected non-uniform heuristics for Drude particles based on mass.
  A lower mass cutoff of 0.001, 0.01, or 0.1 was used in multiple
  places. These now match the value of 0.05 used alongside other mass
  heuristics in Molecule.
  TODO - maybe these heuristics should be centralized?
- Made a small update in the quartic bonded potential to avoid a second
  if condition when drudeHardWall on.

Change-Id: Ief2d0a43f3e12af53867b3b9ddfe8a0be409abb5

src/ComputeBonds.C
src/HomePatch.C
src/Sequencer.C
src/SimParameters.C
ug/ug_forcefield.tex

index 8692634..0ab4a63 100644 (file)
@@ -89,11 +89,10 @@ void BondElem::computeForce(BondElem *tuples, int ntuple, BigReal *reduction,
     // TODO: Instead flag by mass for Drude particles, otherwise mixing and 
     // matching Drude and alchemical "twin" particles will likely not work as 
     // expected.
     // TODO: Instead flag by mass for Drude particles, otherwise mixing and 
     // matching Drude and alchemical "twin" particles will likely not work as 
     // expected.
-    if(simParams->drudeOn) { // for Drude bonds
+    if( simParams->drudeOn && !simParams->drudeHardWallOn ) { // for Drude bonds
       BigReal drudeBondLen = simParams->drudeBondLen;
       BigReal drudeBondConst = simParams->drudeBondConst;
       BigReal drudeBondLen = simParams->drudeBondLen;
       BigReal drudeBondConst = simParams->drudeBondConst;
-      if ( (drudeBondConst > 0) && ( ! simParams->drudeHardWallOn ) &&
-          (r2 > drudeBondLen*drudeBondLen) ) {
+      if ( drudeBondConst > 0 && r2 > drudeBondLen*drudeBondLen ) {
         // add a quartic restraining potential to keep Drude bond short
         BigReal r = sqrt(r2);
         BigReal diff = r - drudeBondLen;
         // add a quartic restraining potential to keep Drude bond short
         BigReal r = sqrt(r2);
         BigReal diff = r - drudeBondLen;
index 87ffd60..f794d6e 100644 (file)
@@ -1929,7 +1929,7 @@ int HomePatch::hardWallDrude(const BigReal timestep, Tensor *virial,
     r_wall_SQ = r_wall*r_wall;
     // Count++;
     for (i=1; i<numAtoms; i++) {
     r_wall_SQ = r_wall*r_wall;
     // Count++;
     for (i=1; i<numAtoms; i++) {
-      if ( (atom[i].mass > 0.01) && ((atom[i].mass < 1.0)) ) { // drude particle
+      if ( (0.05 < atom[i].mass) && ((atom[i].mass < 1.0)) ) { // drude particle
         ia = i-1;
         ib = i;
 
         ia = i-1;
         ib = i;
 
index 67f0619..cbe94ac 100644 (file)
@@ -748,7 +748,7 @@ void Sequencer::newMinimizeDirection(BigReal c) {
   for ( int i = 0; i < numAtoms; ++i ) {
     a[i].velocity *= c;
     a[i].velocity += f1[i];
   for ( int i = 0; i < numAtoms; ++i ) {
     a[i].velocity *= c;
     a[i].velocity += f1[i];
-    if ( drudeHardWallOn && i && (a[i].mass > 0.01) && ((a[i].mass < 1.0)) ) { // drude particle
+    if ( drudeHardWallOn && i && (0.05 < a[i].mass) && ((a[i].mass < 1.0)) ) { // drude particle
       a[i].velocity = a[i-1].velocity;
     }
     if ( fixedAtomsOn && a[i].atomFixed ) a[i].velocity = 0;
       a[i].velocity = a[i-1].velocity;
     }
     if ( fixedAtomsOn && a[i].atomFixed ) a[i].velocity = 0;
@@ -760,10 +760,10 @@ void Sequencer::newMinimizeDirection(BigReal c) {
 
   maxv2 = 0.;
   for ( int i = 0; i < numAtoms; ++i ) {
 
   maxv2 = 0.;
   for ( int i = 0; i < numAtoms; ++i ) {
-    if ( drudeHardWallOn && i && (a[i].mass > 0.01) && ((a[i].mass < 1.0)) ) { // drude particle
+    if ( drudeHardWallOn && i && (0.05 < a[i].mass) && ((a[i].mass < 1.0)) ) { // drude particle
       a[i].velocity = a[i-1].velocity;
       a[i].velocity = a[i-1].velocity;
-      if ( fixedAtomsOn && a[i].atomFixed ) a[i].velocity = 0;
     }
     }
+    if ( fixedAtomsOn && a[i].atomFixed ) a[i].velocity = 0;
     BigReal v2 = a[i].velocity.length2();
     if ( v2 > maxv2 ) maxv2 = v2;
   }
     BigReal v2 = a[i].velocity.length2();
     if ( v2 > maxv2 ) maxv2 = v2;
   }
@@ -806,7 +806,7 @@ void Sequencer::newMinimizePosition(BigReal c) {
 
   if ( simParams->drudeHardWallOn ) {
     for ( int i = 1; i < numAtoms; ++i ) {
 
   if ( simParams->drudeHardWallOn ) {
     for ( int i = 1; i < numAtoms; ++i ) {
-      if ( (a[i].mass > 0.01) && ((a[i].mass < 1.0)) ) { // drude particle
+      if ( (0.05 < a[i].mass) && ((a[i].mass < 1.0)) ) { // drude particle
         a[i].position -= a[i-1].position;
       }
     }
         a[i].position -= a[i-1].position;
       }
     }
@@ -816,7 +816,7 @@ void Sequencer::newMinimizePosition(BigReal c) {
 
   if ( simParams->drudeHardWallOn ) {
     for ( int i = 1; i < numAtoms; ++i ) {
 
   if ( simParams->drudeHardWallOn ) {
     for ( int i = 1; i < numAtoms; ++i ) {
-      if ( (a[i].mass > 0.01) && ((a[i].mass < 1.0)) ) { // drude particle
+      if ( (0.05 < a[i].mass) && ((a[i].mass < 1.0)) ) { // drude particle
         a[i].position += a[i-1].position;
       }
     }
         a[i].position += a[i-1].position;
       }
     }
@@ -1208,7 +1208,7 @@ void Sequencer::langevinVelocitiesBBK1(BigReal dt_fs)
       for (i = 0;  i < numAtoms;  i++) {
 
         if (i < numAtoms-1 &&
       for (i = 0;  i < numAtoms;  i++) {
 
         if (i < numAtoms-1 &&
-            a[i+1].mass < 1.0 && a[i+1].mass >= 0.001) {
+            a[i+1].mass < 1.0 && a[i+1].mass > 0.05) {
           //printf("*** Found Drude particle %d\n", a[i+1].id);
           // i+1 is a Drude particle with parent i
 
           //printf("*** Found Drude particle %d\n", a[i+1].id);
           // i+1 is a Drude particle with parent i
 
@@ -1295,7 +1295,7 @@ void Sequencer::langevinVelocitiesBBK2(BigReal dt_fs)
       for (i = 0;  i < numAtoms;  i++) {
 
         if (i < numAtoms-1 &&
       for (i = 0;  i < numAtoms;  i++) {
 
         if (i < numAtoms-1 &&
-            a[i+1].mass < 1.0 && a[i+1].mass >= 0.001) {
+            a[i+1].mass < 1.0 && a[i+1].mass > 0.05) {
           //printf("*** Found Drude particle %d\n", a[i+1].id);
           // i+1 is a Drude particle with parent i
 
           //printf("*** Found Drude particle %d\n", a[i+1].id);
           // i+1 is a Drude particle with parent i
 
@@ -1637,7 +1637,7 @@ void Sequencer::reinitVelocities(void)
           a[i].mass <= 0.) ? Vector(0,0,0) :
         sqrt(kbT * (a[i].partition ? tempFactor : 1.0) * a[i].recipMass) * 
         random->gaussian_vector() );
           a[i].mass <= 0.) ? Vector(0,0,0) :
         sqrt(kbT * (a[i].partition ? tempFactor : 1.0) * a[i].recipMass) * 
         random->gaussian_vector() );
-    if ( simParams->drudeOn && i+1 < numAtoms && a[i+1].mass < 1.0 && a[i+1].mass >= 0.001 ) {
+    if ( simParams->drudeOn && i+1 < numAtoms && a[i+1].mass < 1.0 && a[i+1].mass > 0.05 ) {
       a[i+1].velocity = a[i].velocity;  // zero is good enough
       ++i;
     }
       a[i+1].velocity = a[i].velocity;  // zero is good enough
       ++i;
     }
@@ -2105,7 +2105,7 @@ void Sequencer::submitReductions(int step)
 
         for (i = 0;  i < numAtoms;  i++) {
           if (i < numAtoms-1 &&
 
         for (i = 0;  i < numAtoms;  i++) {
           if (i < numAtoms-1 &&
-              a[i+1].mass < 1.0 && a[i+1].mass >= 0.001) {
+              a[i+1].mass < 1.0 && a[i+1].mass > 0.05) {
             // i+1 is a Drude particle with parent i
 
             // convert from Cartesian coordinates to (COM,bond) coordinates
             // i+1 is a Drude particle with parent i
 
             // convert from Cartesian coordinates to (COM,bond) coordinates
@@ -2319,7 +2319,7 @@ void Sequencer::submitMinimizeReductions(int step, BigReal fmax2)
 
   for ( int i = 0; i < numAtoms; ++i ) {
     f1[i] += f2[i] + f3[i];  // add all forces
 
   for ( int i = 0; i < numAtoms; ++i ) {
     f1[i] += f2[i] + f3[i];  // add all forces
-    if ( drudeHardWallOn && i && (a[i].mass > 0.01) && ((a[i].mass < 1.0)) ) { // drude particle
+    if ( drudeHardWallOn && i && (a[i].mass > 0.05) && ((a[i].mass < 1.0)) ) { // drude particle
       if ( ! fixedAtomsOn || ! a[i].atomFixed ) {
         if ( drudeStep2 * f1[i].length2() > drudeMove2 ) {
           a[i].position += drudeMove * f1[i].unit();
       if ( ! fixedAtomsOn || ! a[i].atomFixed ) {
         if ( drudeStep2 * f1[i].length2() > drudeMove2 ) {
           a[i].position += drudeMove * f1[i].unit();
@@ -2361,7 +2361,7 @@ void Sequencer::submitMinimizeReductions(int step, BigReal fmax2)
   int numHuge = 0;
   for ( int i = 0; i < numAtoms; ++i ) {
     if ( simParams->fixedAtomsOn && a[i].atomFixed ) continue;
   int numHuge = 0;
   for ( int i = 0; i < numAtoms; ++i ) {
     if ( simParams->fixedAtomsOn && a[i].atomFixed ) continue;
-    if ( drudeHardWallOn && (a[i].mass > 0.01) && ((a[i].mass < 1.0)) ) continue; // drude particle
+    if ( drudeHardWallOn && (a[i].mass > 0.05) && ((a[i].mass < 1.0)) ) continue; // drude particle
     Force f = f1[i];
     BigReal ff = f * f;
     if ( ff > fmax2 ) {
     Force f = f1[i];
     BigReal ff = f * f;
     if ( ff > fmax2 ) {
index d145787..1b3360f 100644 (file)
@@ -1224,17 +1224,17 @@ void SimParameters::config_parser_methods(ParseOptions &opts) {
    opts.optional("drude", "drudeDamping", "Damping coefficient (1/ps) for "
        "Drude oscillators", &drudeDamping);
    opts.range("drudeDamping", POSITIVE);
    opts.optional("drude", "drudeDamping", "Damping coefficient (1/ps) for "
        "Drude oscillators", &drudeDamping);
    opts.range("drudeDamping", POSITIVE);
+   opts.optional("drude", "drudeNbtholeCut", "Nonbonded Thole interactions "
+       "interaction radius", &drudeNbtholeCut, 5.0);
+   opts.range("drudeNbtholeCut", POSITIVE);
+   opts.optionalB("drude", "drudeHardWall", "Apply maximum Drude bond length "
+       "restriction?", &drudeHardWallOn, TRUE);
    opts.optional("drude", "drudeBondLen", "Drude oscillator bond length "
    opts.optional("drude", "drudeBondLen", "Drude oscillator bond length "
-       "beyond which to apply restraint", &drudeBondLen);
+       "beyond which to apply restraint", &drudeBondLen, 0.25);
    opts.range("drudeBondLen", POSITIVE);
    opts.optional("drude", "drudeBondConst", "Drude oscillator restraining "
    opts.range("drudeBondLen", POSITIVE);
    opts.optional("drude", "drudeBondConst", "Drude oscillator restraining "
-       "force constant", &drudeBondConst);
+       "force constant", &drudeBondConst, 40000.0);
    opts.range("drudeBondConst", POSITIVE);
    opts.range("drudeBondConst", POSITIVE);
-   opts.optional("drude", "drudeNbtholeCut", "Nonbonded Thole interactions "
-       "interaction radius", &drudeNbtholeCut);
-   opts.range("drudeNbtholeCut", POSITIVE);
-   opts.optionalB("drude", "drudeHardWall", "Apply maximum Drude bond length "
-       "restriction?", &drudeHardWallOn, FALSE);
 
    // Pair interaction calculations
     opts.optionalB("main", "pairInteraction", 
 
    // Pair interaction calculations
     opts.optionalB("main", "pairInteraction", 
@@ -3678,21 +3678,6 @@ void SimParameters::check_config(ParseOptions &opts, ConfigList *config, char *&
        iout << iWARN << "Undefined 'drudeDamping' will be set to "
          "value of 'langevinDamping'\n" << endi;
      }
        iout << iWARN << "Undefined 'drudeDamping' will be set to "
          "value of 'langevinDamping'\n" << endi;
      }
-     if ( ! opts.defined("drudeBondConst")) {
-       drudeBondConst = 0;
-       if (! drudeHardWallOn) {
-         drudeBondLen = 0;
-         if (opts.defined("drudeBondLen")) {
-           iout << iWARN << "Resetting 'drudeBondLen' to 0 "
-             "since 'drudeBondConst' and 'drudeHardWall' are unset\n" << endi;
-         }
-       }
-     }
-     if ( ! opts.defined("drudeNbtholeCut")) {
-       drudeNbtholeCut = 0;
-         iout << iWARN << "Resetting 'drudeNbtholeCut' to 0 "
-           "since 'drudeNbtholeCut' is unset\n" << endi;
-     }
    }
 
    //  Set up load balancing variables
    }
 
    //  Set up load balancing variables
@@ -5669,18 +5654,17 @@ if ( openatomOn )
         iout << iINFO << "DRUDE DAMPING COEFFICIENT IS "
              << drudeDamping << " INVERSE PS\n";
       }
         iout << iINFO << "DRUDE DAMPING COEFFICIENT IS "
              << drudeDamping << " INVERSE PS\n";
       }
-      if (drudeBondConst > 0.0) {
+      if (drudeHardWallOn) {
+        iout << iINFO << "DRUDE HARD WALL RESTRAINT IS ACTIVE FOR DRUDE BONDS\n";
+        iout << iINFO << "DRUDE MAXIMUM BOND LENGTH BEFORE RESTRAINT IS   "
+             << drudeBondLen << "\n";
+      } else if (drudeBondConst > 0.0) {
         iout << iINFO << "DRUDE QUARTIC RESTRAINT IS ACTIVE FOR DRUDE BONDS\n";
         iout << iINFO << "DRUDE MAXIMUM BOND LENGTH BEFORE RESTRAINT IS   "
              << drudeBondLen << "\n";
         iout << iINFO << "DRUDE BOND RESTRAINT CONSTANT IS                "
              << drudeBondConst << "\n";
       }
         iout << iINFO << "DRUDE QUARTIC RESTRAINT IS ACTIVE FOR DRUDE BONDS\n";
         iout << iINFO << "DRUDE MAXIMUM BOND LENGTH BEFORE RESTRAINT IS   "
              << drudeBondLen << "\n";
         iout << iINFO << "DRUDE BOND RESTRAINT CONSTANT IS                "
              << drudeBondConst << "\n";
       }
-      if (drudeHardWallOn) {
-        iout << iINFO << "DRUDE HARD WALL RESTRAINT IS ACTIVE FOR DRUDE BONDS\n";
-        iout << iINFO << "DRUDE MAXIMUM BOND LENGTH BEFORE RESTRAINT IS   "
-             << drudeBondLen << "\n";
-      }
       if (drudeNbtholeCut > 0.0) {
         iout << iINFO << "DRUDE NBTHOLE IS ACTIVE\n";
         iout << iINFO << "DRUDE NBTHOLE RADIUS IS   "
       if (drudeNbtholeCut > 0.0) {
         iout << iINFO << "DRUDE NBTHOLE IS ACTIVE\n";
         iout << iINFO << "DRUDE NBTHOLE RADIUS IS   "
index e5dd11e..373d4f3 100644 (file)
@@ -865,133 +865,119 @@ Alternative orderings will fail. }
 \subsection{Drude polarizable force field}
 \label{section:drude_forcefield}
 
 \subsection{Drude polarizable force field}
 \label{section:drude_forcefield}
 
-The Drude oscillator model represents induced electronic polarization 
-by introducing an auxiliary particle attached to each polarizable 
-atom via a harmonic spring.  The advantage with the Drude model is 
-that it preserves the simple particle-particle Coulomb electrostatic 
-interaction employed in nonpolarizable force fields, 
-thus its implementation in \NAMD\ is more straightforward than 
-alternative models for polarization.  
-\NAMD\ performs the integration of Drude oscillators 
-by employing a novel dual Langevin thermostat to freeze the Drude 
-oscillators while maintaining the warm degrees of freedom 
-at the desired temperature~\cite{JIAN2011}.  
-Use of the Langevin thermostat enables better parallel scalability 
-than the earlier reported implementation which made use of 
-a dual Nos\'e-Hoover thermostat acting on, and within,
-each nucleus-Drude pair~\cite{Lamoureux-2003a}.  
-Performance results 
-show that the \NAMD\ implementation of the Drude model maintains good 
-parallel scalability, 
-with an increase in computational cost by not more than twice 
-that of using a nonpolarizable force field~\cite{JIAN2011}.  
-
-The Drude polarizable force field requires some extensions to the
-CHARMM force field.  
-The Drude oscillators differ from typical spring bonds only in that they 
-have an equilibrium length of zero.  
-The Drude oscillators are optionally supplemented by a maximal 
-bond length parameter, beyond which a quartic restraining potential 
-is also applied.  
-The force field is also extended by an anisotropic spring term 
-that accounts for out-of-plane forces from a polarized atom and its 
-covalently bonded neighbor with two more covalently bonded 
-neighbors (similar in structure to an improper bond).  
-The screened Coulomb correction of Thole is calculated between 
-pairs of Drude oscillators that would otherwise be excluded from 
-nonbonded interaction and 
-optionally between non-excluded, nonbonded pairs of Drude oscillators 
-that are within a prescribed cutoff distance~\cite{Thole81,van1998molecular}.  
-
-Also included in the Drude force field 
-are the use of off-centered massless interaction sites, 
-so called ``lone pairs'' (LPs),
-to avoid the limitations 
-of centrosymmetric-based Coulomb interactions~\cite{Harder2006}.  
-The coordinate of each LP site is constructed based on three ``guide'' atoms.
-The calculated forces on the massless LP must be transferred 
-to the guide atoms, 
-preserving total force and torque.  
-After an integration step of velocities and positions, 
-the position of the LP is updated based on the three guide atoms, 
-along with additional geometry parameters that give displacement 
-and in-plane and out-of-plane angles.  
-See our research web page
-(\url{http://www.ks.uiuc.edu/Research/Drude/})
-for additional details and parallel performance results.  
-
+The Drude oscillator model represents induced electronic polarization by
+  introducing an auxiliary particle attached to each polarizable atom via a
+  zero-length harmonic spring.
+The advantage with the Drude model is that it preserves the simple
+  particle-particle Coulomb electrostatic interaction employed in
+  nonpolarizable force fields, thus its implementation in \NAMD\ is more
+  straightforward than alternative models for polarization.
+\NAMD\ performs the integration of Drude oscillators by employing a novel dual
+  Langevin thermostat to ``freeze'' the Drude oscillators while maintaining the
+  normal ``warm'' degrees of freedom at the desired
+  temperature~\cite{JIAN2011}.
+Use of the Langevin thermostat enables better parallel scalability than the
+  earlier reported implementation which made use of a dual Nos\'e-Hoover
+  thermostat acting on, and within, each
+  nucleus-Drude pair~\cite{Lamoureux-2003a}.
+Performance results show that the \NAMD\ implementation of the Drude model
+  maintains good parallel scalability with an increase in computational cost by
+  not more than twice that of using a nonpolarizable force
+  field~\cite{JIAN2011}.
+
+Excessive ``hyperpolarization'' of Drude oscillators can be prevented by two
+  different schemes.
+The default ``hard wall'' option reflects elongated springs back towards the
+  nucleus using a simple collision model.
+Alternatively, the Drude oscillators can be supplemented by a flat-bottom
+  quartic restraining potential (usually with a large force constant).
+
+The Drude polarizable force field requires some extensions to the CHARMM force
+  field.
+An anisotropic spring term is added to account for out-of-plane forces from a
+  polarized atom and its covalently bonded neighbor with two more covalently
+  bonded neighbors (similar in structure to an improper bond).
+The screened Coulomb correction of Thole is calculated between pairs of Drude
+  oscillators that would otherwise be excluded from nonbonded interaction and
+  optionally between non-excluded, nonbonded pairs of Drude oscillators
+  that are within a prescribed cutoff distance~\cite{Thole81,van1998molecular}.
+Also included in the Drude force field are the use of off-centered massless
+  interaction sites, so called ``lone pairs'' (LPs), to avoid the limitations
+  of centrosymmetric-based Coulomb interactions~\cite{Harder2006}.
+The coordinate of each LP site is constructed based on three ``host'' atoms.
+The calculated forces on the massless LP must be transferred to the host atoms,
+  preserving total force and torque.
+After an integration step of velocities and positions, the position of the LP
+  is updated based on the three host atoms, along with additional geometry
+  parameters that give displacement and in-plane and out-of-plane angles.
+See our research web page (\url{http://www.ks.uiuc.edu/Research/Drude/}) for
+  additional details and parallel performance results.
 
 \subsubsection{Required input files}
 
 
 \subsubsection{Required input files}
 
-No additional files are required by \NAMD\ to use the 
-Drude polarizable force field.  
-However, it is presently beyond the capability of 
-the Psfgen tool to generate the PSF file needed to perform 
-a simulation using the Drude model.  
-For now, CHARMM is needed to generate correct input files.  
-
-The CHARMM force field parameter files
-specific to the Drude model are required.  
-The PDB file must also include 
-the Drude particles (mass between 0.1 and 1.0) 
-and the LPs (mass 0).  
-The Drude particles always immediately follow their parent atom.  
-The PSF file augments the ``atom'' section with 
-additional columns that include the 
-``Thole'' and ``alpha'' parameters for the screened Coulomb 
-interactions of Thole.  
-The PSF file also requires additional sections that 
-list the LPs, including their guide atoms and geometry parameters, 
-and list the anisotropic interaction terms, including their parameters.  
-A Drude-compatible PSF file is denoted by the 
-keyword ``DRUDE'' given along the top line.  
-
+No additional files are required by \NAMD\ to use the Drude polarizable force
+  field.
+However, it is presently beyond the capability of the \PSFGEN{} tool to
+  generate the PSF file needed to perform a simulation using the Drude model.
+For now, CHARMM is needed to generate correct input files.
+
+The CHARMM force field parameter files specific to the Drude model are
+  required.
+The PDB file must also include the Drude particles (mass between 0.05 and 1.0)
+  and the LPs (mass 0).
+The Drude particles always immediately follow their parent atom.
+The PSF file augments the ``atom'' section with additional columns that include
+  the ``Thole'' and ``alpha'' parameters for the screened Coulomb interactions
+  of Thole.
+The PSF file also requires additional sections that list the LPs, including
+  their host atoms and geometry parameters, and list the anisotropic
+  interaction terms, including their parameters.
+A Drude-compatible PSF file is denoted by the keyword ``DRUDE'' given along the
+  top line.
 
 \subsubsection{Standard output}
 
 
 \subsubsection{Standard output}
 
-The \NAMD\ logging to standard output is extended to provide additional 
-temperature data on the cold and warm degrees of freedom.  
-Four additional quantities are listed 
-on the {\tt ETITLE} and {\tt ENERGY} lines:
+The \NAMD\ logging to standard output is extended to provide additional
+  temperature data on the cold and warm degrees of freedom.
+Four additional quantities are listed on the {\tt ETITLE} and {\tt ENERGY}
+  lines:
 \begin{description}
 \begin{description}
-\item[{\tt DRUDECOM}] gives the temperature for the warm center-of-mass 
-degrees of freedom,
-\item[{\tt DRUDEBOND}] gives the temperature for the cold Drude oscillator 
-degrees of freedom,
-\item[{\tt DRCOMAVG}] gives the average temperature 
-(averaged since the previously reported temperature) 
-for the warm center-of-mass degrees of freedom, 
-\item[{\tt DRBONDAVG}] gives the average temperature 
-(averaged since the previously reported temperature) 
-for the cold Drude oscillator degrees of freedom.
+\item[{\tt DRUDECOM}] gives the temperature for the warm center-of-mass
+  degrees of freedom,
+\item[{\tt DRUDEBOND}] gives the temperature for the cold Drude oscillator
+  degrees of freedom,
+\item[{\tt DRCOMAVG}] gives the average temperature (averaged since the
+  previously reported temperature) for the warm center-of-mass degrees of
+  freedom,
+\item[{\tt DRBONDAVG}] gives the average temperature (averaged since the
+  previously reported temperature) for the cold Drude oscillator degrees of
+  freedom.
 \end{description}
 \end{description}
-The energies resulting from the Drude oscillators and the 
-anisotropic interactions are summed into the {\tt BOND} energy.  
-The energies resulting from the LPs and the screened Coulomb 
-interactions of Thole are summed into the {\tt ELECT} energy.  
-
+The energies resulting from the Drude oscillators and the anisotropic
+  interactions are summed into the {\tt BOND} energy.
+The energies resulting from the LPs and the screened Coulomb interactions of
+  Thole are summed into the {\tt ELECT} energy.
 
 \subsubsection{Drude force field parameters}
 
 The Drude model should be used with the Langevin thermostat enabled
 
 \subsubsection{Drude force field parameters}
 
 The Drude model should be used with the Langevin thermostat enabled
-({\tt Langevin=on}).  
-Doing so permits the use of normal sized time steps (e.g., 1 fs).  
-The Drude model is also compatible with constant pressure simulation 
-using the Langevin piston.  
-Long-range electrostatics may be calculated using PME.  
-The nonbonded exclusions should generally be set to use either the
-1-3 exclusion policy ({\tt exclude=1-3})
-or the scaled 1-4 exclusion policy ({\tt exclude=scaled1-4}).  
-
-The Drude water model (SWM4-NDP) is a 5-site model
-with four charge sites and 
-a negatively charged Drude particle~\cite{Lamoureux-2006a}, 
-with the particles ordered in the input files as 
-oxygen, Drude particle, LP, hydrogen, hydrogen.  
-The atoms in the water molecules should be 
-constrained ({\tt rigidBonds=water}),
-with use of the SETTLE algorithm recommended ({\tt useSettle=on}).  
-Explicitly setting the water model ({\tt waterModel=swm4}) is optional.  
+  ({\tt Langevin=on}).
+Doing so permits the use of normal sized time steps (e.g., 1 fs).
+The Drude model is also compatible with constant pressure simulation using the
+  Langevin piston.
+Long-range electrostatics may be calculated using PME.
+The nonbonded exclusions should generally be set to use either the 1-3
+  exclusion policy ({\tt exclude=1-3}) or the scaled 1-4 exclusion policy
+  ({\tt exclude=scaled1-4}).
+
+The Drude water model (SWM4-NDP) is a 5-site model with four charge sites and
+  a negatively charged Drude particle~\cite{Lamoureux-2006a}, with the
+  particles ordered in the input files as oxygen, Drude particle, LP, hydrogen,
+  hydrogen.
+The atoms in the water molecules should be constrained
+  ({\tt rigidBonds=water}), with use of the SETTLE algorithm recommended
+  ({\tt useSettle=on}).
+Explicitly setting the water model ({\tt waterModel=swm4}) is optional.
 
 \begin{itemize}
 
 
 \begin{itemize}
 
@@ -999,46 +985,74 @@ Explicitly setting the water model ({\tt waterModel=swm4}) is optional.
 {Perform integration of Drude oscillators?}
 {{\tt on} or {\tt off}}
 {{\tt off}}
 {Perform integration of Drude oscillators?}
 {{\tt on} or {\tt off}}
 {{\tt off}}
-{The integration uses a dual Langevin thermostat to freeze the Drude 
-oscillators while maintaining the warm degrees of freedom at the
-desired temperature.  Must also enable the Langevin thermostat.
-If {\tt drude} is set to {\tt on}, then {\tt drudeTemp} must also be defined.}
+{
+The integration uses a dual Langevin thermostat to freeze the Drude
+  oscillators while maintaining the warm degrees of freedom at the desired
+  temperature.
+Must also enable the Langevin thermostat.
+If {\tt drude} is set to {\tt on}, then {\tt drudeTemp} must also be defined.
+}
 
 \item\NAMDCONF{drudeTemp}
 {temperature for freezing the Drude oscillators (K)}
 {non-negative decimal}
 
 \item\NAMDCONF{drudeTemp}
 {temperature for freezing the Drude oscillators (K)}
 {non-negative decimal}
-{For stability, the Drude oscillators must be kept at a very cold termpature.
-Using a Langevin thermostat, it is possible to set this temperature to 0 K.}
+{
+For stability, the Drude oscillators must be kept at a very cold termpature.
+Using a Langevin thermostat, it is possible to set this temperature to 0 K.
+}
 
 \item\NAMDCONF{drudeDamping}
 {damping coefficient for Drude oscillators (1/ps)}
 {positive decimal}
 
 \item\NAMDCONF{drudeDamping}
 {damping coefficient for Drude oscillators (1/ps)}
 {positive decimal}
-{The Langevin coupling coefficient to be applied to the Drude oscillators.  
-If not given, {\tt drudeDamping} is set to the value of {\tt langevinDamping}.}
+{
+The Langevin coupling coefficient to be applied to the Drude oscillators.
+If not given, {\tt drudeDamping} is set to the value of {\tt langevinDamping},
+  but values of as much as an order of magnitude greater may be appropriate.
+}
 
 
-\item\NAMDCONF{drudeBondLen}
-{Drude oscillator bond length, beyond which to apply restraint (\AA)}
+\item\NAMDCONFWDEF{drudeNbTholeCut}
+{nonbonded Thole interaction radius (\AA)}
 {positive decimal}
 {positive decimal}
-{An additional quartic restraining potential is applied to a Drude 
-oscillator if its length exceeds {\tt drudeBondLen}.  
-The recommended value is 0.2 \AA, fitted from QM calculations.}
+{5.0}
+{
+If {\tt drudeNbTholeCut} is non-zero, the screened Coulomb correction of Thole
+  is also calculated for non-excluded, nonbonded pairs of Drude oscillators
+  that are within this radius of interaction.
+}
 
 
-\item\NAMDCONF{drudeBondConst}
-{Drude oscillator restraining force constant}
+\item\NAMDCONFWDEF{drudeHardWall}
+{use collisions to correct hyperpolarization?}
+{on or off}
+{on}
+{
+Excessively elongated Drude oscillator bonds are avoided by reflective
+  collisions induced at a fixed cutoff, {\tt drudeBondLen}.
+A large number of such events is usually indicative of unstable/unphysical
+  dynamics and a simulation will stop if double the cutoff is exceeded.
+}
+
+\item\NAMDCONFWDEF{drudeBondLen}
+{hyperpolarization cutoff (\AA)}
 {positive decimal}
 {positive decimal}
-{If {\tt drudeBondConst} is defined,
-an additional quartic restraining potential is applied to a Drude 
-oscillator if its length exceeds {\tt drudeBondLen}.
-The recommended value is 40000, fitted from QM calculations.}
+{0.25}
+{
+If using {\tt drudeHardWall on}, this is the distance at which collisions
+  occur.
+Otherwise, this is the distance at which an additional quartic restraining
+  potential is applied to each Drude oscillator.
+In this latter case, a value of 0.2~\AA{} (slightly smaller than default) is
+  recommended.
+}
 
 
-\item\NAMDCONF{drudeNbTholeCut}
-{nonbonded Thole interaction radius (\AA)}
+\item\NAMDCONFWDEF{drudeBondConst}
+{Drude oscillator restraining force constant}
 {positive decimal}
 {positive decimal}
-{If {\tt drudeNbTholeCut} is defined, 
-the screened Coulomb correction of Thole is also calculated
-for non-excluded, nonbonded pairs of Drude oscillators that are
-within this radius of interaction.
-The recommended value is 5.0 \AA.}
+{40000.0}
+{
+If {\tt drudeHardWall off} and {\tt drudeBondConst} is non-zero, an additional
+  quartic restraining potential is applied to a Drude oscillator if its length
+  exceeds {\tt drudeBondLen}.
+}
 
 \end{itemize}
 
 
 \end{itemize}