Multipole Module
by Nuria Plattner (nuria_plattner@brown.edu)
and Myung Won Lee (mw.lee@unibas.ch)
and Markus Meuwly (m.meuwly@unibas.ch)
Questions and comments regarding MTP should be directed to
----------------------------------------------------------
Nuria Plattner or Myung Won Lee
References:
N. Plattner and M. Meuwly, Biophys. J., 94, 2505 (2008)
N. Plattner, Distributed multipole moments in atomistic force fields:
implementation and applications, Ph.D. Thesis, University of Basel,
Basel, Switzerland, (2009).
A. J. Stone, The Theory of Intermolecular Forces, Oxford University
Press (1996)
The multipole (MTP) module allows to include electrostatic
interactions between multipolar distributions and point charges
or between several multipolar distributions.
* Menu:
* Interaction:: INTERACTION TYPES INCLUDED
* Input:: MTP INPUT FILES
* Parameters:: HOW TO OBTAIN ATOMIC MULTIPOLE PARAMETERS
* Gradients:: GRADIENTS
* Potentials:: COMBINATION WITH ANHARMONIC BOND POTENTIALS
INTERACTION TYPES INCLUDED -------------------------- Currently interactions up to Rank 2 (quadrupole) are fully implemented and interactions up to Rank 3 (octopole) are included for linear molecules. In detail, the following components are available: 3 components for dipole moments 5 components for quadrupole moments 1 component for octopole moments The interactions of these components include interactions of all components with point charges as well as the interactions of the components with each other. The number of components needed for a given intermolecular interaction depends on the highest multipole rank on each atom as well as on the molecular geometry.
MTP INPUT FILES --------------- A typical input sequence to activate MTP looks as follows: open unit 40 card read name mtp.dma MTP MTPUNIT 40 close unit 40 where mtp.dma is the name of a file containing parameters necessary for MTP (see below). MTP parameters are read initially by the MTP module and used in calculating energies and gradients. A typical MTP parameter file looks like the following: 1 ! number of molecules (a) 1 ! number of types of molecules (a) 2 ! number of sites (b-1) 533 534 0 0 0 ! reference atoms (b-2) QFLC ! and QFLC (b-3) 1.151 ! Re for QFLC (in angstrom) (b-3) NoLINT ! No linear triatomic (b-4) 0 0 ! reference range (b-5) 533 ! reference atom (c) 2 ! rank (Limit 2) 0.0648 ! charge (Q0) 0.3686 ! charge (Q') Q = Q0 + Q' ( R - Re ) 0.5139 ! qu10 (Q0) -1.1332 ! qu10 (Q') 0.3787 ! qu20 (Q0) -0.5204 ! qu20 (Q') 534 ! reference atom 2 ! rank (Limit 2) -0.0648 ! charge (Q0) -0.3686 ! charge (Q') -0.3230 ! qu10 (Q0) 0.9811 ! qu10 (Q') 0.5716 ! qu20 (Q0) -0.1324 ! qu20 (Q') (a) The 1st line specifies the total number of molecules containing MTP, and the 2nd line the number of types of MTP molecules. (b) Parameters for each type of molecule should be given. (b-1) The number of sites is the number of points (atoms) in a molecule that are described by MTP. (b-2) A line with the MTP atom indices should follow, which are atom numbers used in the CHARMM coordinate file. In the current implementation, the total number of MTP atoms in a single molecule cannot exceed 3. Add zeros after the atom numbers on this line so that the total number of integer numbers on this line is 5. (b-3) If the following line contains QFLC, the equilibrium bond length Re should be given in angstrom on the next line and each moment is described by Q = Q0 + Q' ( R - Re ), where R is the interatomic distance, Re is the equilibrium bond length (given after QFLC), Q0 is the moment at R = Re, and Q' is the slope of Q with respect to the change in bond distance. (N.B. Q0 should be in atomic unit, while Re is in angstrom. Therefore, Q' has the unit of Q0 per angstrom.) QFLC is allowed only for diatomic molecules at present. If static moments are to be used, put NoQFLC instead of QFLC. In the case of NoQFLC, equilibrium bond length Re should not be given on the next line. (b-4) Following QFLC/NoQFLC line, LINT/NoLINT should be given, where LINT is used for linear triatomic molecule and NoLINT otherwise. (b-5) Then, the range of atoms should be given. This is used when there are more than one MTP molecule of the same type. The order of atoms in each molecule should be identical. The initial and final atom indices for the molecules of the same type, excluding the first molecule, are provided on this line. The atom indices for the first molecule of this type are given in (c). If there is only a single molecule of this type, just use '0 0'. (c) The parameters for each MTP site in a molecule follow. The atom index is given first, and then the rank is given. Rank 0 is used for charge, rank 1 up to dipole, and rank 2 up to quadrupole moments. In the case of QFLC, Q0 and Q' values for the atom are given for each of the following components: Q00 / Q1Z / Q20 (Q22C) / Q30 for rank 0 / 1 / 2 / 3. Charge for this atom should be set to zero in the psf file. If necessary, Q0 and Q' values for Q22C can be added after Q20. If Q22C is used for a diatomic molecule, it is required to change the 'number of sites' from 2 to 3 and set up a dummy atom as described in (d') below. In the case of NoQFLC, only Q0 values are given for the following components: Q1Z Q1X Q1Y / Q20 Q21C Q21S Q22C Q22S for rank 1 / 2. Charge in the psf file should not be changed for NoQFLC. Here, Q00 means charge, Q1[XYZ] X, Y, and Z components of dipole moments, Q2* quadrupole moments, etc. (d) Part (c) is repeated for each MTP site (atom) in a molecule. (d') When a nonzero Q22C quadrupole moment component is used for a diatomic molecule, a dummy atom is used to set up local coordinate system and information on the dummy atom should be provided. As dummy atom in the MTP module is used only to set up the local axis system of diatomic molecule, the MTPs on the dummy atom are not considered in the calculation of energy. It should be combined with QFLC, and thus Q' values for real atoms are set to 0.0 if fluctuation of MTPs is not desired. Information on dummy atom should be provided after specifying MTPs of all real atoms in a molecule. In the example below, a dummy atom is attached to atom 533, at a distance of 1.5 angstrom from atom 533 with the angle formed by (dummy atom)--(atom 533)--(atom 534) to be 120 degrees. Scan interval should be given in degree. Only integer scan interval is allowed. -533 ! dummy atom attached to atom 533 (N.B. negative sign) 534 1.5 120. 1 ! atom 534, r (in A), theta (in deg), scan_int (in deg) 0 ! rank 0.0 ! charge Q0 for dummy atom (not used) 0.0 ! charge gradient Q' for dummy atom (not used) (e) Parts (b)-(d) are repeated for each type of MTP molecule. More example inputs are provided in the MTP test cases mtp-no-h2o.inp, mtp-h2o-trimer.inp and mtp-cluster.inp
HOW TO OBTAIN ATOMIC MULTIPOLE PARAMETERS
-----------------------------------------
The parameters that should be provided in part (c) can be obtained
from quantum chemical calculations and the GDMA program. A typical
procedure is as follows:
(1) Prepare an input file for Gaussian 03 (G03). Checkpoint file
should be specified. HF, DFT, or MP2 can be used to generate MTP
parameters.
(2) Run the G03 job.
(3) Convert the checkpoint file into a formatted file by 'formchk'
command of G03, which will be used by GDMA.
(4) Prepare GDMA input file. For details, refer to the following web
site:
http://www-stone.ch.cam.ac.uk/documentation/gdma/README.html
(5) Execute GDMA and obtain a punch file. Punch file contains
multipole moment parameters, which can be used in MD simulation.
The multipole scheme used by GDMA and in the MTP module is spherical
tensor notation. If multipoles of another program are to be used or
compared, it should be considered that they may be in Cartesian tensor
notation and have to be transformed to spherical tensor notation
first.
In order to include a new parameter set, it is important to check
whether the order of atom numbers is in agreement with the sign
convention used for the reference axis system. E.g. for linear
molecules, the sign of the dipole and octopole moments taken from
GDMA needs to be adapted to the order of the atoms. For CO, if the C
atom comes first in the atom ordering, the signs are reversed with
respect to the standard orientation in a GAUSSIAN calculation which
puts the oxygen in the positive side of the coordinate system and the
carbon on the negative side of the coordinate system.
GRADIENTS --------- The gradients of static atomic multipole moments are composed of two components: 1) A component depending on the atom position on which the multipole moment is placed. This component is included generally for all multipole interactions in the code. 2) A component depending on the reference axis system in which the orientation of the multipole moment is defined, i.e. on the positions of all atoms defining the corresponding reference axis system. These components have to be coded separately for each new reference axis system. Currently, they are included for linear molecules, linear triatomic and water. For details, see the comments in the code for each multipole interaction subroutine. For the gradients of fluctuating atomic multipole moments (QFLC option), complete gradients include an additional term accounting for the changes of moments with molecular geometry. This contribution has not been incorporated in the MTP module yet.
COMBINATION WITH ANHARMONIC BOND POTENTIALS ------------------------------------------- Besides harmonic potentials for bond stretching and bending, other types of potential may be used, such as Morse potential, RRKR potential, KKY water potential, etc. The potentials mentioned above have been tested (and results from their use are in the literature), but have not been included in the current distribution of CHARMM.
PARALLELIZATION --------------- For each process, separate pair list is built outside MTP module. Energies and gradients due to MTPs are computed for the pair list passed to MTP module in each process, and then are combined outside this module. For a system containing diatomic molecule(s) with non- zero Q22C quadrupole component, only a single processor is supported. Test for the parallel MTP has been made with and without PBCs using an executable built with OpenMPI and PGI compilers.
CHARMM Documentation / Rick_Venable@nih.gov