$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