drude (c44b1)
Polarizable Drude Oscillator Format
by Benoit Roux (roux@uchicago.edu)
and Guillaume Lamoureux (Guillaume.Lamoureux@umontreal.ca)
and Alex MacKerell Jr. (alex@outerbanks.umaryland.edu)
Classical drude particles are generated for all atoms for which
polarizabilities (ie. via the ALPHA keyword) are specified in the RTF
file. This allows for the drudes to be generated automatically when a
molecule is generated in CHARMM (only the heavy atoms carry a
polarizability in the current force field but the code is general and
allows assigning a polarizability to any atom including hydrogens). In
addition, code has been developed to allow for inclusion of atom-based
Thole scale factors, atom-based anisotropic polarizabilities and the
addition of lone pairs to selected atoms at the RTF level. These
capabilities allow for all the information for the polarizable Drude
force field to be included in the RTF and parameter files. This
implementation replaces the old "Drude Oscillator Command" and
"ANISOTROPY keyword" (see below) used to generate the Drude
polarizable model prior to July 2007.
For each atom with a specified polarizability, a "Drude oscillator" is
created by attaching to the atom an additional particle (using a
fictitious chemical bond of length zero and of force constant 'KDRUDE
= k/2'). Each Drude particle is given a mass and a charge, taken from
the mass and the charge of its atom (so that the total mass and charge
are conserved for the "atom-Drude" pair).
As a whole, each "atom-Drude" pair carries a total charge 'q',
unchanged From the partial charge the non-polarizable atom had prior
to calling the DRUDE command during system generation. The
"atom-Drude" pair forms a dipole 'q*d', where 'q' is the charge on the
Drude particle and 'd' is the displacement vector going from the atom
to its Drude particle. Any external electrostatic field 'E' creates a
net displacement 'd = q*E/k', and thus the "atom-Drude" pair behaves
as a point charge 'q' with a polarizability 'alpha = q**2/k'. The
polarizabilities (in Angst**3) are specified with the ALPHA command in
the RTF, along with the THOLE factors (see below), and converted into
charges 'q', assuming a force constant 'k = 2*KDRUDE', where KDRUDE is
specified in the parameter file (see below); a value of 500 should
always be used.
See J. Chem. Phys. 119, 3025-3039 (2003) for more details.
In the Drude treatment of polarization the force constant of the
spring between the atom-Drude pair is a diagonal rank-2 tensor with
components KDRUDE. This leads to an isotropic atomic polarizability,
'alpha = alpha_{11} = alpha_{22} = alpha_{33} = q**2/k'. The
ANISOTROPY command, which is typically included in the RTF, modifies
the components of the Drude force constant tensor allowing for
anisotropic atomic polarizabilites.
See J. Chem. Theory Comput. 2, 1587-1597 (2006) for more details.
The bonded lists are modified so that, if the "real" atoms are in a
1-2, 1-3, or 1-4 relationship, their corresponding Drude particles
will also be in a 1-2, 1-3, or 1-4 relationship, respectively.
For a single atom (charges in parenthesis):
DRUDE
A --------> A~DA
For a diatomic molecule:
A1 DRUDE A1~DA1
| --------> |
A2 A2~DA2
1-2 pairs: A1-A2, A1-DA2, DA1-A2, DA1-DA2
For a triatomic molecule:
A1 A1~DA1
\ DRUDE \
A2 --------> A2~DA2
/ /
A3 A3~DA3
1-2 pairs: A1-A2, A1-DA2, DA1-A2, DA1-DA2,
A2-A3, A2-DA3, DA2-A3, DA2-DA3
1-3 pairs: A1-A3, A1-DA3, DA1-A3, DA1-DA3
Thus, the bonded interactions are modified to allow 1-2 and 1-3
screened dipole-dipole interactions, as proposed by Thole [see
Chem.Phys. 59, 341 (1981)]. If two atoms were not "seeing" each others
in the non-polarizable force field, their dipoles (and only their
dipoles, not their net partial charges) will "see" each others in the
polarizable force field. The screening function strength is modulated
by a parameter "a" that has been generalized to depend on the atom
type associated with the interacting pair of Drude oscillators and are
then generated via a standard arithmetic combination rule (a = a_i +
a_j). These variable atomic THOLE parameters "a_i" are specified in
the RTF through the THOLE keyword.
* Description | Description of new DRUDE model RTF
Description of new DRUDE model RTF
1) Items in RTF file
A) The AUTOGENERATE command now includes a flag for the drude
particles. This should be specified at the beginning of the RTF. Note
that the use of the PATCh keyword assures that autogenerate is applied
after all patches. This should be used.
AUTO ANGLES DIHE DRUDE PATCH
B) MASS statements must be included for Drude particle types.
MASS 104 DRUD 0.00000 DD ! drude particle
The default drude type is DRUD. However, alternate drude types may be
specified. These additional drude types must be i) specified with MASS
statements, ii) the ATOMs to which they are applied must include a
"TYPE drudetype" statement and iii) the NONBond parameters must be
specified. For example, to include a special drude type on water the
following would be needed in the RTF
i) MASS 153 DOH2 0.00000 DD ! water Drude
ii) ATOM OH2 ODW 0.00000 ALPHA -0.97825258 TYPE DOH2
iii) (NONBOND section) D* 0.0 -0.0000 0.0000
! Wildcard for Drudes and dummy atom
iv) (via NBFIX) MAGD DOH2 -0.01000 3.10000
Term iv) is an NBFIX that allows for the interactions between, for
example, magnesium and the water drude to be treated with a special LJ
term. Indeed, the motivation for additional drude types are such
special terms may be applied in special cases. Note that NBFIX terms
have been used between select real atoms in the Drude force field to
improve the free energies of hydration.
See J. Chem. Theory Comput. 6. 1181-1198 (2010) for more details.
C) ALPHA and THOLE parameters are set in the RTF after the charge.
The presence of the "ALPHA" keyword flags that atom as being
polarizable and a drude particle is created and assigned to it when
calling the GENERATE command. An atom based Thole scale factor may
also be specified via the "THOLE" keyword. If the THOLE term is
missing, a default value of 1.3 is assigned to that atom. This
corresponds to a parameter a = a_i + a_j = 2.6 which is the THOLE
parameter in the old Drude command syntax.
ATOM C C 0.630 ALPHA -1.104 THOLE 1.073
Note that a negative sign on alpha leads to the charge on the Drude particle
to be negative. This convention was chosen because the Drudes are meant to
represent electronic degrees of freedom of the system.
D) Specific drude types for a polarizable atom may be specified as
follows using the "TYPE" keyword. See 1A) above, sections i to iv for
more details.
ATOM C C 0.630 ALPHA -1.104 TYPE DC
E) Lone pairs may now be set directly in the RTF using a syntax
similar to the standard lonepair command » lonepair . All options
for the generation of lonepairs (FIXed, CENTer, COLOcate etc.) as
specified in the lonepair documentation may be used although the atom
selection specification is simplified as shown in the following
example.
LONEPAIR relative LPA O C CL distance 0.3 angle 91.0 dihe 0.0
LONEPAIR relative LPB O C N distance 0.3 angle 91.0 dihe 0.0
F) Atom based anisotropic polarizabilities may be assigned to selected
atoms via the ANISOTROPY keyword in the RTF. The anisotropic
polarizability tensor is defined as ALPHA_iso (the isotropic
polarizability) times a tensor A whith trace{A}=3. The tensor A is
diagonal in the local reference frame and is completely specified by
{A11, A22, A33}. The first atom selects the atom on which the
polarizability is to be made anisotropic. The second atom along with
the first defines the "11" direction of the {A} tensor and the 3rd and
4th atom selections define the "22" direction of the tensor, with the
anisotropy defined by fractional polarizability components.
Accordingly, only the 11 and 22 direction are needed since the trace
of the {A} tensor is set to 3. An example follows
ANISOTROPY O C CA N A11 0.6968 A22 1.2194
In the code, the variables A11 and A22 are read into the variables
A11ANIS(IPT) and A22ANIS(IPT). We define the components of the tensor
relative to the isotropic drude, so ALPHA11+ALPHA22+ALPHA33 =
3*ALPHA_iso. The factors A11 and A22 (and A33) represent the factors
modifying the components of the isotropic polarizability tensor
q**2/KDRUDE. For example, setting A11=A22=1 would keep the Drude
completely isotropic. By virtue of the trace we know that
A11+A22+A33=3, so we only need to set two of these terms. Thus A33 =
THREE-A11ANIS(IPT)-A22ANIS(IPT). This sets the relative magnitude of
the polarizability tensor in all 3 directions.
Once the polarizability along a given direction is known, we need to
get the harmonic force constant in this direction. As the
polarizability in a given direction is defined as
QDRUDE**2/K_direction, the 3 harmonic force constants are thus
K11(NANISO) = ONE/A11ANIS(IPT)
K22(NANISO) = ONE/A22ANIS(IPT)
K33(NANISO) = ONE/A33
Now, in principle this would mean that you have to calculate the
forces from 3 different springs (K11, K22, K33) in 3 different
directions. But since we know the overall trace which is related to
the isotropic polarizability, we are going to re-write this so that
there are only 2 springs with force constants KPAR in the parallel and
KPERP in the perpendicular direction, and add an isotropic KDRUDE
spring on top of that.
2) Items in the Parameter file
A) The KDRUDE force constant for all the atoms is set by a wildcard in
the bond parameters. The term must be included for all drude types
included in the model. The wildcard terms can be overwritten by
putting chemical type specific bond parameters following the wild card
term. Note that 500.00 is the default force constant and should not be
changed.
DRUD X 500.000 0.0000
DRUD O 487.740 0.0000
B) NONBOND parameters must be included for all drude types. In
addition, NBFIX terms for the drudes may be included as needed (see 1A
above).
C) A pair-specific Thole shielding can be introduced in the parameter
file for any nonbond pair (following a syntax similar to NBFIX) using
the keyword THOLe. This feature was introduced to permit accurate
models of divalent cations with polarizable atoms. The TCUT assigns a
radius within which all Thole interaction are calculated between the
specificed atoms. The 'aij' typically is assigned by the user. A
default value 1.2 can also be used.
THOLE TCUT 5.0 MAXNBTHOLE 50000
CL SODD @aij
3) Command Line Items
A) A "DRUDE" flag must be included with the GENErate command to
specify the inclusion of drudes on the molecule. The drude mass is set
in this command using the DMASS keyword followed by the mass in AMU; a
value of 0.4 is highly suggested.
generate NMA first none last none setup warn DRUDE DMASS 0.4
An additional command that may be included with the generate command
is the HYPErpolarizability term. This command invokes a higher
exponent half-wall potential on the Drude particles to minimize the
potential for polarization catastrophe. In the HYPEpolarzation command
HORD is the order of the exponent, KHYP is the force constant and RHYP
represent atom-Drude distance at which the hypepolarization term is
invoked. An example of the command follows.
Note that if DrudeHardWall, as described in the following section, use
of HYPErpolarization is not required. Experience with several systems
suggest that the DrudeHardWall feature is a more effective way than
HYPER to generate stable simulations with a 1 fs timesteps.
generate NMA first none last none setup warn DRUDE DMASS 0.4 HYPE HORD 4 KHYP 40000 RHYP 0.2
B) DrudeHardWall command
When the Drude particles move far away from the nucleus (eg. > 0.2
Ang.), simulations become unstable and polarization catastrophe may
occur. Introducing hard wall constraints on the length of Drude-nuclei
bonds is a way to avoid such instability. In each MD step, the length
of Drude bonds are checked and compared with "L_WALL". When the length
of a Drude bond is larger than "L_WALL", the relative velocities along
bond vector are flipped and scaled down according to the temperature
of the drude bond. In addition, the displacement relative to "L_WALL"
is flipped and scaled down according to new velocities along the bond
vector to make sure the Drude bond is not longer than "L_WALL".
DrudeHardWall L_WALL real
Example:
DrudeHardWall L_WALL 0.2
! set the hard wall at 0.2 Angstrom
DrudeHardWall L_WALL 0.0
! turn off the hard wall when L_WALL equals 0.0
Description of Extended PSF format for CHARMM and NAMD
The PSF has been extended for the Drude model. The default format is
for CHARMM and the XPLOR format is for NAMD. The PSF following
generation of the Drude model explicitly provides the THOLE, ALPHA,
lone pair and anisotropic polarization information for CHARMM/NAMD PSF
reader.
* Syntax | Syntax of the DRUDE command
* Description | Description of the DRUDE command
* Toppar | Effect on the topology and parameters
* Warnings | To be aware of when using the DRUDE command
* Examples | Usage examples of the DRUDE command
by Benoit Roux (roux@uchicago.edu)
and Guillaume Lamoureux (Guillaume.Lamoureux@umontreal.ca)
and Alex MacKerell Jr. (alex@outerbanks.umaryland.edu)
Classical drude particles are generated for all atoms for which
polarizabilities (ie. via the ALPHA keyword) are specified in the RTF
file. This allows for the drudes to be generated automatically when a
molecule is generated in CHARMM (only the heavy atoms carry a
polarizability in the current force field but the code is general and
allows assigning a polarizability to any atom including hydrogens). In
addition, code has been developed to allow for inclusion of atom-based
Thole scale factors, atom-based anisotropic polarizabilities and the
addition of lone pairs to selected atoms at the RTF level. These
capabilities allow for all the information for the polarizable Drude
force field to be included in the RTF and parameter files. This
implementation replaces the old "Drude Oscillator Command" and
"ANISOTROPY keyword" (see below) used to generate the Drude
polarizable model prior to July 2007.
For each atom with a specified polarizability, a "Drude oscillator" is
created by attaching to the atom an additional particle (using a
fictitious chemical bond of length zero and of force constant 'KDRUDE
= k/2'). Each Drude particle is given a mass and a charge, taken from
the mass and the charge of its atom (so that the total mass and charge
are conserved for the "atom-Drude" pair).
As a whole, each "atom-Drude" pair carries a total charge 'q',
unchanged From the partial charge the non-polarizable atom had prior
to calling the DRUDE command during system generation. The
"atom-Drude" pair forms a dipole 'q*d', where 'q' is the charge on the
Drude particle and 'd' is the displacement vector going from the atom
to its Drude particle. Any external electrostatic field 'E' creates a
net displacement 'd = q*E/k', and thus the "atom-Drude" pair behaves
as a point charge 'q' with a polarizability 'alpha = q**2/k'. The
polarizabilities (in Angst**3) are specified with the ALPHA command in
the RTF, along with the THOLE factors (see below), and converted into
charges 'q', assuming a force constant 'k = 2*KDRUDE', where KDRUDE is
specified in the parameter file (see below); a value of 500 should
always be used.
See J. Chem. Phys. 119, 3025-3039 (2003) for more details.
In the Drude treatment of polarization the force constant of the
spring between the atom-Drude pair is a diagonal rank-2 tensor with
components KDRUDE. This leads to an isotropic atomic polarizability,
'alpha = alpha_{11} = alpha_{22} = alpha_{33} = q**2/k'. The
ANISOTROPY command, which is typically included in the RTF, modifies
the components of the Drude force constant tensor allowing for
anisotropic atomic polarizabilites.
See J. Chem. Theory Comput. 2, 1587-1597 (2006) for more details.
The bonded lists are modified so that, if the "real" atoms are in a
1-2, 1-3, or 1-4 relationship, their corresponding Drude particles
will also be in a 1-2, 1-3, or 1-4 relationship, respectively.
For a single atom (charges in parenthesis):
DRUDE
A --------> A~DA
For a diatomic molecule:
A1 DRUDE A1~DA1
| --------> |
A2 A2~DA2
1-2 pairs: A1-A2, A1-DA2, DA1-A2, DA1-DA2
For a triatomic molecule:
A1 A1~DA1
\ DRUDE \
A2 --------> A2~DA2
/ /
A3 A3~DA3
1-2 pairs: A1-A2, A1-DA2, DA1-A2, DA1-DA2,
A2-A3, A2-DA3, DA2-A3, DA2-DA3
1-3 pairs: A1-A3, A1-DA3, DA1-A3, DA1-DA3
Thus, the bonded interactions are modified to allow 1-2 and 1-3
screened dipole-dipole interactions, as proposed by Thole [see
Chem.Phys. 59, 341 (1981)]. If two atoms were not "seeing" each others
in the non-polarizable force field, their dipoles (and only their
dipoles, not their net partial charges) will "see" each others in the
polarizable force field. The screening function strength is modulated
by a parameter "a" that has been generalized to depend on the atom
type associated with the interacting pair of Drude oscillators and are
then generated via a standard arithmetic combination rule (a = a_i +
a_j). These variable atomic THOLE parameters "a_i" are specified in
the RTF through the THOLE keyword.
* Description | Description of new DRUDE model RTF
Description of new DRUDE model RTF
1) Items in RTF file
A) The AUTOGENERATE command now includes a flag for the drude
particles. This should be specified at the beginning of the RTF. Note
that the use of the PATCh keyword assures that autogenerate is applied
after all patches. This should be used.
AUTO ANGLES DIHE DRUDE PATCH
B) MASS statements must be included for Drude particle types.
MASS 104 DRUD 0.00000 DD ! drude particle
The default drude type is DRUD. However, alternate drude types may be
specified. These additional drude types must be i) specified with MASS
statements, ii) the ATOMs to which they are applied must include a
"TYPE drudetype" statement and iii) the NONBond parameters must be
specified. For example, to include a special drude type on water the
following would be needed in the RTF
i) MASS 153 DOH2 0.00000 DD ! water Drude
ii) ATOM OH2 ODW 0.00000 ALPHA -0.97825258 TYPE DOH2
iii) (NONBOND section) D* 0.0 -0.0000 0.0000
! Wildcard for Drudes and dummy atom
iv) (via NBFIX) MAGD DOH2 -0.01000 3.10000
Term iv) is an NBFIX that allows for the interactions between, for
example, magnesium and the water drude to be treated with a special LJ
term. Indeed, the motivation for additional drude types are such
special terms may be applied in special cases. Note that NBFIX terms
have been used between select real atoms in the Drude force field to
improve the free energies of hydration.
See J. Chem. Theory Comput. 6. 1181-1198 (2010) for more details.
C) ALPHA and THOLE parameters are set in the RTF after the charge.
The presence of the "ALPHA" keyword flags that atom as being
polarizable and a drude particle is created and assigned to it when
calling the GENERATE command. An atom based Thole scale factor may
also be specified via the "THOLE" keyword. If the THOLE term is
missing, a default value of 1.3 is assigned to that atom. This
corresponds to a parameter a = a_i + a_j = 2.6 which is the THOLE
parameter in the old Drude command syntax.
ATOM C C 0.630 ALPHA -1.104 THOLE 1.073
Note that a negative sign on alpha leads to the charge on the Drude particle
to be negative. This convention was chosen because the Drudes are meant to
represent electronic degrees of freedom of the system.
D) Specific drude types for a polarizable atom may be specified as
follows using the "TYPE" keyword. See 1A) above, sections i to iv for
more details.
ATOM C C 0.630 ALPHA -1.104 TYPE DC
E) Lone pairs may now be set directly in the RTF using a syntax
similar to the standard lonepair command » lonepair . All options
for the generation of lonepairs (FIXed, CENTer, COLOcate etc.) as
specified in the lonepair documentation may be used although the atom
selection specification is simplified as shown in the following
example.
LONEPAIR relative LPA O C CL distance 0.3 angle 91.0 dihe 0.0
LONEPAIR relative LPB O C N distance 0.3 angle 91.0 dihe 0.0
F) Atom based anisotropic polarizabilities may be assigned to selected
atoms via the ANISOTROPY keyword in the RTF. The anisotropic
polarizability tensor is defined as ALPHA_iso (the isotropic
polarizability) times a tensor A whith trace{A}=3. The tensor A is
diagonal in the local reference frame and is completely specified by
{A11, A22, A33}. The first atom selects the atom on which the
polarizability is to be made anisotropic. The second atom along with
the first defines the "11" direction of the {A} tensor and the 3rd and
4th atom selections define the "22" direction of the tensor, with the
anisotropy defined by fractional polarizability components.
Accordingly, only the 11 and 22 direction are needed since the trace
of the {A} tensor is set to 3. An example follows
ANISOTROPY O C CA N A11 0.6968 A22 1.2194
In the code, the variables A11 and A22 are read into the variables
A11ANIS(IPT) and A22ANIS(IPT). We define the components of the tensor
relative to the isotropic drude, so ALPHA11+ALPHA22+ALPHA33 =
3*ALPHA_iso. The factors A11 and A22 (and A33) represent the factors
modifying the components of the isotropic polarizability tensor
q**2/KDRUDE. For example, setting A11=A22=1 would keep the Drude
completely isotropic. By virtue of the trace we know that
A11+A22+A33=3, so we only need to set two of these terms. Thus A33 =
THREE-A11ANIS(IPT)-A22ANIS(IPT). This sets the relative magnitude of
the polarizability tensor in all 3 directions.
Once the polarizability along a given direction is known, we need to
get the harmonic force constant in this direction. As the
polarizability in a given direction is defined as
QDRUDE**2/K_direction, the 3 harmonic force constants are thus
K11(NANISO) = ONE/A11ANIS(IPT)
K22(NANISO) = ONE/A22ANIS(IPT)
K33(NANISO) = ONE/A33
Now, in principle this would mean that you have to calculate the
forces from 3 different springs (K11, K22, K33) in 3 different
directions. But since we know the overall trace which is related to
the isotropic polarizability, we are going to re-write this so that
there are only 2 springs with force constants KPAR in the parallel and
KPERP in the perpendicular direction, and add an isotropic KDRUDE
spring on top of that.
2) Items in the Parameter file
A) The KDRUDE force constant for all the atoms is set by a wildcard in
the bond parameters. The term must be included for all drude types
included in the model. The wildcard terms can be overwritten by
putting chemical type specific bond parameters following the wild card
term. Note that 500.00 is the default force constant and should not be
changed.
DRUD X 500.000 0.0000
DRUD O 487.740 0.0000
B) NONBOND parameters must be included for all drude types. In
addition, NBFIX terms for the drudes may be included as needed (see 1A
above).
C) A pair-specific Thole shielding can be introduced in the parameter
file for any nonbond pair (following a syntax similar to NBFIX) using
the keyword THOLe. This feature was introduced to permit accurate
models of divalent cations with polarizable atoms. The TCUT assigns a
radius within which all Thole interaction are calculated between the
specificed atoms. The 'aij' typically is assigned by the user. A
default value 1.2 can also be used.
THOLE TCUT 5.0 MAXNBTHOLE 50000
CL SODD @aij
3) Command Line Items
A) A "DRUDE" flag must be included with the GENErate command to
specify the inclusion of drudes on the molecule. The drude mass is set
in this command using the DMASS keyword followed by the mass in AMU; a
value of 0.4 is highly suggested.
generate NMA first none last none setup warn DRUDE DMASS 0.4
An additional command that may be included with the generate command
is the HYPErpolarizability term. This command invokes a higher
exponent half-wall potential on the Drude particles to minimize the
potential for polarization catastrophe. In the HYPEpolarzation command
HORD is the order of the exponent, KHYP is the force constant and RHYP
represent atom-Drude distance at which the hypepolarization term is
invoked. An example of the command follows.
Note that if DrudeHardWall, as described in the following section, use
of HYPErpolarization is not required. Experience with several systems
suggest that the DrudeHardWall feature is a more effective way than
HYPER to generate stable simulations with a 1 fs timesteps.
generate NMA first none last none setup warn DRUDE DMASS 0.4 HYPE HORD 4 KHYP 40000 RHYP 0.2
B) DrudeHardWall command
When the Drude particles move far away from the nucleus (eg. > 0.2
Ang.), simulations become unstable and polarization catastrophe may
occur. Introducing hard wall constraints on the length of Drude-nuclei
bonds is a way to avoid such instability. In each MD step, the length
of Drude bonds are checked and compared with "L_WALL". When the length
of a Drude bond is larger than "L_WALL", the relative velocities along
bond vector are flipped and scaled down according to the temperature
of the drude bond. In addition, the displacement relative to "L_WALL"
is flipped and scaled down according to new velocities along the bond
vector to make sure the Drude bond is not longer than "L_WALL".
DrudeHardWall L_WALL real
Example:
DrudeHardWall L_WALL 0.2
! set the hard wall at 0.2 Angstrom
DrudeHardWall L_WALL 0.0
! turn off the hard wall when L_WALL equals 0.0
Description of Extended PSF format for CHARMM and NAMD
The PSF has been extended for the Drude model. The default format is
for CHARMM and the XPLOR format is for NAMD. The PSF following
generation of the Drude model explicitly provides the THOLE, ALPHA,
lone pair and anisotropic polarization information for CHARMM/NAMD PSF
reader.
* Syntax | Syntax of the DRUDE command
* Description | Description of the DRUDE command
* Toppar | Effect on the topology and parameters
* Warnings | To be aware of when using the DRUDE command
* Examples | Usage examples of the DRUDE command
Top
Syntax of the DRUDE command
DRUDe RESEt
GENErate { generate-spec } DRUDe [DMASS real] [KDRUde real] [HYPE] { hype-spec }
generate-spec::= [» struct ]
Syntax of the DRUDE command
DRUDe RESEt
GENErate { generate-spec } DRUDe [DMASS real] [KDRUde real] [HYPE] { hype-spec }
generate-spec::= [» struct ]
Top
Description of the DRUDE command
-------------------------------------------------------------------
Keyword Default Purpose
-------------------------------------------------------------------
RESET Deactivates the Drude particles. After a reset,
the user should delete all Drude particles.
DMASS 0.0 Mass of the Drude oscillators (in amu). The
default zero value causes the masses to be read
from the topology file. Any nonzero value will
override the topology file. A value of 0.4 is
suggested.
KDRUDE 0.0 Force constant of the atom-Drude bonds
(in kcal/mol/Angst**2). The default zero value
causes the bond force constants to be read from
the parameter file. Any nonzero value will
override the parameter file. A value of 500.
is suggested.
VTHOLE Uses variable thole parameters for dipole-dipole
interactions. This term is set by default and is only
used in the fitcharge routine.
ALPHA 0.0 Atomic polarizability specificed in the parameter file.
Zero indicates no polarizability on that atom.
THOLE 2.6 The screening factor for dipole-dipole interactions
between atoms excluded from the non-bonded
interactions. To have no dipole-dipole
interactions between these bonded atoms, use
THOLE = 0. THOLE is specified in the parameter file.
SHAPE 1 Specifies the shape of the dipole-dipole
screening function. Only 1 option is
currently available.
HYPE Off Specify the inclusion of hyperpolarization term on
the Drude oscillators to minimize potential of
polarization catastophe
hype-spec::= [HORD int] [KHYP int] [RHYP real]
(XXXXX defaults have not been set for the following, may be a good idea XXXXXX)
HORD xx Order of the exponent the hyperpolarization. Typically 4.
KHYP xx Force constant. Typically 40,000.
RHYP xx Atom-Drude distance of the flat-well potential. The
hyperpolarization is not activited until the atom-Drude
distance is greater than RHYP. Typically 0.2 Angst.
-------------------------------------------------------------------
1) KDRUDE
KDRUDE is the force constant (in kcal/mol/Angst**2) for the bond
between each atom and its Drude particle the user wishes to use. It
overrides any bond constant found in the parameter file. The default
value for KDRUDE is 500 kcal/mol/Angst**2. This should NOT be
changed!
The atomic polarizabilities (in Angst**3) may be read from the WMAIN
array:
alpha = abs(WMAIN)
The charge on every Drude particle is computed using the following
formula:
q = sqrt( 2*KDRUDE * alpha / CCELEC ) * sign(WMAIN)
The charges are given the signs of the WMAIN values. As long as KDRUDE
is large enough, the Drude particles will stay very close to their
atoms, and the sign of 'q' is irrelevant. However, as the Drude
represents the electronic degrees of freedom, the sign of alpha is set
to a negative value, thereby yielding a negative charge on the Drude
particles.
-------------------------------------------------------------------
2) THOLE, SHAPE
For 1-2 and 1-3 atoms only. All other interactions are regular
Coulombic. See THOLE keyword in the parameter file to perform Thole
scaling between nonbonded atoms.
The THOLE parameter is a dimensionless factor that specifies the
extent of the smearing of the charge 'q' on the Drude oscillators and
of a contribution '-q' on the "real" atom. The default value of THOLE
is 2.6, that is, the 1-2 and 1-3 dipole-dipole interactions are turned
on by default. To turn the interactions off, set THOLE to zero.
The default constant THOLE of 2.6 can be replaced by variable thole by
specifying an atom specific Thole value in the topology file. The
THOLE parameter between oscillators I and J is given by THOLE = THOLEI
+ THOLEJ. Values of THOLEI may also be given in the WCOMP array for
all Drude oscillator containing atoms prior to the DRUDE command.
Because the dipole are explicitly made of two charges, the screened
dipole-dipole interaction between two polarizable atoms (that is, two
"atom-Drude" pairs) is actually the sum of the following four screened
charge-charge interactions:
('q1' on Drude 1) - ('q2' on Drude 2)
('q1' on Drude 1) - ('-q2' on atom 2)
('-q1' on atom 1) - ('q2' on Drude 2)
('-q1' on atom 1) - ('-q2' on atom 2)
The screened charge-charge interaction has the form:
U(r12) = CCELEC * q1*q2 * S(u12)/r12
where 'u12' is the normalized distance:
u12 = r12 * THOLE / (alpha1*alpha2)**(1/6)
'S' is a screening function defined by the SHAPE parameter:
SHAPE Screening function Charge distributions
1 S(u) = 1 - (1+u/2)*exp(-u) Slater-Delta
2 S(u) = erf(u) Gaussian
The default value of SHAPE is 1, which is also the only shape
currently implemented. SHAPE=2 is reserved for Gaussian-Gaussian
distributions.
Two "atom-Drude" pairs have dipole-dipole interactions if the
following conditions are met:
1. The THOLE parameter is nonzero.
2. In the non-polarizable force field, the two atoms were
in the nonbonded exclusion list.
To see if all the desired atoms have dipole-dipole interactions, use
PRNLEV > 7. Each call to the energy will print the atom numbers,
polarizabilities and Drude particles's charges of each interacting
pair.
The energy from the dipole-dipole interactions is added to the ELEC
energy term, and "SKIP ELEC" will skip the Thole interactions as well.
Description of the DRUDE command
-------------------------------------------------------------------
Keyword Default Purpose
-------------------------------------------------------------------
RESET Deactivates the Drude particles. After a reset,
the user should delete all Drude particles.
DMASS 0.0 Mass of the Drude oscillators (in amu). The
default zero value causes the masses to be read
from the topology file. Any nonzero value will
override the topology file. A value of 0.4 is
suggested.
KDRUDE 0.0 Force constant of the atom-Drude bonds
(in kcal/mol/Angst**2). The default zero value
causes the bond force constants to be read from
the parameter file. Any nonzero value will
override the parameter file. A value of 500.
is suggested.
VTHOLE Uses variable thole parameters for dipole-dipole
interactions. This term is set by default and is only
used in the fitcharge routine.
ALPHA 0.0 Atomic polarizability specificed in the parameter file.
Zero indicates no polarizability on that atom.
THOLE 2.6 The screening factor for dipole-dipole interactions
between atoms excluded from the non-bonded
interactions. To have no dipole-dipole
interactions between these bonded atoms, use
THOLE = 0. THOLE is specified in the parameter file.
SHAPE 1 Specifies the shape of the dipole-dipole
screening function. Only 1 option is
currently available.
HYPE Off Specify the inclusion of hyperpolarization term on
the Drude oscillators to minimize potential of
polarization catastophe
hype-spec::= [HORD int] [KHYP int] [RHYP real]
(XXXXX defaults have not been set for the following, may be a good idea XXXXXX)
HORD xx Order of the exponent the hyperpolarization. Typically 4.
KHYP xx Force constant. Typically 40,000.
RHYP xx Atom-Drude distance of the flat-well potential. The
hyperpolarization is not activited until the atom-Drude
distance is greater than RHYP. Typically 0.2 Angst.
-------------------------------------------------------------------
1) KDRUDE
KDRUDE is the force constant (in kcal/mol/Angst**2) for the bond
between each atom and its Drude particle the user wishes to use. It
overrides any bond constant found in the parameter file. The default
value for KDRUDE is 500 kcal/mol/Angst**2. This should NOT be
changed!
The atomic polarizabilities (in Angst**3) may be read from the WMAIN
array:
alpha = abs(WMAIN)
The charge on every Drude particle is computed using the following
formula:
q = sqrt( 2*KDRUDE * alpha / CCELEC ) * sign(WMAIN)
The charges are given the signs of the WMAIN values. As long as KDRUDE
is large enough, the Drude particles will stay very close to their
atoms, and the sign of 'q' is irrelevant. However, as the Drude
represents the electronic degrees of freedom, the sign of alpha is set
to a negative value, thereby yielding a negative charge on the Drude
particles.
-------------------------------------------------------------------
2) THOLE, SHAPE
For 1-2 and 1-3 atoms only. All other interactions are regular
Coulombic. See THOLE keyword in the parameter file to perform Thole
scaling between nonbonded atoms.
The THOLE parameter is a dimensionless factor that specifies the
extent of the smearing of the charge 'q' on the Drude oscillators and
of a contribution '-q' on the "real" atom. The default value of THOLE
is 2.6, that is, the 1-2 and 1-3 dipole-dipole interactions are turned
on by default. To turn the interactions off, set THOLE to zero.
The default constant THOLE of 2.6 can be replaced by variable thole by
specifying an atom specific Thole value in the topology file. The
THOLE parameter between oscillators I and J is given by THOLE = THOLEI
+ THOLEJ. Values of THOLEI may also be given in the WCOMP array for
all Drude oscillator containing atoms prior to the DRUDE command.
Because the dipole are explicitly made of two charges, the screened
dipole-dipole interaction between two polarizable atoms (that is, two
"atom-Drude" pairs) is actually the sum of the following four screened
charge-charge interactions:
('q1' on Drude 1) - ('q2' on Drude 2)
('q1' on Drude 1) - ('-q2' on atom 2)
('-q1' on atom 1) - ('q2' on Drude 2)
('-q1' on atom 1) - ('-q2' on atom 2)
The screened charge-charge interaction has the form:
U(r12) = CCELEC * q1*q2 * S(u12)/r12
where 'u12' is the normalized distance:
u12 = r12 * THOLE / (alpha1*alpha2)**(1/6)
'S' is a screening function defined by the SHAPE parameter:
SHAPE Screening function Charge distributions
1 S(u) = 1 - (1+u/2)*exp(-u) Slater-Delta
2 S(u) = erf(u) Gaussian
The default value of SHAPE is 1, which is also the only shape
currently implemented. SHAPE=2 is reserved for Gaussian-Gaussian
distributions.
Two "atom-Drude" pairs have dipole-dipole interactions if the
following conditions are met:
1. The THOLE parameter is nonzero.
2. In the non-polarizable force field, the two atoms were
in the nonbonded exclusion list.
To see if all the desired atoms have dipole-dipole interactions, use
PRNLEV > 7. Each call to the energy will print the atom numbers,
polarizabilities and Drude particles's charges of each interacting
pair.
The energy from the dipole-dipole interactions is added to the ELEC
energy term, and "SKIP ELEC" will skip the Thole interactions as well.
Top
Effect on the topology and parameters
-------------------------------------------------------------------
1) New atoms
The Drude particles are inserted immediately after their corresponding
atoms in the PSF and coordinate file. For an atom type 'CA1', the
DRUDE command will assign the atom type 'DCA1' for the Drude particle.
Since no regular atoms have names starting with a 'D', the Drude
oscillators can be selected with "SELECT TYPE D* END".
-------------------------------------------------------------------
2) Masses
The masses for the selected atoms are modified so that the total mass
of the atom-Drude pair corresponds to the atomic mass. Try "SCALAR
MASS SHOW" before and after calling the DRUDE command.
-------------------------------------------------------------------
3) Charges
The charges for the selected atoms are modified so that the total
charge of the atom-Drude pair corresponds to the atomic partial
charge. Try "SCALAR CHARGE SHOW" before and after calling the DRUDE
command.
-------------------------------------------------------------------
4) Bonded interactions
In addition to the atom-Drude bonds, zero force bonds are created
between the atom-Drude pairs to maintain the same 1-2, 1-3, and 1-4
relationships that existed prior to the DRUDE call. For example, for
two bonded atoms A1 and A2, with Drude particles DA1 and DA2,
DA1 DA2
| |
A1 ----- A2
zero force bonds are created between DA1 and DA2, between DA1 and A2,
and between A1 and DA2, so that any particle of the A1-DA1 pair is 1-2
bonded to any particle of the A2-DA2 pair. Since the force constants
of these fictitious bonds are zero, the computational overhead is
minimal.
-------------------------------------------------------------------
5) Non-bonded interactions
Whether the Drude particles have Lennard-Jones parameters or not, the
Lennard-Jones parameters of the selected atoms are kept unchanged.
Since the polarizable force field is built from the same "ingredients"
as the non-polarizable force field, all the NBONDS options can be used
as before (notably the PMEWALD method).
Effect on the topology and parameters
-------------------------------------------------------------------
1) New atoms
The Drude particles are inserted immediately after their corresponding
atoms in the PSF and coordinate file. For an atom type 'CA1', the
DRUDE command will assign the atom type 'DCA1' for the Drude particle.
Since no regular atoms have names starting with a 'D', the Drude
oscillators can be selected with "SELECT TYPE D* END".
-------------------------------------------------------------------
2) Masses
The masses for the selected atoms are modified so that the total mass
of the atom-Drude pair corresponds to the atomic mass. Try "SCALAR
MASS SHOW" before and after calling the DRUDE command.
-------------------------------------------------------------------
3) Charges
The charges for the selected atoms are modified so that the total
charge of the atom-Drude pair corresponds to the atomic partial
charge. Try "SCALAR CHARGE SHOW" before and after calling the DRUDE
command.
-------------------------------------------------------------------
4) Bonded interactions
In addition to the atom-Drude bonds, zero force bonds are created
between the atom-Drude pairs to maintain the same 1-2, 1-3, and 1-4
relationships that existed prior to the DRUDE call. For example, for
two bonded atoms A1 and A2, with Drude particles DA1 and DA2,
DA1 DA2
| |
A1 ----- A2
zero force bonds are created between DA1 and DA2, between DA1 and A2,
and between A1 and DA2, so that any particle of the A1-DA1 pair is 1-2
bonded to any particle of the A2-DA2 pair. Since the force constants
of these fictitious bonds are zero, the computational overhead is
minimal.
-------------------------------------------------------------------
5) Non-bonded interactions
Whether the Drude particles have Lennard-Jones parameters or not, the
Lennard-Jones parameters of the selected atoms are kept unchanged.
Since the polarizable force field is built from the same "ingredients"
as the non-polarizable force field, all the NBONDS options can be used
as before (notably the PMEWALD method).
Top
Issues to be aware of when using the DRUDE command
-------------------------------------------------------------------
1) Once all molecules in the system have been generated and any patching
performed, invoke the following commands to create the drude and lone pair
coordinates and set up the needed arrays.
COOR SDRUDE
COOR SHAKE
-------------------------------------------------------------------
2) The preferred call to SHAKE is:
COOR COPY COMP
SHAKE BONH PARAM TOLERANCE 10E-9 -
NOFAST -
SELECT .NOT. TYPE D* END -
SELECT .NOT. TYPE D* END
-------------------------------------------------------------------
3) Always delete all the Drude particles after a RESET
Treat this sequence of commands a single command:
DRUDE RESET
DELETE ATOMS SELECT TYPE D* END
The "DRUDE RESET" command puts back the mass and the charge of the
Drude particles on the heavy atoms, and erases the distinction between
a Drude particle and a regular atom.
Note to fully clear all the relevant arrays to allow for a new molecule
to be generated all the following commands must be performed.
SHAKE OFF
DRUDE RESET
LONE CLEAR
DELETE ATOM SELE ALL END
-------------------------------------------------------------------
4) DMASS and KDRUDE are overriding the toppar files
Any nonzero value for DMASS and KDRUDE specified by the user in a
GENErate statement overrides the values from the topology and
parameter files.
-------------------------------------------------------------------
5) Syntax and Description of the DrudeHardWall command
DrudeHardWall L_WALL real
Example:
DrudeHardWall L_WALL 0.2
! set the hard wall at 0.2 Angstrom
DrudeHardWall L_WALL 0.0
! turn off the hard wall when L_WALL equals 0.0
When the Drudes move too far away from the nucleus, the simulations
become unstable leading to polarization catastrophe. To avoid this
instability and "hard wall constraint" on the length of Drude-nuclei
bonds is introduced. In each MD step, the length of Drude bonds are
checked and compared with "L_WALL". When the length of Drude bond is
larger than "L_WALL", the relative velocities along the bond vector
are reversed and scaled down according to the temperature of the Drude
bond. In addition, the displacement relative to "L_WALL" is reversed
and scaled down according to new the velocities along the bond vector
to make sure the Drude bond is not longer than "L_WALL".
Issues to be aware of when using the DRUDE command
-------------------------------------------------------------------
1) Once all molecules in the system have been generated and any patching
performed, invoke the following commands to create the drude and lone pair
coordinates and set up the needed arrays.
COOR SDRUDE
COOR SHAKE
-------------------------------------------------------------------
2) The preferred call to SHAKE is:
COOR COPY COMP
SHAKE BONH PARAM TOLERANCE 10E-9 -
NOFAST -
SELECT .NOT. TYPE D* END -
SELECT .NOT. TYPE D* END
-------------------------------------------------------------------
3) Always delete all the Drude particles after a RESET
Treat this sequence of commands a single command:
DRUDE RESET
DELETE ATOMS SELECT TYPE D* END
The "DRUDE RESET" command puts back the mass and the charge of the
Drude particles on the heavy atoms, and erases the distinction between
a Drude particle and a regular atom.
Note to fully clear all the relevant arrays to allow for a new molecule
to be generated all the following commands must be performed.
SHAKE OFF
DRUDE RESET
LONE CLEAR
DELETE ATOM SELE ALL END
-------------------------------------------------------------------
4) DMASS and KDRUDE are overriding the toppar files
Any nonzero value for DMASS and KDRUDE specified by the user in a
GENErate statement overrides the values from the topology and
parameter files.
-------------------------------------------------------------------
5) Syntax and Description of the DrudeHardWall command
DrudeHardWall L_WALL real
Example:
DrudeHardWall L_WALL 0.2
! set the hard wall at 0.2 Angstrom
DrudeHardWall L_WALL 0.0
! turn off the hard wall when L_WALL equals 0.0
When the Drudes move too far away from the nucleus, the simulations
become unstable leading to polarization catastrophe. To avoid this
instability and "hard wall constraint" on the length of Drude-nuclei
bonds is introduced. In each MD step, the length of Drude bonds are
checked and compared with "L_WALL". When the length of Drude bond is
larger than "L_WALL", the relative velocities along the bond vector
are reversed and scaled down according to the temperature of the Drude
bond. In addition, the displacement relative to "L_WALL" is reversed
and scaled down according to new the velocities along the bond vector
to make sure the Drude bond is not longer than "L_WALL".
Top
Usage examples of the DRUDE command
-------------------------------------------------------------------
1) Polarizable benzene
(see test/c30test/drude_benzene.inp)
After reading the standard, non-polarizable topology and parameter
files, a standard benzene molecule is generated:
READ SEQUENCE BENZ 1
GENERATE BENZ SETUP FIRST NONE LAST NONE DRUDE DMASS 0.4
Obtain the benzene coordinates using IC build or reading in the
coordinates. READ COOR CARD....
Then generate drude and lone pair coordinates.
COOR SDRUDE
COOR SHAKE
The structure is minimized:
MINI SD STEP 0.001 NSTEP 1000 NPRINT 100
Since benzene is a nonpolar molecule, the Drude particles are not
significantly moved from their heavy atoms. To find the induced
atomic dipoles for a given structure, one should use "CONS FIX SELECT
.NOT. TYPE D* END" before calling MINI.
The molecular polarizability is obtained using the VIBRAN command with
the DIPOLES keyword:
VIBRAN
DIAGONALIZE
PRINT NORMAL VECTORS DIPOLES SELECT ALL END
END
The total polarizability is an anisotropic tensor similar to the
experimental results for benzene [J.Chem.Phys. 95, 5873 (1991)]. The
strong anisotropy comes from the 1-2 and 1-3 dipole-dipole
interactions. Deactivating these interactions by using THOLE=0 leads
to the polarizability tensor being almost isotropic.
-------------------------------------------------------------------
2) SWM4-NDP water
(see test/c30test/swm4.inp)
See J. Chem. Phys. 119, 5185-5197 (2003) and Chem. Phys. Lett. 418,
245-249 (2005) for a complete description of the model. After reading
the topology and parameter files, the model is built as following:
READ SEQUENCE SWM4 ...
GENERATE WAT SETUP NOANGLE NODIHEDRAL DRUDE DMASS 0.4
READ COOR CARD ...
Then generate drude and lone pair coordinates.
COOR SDRUDE
COOR SHAKE
SHAKE the waters
COOR COPY COMP
SHAKE BONH PARAM TOLERANCE 1.0E-9 -
NOFAST -
SELECT ( SEGID WAT .AND. .NOT. TYPE D* ) END -
SELECT ( SEGID WAT .AND. .NOT. TYPE D* ) END
-------------------------------------------------------------------
3) Protein, with disulfide patches, SWM4-NDP water and ions
Note that for complex systems it is suggested that the system
initially be equilibrated using the additive force field, with the
final coordinates used to initiate the Drude simulations. Tools to
perform this may be found on the CHARMM Gui or in the MacKerell web
site.
! Read protein
open read card unit 10 name protein.crd
read sequence coor unit 10
generate PROA setup warn first NTER last CTER DRUDE DMASS 0.4
! NOTE: the DRUDE statement after 'patch' is not necessary Disulfide
! bonds and if PATCH is included in the AUTOgenerate command it is not
! needed following patches.
patch disu2 PROA 26 PROA 84 setup warn DRUDE DMASS 0.4
autogenerate angles dihedrals
patch disu2 PROA 40 PROA 95 setup warn DRUDE DMASS 0.4
autogenerate angles dihedrals
patch disu2 PROA 58 PROA 110 setup warn DRUDE DMASS 0.4
autogenerate angles dihedrals
patch disu2 PROA 65 PROA 72 setup warn DRUDE DMASS 0.4
autogenerate angles dihedrals
! Read xtal waters
open read card unit 10 name xtal_water.crd
read sequence coor unit 10
generate WATA setup warn noangle nodihedral DRUDE DMASS 0.4
read coor card unit 10 resid
! Read solvent waters
open read card unit 10 name solv_water.crd
read sequence coor unit 10
generate SOLV setup warn noangle nodihedral DRUDE DMASS 0.4
read coor card unit 10 resid
! Read IONS
open read card unit 10 name ions.crd
read sequence coor unit 10
generate CL setup warn noangle nodihedral DRUDE DMASS 0.4
read coor card unit 10 resid
COOR SDRUDE
COOR SHAKE
! cell parameters
SET BOXTYPE = RECT
SET XTLTYPE = CUBIC
SET A = 62.75993
SET B = 62.75993
SET C = 62.75993
SET ALPHA = 90.0
SET BETA = 90.0
SET GAMMA = 90.0
SET XCEN = 0
SET YCEN = 0
SET ZCEN = 0
! Setup PBC (Periodic Boundary Condition)
CRYSTAL DEFINE @XTLtype @A @B @C @alpha @beta @gamma
crystal build noper 0 cutoff @3
!Image centering by residue
IMAGE BYRESID XCEN @xcen YCEN @ycen ZCEN @zcen sele resname SWM4 end
IMAGE BYRESID XCEN @xcen YCEN @ycen ZCEN @zcen sele ( segid pot .or. segid cl ) end
! Nonbonded Options
set 3 16.0 ! cutim
set 4 16.0 ! cutnb
set 5 10.0 ! ctonnb
set 6 12.0 ! ctofnb
set 7 switch
set 8 atom
set 9 vatom
set 10 vswitch !Note use of vswitch
set kap 0.34 ! epsilon for the PME
set ord 6 ! spline order for PME
! fftx dimensions should correspond to 1 A or less; must
! be power of 2, 3 or 5
set fftx 64
set ffty 64
set fftz 64
!invoke drudehardwall to avoid polarization catastrophe.
DrudeHardWall L_WALL 0.2
update inbfrq -1 imgfrq -1 ihbfrq 0 -
ewald pmewald kappa @kap fftx @fftx ffty @ffty fftz @fftz order @ord -
@7 @8 @9 @10 cutimg @3 cutnb @4 ctofnb @6 ctonnb @5
energy
!minimize only the Drudes
cons harm force 100000.0 sele .not. type D* end
mini SD nstep 50
cons harm force 0.0 sele .not. type D* end
coor copy comp
SHAKE bonh param tolerance 1.0e-9 -
nofast -
select ( .not. (type D*)) end -
select ( .not. (type D*)) end noreset
!minimize full system
mini SD nstep 50
! with proper equilibration, Hyperpolarization and/or Hardwall
! 1 fs timestep should be possible
set timestep 0.0005 !0.001
set nstep 200000
set seed 21320
set ii 0
set jj 1
open read unit 11 card name @tmp/hsd12.psf@psfid.@ii.res
open write unit 12 card name @tmp/hsd12.psf@psfid.@jj.res
open write unit 13 file name @tmp/hsd12.psf@psfid.@jj.dcd
! Set up temperature control for NVT simulation
!scf THER 2 TAU 0.005 TREF 0.00 SELECT TYPE D* END -
TPCONTROL NTHER 4 CMDAMP 10.0 NSTEP 20 -
THER 1 TAU 0.1 TREF @temp SELECT .NOT. (resname swm4 .or. TYPE D* ) END -
THER 2 TAU 0.1 TREF @temp SELECT resname swm4 .and. .NOT. TYPE D* END -
THER 3 TAU 0.005 TREF 1.00 SELECT TYPE D* .and. .not. resname swm4 END -
THER 4 TAU 0.005 TREF 1.00 SELECT resname swm4 .and. TYPE D* END !-
! BARO BTAU 0.1 PREF 1.00 DSCY !Activate BARO for NPT simulation
! Setup and run the Thermodynamic Integration loop using Langevine
! Dynamics.
DYNAMICS vv2 start nstep @nstep timestp @timestep -
ntrfrq @nstep iprfrq -1 -
nprint 1000 iseed @seed -
iasvel 1 firstt @temp finalt @temp -
inbfrq -1 imgfrq -1 ihbfrq 0 ilbfrq 0 -
iunread 11 -
iunwrite 12 -
iuncrd 13 nsavcrd 2000
Details of the command to perform molecular dynamics for the polarizable Drude mdoel is explained in
» vv2
Usage examples of the DRUDE command
-------------------------------------------------------------------
1) Polarizable benzene
(see test/c30test/drude_benzene.inp)
After reading the standard, non-polarizable topology and parameter
files, a standard benzene molecule is generated:
READ SEQUENCE BENZ 1
GENERATE BENZ SETUP FIRST NONE LAST NONE DRUDE DMASS 0.4
Obtain the benzene coordinates using IC build or reading in the
coordinates. READ COOR CARD....
Then generate drude and lone pair coordinates.
COOR SDRUDE
COOR SHAKE
The structure is minimized:
MINI SD STEP 0.001 NSTEP 1000 NPRINT 100
Since benzene is a nonpolar molecule, the Drude particles are not
significantly moved from their heavy atoms. To find the induced
atomic dipoles for a given structure, one should use "CONS FIX SELECT
.NOT. TYPE D* END" before calling MINI.
The molecular polarizability is obtained using the VIBRAN command with
the DIPOLES keyword:
VIBRAN
DIAGONALIZE
PRINT NORMAL VECTORS DIPOLES SELECT ALL END
END
The total polarizability is an anisotropic tensor similar to the
experimental results for benzene [J.Chem.Phys. 95, 5873 (1991)]. The
strong anisotropy comes from the 1-2 and 1-3 dipole-dipole
interactions. Deactivating these interactions by using THOLE=0 leads
to the polarizability tensor being almost isotropic.
-------------------------------------------------------------------
2) SWM4-NDP water
(see test/c30test/swm4.inp)
See J. Chem. Phys. 119, 5185-5197 (2003) and Chem. Phys. Lett. 418,
245-249 (2005) for a complete description of the model. After reading
the topology and parameter files, the model is built as following:
READ SEQUENCE SWM4 ...
GENERATE WAT SETUP NOANGLE NODIHEDRAL DRUDE DMASS 0.4
READ COOR CARD ...
Then generate drude and lone pair coordinates.
COOR SDRUDE
COOR SHAKE
SHAKE the waters
COOR COPY COMP
SHAKE BONH PARAM TOLERANCE 1.0E-9 -
NOFAST -
SELECT ( SEGID WAT .AND. .NOT. TYPE D* ) END -
SELECT ( SEGID WAT .AND. .NOT. TYPE D* ) END
-------------------------------------------------------------------
3) Protein, with disulfide patches, SWM4-NDP water and ions
Note that for complex systems it is suggested that the system
initially be equilibrated using the additive force field, with the
final coordinates used to initiate the Drude simulations. Tools to
perform this may be found on the CHARMM Gui or in the MacKerell web
site.
! Read protein
open read card unit 10 name protein.crd
read sequence coor unit 10
generate PROA setup warn first NTER last CTER DRUDE DMASS 0.4
! NOTE: the DRUDE statement after 'patch' is not necessary Disulfide
! bonds and if PATCH is included in the AUTOgenerate command it is not
! needed following patches.
patch disu2 PROA 26 PROA 84 setup warn DRUDE DMASS 0.4
autogenerate angles dihedrals
patch disu2 PROA 40 PROA 95 setup warn DRUDE DMASS 0.4
autogenerate angles dihedrals
patch disu2 PROA 58 PROA 110 setup warn DRUDE DMASS 0.4
autogenerate angles dihedrals
patch disu2 PROA 65 PROA 72 setup warn DRUDE DMASS 0.4
autogenerate angles dihedrals
! Read xtal waters
open read card unit 10 name xtal_water.crd
read sequence coor unit 10
generate WATA setup warn noangle nodihedral DRUDE DMASS 0.4
read coor card unit 10 resid
! Read solvent waters
open read card unit 10 name solv_water.crd
read sequence coor unit 10
generate SOLV setup warn noangle nodihedral DRUDE DMASS 0.4
read coor card unit 10 resid
! Read IONS
open read card unit 10 name ions.crd
read sequence coor unit 10
generate CL setup warn noangle nodihedral DRUDE DMASS 0.4
read coor card unit 10 resid
COOR SDRUDE
COOR SHAKE
! cell parameters
SET BOXTYPE = RECT
SET XTLTYPE = CUBIC
SET A = 62.75993
SET B = 62.75993
SET C = 62.75993
SET ALPHA = 90.0
SET BETA = 90.0
SET GAMMA = 90.0
SET XCEN = 0
SET YCEN = 0
SET ZCEN = 0
! Setup PBC (Periodic Boundary Condition)
CRYSTAL DEFINE @XTLtype @A @B @C @alpha @beta @gamma
crystal build noper 0 cutoff @3
!Image centering by residue
IMAGE BYRESID XCEN @xcen YCEN @ycen ZCEN @zcen sele resname SWM4 end
IMAGE BYRESID XCEN @xcen YCEN @ycen ZCEN @zcen sele ( segid pot .or. segid cl ) end
! Nonbonded Options
set 3 16.0 ! cutim
set 4 16.0 ! cutnb
set 5 10.0 ! ctonnb
set 6 12.0 ! ctofnb
set 7 switch
set 8 atom
set 9 vatom
set 10 vswitch !Note use of vswitch
set kap 0.34 ! epsilon for the PME
set ord 6 ! spline order for PME
! fftx dimensions should correspond to 1 A or less; must
! be power of 2, 3 or 5
set fftx 64
set ffty 64
set fftz 64
!invoke drudehardwall to avoid polarization catastrophe.
DrudeHardWall L_WALL 0.2
update inbfrq -1 imgfrq -1 ihbfrq 0 -
ewald pmewald kappa @kap fftx @fftx ffty @ffty fftz @fftz order @ord -
@7 @8 @9 @10 cutimg @3 cutnb @4 ctofnb @6 ctonnb @5
energy
!minimize only the Drudes
cons harm force 100000.0 sele .not. type D* end
mini SD nstep 50
cons harm force 0.0 sele .not. type D* end
coor copy comp
SHAKE bonh param tolerance 1.0e-9 -
nofast -
select ( .not. (type D*)) end -
select ( .not. (type D*)) end noreset
!minimize full system
mini SD nstep 50
! with proper equilibration, Hyperpolarization and/or Hardwall
! 1 fs timestep should be possible
set timestep 0.0005 !0.001
set nstep 200000
set seed 21320
set ii 0
set jj 1
open read unit 11 card name @tmp/hsd12.psf@psfid.@ii.res
open write unit 12 card name @tmp/hsd12.psf@psfid.@jj.res
open write unit 13 file name @tmp/hsd12.psf@psfid.@jj.dcd
! Set up temperature control for NVT simulation
!scf THER 2 TAU 0.005 TREF 0.00 SELECT TYPE D* END -
TPCONTROL NTHER 4 CMDAMP 10.0 NSTEP 20 -
THER 1 TAU 0.1 TREF @temp SELECT .NOT. (resname swm4 .or. TYPE D* ) END -
THER 2 TAU 0.1 TREF @temp SELECT resname swm4 .and. .NOT. TYPE D* END -
THER 3 TAU 0.005 TREF 1.00 SELECT TYPE D* .and. .not. resname swm4 END -
THER 4 TAU 0.005 TREF 1.00 SELECT resname swm4 .and. TYPE D* END !-
! BARO BTAU 0.1 PREF 1.00 DSCY !Activate BARO for NPT simulation
! Setup and run the Thermodynamic Integration loop using Langevine
! Dynamics.
DYNAMICS vv2 start nstep @nstep timestp @timestep -
ntrfrq @nstep iprfrq -1 -
nprint 1000 iseed @seed -
iasvel 1 firstt @temp finalt @temp -
inbfrq -1 imgfrq -1 ihbfrq 0 ilbfrq 0 -
iunread 11 -
iunwrite 12 -
iuncrd 13 nsavcrd 2000
Details of the command to perform molecular dynamics for the polarizable Drude mdoel is explained in
» vv2