 |
|
|
$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 amove 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.
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=1, 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 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, or ISOFTCR>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.
= 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.
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.
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 IFEPTOP=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 IFEPTOP=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.
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.
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.
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.
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.
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)
= 1 use MM Lennard-Jones terms (default)
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.
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, HESSIAN and MD to run
QuanPol geometry optimization (including 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).
For Hessian jobs, the $FORCE options are not used. 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. For QM/MM systems,
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.
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.
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 =< 10)
= 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=<10) 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 > 10)
= 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 =< 10 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 =< 10)
= 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 < 17.0 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=<10) 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 > 10)
= explicitly define n QM atoms in $DATA to be
active site atoms. This rather lengthy array
should be generated by using LACTQM with
n =< 10 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.
= 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.
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.
KOUTPBC= request coordinates be printed in PBC master box.
= 0 no PBC print out (default)
= 1 $PBCDATA and $PBCFFDATA and $PBCFFDATB are
printed to the dat file every KOUT steps.
These coordinates should not be used to
restart jobs.
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.
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.
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,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, QM/MM and QM/ systems, but
QM atoms will not be rattled unless IRATQM=1.
RATTLE does not affect atoms fixed by NFIXMM,
NFIXMMB, NFIXQM, NFIXQMB, 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, NFIXMMB, NFIXQM, NFIXQMB. 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.
Usuaully 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 (defautl=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.
NFIXQMB= similar to NFIXQM, but for atoms in QM state B,
which are input in $DATA and specified by
MATOMB, MCHARGB and MULTB.
At most 200 QM atoms can be fixed. Default=0.
If both NFIXQM and NFIXQMB are used, the total
number cannot exceed 200.
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.
NFIXMMB= similar to NFIXMM, but for atoms in $FFDATB.
At most 200 MM atoms can be fixed. Default=0.
If both NFIXMM and NFIXMMB are used, the total
number cannot exceed 200.
NFIXQM, NFIXQMB, NFIXMM, NFIXMMB are absolutely enforced,
even IRATTLE and JRATTLE cannot affect them. They also
fix active site atoms defined using LACTMM and LACTQM.
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 and QM/ style 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.
IFEPTOP= specify single or dual topology MD simulation.
Single topology QM/MM system can also run OPT
for IFEPTOP=1 and IFEPTYP=1.
Both IFEPTOP=1 and 2 need $FFDATB to define
state B (the second or target state in FEP), and
need KFREEA and KFREEB to specify the alchemical
atoms in $FFDATA and $FFDATB.
When QM atoms are involved, only IFEPTOP=1 is
allowed, and MATOMA and MATOMB are required to
define QM states A and B.
For IFEPTYP=1, the coordinates of MATOMB atoms
must be the same as those in MATOMA, but may be
less in number (e.g., deprotonated) or smaller
in atomic numbers (e.g., F atoms become H atoms)
to define an alchemical perturbation.
For IFEPTYP=2, the coordinates of MATOMB atoms
may be different from those of MATOMA 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.
= 0 use single topology in $FFDATA (default).
= 1 use single topology scheme, in which only one
set of atoms is used to run MD and sample the
phase space, but the solute atoms in the
KFREEA and KFREEB lists can have different
force field parameters or different (fixed)
coordinates for states A and B.
= 2 use dual topology scheme, in which two sets
of solute atoms coexist in the MD (but those
of A not seeing those of B) and sample
different phase spaces. Soft-core charge
and LJ potentials are usually used to avoid
sampling difficulty arising from singularity.
Dual topology is only implemented for pure MM
systems (no induced dipole) and IFEPTYP=1.
IFEPTYP= specify the type of free energy perturbation
(FEP) calculation (i.e. what free energy or
free energy change to be calculated).
When NFIXMM and NFIXMMB are used to fix two,
three, or more atoms in states A and B, state B
can differ from A in coordinates of the atoms in
NFIXMM and NFIXMMB (all of the other atoms must
have the same coordinates).
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.
= 1 For pure MM, use solvation potential energy
of the solute molecules (i.e. the KFREEA and
KFREEB 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,
and both IFEPTOP=1 and 2 can be used.
= 1 For regular QM/MM, the total QM/MM potential
energy (including solvation and QM internal
energy) is used. The coordinates of MATOMB
atoms must be the same as those of MATOMA.
= 1 For mean-field QM/, the solvation free
energy change is derived from MM simulation,
and corrected by the QM internal energy and
QM/ mean-field electrostatic potential.
The coordinates of MATOMB atoms must be the
same as those of MATOMA.
= 2 For pure MM, use solvation potential energy
of the solute molecules (i.e. the KFREEA and
KFREEB atoms), the covalent potential energy
within the KFREEA atoms (and KFREEB atoms for
B), and the noncovalent potential energy
within the KFREEA atoms (and KFREEB atoms for
B). Require NFIXMM and NFIXMMB.
The following settings are enforced:
WSIMUL=0.0,WPERT1=1.0,WPERT2=1.0,ISOFTCR=0
Only IFEPTOP=1 can be used for IFEPTYP=2.
For pure MM systems, this option is usually
called potential of mean force (PMF).
= 2 For regular QM/MM, the total QM/MM potential
energy (including solvation and QM internal
energy) is used. The coordinates of MATOMB
atoms may be different from those of MATOMA.
= 2 For mean-field QM/, the solvation free
energy change is derived from MM simulation,
and corrected by the QM internal energy and
QM/ mean-field electrostatic potential.
The coordinates of MATOMB atoms may be
different from those of MATOMA.
KFREEA = n, K1, K2, ..., Kn
n specifies the number of atoms in $FFDATA
included in FEP calculation,
K1 - Kn are the sequencial number of the n atoms
in $FFDATA. The limit of n is 500.
If this input is lengthy, use multiple lines
and '>' at the end of each line to glue them
together.
KFREEB = m, K1, K2, ..., Km
m specifies the number of atoms in $FFDATB
included in FEP calculation,
K1 - Km are the sequencial number of the m atoms
in $FFDATB. The limit of m is 500.
For IFEPTOP=1 KFREEB is the same
as KFREEA (it is not necessary to input KFREEB).
If this input is lengthy, use multiple lines
and '>' at the end of each line to glue them
together.
WSIMUL =, WPERT1=, WPERT2=
= three weights (i.e., coupling coefficient
lambda) of the target state B in an FEP
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
ISOFTCR= specify the use of soft-core potentials for
Lennard-Jones, charge-charge interactions in
alchemical free energy perturbation simulation
(only for IFEPTOP=2).
= 0 do not use soft-core (default).
= 1 use soft-core LJ and CH potentials:
LJ=lambda*4*EPS*
{1/[SOFTALJ*(1-lambda)**2+(R/SIG)**6]**2
+1/[SOFTALJ*(1-lambda)**2+(R/SIG)**6]}
CH=lambda*QI*QJ/
SQRT[SOFTACH*(1-lambda)**2+R**2]
lambda = WSIMUL, WPERT1, or WPERT2
SOFTALJ= soft-core parameter alpha for Lennard-Jones.
Default=0.3.
SOFTACH= soft-core parameter alpha for charge-charge.
Default=2.8 A**2.
For IFEPTOP=1, KFREEB must be the same as KFREEA, and
$FFDATB must have the same topology as $FFDATA: same
atoms with the same coordinates (except for atoms fixed
by NFIXMM and NFIXMMB for IFEPTYP=2 jobs), same number
of bonds, angles, dihedrals and other covalent terms.
The parameters (e.g., mass, charge, LJ potential, bond
constant, angle bending constant and others) associated
with alchemical atoms in the KFREEA (=KFREEB) list can be
different in order to define two different states.
IFEPTOP=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).
The best way to create the $FFDATB for IFEPTOP=1 is to
modify the $FFDATA by changing the names, masses, bond
constants, charges and LJ parameters (but never the
coordinates) of the solute atoms in the KFREEA list.
For IFEPTOP=2, KFREEB can be different from KFREEA.
$FFDATB must be very similar to $FFDATA: only the atoms
in KFREEA and KFREEB can be different, and all other
atoms must be the same with the same coordinates.
IFEPTOP=2 MD is performed on the potential energy
surface (PES) for that the covalent potential parameters
(e.g. bond lengths and constants) are in full strength
for both A and B (i.e. ideal-gas-molecule end states),
while charge and LJ potential energies are combined using
E(W)=(1-W)*E(A)+W*E(B). Soft-core charge/LJ potentials
are generally required to achieve better sampling.
The best way to create the $FFDATB for IFEPTOP=2 is to
modify the $FFDATA by changing and/or inserting atoms and
their covalent and noncovalent potentials. The KFREEB
atoms and all of their covalent/noncovalent potentials
specified in $FFDATB are identified and added to
$FFDATA. Since other parts in $FFDATB are not used,
it is not necessary to change them (even though they may
look wrong). One may also use two similar $FFPDB to
create two similar $FFDATA and rename one as $FFDATB.
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.
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, QM/MM and QM/ calculations (MD, OPT).
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.
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.
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 for 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 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.
**** mean field QM/ MD simulation ****
MEANFLD= average position mean field QM/ calculation:
Cui and Li, J. Chem. Phys. 138, 174114 (2013)
= 0 normal QM/MM, no mean field (default)
= n run MM MD simulation for n steps in the
presence of a rigid MM image of the QM
region, store and use the n sets of the MM
coordinates to run a mean field QM/
calculation to obtain QM wavefunction and
energy. For polarizable MM, the (n+1)/2 step
coordinates from the n sets are used to run a
QM/MMpol calculation to obtain polarization
energy. Only energy can be run for QM
atoms in the MM mean field.
MEANFLD can vary from 1 to 20,000, or even
larger if there is enough computer memory.
MEANFLD=10000 and MFMERGE=20 work very well.
MFMERGE= specify how the n (n = MEANFLD) sets of MM
coordinates are merged (averaged) to reduce
the computing time in evaluating QM 1-electron
integrals of the MM charges that are more than
SWRAQ angstrom away from the QM center point.
For MM charges within SWRAQ angstrom, MIN(n,10)
will be used to replace MFMERGE.
= 1 no mergence, so the n sets of MM coordinates
are used explicitly (but very slow).
= m merge every m sets (m.
= 1 use the force field atomic charges.
This method can be very inaccurate.
= 10 use multipole points at each atom (default).
The multipole points are generated with a
density based 3D grid point expansion
method. The solution phase wavefunction is
optimized using the FixSol model.
This option does not need IFIXSOL=1.
QuanPol QM/ MD simulation using MEANFLD=10000:
** QM/ MD step 0 **
1. Initial: (x0_QM, v0_QM, x0_MM, v0_MM)
2. IMMM MD 0 and 1-10000
a. Create an MM image (e.g. charges and LJ potential)
of the QM region using pure QM method (no MM). The
image does not contain polarizability.
b. Using the rigid MM image of the QM, obtain E_IMMM
and forces on all MM atoms. If induced dipoles are
used, they are used only for the MM region, not
for the QM region (i.e. no polarizability for the
image of the QM).
c. MM atoms move. Record 10000 sets of MM coordinates.
Record LJ interaction energy and forces between QM
and MM atoms as E0_. For LJQMMM=1,
E0_ is part of E0_.
d. Report IMMM average energies:
E0_ = E0_ + E0_
+ E0_ + K0_
+ E0_
Tn = K0_QM + Kn_MM (for T scaling)
f. T and P scaling, but QM atoms are not scaled.
3. Run a QM/ calculation to get E0_QM, E0_QM
and E0_QMMMpol. If induced dipoles are used, the
middle step MM coordinates are used to calculate
polarization energy. Now polarization is described
for the QM region using QM method.
4. Report QM/ MD energies:
E0_QM = E0_QM + E0_QM + E0_
+ K0_QM + E0_ + K0_
+ E0_QMMMpol
5. Print out all coordinates and velocity for restart:
x0_QM, v0_QM, x0_MM, v0_MM
** QM/ MD step 1 **
1. No change: (x0_QM, v0_QM, x10000_MM, v10000_MM)
2. IMMM MD 10001-20000
a. Using the rigid MM image of the QM, obtain E_IMMM
and forces on all MM atoms.
b. MM atoms move. Record 10000 sets of MM coordinates.
Record LJ interaction energy and forces between QM
and MM atoms as E1_. For LJQMMM=1,
E1_ is part of E1_.
c. Report IMMM average energies:
E1_ = E1_ + E1_
+ E1_ + K1_
+ E1_
Tn = K0_QM + Kn_MM (for T scaling)
d. TP scaling. QM atoms are not scaled.
3. Run a QM/ calculation to get E1_QM, E1_QM
and E1_QMMMpol.
4. Report QM/ MD energies:
E1_QM = E1_QM + E1_QM + E1_
+ K0_QM + E1_ + K1_
+ E1_QMMMpol
5. Print out all coordinates and velocity for restart:
x0_QM, v0_QM, x10000_MM, v10000_MM
6. T and P scaling, but QM atoms are not scaled.
**** 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.
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 (default=1).
Switching functions work in the tail region,
from SWRA to SWRB.
= 0 no switching function
= 1 atom-atom switching function for LJ;
QMcenter-MMatom switching function for
QM-rep, QM-charge and QM-pol interactions;
If IPOLSHF=0 is specified, atom-atom
switching function is also used for
charge-pol and pol-pol interactions.
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)
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 charge-charge interaction. If IPOLSHF=1
is specified, shifting function is also used for
charge-pol and pol-pol interactions.
= 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 (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.
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
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.
For MM atoms only. 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 interaction.
SWRAQ and SWRBQ should be as large as possible.
The defaults are 10 A and 12 A. For protein
calculations, 22 A and 32 A are good.
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,2 (with IFEPTOP=1).
= 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, GVB, MCSCF, TDDFT, and MP2.
= 0 skip (default)
= 1 perform FixSol calculation
When FixSol is used, PBC and switching/shifting
functions are turned off automatically.
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.
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.
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*( - **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.
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 60 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/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.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 1.00 should be used.
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.
Use $END or a $END line to end $QUANPO.
==========================================================
generated on 7/7/2017