From 229af5b048b360670f50726a0ff8fdc3e30e7252 Mon Sep 17 00:00:00 2001
From: radakb
Date: Thu, 28 Mar 2019 15:21:09 0400
Subject: [PATCH] Update default behavior and documentation for Drude
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 nonuniform 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.
ChangeId: Ief2d0a43f3e12af53867b3b9ddfe8a0be409abb5

src/ComputeBonds.C  5 +
src/HomePatch.C  2 +
src/Sequencer.C  22 ++
src/SimParameters.C  40 ++
ug/ug_forcefield.tex  296 ++++++++++++++++++++++
5 files changed, 181 insertions(+), 184 deletions()
diff git a/src/ComputeBonds.C b/src/ComputeBonds.C
index 86926343..0ab4a633 100644
 a/src/ComputeBonds.C
+++ b/src/ComputeBonds.C
@@ 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.
 if(simParams>drudeOn) { // for Drude bonds
+ if( simParams>drudeOn && !simParams>drudeHardWallOn ) { // for Drude bonds
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;
diff git a/src/HomePatch.C b/src/HomePatch.C
index 87ffd600..f794d6e3 100644
 a/src/HomePatch.C
+++ b/src/HomePatch.C
@@ 1929,7 +1929,7 @@ int HomePatch::hardWallDrude(const BigReal timestep, Tensor *virial,
r_wall_SQ = r_wall*r_wall;
// Count++;
for (i=1; i 0.01) && ((atom[i].mass < 1.0)) ) { // drude particle
+ if ( (0.05 < atom[i].mass) && ((atom[i].mass < 1.0)) ) { // drude particle
ia = i1;
ib = i;
diff git a/src/Sequencer.C b/src/Sequencer.C
index 67f06190..cbe94ac7 100644
 a/src/Sequencer.C
+++ b/src/Sequencer.C
@@ 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];
 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[i1].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 ) {
 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[i1].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;
}
@@ 806,7 +806,7 @@ void Sequencer::newMinimizePosition(BigReal c) {
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[i1].position;
}
}
@@ 816,7 +816,7 @@ void Sequencer::newMinimizePosition(BigReal c) {
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[i1].position;
}
}
@@ 1208,7 +1208,7 @@ void Sequencer::langevinVelocitiesBBK1(BigReal dt_fs)
for (i = 0; i < numAtoms; i++) {
if (i < numAtoms1 &&
 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
@@ 1295,7 +1295,7 @@ void Sequencer::langevinVelocitiesBBK2(BigReal dt_fs)
for (i = 0; i < numAtoms; i++) {
if (i < numAtoms1 &&
 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
@@ 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() );
 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;
}
@@ 2105,7 +2105,7 @@ void Sequencer::submitReductions(int step)
for (i = 0; i < numAtoms; i++) {
if (i < numAtoms1 &&
 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
@@ 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
 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();
@@ 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;
 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 ) {
diff git a/src/SimParameters.C b/src/SimParameters.C
index d1457871..1b3360f5 100644
 a/src/SimParameters.C
+++ b/src/SimParameters.C
@@ 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", "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 "
 "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 "
 "force constant", &drudeBondConst);
+ "force constant", &drudeBondConst, 40000.0);
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",
@@ 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;
}
 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
@@ 5669,18 +5654,17 @@ if ( openatomOn )
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";
}
 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 "
diff git a/ug/ug_forcefield.tex b/ug/ug_forcefield.tex
index e5dd11e9..373d4f38 100644
 a/ug/ug_forcefield.tex
+++ b/ug/ug_forcefield.tex
@@ 865,133 +865,119 @@ Alternative orderings will fail. }
\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 particleparticle 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\'eHoover thermostat acting on, and within,
each nucleusDrude pair~\cite{Lamoureux2003a}.
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 outofplane 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 nonexcluded, 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 offcentered massless interaction sites,
so called ``lone pairs'' (LPs),
to avoid the limitations
of centrosymmetricbased 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 inplane and outofplane 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
+ zerolength harmonic spring.
+The advantage with the Drude model is that it preserves the simple
+ particleparticle 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\'eHoover
+ thermostat acting on, and within, each
+ nucleusDrude pair~\cite{Lamoureux2003a}.
+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 flatbottom
+ 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 outofplane 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 nonexcluded, 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 offcentered massless
+ interaction sites, so called ``lone pairs'' (LPs), to avoid the limitations
+ of centrosymmetricbased 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 inplane and outofplane 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}
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 Drudecompatible 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 Drudecompatible PSF file is denoted by the keyword ``DRUDE'' given along the
+ top line.
\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}
\item[{\tt DRUDECOM}] gives the temperature for the warm centerofmass
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 centerofmass 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 centerofmass
+ 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 centerofmass 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}
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
({\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.
Longrange electrostatics may be calculated using PME.
The nonbonded exclusions should generally be set to use either the
13 exclusion policy ({\tt exclude=13})
or the scaled 14 exclusion policy ({\tt exclude=scaled14}).

The Drude water model (SWM4NDP) is a 5site model
with four charge sites and
a negatively charged Drude particle~\cite{Lamoureux2006a},
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.
+Longrange electrostatics may be calculated using PME.
+The nonbonded exclusions should generally be set to use either the 13
+ exclusion policy ({\tt exclude=13}) or the scaled 14 exclusion policy
+ ({\tt exclude=scaled14}).
+
+The Drude water model (SWM4NDP) is a 5site model with four charge sites and
+ a negatively charged Drude particle~\cite{Lamoureux2006a}, 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}
@@ 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}}
{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)}
{nonnegative 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}
{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}
{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 nonzero, the screened Coulomb correction of Thole
+ is also calculated for nonexcluded, 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}
{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}
{If {\tt drudeNbTholeCut} is defined,
the screened Coulomb correction of Thole is also calculated
for nonexcluded, 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 nonzero, an additional
+ quartic restraining potential is applied to a Drude oscillator if its length
+ exceeds {\tt drudeBondLen}.
+}
\end{itemize}

2.20.1