dcm (c46b2)

Distributed Charge Model (DCM)

By Mike Devereux, Oliver Unke, Eric Michoulier, Markus Meuwly
(Michael.Devereux@unibas.ch, m.meuwly@unibas.ch )


The DCM command allows the use of the Distributed Charge Model to
compute the electrostatic energy and force. The method is documented
in the article 'A Novel, Computationally Efficient Multipolar Model
Employing Distributed Charges for Molecular Dynamics Simulations'
(J. Chem. Theory Comput., 2014, 10, pp
4229-4241. http://dx.doi.org/10.1021/ct500511t). If called, the DCM
module zeroes all atomic charges in the standard CHARMM 'CG' array,
copies their original values to a new 'DCG' array, adds off-center DCM
charges and then handles all electrostatics internally. The
electrostatic term in the output file therefore, arises purely from
the DCM module. The module is also compatible with charges fitted in
line with the "minimal distributed charge model" (MDCM) extension to
the original method, conceived to reproduce the electrostatic
potential (ESP) generated by a molecule with similar accuracy using
fewer charges. This extension is documented in (J. Chem. Phys., 2017,
147, 161712, http://dx.doi.org/10.1063/1.4993424).


* Syntax | Syntax of DCM Commands
* Description | Description of DCM Commands
* Requirements | Requirements to run DCM calculations
* Capabilities | Description of capabilities
* DCM format | Description of the .dcm topology file
* Examples | Input Files for various examples


Top
Syntax of DCM Commands

[Syntax DCM]

DCM { IUDCm elec-opt-spec [ XYZ IUXYZ ] [ MUL IUMUL ] [ LAMBda real ] }

DCM { S1 IUS1 }


elec-opt-spec::= [ SWITCH ] [ SHIFT ] [ FSHIFT ] [ TSWITCH ] [ TSHIFT ]


Top
Keyword Default Purpose

IUDCm - Fortran unit of .dcm parameter file
(should be opened in input file before using)

LAMBda - Scaling factor to be applied to DCM energy for each energy
evaluation. This option is used for evaluating free energy
contributions in line with the approach of Bereau et al. (Bereau, T.,
Kramer, C., Meuwly, M., J.Chem.Theor.Comput., 2013, 9, 5450-5459)

XYZ - Write Cartesian coordinates of all DCM charges along with
nuclear coordinates to a file for diagnostic purposes or for use in
QM/MM calculations

IUXYZ - Fortran unit of file where all coordinates and charges should
be dumped at every timestep

MUL - Calculate and write atomic multipoles from DCM charge
arrangements of each atom in global reference axis system in GDMA's
'.lpun' format (diagnostic)

IUMUL - Fortran unit of file where total atomic multipole moments
should be dumped

S1 - For thermodynamic integration calculations in
combination with CHARMM's "PERT" keyword only. Defines the "Lambda=1"
product electrostatic state, note that this keyword assumes that the
"Lambda=0" or "initial" TI electrostatic state has already been
defined and stores all existing DCM charges before the "S1" option is
called in a separate array. See test case.

IUS1 - Fortran unit of .dcm parameter file containing DCM charges for
final (product) state


Electrostatic cutoff options:

The standard SWITCH, SHIFT and FSHIFT cutoffs have been implemented
but need to be declared separately in the DCM section, overriding
whatever was requested by NBOND. The cutoff values (inner and outer
cutoff, nonbond list cutoff and image cutoff) will be taken from the
NBOND section of the standard CHARMM input file, so only the desired
method needs to be declared in the DCM section.

Note that "SHIFT", "FSHIFT" and "SWITCH" cutoffs have been implemented
such that all DCM charges belonging to the same atom experience the
same damping (the charge-charge distance evaluated is actually the
nucleus-nucleus distance). This can be important for true DCM models
derived directly from multipole expansions, but is not appropriate for
MDCM, TIP4P or TIP5P, where charges are further from the nucleus and
the models were fitted by damping based on the true charge-charge
separation, rather than the nucleus-nucleus separation. For such
models TSHIFT (recommended) or TSWITCH (not recommended) can be used.

SWITCH - "SWITched" cutoffs will be used, nuclear coordinates are
used to approximate charge-charge distance. See Nbonds.

SHIFT - "SHIFted" cutoffs will be used, nuclear coordinates are used
to approximate charge-charge distance. See Nbonds.

FSHIFT - "FSHIfted" cutoffs will be used, nuclear coordinates are
used to approximate charge-charge distance. See Nbonds.

TSWITCH - Not recommended. "SWITched" electrostatic cutoffs will be
used, charge-charge distances are used to evaluate cutoffs. See
Nbonds.

TSHIFT - Recommended. CHARMM "SHIFted" electrostatic cutoffs will be
used, charge-charge distances are used to evaluate cutoffs. See
Nbonds.


Top
Requirements

A DCM calculation requires a 'DCM' section in the CHARMM input file
and an external file (typically with extension '.dcm') containing the
DCM charges for each DCM residue.

The calculation proceeds by first setting up the simulation, defining
atom types, force field parameters, sequence, atomic coordinates
etc. Periodic boundaries and nonbonded cutoffs can also be configured
before calling DCM. Once simulation setup is otherwise complete a call
is made in the input file to the DCM module to add DCM charges to all
atoms with a residue type mentioned in the .dcm parameter file. A call
to 'ENERgy' before and after calling 'DCM' in the CHARMM input file
should reveal a change in the electrostatic interaction energy after
DCM initialization. After initializing the DCM module it will be
called each time the energy is updated.

If DCM is to be used in conjunction with CHARMM's "PERT" module for
thermodynamic integration, then the initial lambda=0 state should be
set up before calling DCM with the corresponding lambda=0 .dcm
parameter file, and then "PERT" to store the initial state. The system
should then be modified as appropriate to form the lambda=1 final
state and the DCM command should be called again with the "S1" option
to read in a second "lambda=1" .dcm parameter file prepared by the
user to define the final state. Dynamics for the given TI window can
then begin. See test case.


Top
Capabilities

The DCM code has been integrated with standard CHARMM machinery for
periodic boundaries, non-bonded cutoffs and constant-pressure
simulations. Ewald summation is not yet supported.

There is no requirement for all moieties in a system to carry DCM
charges, akin to QM/MM simulations. A standard CHARMM simulation can
be setup with only certain residues of interest assigned DCM charges,
any atom-centered charges will be passed to the DCM module and the
split-level atom-centered/DCM simulation will proceed automatically.

The DCM charges and Cartesian coordinates can be dumped to a file at
every timestep using the 'XYZ' command in the DCM block of the input
file. This file can be used diagnostically or parsed and used to pass
the DCM charges as external point charges in a QM/MM calculation. It
is not yet possible to specify how frequently the file should be
updated, or to specify a subgroup of the atoms to be dumped.

The code works for diatomics so long as all DCM charges lie along the
bond axis, as the bond is used to define the z-axis while the x- and
y-axes are undefined for diatomics.

Free Energy Calculations have been implemented using CHARMM's PERT
routines. To use this feature, add a standard DCM section before the
PERT command is called that reads a .dcm parameter file for the
LAMBDA=0 state. Call PERT as usual for a standard CHARMM TI run. After
PERT you should define your LAMBDA=1 state by calling the DCM routine
again, but this time read in a second .dcm file with the keyword 'S1'

A simple polarization routine is included in the DCM module, where an
isotropic polarizability can be assigned to each atom in the dcm
charge parameter file. To maintain computational efficiency a
single-iteration approximation is used, i.e. rather than the induced
dipole on an atom modifying the electric field of other atoms and
inducing a new dipole, iterating until self-consistency, we stop after
a single iteration. This is equivalent to assuming that induced
dipoles are caused only by the permanent electric field of surrounding
moieties, and are not affected by the surrounding induced dipoles. As
such, the polarizability should be treated as a parameter to be
fitted, rather than a physically measurable quantity. The polarization
catastrophe is avoided using a damping function employed in the AMOEBA
force field (J Phys Chem B. 2010 March 4; 114(8):
2549–2564. doi:10.1021/jp910674d ).

The maximum number of charge sites per atom is 7.


SHAKE

Use of SHAKE is supported, but is not yet compatible with the "FAST
WATER" option. Use instead the older (and more general) "SHAKE BOND"
approach.


Top
DCM file format

The DCM parameter file (extension .dcm ) will load the following
information / parameters:

NR: The number of residues with DCM charges and polarizabilities. For
solute in solvent NR=2.

NAME: The name of the residue as in the psf file.

NFRAMES: The total number of axis system frames required for the
residue, see J. Chem. Theory Comput., 2014, 10 (10), pp 4229–4241 for
a description of the local axis system and how to define them.

Na Nb Nc SAXIS: Atom indices to define the axis system frames and the
type of axis system to use for atoms Na, Nb, Nc. SAXIS is the type of
axis system. Two keywords exist BO and BI corresponding to bond or
bisector z-axis. BO will use atom a and b to define the local z axis,
BI uses the bisector of a b c but is not supported in this release.
Note that in the line 'Na Nb Nc', if the third atom index Nc is set to
'0' the molecule is taken to be a diatomic and only two atoms
(i.e. information for atoms a and b) are expected to follow.

NCa NPa: the number of charges and polarizable sites for atom a

POSX POSY POSZ CV: x, y, z positions of the charge sites with the atom
at (0 0 0) and the magnitude of the charge

PVa: value of isotropic polarizability for atom a

NOTES:

Identical blocks follow for atoms b and c. If an atom index occurs in
multiple frames only one such block is given and in all other frames
the number of charges and polarizable sites is set to "0 0", see
example below.

For each axis system frame defined in a DCM parameter file, three
local axis systems (one for each atom) are constructed with the
origins taken at the position of each atomic nucleus. The number of
charges and polarizable sites per atom (local axis system) and the
position (in Angstrom) of each charge in its local axis system are
also defined, along with the corresponding values of charges and
polarizabilities (optional) in (a.u).

To summarize, once an axis system frame for atoms a,b,c is defined,
the number, (x,y,z)-position and magnitude of charges (and polarizable
sites if requested) for atom a in its local axis system, then for atom
b and then the same for atom c are defined next.

The .dcm file has a fixed format, see examples.

If an atom appears in more than one axis system frame, then in the
.dcm file, the number of charges and charge values for this same atom
have to be declared multiple times, but the charges will be non zero
in only one axis system frame. See also the explanation in the first
paragraph, page 4232 of J. Chem. Theory Comput., 2014, 10 (10), pp
4229–4241 . See example below.


Top
Examples

!----------------------------------EXAMPLE H2O--------------------------------!
The following is a .dcm file for water, which provides a simple example as
there is only one frame and one off-center charge per atom. All atoms are
polarizable.
!-----------Next line is begining of the .dcm file for TIP4P------------------!
1 ! Number of residue types defined here

WAT ! residue name
1 ! Number of axis system frames
1 2 3 BO ! atom indices involved in frame 1
1 1 ! number of charges and polarizable sites for atom 1
-0.0288081577 -0.0000000000 -0.1785141205 0.9378318804
0.3047450322
1 1 ! number of charges and polarizable sites for atom 2
-0.2071277011 -0.0000000000 0.1603893325 -1.8756637608
0.7819236617
1 1 ! number of charges and polarizable sites for atom 3
0.0288081601 0.0000000000 -0.1785141067 0.9378318804
0.3047450322

! Note that in the line '1 2 3 BO', if the third atom index is set to
! '0' the molecule is taken to be a diatomic and only two atoms (with
! up to 7 DCM charges each) are expected to follow.
!-----------------------end of the .dcm file for TIP4P------------------------!


!----------------Example .dcm file for Chlorobenzene and water--------!
2 ! Number of residue types defined here

DCLB
6 ! Number of axis system frames
12 11 9 BO ! Atom indices involved in this frame and axis type
6 0
0.1322883528 0.0000000000 -0.0012546326 -2.4350397369
-0.1322883528 0.0000000000 0.0012546326 -2.3922930817
0.0000000000 -0.1322943021 0.0000000000 0.0000000000
0.0012546326 0.0000000000 0.1322883528 14.5389393293
-0.0012546326 0.0000000000 -0.1322883528 15.9598702605
0.0000000000 0.0000000000 0.0000000000 -25.7497399110
6 0
0.1322666498 0.0000000000 -0.0027047604 -3.3405778203
-0.1322666498 0.0000000000 0.0027047604 -3.2891731232
0.0000000000 -0.1322943021 0.0000000000 0.0000000000
0.0027047604 0.0000000000 0.1322666498 4.7843285906
-0.0027047604 0.0000000000 -0.1322666498 5.0234232554
0.0000000000 0.0000000000 0.0000000000 -3.2386563845
6 0
-0.0670979693 0.0000000000 0.1140159853 -1.7871855727
0.0670979693 0.0000000000 -0.1140159853 -1.0176053084
0.0000000000 -0.1322943021 0.0000000000 0.0000000000
-0.1140159853 0.0000000000 -0.0670979693 6.1718754315
0.1140159853 0.0000000000 0.0670979693 4.8658246976
0.0000000000 0.0000000000 0.0000000000 -8.6473962597
10 9 7 BO
6 0 ! charge positions and magnitudes for atom 10
0.0000000000 -0.1322943021 0.0000000000 -2.0665148368
0.0000000000 0.1322943021 0.0000000000 -2.0665148368
-0.1314857660 0.0000000000 -0.0146039624 0.0267284440
0.0146039624 0.0000000000 -0.1314857660 4.8178288032
-0.0146039624 0.0000000000 0.1314857660 2.3956766186
0.0000000000 0.0000000000 0.0000000000 -2.5726735323
0 0 ! charge positions and magnitudes for atom 9 already given above!
6 0 ! charge positions and magnitudes for atom 7
-0.0555789916 0.0000000000 0.1200531468 -1.6473020243
0.0555789916 0.0000000000 -0.1200531468 -2.1338069610
0.0000000000 -0.1322943021 0.0000000000 0.0000000000
-0.1200531468 0.0000000000 -0.0555789916 9.8311970287
0.1200531468 0.0000000000 0.0555789916 13.1161265517
0.0000000000 0.0000000000 0.0000000000 -18.9810137226
8 7 5 BO
6 0
0.0000000000 -0.1322943021 0.0000000000 -0.7955207897
0.0000000000 0.1322943021 0.0000000000 -0.7955207897
-0.1322872437 0.0000000000 0.0013665776 0.0524478167
-0.0013665776 0.0000000000 -0.1322872437 7.5205576582
0.0013665776 0.0000000000 0.1322872437 6.5326075391
0.0000000000 0.0000000000 0.0000000000 -12.6624181512
0 0
6 0
0.0000000000 -0.1322943021 0.0000000000 -1.2344997536
0.0000000000 0.1322943021 0.0000000000 -1.2344997536
0.1026663708 0.0000000000 0.0834349968 -1.2668789479
-0.0834349968 0.0000000000 0.1026663708 1.9634869222
0.0834349968 0.0000000000 -0.1026663708 1.8207392319
0.0000000000 0.0000000000 0.0000000000 -0.7470696934
6 5 3 BO
6 0
-0.1321728140 0.0000000000 0.0056683007 -1.5844020995
0.1321728140 0.0000000000 -0.0056683007 -1.4546973784
0.0000000000 -0.1322943021 0.0000000000 0.0000000000
-0.0056683007 0.0000000000 -0.1321728140 3.2681447627
0.0056683007 0.0000000000 0.1321728140 1.5583664720
0.0000000000 0.0000000000 0.0000000000 -1.2641268969
0 0
6 0
0.0000000000 -0.1322943021 0.0000000000 -1.0643317857
0.0000000000 0.1322943021 0.0000000000 -1.0643317857
0.0708638376 0.0000000000 -0.1117143630 0.0770603839
0.1117143630 0.0000000000 0.0708638376 5.9321977891
-0.1117143630 0.0000000000 -0.0708638376 7.0211108568
0.0000000000 0.0000000000 0.0000000000 -11.0495521750
4 3 1 BO
6 0
0.0000000000 -0.1322943021 0.0000000000 -0.8929191626
0.0000000000 0.1322943021 0.0000000000 -0.8929191626
0.1322222640 0.0000000000 -0.0043652355 0.1409144878
0.0043652355 0.0000000000 0.1322222640 6.6091940574
-0.0043652355 0.0000000000 -0.1322222640 9.5876195788
0.0000000000 0.0000000000 0.0000000000 -14.0828937950
0 0
6 0
0.0000000000 -0.1322943021 0.0000000000 -1.2972793586
0.0000000000 0.1322943021 0.0000000000 -1.2972793586
-0.0021140442 0.0000000000 -0.1322774100 -1.2721486192
0.1322774100 0.0000000000 -0.0021140442 1.7273644808
-0.1322774100 0.0000000000 0.0021140442 1.8110129248
0.0000000000 0.0000000000 0.0000000000 -0.3703920636
2 1 11 BO
6 0
0.0000000000 -0.1322943021 0.0000000000 -1.8934057808
0.0000000000 0.1322943021 0.0000000000 -1.8934057808
0.1320838874 0.0000000000 -0.0074584891 -0.1562765923
0.0074584891 0.0000000000 0.1320838874 2.7844644029
-0.0074584891 0.0000000000 -0.1320838874 5.3044822791
0.0000000000 0.0000000000 0.0000000000 -3.6113278681
0 0
0 0

DCWA
1 ! The number of axis system frames
1 2 3 BO ! Atom indices involved in frame 1
6 0 ! Number of charges for atom 1
0.0000000000 -0.1322943021 0.0000000000 -0.0060093533
0.0000000000 0.1322943021 0.0000000000 -0.0060093533
-0.1040648542 0.0000000000 0.0816840774 0.0325927329
-0.0816840774 0.0000000000 -0.1040648542 0.5076023251
0.0816840774 0.0000000000 0.1040648542 0.5488578071
0.0000000000 0.0000000000 0.0000000000 -0.6978266885
6 0
0.0000000000 -0.1322943021 0.0000000000 -2.8364333031
0.0000000000 0.1322943021 0.0000000000 -2.8364333031
-0.1046038419 0.0000000000 0.0809927073 -0.5067793796
-0.0809927073 0.0000000000 -0.1046038419 2.7861557042
0.0809927073 0.0000000000 0.1046038419 2.7861544758
0.0000000000 0.0000000000 0.0000000000 -0.1510791340
6 0
0.0000000000 -0.1322943021 0.0000000000 -0.0060085640
0.0000000000 0.1322943021 0.0000000000 -0.0060085640
0.1040623245 0.0000000000 0.0816873001 0.0325939203
-0.0816873001 0.0000000000 0.1040623245 0.5488583230
0.0816873001 0.0000000000 -0.1040623245 0.5076037791
0.0000000000 0.0000000000 0.0000000000 -0.6978314246


!---------------------END of file .dcm---------------------------------------!