$CIMFRG group
(required when CIMTYP=GSECIM in $CONTRL;
relevant to single- and multi-level GSECIM calculations)
In free-format, provide an array of the comma or dash
separated integers labeling the atoms. Integers in each
row, which must be listed in parentheses, represent the
user-specified group of atoms that are used to identify
central occupied LMOs in the initial steps of the GSECIM
subsystem design. When generating subsystems, the atoms
specified in a given row will always be kept together.
Optionally, one can provide the name of the method to be
applied to the group of atoms listed in a given row. The
hydrogen atoms do not have to be included, since each
hydrogen atom is assigned to the nearest non-hydrogen atom,
although there may be cases where listing hydrogens is
useful if there is ambiguity in assigning hydrogens to non-
hydrogen atoms. The subsystems involving atoms that are not
represented in $CIMFRG will be determined using the
standard SECIM algorithm and treated with the method
defined by SUBMTD in $CIMINP. For example,
$CIMFRG
1 (3-4,6) CR-CCL
2 (1-2,11) CCSD
3 (14,34,50)
4 (16,18,54)
5 (21-22,26,30)
6 (37-38,42,46)
$END
implies that atoms 3, 4, and 6 (and adjacent hydrogens, if
any) will be used to identify central LMOs that define the
1st subsystem in the initial steps of the GSECIM algorithm,
which will be treated by CR-CC(2,3). Similarly, atoms 1, 2,
and 11 (and adjacent hydrogens, if any) will be used to
identify central LMOs that define the 2nd subsystem in the
initial steps of the GSECIM algorithm, which will be
treated by CCSD, whereas atoms 14, 34, and 50 (and adjacent
hydrogens, if any) will be used to identify central LMOs
that define the 3rd subsystem in the initial steps of the
GSECIM algorithm which will be treated by the method
defined by SUBMTD in $CIMINP, etc. The remaining atoms will
be treated by the standard SECIM algorithm and the method
defined by SUBMTD in $CIMINP.
$QUANPO group (optional, relevant if QuanPol is used)
QuanPol (quantum chemistry polarizable force field
program) can perform MM, QM/MM or QM/MM/Continuum solvent
MD simulation, geometry optimization, saddle point search
and Hessian vibration frequency calculation, using HF, DFT,
DFTB, GVB, MCSCF, MP2, TDDFT or TDDFTB wavefunctions.
Semi-empirical metods (AM1, PM3) are implemented using LJ
and charge interactions between QM and MM atoms.
"QuanPol: A full spectrum and seamless QM/MM program",
Journal of Computational Chemistry, 2013, 34, 2816-2833.
Either $FFDATA (explicit input of coordinates and force
field parameters) or $FFPDB (Protein Data Bank coordinates
and built-in force fields) needs to be present, to define
the MM atoms. Quantum atoms, if any, are given in $DATA as
usual. If no quantum atom is used, use COORD=FRAGONLY in
$CONTRL. For free energy perturbation calculations,
$FFDATA and $FFDATB are required to define the reference
and target MM states, and $QUANPO keywords MATOMA and
MATOMB are required to define the reference and target QM
states.
Force field data sets are located by the environment
variable QUANPOL. Some parameter and topology files from
the CHARMM, AMBER and MMFF94 programs are included in
QUANPOL and can be read by QuanPol.
When some of the $QUANPO inputs become lengthy, use
multiple lines and '>' at the end of each line to glue them
together.
**** set up MM force field ****
MXFFAT = maximum number of MM atoms.
MXBOND = maximum number of MM bonds.
MXANGL = maximum number of MM bond angles.
MXDIHR = maximum number of MM dihedral rotation angles.
MXDIHB = maximum number of MM dihedral bending angles.
(i.e. improper torsion in CHARMM).
MXWAGG = maximum number of MM wagging angles.
MXCMAP = maximum number of CHARMM correction map cases.
All of the above are for memory allocation purposes.
Defaults are automatically determined and are almost always
good.
NFFTYP = select the force field type (no default)
= 0 user defined force field
= 20000-29999 CHARMM
= 30000-39999 AMBER (including GAFF)
= 40000-49999 OPLSAA
= 50000-59999 MMFF94
User defined force field can be input from $FFDATA, and/or
by using IADDWAT and supplying water potential in the
path/gamess/auxdata/QUANPOL/WATERIONS.DAT file. For
NFFTYP=0 the default WT14CH and WT14LJ are both 1.0. For
any molecule, it is a good idea to use LOUT=1 and NFFTYP=0
to generate a template $FFDATA, and then modify it.
QuanPol can also use combined force field, such as
AMBER/MMFF94. See ICOMBIN options.
CHARMM force field can be established for amino acids,
nucleic acids and simple ions by using the top/par files
from CHARMM developers. It is not made available for
general molecules. WT14CH=1.0, WT14LJ=1.0 (but uses a
second set of LJ potential). Note CHARMM typically
requires the use of ISHIFT=4, ISWITCH=2, SWRA=10, SWRB=12.
These must be input in $QUANPO.
* To use CHARMM36, give the following in $QUANPO:
NFFTYP=20000 NFFFILE=2
TOPFILE='path/gamess/auxdata/QUANPOL/top_all36_prot.rtf'
PARFILE='path/gamess/auxdata/QUANPOL/par_all36_prot.prm'
* To use CHARMM27, give the following in $QUANPO:
NFFTYP=20000 NFFFILE=2
TOPFILE='path/top_all27_prot_na.rtf'
PARFILE='path/par_all27_prot_na.prm'
AMBER force field can be established for amino acids,
nucleic acids and simple ions by using the top/par files
from AMBER developers. It is not made available for
general molecules. However, QuanPol LOUT=1 is able to
read AMBER GAFF files in the mol2 format (generated by
AmberTools), and read the AMBER gaff.dat parameter file
to establish the force field. WT14CH=1/1.2, WT14LJ=0.5
(if a second set of LJ potential is used, WT14LJ is set
to be 1.0). AMBER typically requires the use of
ISHIFT=0, ISWITCH=0, together with SWRB=12 or a larger
value. These must be input in $QUANPO.
* To use AMBER14 nonpolarizable protein force field:
NFFTYP=30000 NFFFILE=3 IDOPOL=0
TOPAMIA='path/gamess/auxdata/QUANPOL/amino12.in'
TOPNTER='path/gamess/auxdata/QUANPOL/aminont12.in'
TOPCTER='path/gamess/auxdata/QUANPOL/aminoct12.in'
PARFILE='path/gamess/auxdata/QUANPOL/parm10.dat'
PARFIL2='path/gamess/auxdata/QUANPOL/frcmod.ff14SB'
* To use AMBER14 ipq nonpolarizable protein force field:
NFFTYP=30000 NFFFILE=3 IDOPOL=0
TOPAMIA='path/gamess/auxdata/QUANPOL/amino12.in'
TOPNTER='path/gamess/auxdata/QUANPOL/aminont12.in'
TOPCTER='path/gamess/auxdata/QUANPOL/aminoct12.in'
PARFILE='path/gamess/auxdata/QUANPOL/parm14ipq.dat'
* To use AMBER12 polarizable protein force field:
NFFTYP=30000 NFFFILE=3 IDOPOL=100
TOPAMIA='path/all_amino12pol*.in'
TOPNTER='path/all_aminont12pol*.in'
TOPCTER='path/all_aminoct12pol*.in'
PARFIL2='path/frcmod.ff12pol*'
PARFILE='path/gamess/auxdata/QUANPOL/parm99.dat'
* To use AMBER12 nonpolarizable protein force field:
NFFTYP=30000 NFFFILE=3 IDOPOL=0
TOPAMIA='path/gamess/auxdata/QUANPOL/amino12.in'
TOPNTER='path/gamess/auxdata/QUANPOL/aminont12.in'
TOPCTER='path/gamess/auxdata/QUANPOL/aminoct12.in'
PARFILE='path/gamess/auxdata/QUANPOL/parm10.dat'
PARFIL2='path/gamess/auxdata/QUANPOL/frcmod.ff12SB'
* To use AMBER10 nonpolarizable protein/na force field:
NFFTYP=30000 NFFFILE=3 IDOPOL=0
TOPAMIA='path/gamess/auxdata/QUANPOL/amino10.in'
TOPNTER='path/gamess/auxdata/QUANPOL/aminont10.in'
TOPCTER='path/gamess/auxdata/QUANPOL/aminoct10.in'
TOPNUCA='path/gamess/auxdata/QUANPOL/nucleic10.in'
PARFILE='path/gamess/auxdata/QUANPOL/parm10.dat'
* To use AMBER02 polarizable protein/na force field:
NFFTYP=30000 NFFFILE=3 IDOPOL=100
TOPAMIA='path/gamess/auxdata/QUANPOL/all_amino02.in'
TOPNTER='path/gamess/auxdata/QUANPOL/all_aminont02.in'
TOPCTER='path/gamess/auxdata/QUANPOL/all_aminoct02.in'
TOPNUCA='path/gamess/auxdata/QUANPOL/all_nuc02.in'
PARFILE='path/gamess/auxdata/QUANPOL/parm99.dat'
PARFIL2='path/gamess/auxdata/QUANPOL/frcmod.ff02pol.r1'
* To use AMBER94 nonpolarizable protein/na force field:
NFFTYP=30000 NFFFILE=3 IDOPOL=0
TOPAMIA='path/gamess/auxdata/QUANPOL/all_amino94.in'
TOPNTER='path/gamess/auxdata/QUANPOL/all_aminont94.in'
TOPCTER='path/gamess/auxdata/QUANPOL/all_aminoct94.in'
TOPNUCA='path/gamess/auxdata/QUANPOL/all_nuc94.in'
PARFILE='path/gamess/auxdata/QUANPOL/parm94.dat'
* To use AMBER94 nonpolarizable protein/na force field:
NFFTYP=30000 NFFFILE=2 IDOPOL=0
TOPFILE=
'path/gamess/auxdata/QUANPOL/top_amber_cornell.inp'
PARFILE=
'path/gamess/auxdata/QUANPOL/par_amber_cornell.inp'
OPLSAA force field can be established for amino acids
and simple ions by using the top/par files from the
CHARMM package. It is not made available for general
molecules. WT14CH=0.5, WT14LJ=0.5 (if a second set of
LJ potential is used, WT14LJ is set to be 1.0). QuanPol
does not have the original OPLS switching functions for
charge-charge potential. The original OPLS uses a
switching function in 10.5-11.0 A (but in some cases
12.5-13.0 and 14.5-15.0 A) for both charge-charge and
LJ potentials.
* To use OPLSAA-96 nonpolarizable protein force field:
NFFTYP=40000 NFFFILE=2
TOPFILE='path/gamess/auxdata/QUANPOL/top_opls_aa.inp'
PARFILE='path/gamess/auxdata/QUANPOL/par_opls_aa.inp'
MMFF94 is implemented in QuanPol for general organic
molecules and some metal ions as described in the
original MMFF94 papers, and tested with the validation
suit (http://server.ccl.net/cca/data/MMFF94/) of 761
tests. All 761 tests can pass, with the largest positive
difference of +0.011 kcal/mol, the largest negative
difference of -0.007 kcal/mol, and a root mean square
difference of 0.002 kcal/mol. Therefore, this is a
complete implementation of MMFF94. It also works for
proteins and DNA/RNA molecules. When using LOUT=1 and
NFFTYP=50000 to generate MMFF94 force field for a
molecule in $FFDATA or $FFPDB, the parameter file
'MMFF-I_AppendixB.ascii' must be used. This file can be
downloaded from a JCC ftp server:
'ftp.wiley.com/public/journals/jcc/suppmat/17/490/'.
The keyword MMFF94Q is required for some cases.
The charge-charge interaction has a buffer distance of
0.05 A in MMFF94. However, in QuanPol implementation
when ISHIFT>0 or IEWALD>0 is used, this buffer
distance is not used (i.e. set to be zero).
There are a dielectric constant D and an index n in the
MMFF94 formula. Only D=1.0 and n=1 are implemented in
QuanPol. WT14CH=0.75, WT14LJ=1.0 (MMFF94 uses a 14-7
potential). It is not clear what shifting function,
switching function and cutoff distances should be used
in MMFF94 bulk MD simulations.
* To use MMFF94 nonpolarizable force field:
LOUT=1 NFFTYP=50000
PARFILE='path/MMFF-I_AppendixB.ascii'
MMFF94s is a variant of MMFF94. To use MMFF94s, one
needs to download two parameter files from:
ftp://ftp.wiley.com/public/journals/jcc/suppmat/20/720
These two files are named 'Table1.txt' and 'Table2.txt'.
One should replace the corresponding sections in
'MMFF-I_AppendixB.ascii' with 'Table1.txt' and
'Table2.txt', and save it as a new file such as
'MMFF94s.par'. QuanPol can reproduce MMFF94s results
available in the validation suit of 265 tests
(http://server.ccl.net/cca/data/MMFF94s/index.shtml).
* To use MMFF94s nonpolarizable force field:
LOUT=1 NFFTYP=50000
PARFILE='path/MMFF94s.par'
LOUT = create a $FFDATA group containing force field
parameters for a molecule. Require inputting
a $FFDATA group containing only COORDINATES.
Other sections in $FFDATA are not used.
= 0 no action (default)
= 1 create the $FFDATA group for a molecule
using force fields defined by NFFTYP.
For NFFTYP=0:
Can use JRATTLE to define coordination bonds that
are usually not considered as covalent bonds.
The equilibrium bond lengths and angles are
set to be those in the input geometry. Trivial
force constants and potential parameters are
assigned to covalent and noncovalent terms.
Users are supposed to know and edit the output
$FFDATA to assign formal charges to ions and
ionized groups. This option is most useful for
preparing a simplified force field for the QM
molecule to be used in QuanPol QM/MM calculation.
For NFFTYP=30000:
Need
PARFILE='path/gamess/
auxdata/QUANPOL/gaff.dat'
TOPFILE='path/xxx.mol2'
The xxx.mol2 file is generated by AmberTools.
For NFFTYP=50000:
Need PARFILE='path/MMFF-I_AppendixB.ascii'.
Formal charges on some atoms/ions must be
specified via keyword MMFF94Q.
= 910 to calculate the averaged projection area
of the molecule in $FFDATA, with relevance
to ion mobility spectrometry (IMS). This
also needs RALLMM to define the probing
gas (e.g. He, N2) radius in angstrom, and
NRADMM to define the atomic radii
(angstrom) in $FFDATA. Here for LOUT=910,
NRADMM is to define the atomic radii for
different elements. For example,
NRADMM=4 1 0.40 6 1.641 7 1.585 8 1.496
is to define H, C, N, O radii.
The average projection area is obtained
from 960 directions from the COM of the
$FFDATA molecule.
MMFF94Q= n, I1, Q1, I2, Q2, ... In, Qn
= specify the formal charge for up to 50 atoms
when LOUT=1 and NFFTYP=50000 are selected.
Note the default formal charge is zero, and
does not need user specification. QuanPol
is able to determine almost all formal charges
for H, C, N, O, F, Si, P, S, Cl, Br, I atoms
and simple ions. Therefore, MMFF94Q
is required only for multivalent metal
ions (e.g. Fe, Cu) and some special molecules.
MMFF94 atom types depend on formal charges.
n = number of atoms (default=0)
In = atom sequential number in $FFDATA
Qn = formal charge (e.g. -0.50, +2.0, +3.0)
For example, MMFF94Q=1 30 +3.0 is to assign
the 30th atom with +3.0 e charge.
**** set up QM/MM ****
QuanPol automatically performs QM/MM calculation when
$FFDATA (or $FFPDB) and $DATA are both detected in the
input deck and the numbers of QM atoms and MM atoms are
both greater than zero. Any MM atom that has virtually
the same Cartesian coordinates of a QM atom will be
identified and labeled, and vice versa. QM atoms will
be enforced to have the coordinates and masses of their
matching MM atoms.
QuanPol can also automatically generate a full set of
QM atoms (i.e., the $DATA group) from the $FFDATA based
on an input $DATA that contains one and only one atom
(e.g., a Zn2+ ion at the catalytic site). This function
is invoked by using IDIMER=3 with NMOLE=-1. The atoms
within 4.0 angstrom to the atom in $DATA are selected
to be QM atoms, as well as those forming covalent bonds
to these atoms. QM/MM boundaries are automatically
determined, and QM/MM link capping H atoms are set.
In addition, the generated QM atoms are automatically
grouped into fragments (each with roughly 5~25 atoms) to
facilitate the IDIMER=3 calculation. The fragmentation
at covalent bonds is also automatically determined for
large pieces of molecules in the QM region.
If there is no covalent bond between the QM and MM
regions, QM atoms will only feel the noncovalent
interaction potentials, such as charge, induced dipole
and LJ (or QMMMREP potential if LJQMMM=0), of the MM
atoms. There will be no covalent interactions between
QM and MM atoms. In general, there are no covalent and
noncovalent interactions between QM atoms, with one
exception for IFEPTYP=1 with MATOMA > MATOMB. In this
case, the bond, angle, dihedral rotation/bending, and
wagging terms involving any QM dummy atoms (i.e. those
appear in QM state A but not in QM state B) are retained
and scaled by WSIMUL. When WSIMUL=1.0, the pure state B
has a full strength of these MM covalent terms to ensure
that the QM dummy atoms stay in their positions. Using
WSIMUL is good because (1-WSIMUL)*QM + WSIMUL*MM is just
right.
If there are covalent bonds between QM atoms and MM atoms
all covalent force field interactions will remain in full
strengths if they involve at least one MM atom. If the
QM atoms are capping QM H atoms, the bond, angle and
dihedral rotation (only these three) terms involving the
QM H atoms and other QM atoms (but no MM atoms) will be
retained and scaled by RETAIN, which is defaulted as 0.5.
This is to compensate for the weakening of the covalent
terms due to the elongation of the capping H atom bond
length. If the QM atoms are not capping QM H atoms with
elongated bond lengths, it is not necessary to compensate
for the covalent terms, and RETAIN should be 0. All
covalent terms involving only QM atoms are excluded,
but may be retained for IFEPTYP=1 with MATOMA > MATOMB
(see the above paragraph).
To prepare an input file for QM/MM with covalent bonds:
a. Prepare a good input for pure MM calculation.
b. Copy some MM atoms from $FFDATA to $DATA to be QM
atoms. It is not necessary to delete these atoms
from $FFDATA. $DATA can have atoms not in $FFDATA
if semi-empirical (AM1, PM3) methods are not used.
The order of the atoms in $DATA does not matter for
their mapping to atoms in $FFDATA is solely based
on coordinates (but see MATOMA and MATOMB).
c. If necessary, change one or several of the QM atoms
in $DATA to be capping H atoms. The simplest way
is to only change the nuclear charge to 1.0 because
GAMESS recognizes atoms using their nuclear charges
rather than their names. The atom names can also be
changed to enhance readibility. For example, for
an alpha carbon in a protein, the following can
be seen in $DATA:
CX 1.0 X Y Z
The Cartesian coordinates X, Y, Z cannot be changed.
For DFTB, however, the correct chemical symbol must
be used in $DATA:
H 1.0 X Y Z
QuanPol will automatically match QM and MM atoms and
a. Zero off all covalent force field terms that involve
only QM atoms. Scale the covalent force field terms
(only bond, angle and dihedral rotation) that
involve QM and capping QM H atoms by RETAIN, which
is typically 0.5.
b. Zero off force field charges and polarizabilities
for all QM atoms. However, charges are saved for
QM atoms when semi-empirical methods (AM1, PM3)
are used.
c. Zero off force field LJ potentials for all QM atoms
if LJQMMM=0 is used. Retain LJ terms for QM atoms
if LJQMMM=1 is used, but exclude any LJ terms
between QM atoms. Semi-empirical methods must use
LJQMMM=1.
d. Zero off QMMMREP terms for all QM atoms, and also
for MM atoms that form covalent bonds to QM atoms.
e. Apply special QMREP to capping QM H atoms for
standard QM methods, or adjust the corresponding
force field equilibrium bond lengths for capping QM
H atoms for semi-empirical AM1 and PM3 methods.
For DFTB, the corresponding H, S and repulsion
values in the Slater-Koster tables will be
automatically adjusted for capping QM H atoms to
reproduce the correct bond lengths and electronic
structure.
RETAIN = retaining factor (0.0 - 1.0, default=0.5) for
force field covalent terms that involve only QM
atoms, one of which is a capping H atom with
a repulsion potential. The purpose is to
strengthen the weakened QM covalent terms
involving the capping H atom.
QMREP = NQMREP, IATOM1, NTERM1, C11, Z11, C12, Z12 ...,
IATOM2, NTERM2, C21, Z21, C22, Z22 ...
= specify effective Gaussian repulsion potentials
at capping QM atoms (typically H atoms in
place of C and N atoms of a peptide) to produce
the desired longer bond lengths.
NQMREP = number of QM atoms with Gaussian
potentials. Up to 200 atoms.
IATOMn = QM atom sequential number in $DATA.
When MATOMB > 0, IATOMn is only for
state A.
NTERMn = number of Gaussians at IATOMn, Up to 4
C11,.. = strength part of the Gaussian, e/bohr
Z11,.. = radial part of the Gaussian, 1/bohr**2
Must enter NTERMn pairs of C and Z for
atom IATOMn.
For example, to define 1 Gaussian for QM atom 1
and 3 Gaussians for QM atom 7, give
QMREP=2,1,1,3.0,3.0,7,3,8.0,6.0,3.0,3.0,0.3,1.0
For H atom forming C-H bond, a single Gaussian
with C=3.0 e/bohr and Z=3.0 bohr**(-2) is a good
option because it can create an equilibrium bond
length of ~1.50 angstrom for a C-C bond in
proteins.
** QMREP default **
QuanPol always automatically applies a single
Gaussian potentials (C=3.0, Z=3.0) to all QM
capping H atoms. Explicit input of QMREP can
override the QuanPol default values. If no
QMREP is wanted for a capping H atom, simply
input C=0.0 and Z=0.0 for the capping H atom.
LJQMMM = specify how the QM-MM repulsion and dispersion
are handled in a QM/MM system.
= 0 use QMMMREP (Gaussian potentials on MM).
Not available for semi-empirical methods.
= 1 use MM Lennard-Jones terms (default).
Only choice for semi-empirical methods.
LJQM = use LJ dispersion correction between QM atoms
when DFTB, HF, DFT methods are used. LJQM will
be set to 0 for MP2 methods. When LJQM=1 is
selected, other dispersion corrections for DFTB
or DFT methods should be turned off to avoid
over-corrections.
= 0 no use (default)
= 1 use LJ dispersion correction in a way that
only the attraction region (R>Rm in 12-6 LJ)
is used, the repulsion region is replaced by
a constant (negative epsilon in 12-6 LJ).
For MMFF94 style 'buffered 14-7 LJ', the
minimum energy is actually -1.0005648472872
times of epsilon at the distance Rm =
0.996182687569541*R_ij for a pair of atoms.
LJ parameters are from the input $FFDATA.
RDAMP = specify the effective distance (A) in the damping
function for polarizability. Default=3.0 A.
QuanPol scales MM polarizability using a Gaussian
function of interatomic distances R (in bohr):
S(R) = EXP[-0.0863*(R-RDAMP)**2]
This method works only for QM/MM. R are QM-MM
distances. A MM polarizability point is scaled
by all QM atoms within RDAMP to the MM point.
DFTBMM = specify the damping distance (A) in the damping
function for DFTB/MM charge-charge interaction.
Default=0.7 A for DFTB3. Reasonable values are
0.6 - 0.8 A. Values smaller than 0.5 A should
not be used because the QM-MM coupling is too
strong and the DFTB SCF tends to diverge.
QuanPol uses the following distance function
for charge-charge interaction in DFTB/MM:
S(R) = 1/SQRT(R**2+DFTBMM**2)
R is the QM-MM internuclear distances. This
applies to both DFTB nuclear charge and electron
charge.
For other QM/MM methods, DFTBMM is set to 0.0.
INTCHG = specify how noninteger charge is treated in
QM/MM. When covalent bonds are cut, the MM
region is often left with a noninteger charge.
= 0 no action (not recommended).
= 1 add missing charge to the MM atoms that form
covalent bonds to QM atoms.
For multiple such MM atoms, the missing
charge is evenly distributed. (default)
**** set up optimization and MD simulation ****
Use $CONTRL RUNTYP=OPTIMIZE, SADPOINT, HESSIAN and MD to
run QuanPol geometry optimization, saddle point or
transition state search, vibration frequency calculation
and MD simulation. For optimization jobs, most of the
$STATPT options can be used, but not all. $STATPT
PROJCT=.T. makes little sense to QM/MM jobs so it is
always false. The Hessian from pure MM calculation can
be used with MMHESS=1 and IHESS=1 to guide QM/MM geometry
search (RUNTYP=OPTIMIZE).
QM/MM saddle point search can be run in three steps. In
the first step the Hessian from a prior QM/MM reactant
optimization (RUNTYP=OPTIMIZE) can be used together with
JUMBUP=1 and JUMBPOT to run an optimization (RUNTYP=
OPTIMIZE) to obtain an approximate saddle point structure
and a BFGS updated positive definite Hessian. In the
second step, RUNTYP=HESSIAN, the Hessian from the first
step is read in and partially updated only for elements
associated with 3~6 selected QM atoms defined in LACTQM
(LACTMM must also be given to define MM atoms for the
Hessian). The QM atoms constrained by JUMBPOT in the
first step must be included in LACTQM, plus a few more
QM atoms. The partially updated Hessian should have only
one negative mode corresponding to the TS. In the third
step, RUNTYP=SADPOINT and the Hessian from the second
step is used (without JUMBPOT) to search for the saddle
point.
For Hessian jobs, the $FORCE options are not used.
(1). For pure MM systems, all MM atoms (if less than 2000
and LACTMM=0) or active site MM atoms (when LACTMM > 0)
are included in Hessian calculation. Pure MM Hessian
calculations are not restartable because they are cheap.
(2). For QM/MM systems without a LACTMM list, only QM
atoms are used to construct the force constant matrix,
but MM effects are included. When LACTQM > 0 is used,
only the active site QM atoms are displaced to calculate
force constants, nonactive site QM atoms are not displaced
(time saving!) and some trivial negative force constants
are used to produce imaginary vibration frequencies of
1.000*i cm-1. To get isotope effects, an already computed
$HESS group can be supplied in the input file, together
with a modified atomic mass in the PARAMETER section of
$FFDATA ($MASS will not work). Of course, Hessian jobs
do not work with RATTLE at all. QuanPol QM/MM Hessian
jobs are restartable by using the saved $VIB group in the
new input (the LACTQM list of QM atoms can be different
in a restart job).
(3). For QM/MM systems with a LACTMM list and an input
$HESS group, the input Hessian is partially updated only
for the elements associated with the QM atoms in LACTQM.
The updated Hessian is saved in the .HS1 file. This is
typically done to obtain a good initial Hessian for TS
search. At the end frequencies will be generated from
the QM part of the QM/MM Hessian for inspection. No
RATTLE can be used. By default single displacment will
be used with a size of 0.01 bohr, but double displacement
will be used with a size of 0.01 bohr if $QUANPO NSTEP=2
is specified. If $QUANPO NSTEP=0 is given, the QM/MM
$HESS will be read in and the frequencies of the QM part
will be generated. Restartable by using the saved $VIB
group in the new input (the LACTQM list of QM atoms can
be different in a restart job, and NSTEP=2 can be added
to request double displacement).
DT = MD time step size. Default=1.0D-15 second is
usually good, especially when RATTLE is used.
QuanPol sets a maximum value of 0.2 bohr/DT
(i.e., 10583.5 m/s for DT=1 fs) for each
velocity component when ITSTAT is applied.
Such a limit can prevent chaotic behavior,
which may occur frequently for large systems.
NSTEP = number of MD or OPTIMIZE steps (default=1000).
This NSTEP overrides the $STATPT NSTEP.
MDOPT = controls global geometry optimization on an MD
trajectory. At every MDOPT (e.g., 1000) steps,
the MD is paused for a geometry optimization
(maximum 20000 steepest descent steps). The
optimized coordinates are printed out in the log
file. The MD can be run at any TEMP0, like 300
K or 900 K. This method can be used to find the
gas phase global minimum of peptides.
Default=0, meaning do not do it.
ISWAP = n, R
= swapping QM and MM water molecules in QM/MM
MD simulation on-the-fly at step 0 and every n
steps.
A QM water molecule farthest from the QM
center will be identified, and swapped with
a MM water molecule closest to the QM center
(and closer than the farthest QM one).
Only one QM water will be swapped with one
MM water at every n steps, and only those
MM water molecules within R angstrom
to the QM center will be considered. This
works only for 3-point water models, with
the purpose of keeping QM water molecules
in the QM region. Energy and forces are
discontinuous upon the swapping.
The default is ISWAP=0, 12.0, meaning no
swapping. For example, ISWAP=1000, 12.0 will
try to swap at step 0 and every 1000 MD
steps for water molecules within 12.0
angstrom to the QM center. It is found that
for DT=1.0D-15 s, ISWAP=1000 works very well.
ISWAP<500 may be unnecessarily too frequent,
while ISWAP>2000 may be too infrequent.
MMHESS = specify if the supplied $HESS group for QM/MM
geometry optimization is generated from pure
MM Hessian calculation. If yes, the force
constant matrix will be re-ordered to match
the atom sequence in the QM/MM job, in which
the Hessian always has QM atoms placed ahead
of MM atoms. Default = 0.
= 0 no, the $HESS does not need re-ordering.
This means that the $HESS is from a
previous QM/MM calculation.
= 1 yes, the $HESS needs re-ordering to match
QM/MM atoms. This means that the $HESS
is from a pure MM Hessian calculation.
LACTMM = n, K1, R1, ..., Kn, Rn (case 1, when n =< 2)
= define sphere radii in angstrom for n MM atoms
in $FFDATA. All MM and QM atoms included in
these spheres are defined as active site atoms.
At most 2000 QM and MM atoms can be active site
atoms so R should be typically < 16.5 A.
Default n=0 means that all MM atoms are active
site atoms but note that LACTQM may also invoke
some active site MM atoms. Apparently, atomic
position changes will affect the number of
active site atoms. To be consistent, use the
automatically generated LACTMM for restart jobs.
For example, LACTMM=2 100 16.5 106 15.0 is to
define all QM and MM atoms within 16.5 A to
the 100th MM atom and those within 15.0 A to
the 106th MM atom as active site atoms.
To explicitly define n (n=<2) MM atoms as
active site atoms, one can give a series of
small radii such as 0.1 A after each atom.
= n, I1, I2, ..., In (case 2, when n > 2)
= explicitly define n MM atoms in $FFDATA to be
active site atoms. This rather lengthy array
is always generated (and printed into the .dat
file) by using LACTMM with n =< 2 for restart
jobs. Default n=0 means that all MM atoms are
active site atoms. Maximum n=2000.
LACTQM = n, K1, R1, ..., Kn, Rn (case 1, when n =< 2)
= define sphere radii in angstrom for n QM atoms
in $DATA. All MM and QM atoms sitting in
these spheres are defined as active site atoms.
At most 2000 QM and MM atoms can be active site
atoms so R should be typically < 16.5 A.
Default n=0 means that all QM atoms are active
site atoms but note that LACTMM may also invoke
some active site QM atoms. Use the
automatically generated LACTQM for restart jobs.
For example, LACTQM=2, 5 15.0, 8 15.0 is to
define all QM and MM atoms within 15.0 A to
the 5th QM atom and those within 15.0 A to
the 8th QM atom as active site atoms.
To explicitly define n (n=<2) QM atoms as
active site atoms, one can give a series of
small radii such as 0.1 A after each atom.
= n, I1, I2, ..., In (case 2, when n > 2)
= explicitly define n QM atoms in $DATA to be
active site atoms. This rather lengthy array
should be generated by using LACTQM with
n =< 2 for restart jobs. Default n=0 means
that all QM atoms are active site atoms.
Maximum n=2000.
Active site QM and MM atoms defined using LACTMM and
LACTQM will move in geometry optimization, Hessian
calculation, and MD simulation. Nonactive site atoms
will be absolutely fixed in their input Cartesian
coordinates (even IRATLLE and JRATTLE cannot move them).
However, active site QM and MM atoms can still be
explicitly fixed by using NFIXMM and NFIXQM.
IHESS = request Hessian guided geometry optimization.
Works for all MM and QM/MM cases but is limited
to 2000 movable atoms. This keyword is
irrelevant to RUNTYP=HESSIAN.
Default = 1 for QM/MM systems. Otherwise 0.
= 0 no Hessian, use steepest descent. It can
converge, but very slow and may require
$QUANPO NSTEP=30000 or more. This is
recommended for pure MM systems because
Hessian methods are very time consuming.
In addition, for large pure MM systems, we
suggest the use of $STATPT OPTTOL=1.0D-05.
It is the default for QM=AM1, PM3, DFTB.
= 1 use Hessian. Hessian is usually
good for ~100 optimization steps. After
that, the Hessian may become wrong and
lose its guiding power. Use $STATPT to
input Hessian options. If $HESS is from
a previous pure MM run, MMHESS=1 must be
used. IHESS=1 is recommended for QM/MM
because steepest descent method converges
very slowly. Typically the default $STATPT
OPTTOL=1.0D-04 is good, but occasionary
for very flat potential energy surfaces,
such as H transfer from one O to another O,
1.0D-05 should be used to avoid false
identifications of minima. When IHESS=1
is selected, the $HESS at each optimization
step is alternatively printed out in the
.hs1 and .hs2 files.
JOUT = report simulation information such as energies
and temperature to the log file every JOUT steps.
Default=1.
KOUT = in RUNTYP=MD write coordinates and velocity to
the trj file every KOUT steps (default=100).
= in RUNTYP=OPTIMIZE coordinates are written to
the trj file every KOUT steps (default=1).
CENTX =
CENTY =
CENTZ = define the center of the PBC master box or the
sphere center. If not found, it is
automatically calculated. Use the same CENTX,
CENTY, CENTZ for restart jobs.
XBOX =
YBOX =
ZBOX = size of the periodic box in angstrom.
Default 1.0D+30 means no PBC is imposed.
EFIELDX=
EFIELDY=
EFIELDZ= external electric field in x, y, z directions
applied to MM charges. Defaults 0.0 V/m.
Note: 1 V/m = 1.9446905D-12 atomic unit.
IADDWAT= add water molecules to the system
= 0 no adding water (default)
= 1 add water in PBC master box
= 2 add water in a sphere
When MMFF94 is used, it is better to add water
with ITYPWAT=301 and obtain a good $FFDATA,
then use LOUT=1 and NFFTYP=50000 to apply the
MMFF94 force field to all atoms in $FFDATA.
ITYPWAT= select the water model. Rigid water models
should use IRATTLE=10 or 20 in MD. For 5-point
water models the following (real) masses are
implemented:
O=14.000, H=1.008, M=1.000 (amu)
Users can also define their own water models.
One way is to add new water models to the library
file path/gamess/auxdata/QUANPOL/WATERIONS.LIB,
the other is to directly modify the parameters in
$FFDATA and $FFDATB groups already generated by
QuanPol for one of the built-in models (use a
similar one). ITYPWAT=300 and 500 should be used
for user defined 3-point and 5-point models.
4-point and 6-point models are not supported.
= 0 no specification (default)
= 301 flexible nonpolarizable 3-point model
= 302 flexible polarizable 3-point model
= 303 flexible SPC/Fw model by Wu/Tepper/Voth,
J.Chem.Phys. 124, 024503 (2006).
= 304 rigid TIP3P (LJ terms for H atoms, CHARMM)
= 305 rigid TIP3P (no LJ term for H atoms, AMBER)
= 306 rigid SPC
= 307 rigid SPC/E (extended SPC)
= 308 rigid POL3, J.Phys.Chem.99,6208 (1995)
= 504 rigid TIP5P-E, J.Chem.Phys.120,6085 (2004)
= 505 rigid TIP5P, J.Chem.Phys.112,8910 (2000)
= 300 user defined 3-point model. See IRATTLE.
= 500 user defined 5-point model. See IRATTLE.
JADDNA1= add Na+ ions to DNA/RNA PO4- sites.
JADDK1 = add K+ ions to DNA/RNA PO4- sites.
= 0 skip (default)
= 1 add NA+/K+ ions to all possible PO4- sites
IADDNA1= number of Na+ ions randomly added. Default=0.
IADDK1 = number of K+ ions randomly added. Default=0.
IADDCA2= number of Ca2+ ions randomly added. Default=0.
IADDMG2= number of Mg2+ ions randomly added. Default=0.
IADDCL1= number of Cl- ions randomly added. Default=0.
KMASTER= request relocating water molecules and ions to
the PBC master box (centered at CENTX, CENTY,
CENTZ) in MM and QM/MM MD simulation. The
default KMASTER is n=100. No effect on none
MD or none PBC jobs.
= 0 no relocation
= n do relocation every n MD steps. This
will not cause problem for NDFS or
IFEPTYP>0. This will affect NDIEL if the
PBC system has mobile ions, which will be
relocated and lead to large fluctuations
in dipole moments.
KOUTACT= n, R
n specifies the n-th atom in $FFDATA;
R is a radius, typically 10 A.
Active site atoms within R to the n-th atom are
printed out in log file for visualization.
Periodic boundary condition is considered.
default n=0 and R=0.0.
ITSTAT = enable the thermostat (velocity scaling)
= 0 no thermostat (default)
= 1 Berendsen. Scale all velocities at every MD
step so that
T' = (1-(DT/TT))*T + (DT/TT)*T0
Eq (11) in J.Chem.Phys. 81, 3684 (1984).
T = temperature
T0 = target or bath temperature TEMP0
TT = BERENDT (default 200 fs).
If TT>>DT, virtually no scaling.
If TT =DT, complete scaling to T0.
If T-T0 > 100 K, TT=10*DT is used.
If T-T0 > 200 K, TT= DT is used.
Berendsen thermostat tends to give an average
T slightly (0.01~0.3 K) lower than T0.
= 2 Andersen. Reassign Maxwell-Boltzmann
velocities to 20% randomly selected atoms
at every MD step. The velocity components
of all selected atoms are assigned as:
v = sigma*SQRT(-2*Ln(u1))*Cos(2*Pi*u2)
u1 = random number in (0,1)
u2 = random number in (0,1)
sigma = SQRT(k*T0/m)
k = Boltzmann constant
T0 = target or bath temperature TEMP0
m = mass of the atom
Andersen thermostat is not good for time-
dependent properties such as diffusion
coefficient and vibrational spectrum.
IPSTAT = enable the barostat (volume scaling)
= 0 no barostat (default)
= 1 Berendsen barostat at every MD step:
mu = [1 - (BETA*DT/TP)*(P0-P)]**(1/3)
Eq (21) in J.Chem.Phys. 81, 3684 (1984)
misses the BETA.
P = pressure, could be 1000 bar.
P0 = target or bath pressure PRES0
BETA = 4.9D-05 bar-1
TP = BERENDP (default 200 fs)
To enhance MD stability, mu larger than
1.0001 is set to be 1.0001, smaller than
0.9999 is set to be 0.9999.
= 3 Berendsen barostat at every MD step, but
separately in x, y, z directions. This
is necessary for anisotropic simulation
systems such as lipids in water. This
has little effect on isotropic systems.
A barostat is meaningful only for PBC system.
BEREND = BERENDT, BERENDP
= Berendsen thermostat coupling time for ITSTAT=1,
Berendsen barostat coupling time for IPSTAT=1,3.
Default = 200.0D-15 second (200 fs) for both.
Example: BEREND=200.0D-15, 200.0D-15.
TEMP0 = bath temperature in K. For MD jobs, there is no
default, and must be input by user. For other
jobs, the default=298.15 K.
This is the temperature for Hessian
thermochemistry calculation.
PRES0 = bath pressure in bar. Default=1.0 bar.
A pre-equilibrium system may show huge positive
or negative pressures like 100,000 bar.
An equilibrium system may show pressures
fluctuating by several tens or hundreds bar.
When IPSTAT=1 or 3 is used, the average pressure
converges slowly to PRES0 in ~100,000 fs.
INTALG = MD integrator algorithm.
= 1 Beeman algorithm (default)
This Beeman algorithm does not require a
step back. Instead, its first step is simply
a velocity verlet step.
= 2 velocity verlet algorithm
NRANDOM= selects the seed for QuanPol's random number
generator:
= 0 use fixed seeds (default)
= 1 use time/date to generate seeds
QuanPol uses a 16-bit pseudo-random integer
generator with a cycle length 6953607871644. See
Wichmann & Hill, Appl.Statist. 31, 188-190 (1982)
Fixed and time/date seeds should give the same
randomness.
IRATTLE= apply RATTLE to constrain bond lengths in MD.
It is irrelevant to geometry optimization or
Hessian vibration calculation.
Good for all MM and QM/MM systems, but
QM atoms will not be rattled unless IRATQM=1.
RATTLE does not affect atoms fixed by NFIXMM,
NFIXQM, or nonactive site
atoms defined by LACTMM.
IRATTLE=1 is recommended because it is fast.
Others are relatively slow. See SCALRAT.
QuanPol does not distinguish between 1-3
Urey-Bradley terms and regular 1-2 bond terms:
terms with force constants < 100 kcal/mol/A**2
are not constrained by RATTLE (except for zero-
strength bond between two H atoms).
QuanPol RATTLE recognizes H atoms by the
nuclear charge (NUC=1.0), which does not affect
any interaction potential. So, any point can
be given NUC=1.0 (and also a mass), and treated
like an H atom by RATTLE. Five-point water
models have 4 points given NUC=1.0 in order to
use IRATTLE options 10 and 20. It is not
recommended to use NUC=1.0 for other atoms
because some calculations, such as FIXSOL, also
rely on NUC to recognize atoms.
= 0 skip (default).
= 1 constrain bonds that involve H atoms and have
bond constants larger than 100 kcal/mol/A**2.
If water models are involved, this will
constrain the O-H and O-M bonds in 3-point
and 5-point models because they are usually
500 kcal/mol/A**2 strong. The H-H, H-M and
M-M distances will not be constrained if
their strengths are 0.0, but will be
constained if their strengths are > 100
kcal/mol/A**2.
= 10 constrain bonds that involve H atoms and have
bond constants larger than 100 kcal/mol/A**2,
and bonds that involve two H atoms and have
zero bond constants. This option is designed
for systems involving rigid 3-point and
5-point water models: it will constrain the
O-H and O-M bonds, as well as the H-H, H-M
and M-M distances.
= 2 constrain all bonds that have bond constants
larger than 100 kcal/mol/A**2. This option
can be used for any rigid molecules when
3N-6 (3N-5 for linear molecule, N=number of
mass points) independent bonds are defined.
Rigid solvent models must use RATTLE.
If water models are involved, this will
constrain the O-H and O-M bonds in 3-point
and 5-point models because they are usually
500 kcal/mol/A**2 strong. The H-H, H-M and
M-M distances will not be constrained if
their strengths are 0.0, but will be
constained if their strengths are > 100
kcal/mol/A**2.
= 20 constrain all bonds that have bond constants
larger than 100 kcal/mol/A**2, and bonds that
involve two H atoms and have zero bond
constants. This option is designed
for systems involving rigid 3-point and
5-point water models: it will constrain the
O-H and O-M bonds, as well as the H-H, H-M
and M-M distances.
JRATTLE= n, I1, J1, R1, ..., In, Jn, Rn
n specifies the number of atom pairs to be
constrained using the RATTLE scheme, or some
additional bonds (especially some coordinate
bonds that are significantly longer than normal
covalent bonds) to be used by LOUT=1 in
generating a force field (in this case, Rn must
be given but will not be used).
I1, J1, R1 are the sequential numbers and target
distance (A) of the 1st pair of atoms. Must give
n pairs. n can be 0 - 10.
When both IRATTLE and JRATTLE are used, JRATTLE
pairs (if new after check) are added to IRATTLE
pairs. JRATTLE does not affect atoms fixed
by NFIXMM, NFIXQM. Default=0.
IRATQM = specify the rattle of QM atoms in a QM/MM system
when IRATTLE and/or JRATTLE are used.
= 0 use no rattle for QM atoms, and no rattle
between QM atoms and MM atoms.
= 1 use rattle for QM atoms. QM atoms will be
rattled if and only if they appear as
rattled MM atoms in $FFDATA. This option
must be used (thus the default) for QM/MM
jobs when the QM part contains TIP5P style
water molecules.
The defaults should not be changed unless one
really knows the working mechanism of rattle
in QM/MM MD.
RATOLC =
RATOLV = convergence criteria in RATTLE step 1 for
coordinate and step 2 for velocity.
Default RATOLC=1.0D-05 and RATOLV=1.0D-08.
Loose criteria may destroy energy conservation
while tight criteria are costly.
MXRATT = maximum iterations in RATTLE steps 1 and 2.
Usually 4 iterations are enough for IRATTLE=1.
For rigid 5-point water models (IRATTLE=10 or 20)
it requires ~30 iterations when SCALRAT=1.5 is
used (default=200).
SCALRAT= scaling factor in RATTLE correction to coordinate
and velocity. Default is:
1.0 for IRATTLE=1
1.3 for IRATTLE=10 or 20 and 3-point water models
1.5 for IRATTLE=10 or 20 and 5-point water models
NFIXQM = specifies QM atoms in $DATA to be fixed in MD
simulation and geometry optimization. If any of
these fixed QM atoms appear in $FFDATA as MM
atoms, these MM atoms will also be fixed.
To fix 2 QM atoms, 5 and 18, give:
NFIXQM = 2, 5, 18
At most 200 QM atoms can be fixed. Default=0.
If this input is lengthy, use multiple lines
and '>' at the end of each line to glue them
together.
NFIXMM = specifies MM atoms in $FFDATA to be fixed in MD
simulation and geometry optimization. If any of
these fixed MM atoms appear in $DATA as QM atoms,
these QM atoms will also be fixed. To fix 4
MM atoms, 100, 1234, 9999 and 70012, give:
NFIXMM = 4, 100, 1234, 9999, 70012
At most 200 MM atoms can be fixed. Default=0.
NFIXQM and NFIXMM are absolutely enforced, even
IRATTLE and JRATTLE cannot affect them. They also
fix active site atoms defined using LACTMM and LACTQM.
NFXIQM and NFIXMM will be automatically set in potential
of mean force (PMF) free energy perturbation (FEP) MD
simulation (i.e., IFEPTYP=2) by the KFREEAB keyword.
SCFTYP2, TDDFT2, CITYP2, MPLEVL2, MULT2, ICHARG2
= similar and default to the keywords SCFTYP,
TDDFTYP, CITYP, MPLEVL, MULT and ICHARG in group
$CONTRL, but to define a different QM calculation
on the MD trajectory. For example, when DFT/MM MD
is performed, one can use TDDFT2=EXCITE to
request a TDDFT calculation and obtain vertical
excitation energies. Apply only to MD.
The results are printed out as '2ND' potential
energies in the log file.
Defaults are their counterparts in $CONTRL,
meaning no additional QM calculation.
**** MD free energy simulation ****
MATOMA =,MATOMB=,MCHARGA=,MCHARGB=,MULTA=,MULTB=
= specify the numbers of atoms, the charges, and
the multiplicities of QM state A and state B in
QM/MM free energy calculations.
MATOMB can be smaller but never be greater than
MATOMA. Defaults:
MATOMA = the number of atoms in $DATA
MCHARGA = ICHARG in $CONTRL
MULTA = MULT in $CONTRL
MATOMB = 0
MCHARGB = 0
MULTB = 1
The input in $DATA must have all the atoms for
state A defined before the atoms for state B.
If QM state A is listed in $FFDATA, QM state B
must also be listed in $FFDATB. For IFEPTYP=2,
the coordinates of MATOMB atoms may be
different from those of MATOMA.
(1) The sum of MATOMA and MATOMB must equal
to the total number of atoms in the $DATA.
(2) The sum of MCHARGA and MCHARGB must
equal to the total charge defined by
ICHARG in $CONTRL.
(3) MULTA and MULTB must be reasonable for
state A and state B, respectively.
IFEPTYP= specify the type of free energy perturbation
(FEP) calculation (i.e. what free energy or
free energy change to be calculated).
In FEP only one set of atoms are used to run MD
and sample the phase space, but the solute atoms
in the KFREEAB list can have different force
field parameters (IFEPTYP=1), or different
(fixed) coordinates (IFEPTYP=2), for states A
(in $FFDATA) and B (in $FFDATB). Both IFEPTYP=1
and 2 need KFREEAB to specify the perturbed
atoms in $FFDATA and $FFDATB. When QM atoms are
involved, MATOMA and MATOMB are also required to
define QM states A and B (both in $DATA). Note
it is almost meaningless to use different
settings in switching and shifting functions in
relative free energy calculation.
= 0 use no potential energy (default). This null
option gives zero free energy change, so no
FEP will be performed, and no need for
$FFDATB, KFREEAB, MATOMA, MATOMB...
= 1 For pure MM, use solvation potential energy
of the solute molecules (i.e. the KFREEAB
atoms). This is equivalent to
excluding the internal potential energy of
the solute molecules from the total energy.
For pure MM systems, this option is usually
called solvation free energy perturbation.
For QM/MM, the total QM/MM potential
energy (including solvation and QM internal
energy) is used. The coordinates of MATOMB
atoms (i.e., KFREEAB atoms in $FFDATB) must
be the same as those of MATOMA atoms
(i.e., KFREEAB atoms in $FFDATA) but may be
less in number (e.g., deprotonated) or
smaller in atomic numbers (e.g., F atoms
become H atoms) as to define an alchemical
perturbation.
= 2 Potential of mean force (PMF) FEP MD.
For pure MM, use solvation potential energy
of the solute molecules (i.e. the KFREEAB
atoms), the covalent potential energy
within the KFREEAB atoms, and noncovalent
potential energy within the KFREEAB atoms.
KFREEAB atoms are fixed in states A and B.
The only differences between $FFDATA and
$FFDATB are the coordinates of the KFREEAB
atoms.
For QM/MM, the total QM/MM potential
energy (including solvation and QM internal
energy) is used. The coordinates of MATOMB
atoms (i.e., KFREEAB atoms in $FFDATB) may
be different from those of MATOMA atoms
(i.e., KFREEAB atoms in $FFDATA) as to form
a geometric perturbation, but A and B must
be the same molecule with the same number
and types of atoms. In any case, MATOMB
cannot be greater than MATOMA.
KFREEAB= n, K1, K2, ..., Kn
n specifies the number of atoms in $FFDATA and
$FFDATB to be included in FEP calculation,
K1 - Kn are the sequencial number of the n atoms
in $FFDATA and $FFDATB. The limit of n is 200.
If this input is lengthy, use multiple lines
and '>' at the end of each line to glue them
together. Default is 0.
WSIMUL =, WPERT1=, WPERT2=
= three weights (i.e., coupling coefficient
lambda) of the target state B in IFEPTYP=1
calculation. All values must be from 0.0 to 1.0.
Defaults=0.0, 1.0, 1.0.
Usually a series of values are used to build a
perturbation route forth and back. The system is
simulated in the state defined with WSIMUL, and
the free energy differences are calculated for
the two states defined with WPERT1 and WPERT2.
The default value of WPERT2 = WPERT1. Examples:
WSIMUL=0.5 WPERT1=0.0 WPERT2=1.0
WSIMUL=0.5 WPERT1=0.6 WPERT2=0.7
WSIMUL=1.0 WPERT1=0.9 WPERT2=0.8
For IFEPTYP>0, $FFDATB must have the same topology as
$FFDATA: same number of atoms, same coordinates except
for atoms fixed by KFREEAB for IFEPTYP=2, same number
of bonds, angles, dihedrals and other covalent terms.
For IFEPTYP=1, the parameters (e.g., mass, charge, LJ
potential, bond constant, angle bending constant and
others) associated with alchemical atoms in the KFREEAB
list can be different in order to define two different
states. IFEPTYP=1 MD is performed on the potential energy
surface (PES) for that the covalent potential parameters
(e.g. bond lengths and constants) are combined using
F(W)=(1-W)*F(A)+W*F(B), charge and LJ potential energies
are combined using E(W)=(1-W)*E(A)+W*E(B).
IFEPTYP=2 MD is performed on the potential energy
surface (PES) of A (A equals to B for IFEPTYP=2), the
energy change is caused by the perturbation in the
coordinates from A to B (and only for atoms in KFREEAB).
The best way to create the $FFDATB for IFEPTYP>0 is to
modify the $FFDATA by changing the names, masses, bond
constants, charges and LJ parameters (or coordinates if
IFEPTYP=2) of the solute atoms in the KFREEAB list.
JUMBUP = 0 no action. (default)
= -1 adjust the R0 value on the fly for
JUMBPOT=12 when RUNTYP=OPTIMIZE so that the
energy is minimized. This is useful when
JUMBPOT=12 is used to locate a minimum point
on the potential energy surface. If the
adjustment and optimization are not
converging, there is unlikely a minimum
point in the given region of the potential
energy surface. The convergence criterion
of R0 is that the bias potential energy be
less than 1.6D-6 hartree.
= 1 adjust the R0 value on the fly for
JUMBPOT=12 when RUNTYP=OPTIMIZE so that the
energy is maximized. This is useful when
JUMBPOT=12 is used to locate a saddle point
on the potential energy surface. If the
adjustment and optimization are not
converging, there is unlikely a saddle
point in the given region of the potential
energy surface. Saddle points sometimes
are difficult to locate so a few trials with
the bias potential on different bonds may
be required. The convergence criterion
of R0 is that the bias potential energy be
less than 1.6D-6 hartree. The Hessian,
after a necessary modification, can be used
for RUNTYP=SADPOINT, which can find the true
saddle point.
JUMBPOT= NTYP, I1, I2, I3, I4, FC, R0
= apply umbrella sampling bias (harmonic) potential
to a reduced or combined MM internal coordinate:
V_bias = 0.5*FC*(R - R0)**2
1D histograms are printed out to the .log file
every JOUT steps, with 61 bins and bin size of
either 0.01 A or 1.0 degree.
JUM2POT= NTYP, I1, I2, I3, I4, FC, R0
= apply a second umbrella sampling bias potential
to a reduced or combined MM internal coordinate.
2D histograms are printed out to the .trj file
every KOUT steps, with 3721 bins and bin size of
either 0.01 A or 1.0 degree.
If selected, these bias potentials are added to
all MM and QM/MM MD/OPT calculations.
So, they can also be used for transition state
search. A transiton state can often be located
by using RUNTYP=OPTIMIZE and a single bias
potential JUMBPOT=12 with a FC value such
as 300 to 3000 kcal/mol/A**2. JUMBUP=1 can be
used to automatically adjust R0 values on-the-
fly to precisily determine the transition state
geometry. The FC value may heavily affect the
convergency of the R0 value and the optimization
process.
RUNTYP=HESSIAN/SADPOINT cannot have JUMBPOT or
JUM2POT.
The QuanPol Weighted Histogram Analysis (QPWHA)
program can be used to obtain 1D and 2D PMF
profiles. For NTYP=12 (e.g. for Na+ and Cl-
ions), the 1D PMF obtained from QPWHA program
must be corrected by a relative volume-entropy
term, which is kT*Ln((R/R0)**2). Here k is
Boltzmann constant, T is temperature, R is the
distance, and R0 is the reference distance at
which the PMF is set to be zero.
To obtain good 1D PMF, at least 100,000 MD steps
are required for each window (i.e. each R0). To
obtain good 2D PMF, at least 1000,000 MD steps
are required for each 2D window. Therefore, 2D
umbrella sampling is very expensive.
NTYP= define the internal coordinate R
= 0 nothing (default)
= 12 R = R12 (needs kT*Ln((R/R0)**2) )
= 1212 R = R12 - R'12
= 123 R = angle 123 (0-180 deg)
= 1234 R = dihedral angle 1234 (0-360 deg)
Ii = atoms in $FFDATA. Must give four integers,
but some or all can be 0.
FC = force constant, either in kcal/mol/A**2 or
kcal/mol/deg**2, depending on NTYP
R0 = equilibrium R, either in A or degree
Six examples for setting up 1D umbrella sampling:
JUMBPOT= 12 8 5 0 0 120.000 3.00
JUMBPOT= 1212 3 4 4 5 80.000 -0.20
JUMBPOT= 1212 3 5 7 8 100.000 1.50
JUMBPOT= 123 6 2 7 0 0.010 120.00
JUMBPOT= 1234 2 6 9 4 0.010 340.00
An example for setting up 2D umbrella sampling:
JUMBPOT= 12 8 5 0 0 120.000 2.20
JUM2POT= 12 10 25 0 0 120.000 1.50
An example for setting up 2D umbrella sampling:
JUMBPOT= 12 8 5 0 0 120.000 1.20
JUM2POT= 1234 10 11 19 20 0.010 120.00
IRMDF = I1, I2, R1, R2, N
= apply thermodynamic integration in MD simulation
to evaluate the mean force between two atoms,
the distance between which is constrained via a
RATTLE-like scheme. Can coexist with RATTLE, but
the IRMDF atoms will not be affected by RATTLE.
Works for MM and QM/MM MD simulations. The two
atoms can be MM or QM, both or either.
Default is 0, 0, 0, 0, 0.
I1 = MM atom in $FFDATA.
I2 = MM atom in $FFDATA.
R1 = starting distance between I1 and I2 in A,
must be between 0-100 A.
R2 = ending distance between I1 and I2 in A,
must be between 0-100 A. R2 can be larger
or smaller than R1.
N = the number of evenly distributed distances
in between R1 and R2 for that the MD
simulation will be consecutively run,
NSTEP/N steps for each distance.
N must be an integer between 1-100. If N=1,
the simulation will be run for (R1+R2)/2.
For example, inputing
IRMDF= 98, 100, 2.0, 3.0, 10
will evaluate the mean force between MM atoms
98 and 100 at 10 distances from 2.05 to 2.95 A:
2.05, 2.15, 2.25, ..., 2.85, 2.95
If NSTEP=1000000, for each distance 100000 MD
steps will be run to obtain the mean force.
The mean force will be multiplied by 0.10 A,
which is (R2-R1)/N, to produce the free
energy change for the 10 distance bins:
delta G from 2.00 to 2.10
delta G from 2.10 to 2.20
...
delta G from 2.90 to 3.00
Adding these 10 values up will give the
free energy change from 2.00 to 3.00 A.
The system should have been well equilibrated
with the distance restricted at 2.05 A via
RATTLE or IRMDF.
If MM atoms 98 and 100 are defined as QM atoms
in $DATA, the mean force is for two QM atoms.
**** cell-list and fast-list ****
QuanPol uses a standard cell-list scheme to generate a
large neighbor list, which is typically 2.0 times larger
than the small neighbor list and has a long updating
cycle like 55 fs. The small list can be efficiently and
frequently (e.g. every 11 fs) generated from the large
list. QuanPol uses an automatic method to determine when
to update a neighbor list. The atoms displace more than
0.2 and less than 0.9 of the buffer width are stored in
7 lists called fast-lists. When there are ~100 atoms in
the 4th fast-list, which stores atoms that have displaces
between 0.5 and 0.6 of the buffer width, it is fairly
quick to check the pair distances between the atoms in all
fast-lists. New atom pairs are added to the current list
to avoid an immediate update, unless the number of atoms
in the 4th fast-list exceeds MXCHECK (typically 300).
For an equilibrium system, the frequencies of updating
the large and small neighbor lists are almost constants.
QuanPol identifies the frequencies and skips unnecessary
checking of the fast-lists. For example, when BUFWID1
=1.0 A and BUFWID2=4.0 A are used, the lists update every
~55 and ~11 MD steps (DT=1 fs) for a PBC system with 9121
protein atoms, 45 ions, and 60759 water atoms at T=310 K,
P=1 bar and V=88.77**3 A**3. In this case, it is safe to
skip the first 48 steps [estimated as NINT(55-SQRT(55))]
for the large list and the first 8 steps for the small
list.
Fast-list updating information is printed in the dat file
for the first 10,000 MD steps.
It is rare to use the following keywords to change the
default values.
MXCHECK= maximum number of atoms to be checked for the
4th fast-list, which stores atoms that have
displaces between 0.5 and 0.6 of the buffer
width. Default=100, maximum=300.
MXCHECK=1 is essentially the CHARMM heuristic
method.
MXLIST2= maximum number of neighbor MM atoms around a
given MM atom in the large neighbor list.
Default=3400 is good for SWRB=12.0 and BUFWID2
=4.0 A. MXLIST2 can be estimated as
((SWRB+BUFWID2)**3)*3/4.
BUFWID2= The width of the buffer region for the large
neighbor list. This width is added to SWRB to
define the sphere. Default=4.0 A is good for
water and biological systems consisting of water.
3.0-6.0 A are reasonable values for SWRB=12.0 A.
It is good to have BUFWID2 > BUFWID1 + 3.0 A.
If BUFWID2 equals to BUFWID1, only one neighbor
list (with MXLIST2) will be used.
MXLIST1= maximum number of neighbor MM atoms around a
given MM atom in the small neighbor list.
Default=1700 is good for SWRB=12.0 and BUFWID1
=1.0 A. MXLIST1 can be estimated as
((SWRB+BUFWID1)**3)*3/4.
BUFWID1= The width of the buffer region for the small
neighbor list. This width is added to SWRB to
define the sphere. Default=1.0 A is good for
water and biological systems consisting of water.
1.0-2.0 A are reasonable values for SWRB=12.0 A.
**** long-range interactions ****
ISWITCH= selects switching function, which works in the
tail region from SWRA to SWRB (or SWRAQ to
SWARBQ). For standard QM methods (HF, GVB,
MCSCF, DFT, TDDFT, MP2) ISWITCH=0,1,2 can be
used and the default is 1. For semi-empirical
QM methods (AM1, PM3 and DFTB), ISWITCH=0,2 can
be used and the default is 2 (even if ISWITCH=1
is input it will be set to 2).
ISWITCH=1 should be used if the MM perturbation
on QM electronic property is important.
ISWITCH=2 (QM-MM coupling is affected by ISHIFT)
may result in unnatural MM perturbation on QM,
especially when QM (or a QM group) has charges.
ISWITCH=0 is set for continuum solvation models.
= 0 no switching function
= 1 atom-atom switching function for LJ (for both
MM-MM and QM-MM, from SWRA to SWRB);
atom-atom shifting (ISHIFT) function for MM
charge-charge interactions;
QMcenter-MMatom switching function for
QM-charge interaction (from SWRAQ to SWRBQ);
If IPOLSHF=0, atom-atom switching function
is also used for QM-pol, charge-pol and pol-
pol interactions (from SWRA to SWRB).
= 2 atom-atom switching function for LJ (for both
MM-MM and QM-MM, from SWRA to SWRB);
atom-atom shifting (ISHIFT) function for MM
charge-charge interactions;
atom-atom shifting (ISHIFT) function for
QM-charge interaction (from 0 to SWRB);
If IPOLSHF=0, atom-atom switching function
is also used for QM-pol, charge-pol and pol-
pol interactions (from SWRA to SWRB).
The switching function implemented in QuanPol is
W(r) = 1 - 10*D**3 + 15*D**4 - 6*D**5
with
D=(r**2 - SWRA**2)/(SWRB**2 - SWRA**2)
or
D=(r**2 - SWRAQ**2)/(SWRBQ**2 - SWRAQ**2)
ISHIFT = selects shifting function (default=4). The order
of aggressiveness in shifting is 1 > 2 > 3 > 4.
Shifting functions operate on the range zero to
SWRB for MM charge-charge interaction, but may
also affect QM-MM charge interactions (see
ISWITCH). If IPOLSHF=1 is specified, shifting
function is also used for charge-pol and pol-pol
interactions.
ISHIFT=0 is set for continuum solvation models.
= 0 no shifting function
= 1 use the atom-atom shifting function
S(r)=(1-r/SWRB)**2
This shifting function is used by the ENCAD
and ilMM codes.
= 2 use the atom-atom shifting function
S(r)=1-[3*RXNEPS/(2*RXNEPS+1)]*(r/SWRB)+
[(RXNEPS-1)/(2*RXNEPS+1)]*[(r/SWRB)**3]
to mimic a dielectric reaction field.
RXNEPS is required. This shifting function is
Eq (5) in Rick, J.Chem.Phys. 120, 6085 (2004)
= 3 use the simple atom-atom level shifting:
S(r)=(1-r/SWRB)
This is not a smooth function.
= 4 use the atom-atom shifting function
S(r)=[1-(r/SWRB)**2]**2
This is one of the CHARMM shifting functions.
For dipolar bulk systems, if Ewald summation is not used,
a shifting function (rather than a switching function)
should be used for charge-charge interaction (otherwise
structures and energies may be wrong). Many force fields,
especially water models, are optimized for particular
shifting functions, switching functions, and cutoff
distances. Very different results may be obtained when
different shifting and switching functions are used.
Note in QuanPol only ISHIFT is implemented for MM charge
interactions (ISWITCH has no effect on them).
For relative energy or free energy calculations, it is
almost meaningless to use different settings in switching
and shifting functions.
IPOLSHF= select atom-atom shifting function for
MM charge-pol and pol-pol interactions.
= 0 use switching function (default)
= 1 use shifting function. This is not
recommended because induced dipole energy is
sensitive to shifting functions. Induced
dipole energy is much less sensitive to
switching functions because they only work
at far distances.
SWRA =
SWRB = distance cutoffs for the switching function
that gradually drops the interactions from full
strength at SWRA to zero at SWRB. In angstrom.
SWRB is also the cutoff for shifting functions.
Default SWRA=10 A, SWRB=12 A when PBC is used.
Defaults are huge values when PBC is not used.
SWRAQ =
SWRBQ = same as SWRA and SWRB, but for QM-MM switching
function (QMcenter-MMatom).
SWRAQ and SWRBQ should be as large as possible.
The defaults are 10 A and 12 A when PBC is used,
and are huge values when PBC is not used. For
protein calculations with a large QM region, 22
A and 32 A may be good for ISWITCH=1, 22 A and
24 A may be good for ISWITCH=2.
IEWALD = request Ewald summation for PBC charge-charge
interaction. Only charge-charge is implemented,
with the tin-foil conductor boundary condition.
Works only for neutral and pure MM systems.
Also works for MM IFEPTYP=1 and 2.
= 0 no Ewald summation (default)
= 1 use cubic Ewald summation
= 2 use near-spherical Ewald summation, 2~3
times faster than IEWALD=1 (recommended)
KEWALD = the number of cubic or spherical shells in Ewald.
Often called K-vector in the literature.
Default=10 (should increase for XBOX > 30 A).
Maximum 100.
When 10 shells are used, there are 9261 boxes
for IEWALD=1 and 5833 boxes for IEWALD=2,
including the master box. The direct charge-
charge interaction (i.e. real space sum) is
calculated within the master box (i.e. minimum
image convention) and with a cutoff = SWRB,
which is typically 12.0 A (22.68 bohr). See:
KEWALD = 5 10 20 40
IEWALD 1 # boxes= 1331 9261 68921 531441
IEWALD 2 # boxes= 967 5833 39913 293621
SPLIT = the Ewald splitting parameter in the Gauss error
function ERF(SPLIT*R). Default 0.15 bohr**(-1)
is good for SWRB = 12.0 A (22.68 bohr) because
ERFC(0.15*22.68) = 1.5D-06.
Larger SPLIT, smaller SWRB, larger KEWALD.
Smaller SPLIT, larger SWRB, smaller KEWALD.
For bulk water, when IEWALD=2 SWRB=12 SPLIT=0.15 are used,
the following settings can likely converge the Ewald
energy to within 0.1 kcal/mol:
XBOX = 25 50 75 100 125 150 (in Ang)
KEWALD = 6 14 22 31 41 51
For a given system, inclreasing KEWALD by 2 can typically
decrease the error of its Ewald energy by 10 times.
**** continuum solvation models ****
RXNEPS = dielectric constant in ISPHSOL, IFIXSOL and
ISHIFT=2 calculations. Default=78.39.
IFIXSOL= enable the FixSol solvation model
calculation, which is available for QM/MM and
pure MM systems. FixSol paper:
Thellamurege and Li, JCP 137, 246101 (2012)
FixSol is equivalent to CPCM or COSMO, but uses
the FIXPVA2 tessellation scheme. FixSol works
for HF, DFT, DFTB, GVB, MCSCF, TDDFT, TDDFTB,
and MP2.
= 0 skip (default)
= 1 perform FixSol calculation
When FixSol is used, PBC and switching/shifting
functions are turned off automatically.
See J. Chem. Theo. & Comp., 2015, 11(9)
4220-4225 for better choices of the scaling
factor used in COSMO and CPCM. It may also work
well for FixSol.
FIXTOL = convergency criterion in FixSol iterative
calculation of surface charges.
Default=1.0D-10 e is almost always good.
For large systems, FIXTOL=1.0D-06 e may be used.
MXFFTS = maximum number of surface tesserae to be used in
FixSol calculation. Default is usually enough.
NTSATM = number of surface tesserae per atom to be used in
FixSol calculation. Only 60, 240 and 960 are
allowed. Default=60. FixSol uses the FIXPVA2
tessellation method.
By default, FixSol uses a set of simplified
united atomic radii (SUAR):
H 0.000 A
Li - B 1.400 A
C 2.100 A
N 2.000 A
O 1.900 A
F - Al 1.800 A
Si 2.000 A
P 2.200 A
S 2.400 A
Cl 2.760 A
Ar 3.000 A
All others 2.400 A
NRADQM and NRADMM values override the RALLQM,
RALLMM or SUAR defaults, and NRADQM overrides
NRADMM. The override order is:
NRADQM > NRADMM > RALLQM > RALLMM > SUAR
RALLMM = FixSol radii for all heavy MM atoms in $FFDATA.
Default = 0.0 A, use SUAR.
RALLQM = FixSol radii for all heavy QM atoms in $DATA.
The capping QM H atoms in QM/MM systems will be
treated as heavy atoms. Default = 0.0 A, use
SUAR.
NRADMM = n, I1, R1, I2, R2, ... In, Rn
= specify the FixSol radii (in angstrom) for up
to 200 MM atoms in $FFDATA.
n = number of atoms (default=0)
In = MM atom sequential number in $FFDATA
Rn = radius (e.g. 0.001, 1.80, 500.0)
For example, NRADMM=2 500 2.0 502 2.5 is to
assign the 500th MM atom with 2.0 A radius and
the 502nd MM atom with 2.5 A radius.
NRADQM = n, I1, R1, I2, R2, ... In, Rn
= specify the FixSol radii (in angstrom) for up
to 200 QM atoms in $DATA.
n = number of atoms (default=0)
In = QM atom sequential number in $DATA
Rn = radius (e.g. 0.001, 1.80, 500.0)
For example, NRADQM=2 5 1.7 6 1.9 is to
assign the 5th QM atom with 1.7 A radius and
the 6th QM atom with 1.9 A radius.
** Spherical boundary condition and **
** SphSol have strong surface effect **
** Do not use them **
SPHRAD = radius of the sphere containing the QM/MM system.
Default is a huge value, meaning no sphere.
A Lennard-Jones type potential is applied to keep
the heavy atoms in the sphere. For each atom:
V=4*SPHEPS*{[SPHSIG/(r-R)]**12 - [SPHSIG/(r-R)]**6}
R= SPHRAD + [2**(1/6)-1]*SPHSIG
V= -SPHEPS when r = SPHRAD - SPHSIG
SPHEPS = Lennard-Jones epsilon parameter for SPHRAD.
Default=0.15 kcal/mol is good for water.
Proper values should be determined empirically.
SPHSIG = Lennard-Jones sigma parameter for SPHRAD.
Default=1.5 A is good for water.
Proper values should be close to the radii of
the solvent atoms, which are usually around 1.5.
ISPHSOL= enable spherical solvation model (SphSol)
= 0 no SphSol (default)
= 1 image charge method, currently only for
pure MM system
= 60, 240, 960, 3840 to choose surface charge
method and define the number of surface
elements. Available for MM and QM/MM.
When SphSol is used, PBC and switching/shifting
functions are turned off automatically.
RSPHSOL= radius of sphere in angstrom (default=1.0D+30)
used in the SphSol calculation.
SPHRAD is also required. For water solvent,
RSPHSOL = SPHRAD + 0.60 A
RXNEPS = 78.39
SPHEPS = 0.15
SPHSIG = 1.50
are strongly suggested.
**** MD properties ****
NRDF = n, NAME1, NAME2, ...
= specifies the number of pairs for the radial
distribution function calculation, and the names
of the atoms. Must give n pairs of names.
This option works for both periodic and spherical
boundaries (defined by XBOX and SPHRAD).
Default n = 0.
The RDF is calculated at every MD step but
printed out every JOUT steps.
NRDEN = n, NAME1, NAME2, ...
= specifies the number of atoms for the radial
density profile calculation, and the names of
the atoms. Must give n names (default n = 0).
The profile is calculated at every MD step but
printed out every JOUT steps. Works only for
spheric boundary condition (SPHRAD).
DELRDF = specifies the radial increment in the radial
distribution function calculation (NRDF) and the
radial density profile (NRDEN) calculation.
Default=0.05 angstrom.
DIFFUSE= n, NAME1, NAME2, ...
= specifies the number of atoms for diffusion
coefficient calculation, and the names of the
atoms. Must give n names. Default n=0.
If DIFFUSE=1 COM are given, the center of mass
of the first NATPDB atoms in $FFDATA will be
used to determine the diffusion coefficient
and gas phase ion mobiliy, as well as the ion-
gas collision cross section (CCS relevant to ion
mobility spectrometry, IMS) in MD simulation.
In this case, NATPDB must be given by user
(cannot be determined automatically). In
addition, EFIELDX, EFIELDY, EFIELDZ should be
used to specify the external electric field.
Ref: J. Chem. Phys. 148, 064109, 2018.
Currently only He, Ne, Ar, Kr, H2, N2, O2,
CO2, CH4 can be used as 'single-atom' gas
molecules, and only H and H, N and N, O and O
can be used as 'two-atom' gas molecules, and
only C,O,O and O1,H2,H3 can be used as 'three
-atom' CO2 and H2O gas molecules. A typical IMS
CCS simulation has the following keywords in
$QUANPO:
DIFFUSE=1 COM TIMDFS=2.0D-11 NATPDB=XX
EFIELDX=-1.5D+07
EFIELDY=-1.5D+07
EFIELDZ=-1.5D+07
PRES0=10.0 DT=1.0D-15 TEMP0=290
NSTEP=100000000 JOUT=20000 KOUT=100000
NFFTYP=0 SWRA=10 SWRB=12
ITSTAT=1 IPSTAT=0 IRATTLE=0
ISHIFT=4 IPOLSHF=0 NDIEL=0 IDOPOL=1
XBOX=100.0 YBOX=100.0 ZBOX=100.0
TIMDFS = time interval for diffusion coefficient
calculation.
Default=3.0D-12 second is good for water.
Can be larger, but should not be smaller.
There must be sufficient displacement in order
to apply the statistical formula.
NATPDB = number of atoms in the PDB file (but waters in
PDB are excluded). If $FFPDB is used, NATPDB
will be automatically determined. The main use
is for restart jobs in which only $FFDATA is
provided.
NRIJMM = NRIJMM, I1, J1, I2, J2, ...
= specifies up to 100 pairs of MM atoms to print
out their distances at every JOUT steps.
Works for both MD and OPTIMIZE. Useful when one
wants to monitor H-bond distances. Default
NRIJMM = 0.
NRIJQM = NRIJQM, I1, J1, I2, J2, ...
= specifies up to 100 pairs of QM atoms to print
out their distances at every JOUT steps.
Default NRIJQM = 0.
NAIJKMM= NAIJKMM, I1, J1, K1, I2, J2, K2, ...
= specifies up to 100 sets of MM atoms to print
out their angles (IJK) at every JOUT steps.
Default NAIJKMM = 0.
NAIJKQM= NAIJKQM, I1, J1, K1, I2, J2, K2, ...
= specifies up to 100 sets of QM atoms to print
out their angles (IJK) at every JOUT steps.
Default NAIJKQM = 0.
NRMSD = root-mean-square-displacement calculation
for all NATPDB atoms in $FFPDB or $FFDATA.
Works for both MD and OPTIMIZE.
= 0 skip (default)
= 1 calculate RMSD from the initial coordinates
at every JOUT steps. The average unsigned
displacement is also printed out.
NGYRA = radius of gyration calculation for all
NATPDB and non-hydrogen NATPDB atoms in $FFPDB
or $FFDATA(see TIMGYRA).
Works for both MD and OPTIMIZE.
= 0 skip (default)
= 1 calculate radius of gyration using formula:
R=SQRT[sum(m*r*r)/sum(m)]
r: distance from COM
m: atom mass
So R is mass-weighted RMS distance from COM.
TIMGYRA= time interval for radius of gyration calculation.
Default=1.0D-12 s. Can be larger or smaller.
For OPTIMIZE, it is every JOUT steps.
NRALL = activate internuclear distance calculation
for all NATPDB atoms in $FFPDB or $FFDATA
(see TIMRALL). Works for both MD and OPTIMIZE.
= 0 skip (default)
= 1 calculate internuclear distances and compare
to those in the initial structure. RMS
deviation is printed out at every JOUT steps.
TIMRALL= time interval for internuclear distance
calculation. Default=1.0D-12 second. Can be
larger or smaller, but frequent calculation
slows down the MD. For OPTIMIZE, it is every
JOUT steps.
NDIEL = MD simulation of dielectric constant.
= 0 skip
= 1 calculate dielectric constant for the whole
system, including all QM and MM atoms
(default).
If NATPDB>0, it also calculates dielectric
constant for the subsystem defined by NATPDB
(i.e. a protein or DNA/RNA molecule).
The following formula in atomic units is used:
Eps = 1 + 4*Pi*( - (M)**2)/(3kTV)
M = total dipole moment of the system or
the subsystem (including induced atomic
dipoles) at the center of mass.
k = Boltzmann constant
T = average temperature
V = average volume. For NATPDB atoms, V is
estimated as 6.72 A**3 per atom.
For open systems, the volume is infinite, so the
dielectric constant is 1.
IVIBMM = n, I1, I2, I3, ... In
= specifies up to 200 atoms in $FFDATA to calculate
their center of mass and dipole moment at each MD
step. In addition, the velocities of these MM
atoms and the velocity sum will be printed out
at every MD step.
Default n=0.
If this input is lengthy, use multiple lines
and '>' at the end of each line to glue them
together.
Note that in any case, the dipole moment of all
MM or QM or QM/MM atoms, the velocities of all
QM and IVIBMM atoms, the velocity sums of all MM,
QM, QM/MM, and IVIBMM atoms are always printed
out at every MD step.
The QuanPol Vibrational Spectrum Program can be
used to analyze the time dependence of the dipole
moment and velocities, and generate IR and
vibrational spectra.
**** preparation tools ****
NFOLD = this is used only for $FFDATA to duplicate the
input molecule in 3D space NFOLD times.
Reasonable values are 0, 3, 6, 9, 12 and 15,
which leads to 1, 8, 64, 512, 4096 and 32768
copies. 0, 1, 2, 3, ..., 14, 15 can be used.
Default=0, no action.
RFOLD = specifies the spacing when NFOLD is active.
The value should be typically the cubic root of
the volume of the duplicated molecule, and should
be calculated using density. For example, 3.1,
4.7 and 4.9 A are good for H2O, CH2Cl2 and
CH3COCH3, respectively. Default=0.0 A.
ICOMBIN= combine $FFDATA and $FFDATB to be a new $FFDATA.
This can be used to combine solutes with a box
of solvent molecules prepared using NFOLD, or
to combine two molecules with a covalent bond
between them.
ICOMBIN can combine $FFDATA and $FFDATB from
different force fields (e.g., AMBER + MMFF94),
the combined $FFDATA can be used by QuanPol to
run MM and QM/MM calculations (NFFTYP=0, 20000,
30000, or 40000 must be used). The MMFF94
force field will be 'reduced' to the AMBER or
CHARMM style (12-6 LJ will be used for MMFF94).
This is particularly useful when some small
molecules need an ad hoc force field.
See IDELETE if overlap atoms need to be deleted.
= 0 skip (default)
= 1 combine $FFDATA and $FFDATB, both remain
in their original Cartesian coordinates.
= 2 combine $FFDATA and $FFDATB, and translate
$FFDATB so its geometric center coincides
with that of $FFDATA (move B to match A).
= 3 combine $FFDATA and $FFDATB, between that
there is one covalent bond specified via
the keyword MATCHAB.
MATCHAB= IA1, IA2, IB1, IB2
= specify the sequence numbers of a pair of atoms
forming a covalent bond in $FFDATA and $FFDATB
when ICOMBIN=3 is used.
IA1 and IA2 for the two bonded atoms in $FFDATA.
IB1 and IB2 for the two bonded atoms in $FFDATB.
Atoms IA1 and IB1 should have the same Cartesian
coordinates, so do atoms IA2 and IB2.
When ICOMBIN=3 is used, atoms in $FFDATA are all
deleted if they are on the IA2 side, atoms in
$FFDATB are all deleted if they are on the IB1
side. Covalent terms across this bond is
estimated using existing values in $FFDATA and
$FFDATB.
IDELETE= check the atoms in $FFDATA and delete those
are within 1.0 A to any one of the first n atoms
(n=IDELETE). The atoms forming covalent bonds
with any deleted atoms will also be deleted
(molecule deletion). Default=0, no action.
This can be used to remove overlaping atoms in a
$FFDATA generated from ICOMBIN=1 and 2 (not 3).
ISCOOP = scoop out a subset of atoms/molecules from
a given $FFDATA. The scooped-out atoms are
centered at CENTX, CENTY, CENTZ, which are either
given or determined from the input $FFDATA.
= 0 skip (default)
= 1 scoop out a rectangular box with side lengths
XBOX, YBOX, ZBOX.
= 2 scoop out a sphere with radius = SPHRAD
**** force field files ****
NFFFILE= select force field parameter and topology files
= 0 use no such files (default)
= 2 use parameter and topology files from CHARMM
= 3 use parameter and topology files from AMBER
TOPFILE= path/name of a CHARMM or AMBER GAFF topology
file, in single quotes. For example, if
yyy=/home/user,
'yyy/gamess/auxdata/QUANPOL/top_all27_prot_na.rtf'
'yyy/gamess/auxdata/QUANPOL/top_all36_na.rtf'
'yyy/gamess/auxdata/QUANPOL/top_all36_prot.rtf'
'yyy/gamess/auxdata/QUANPOL/top_amber_cornell.inp'
'yyy/gamess/auxdata/QUANPOL/top_opls_aa.inp'
'yyy/amber-gaff.mol2'
The amber-gaff.mol2 file must be generated by
using AmberTools (http://ambermd.org/), and in
the mol2 format.
There must be no space between 'TOPFILE' & '=',
and the path/name must be in single quotes,
and less than 90 characters.
* Correct examples:
TOPFILE='yyy/gamess/auxdata/QUANPOL/xxx'
TOPFILE= 'yyy/xxx'
* Wrong examples:
TOPFILE ='yyy/gamess/auxdata/QUANPOL/xxx'
TOPFILE='~/gamess/auxdata/QUANPOL/xxx'
TOPFILE=yyy/gamess/auxdata/QUANPOL/xxx
TOPAMIA= path/name of an AMBER topology file for amino
acids, in single quotes. For example, if
yyy=/home/user,
'yyy/gamess/auxdata/QUANPOL/all_amino94.in'
'yyy/gamess/auxdata/QUANPOL/all_amino02.in'
'yyy/gamess/auxdata/QUANPOL/amino10.in'
'yyy/gamess/auxdata/QUANPOL/amino12.in'
See TOPFILE for correct input format.
TOPCTER= path/name of an AMBER topology file for
C-terminal amino acids, in single quotes.
For example, if yyy=/home/user,
'yyy/gamess/auxdata/QUANPOL/all_aminoct94.in'
'yyy/gamess/auxdata/QUANPOL/all_aminoct02.in'
'yyy/gamess/auxdata/QUANPOL/aminoct10.in'
'yyy/gamess/auxdata/QUANPOL/aminoct12.in'
See TOPFILE for correct input format.
TOPNTER= path/name of an AMBER topology file for
N-terminal amino acids, in single quotes.
For example, if yyy=/home/user,
'yyy/gamess/auxdata/QUANPOL/all_aminont94.in'
'yyy/gamess/auxdata/QUANPOL/all_aminont02.in'
'yyy/gamess/auxdata/QUANPOL/aminont10.in'
'yyy/gamess/auxdata/QUANPOL/aminont12.in'
See TOPFILE for correct input format.
TOPNUCA= path/name of an AMBER topology file for nucleic
acids, in single quotes. For example, if
yyy=/home/user,
'yyy/gamess/auxdata/QUANPOL/all_nuc94.in'
'yyy/gamess/auxdata/QUANPOL/all_nuc02.in'
'yyy/gamess/auxdata/QUANPOL/nucleic10.in'
See TOPFILE for correct input format.
PARFILE= path/name of a CHARMM/AMBER/MMFF parameter file,
in single quotes. For example, if
yyy=/home/user,
'yyy/gamess/auxdata/QUANPOL/par_all27_prot_na.prm'
'yyy/gamess/auxdata/QUANPOL/par_all36_prot.prm'
'yyy/gamess/auxdata/QUANPOL/par_all36_na.prm'
'yyy/gamess/auxdata/QUANPOL/par_amber_cornell.inp'
'yyy/gamess/auxdata/QUANPOL/par_amber_98.inp'
'yyy/gamess/auxdata/QUANPOL/par_opls_aa.inp'
'yyy/gamess/auxdata/QUANPOL/parm91.dat'
'yyy/gamess/auxdata/QUANPOL/parm94.dat'
'yyy/gamess/auxdata/QUANPOL/parm96.dat'
'yyy/gamess/auxdata/QUANPOL/parm98.dat'
'yyy/gamess/auxdata/QUANPOL/parm99.dat'
'yyy/gamess/auxdata/QUANPOL/parm10.dat'
'yyy/gamess/auxdata/QUANPOL/parm14ipq.dat'
'yyy/gamess/auxdata/QUANPOL/gaff.dat'
'yyy/gamess/auxdata/QUANPOL/MMFF-I_AppendixB.ascii'
See TOPFILE for correct input format.
PARFIL2= path/name of a second AMBER parameter file
frcmod.* that is to add and replace parameters
in regular parameter file parm**.dat.
For example, if yyy=/home/user,
'yyy/gamess/auxdata/QUANPOL/frcmod.ff99SB'
'yyy/gamess/auxdata/QUANPOL/frcmod.ff12SB'
'yyy/gamess/auxdata/QUANPOL/frcmod.ff14SB'
'yyy/gamess/auxdata/QUANPOL/frcmod.ff02pol.r1'
'yyy/gamess/auxdata/QUANPOL/frcmod.parmbsc0'
See TOPFILE for correct input format.
PARFIL3= path/name of a third AMBER parameter file
frcmod.* that is to add or replace parameters
in regular parameter file parm**.dat and
PARFIL2. This is seldom used.
LJSIGMA= select the use of sigma or Rmin/2 for LJ in the
input and output $FFDATA (and $FFDATB). This is
only for I/O purposes.
= 0 use Rmin/2 (default)
= 1 use sigma, which is 1.781797436280679*Rmin/2
WT14LJ = scaling factor for 1-4 Lennard-Jones interaction.
Default=1.00.
For CHARMM, QuanPol uses an additional set of
parameters for 1-4 LJ interaction. In this case
WT14LJ must be 1.00. If not, only the first set
of LJ parameters will be used, and scaled by
WT14LJ for 1-4 cases.
For AMBER and OPLSAA, QuanPol has two ways to
scale the 1-4 LJ interaction by 0.50:
1. Use WT14LJ = 1.00 but an additional set of
pre-scaled LJ parameters. (default)
2. Use WT14LJ = 0.50. The additional set of LJ
parameters is not used.
For MMFF94, the default and only option is 1.00.
Users can input WT14LJ to override the defaults.
WT14CH = scaling factor for 1-4 charge-charge interaction.
Default=1, 1/1.2, 1/2, 3/4 for CHARMM, AMBER,
OPLSAA, MMFF94, respectively, and = 1 for other
cases.
Users can input WT14CH to override the defaults.
**** others ****
IDOCHG = include MM charges
IDOLJ = include MM Lennard-Jones
IDOCMAP= include CHARMM CMAP for proteins
For all of these,
= 1 include (default)
= 0 exclude
IDOPOL = specify how to include induced dipoles. For
large systems, IDOPOL=1 is ~2 times slower than
IDOPOL=0, and IDOPOL=100 is ~10 times slower
than IDOPOL=0. Most induced dipole models are
parameterized using IDOPOL=100, and must use
IDOPOL=100. Only those parameterized using
IDOPOL=1 can use IDOPOL=1. For the same system,
IDOPOL=1 gives 85%~90% of the polarization
energy as compared to IDOPOL=100.
= 0 do not include
= 1 dipoles are induced by external field due
to MM charges, QM nuclei and electrons, and
induced surface charges. No interaction
between induced dipoles are considered and
no iteration is required, thus very fast.
= 100 Dipoles are induced by external field and
the field due to other induced dipoles.
It requires many iterations (maximum=100).
ITYPWAT=302 and NFFTYP=30000 (polarizable
version from 2002) should use IDOPOL=100
(default).
POLTOL = convergency criterion in induced dipole iterative
calculation when IDOPOL=100. Default=1.0D-09
e*bohr.
IPODAMP= specify methods for damping interactions between
induced dipoles at short distances. Damping is
necessary only for IPO1213=1.
= 0 no damping (default)
= 1 linear Thole model (energy not conserved)
= 2 exponential Thole model
= 3 Tinker-exponential model (Thole-Amoeba)
For details of these methods, see Eq 41, 42, 43
in J. Phys.: Condens. Matter 21 (2009) 333102
APODAMP= the unitless factor a in the damping formulas
for IPODAMP=1, 2, and 3. Defaults are 2.500,
2.000, and 0.300, respectively.
IPO1213= specify inclusion of 1-2 and 1-3 interactions
of induced dipoles.
= 0 exclude 1-2 and 1-3 pairs (default)
= 1 include 1-2 and 1-3 pairs
Inclusion of 1-2 and 1-3 interactions often
requires the use of IPODAMP=1,2,3 and is
typically 2~3 times slower than excluding them,
due to the stronger couplings between induced
dipoles. Induced dipoles may have difficulty
to converge if the factor a (see APODAMP) is too
small for IPODAMP=1 or too large for IPODAMP=
2 and 3.
IDIMER = request a fragmentation calculation of the QM
region. Works for RHF, ROHF, UHF types of
wavefunction (and RMP2, ROMP2, UMP2, DFT and
TDDFT) in QM/MMpol/FixSol, and can run geometry
optimization, saddle point search, Hessian, MD.
GVB and MCSCF wavefunctions are not implemented.
AM1, PM3, DFTB/TDDFTB are not implemented.
Free energy perturbation (IFEPTYP=1,2) is not
implemented. LJQMMM=1 must be used.
= 0 no fragmentation of the QM region (default)
= 1 use fragmentation and the sum of the
monomer energy as the total energy. This
is only for speacial tests.
= 2 use fragmentation and the dimer-level
energy as the total energy. This lacks of
many-body polarization correction, thus
not recommended.
= 3 use fragmentation and the dimer-level
energy with the many-body polarization
correction as the total energy, thus the
right choice.
If IDIMER=1, 2 or 3 is specified, the following keywords
(NMOLE, MCHARG, MMULT, IBREAK, RDIMER) are required. In
addition, the MMFF94 force field parameter file is also
required in the $QUANPO group:
PARFILE='your-path/MMFF-I_AppendixB.ascii'
This file is not included in ~/GAMESS/auxdata/QUANPOL/
NMOLE = n, I1, I2, ..., In
= specify how many fragments the QM region is
divided for IDIMER=1,2,3 calculation. I1, I2
tells how many QM atoms (in the $DATA group)
are in each fragment. For example,
NMOLE=5, 10 7 13 3 6
is to use 5 fragments (monomers), the first
one is the first 10 atoms in $DATA, the second
one is the next 7 atoms, the third one is the
next 13 atoms, the fourth one is the next 3
atoms, and the fifth one is the last 6 atoms.
Atoms in $DATA must be in this order, and all
atoms in $DATA must be specified.
The maximum number of fragment, n, is 1024.
The default=0. When n=0 or -1 is used, I1,
I2, ... should not present.
Use IDIMER=3 and NMOLE=0 and IBREAK to
automatically generate the NMOLE list for a
given $DATA group of QM atoms. Using IDIMER=3
and NMOLE=-1, together with a $DATA containing
one and only one QM atom (roughly the central
atom of the desired QM region, such as a Zn2+
ion at the catalytic site) will generate a full
set of QM atoms (printed out as a $DATA group
in the *.log file) that are grouped into small
fragments. In general, atoms in covalent bonds
are grouped into the same fragment. Use IBREAK
to specify covalent bond cuts, or to link non-
bonded groups (e.g., several water molecules).
The charges and multiplicity of each fragment
are also automatically determined, but may need
manual edits because no one knows the exact
values except the user.
For TDDFT, only the first monomer in $DATA
group will actually run TDDFT, all the other
monomers are reduced to DFT (no excitation).
When use IDIMER=3 and NMOLE=0 to group the
$DATA atoms, manually move one atom in the
TDDFT monomer to be the first atom in $DATA
will guarantee the monomer be the first one.
MCHARG = n, I1, I2, ..., In
= specify the charge of each fragment for IDIMER.
MCHARG=5, 0 0 1 0 -2
is to tell that the first, second and fourth
fragments are neutral, the third +1 e charged,
and the last one -2 e charged. These charges
determine the number of electrons in the QM
calculation of each fragment.
Use IDIMER=3 and NMOLE=0 or -1 to automatically
generate the MCHARG list (may need IBREAK
and manual edits for some metal ions).
MMULT = n, M1, M2, ..., Mn
= specify the electronic multiplicity of each
fragment for IDIMER.
MMULT=5, 1 2 1 1 1
is to tell that all the five fragments have S=0
except for the second one.
For fragments with open bonds, capping H atoms
will be automatically added so these open bonds
are closed-shell rather than open-shell. The
MMULT values should be for the fragments with
capping H atoms (thus typically 1).
Use IDIMER=3 and NMOLE=0 or -1 to automatically
generate the MMULT list (may need IBREAK
and manual edits for some metal ions).
IBREAK = m, I1 J1, I2 J2, ..., Im Jm
= provide m pairs of atoms in the $DATA group
if fragmentation is intended for the covalent
bonds (e.g. C-C, C-N, C-O) between these pairs
of atoms in the IDIMER=1,2,3 calculation.
Of course, these atoms must be in bonds.
Note that metal ions (e.g., Li+, Na+, K+, Fe2+
Cu+, Zn2+) are considered as 'isolated' groups
and their coordinate bonds are ignored. To
enforce metal ion coordinate bonds, specify
them in the IBREAK list using negative sign.
IBREAK=2, 30 28 -60 -62
is to force a covalent bond cut at atoms 30
and 28, and force a covalent bond between atoms
60 and 62 (not necessarily involving metal
ions, can be two water molecules).
Up to 40 pairs of atoms can be specified.
Use IDIMER=3 NMOLE=0 and a IBREAK list to
automatically generate the fragments and a new
$DATA group with atoms re-ordered.
Capping H atoms will be used to fill the open
bonds via IBREAK to facilitate the QM and MM
calculation. However, a user will not notice
any thing behind the screen, unless NSTEP=0
is specified to print out some details to the
*.log and *.dat files.
RDIMER = the cutoff distance (A) in IDIMER QM calculation
of pairs of fragments (monomers). No real QM
calculation will be performed for a pair of
monomers if their shortest atomic distance is
beyond RDIMER. Default = 4.0 angstrom. Using a
large value (e.g., 10000.0 A) will enforce all
dimer pairs be calculated via QM. The
interaction between distant pairs are estimated
with an ad hoc force field (only LJ, charge and
induced dipole).
Use $END or a $END line to end $QUANPO.
2315 lines are written.
Edited by Shiro KOSEKI on Tue May 17 15:19:38 2022.