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 [51]. 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 [62]. 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 [51].

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 [110,111]. 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 [43]. 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.

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.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.

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 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 energies resulting from the Drude oscillators and the anisotropic interactions are summed into the BOND energy. The energies resulting from the LPs and the screened Coulomb interactions of Thole are summed into the ELECT energy.

Drude force field parameters

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 [61], 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.

• drude $<$ Perform integration of Drude oscillators? $>$
  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.

• drudeTemp $<$ temperature for freezing the Drude oscillators (K) $>$
  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.

• drudeDamping $<$ damping coefficient for Drude oscillators (1/ps) $>$
  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.

• drudeNbTholeCut $<$ nonbonded Thole interaction radius (Å) $>$
  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.

• drudeHardWall $<$ use collisions to correct hyperpolarization? $>$
  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.

• drudeBondLen $<$ hyperpolarization cutoff (Å) $>$
  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.

• drudeBondConst $<$ Drude oscillator restraining force constant $>$
  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.