 |
|
|
$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