Drude polarizable force field

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 [30]. Use of the Langevin thermostat enables better parallel scalability than the earlier reported implementation which made use of a dual Nosé-Hoover thermostat acting on, and within, each nucleus-Drude pair [38]. 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 [30].

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 [75,76].
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 [26].
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 (`http://www.ks.uiuc.edu/Research/Drude/`) for
additional details and parallel performance results.

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.

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 `ETITLE` and `ENERGY`
lines:

`DRUDECOM`- gives the temperature for the warm center-of-mass degrees of freedom,
`DRUDEBOND`- gives the temperature for the cold Drude oscillator degrees of freedom,
`DRCOMAVG`- gives the average temperature (averaged since the previously reported temperature) for the warm center-of-mass degrees of freedom,
`DRBONDAVG`- gives the average temperature (averaged since the previously reported temperature) for the cold Drude oscillator degrees of freedom.

The Drude model should be used with the Langevin thermostat enabled
(`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 (`exclude=1-3`) or the scaled 1-4 exclusion policy
(`exclude=scaled1-4`).

The Drude water model (SWM4-NDP) is a 5-site model with four charge sites and
a negatively charged Drude particle [37], with the
particles ordered in the input files as oxygen, Drude particle, LP, hydrogen,
hydrogen.
The atoms in the water molecules should be constrained
(`rigidBonds=water`), with use of the SETTLE algorithm recommended
(`useSettle=on`).
Explicitly setting the water model (`waterModel=swm4`) is optional.

Perform integration of Drude oscillators?`drude`**Acceptable Values:**`on`or`off`**Default Value:**`off`**Description:**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`drude`is set to`on`, then`drudeTemp`must also be defined.temperature for freezing the Drude oscillators (K)`drudeTemp`**Acceptable Values:**non-negative decimal**Description:**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.damping coefficient for Drude oscillators (1/ps)`drudeDamping`**Acceptable Values:**positive decimal**Description:**The Langevin coupling coefficient to be applied to the Drude oscillators. If not given,`drudeDamping`is set to the value of`langevinDamping`, but values of as much as an order of magnitude greater may be appropriate.nonbonded Thole interaction radius (Å)`drudeNbTholeCut`**Acceptable Values:**positive decimal**Default Value:**5.0**Description:**If`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.use collisions to correct hyperpolarization?`drudeHardWall`**Acceptable Values:**on or off**Default Value:**on**Description:**Excessively elongated Drude oscillator bonds are avoided by reflective collisions induced at a fixed cutoff,`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.hyperpolarization cutoff (Å)`drudeBondLen`**Acceptable Values:**positive decimal**Default Value:**0.25**Description:**If using`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 Å (slightly smaller than default) is recommended.Drude oscillator restraining force constant`drudeBondConst`**Acceptable Values:**positive decimal**Default Value:**40000.0**Description:**If`drudeHardWall off`and`drudeBondConst`is non-zero, an additional quartic restraining potential is applied to a Drude oscillator if its length exceeds`drudeBondLen`.