Spherical Multipole Electrostatic Module using Local
Reference Axis Systems
by Tristan Bereau (bereau@alumni.cmu.edu)
Christian Kramer (christian.kramer@uibk.ac.at)
Markus Meuwly (m.meuwly@unibas.ch)
References:
T. Bereau, C. Kramer, M. Meuwly, submitted
T. Bereau et al., J. Phys. Chem. B _117_ 5460-5471 (2013)
C. Kramer, P. Gedeck, M. Meuwly, J. Comp. Chem. _33_ 1673--1688 (2012)
The MTPL module represents the charge distribution of small (diatomics) to
arbitrarily large molecules using a multipole (MTP) expansion in
spherical harmonics. MTP interactions are computed in the atoms' local
reference axis systems. This reduces the number of MTP interactions to be
evaluated by allowing to set coefficients to zero on the basis of symmetry. The
module computes both interaction energies and forces/torques for molecular
dynamics simulations.
* Menu:
* Coefficients:: MTP coefficients
* Syntax:: Syntax of the MTPL command
* Description:: Description of the keywords and options
* Notes:: General notes
MTP coefficients
MTP interactions up to, and including, rank 2 (quadrupole) are implemented.
All MTP coefficients are expressed in spherical coordinates:
* 1 component for the monopole (i.e., partial charge): Q_00
* 3 components for the dipole: Q_10, Q_11c, Q_11s
* 5 components for the quadrupole: Q_20, Q_21c, Q_21s, Q_22c, Q_22s
The partial charges *replace* the charges read from a topology/PSF file.
The MTPL module computes all interactions *except* the charge-charge
interactions, which need to be taken care of by other modules (e.g., PME).
Any atom can be given an MTP rank between 0 (i.e., monopole) and 2 (i.e.,
quadrupole). Different ranks on different atoms can be used within the same
simulation.
MTP coefficients must be expressed in a modified PUN file that describes both
the MTP parameters on every atomic site and the associated local axis
system. All atoms must be specified in the same file, using the same atom ID
(and order) as in the psf file. The structure of the file consists is:
1. 3 lines for title/comments (unused)
2. atom ID; position x; position y; position z; "Rank"; rank
3. "LRA:"; lra; neighbor 1; neighbor 2; neighbor 3; neighbor 4
4. Q_00
5. Q_10; Q_11c; Q_11s
6. Q_20; Q_21c; Q_21s; Q_22c; Q_22s
7. *empty line*
8. *Repeat from step 2. to include more atoms*
where ";" corresponds to a blank space; 'lra' is the type of local reference
axis system, which can be one of the following: 'c3v', 'int', 'ter' (see
[Kramer, Gedeck, Meuwly, JCC _33_ (2012)]) as well as 'lin' for linear/diatomic
molecules (all of them without the single quotes); neighbor # is the atom ID of
a neighbor used to define the local axis system ('0' means no neighbor); and
Q_xx are the spherical MTP coefficients expressed in the local axis system.
The order and number of neighbors for each type of axis system is described in
Kramer et al., except for 'lin', where an atom has one neighbor if it's at a
terminal position and two neighbors (with the same priority rules as in Kramer
et al.) if it's an internal atom. The positions x, y, and z are not used in the
present implementation. The atom ID starts at 1 for the first particle, not 0.
The following is an example of such a file for a single water molecule
# LPUN file for a single water molecule
# MTP interactions expressed in the local axis system
#
1 O2HH 1.418 2.967 2.166 Rank 2
LRA: int 3 2 0 0
-0.4188
0.00 -0.00 -0.4358
-0.9632 0.00 0.00 0.4747 0.00
2 HO2H 2.185 3.54 2.157 Rank 2
LRA: ter 1 3 0 0
0.2094
-0.0279 0.00 0.00
0.17852 -0.0132 0.00 0.0278 0.00
3 HO2H 0.677 3.56 2.287 Rank 2
LRA: ter 1 2 0 0
0.2094
-0.0279 -0.00 0.00
0.17852 -0.0132 0.00 0.0278 0.00
Further assistance on the generation of this file will be provided in an
upcoming publication [Kramer, Bereau, Spinn, Liedle, Gedeck, Meuwly, in
preparation].
Syntax of the MTPL command
MTPL MTPUNIT integer [PREF real] [cutoff-spec]
cutoff-spec::= [RON2 real] [ROFF2 real] [RON3 real] [ROFF3 real]
[RON4 real] [ROFF4 real] [RON5 real] [ROFF5 real]
Description of the keywords and options
Keyword default Purpose
MTPUNIT read above-mentioned local PUN file (card format)
PREF 1.0 prefactor that scales all interaction energies and
forces/torques by PREF (needs to be between 0.0 and 1.0).
RON2 CTONNB Threshold value of the switching function for all
interactions with power law ~1/R^2.
RON3 CTONNB id. for ~1/R^3.
RON4 CTONNB id. for ~1/R^4.
RON5 CTONNB id. for ~1/R^5.
ROFF2 CTOFNB Cutoff value of the switching function for all interactions
with power law ~1/R^2.
ROFF3 CTOFNB id. for ~1/R^3.
ROFF4 CTOFNB id. for ~1/R^4.
ROFF5 CTOFNB id. for ~1/R^5.
General notes
* The use of the present module requires CHARMM to be compiled with the MTPL
flag.
* Parallelization (through MPI) is supported.
* The MTPL module converts the torques generated by the interactions into forces
on the neighboring atoms (the ones involved in an atom's local axis system).
For stability reasons, we advise the use of bond-constraint algorithms (e.g.,
SHAKE) on all bonds to hydrogens.
* The PREF option can show useful when one wishes to perform thermodynamic
integration. The interaction potential can simply be scaled to the coupling
constant linearly:
U_TI(q; lambda) = lambda * U_MTP(q)
where U_MTP(q) is the original interaction energy, q is the set of
coordinates, and lambda is the coupling parameter. In that case, the
simulation is run with PREF=lambda, and a later postprocessing read of the
trajectory can extract the energies with the original Hamiltonian (i.e.,
PREF=1.0), since the derivative of the interaction energy with respect to
lambda will yield the original interaction energy U_MTP(q).
* Increasing the PRNLEV to at least 5 will output the total MTP energy computed
by the module at every MD step.
* Compiling CHARMM with the MTPL_DEBUG flag will show detailed information about
the calculation of every energy and force/torque.
CHARMM Documentation / Rick_Venable@nih.gov