$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