$QUANPO group (optional, relevant if QuanPol is used) QuanPol (quantum chemistry polarizable force field program) can perform MM, QM/MM or QM/MM/Continuum solvent MD simulation, geometry optimization, saddle point search and Hessian vibration frequency calculation, using HF, DFT, GVB, MCSCF, MP2 or TDDFT wavefunctions. "QuanPol: A full spectrum and seamless QM/MM program", Journal of Computational Chemistry, 2013, 34, 2816-2833. Either $FFDATA (explicit input of coordinates and force field parameters) or $FFPDB (Protein Data Bank coordinates and built-in force fields) needs to be present, to define the MM atoms. Quantum atoms, if any, are given in $DATA as usual. If no quantum atom is used, use COORD=FRAGONLY in $CONTRL. For free energy perturbation calculations, $FFDATA and $FFDATB are required to define the reference and target MM states, and $QUANPO keywords MATOMA and MATOMB are required to define the reference and target QM states. Force field data sets are located by the environment variable QUANPOL. Some parameter and topology files from the CHARMM, AMBER and MMFF94 programs are included in QUANPOL and can be read by QuanPol. When some of the $QUANPO inputs become lengthy, use multiple lines and '>' at the end of each line to glue them together. **** set up MM force field **** MXFFAT = maximum number of MM atoms. MXBOND = maximum number of MM bonds. MXANGL = maximum number of MM bond angles. MXDIHR = maximum number of MM dihedral rotation angles. MXDIHB = maximum number of MM dihedral bending angles. (i.e. improper torsion in CHARMM). MXWAGG = maximum number of MM wagging angles. MXCMAP = maximum number of CHARMM correction map cases. All of the 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