Drude Oscillator Command
by Benoit Roux (Benoit.Roux@med.cornell.edu)
and Guillaume Lamoureux (Guillaume.Lamoureux@umontreal.ca)
The DRUDE command generates a polarizable system by modifying the
topology and parameters of an existing non-polarizable system. For
each selected atom, it creates a "Drude oscillator" 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 has a charge 'Q', unchanged from
the partial charge the non-polarizable atom had prior to calling the
DRUDE command. 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 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 read from WMAIN, and converted
into charges 'q', assuming a force constant 'k = 2*KDRUDE'.
See J. Chem. Phys. 119, ??? (2003) 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. (This
is done by creating additional fictitious bonds of force constant zero
between the particles.)
For a single atom (charges in parenthesis):
DRUDE (q)
A --------> A~DA
(Q) (Q-q)
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
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.
* Menu:
* 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
Syntax of the DRUDE command
DRUDe RESEt
DRUDe [MASS real] [KDRUde real] -
[THOLe real [SHAPe integer]] -
atom-selection
atom-selection::= (see *note select:(select.doc).)
Description of the DRUDE command
-------------------------------------------------------------------
Keyword Default Purpose
-------------------------------------------------------------------
RESET Desactivates the Drude particles. After a reset,
the user should delete all Drude particles.
MASS 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.
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.
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.
SHAPE 1 Specifies the shape of the dipole-dipole
screening function.
-------------------------------------------------------------------
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
is overriding any bond constant found in the parameter file. For
highly polar molecules like water, the recommended value for KDRUDE is
500 kcal/mol/Angst**2.
The atomic polarizabilities (in Angst**3) are 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.
-------------------------------------------------------------------
2) THOLE, SHAPE
(For 1-2 and 1-3 atoms only. All other interactions are regular
Coulomb.)
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.
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 where
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.
Effect on the topology and parameters
-------------------------------------------------------------------
1) New atoms
The Drude particles are inserted immediately after their corresponding
atoms. 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 to
maintain between the atom-Drude pairs the same 1-2, 1-3, and 1-4
relationships that were existing previously to the DRUDE call. 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
Weither 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).
To be aware of when using the DRUDE command
-------------------------------------------------------------------
1) Call the DRUDE command after all the atoms are built
Otherwise, some zero-force bonds between the Drude particles and
neighboring atoms (as discussed in the previous section) may be
missing. And since bad contacts are difficult to resolve with a
polarizable force field, it is probably safer to minimize/equilibrate
the system using first a non-polarizable force field.
-------------------------------------------------------------------
2) Call the DRUDE command before SHAKE and LONEPAIR
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
-------------------------------------------------------------------
2) 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.
-------------------------------------------------------------------
3) Beware of atom names conflicts
Since the atom names don't have more than four letters, atoms with
different names may end up having Drude particles with the same
name:
C210 --> DC21
C211 --> DC21
-------------------------------------------------------------------
4) MASS and KDRUDE are overriding the toppar files
Any nonzero value for MASS and KDRUDE specified by the user is
overriding the values from the topology and parameter files.
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
The polarizabilities on all carbon atoms are set to 1.5 Angst**3:
SCALAR WMAIN SET +1.5 SELECT .NOT. TYPE H* END
DRUDE SELECT .NOT. TYPE H* END
The selection contains the atoms CG, CD1, CD2, CE1, CE2, and CZ. The
DRUDE command will look for atom types DCG, DCD1, DCD2, DCE1, DCE2,
and DCZ. If these types are unknown, the program will crash. For
this reason, it is necessary to append the atom types of the Drude
particles when reading the topology:
OPEN READ CARD UNIT 1 NAME @TOPPAR/top_all22_drude.inp
READ RTF CARD APPEND UNIT 1
Similarly, the program will crash if bond parameters are missing, and
the additional bond parameters should be appended to the parameters:
OPEN READ CARD UNIT 1 NAME @TOPPAR/par_all22_drude.inp
READ PARAM CARD APPEND UNIT 1
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. Desactivating these interactions by using THOLE=0, the
polarizability tensor is almost isotropic.
-------------------------------------------------------------------
2) SWM4-DP water
(see test/c30test/swm4.inp)
See J. Chem. Phys. ??? (2003) 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
READ COOR CARD ...
SET ALPHAO = 1.042520
SET DOM = 0.238080
SCALAR WMAIN SET @ALPHAO SELECT ( SEGID WAT .AND. TYPE OH2 ) END
DRUDE SELECT ( SEGID WAT .AND. TYPE OH2 ) END
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
LONEPAIR BISECTOR DIST @DOM ANGLE 0.0 DIHE 0.0 -
SELECT ATOM WAT * OM END SELECT ATOM WAT * OH2 END -
SELECT ATOM WAT * H1 END SELECT ATOM WAT * H2 END
The molecular dynamics for polarizable water is explained in
*note vv2:(chmdoc/vv2.doc).
NIH/DCRT/Laboratory for Structural Biology
FDA/CBER/OVRR Biophysics Laboratory
Modified, updated and generalized by C.L. Brooks, III
The Scripps Research Institute