Multi-Surface Adiabatic Reactive Molecular Dynamics (MS-ARMD)
by
Tibor Nagy (tibornagy@chem.elte.hu) [coder]
Juvenal Yosa Reyes (juvenal.yosa@unibas.ch)
Markus Meuwly. (m.meuwly@unibas.ch)
RKHS extension by
Oliver Unke (oliver.unke@unibas.ch)
Marco Pezzella (marco.pezzella@unibas.ch)
References:
[1] Multi-Surface Adiabatic Reactive Molecular Dynamics
Tibor Nagy*, Juvenal Yosa Reyes, Markus Meuwly*
submitted to JCTC (Oct, 2013).
[2] State-selected ion-molecule reactions with Coulomb-crystallized
molecular ions in traps
Xin Tong, Tibor Nagy, Juvenal Yosa Reyes, Matthias Germann,
Markus Meuwly*, Stefan Willitsch*
Chemical Physics Letters, Volume 547, 21 September 2012, Pages 1–8
[3] Constructing multidimensional molecular potential energy surfaces from
ab initio data
Timothy Hollebeek, Tak-San Ho, Herschel Rabitz*
The Multi-Surface Adiabatic Reactive Molecular Dynamics (MS-ARMD) method allows
the construction of global reactive potential energy surface from standard
force fields and running molecular dynamics or Monte Carlo simulations on it.
The effective surface is always the lowest energy surface, except for the region
where several surfaces have the same low energy, where it switches smoothly
between them by changing their weights. The algorithm is based on an energy
difference-based switching method, which conserves total energy during dynamics.
The code can be run also with a single state allowing the easy usage of some
advanced parametrization (Morse potential, MIE potential) also for non-reactive
simulations. The CROSS module of CHARMM is based on the ARMD method, which uses
an alternative, time-dependent switching function. The ARMD method has been used
successfully for the simulation of reactions in large molecules (see cross.doc
for references). The comparison of MS-ARMD and ARMD methods and their advantages
features are discussed in detail in reference [1]. Looking at the example
provided as a test case with CHARMM (mrmd_h2so4.inp with mrmd_h2so4.par MRMD
parameter file, see Example Section) can help in understanding this
documentation greatly.
* Menu:
* Syntax::
Syntax of the MRMD command
* Description::
Description on the MRMD command
* Parameter file::
Description of multiple states and their connection
* States::
Declaration of states: leveling and their connection
* Nonbonded::
(Re)parametrization of nonbonded interactions for
each state
* Bonded::
(Re)parametrization of bonded interactions for each
state
* Output::
Detailed description of output generated during normal
termination
* Example::
Water elimination from vibrationally highly
excited sulfuric acid molecule
* Troubleshooting::
Hint for troubleshooting
* Code::
Essential code related notes
Syntax for the Reactive Molecular Dynamics commands
|---------------optional arguments---------------|
MRMD UPAR integer [ UCRG integer ] [ PRDY integer ] [ PRCA integer ]
MRMD RSET
UPAR Int Unit containing the mrmd parametrization for all surfaces.
This unit must be opened (FORMatted) for reading before MRMD
is called.
UCRG -1 Unit to which a geometry will be written (FORMatted write) in PDB
format at each crossing point (default=-1 => no writing).
PRCA -1 Period in module calls at which MRMD prints out the energy
and the weight of each surface and the energy of the effective
surface when dynamics is not active (default=-1 => no writing).
PRDY -1 Period in time-steps at which MRMD prints out the energies
and the weights of each surface and the energy of the effective
surface during dynamics. (default=-1 => no writing)
RSET Deactivate currently active MRMD parametrization. Should not be
used together with other MRMD keywords.
Description of the Reactive Molecular Dynamics command
-----------------------------
MRMD command in CHARMM input:
-----------------------------
MS-ARMD force field is invoked by the MRMD command. On the first call the
content of MRMD parameter file is fully processed. This contains information on
one or more force fields by providing the changes from the PSF/parameter file
by means of modifying/removing/adding FF potential terms. The PSF and parameter
arrays in CHARMM are never modified! MRMD will affect only the total energy
and the forces on individual atoms. The parameter file also contains information
on the switching function, that is how the surfaces are connected.
----------------------------------------------------------------------------
Position of the MRMD command during input relative to other CHARMM commands:
----------------------------------------------------------------------------
Before MRMD is called, it is important that the bonded parameter list
is up-to-date, therefore it is recommended to execute the UPDAte
command before calling MRMD. It requires that the PSF is already
built.
Before executing the DYNA command it is important that crossing
geometry (UCRG) file is opened.
All commands (MINI, VIBR, DYNA) which are supposed to use the force
field parametrization provided by MRMD should be preceded by
calling MRMD.
A typical input sequence for an MRMD simulation is as follows:
"
...read or generate PSF...
UPDATE
OPEN UNIT 9 READ FORMATTED NAME mrmd.par
OPEN UNIT 14 WRITe FORMatted NAME crossings.pdb
MRMD UPAR 9 UCRG 14 PRCA 10 PRDY 1000
MINI ...
OPEN UNIT 12 WRITe FORMatted NAME new.res
OPEN UNIT 13 WRITe UNFORMatted NAME mrmd_traj.dcd
DYNA ....
"
----------------------------------------------
Compatibility with other major CHARMM commands
----------------------------------------------
- PSF modifying routines -
Calling MRMD assumes that the PSF will not be modified later. If it is
modified then MRMD has to be called again either with the same or a
new parameter file (UPAR).
- VIBRan -
MRMD does not provide (analytical) second derivatives, therefore all
routines which requires second derivatives of the energy can work only
if they have an option for numerical second derivatives (i.e. VIBR
with FINITE).
- IMAGe -
The nonbonded interaction between atoms with modified non-bonded
parameters and image atoms are not updated. If the atoms with modified
nonbonded parameters remain far from the edges of the center box then
this will not cause a problem at all. It is planned to provide this
functionality in the future.
- NBONded -
compatible with NBONded for the following settings:
ELEC : electrostatic interaction
CDIElec: constant dielectric
ATOM : atom-based electrostatic interaction
SHIFt : shifting cut-off method for electrostatic interaction
VDW : van der Waals interaction by Lennard-Jones potential
VSWItch: switching cut-off method for van der Waals interaction
VATOm : atom-based van der Waals interaction
NBXMOD : +1,+2,+3,+4,+5 are all supported
not compatible with NBONded for the following settings:
GROUp : group-based electrostatic interaction
SWITch : switching cut-off method for electrostatic interaction
FSWItch : force-switching cut-off method for electrostatic interaction
VGROup : group-based van der Waals interaction
VSHIft : shifting cut-off method for van der Waals interaction
FVSWitch: force-switching cut-off method for van der Waals interaction
- SKIP -
ARMD is compatible with SKIPe commands
SKIPe BOND : skips bond energy correction defined in BOND HARM
skips bond energy correction defined in BOND MORS
SKIPe ANGL : skips energy correction defined in angle part of ANGL HARM
SKIPe UREY : skips energy correction defined in Urey-Bradley part of ANGL HARM
SKIPe DIHE : skips energy correction defined in DIHE FOUR
SKIPe IMPR : skips energy correction defined in IMPR HARM
SKIPe ELEC : skips energy correction defined in point charge part of NBON ATOM
SKIPe VDW : skips energy correction defined in Lennard-Jones part of NBON ATOM
skips energy correction defined in NBON GVDW
The user has to provide an additional external parameter file beyond the usual CHARMM parameter file to describe one or more potential energy surfaces, which together define a global surface. This (re)parametrization is based on atom indices instead of atom types, therefore allows specific reparametrization of any part of the system. ------------------------------------------------------- Determining effective reactive potential energy surface ------------------------------------------------------- First the energy correction to the CHARMM energy for each surface is calculated based on the (re)parametrization given in this file. Then, the energy shift is added to each surface. Then a smooth connection of the low-lying energy surfaces is calculated using the energy difference based switching functions. Finally, the Gaussian*polynomial functions of pairwise energy differences are added optionally to adjust the shape of the crossing regions between surfaces. The MRMD module in CHARMM corrects the total energy, and does not modify the individual energy terms (BOND, VDW, etc). MRMD energy correction will be added only to the total energy. The parameter file should always have all the following blocks even if the parameters are not changed/used: SURF: defines the number of surfaces and their level shifts SWCH: defines switching function and switching related parameters SHAP GAPO: defines shaping functions to crossing regions between surfaces NBON ATOM: (re)parametrization of charges and Lennard-Jones potentials of atoms NBON GVDW: parametrization of pair-wise generalized Lennard-Jones potentials BOND HARM: (re)parametrization of harmonic bond potentials BOND MORS: parametrization of Morse potentials ANGL HARM: (re)parametrization of harmonic angle potentials DIHE FOUR: (re)parametrization of Fourier dihedral angle potentials IMPR HARM: (re)parametrization of harmonic improper dihedral angle potentials
Description of the individual blocks:
================================================================================
1. Block SURF
================================================================================
ROLE:
Defines the number of surfaces (nsurf) and the energy shifts (VshiftI). At each
call to MRMD, the energy correction of all surfaces are determined, and based
on the energy difference based switching formula, corresponding weights are
determined and the effective correction is calculated.
--------------------------Structure of block SURF-------------------------------
SURF nsurf
1 Vshift1
2 Vshift2
...
I VshiftI
...
nsurf Vshiftnsurf
--------------------------------------------------------------------------------
SHIFTING THE SURFACES (kcal/mol):
dVI(r)=dV0I(r) +VshiftI
VI(r) =VCHARMM(r)+dVI(r)
dVI0:
Energy correction of surface I to the CHARMM energy based on the
(re)parametrization of force field terms given in the MRMD parameter.
dVI:
Total energy correction of surface I to CHARMM energy. Their weighted sum plus
the shaping functions determine the effective correction to the CHARMM energy.
VCHARMM:
CHARMM energy calculated without MRMD module
VI:
Energy of MRMD surface I. Their weighted sum plus the shaping functions
determine the effective global surface.
--------------------------------------------------------------------------------
nsurf ({integer,>0})
Number of surfaces. At least one is required. If only one surface is given then
no surfaces crossing possible.
VshiftI ({real,kcal/mol}):
Energy shift for the I-th surface. It is defined in separate lines for each
surface in increasing order of the surface index. Energy shifting is used to
adjust the level of each surface to a global reference scale (eg. ab initio
energies), which is necessary to reproduce reaction energetics and predict
dividing surfaces between them.
================================================================================
2. Block SWCH
================================================================================
ROLE:
Defines the switching function and switching related parameters.
MS-ARMD switching concept:
1. Normalized energy difference from the lowest-energy surface is determined.
2. Raw weights(wi0) are defined using simple mathematical switching functions.
3. Raw weights are renormalized to give mixing weights (wi).
4. The weighted linear combination of surfaces plus the shaping functions
with the weights give the MS-ARMD surface.
--------------------------Structure of block SURF-------------------------------
SWCH
switching_function deltaV
--------------------------------------------------------------------------------
switching_function (string): mathematical switching function
deltaV({real,>0,kcal/mol}) : switching parameter
Note, input for SWCH has to be provided even if there is only one surface!
MS-ARMD switching concept:
1. Normalized energy difference from the lowest-energy surface is determined:
Vi ! energy of the surface i
Vmin(r) = min(Vi(r),i=1...n) ! energy of the lowest-energy surface
deltai(r) = (Vi(r)-Vmin(r))/deltaV ! normalized energy difference of surface
! i from the lowest-energy surface
2. Raw weights wi0 are defined using switching functions (f).
Currently available switching functions:
A, JOHNSON7: 7th order switching function by B.R.Johnson
B, EXPDECAY: exponential decay switching function (recommended)
w0i(r)=f(deltai(r))
3. Raw weights are renormalized and the linear combination of surface
energies are formed using the renormalized weights.
wi(r) = w0i(r)/(w01(r)+w02(r)+...) ! normalized weights
Veff0(r) = sum(wi(r)*Vi(r),i=1...n) ! effective, smooth surface
Veff0(r) = sum(wi(r)*Vi(r),i=1...n)=VCHARMM+sum(wi(r)*dVi(r),i=1...n)
Veff(r) = Veff0(r)+shaping functions
--------------------------------------------------------------------------------
7th order switching function by B.R. Johnson
B.R. Johnson: J. Chem. Phys. 83, 1204 (1985)
--------------------------------------------------------------------------------
JOHNSON7 delta
--------------------------------------------------------------------------------
w0i(r) = f(deltai(r))=fJ(1-deltai(r))
where:
{ 0 if x<0
fJ(x)={ -x**4*(((20*x-70)*x+84)*x-35) if 0<x<1
{ 1 if x>1
deltaV ({real,>0,kcal/mol}):
if the energy of a surface is less than the energy of the lowest-lying surface
+delta, then it will have a non-zero weight, and start contributing to the
effective surface. When two or more lowest lying surfaces are crossing, their
weights change between 0-1.
NOTE: while the Johnson's switching function is 3 times continuously
differentiable, the final MS-ARMD surface is not differentiable if the
three lowest surfaces are within the deltaV energy of each other.
--------------------------------------------------------------------------------
EXPDECAY switching function
(See JCTC 2013 MS-ARMD paper)
--------------------------------------------------------------------------------
EXPDECAY delta
--------------------------------------------------------------------------------
w0i(r) = exp(-deltai(r))
--------------------------------------------------------------------------------
deltaV ({real,>0,kcal/mol}):
The raw weight of surfaces drops exponentially with increasing energy.
This exponential decay has a characteristic energy of deltaV, that is
an increase in energy by deltaV causes a drop by a factor of
e~=2.7 in the weight. EXPDECAY switching function is recommended over JOHNSON7
switching as it always provides a smooth (analytical, infinite times
differentiable) effective surface.
================================================================================
3. Block SHAP GAPO
================================================================================
ROLE:
Defines the parameters of Gaussian*Polynomial (GAPO) shaping functions to adjust
the crossing region for pairs of surfaces.
--------------------------Structure of block SURF-------------------------------
SHAP GAPO ngapo
surf11 surf12 n1 x10 x11 a10 a11 a12 ... a1n1
surf21 surf22 n2 x20 x21 a20 a21 a22 ... a2n2
...
surfI1 surfI2 nI xI0 xI1 aI0 aI1 aI2 ... aI2nI
...
ngapo...
--------------------------------------------------------------------------------
Energy correction:
Veff(r)=Veff0(r)+sum_I[ {wsurfI1+wsurfI2}*VGAPOI(VsurfI2(r)-VsurfI1(r)) ]
--------------------------------------------------------------------------------
GAUSSIAN*POLYNOMIAL FUNCTION (kcal/mol):
xI=VsurfI2(r)-VsurfI1(r)
VGAPOI(xI) = exp(-(xI-xI0)^2/(2 xI1^2))*
*[aI0+aI1*(xI-xI0)+aI2*(xI-xI0)^2+...+aInI*(xI-xI0)^nI]
where r(angstroem) is the coordinates of all atom.
--------------------------------------------------------------------------------
ngapo ({integer,>=0}):
Total number of GAPO functions used for adjusting the crossing regions.
surfI1, surfI2 ({integer,nsurf>=,>=1}):
The indices of the two surfaces of the I-th GAPO function.
nI ({integer,>=0})
Polynomial order of the I-th GAPO function.
xI0 ({real,1/(kcal/mol)})
Shift of the I-th GAPO function.
xI1 ({real,>0,(kcal/mol)})
Standard deviation parameter of the Gaussian of the I-th GAPO function.
aI0,aI1,...,aInI ({real,(kcal/mol)^(1-j) where j=0,1,2,3})
0-th, 1-st, ..nI-th order polynomial coefficients of the I-th GAPO function.
================================================================================
4. Block NBON ATOM
================================================================================
ROLE:
Reparametrizes LJ interaction and point-charge electrostatic interaction on
each surface.
--------------------------Structure of block ATOM-------------------------------
NBON ATOM natom
atom1 q11 eps11 rmh11 xeps11 xrmh11 q12 eps12 rmh12 xeps12 xrmh12 ...
atom2 q21 eps21 rmh21 xeps21 xrmh21 q22 eps22 rmh22 xeps22 xrmh22 ...
...
atomI qI1 epsI1 rmhI1 xepsI1 xrmhI1 qI2 epsI2 rmhI2 xepsI2 xrmhI2 ...
...
atomnatom ...
--------------------------------------------------------------------------------
ELECTROSTATIC INTERACTION: COULOMB POTENTIAL WITH SHIFTING (kcal/mol):
V(r)=fshift(r)/4pi/EPS0/EPS*q*q/r
where r is the distance of two atoms.
S1(r)=fshift(r)={ 0 if r >roff }
{ (1-(r/roff)^2)^2 if r<=roff }
V(r)=V(r)*E14FAC for atoms in 1-4 position if NBXMOD=1,2,3,5
where roff=CTOFNB
NBXMOD,CTOFNB,EPS,E14FAC are described in nbond.doc
--------------------------------------------------------------------------------
VAN DER WAALS INTERACTION: LENNARD-JONES POTENTIAL WITH SWITCHING (kcal/mol):
eps=sqrt(eps1*eps2)
rmin=rmh1+rmh2
V(r;eps,rmin)=fswitch(r)*eps*[(rmin/r)^12-2*(rmin/r)^6]
where r(angstroem) is the distance of atomI1 and atomI2.
{1 if r<ron }
fswitch(r)={(roff^2-r^2)^2*(roff^2+2r^2-3ron^2)/(roff^2-ron^2)^3 if ron<r<roff}
{0 if roff<r }
where:roff=CTOFNB,ron=CTONNB
CTONNB,CTOFNB are described in nbond.doc
--------------------------------------------------------------------------------
natom ({integer,>=0})
Number of atoms with reparametrized nonbonded interaction.
atomI ({integer,>0}):
PSF index of atom I to be reparametrized, no duplicate definition is
allowed for the same atom.
qIJ ({real, elementary charge unit}):
Charge of atom I on surface J
epsIJ ({real,>=0,kcal/mol}):
Well depth of Lennard-Jones potential for atom I on surface J.
rmhIJ ({real,>0,angstroem}):
Half of the separation corresponding to the minimum of the Lennard-Jones
potential well ("Rmin/2") for atom I on surface J (usual CHARMM convention).
xepsIJ ({x} or {X} or {real,>=0,kcal/mol}):
epsIJ value for "special van der Waals interaction" (CHARMM jargon), which acts
between atoms in 1-4 position and described by Lennard-Jones potential, which is
active in the case of NBXMOD=+5. Opposite to CHARMM convention, this is always
given as non-negative value. If letter x or X is given, then the value copied
from epsIJ.
xrmhIJ ({x} or {X} or {real,>0,angstroem}):
rmhIJ value for "special van der Waals interaction" (CHARMM jargon), which acts
between atoms in 1-4 position and described by Lennard-Jones potential, which is
active in the case of NBXMOD=+5. If letter x or X is given, then the value
copied from rmhIJ.
================================================================================
5. Block NBON GVDW
================================================================================
ROLE:
Adds general-exponent Lennard-Jones potential (equivalent to the Mie-potential)
and removes 6-12 Lennard-Jones potential for pair of atoms on each surface.
General-exponent Lennard-Jones potential can improve the description of the van
der Waals interaction between atoms in the reaction region and it also better
captures the energetics for the reactant and product complexes and around the
transition state.
--------------------------------------------------------------------------------
VAN DER WAALS INTERACTION: GENERAL-EXPONENT LENNARD-JONES POTENTIAL (kcal/mol)
V(r;eps,rmin,rep,att)=att/(rep-att)*eps*[(rmin/r)^rep-rep/att*(rmin/r)^att]
--------------------------Structure of block GVDW-------------------------------
NBON GVDW ngvdw
atom11 atom12 eps11 rmin11 att11 rep11 eps12 rmin12 att12 rep12 ....
atom21 atom22 eps21 rmin21 att21 rep21 eps22 rmin22 att22 rep22 ....
...
atomI1 atomI2 epsI1 rminI1 attI1 repI1 epsI2 rminI2 attI2 repI2 ....
...
atomngvdw1 atomngvdw2 ...
--------------------------------------------------------------------------------
ngvdw ({integer,>=0})
Number of atom pairs with reparametrized van der Waals interaction.
atomI1,atomI2 ({integer,>0}):
PSF indices of atom pair I.
No duplicate definition is allowed for the same pair of atoms (ij=ji).
epsIJ ({x} or {X} or {real,>=0,kcal/mol}):
Well depth of general-exponent Lennard-Jones potential for atom pair I on
surface J. If letter x or X is given for epsIJ, then following values of rminIJ,
attIJ and repIJ are not used, and the 6-12 Lennard-Jones potential will be used
for the corresponding atom pair on surface J.
rminIJ ({real,>0,angstroem}):
Distance of atoms at the minimum of general-exponent Lennard-Jones potential I
on surface J. If letter x or X is given for epsIJ, then following values of
rminIJ, attIJ and repIJ are not used, and the 6-12 Lennard-Jones potential from
the CHARMM parameter file will be used for the corresponding atom pair on
surface J.
attIJ ({real,>0,repIJ>}):
Attractive exponent of general-exponent Lennard-Jones potential I on surface J.
If letter x or X is given for epsIJ, then following values of rminIJ, attIJ and
repIJ are not used, and the 6-12 Lennard-Jones potential will be used for the
corresponding atom pair on surface J.
repIJ ({real,>0,>attIJ}):
Repulsive exponent of general-exponent Lennard-Jones potential I on surface J.
If letter x or X is given for epsIJ, then values of rminIJ, attIJ and repIJ are
not used, and the 6-12 Lennard-Jones potential will be used for the
corresponding atom pair on surface J.
================================================================================
6. Block BOND HARM
================================================================================
ROLE:
Reparametrizes, adds or removes harmonic bond potentials on surfaces.
--------------------------Structure of block HARM-------------------------------
BOND HARM nharm
atom11 atom12 fch11 req11 fch12 req12 ...
atom21 atom22 fch21 req21 fch22 req22 ...
...
atomI1 atomI2 fchI1 reqI1 fchI2 reqI2 ...
...
atomnharm1 atomnharm2 ...
--------------------------------------------------------------------------------
HARMONIC BOND POTENTIAL (kcal/mol):
V(r;fch,req)=fch*(r-req)^2
where r(angstroem) is the length of bond atomI1-atomI2.
--------------------------------------------------------------------------------
nharm({integer,>=0}):
Number of those atom pairs between which the bond is (re)defined.
atomI1, atomI2 ({integer,>0}):
PSF index of the two atoms forming the I-th bond. No duplicate definition is
allowed for the same pair of atoms (ij=ji). The definition for the same pair
of atoms is not allowed in the MORS block.
fchIJ ({x} or {X} or {real,>=0, kcal/mol/angstroem^2}):
Half force constant of bond potential I on surface J (usual CHARMM convention)
If letter x or X is given then the bond is considered as broken in this surface
(even if it is in the PSF file) If letter x or X is given for fchIJ, then the
following value of reqIJ is not used.
reqIJ ({>0,real,angstroem}):
Equilibrium bond length of bond I on surface J. If letter x or X is given for
fchIJ, then the following value of reqIJ is not used.
================================================================================
7. Block BOND MORS
================================================================================
ROLE:
Adds Morse bond potential, removes harmonic bond potential on surfaces.
If the same bond exists in PSF, then it is automatically replaced with this.
--------------------------Structure of block MORS-------------------------------
BOND MORS nmors
atom11 atom12 de11 req11 beta11 de12 req12 beta12 ...
atom21 atom22 de21 req21 beta21 de22 req22 beta22 ...
...
atomI1 atomI2 deI1 reqI1 betaI1 deI2 reqI2 betaI2 ...
...
atomnmors1 atomnmors2 ...
--------------------------------------------------------------------------------
MORSE-POTENTIAL (kcal/mol):
V(r)=de*[1-exp(-beta*(r-req))]^2
where r(angstroem) is the length of bond atomI1-atomI2.
--------------------------------------------------------------------------------
nmors({integer,>=0}):
Number of redefined Morse bonds.
atomI1, atomI2 ({integer,>0}):
PSF indices of atom pair I.
No duplicate definition is allowed for the same pair of atoms (ij=ji).
The definition for the same pair cannot be used in the HARM block.
deIJ ({x} or {X} or {real,>=0, kcal/mol}):
Well depth of Morse bond I on surface J.
If letter x or X is given the bond is broken in this surface (even if it is a
PSF harmonic bond). If letter x or X is given for deIJ, then following values
of reqIJ and betaIJ are not used.
reqIJ ({>0,real,angstroem}):
Equilibrium bond length of Morse bond I on surface J. If letter x or X is given
for deIJ, then following values of reqIJ and betaIJ are not used.
betaIJ ({>0,real, 1/angstroem}):
Beta parameter of Morse bond I on surface J.
If letter x or X is given for deIJ, then folowing values of reqIJ and betaIJ are
not used.
================================================================================
8. Block BOND RKHS (optional if BOMBlev < 5)
================================================================================
ROLE:
Adds RKHS bond potential, removes harmonic bond potential on surfaces.
If the same bond exists in PSF, then it is automatically replaced with this.
For backwards compatibility, this block can be omitted if BOMBlev is smaller than 5.
The Reproducing Kernel Hilber Space (RKHS) is an interpolation scheme which can interpolate a grid of energies. In comparison to for example spline based methods, the gridpoints are necessarily reproduced exactly. An RKHS bond should be used for example if a grid of ab initio energies is available, which is poorly fitted by a Morse potential. Grid data needs to be provided in a simple text file in the following format (without header):
r(angstroem) E(kcal/mol) <-- omit this header line!
r1 E1
r2 E2
. .
. .
. .
ri Ei
. .
. .
. .
rngrid Engrid
Where in each line, ri and Ei are separated by whitespace. It is very important that the energy values of the grid are given such that the energy asymptotically approaches 0 for large values of r. This can be achieved for any kind of grid data by addition of a constant to the energy values.
--------------------------Structure of block RKHS-------------------------------
BOND RKHS nrkhs
atom11 atom12 ngrid11 ugrid11 ucoef11 ngrid12 ugrid12 ucoef12 ...
atom21 atom22 ngrid21 ugrid21 ucoef21 ngrid22 ugrid22 ucoef22 ...
...
atomI1 atomI2 ngridI1 ugridI1 ucoefI1 ngridI2 ugridI2 ucoefI2 ...
...
atomnrkhs1 atomnrkhs2 ...
--------------------------------------------------------------------------------
RKHS-POTENTIAL (kcal/mol):
V(r)=sum[i=1 to Ngrid](ci*k(r,ri))
where r(angstroem) is the length of bond atomI1-atomI2,
Ngrid the number of grid points in the energy grid, ci the
kernel coefficient for grid point i and ri the length of bond
atomI1-atomI2 for grid point i. The kernel function
k(r,r') is described elsewhere [3].
--------------------------------------------------------------------------------
nrkhs({integer,>=0}):
Number of redefined RKHS bonds.
atomI1, atomI2 ({integer,>0}):
PSF indices of atom pair I.
No duplicate definition is allowed for the same pair of atoms (ij=ji).
The definition for the same pair cannot be used in the HARM or MORS block.
ngridIJ ({integer,>0}):
Number of grid points for bond I on surface J. If letter x or X is given for ngridIJ, then following values ugridIJ and ucoefIJ are not used.
ugridIJ ({integer,>0}):
Unit number for the file containing the grid of distance and energy values for bond I on surface J. The file also needs to be opened in the CHARMM input file prior to initializing the MRMD module.
ucoefIJ ({integer}):
Unit number for the file containing the kernel coefficients for bond I on surface J. If a positive value is provided, the coefficients are assumed to be precomputed and read out from unit ucoefIJ. If a negative value is provided, the coefficients are calculated and written to unit -ucoefIJ. If the value 0 is provided, coefficients are calculated and not written to any file. Note that for moderate grid sizes of ngridIJ (< 500) the computation of coefficients is reasonably fast. Reading precomputed values is only interesting if a large grid (ngridIJ >> 1000) is used for many individual CHARMM runs. When in doubt, the user should always use a value ucoefIJ=0. In case a coefficient file should be read or written, the corresponding file also needs to be opened in the CHARMM input file prior to initializing the MRMD module.
================================================================================
9. Block ANGL HARM
================================================================================
ROLE:
Reparameterizes, adds or removes harmonic angle potentials and Urey-Bradley
potentials on surfaces.
--------------------------Structure of block ANGL-------------------------------
ANGL HARM nangl
atom11 atom12 atom13 fch11 phieq11 ufch11 ureq11 fch12 phieq12 ufch12 ureq12 ...
atom21 atom22 atom23 fch21 phieq21 ufch21 ureq21 fch22 phieq22 ufch22 ureq22 ...
...
atomI1 atomI2 atomI3 fchI1 phieqI1 ufchI1 ureqI1 fchI2 phieqI2 ufchI2 ureqI2 ...
...
atomnangl1 atomnangl2 ...
--------------------------------------------------------------------------------
HARMONIC ANGLE POTENTIAL (kcal/mol):
V(phi;fch.phieq)=fch*(phi-phieq)^2
where phi(radian) is the angle formed by atomI1-atomI2-atomI3.
--------------------------------------------------------------------------------
UREY-BRADLEY POTENTIAL (kcal/mol):
V(r;ufch,ureq)=ufch*(r-ureq)^2
where r(angstroem) is the distance between atomI1 and atomI3.
--------------------------------------------------------------------------------
nangl({integer,>=0}):
Number of (re)defined harmonic angle and Urey-Bradley potentials.
atomI1, atomI2, atomI3 ({integer,>0}):
PSF indices of the three atoms forming the I-th angle (atomI1-atomI2-atomI3).
No duplicate definition is allowed for the same three atoms in the give order
(ijk=kji).
fchIJ ({x} or {X} or {real,>=0,kcal/mol/radian^2}):
Half force constant of angle potential I on surface J (usual CHARMM convention).
If letter x or X is given the angle potential is not present on this surface due
to missing bond. If letter x or X is given for fchIJ, then following values of
phieqIJ, ufchIJ, ureqIJ are not used.
phieqIJ ({>0,real,degree}):
After reading it in, it is immediately converted to radian. Equilibrium angle of
angle potential I on surface J. If letter x or X is given for fchIJ, then
following values of phieqIJ, ufchIJ, ureqIJ are not used.
ufchIJ ({x} or {X} or {real,>=0,kcal/mol/angstroem^2}):
Half force constant of Urey-Bradley potential I on surface J
(usual CHARMM convention).
If letter x or X is given the Urey-Bradley potential is not present on this
surface and the following values of phieqIJ, ufchIJ, ureqIJ are not used.
ureqIJ ({real,>0,degree}):
Equilibrium distance of Urey-Bradley potential I on surface J. If letter x or X
is given for fchIJ, then following values of phieqIJ, ufchIJ, ureqIJ are not
used.
================================================================================
10. Block DIHE FOUR
================================================================================
ROLE:
Reparameterizes, adds or removes proper dihedral potentials on surfaces.
--------------------------Structure of block DIHE-------------------------------
DIHE FOUR ndihe
atom11 atom12 atom13 atom14 per1 amp11 phi011 amp12 phi012 ...
atom21 atom22 atom23 atom24 per2 amp21 phi021 amp22 phi022 ...
...
atomI1 atomI2 atomI3 atomI4 perI ampI1 phi0I1 ampI2 phi0I2 ...
...
atomndihe1 atomndihe2 atomndihe3 atomndihe4 ...
--------------------------------------------------------------------------------
PROPER DIHEDRAL POTENTIAL (kcal/mol):
per=1,2,..,6 => V(phi)=amp*[1+cos(per*phi-phi0)]
per=0 => V(phi)=amp*(phi-phi0)^2
where phi(radian) is the dihedral angle formed by atomI1-atomI2-atomI3-atomI4.
Expected bond connectivity of atoms: atomI1-atomI2-atomI3-atomI4
--------------------------------------------------------------------------------
ndihe({integer,>=0}):
Number of (re)defined proper dihedral potentials.
atomI1, atomI2, atomI3, atomI4 ({integer,>0}):
PSF indices of the four atoms forming the I-th dihedral angle. No duplicate
definition is allowed for the same 4 atoms in the given order (ijkl=lkji)
with the same periodicity.
perI ({integer:1,2,3,4,5,6}):
Periodicity of cosine dihedral potential I. For the same group of atoms multiple
declarations with different periodicity are allowed (Fourier series:1,2,3,4,5,6)
No duplicate definition is allowed for the same 4 atoms in the given order
(ijkl=lkji) with the same periodicity.
ampIJ ({x} or {X} or {real,kcal/mol},{real,>0,kcal/mol}):
Amplitude of dihedral potential I on surface J. If letter x or X is given for
ampIJ, then the dihedral potential does not exist on the corresponding surface
and the following value of phi0IJ is not used.
phi0IJ ({x} or {X} or {real,degree}):
Zero phase for dihedral potential. After reading it in, it is immediately
converted to radian within the code.
The location of extrema are:
if ampI>0:
cos(n*phi_max-phi0) maximal <=> n*phi_max-phi0=2*k*pi k is integer
=>maxima: phi_max=(2*k*pi+phi0)/n k is integer
cos(n*phi_min-phi0) minimal <=> n*phi_min-phi0=(2*k+1)*pi k is integer
=>mixima: phi_min=((2*k+1)*pi+phi0)/n k is integer
if ampI<0:
-cos(n*phi_min-phi0) minimal <=> n*phi_min-phi0=2*k*pi k is integer
=>minima: phi_min=(2*k*pi+phi0)/n k is integer
-cos(n*phi_max-phi0) maximal <=> n*phi_max-phi0=(2*k+1)*pi k is integer
=>maxima: phi_max=((2*k+1)*pi+phi0)/n k is integer
================================================================================
11. Block IMPR HARM
================================================================================
ROLE:
Reparameterizes, adds or removes improper dihedral potentials on surfaces.
--------------------------Structure of block IMPR-------------------------------
IMPR HARM nimpr
atom11 atom12 atom13 atom14 per1 fch11 phi011 fch12 phi012 ...
atom21 atom22 atom23 atom24 per2 fch21 phi021 fch22 phi022 ...
...
atomI1 atomI2 atomI3 atomI4 perI fchI1 phi0I1 fchI2 phi0I2 ...
...
atomnimpr1 atomnimpr2 atomnimpr3 atomnimpr4 ...
--------------------------------------------------------------------------------
IMPROPER DIHEDRAL POTENTIAL (kcal/mol):
perI=0 V(phi)=fch*(phi-phi0)^2
where phi(radian) is the dihedral angle formed by atomI1-atomI2-atomI3-atomI4.
Expected bond connectivity of atoms:
atomI1 - atomI4
| \
atomI2 atomI3
IMPORTANT:
CHARMM topology should use the same topology definition for improper dihedrals
in the reaction center otherwise dihedrals will not be successfully removed
on MRMD surface if a bond gets broken. Improper dihedrals to be removed or
reparameterized have to have the same or reverse atom order.
------------------------------------------------------------------------------
nimpr({integer,>=0}):
Number of (re)defined improper dihedral potentials.
atomI1, atomI2, atomI3, atomI4 ({integer,>0}):
PSF indices of the four atoms forming the I-th improper dihedral angle.
No duplicate definition is allowed for the same group of atoms in the
given order (ijkl=lkji)
perI ({integer:0}):
Periodicity of dihedral potential I. Only perI=0 is allowed and it implies
quadratic angle potential.
fchIJ ({x} or {X} or {real,kcal/mol/radian^2}):
Half force constant of quadratic dihedral potential I on surface J (usual CHARMM
convention). If letter x or X is given for fchIJ, then the improper potential
does not exist on the corresponding surface and the following value of phi0IJ is
not used.
phi0IJ ({x} or {X} or {real,degree}):
Equilibrium angle of quadratic dihedral potential I on surface J. After reading
it in, it is immediately converted to radian within the code.
Output to standard output: ---------------------------------------------- 1. At the call of MRMD command in CHARMM input ---------------------------------------------- Subroutine MRMD_INIT will be executed, which will print: Printed lines start with: 'MRMD_INIT>...' if PRNLEV>=5 the following will be printed: - the arguments with which the MRMD command was called - progress of reading various blocks of MRMD parameter file - after reading the parameter file, the number of records in each block if PRNLEV>=6 the following will be printed: - all the parameters that are read from the MRMD parameter file. - changes in 1-2, 1-3, 1-4, 1-(more than 4) neighbourship lists for each surface ----------------------------------------------------- 2. At energy call if certain conditions are fulfilled ----------------------------------------------------- Subroutine EMRMD will be executed, which calls several other subroutines. Printed lines start with: 'EMRMD>...' 'MRMD_ENERG>...' 'MRMD_SURF>...' ------------------------------- 2A.Printing of surface crossing: ------------------------------- - PRNLEV>=4: time of crossing and the corresponding two surfaces --------------------------------------------- 2B.Printing energetics at various detail level --------------------------------------------- - During dynamics after every integration step defined after PRDY argument. - When dynamics is not active after every PRCA calls. - Initial call of EMRMD routine. - When surface crossing is detected. PRNLEV>=4 MRMD energy correction. PRNLEV>=5 MRMD energy correction for each surface. PRNLEV>=6 Summed energy correction of each type (NBON ELEC,NBON VDW,NBON GVDW,BOND HARM,BOND MORS, ANGL HARM, ANGL UREY,DIHE FOUR,IMPR HARM) for each surface. PRNLEV>=7 Each individual MRMD energy term, except for nonbonded interactions. Each GAPO energies. PRNLEV>=8 Each individual nonbonded energy term is also printed. ------------------------------------------- 2C. Printing forces at various detail level ------------------------------------------- PRNLEV>=9 All nonzero MRMD force corrections for each atom are printed. PRNLEV>=10 All forces of all potential energy surfaces are printed.
The mrmd_h2so4.inp test case with the mrmd_h2so4.par parameter file are distributed with the package. This example simulation describes the water elimination from a vibrationally highly excited sulfuric acid molecule. Two surfaces are defined: one for sulfuric acid molecule and one for water and sulfur-trioxid molecules. The initial conditions are prepared so that water elimination takes place within a few timesteps. All the force field parameters of the H2SO4 molecule are redefined in the parameter file (mrmd_h2so4.par) in order to demonstrate most of the features of the MRMD module.
Before running production calculations, it is highly recommended to check whether all the expected individual energy terms are cancelled out or/and added for each surfaces. Increasing the PRNLEV parameter up to 8, a detailed listings is done, showing also the values of all used parameters and geometric variables (ie. bond length) used for the calculation of the removed and added terms. If total energy seems to be wrong, this listing can help in identifing the energy terms which are responsible for it.
----------------------------------
Interfaces, subroutines, functions
----------------------------------
memory (de)(re)allocation interface
MRMD_ALLOC using subroutines:
MRMD_ALLOC_REAL1,MRMD_ALLOC_REAL2,MRMD_ALLOC_REAL3, MRMD_ALLOC_INTG1,MRMD_ALLOC_INTG2,
MRMD_ALLOC_INTG3,MRMD_ALLOC_LOGI1,MRMD_ALLOC_LOGI2,MRMD_ALLOC_CH16_1
reads and processes parameter file and sets up all global variables
MRMD_INIT
energies, forces for effective and individual surfaces
EMRMD,MRMD_ENERG,MRMD_SURF
switching related
MRMD_SWITCH_FUNCT,MRMD_SWITCH_WEIGHTS
energy terms:
MRMD_GAPO,MRMD_ELEC,MRMD_VDW,MRMD_GVDW,MRMD_HARM,MRMD_MORS,MRMD_RKHS,MRMD_ANGL,MRMD_DIHE
utility functions and subroutines for RKHS interpolation:
MRMD_RKHS_CALC_ASYM,MRMD_RKHS_DKERNEL,MRMD_RKHS_KERNEL,MRMD_CHOL_DECOMP,MRMD_CHOL_SOLVE
write out pdb file
MRMD_PDB
converting strings to uppercase:
function MRMD_UPPERCASE
complete deallocation
MRMD_DEALLOC
----------------------------------
Global variables in other routines
----------------------------------
variable for the preprocessor:
RMD: logical variable to include MRMD module or not into the code to be compiled
Fortran variables, constants:
logical MRMD_ACTIVE = whether MRMD module is active (reactive surface loaded)
integer MRMD = 103
string CETERM(MRMD) = 'ERMD'
=> use ?ERMD to request ERMD energy correction
logical QETERM(MRMD) = whether MRMD modul is compiled
=> use ?MRMD to request whether MRMD module is compiled
?MRMD .eq. 1 => yes
?MRMD .ne. 1 => no
energy ETERM(MRMD) = energy correction to PSF
--------------------------------------
Other modules refering to MRMD module:
--------------------------------------
charmm/miscom.src : calls MRMD_INIT
energy/energy.src : calls EMRMD
energy/energym.src: defines MRMD=103
energy/eutil.src : defines CETERM(MRMD)='ERMD'
misc/genetic.src : calls EMRMD
pert/epert.src : calls EMRMD
pert/icpert.src : calls EMRMD
CHARMM Documentation / Rick_Venable@nih.gov