(15 Jul 2024) ********************************* * * * Section 2 - Input Description * * * ********************************* This section of the manual describes the input to GAMESS. The section is written in a reference, rather than tutorial fashion. However, there are frequent reminders that more information can be found on a particular input group, or type of calculation, in the 'Further Information' section of this manual. Numerous complete input files are shown in the 'Input Examples' section. Note that this chapter of the manual can be searched online by means of the "gmshelp" command, if your computer runs Unix. A command such as gmshelp scf will display the $SCF input group. With no arguments, the gmshelp command will show you all of the input group names. Type "" to see the next screen, "b" to back up to the previous screen, and "q" to exit the pager. If gmshelp does not work, ask the person who installed GAMESS to fix the 'gmshelp' script, as it is extremely useful. The order of this section is chosen to approximate the order in which most people prepare their input ($CONTRL, $BASIS/$DATA, $GUESS, and so on). The next few pages contain a list of all possible input groups, grouped in this way. The PDF version of this file contains an index of all group names in alphabetical order. * name function module:routine ---- -------- -------------- Molecule, basis set, wavefunction specification: $CONTRL chemical control data INPUTA:START $SYSTEM computer related options INPUTA:START $BASIS basis set INPUTB:BASISS $DATA molecule, geometry, basis set INPUTB:MOLE $ZMAT internal coordinates ZMATRX:ZMATIN $LIBE linear bend coordinates ZMATRX:LIBE $SCF HF-SCF wavefunction control SCFLIB:SCFIN $SCFMI SCF-MI input control data SCFMI :MIINP $DFT density functional theory DFT :DFTINP $TDDFT time-dependent DFT TDDFT :TDDINP $CIS singly excited CI CISGRD:CISINP $CISVEC vectors for CIS CISGRD:CISVRD $MP2 2nd order Moller-Plesset MP2 :MP2INP $RIMP2 resolution of the identity MP2 RIMP2 :RIDRVR $AUXBAS RI-MP2's basis set specification RIMP2 :RIDRVR $CCINP coupled-cluster input CCSDT :CCINP $RICC resolution-of-the-identity CC RICCRHF:DRRICC_RHF $EOMINP equation-of-motion CC EOMCC :EOMINP $MOPAC semi-empirical specification MPCMOL:MOLDAT $GUESS initial orbital selection GUESS :GUESMO $VEC orbitals (formatted) GUESS :READMO $MOFRZ freezes MOs during SCF runs EFPCOV:MFRZIN $DFTB DFTB input DFTBLB:INPUT $DFTBSK Slater-Koster table input DFTBSK:SKTAB $REKS REKS input REKS :REXINPUT Note that MCSCF and CI input is listed below. Potential energy surface options: $STATPT geometry search control STATPT:SETSIG $TRUDGE nongradient optimization TRUDGE:TRUINP $TRURST restart data for TRUDGE TRUDGE:TRUDGX $FORCE hessian, normal coordinates HESS :HESSX $CPHF coupled-Hartree-Fock options CPHF :CPINP $CPMCHF coupled-MR-Hartree-Fock options MCPCGX:MCPCGX $MASS isotope selection VIBANL:RAMS $HESS force constant matrix (formatted) HESS :FCMIN $GRAD gradient vector (formatted) HESS :EGIN $DIPDR dipole deriv. matrix (formatted) HESS :DDMIN $ALPDR alpha polar. der. (formatted) RAMAN :ADMIN $VIB HESSIAN restart data (formatted) HESS :HSSNUM $VIB2 num GRAD/HESS restart (formatted) HESS :HSSFUL $VSCF vibrational anharmonicity VSCF :VSCFIN $VIBSCF VSCF restart data (formatted) VSCF :VGRID $GAMMA 3rd nuclear derivatives HESS :GAMMXX $EQGEOM equilibrium geometry data HESS :FFCARX $HLOWT hessian data from equilibrium HESS :FFCARX $GLOWT 3rd derivatives at equilibrium HESS :FFCARX $IRC intrinsic reaction coordinate RXNCRD:IRCX $DRC dynamic reaction path DRC :DRCDRV $MEX minimum energy crossing point MEXING:MEXINP $CONICL conical intersection search $MD molecular dynamics trajectory MDEFP :MDX $RDF radial dist. functions for MD MDEFP :RDFX $GLOBOP Monte Carlo global optimization GLOBOP:GLOPDR $GLBFRG Monte Carlo atom groups GLOBOP:GLOPDR $GRADEX gradient extremal path GRADEX:GRXSET $SURF potential surface scan SURF :SRFINP $NEB nudged elastic band NEBPATH:NEBPX Interpretation, properties: $LOCAL localized molecular orbitals LOCAL :LMOINP $TRUNCN localized orbital truncations EFPCOV:TRNCIN $ELMOM electrostatic moments PRPLIB:INPELM $ELPOT electrostatic potential PRPLIB:INPELP $ELDENS electron density PRPLIB:INPELD $ELFLDG electric field/gradient PRPLIB:INPELF $POINTS property calculation points PRPLIB:INPPGS $GRID property calculation mesh PRPLIB:INPPGS $PDC MEP fitting mesh PRPLIB:INPPDC $MGC mean gradient charges $RADIAL atomic orbital radial data PRPPOP:RADWFN $MOLGRF orbital plots PARLEY:PLTMEM $STONE distributed multipole analysis PRPPOP:STNRD $COMP thermochemical calculation COMP :COMPX $RAMAN Raman intensity RAMAN :RAMANX $NMR NMR shielding tensors NMR :NMRX $MOROKM Morokuma energy decomposition MOROKM:MOROIN $LMOEDA LMO-based energy decomposition MOROKM:MMOEDIN $QMEFP QM/EFP energy decomposition EFINP :QMEFPAX $FFCALC finite field polarizabilities FFIELD:FFLDX $TDHF time dependent HF of NLO props TDHF :TDHFX $TDHFX TDHF for NLO, Raman, hyperRaman TDX:FINDTDHFX Solvation models: $EFRAG use effective fragment potential EFINP :EFINP $EFPARM additional parameters for EFP EFINP :EFPARMINP $FRAGNAME specifically named fragment pot. EFINP :RDSTFR $FRGRPL inter-fragment repulsion EFINP :RDDFRL $EWALD Ewald sums for EFP electrostatics EWALD :EWALDX $MAKEFP generate effective fragment pot. EFINP :EFPX $PRTEFP simplified EFP generation EFINP :PREFIN $DAMP EFP multipole screening fit CHGPEN:CGPINP $DAMPGS initial guess screening params CHGPEN:CGPINP $PCM polarizable continuum model PCM :PCMINP $PCMGRD PCM gradient control PCMCV2:PCMGIN $PCMCAV PCM cavity generation PCM :MAKCAV $TESCAV PCM cavity tesselation PCMCV2:TESIN $REORG solvent reorganization in IEF-PCM REORG :RORGIN $NEWCAV PCM escaped charge cavity PCM :DISREP $IEFPCM PCM integral equation form. data PCM :IEFDAT $PCMITR PCM iterative IEF input PCMIEF:ITIEFIN $DISBS PCM dispersion basis set PCMDIS:ENLBS $DISREP PCM dispersion/repulsion PCMVCH:MORETS $SVP Surface Volume Polarization model SVPINP:SVPINP $SVPIRF reaction field points (formatted) SVPINP:SVPIRF $COSGMS conductor-like screening model COSMO :COSMIN $SCRF self consistent reaction field SCRF :ZRFINP Integral, and integral modification options: $ECP effective core potentials ECPLIB:ECPPAR $MCP model core potentials MCPINP:MMPRED $RELWFN scalar relativistic integrals INPUTB:RWFINP $EFIELD external electric field PRPLIB:INPEF $INTGRL 2e- integrals INT2A :INTIN $FMM fast multipole method QMFM :QFMMIN $TRANS integral transformation TRANS :TRFIN Fragment Molecular Orbital method: $FMO define FMO fragments FMOIO :FMOMIN $FMOPRP FMO properties and convergers FMOIO :FMOPIN $FMOXYZ atomic coordinates for FMO FMOIO :FMOXYZ $AFOMOD capping atom input for FMO FMOLIB:MAKELMO $OPTFMO input for special FMO optimizer FMOGRD:OPTFMO $FMOHYB localized MO for FMO boundaries FMOIO :FMOLMO $FMOBND FMO bond cleavage definition FMOIO :FMOBON $FMOENM monomer energies for FMO restart FMOIO :EMINOU $FMOEND dimer energies for FMO restart FMOIO :EDIN $OPTRST OPTFMO restart data FMOGRD:RSTOPT $GDDI group DDI definition INPUTA:GDDINP Polymer model: $ELG polymer elongation method ELGLIB:ELGINP Divide and conquer model: $DANDC DC SCF input DCLIB :DCINP $DCCORR DC correlation method input DCLIB :DCCRIN $SUBSCF subsystem definition for SCF DCLIB :DFLCST $SUBCOR subsystem definition for MP2/CC DCLIB :DFLCST $MP2RES restart data for DC-MP2 DCMP2 :RDMPDC $CCRES restart data for DC-CC DCCC :RDCCDC clusters in molecules $CIMINP controls clusters in molecules CIMINF:CIMINP $CIMATM fine tune calculation level CIMINF:CIMINP $CIMFRG fine tune atomic fragmentation CIMINF:CIMPRT Quantum mechanics/molecular mechanics model: $QUANPO QuanPol calculation QUANPO:QUANPOL $FFDATA QuanPol coordinates for molecules QUANPO:QUANPOL $FFDATB QuanPol coordinates for molecules QUANPO:QUANPOL $FFPDB QuanPol coordinates for proteins QUANPO:QUANPOL MCSCF and CI wavefunctions, and their properties: $CIINP control over CI calculation GAMESS:WFNCI $DET determinant full CI for MCSCF ALDECI:DETINP $CIDET determinant full CI ALDECI:DETINP $GEN determinant general CI for MCSCF ALGNCI:GCIINP $CIGEN determinant general CI ALGNCI:GCIINP $ORMAS determinant multiple active space ORMAS :FCINPT $CEEIS CI energy extrapolation CEEIS :CEEISIN $CEDATA restart data for CEEIS CEEIS :RDCEEIS $GCILST general MCSCF/CI determinant list ALGNCI:GCIGEN $GMCPT general MCSCF/CI determinant list GMCPT :OSRDDAT $PDET parent determinant list GMCPT :OSMKREF $ADDDET add determinants to reference GMCPT :OSMKREF $REMDET remove determinants from ref. GMCPT :OSMKREF $SODET determinant second order CI FSODCI:SOCINP $DRT GUGA distinct row table for MCSCF GUGDRT:ORDORB $CIDRT GUGA CI (CSF) distinct row table GUGDRT:ORDORB $MCSCF control over MCSCF calculation MCSCF :MCSCF $MRMP MRPT selection MP2 :MRMPIN $DETPT det. multireference pert. theory DEMRPT:DMRINP $MCQDPT CSF multireference pert. theory MCQDPT:MQREAD $EXCORR interface to MPQC's R12 program EXCORR:GETEXC $CASCI IVO-CASCI input IVOCAS:IVODRV $IVOORB fine tuning of IVO-CASCI IVOCAS:ORBREAD $CISORT GUGA CI integral sorting GUGSRT:GUGSRT $GUGEM GUGA CI Hamiltonian matrix GUGEM :GUGAEM $GUGDIA GUGA CI diagonalization GUGDGA:GUGADG $GUGDM GUGA CI 1e- density matrix GUGDM :GUGADM $GUGDM2 GUGA CI 2e- density matrix GUGDM2:GUG2DM $LAGRAN GUGA CI Lagrangian LAGRAN:CILGRN $TRFDM2 GUGA CI 2e- density backtransform TRFDM2:TRF2DM $DIABAT diabatic states DIAB:DIABINP $TRANST transition moments, spin-orbit TRNSTN:TRNSTX LibXC and related sections: $LIBXC control LibXC variables LIBXC :LIBXC_INPUT $LDA_K LDA kinetic functionals LIBXC :LIBXC_CHOOSE_FUNCTIONAL $LDA_X LDA exchange functionals LIBXC :LIBXC_CHOOSE_FUNCTIONAL $LDA_C LDA correlation functionals LIBXC :LIBXC_CHOOSE_FUNCTIONAL $LDA_XC LDA XC functionals LIBXC :LIBXC_CHOOSE_FUNCTIONAL $HLDA_XC hybrid LDA XC functionals LIBXC :LIBXC_CHOOSE_FUNCTIONAL $GGA_K control GGA kinetic functionals LIBXC :LIBXC_CHOOSE_FUNCTIONAL $GGA_X GGA exchange functionals LIBXC :LIBXC_CHOOSE_FUNCTIONAL $GGA_C GGA correlation functionals LIBXC :LIBXC_CHOOSE_FUNCTIONAL $GGA_XC GGA XC functionals LIBXC :LIBXC_CHOOSE_FUNCTIONAL $HGGA_X hybrid GGA exchange functionals LIBXC :LIBXC_CHOOSE_FUNCTIONAL $HGGA_XC hybrid GGA XC functionals LIBXC :LIBXC_CHOOSE_FUNCTIONAL $MGGA_K m-GGA kinetic functionals LIBXC :LIBXC_CHOOSE_FUNCTIONAL $MGGA_X m-GGA exchange functionals LIBXC :LIBXC_CHOOSE_FUNCTIONAL $MGGA_C m-GGA correlation functionals LIBXC :LIBXC_CHOOSE_FUNCTIONAL $MGGA_XC m-GGA XC functionals LIBXC :LIBXC_CHOOSE_FUNCTIONAL $HMGGA_XC hybrid m-GGA XC functionals LIBXC :LIBXC_CHOOSE_FUNCTIONAL $HMGGA_X hybrid m-GGA exchange functionals LIBXC :LIBXC_CHOOSE_FUNCTIONAL RISM and related sections: $RISMUN solvent effect using RISM RSM_UNINP :RSM_UNCN $FITINP fitting of cSED RSM_SDDUTL:RSM_FITINP $ABSINP specifies the weight in the fitting RSM_SDDUTL:RSM_ABSIN OpenMP offload and accelerator sections: $DEVICE Accelerator device settings RIMP2GPU:DEVICEINP $OFFLOAD OpenMP target offload controls MOD_OFFLOAD:MOD_OFFLOAD * this column is more useful to programmers than to users. ========================================================== $CONTRL group (note: only one "oh"!) This group specifies the type of wavefunction, the type of calculation, use of core potentials, spherical harmonics, coordinate choices, and similar fundamental job options. Because this is a very long input group, here is a short list of its most important keywords: SCFTYP, MPLEVL, CITYP, CCTYP, DFTTYP, TDDFT RUNTYP, ICHARG, MULT, RELWFN/PP, NZVAR, ISPHER SCFTYP specifies the self-consistent field wavefunction. You may choose from = RHF Restricted Hartree Fock calculation (default) = UHF Unrestricted Hartree Fock calculation = ROHF Restricted open shell Hartree-Fock. (high spin, see GVB for low spin) = GVB Generalized valence bond wavefunction, or low spin ROHF. (needs $SCF input) = MCSCF Multiconfigurational SCF wavefunction (this requires $DET or $DRT input) = REKS Spin-restricted ensemble referenced Kohn-Sham method. (this requires $REKS input) = NONE indicates a single point computation, rereading a converged SCF function. This option requires that you select CITYP=ALDET, ORMAS, FSOCI, GENCI, or GUGA, requesting only RUNTYP=ENERGY or TRANSITN, and using GUESS=MOREAD. The treatment of electron correlation for the above SCF wavefunctions is controlled by the keywords DFTTYP, VBTYP, MPLEVL, CITYP, and CCTYP contained in this group. No more than one of these may be chosen in a single run (except as part of RUNTYP=SURFACE, or if CITYP=SFDET). Scalar relativistic effects may be incorporated using RELWFN for any of these wavefunction choices, correlated or not. DFTTYP = NONE ab initio computation (default) = XXXXXX perform density functional theory run, using the functional specified. Many choices for XXXXXX are listed in the $DFT and $TDDFT input groups. TDDFT = NONE no excited states (default) = EXCITE generate time-dependent DFT excitation energies, using the DFTTYP= functional, for RHF or UHF references. Analytic nuclear gradients are available for RHF. See $TDDFT. = SPNFLP spin-flip TD-DFT, for either UHF or ROHF references. Nuclear gradients and solvent effects are programmed. See $TDDFT. = MRSF Mixed-Reference Spin-Flip (MRSF) TD-DFT for ROHF reference. Nuclear gradients and nonadiabatic MD are programmed for ROHF. See $TDDFT. CUHF option is also utilized to improve the convergence of MRSF. = POL (hyper)polarizability calculation, for RHF only. See $TDDFT. * * * * * PDFTYP = NONE ab initio computation (default) = XXXXXX perform multiconfiguration pair-density functional theory run, using the functional specified. Many choices for XXXXXX are TPBE, FTPBE, TBLYP,FTBLYP,TREVPBE and FTREVPBE. SCFTYP=MCSCF is required. * * * * * VBTYP = NONE no valence bond calculation (default) = VB2000 use the VB2000 program to generate VB wavefunctions, for SCFTYP=RHF or ROHF. Analytic nuclear gradients are not available. A $VB2000 input group is required. See ~/gamess/vb2000/DOC/readme.GAMESS for info about $VB2000, and see also http://www.scinetec.com/~vb * * * * * MPLEVL = chooses Moller-Plesset perturbation theory level, after the SCF. See $MP2, or $MRMP for MCSCF. = 0 skip the MP computation (default) = 2 perform second order energy correction. MP2 (a.k.a. MBPT(2)) is implemented for RHF, UHF, ROHF, and MCSCF wavefunctions, but not GVB. Gradients are available for RHF, UHF, or ROHF based MP2, but for MCSCF, you must choose numerical derivatives to use any RUNTYP other than ENERGY, TRUDGE, SURFACE, or FFIELD. MP2 can also be used on top of CITYP=SFDET to get the SF-MPRMP2 energy correction. * * * * * CITYP = chooses CI computation after the SCF, for any SCFTYP except UHF. = NONE skips the CI. (default) = CIS single excitations from a SCFTYP=RHF reference, only. This is for excited states, with analytic nuclear gradients available. See the $CIS input group. = SFCIS spin-flip style CIS, see $CIS input. = ALDET runs the Ames Laboratory determinant full CI package, requiring $CIDET. = ORMAS runs an Occupation Restricted Multiple Active Space determinant CI. The input is $CIDET and $ORMAS. = FSOCI runs a full second order CI using determinants, see $CIDET and $SODET. = GENCI runs a determinant CI program that permits arbitrary specification of the determinants, requiring $CIGEN. = GUGA runs the Unitary Group CI package, which requires $CIDRT input. Analytic gradients are available only for RHF, so for other SCFTYPs, you may choose only RUNTYP=ENERGY, TRUDGE, SURFACE, FFIELD, TRANSITN. = SFORM performs SF-ORMAS calculation. SCFTYP must be ROHF, and the MULSF must be specified in $CIDET. The active space may be defined as usual in $ORMAS. Analytic gradients are available. = SFDET performs SF-CI calculation using the ALDET code. SCFTYP must be ROHF and MULSF must be specified. This is mostly useful when one wants to perform an SF-MRMP2 calculation. PMTD1 = For CITYP=ALDET or ORMAS, or for these two CI steps in MCSCF runs, for EFP solvent calculations, this flag enables use of "polarization method 1" for the effective fragments. See also FSTATE in $CIDET or $DET = .TRUE. The EFP dipoles will not be re-polarized to the CITYP wavefunction (default) = .FALSE. The EFP dipoles will be re-polarized to the CITYP wavefunction * * * * * CCTYP chooses a ground-state coupled-cluster (CC) or, optionally, equation-of-motion coupled-cluster (EOMCC) calculation for excited, electron-attached, and ionized electronic states using either the RHF or ROHF reference. See also $CCINP and $EOMINP. Only CCSD and CCSD(T) for RHF can run in parallel using DDI. All CC and EOMCC methods may run in parallel on a single SMP node if GAMESS is compiled with a parallel math library, such as MKL. Only CCSD and CCSD(T) calculations can run using the RI approximation for two-electron repulsion integrals and a hybrid MPI/OpenMP parallelization model. The ground-state CC options: = NONE skips CC computation (default). = LCCD perform a CC calculation using the linearized coupled-cluster method with double excitations. = CCD perform a CC calculation using the coupled-cluster method with doubles. = CCSD perform a CC calculation with single and double excitations. = CCSD(T) in addition to CCSD, the non-iterative triples corrections are computed, giving standard CCSD[T] and CCSD(T) energies. = R-CC in addition to all CCSD(T) calculations, the renormalized R-CCSD[T] and R-CCSD(T) energies are computed. = CR-CC in addition to all R-CC calculations, the completely renormalized CR-CCSD[T] and CR-CCSD(T) energies are computed. = CR-CCL in addition to a CCSD ground state and the associated left CCSD calculations, the non-iterative triples energy correction defining the rigorously size-extensive completely renormalized CR-CC(2,3) approach (sometimes also called CR-CCSD(T)_L) is computed. Ground state only (zero NSTATE vector) CCTYP=CR-EOM type CR-EOMCCSD(T) energies and CCSD properties are also generated. For further information about accuracy, and A to D CR-CC(2,3) energy types, see Section 4 of this manual. Generally, CR-CC(2,3),D provides the most accurate energetics. = CCSD(TQ) in addition to all R-CC calculations, non-iterative triples and quadruples corrections are determined, to give CCSD(TQ) and various R-CCSD(TQ) energies. = CR-CC(Q) in addition to all CR-CC and CCSD(TQ) calculations, the CR-CCSD(TQ) energies are obtained. = CCSD3A perform CC calculations with singles, doubles, and active-space triples, designated as CCSDt. When all orbitals in the MO basis are set to be active, CCSDt becomes full CCSDT. = CCT3 in addition to CCSDt, compute the non-iterative CC(P;Q)-based CC(t;3) energy correction due to the missing triples not captured by CCSDt. As for CR-CCL, variant D of CC(t;3) is generally most accurate and closest to CCSDT. = ACCSD perform ACCSD calculations, i.e., CCSD-type calculations that select or scale T2^2 diagrams in the equations projected on doubly excited determinants in the spirit of approximate coupled-pair (ACP) approaches. = ACC23 in addition to ACCSD calculations, compute the CR-CC(2,3)-type non-iterative triples correction to the ACCSD energy. = ACCSD3A perform ACCSDt calculations, i.e., CCSDt-type calculations that select or scale T2^2 diagrams in the equations projected on doubly excited determinants in the spirit of ACP approaches. = ACCT3 in addition to ACCSDt calculations, compute the CC(t;3)-type non-iterative energy correction due to the missing triples not captured by ACCSDt. ACCSD, ACC23, ACCSD3A, and ACCT3 options require the ACP array with T2^2 diagram scaling factors other than 1.0 to be provided in $CCINP. CCSD, CR-CCL, CCSD3A, CCT3, ACCSD, ACC23, ACCSD3A, and ACCT3 options work for both RHF and ROHF references; the remaining ground-state CC calculations can only be performed using SCFTYP=RHF. Excited-state EOMCC options: = EOM-CCSD in addition to a CCSD ground state, excited states are calculated using the equation-of-motion coupled-cluster method with singles and doubles (EOMCCSD). = CR-EOM in addition to CCSD and EOMCCSD, non-iterative triples corrections to CCSD ground-state and EOMCCSD excited-state energies are determined using completely renormalized CR-EOMCCSD(T) approaches. In typical applications, delta-CR-EOMCCSD(T) triples corrections to EOMCCSD excitation energies are most useful. = CR-EOML in addition to printing all results that CR-EOM obtains and solving the left eigenstate EOMCCSD problem used for properties other than energy, the CR-EOMCC(2,3) triples corrections analogous to those of the ground-state CR-CCL approach are calculated. In typical applications, delta-CR-EOMCC(2,3) triples corrections to EOMCCSD excitation energies are most useful. CCTYP=EOM-CCSD is available for RHF or ROHF references; the CR-EOM and CR-EOML triples corrections can only be calculated using SCFTYP=RHF. Single and double electron-attachment and ionization potential EOMCC options for open-shell systems with one or two electrons or holes outside closed-shell cores (available for RHF or ROHF references, although the EA- and IP-EOMCC calculations using ROHF reference are discouraged): = EA-EOM2 electron-attachment EOMCC calculations with up to 2p1h excitations treated fully (i.e., EA-EOMCCSD). = EA-EOM3A electron-attachment EOMCC calculations with all 1p and 2p1h and active-space 3p2h excitations [i.e., EA-EOMCCSDt or EA-EOMCC(3p2h){Nu}]. = IP-EOM2 ionization potential EOMCC calculations with up to 2h1p excitations treated fully (i.e., IP-EOMCCSD). = IP-EOM3A ionization potential EOMCC calculations with all 1h and 2h1p and active-space 3h2p excitations [i.e., IP-EOMCCSDt or IP-EOMCC(3h2p){No}]. = DEAEOM2 double electron-attachment EOMCC with up to 3p1h excitations treated fully [i.e., DEA-EOMCC(3p1h) or DEA-EOMCCSD]. = DEAEOM2A double electron-attachment EOMCC with all 2p and active-space 3p1h excitations [i.e., DEA-EOMCC(3p1h){Nu}]. = DEAEOM3 double electron-attachment EOMCC with up to 4p2h excitations treated fully [i.e., DEA-EOMCC(4p2h); most accurate, but very expensive!]. = DEAEOM3A double electron-attachment EOMCC with all 2p and active-space 3p1h and 4p2h excitations [i.e., DEA-EOMCC(3p1h,4p2h){Nu}]. = DEAEOM3B double electron-attachment EOMCC with all 2p and 3p1h and active-space 4p2h excitations [i.e., DEA-EOMCC(4p2h){Nu}]. = DIPEOM2 double ionization potential EOMCC with up to 3h1p excitations treated fully [i.e., DIP-EOMCC(3h1p) or DIP-EOMCCSD]. = DIPEOM3 double ionization potential EOMCC with up to 4h2p excitations treated fully [i.e., DIP-EOMCC(4h2p)]. = DIPEOM3A double ionization potential EOMCC with all 2h and 3h1p and active-space 4h2p excitations [i.e., DIP-EOMCC(4h2p){No}]. Labels "p" and "h" in the description of (D)EA- and (D)IP-EOMCC methods refer to particles (correlated orbitals unoccupied in the underlying reference) and holes (correlated orbitals occupied in the underlying reference), respectively. Nu and No in the description of (D)EA- and (D)IP-EOMCC methods indicate the numbers of active unoccupied and active occupied orbitals used to select the dominant 3p2h, 3h2p, 3p1h, 4p2h, and 4h2p contributions. EA/IP-EOMCC runs produce ground and excited states of systems obtained by attaching one electron to (EA) or removing one electron from (IP) the underlying closed-shell cores described by CCSD. DEA/DIP-EOMCC runs produce ground and excited states of systems obtained by attaching two electrons to (DEA) or removing two electrons from (DIP) the underlying closed-shell cores described by CCSD. As all other EOMCC calculations, EA/IP- and DEA/DIP-EOMCC runs read $CCINP as well as $EOMINP inputs. Any publication describing the results of CC calculations obtained using GAMESS should reference the appropriate papers, which are listed in the output of every run and in chapter 4 of this manual. Analytic gradients are not available, so use CCTYP only for RUNTYP=ENERGY, TRUDGE, SURFACE, or maybe FFIELD, or request numerical derivatives. Generally speaking, the renormalized CC energies are obtained at similar cost to the standard values, while completely renormalized ones cost about twice the time. For usage tips and more information about resources on the various coupled cluster methods, see Section 4, 'Further Information'. CIMTYP chooses a Cluster-In-Molecule (CIM) calculation. = NONE skip CIM computation, i.e., perform a canonical calculation (default). = SECIM perform a single-environment CIM (SECIM) computation. = DECIM perform a dual-environment CIM (DECIM) computation. = GSECIM perform a generalized SECIM (GSECIM) computation. The $CIMFRG must be included as well. See also $CIMINP and, optionally, $CIMFRG and $CIMATM. If CIMTYP is given, SUBMTD in $CIMINP is required. Only RUNTYP=ENERGY and SCFTYP=RHF or ROHF work when CIMTYP is given. See SUBMTD in $CIMINP for more details. * * * * * RELWFN Selects all-electron scalar relativity treatment. See the $RELWFN input group for more information, including nuclear derivative availability. = NONE use the basic Schrodinger equation (default) = LUT-IOTC local unitary transformation modification of IOTC, due to H.Nakai, J.Seino, Y.Nakajima. This is the fastest and most numerically reliable scalar relativity method, so it is preferred over RESC, DK, or IOTC. = LUTIOTC2 local unitary transformation modification of IOTC2E for one- and two-electron terms, due to H. Nakai and J. Seino. This option is not available to range separated functionals in DFT. = IOTC infinite-order two-component method of M. Barysz and A.J. Sadlej. = IOTC2E infinite-order two-component method including two-electron term of J. Seino and M. Hada, in addition to IOTC. This option is not available to range separated functionals in DFT. = DK Douglas-Kroll transformation, available at the 1st, 2nd, or 3rd order. = RESC relativistic elimination of small component, the method of T. Nakajima and K. Hirao, available at 2nd order only. = NESC normalised elimination of small component, the method of K. Dyall, 2nd order only. * * * * * RUNTYP specifies the type of computation, for example at a single geometry point: = ENERGY Molecular energy. (default) = GRADIENT Molecular energy plus gradient. = HESSIAN Molecular energy plus gradient plus second derivatives, including harmonic harmonic vibrational analysis. Analytic FMO Hessians and PAVE should use RUNTYP=FMOHESS. Semi-numerical FMO Hessians except for PAVE should use RUNTYP=HESSIAN. See the $FORCE and $CPHF input groups. = FMOHESS used for FMO Hessians: analytic for DFTB, RHF, DFT, UHF, and ROHF; numerical runs also possible when using PAVE. FMOHESS reads $FORCE as HESSIAN. = GAMMA Evaluate up to 3rd nuclear derivatives, by finite differencing of Hessians. See $GAMMA, and also NFFLVL in $CONTRL. multiple geometry options: = OPTIMIZE Optimize the molecular geometry using analytic energy gradients. See $STATPT. = TRUDGE Non-gradient total energy minimization. See $TRUDGE and $TRURST. = SADPOINT Locate saddle point (transition state). See $STATPT. = MEX Locate minimum energy crossing point on the intersection seam of two potential energy surfaces. See $MEX. = CONICAL Locate conical intersection point on the intersection seam of two potential energy surfaces. See $CONICL. = IRC Follow intrinsic reaction coordinate. See $IRC. = VSCF anharmonic vibrational corrections. See $VSCF. = DRC Follow dynamic reaction coordinate. See $DRC. = MD molecular dynamics trajectory, see $MD. = GLOBOP Monte Carlo-type global optimization. See $GLOBOP. = OPTFMO genuine FMO geometry optimization using nearly analytic gradient. See $OPTFMO. = GRADEXTR Trace gradient extremal. See $GRADEX. = SURFACE Scan linear cross sections of the potential energy surface. See $SURF. single geometry property options: = COMP composite thermochemistry calculation, including G3MP2. See $COMP input. = G3MP2 evaluate heat of formation using the G3(MP2,CCSD(T)) methodology. See test example exam43.inp for more information. = PROP Molecular properties will be calculated. Orbital localization can be requested as well. See $ELPOT, etc. Converged orbitals must be input in a $VEC input, which suffice to reproduce the wavefunction only for simple SCF: RHF, UHF, ROHF, or DFT counterparts. GVB also works (CICOEF may be needed). All other calculations must instead use RUNTYP=ENERGY to regenerate the density matrix. = RAMAN computes Raman intensities, see $RAMAN. = NACME non-adiabatic coupling matrix element between two or more state averaged MCSCF wavefunctions. The calculation has no specific input group, but must use only SCFTYP=MCSCF with CISTEP=ALDET or ORMAS. = NMR NMR shielding tensors for closed shell molecules by the GIAO method. See $NMR. = EDA Perform energy decomposition analysis. Give one of $MOROKM or $LMOEDA inputs. = QMEFPEA QM/EFP solvent energy analysis, see $QMEFP. = TRANSITN Compute radiative transition moment or spin-orbit coupling. See $TRANST. = FFIELD applies finite electric fields, most commonly to extract polarizabilities. See $FFCALC. = TDHF analytic computation of time dependent polarizabilities. See $TDHF. = TDHFX extended TDHF package, including nuclear polarizability derivatives, and Raman and Hyper-Raman spectra. See $TDHFX. = MAKEFP creates an effective fragment potential, for SCFTYP=RHF or ROHF only. See $MAKEFP, $DAMP, $DAMPGS, $STONE, ... = FMO0 performs the free state FMO calculation. See $FMO. = FCIDUMP Performs an ROHF calculation and generates a FCIDUMP file. = GRADNAC $GRADIENT + NonAdibatic Coupling Term. NACT is designed in cases of using external MD drivers and requires coordinates of current and previous steps. * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Note that RUNTYPs which require the nuclear gradient are GRADIENT, HESSIAN, OPTIMIZE, SADPOINT, GLOBOP, IRC, GRADEXTR, DRC, and RAMAN These are efficient with analytic gradients, which are available only for certain CI or MP2 calculations, but no CC calculations, as indicated above. See NUMGRD. * * * * * * * * * * * * * * * * * * * * * * * * * * * * * NUMGRD Flag to allow numerical differentiation of the energy. Each gradient requires the energy be computed twice (forward and backward displacements) along each totally symmetric modes. It is thus recommended only for systems with just a few symmetry unique atoms in $DATA. FMO uses GDDI in the fragment-way, other runs use GDDI for offloading for offsets. The default is .FALSE. EXETYP = RUN Actually do the run. (default) = CHECK Wavefunction and energy will not be evaluated. This lets you speedily check input and memory requirements. See the overview section for details. Note that you must set PARALL=.TRUE. in $SYSTEM to test distributed memory allocations. = EXPERT Declare your expert status to enjoy special priveleges: some runs will be allowed to proceed that normally stop with the view of protecting novice users. In particular, unnormalized basis sets may be allowed. Also some warnings may be suppressed. = FMOHMO Normalize HMOs in $FMOHYB and quit. = DEBUG Massive amounts of output are printed, useful only if you hate trees. = routine Maximum output is generated by the routine named. Check the source for the routines this applies to. * * * * * * * ICHARG = Molecular charge. (default=0, neutral) MULT = Multiplicity of the electronic state = 1 singlet (default) = 2,3,... doublet, triplet, and so on. ICHARG and MULT are used directly for RHF, UHF, ROHF. For GVB, these are implicit in the $SCF input, while for MCSCF or CI, these are implicit in $DRT/$CIDRT or $DET/$CIDET input. You must still give them correctly. * * * the next three control molecular geometry * * * COORD = choice for molecular geometry in $DATA. = UNIQUE only the symmetry unique atoms will be given, in Cartesian coords (default). = HINT only the symmetry unique atoms will be given, in Hilderbrandt style internals. = PRINAXIS Cartesian coordinates will be input, and transformed to principal axes. Please read the warning just below!!! = ZMT GAUSSIAN style internals will be input. = ZMTMPC MOPAC style internals will be input. = FRAGONLY means no part of the system is treated by ab initio means, hence $DATA is not given. The system is defined by $EFRAG. Note: the choices PRINAXIS, ZMT, ZMTMPC require input of all atoms in the molecule. They also orient the molecule, and then determine which atoms are unique. The reorientation is likely to change the order of the atoms from what you input. When the point group contains a 3- fold or higher rotation axis, the degenerate moments of inertia often cause problems choosing correct symmetry unique axes, in which case you must use COORD=UNIQUE rather than Z-matrices. Warning: The reorientation into principal axes is done only for atomic coordinates, and is not applied to the axis dependent data in the following groups: $VEC, $HESS, $GRAD, $DIPDR, $VIB, nor Cartesian coords of effective fragments in $EFRAG. COORD=UNIQUE avoids reorientation, and thus is the safest way to read these. Note: the choices PRINAXIS, ZMT, ZMTMPC require the use of a group named $BASIS to define the basis set. The first two choices might or might not use $BASIS, as you wish. UNITS = distance units, any angles must be in degrees. = ANGS Angstroms (default) = BOHR Bohr atomic units NZVAR = 0 Use Cartesian coordinates (default). = M If COORD=ZMT or ZMTMPC, and $ZMAT is not given: the internal coordinates will be those defining the molecule in $DATA. In this case, $DATA may not contain any dummy atoms. M is usually 3N-6, or 3N-5 for linear. = M For other COORD choices, or if $ZMAT is given: the internal coordinates will be those defined in $ZMAT. This allows more sophisticated internal coordinate choices. M is ordinarily 3N-6 (3N-5), unless $ZMAT has linear bends. NZVAR refers mainly to the coordinates used by OPTIMIZE or SADPOINT runs, but may also print the internal's values for other run types. You can use internals to define the molecule, but Cartesians during optimizations! * * * * * * * Pseudopotentials may be of two types: ECP (effective core potentials) which generate nodeless valence orbitals, and MCP (model core potentials) producing valence orbitals with the correct radial nodal structure. At present, ECPs have analytic nuclear gradients and Hessians, while MCPs have analytic nuclear gradients. PP = pseudopotential selection. = NONE all electron calculation (default). = READ read ECP potentials in the $ECP input. = SBKJC use Stevens, Basch, Krauss, Jasien, Cundari ECP potentials for all heavy atoms (Li-Rn are available). = HW use Hay, Wadt ECP potentials for heavy atoms (Na-Xe are available). = MCP use Huzinaga's Model Core Potentials. The correct MCP potential will be chosen to match the requested MCP valence basis set (see $BASIS). * * * * * * * LOCAL = controls orbital localization. = NONE Skip localization (default). = BOYS Do Foster-Boys-like localization. = RUEDNBRG Do Edmiston-Ruedenberg localization. = POP Do Pipek-Mezey population localization. = SVD Do single value decomposition, to project the molecular orbitals onto atoms. This is available only for SCFTYP=RHF, ROHF, and MCSCF (full space or ORMAS). The ORIENT keyword in $LOCAL is pertinent. See the related $LOCAL input. Localization is not available for SCFTYP=GVB. DFTB only works with LOCAL=POP (and NONE). * * * * * * * ISPHER = Spherical Harmonics option = -1 Use Cartesian basis functions to construct symmetry-adapted linear combination (SALC) of basis functions. The SALC space is the linear variation space used. (default) = 0 Use spherical harmonic functions to create SALC functions, which are then expressed in terms of Cartesian functions. The contaminants are not dropped, hence this option has EXACTLY the same variational space as ISPHER=-1. The only benefit to obtain from this is a population analysis in terms of pure s,p,d,f,g functions. = +1 Same as ISPHER=0, but the function space is truncated to eliminate all contaminant Cartesian functions [3S(D), 3P(F), 4S(G), and 3D(G)] before constructing the SALC functions. The computation corresponds to the use of a spherical harmonic basis. QMTTOL = linear dependence threshhold Any functions in the SALC variational space whose eigenvalue of the overlap matrix is below this tolerence is considered to be linearly dependent. Such functions are dropped from the variational space. What is dropped is not individual basis functions, but rather some linear combination(s) of the entire basis set that represent the linear dependent part of the function space. The default is a reasonable value for most purposes, 1.0E-6. When many diffuse functions are used, it is common to see the program drop some combinations. On occasion, in multi-ring molecules, we have raised QMTTOL to 3.0E-6 to obtain SCF convergence, at the cost of some energy. MAXIT = Maximum number of SCF iteration cycles. This pertains only to RHF, UHF, ROHF, or GVB runs. See also MAXIT in $MCSCF. (default = 30) * * * interfaces to other programs * * * MOLPLT = flag that produces an input deck for a molecule drawing program distributed with GAMESS. (default is .FALSE.) PLTORB = flag that produces an input deck for an orbital plotting program distributed with GAMESS. (default is .FALSE.) AIMPAC = flag to create an input deck for Bader's Atoms In Molecules properties code. (default=.FALSE.) For information about this program, see the URL http://www.chemistry.mcmaster.ca/aimpac DGRID = flag to add extra digits in molecular orbitals to the log file for use by Kohout's DGrid program: http://www2.cpfs.mpg.de/~kohout/dgrid.html This is one of the modern alternatives to the old AIMPAC codes, in the QTAIM/ELF arena. (default .FALSE.) NUMCOR = an array setting up the number of core orbitals for each element (up to Z=103). NUMCOR is only used in MP2 currently. Default: -1 for all atoms, which means use the default values hardwired in the MP2 code. FRIEND = string to prepare input to other quantum programs, choose from = HONDO for HONDO 8.2 = MELDF for MELDF = GAMESSUK for GAMESS (UK Daresbury version) = GAUSSIAN for Gaussian 9x = ALL for all of the above PLTORB, MOLPLT, and AIMPAC decks are written to file PUNCH at the end of the job. Thus all of these correspond to the final geometry encountered during jobs such as OPTIMIZE, SAPDOINT, IRC... In contrast, selecting FRIEND turns the job into a CHECK run only, no matter how you set EXETYP. Thus the geometry is that encountered in $DATA. The input is added to the PUNCH file, and may require some (usually minimal) massaging. PLTORB and MOLPLT are written even for EXETYP=CHECK. AIMPAC requires at least RUNTYP=PROP. * * * NFFLVL used to determine energies and gradients away from equilibrium structures, at the coordinates given in $DATA. The method will use a Taylor expansion of the potential surface around the stationary point. See $EQGEOM, $HLOWT, $GLOWT. This may be used with RUNTYP=ENERGY or GRADIENT. = 2 uses only Hessian information, which gives a reasonable energy, but not such a good gradient. = 3 uses Hessian and 3rd nuclear derivatives in the Taylor expansion, producing more accurate values for the energy and for the gradient. * * * computation control switches * * * For the most part, the default is the only sensible value, and unless you are sure of what you are doing, these probably should not be touched. NPRINT = Print/punch control flag See also EXETYP for debug info. (options -7 to 5 are primarily debug) = -7 Extra printing from Boys localization. = -6 debug for geometry searches = -5 minimal output = -4 print 2e-contribution to gradient. = -3 print 1e-contribution to gradient. = -2 normal printing, no punch file = 1 extra printing for basis,symmetry,ZMAT = 2 extra printing for MO guess routines = 3 print out property and 1e- integrals = 4 print out 2e- integrals = 5 print out SCF data for each cycle. (Fock and density matrices, current MOs = 6 same as 7, but wider 132 columns output. This option isn't perfect. = 7 normal printing and punching (default) = 8 more printout than 7. The extra output is (AO) Mulliken and overlap population analysis, eigenvalues, Lagrangians, ... = 9 everything in 8 plus Lowdin population analysis, final density matrix. NOSYM = 0 the symmetry specified in $DATA is used as much as possible in integrals, SCF, gradients, etc. (this is the default) = 1 the symmetry specified in the $DATA input is used to build the molecule, then symmetry is not used again. Some GVB or MCSCF runs (those without a totally symmetric charge density) require you request no symmetry. ETOLLZ = threshold to label molecular orbitals by Lz values. Small matrices of the Lz operator are diagonalized for the sets of MOs whose orbital energies are degenerate to within ETOLLZ. This option may be used in molecules with distorted linear symmetry for approximate labelling. Default: 1.0d-6 for linear, 0 (disable) if not. INTTYP selects the integral package(s) used, all of which produce equally accurate results. This is therefore used only for debugging purposes. = BEST use the fastest integral code available for any particular shell quartet (default): s,p,L or s,p,d,L rotated axis code first. ERIC s,p,d,f,g precursor transfer equation code second, up to 5 units total ang. mom. Rys quadrature for general s,p,d,f,g,L, or for uncontracted quartets. = ROTAXIS means don't use ERIC at all, e.g. rotated axis codes, or else Rys quadrature. = ERIC means don't use rotated axis codes, e.g. ERIC code, or else Rys quadrature. = RYSQUAD means use Rys quadrature for everything. GRDTYP = BEST use Schlegel routines for spL gradient blocks, and Rys quadrature for all other gradient integrals. (default) = RYSQUAD use Rys quadrature for all gradient integrals. This option is only slightly more accurate, but is rather slower. NORMF = 0 normalize the basis functions (default) = 1 no normalization NORMP = 0 input contraction coefficients refer to normalized Gaussian primitives. (default) = 1 the opposite. ITOL = primitive cutoff factor (default=20) = n products of primitives whose exponential factor is less than 10**(-n) are skipped. ICUT = n integrals less than 10.0**(-n) are not saved on disk. (default = 9). Direct SCF will calculate to a cutoff 1.0d-10 or 5.0d-11 depending on FDIFF=.F. or .T. ISKPRP = 0 proceed as usual 1 skip computation of some properties which are not well parallelised. This includes bond orders and virial theorem, and can help parallel scalability if many CPUs are used. Note that NPRINT=-5 disables most property computations as well, so ISKPRP=1 has no effect in that case. (default: 0) MODDOS = a bit-additive option to calculate density of states (DOS) + 1: calculate total DOS + 2: automatically set emin and emax to the whole MO range (see DOSPAR(5)) + 4: set the Fermi level to 0 + 8: omit lines with no density Default: 0 DOSPAR = an array of 5 parameters for DOS calculations (1) minimum energy (eV) (2) maximum energy (eV) (3) energy grid step (eV) (4) broadening Gaussian sigma (eV) (5) (1) and (2) are set automatically using MO energies padded with this value (MODDOS=2). Default: 0,0,0,0,0 * * * restart options * * * IREST = restart control options (for OPTIMIZE run restarts, see $STATPT) Note that this option is unreliable! = -1 reuse dictionary file from previous run, useful with GEOM=DAF and/or GUESS=MOSAVED. Otherwise, this option is the same as 0. = 0 normal run (default) = 1 2e restart (1-e integrals and MOs saved) = 2 SCF restart (1-,2-e integrals and MOs saved) = 3 1e gradient restart = 4 2e gradient restart GEOM = select where to obtain molecular geometry = INPUT from $DATA input (default for IREST=0) = DAF read from DICTNRY file (default otherwise) As noted in the first chapter, binary file restart is not a well tested option! ========================================================== ========================================================== $SYSTEM group (optional) This group provides global control information for your computer's operation. This is system related input, and will not seem particularly chemical to you! MWORDS = the maximum replicated memory which your job can use, on every core. This is given in units of 1,000,000 words (as opposed to 1024*1024 words), where a word is defined as 64 bits. (default=1) In case finer control over the replicated memory is needed, this value can be given in units of words, with the old keyword MEMORY, instead of MWORDS. MEMDDI = the grand total memory needed for the distributed data interface (DDI) storage, given in units of 1,000,000 words. See Section 5 of this manual for an extended explanation of running with MEMDDI. note: the memory required on each processor core for a run using p cores is therefore MEMDDI/p + MWORDS. The parallel runs that currently require MEMDDI are: SCFTYP=RHF MPLEVL=2 energy or gradient SCFTYP=UHF MPLEVL=2 energy or gradient SCFTYP=ROHF MPLEVL=2 OSPT=ZAPT energy or gradient SCFTYP=MCSCF MPLEVL=2 energy SCFTYP=MCSCF using the FULLNR or JACOBI convergers SCFTYP=MCSCF analytic hessian SCFTYP=any CITYP=ALDET, ORMAS, GUGA SCFTYP=any energy localization SCFTYP=RHF CCTYP=CCSD or CCSD(T) All other parallel runs should enter MEMDDI=0, for they use only replicated memory. Some serial runs execute the parallel code (on just 1 CPU), for there is only a parallel code. These serial runs must give MEMDDI as a result: SCFTYP=ROHF MPLEVL=2 OSPT=ZAPT gradient/property run SCFTYP=MCSCF analytic hessian Two kinds of runs (RI-MP2 and parallel CCSD(T)) use an additional type of memory, for which there is no input keyword. Please read EXETYP=CHECK output carefully to learn the total memory/node requirements for these two! TIMLIM = time limit, in minutes. Set to about 95 percent of the time limit given to the batch job (if you use a queueing system) so that GAMESS can stop itself gently. (default=525600.0 minutes) PARALL = a flag to cause the program to execute the parallel algorithm, in cases where different serial and parallel codes exist, if you happen to be running on only one core. The default is .FALSE. if you are running on one core. The main purpose of this keyword is to allow you to do EXETYP=CHECK runs on only one core, when your intent is perform the actual calculation in parallel. PARALL is ignored for runs on more than one core, when of course parallel algorithms are executed. KDIAG = diagonalization control switch = 0 use a vectorized diagonalization routine if one is available on your machine, else use EVVRSP. (default) = 1 use EVVRSP diagonalization. This may be more accurate than KDIAG=0. = 2 use GIVEIS diagonalization (not as fast or reliable as EVVRSP) = 3 use JACOBI diagonalization (this is the slowest method) = 4 use DSYEV diagonalization from LAPACK = 5 use DSYEVD for DFTB Fock matrices and DSYEV elsewhere (both from LAPACK). To use LAPACK, you should link it. You may be able to edit comp and lked to take advantage of your LAPACK library. Some math libraries such as MKL, OpenBLAS, and Fujitsu SSL come with LAPACK. COREFL = a flag to indicate whether or not GAMESS should produce a "core" file for debugging when subroutine ABRT is called to kill a job. This variable pertains only to UNIX operating systems. (default=.FALSE.) BALTYP = Parallel load balance scheme: = SLB uses static load balancing. = DLB uses dynamic load balancing (default). Dynamic load balancing attempts to spread out possibly unequal work assignments based on the rate at which different nodes complete tasks. For historical reasons, it is permissible to spell SLB as LOOP, and DLB as NXTVAL. MXSEQ2 = 300 (default) MXSEQ3 = 150 (default) Matrix/vector problem size in loops requiring either O(N**2) or O(N**3) work, respectively. Problems below these sizes are run purely serial, to avoid poor communication/computation ratios. NODEXT = array specifying node extensions in GDDI for each file. Non-zero values force no extension. E.g., NODEXT(40)=1 forces file 40 (file numbers are unit numbers used in GAMESS, see "rungms" or PROG.DOC) to have the name of $JOB.F40 on all nodes, rather than $JOB.F40, $JOB.F40.001, $JOB.F40.002 etc. This option should not be used on multicore machines, and is really outdated now. (default: all zeros) IOSMP = Parallelise I/O on SMP machines with multiple hard disks. Two parameters are specified, whose meaning should be clear from the example. iosmp(1)=2,6 2 refers to the number of HDDs per SMP box. 6 is the location of the character in the file names that switches HDDs, i.e. if HDDs are mounted as /work1 and /work2, then 6 refers to the position of the number 1 in /work1. The file system should permit disks attached with directory names differing by one symbol. (default: 0,0, disable the feature) MODIO = Global I/O options (bitwise additive) (default: 0) 1 - forbid flushing files 2 - do not close dictionary file in GDDI (only record indices are reset) 4 - do not print timings on each rank at run end. 8 - forbid grid data saving in DFT (prevent F22 from being opened; however, in CPKS F22 is always opened). 16 - reduce I/O (skip various printing) 32 - do not open file F15 (in SOSCF) on workers/agents and do not close on all. 64 - use in-memory F15 (in SOSCF). This also parallelizes one more step in SOSCF. 128 - always run EVVRSP sequentially. This is useful on mixed CPU type clusters. 256 - reduce timing output. 512 - use XYZ file to store coordinates. 1024 - enforce minimal output 2048 - a single flag to turn on various accceleration options 4096 - redirect OUTPUT and PUNCH output to /dev/null on all parallel ranks except the global master. Use only if OS permits it (a UNIX). Mainly intended for GDDI runs to prevent masters from writing anything. It is assumed that the global master writes all important data. 16384 - force workers/agents to write error messages when aborting a job due to an error; mostly useful in GDDI. MEM10 = words used to store dictionary file F10 in memory. Selecting this option will replace any I/O for F10 by memory access. Default: 0 (disk-based F10) MEM22 = words used to store grid file F22 for DFT in memory. Selecting this option will skip any I/O for F22. Mainly useful for DFT Hessian and analytic FMO-DFT gradient, both using CPKS. MEM22 is not used for F22 in DFT energy/gradient (use MODIO=8 for removing I/O to F22). Default: 0 (disk-based F22) NALIGN = align dynamic memory blocks so that relative addresses are multiples of NALIGN (absolute addresses are given by OS, usually a multiple of 64 bits). Default: 1 (effectively, no allignment) ========================================================== ========================================================== $BASIS group (optional) This group allows certain standard basis sets to be easily requested. Basis sets are specified by: a) GBASIS plus optional supplementations such as NDFUNC, b) BASNAM to read custom basis sets from your input, c) EXTFIL to read custom bases from an external file, d) or omit this group entirely, and give the basis set in the $DATA input, which is completely general. GBASIS requests various Gaussian basis sets. These include options for effective core and model core potentials. Rather oddly, GBASIS also can select semi-empirical models, and in that case requests the Slater-type orbitals for the MOPAC-type calculation. Note: The first two groups of GBASIS keywords below (except G3L and G3LX) define only the basic functions, without any polarization functions and/or diffuse functions. For example, main group elements have the basic functions for their s,p valence orbitals. Polarization and/or diffuse supplements are added separately to these GBASIS values, with keywords NPFUNC, NDFUNC, NFFUNC, DIFFS, DIFFSP, POLAR, SPLIT2, and SPLIT3, which are defined at the end of this input group. GBASIS = STO - Pople's STO-NG minimal basis set. Available H-Xe, for NGAUSS=2,3,4,5,6. = N21 - Pople's N-21G split valence basis set. Available H-Xe, for NGAUSS=3. Available H-Ar, for NGAUSS=6. = N31 - Pople's N-31G split valence basis set. Available H-Ne,P-Cl for NGAUSS=4. Available H-He,C-F for NGAUSS=5. Available H-Kr, for NGAUSS=6, note that the bases for K,Ca,Ga-Kr were changed 9/2006. = N311 - Pople's "triple split" N-311G basis set. Available H-Ne, for NGAUSS=6. Selecting N311 implies MC for Na-Ar. = G3L - Pople's G3MP2Large basis set, for H-Kr. = G3LX - Pople's G3MP2LargeXP basis set, for H-Kr. NGAUSS = the number of Gaussians (N). This parameter pertains to GBASIS=STO, N21, N31, or N311. GBASIS = MINI - Huzinaga's 3 gaussian minimal basis set. Available H-Rn. = MIDI - Huzinaga's 21 split valence basis set. Available H-Rn. = DZV - "double zeta valence" basis set. a synonym for DH for H,Li,Be-Ne,Al-Cl. (14s,9p,3d)/[5s,3p,1d] for K-Ca. (14s,11p,5d/[6s,4p,1d] for Ga-Kr. = DH - Dunning/Hay "double zeta" basis set. (3s)/[2s] for H. (9s,4p)/[3s,2p] for Li. (9s,5p)/[3s,2p] for Be-Ne. (11s,7p)/[6s,4p] for Al-Cl. = TZV - "triple zeta valence" basis set. (5s)/[3s] for H. (10s,3p)/[4s,3p] for Li. (10s,6p)/[5s,3p] for Be-Ne. a synonym for MC for Na-Ar. (14s,9p)/[8s,4p] for K-Ca. (14s,11p,6d)/[10s,8p,3d] for Sc-Zn. = MC - McLean/Chandler "triple split" basis. (12s,9p)/[6s,5p] for Na-Ar. Selecting MC implies 6-311G for H-Ne. = MINIX The basis set for HF-3C. GBASIS=MINIX will not add three energy corrections (to add them, use GBASIS=HF-3C). Supported elements are H-Kr, and MINIX means H-Ne, B-Ne: MINIS (a minimum basis set) Li-Be, Na-Mg, Al-Ar: MINIS+p (plus polarization) K-Zn: SV (a minimum basis set) Ga-Kr: SVP (SV + polarization) = HF-3C The minimal basis set MINIX and three add-on energy corrections (3c). These three corrections are D3(BJ), GCP and SRB. See $DFT's dispersion corrections. * * systematic basis set families * * * These four families provide a hierachy of basis sets approaching the complete basis set limits. These families include relevant polarization and diffuse augmentations, as indicated in their names. GBASIS = CCn - Dunning-type Correlation Consistent basis sets, officially called cc-pVnZ. Use n = D,T,Q,5,6 to indicate the level of polarization. Available for H-He, Li-Ne, Na-Ar, Ca, Ga-Kr and for Sc-Zn for n=T,Q. = ACCn - As CCn, but augmented with a set of diffuse functions, e.g. aug-cc-pVnZ. Availability is the same as CCn. = CCnC - As CCn, but augmented with tight functions for recovering core and core-valence correlation, e.g. cc-pCVnZ. Available H-Ar for n=D,T,Q, also n=5 for H-Ne. = ACCnC- As CCn, augmented with diffuse as well as CCnC's tight functions, e.g. aug-cc-pCVnZ. Availability is the same as CCnC. = CCnWC the omega form of CCnC, e.g. cc-pwCVnZ, for H-Ar, for n=T only. CCnWC's tight functions are considered superior to CCnC's for recovery of core/valence correlation. = ACCnWC augmented form of CCnWC: aug-cc-pwCVnZ. See extended notes below! = PCseg-n - Polarization Consistent basis sets. n = 0,1,2,3,4 indicates the level of polarization. (n=0 is unpolarized, n=1 is ~DZP, n=2 is ~TZ2P, etc.). These provide a hierarchy of basis sets suitable for DFT and HF calculations. Available for H-Kr. = APCseg-n - These are the PCseg-n bases, with diffuse augmentation. See extended notes below! Sapporo valence basis sets: = SPK-nZP - Sapporo family of non-relativistic bases, n=D,T,Q, available H-Xe = SPK-AnZP - diffuse augmentation of the above. = SPKrnZP - Sapporo family of relativistic bases n=D,T,Q, available H-Xe. These should be used only with a relativistic transformation of the integrals, such as RELWFN=LUT-IOTC. = SPKrAnZP - diffuse augmentation of the above. See extended notes below! Sapporo core/valence basis sets: = SPK-nZC - Sapporo family of non-relativistic bases, n=D,T,Q, available H-Xe = SPK-nZCD - diffuse augmentation of the above. = SPKrnZC - Sapporo family of relativistic bases n=D,T,Q, available H-Rn. To be used only with a relativistic transformation of the integrals, such as RELWFN=LUT-IOTC. = SPKrnZCD - diffuse augmentation of the above. See extended notes below! = KTZV - Karlsruhe valence triple zeta basis, as developed by Prof.Ahlrichs, see REFS.DOC. = KTZVP- Karlsruhe valence triple zeta basis with a set of single polarization (P). = KTZVPP-Karlsruhe valence triple zeta basis with a set of double polarization (PP). The Karlsruhe sets are provided for H-Ar. Normally these families are used as spherical harmonics, see ISPHER=1 in $CONTRL. Failure to set ISPHER=1 will result in discrepancies in energy values compared to the literature or other programs, difficulties in converging SCF/DFT, CC, CI, and/or response equation iterations, and longer run times due to retention of unimportant MOs. The calculations will refuse to run without ISPHER being set. Important note about the PCseg basis set family: 1. These should be used only in spherical harmonic form. 2. The PCn basis sets included in GAMESS versions prior to March 2014 were generally contracted, but were replaced by computationally more efficient segmented contractions, and renamed to PCseg-n. The segemented contractions have the same or slightly better accuracy (especially for n=0) as the original PCn bases, which are no longer available. Important notes about the CC basis set family: 1. These should be used only in spherical harmonic form. 2. The CC5 and CC6 basis sets (and corresponding augmented versions) contain h-functions, and CC6 also contains i- functions. As of January 2013, GAMESS integral code can correctly use h & i functions, so these three call up the true basis sets. Prior to January 2013, GAMESS' integral codes was restricted to g-functions, so these three truncated away any h & i functions, to spdfg subsets, and therefore were not the true basis sets. 3. Note that the CC basis sets are generally contracted, which GAMESS can only handle by replicating the primitive basis functions, leading to a less than optimum performance in AO integral evaluation. 4. The implementation of the cc-pVnZ and cc-pCVnZ basis sets for Na-Ar include one additional tight d-function, producing the so-called cc-pV(n+d)Z and cc-pV(n+d)Z sets, which are known to improve results (see J.Chem.Phys. 114, 9244(2001) and Theoret.Chem.Acc. 120, 119(2008)). These tight d versions are invoked by GBASIS=CCn or CCnC (and also their augmented counterparts ACCn or ACC). This means the old (and less accurate) basis sets without the tight d's are not available for Na-Ar. 5. Alkali and alkali earth basis sets (Li,Be,Na,Mg) were changed April 2013 so that regular, diffuse, tight d (for Na/Mg), core/valence, and weighted core/valence sets agree with their official publication: Theor. Chem. Acc. 128, 69(2011). 6. In case you are interested in scalar relativistic effects, the CCT-DK and CCQ-DK sets optimized for use with Douglas/Kroll are available for Sc-Kr. These will be used if you type GBASIS=CCT or CCQ along with RELWFN=RESC, DK, IOTC, IOTC2E, LUT-IOTC, or LUTIOTC2, while using NR sets for elements lighter than Sc. DK versions of ACCD or ACCT are available for Sc- Zn (but not for the rest of the row, Ga-Kr). Notes about the Sapporo basis set family: 1. SPK is the international airport city code for Sapporo. 2. These should be used only in spherical harmonic form. 3. The relativistic core/valence sets are available for all atoms including the 6th row of the periodic table (H-Rn). 4. It is extremely illogical to use any of the all-electron relativistic bases without turning on scalar relativity! So, choose RELWFN=LUT-IOTC (or IOTC, IOTC2E, LUTIOTC2, DK, RESC) in $CONTRL. 5. The core/valence basis sets treat (n-1)s,(n-1)p,ns for s-block elements; (n-1)s,(n-1)p,ns,np for p block elements; (n-1)s,(n-1)p,(n-1)d,ns for d block elements, and the 4s- 4f,5s-5d,6s for f block (lanthanides). This suggests you should change the default number of core orbitals, such as NACORE in $MP2 or NCORE in $CCINP, to correlate the indicated semi-core orbitals (this is not automatic). 6. The relativistic sets ("r") are identical to the non- relativistic choices ("-") for atoms H-Ar, where scalar relativity has almost no effect on orbital shapes. 7. The relativistic bases were optimized at the 3rd order of the Douglas-Kroll transformation, with a Gaussian nuclei model. It should be fine to use them with any RESC, DK, IOTC, IOTC2E, LUT-IOTC, or LUTIOTC2 calculation. 8. Because they are stored in an external file supplied with GAMESS, these can only be accessed via GBASIS in this group, not by using them in-line in $DATA. 9. The SPK basis sets were extracted from the data base of Segmented Gaussian Basis Sets, maintained by Takeshi Noro, University of Hokkaido, Sapporo, Japan: http://setani.sci.hokudai.ac.jp/sapporo/Welcome.do The mapping between the data base names and the keywords used in GAMESS is (for n=D,T,Q): data base name keyword Sapporo-nZP SPK-nZP Sapporo-nZP+diffuse SPK-AnZP Sapporo-DK-nZP SPKrnZP Sapporo-DK-nZP+diffuse SPKrAnZP Sapporo-nZP-2012 SPK-nZC Sapporo-nZP-2012+diffuse SPK-nZCD Sapporo-DK-nZP-2012 SPKrnZC Sapporo-DK-nZP-2012+diffuse SPKrnZCD * * * Effective Core Potential (ECP) bases * * * GBASIS = SBKJC- Stevens/Basch/Krauss/Jasien/Cundari valence basis set, for Li-Rn. This choice implies an unscaled -31G basis for H-He. = HW - Hay/Wadt valence basis. This is a -21 split, available Na-Xe, except for the transition metals. This implies a 3-21G basis for H-Ne. * * * Model Core Potential (MCP) bases * * * Notes: Select PP=MCP in $CONTRL to automatically use the model core potential matching your basis choice below. References for these bases, and other information about MCPs can be found in the REFS.DOC chapter. Another family covering almost all elements is available in $DATA only. GBASIS = MCP-DZP, MCP-TZP, MCP-QZP - a family of double, triple, and quadruple zeta quality valence basis sets, which are akin to the correlation consistent sets, in that these include increasing levels of polarization (and so do not require "supplements" like NDFUNC or DIFFSP) and must be used as spherical harmonics (see ISPHER). Availability: MCP-DZP: 56 elements Z=3-88, except V-Zn, Y-Cd, La, Hf-Hg MCP-TZP, MCP-QZP: 85 elements Z=3-88, except La The basis sets for hydrogen atoms will be the corresponding Dunning's cc-pVNZ (N=D,T,Q). = MCP-ATZP, MCP-AQZP - MCP-TZP and MCP-QZP core potentials whose basis sets were augmented with diffuse functions Availability: same as for MCP-TZP, MCP-QZP = MCPCDZP, MCPCTZP, MCPCQZP - based on MCP-DZP, MCP-TZP, MCP-QZP, with core-valence functions provided for the alkali and alkaline earth atoms Na through Ra. = MCPACDZP, MCPACTZP, MCPACQZP - based on MCPCDZP, MCPCTZP, MCPCQZP, with core-valence functions provided for the alkali and alkaline earth atoms Na through Ra, and augmented with diffuse functions. The basis sets were extracted from the data base Segmented Gaussian Basis Sets, maintained by Takeshi Noro, Quantum Chemistry Group, Sapporo, Japan: http://setani.sci.hokudai.ac.jp/sapporo/Welcome.do The mapping between the data base names and the names used in GAMESS is data base name GAMESS keyword MCP/NOSeC-V-DZP MCP-DZP MCP/NOSeC-V-TZP MCP-TZP MCP/NOSeC-V-QZP MCP-QZP MCP/NOSeC-V-TZP+diffuse MCP-ATZP MCP/NOSeC-V-QZP+diffuse MCP-AQZP MCP/NOSeC-CV-DZP MCPCDZP MCP/NOSeC-CV-TZP MCPCTZP MCP/NOSeC-CV-QZP MCPCQZP MCP/NOSeC-CV-DZP+diffuse MCPACDZP MCP/NOSeC-CV-TZP+diffuse MCPACTZP MCP/NOSeC-CV-QZP+diffuse MCPACQZP GBASIS = IMCP-SR1 and IMCP-SR2 - valence basis sets to be used with the improved MCPs with scalar relativistic effects. These are available for transition metals except La, and the main group elements B-Ne, P-Ar, Ge, Kr, Sb, Xe, Rn. The 1 and 2 refer to addition of first and second polarization shells, so again don't use any of the "supplements" and do use spherical harmonics. = IMCP-NR1 and IMCP-NR2 - closely related valence basis sets, but with nonrelativistic model core potentials. GBASIS = ZFK3-DK3, ZFK4-DK3, ZFK5-DK3, or ZFK3LDK3, ZFK4LDK3, ZFK5LDK3 These are a family of model core potential basis sets developed by Zeng/Fedorov/Klobukowski, for the p-block elements from 2p to 6p. The potentials were paramaterized taking into account both DK3 scalar relativistic and DK-SOC effects. The fundamental basis functions are from the Well-Tempered Basis Sets. The number after ZFK indicates the augmentation levels, e.g. ZFK3 means the diffuse functions from aug-cc-pVTZ are added, ZFK4 means from aug- cc-pVQZ, etc. The difference between ZFKn-DK3 and ZFKnLDK3 is that the common s and p exponents have been contracted as a single L-shell for the outermost s and p valence shells to save time in the "L" case. The s-block elements from 1s to 4s have also been put in the library. For H/He, all-electron aug-cc-pVnZ basis sets are used. For Li/Be, the relativistically contracted atomic natural orbital all- electron basis sets (ANO-RCC) are used. For Na/Mg, and K/Ca, unpublished MCP and basis sets based on ANO-RCC are available, although the potentials have not been extensively tested yet. No d-block elements can be used. * * * semiempirical basis sets * * * GBASIS = MNDO - selects MNDO model Hamiltonian = AM1 - selects AM1 model Hamiltonian = PM3 - selects PM3 model Hamiltonian = RM1 - selects RM1 model Hamiltonian = DFTB - selects tight binding Hamiltonian Note: The elements for which these exist can be found in the 'further information' section of this manual. If you pick one of these, all other data in this group is ignored. Semi-empirical runs actually use valence-only Slater type orbitals (STOs), not Gaussian GTOs, but the keyword remains GBASIS. Except for NGAUSS, all other keywords such as NDFUNC, etc. will be ignored for these. If you add NGAUSS, STO-NG expansions of the valence STO functions in terms of Gaussians will be added to the log file. Plotting programs such as MacMolPlt can pick up this approximation to the STOs used up from the ouput, in order to draw the orbitals. The default NGAUSS=0 suppresses this output, but values up to 6 may be given to control the accuracy of the STO-NG printing. --- supplementary functions --- NDFUNC = number of heavy atom polarization functions to be used. These are usually d functions, except for MINI/MIDI. The term "heavy" means Na on up when GBASIS=STO, HW, or N21, and from Li on up otherwise. The value may not exceed 3. The variable POLAR selects the actual exponents to be used, see also SPLIT2 and SPLIT3. (default=0) NFFUNC = number of heavy atom f type polarization functions to be used on Li-Cl. This may only be input as 0 or 1. (default=0) NPFUNC = number of light atom, p type polarization functions to be used on H-He. This may not exceed 3, see also POLAR. (default=0) DIFFSP = flag to add diffuse sp (L) shell to heavy atoms. Heavy means Li-F, Na-Cl, Ga-Br, In-I, Tl-At. The default is .FALSE. DIFFS = flag to add diffuse s shell to hydrogens. The default is .FALSE. Warning: if you use diffuse functions, please read QMTTOL in the $CONTRL input group for numerical concerns. POLAR = exponent of polarization functions = COMMON (default for GBASIS=STO,N21,HW,SBKJC) = POPN31 (default for GBASIS=N31) = POPN311 (default for GBASIS=N311, MC) = DUNNING (default for GBASIS=DH, DZV) = HUZINAGA (default for GBASIS=MINI, MIDI) = HONDO7 (default for GBASIS=TZV) SPLIT2 = an array of splitting factors used when NDFUNC or NPFUNC is 2. Default=2.0,0.5 SPLIT3 = an array of splitting factors used when NDFUNC or NPFUNC is 3. Default=4.00,1.00,0.25 The splitting factors are from the Pople school, and are probably too far apart. See for example the Binning and Curtiss paper. For example, the SPLIT2 value will usually cause an INCREASE over the 1d energy at the HF level for hydrocarbons. The actual exponents used for polarization functions, as well as for diffuse sp or s shells, are described in the 'Further References' section of this manual. This section also describes the sp part of the basis set chosen by GBASIS fully, with all references cited. Note that GAMESS always punches a full $DATA input group. Thus, if $BASIS does not quite cover the basis you want, you can obtain this full $DATA from EXETYP=CHECK, and then change polarization exponents, add Rydbergs, etc. * * * This may only be used with COORD=UNIQUE or HINT! BASNAM = an array of names of customized basis set input groups. BASNAM should obey the rule of no more than six characters starting with a letter names, and must avoid using any GBASIS string. However, the individual basis inputs can use any of the GBASIS sets by its standard name. Basis supplementations such as DIFFS or NDFUNC may only be given by explicit numerical values. This is best explained by an example where a core potential and valence-only basis set is used on a transition metal, but not its ligands: $contrl scftyp=rohf icharg=+3 mult=4 runtyp=gradient pp=read ispher=1 $end $system mwords=1 $end $guess guess=huckel $end $basis basnam(1)=metal, ligO,ligO,ligO,ligO,ligO,ligO, ligH,ligH,ligH,ligH,ligH,ligH, ligH,ligH,ligH,ligH,ligH,ligH $end $data Cr+3(H2O)6 complex...SBKJC & 6-31G(d) geometry Th CHROMIUM 24.0 .0000000000 .0 .0000000000 OXYGEN 8.0 .0000000000 .0 2.0398916104 HYDROGEN 1.0 .7757887450 .0 2.6122732372 $end ! core potential basis for Chromium $metal sbkjc $end ! normal 6-31G(d) for oxygen ligands $ligO n31 6 d 1 ; 1 0.8 1.0 $end ! unpolarized basis for hydrogens $ligH n31 6 $end $ecp Cr-ecp SBKJC O-ecp none ...snipped... there must be 6 O's given here O-ecp none H-ecp none ...snipped... there must be 12 H's given here H-ecp none $end * * * EXTFIL = a flag to read basis sets from an external file, defined by EXTBAS, rather than from $DATA. (default=.false.) It may be easier to use BASNAM to create custom basis sets! BASNAM has the bonus that your input file contains all information about the calculation, explicitly. Except for MCP basis sets, no external file is provided with GAMESS, thus you must create your own. The GBASIS keyword must give an 8 or less character string, obviously not using any internally stored names. Every atom must be defined in the external file by a line giving the chemical symbol, and this chosen string. Following this header line, give the basis in free format $DATA style, containing only S, P, D, F, G, and L shells, and terminating each atom by the usual blank line. The external file may have several families of bases in the same file, identified by different GBASIS strings. =========================================================== ========================================================== $DATA group (required) $DATAS group (if NESC chosen, for small component basis) $DATAL group (if NESC chosen, for large component basis) This group describes the global molecular data such as point group symmetry, nuclear coordinates, and possibly the basis set. It consists of a series of free format card images. See $RELWFN for more information on large and small component basis sets. The input structure of $DATAS and $DATAL is identical to the COORD=UNIQUE $DATA input. ---------------------------------------------------------- -1- TITLE a single descriptive title card. ---------------------------------------------------------- -2- GROUP, NAXIS GROUP is the Schoenflies symbol of the symmetry group, you may choose from C1, Cs, Ci, Cn, S2n, Cnh, Cnv, Dn, Dnh, Dnd, T, Th, Td, O, Oh. NAXIS is the order of the highest rotation axis, and must be given when the name of the group contains an N. For example, "Cnv 2" is C2v. "S2n 3" means S6. Use of NAXIS up to 8 is supported in each axial groups. For linear molecules, choose either Cnv or Dnh, and enter NAXIS as 4. Enter atoms as Dnh with NAXIS=2. If the electronic state of either is degenerate, check the note about the effect of symmetry in the electronic state in the SCF section of REFS.DOC. ---------------------------------------------------------- In order to use GAMESS effectively, you must be able to recognize the point group name for your molecule. This presupposes a knowledge of group theory at about the level of Cotton's "Group Theory", Section 3. Armed with only the name of the group, GAMESS is able to exploit the molecular symmetry throughout almost all of the program, and thus save a great deal of computer time. GAMESS does not require that you know very much else about group theory, although a deeper knowledge (character tables, irreducible representations, term symbols, and so on) is useful when dealing with the more sophisticated wavefunctions. Cards -3- and -4- are quite complicated, and are rarely given. A *SINGLE* blank card may replace both cards -3- and -4-, to select the 'master frame', which is defined on the next page. If you choose to enter a blank line, skip to one of the -5- input sequences. Note! If the point group is C1 (no symmetry), skip over cards -3- and -4- (which means no blank card). ---------------------------------------------------------- -3- X1, Y1, Z1, X2, Y2, Z2 For C1 group, there is no card -3- or -4-. For CI group, give one point, the center of inversion. For CS group, any two points in the symmetry plane. For axial groups, any two points on the principal axis. For tetrahedral groups, any two points on a two-fold axis. For octahedral groups, any two points on a four-fold axis. ---------------------------------------------------------- -4- X3, Y3, Z3, DIRECT third point, and a directional parameter. For CS group, one point of the symmetry plane, noncollinear with points 1 and 2. For CI group, there is no card -4-. For other groups, a generator sigma-v plane (if any) is the (x,z) plane of the local frame (CNV point groups). A generator sigma-h plane (if any) is the (x,y) plane of the local frame (CNH and dihedral groups). A generator C2 axis (if any) is the x-axis of the local frame (dihedral groups). The perpendicular to the principal axis passing through the third point defines a direction called D1. If DIRECT='PARALLEL', the x-axis of the local frame coincides with the direction D1. If DIRECT='NORMAL', the x-axis of the local frame is the common perpendicular to D1 and the principal axis, passing through the intersection point of these two lines. Thus D1 coincides in this case with the negative y axis. ---------------------------------------------------------- The 'master frame' is just a standard orientation for the molecule. By default, the 'master frame' assumes that 1. z is the principal rotation axis (if any), 2. x is a perpendicular two-fold axis (if any), 3. xz is the sigma-v plane (if any), and 4. xy is the sigma-h plane (if any). Use the lowest number rule that applies to your molecule. Some examples of these rules: Ammonia (C3v): the unique H lies in the XZ plane (R1,R3). Ethane (D3d): the unique H lies in the YZ plane (R1,R2). Methane (Td): the H lies in the XYZ direction (R2). Since there is more than one 3-fold, R1 does not apply. HP=O (Cs): the mirror plane is the XY plane (R4). In general, it is a poor idea to try to reorient the molecule. Certain sections of the program, such as the orbital symmetry assignment, do not know how to deal with cases where the 'master frame' has been changed. Linear molecules (C4v or D4h) must lie along the z axis, so do not try to reorient linear molecules. You can use EXETYP=CHECK to quickly find what atoms are generated, and in what order. This is typically necessary in order to use the general $ZMAT coordinates. * * * * Depending on your choice for COORD in $CONTROL, if COORD=UNIQUE, follow card sequence U if COORD=HINT, follow card sequence U if COORD=CART, follow card sequence C if COORD=ZMT, follow card sequence G if COORD=ZMTMPC, follow card sequence M Card sequence U is the only one which allows you to define a completely general basis here in $DATA. Recall that UNIT in $CONTRL determines the distance units. ---------------------------------------------------------- -5U- Atom input. Only the symmetry unique atoms are input, GAMESS will generate the symmetry equivalent atoms according to the point group selected above. if COORD=UNIQUE NAME, ZNUC, X, Y, Z *************** NAME = 10 character atomic name, used only for printout. Thus you can enter H or Hydrogen, or whatever. ZNUC = nuclear charge. It is the nuclear charge which actually defines the atom's identity. X,Y,Z = Cartesian coordinates. if COORD=HINT ************* NAME,ZNUC,CONX,R,ALPHA,BETA,SIGN,POINT1,POINT2,POINT3 NAME = 10 character atomic name (used only for print out). ZNUC = nuclear charge. CONX = connection type, choose from 'LC' linear conn. 'CCPA' central conn. 'PCC' planar central conn. with polar atom 'NPCC' non-planar central conn. 'TCT' terminal conn. 'PTC' planar terminal conn. with torsion R = connection distance. ALPHA= first connection angle BETA = second connection angle SIGN = connection sign, '+' or '-' POINT1, POINT2, POINT3 = connection points, a serial number of a previously input atom, or one of 4 standard points: O,I,J,K (origin and unit points on axes of master frame). defaults: POINT1='O', POINT2='I', POINT3='J' ref- R.L. Hilderbrandt, J.Chem.Phys. 51, 1654 (1969). You cannot understand HINT input without reading this. Note that if ZNUC is negative, the internally stored basis for ABS(ZNUC) is placed on this center, but the calculation uses ZNUC=0 after this. This is useful for basis set superposition error (BSSE) calculations. ---------------------------------------------------------- * * * If you gave $BASIS, continue entering cards -5U- until all the unique atoms have been specified. When you are done, enter a " $END " card. * * * If you did not, enter cards -6U-, -7U-, -8U-. ---------------------------------------------------------- -6U- GBASIS, NGAUSS, (SCALF(i),i=1,4) GBASIS can have exactly the same meaning as the keyword in $BASIS. You may choose from STO, N21, N31, N311, ACCT, PC4, ... A few of these require NGAUSS below. In addition, GBASIS can be S, P (or L), D, F, G, H, or I to enter an explicit basis set. L means both an S and P shell with the same exponent. See NGAUSS below. In addition, GBASIS may be defined as MCP, to indicate that the current atom is represented by a model core potential, and valence basis set. An internally stored basis and potential will be applied (see REFS.DOC for the details). The MCP basis supplies only the occupied atomic orbitals, e.g. sp for a main group element, so please supplement with any desired polarization. In case the keyword MCP is followed by the keyword READ, everything will be taken from the input file, namely the basis functions are read using the sequence -6U-, -7U-, and -8U-, from lines following the "MCP READ" line. In addition, "MCP READ" implies that the parameters of the model core potentials, together with core basis functions are in the input stream, in a $MCP input group. Other MCP bases are available in the $BASIS input, but note that to locate the MCP, the atom name must be a chemical symbol, that is "P" instead of "Phosphorus". NGAUSS is the number of Gaussians (N) in the Pople style basis, or user input general basis. It has meaning only for GBASIS=STO, N21, N31, or N311, or explicit GTO types such as S,P,D,F... Up to 4 scale factors may be entered. If omitted, standard values are used. They are not documented as every GBASIS treats these differently. Read the source code if you need to know more. They are seldom given. ---------------------------------------------------------- * * * If GBASIS is not S,P,D,F,... either add more shells by repeating card -6U-, or go on to -8U-. * * * If GBASIS=S,P,D,F,... enter NGAUSS cards -7U-. ---------------------------------------------------------- -7U- IG, ZETA, C1, C2 IG = a counter, IG takes values 1, 2, ..., NGAUSS. ZETA = Gaussian exponent of the IG'th primitive. C1 = Contraction coefficient for S,P,D,F,G shells, and for the s function of L shells. C2 = Contraction coefficient for the p in L shells. ---------------------------------------------------------- * * * For more shells on this atom, go back to card -6U-. * * * If there are no more shells, go on to card -8U-. ---------------------------------------------------------- -8U- A blank card ends the basis set for this atom. ---------------------------------------------------------- Continue entering atoms with -5U- through -8U- until all are given, then terminate the group with a " $END " card. --- this is the end of card sequence U --- COORD=CART input: ---------------------------------------------------------- -5C- Atom input. Cartesian coordinates for all atoms must be entered. They may be arbitrarily rotated or translated, but must possess the actual point group symmetry. GAMESS will reorient the molecule into the 'master frame', and determine which atoms are the unique ones. Thus, the final order of the atoms may be different from what you enter here. NAME, ZNUC, X, Y, Z NAME = 10 character atomic name, used only for printout. Thus you can enter H or Hydrogen, or whatever. ZNUC = nuclear charge. It is the nuclear charge which actually defines the atom's identity. X,Y,Z = Cartesian coordinates. ---------------------------------------------------------- Continue entering atoms with card -5C- until all are given, and then terminate the group with a " $END " card. --- this is the end of card sequence C --- COORD=ZMT input: (GAUSSIAN style internals) ---------------------------------------------------------- -5G- ATOM Only the name of the first atom is required. See -8G- for a description of this information. ---------------------------------------------------------- -6G- ATOM i1 BLENGTH Only a name and a bond distance is required for atom 2. See -8G- for a description of this information. ---------------------------------------------------------- -7G- ATOM i1 BLENGTH i2 ALPHA Only a name, distance, and angle are required for atom 3. See -8G- for a description of this information. ---------------------------------------------------------- -8G- ATOM i1 BLENGTH i2 ALPHA i3 BETA i4 ATOM is the chemical symbol of this atom. It can be followed by numbers, if desired, for example Si3. The chemical symbol implies the nuclear charge. i1 defines the connectivity of the following bond. BLENGTH is the bond length "this atom-atom i1". i2 defines the connectivity of the following angle. ALPHA is the angle "this atom-atom i1-atom i2". i3 defines the connectivity of the following angle. BETA is either the dihedral angle "this atom-atom i1- atom i2-atom i3", or perhaps a second bond angle "this atom-atom i1-atom i3". i4 defines the nature of BETA, If BETA is a dihedral angle, i4=0 (default). If BETA is a second bond angle, i4=+/-1. (sign specifies one of two possible directions). ---------------------------------------------------------- o Repeat -8G- for atoms 4, 5, ... o The use of ghost atoms is possible, by using X or BQ for the chemical symbol. Ghost atoms preclude the option of an automatic generation of $ZMAT. o The connectivity i1, i2, i3 may be given as integers, 1, 2, 3, 4, 5,... or as strings which match one of the ATOMs. In this case, numbers must be added to the ATOM strings to ensure uniqueness! o In -6G- to -8G-, symbolic strings may be given in place of numeric values for BLENGTH, ALPHA, and BETA. The same string may be repeated, which is handy in enforcing symmetry. If the string is preceded by a minus sign, the numeric value which will be used is the opposite, of course. Any mixture of numeric data and symbols may be given. If any strings were given in -6G- to -8G-, you must provide cards -9G- and -10G-, otherwise you may terminate the group now with a " $END " card. ---------------------------------------------------------- -9G- A blank line terminates the Z-matrix section. ---------------------------------------------------------- -10G- STRING VALUE STRING is a symbolic string used in the Z-matrix. VALUE is the numeric value to substitute for that string. ---------------------------------------------------------- Continue entering -10G- until all STRINGs are defined. Note that any blank card encountered while reading -10G- will be ignored. GAMESS regards all STRINGs as variables (constraints are sometimes applied in $STATPT). It is not necessary to place constraints to preserve point group symmetry, as GAMESS will never lower the symmetry from that given at -2-. When you have given all STRINGs a VALUE, terminate the group with a " $END " card. --- this is the end of card sequence G --- * * * * The documentation for sequence G above and sequence M below presumes you are reasonably familiar with the input to GAUSSIAN or MOPAC. It is probably too terse to be understood very well if you are unfamiliar with these. A good tutorial on both styles of Z-matrix input can be found in Tim Clark's book "A Handbook of Computational Chemistry", published by John Wiley & Sons, 1985. Both Z-matrix input styles must generate a molecule which possesses the symmetry you requested at -2-. If not, your job will be terminated automatically. COORD=ZMTMPC input: (MOPAC style internals) ---------------------------------------------------------- -5M- ATOM Only the name of the first atom is required. See -8M- for a description of this information. ---------------------------------------------------------- -6M- ATOM BLENGTH Only a name and a bond distance is required for atom 2. See -8M- for a description of this information. ---------------------------------------------------------- -7M- ATOM BLENGTH j1 ALPHA j2 Only a bond distance from atom 2, and an angle with respect to atom 1 is required for atom 3. If you prefer to hook atom 3 to atom 1, you must give connectivity as in -8M-. See -8M- for a description of this information. ---------------------------------------------------------- -8M- ATOM BLENGTH j1 ALPHA j2 BETA j3 i1 i2 i3 ATOM, BLENGTH, ALPHA, BETA, i1, i2 and i3 are as described at -8G-. However, BLENGTH, ALPHA, and BETA must be given as numerical values only. In addition, BETA is always a dihedral angle. i1, i2, i3 must be integers only. The j1, j2 and j3 integers, used in MOPAC to signal optimization of parameters, must be supplied but are ignored here. You may give them as 0, for example. ---------------------------------------------------------- Continue entering atoms 3, 4, 5, ... with -8M- cards until all are given, and then terminate the group by giving a " $END " card. --- this is the end of card sequence M --- ========================================================== This is the end of $DATA! If you have any doubt about what molecule and basis set you are defining, or what order the atoms will be generated in, simply execute an EXETYP=CHECK job to find out! ========================================================== $ZMAT group (required if NZVAR is nonzero in $CONTRL) This group lets you define the internal coordinates in which the gradient geometry search is carried out. These need not be the same as the internal coordinates used in $DATA. The coordinates may be simple Z-matrix types, delocalized coordinates, or natural internal coordinates. You must input a total of M=3N-6 internal coordinates (M=3N-5 for linear molecules). NZVAR in $CONTRL can be less than M IF AND ONLY IF you are using linear bends. It is also possible to input more than M coordinates if they are used to form exactly M linear combinations for new internals. These may be symmetry coordinates or natural internal coordinates. If NZVAR > M, you must input IJS and SIJ below to form M new coordinates. See DECOMP in $FORCE for the only circumstance in which you may enter a larger NZVAR without giving SIJ and IJS. **** IZMAT defines simple internal coordinates **** IZMAT is an array of integers defining each coordinate. The general form for each internal coordinate is code number,I,J,K,L,M,N IZMAT =1 followed by two atom numbers. (I-J bond length) =2 followed by three numbers. (I-J-K bond angle) =3 followed by four numbers. (dihedral angle) Torsion angle between planes I-J-K and J-K-L. =4 followed by four atom numbers. (atom-plane) Out-of-plane angle from bond I-J to plane J-K-L. =5 followed by three numbers. (I-J-K linear bend) Counts as 2 coordinates for the degenerate bend, normally J is the center atom. See $LIBE. =6 followed by five atom numbers. (dihedral angle) Dihedral angle between planes I-J-K and K-L-M. =7 followed by six atom numbers. (ghost torsion) Let A be the midpoint between atoms I and J, and B be the midpoint between atoms M and N. This coordinate is the dihedral angle A-K-L-B. The atoms I,J and/or M,N may be the same atom number. (If I=J AND M=N, this is a conventional torsion). Examples: N2H4, or, with one common pair, H2POH. Example - a nonlinear triatomic, atom 2 in the middle: $ZMAT IZMAT(1)=1,1,2, 2,1,2,3, 1,2,3 $END This sets up two bonds and the angle between them. The blanks between each coordinate definition are not necessary, but improve readability mightily. **** the next define delocalized coordinates **** DLC is a flag to request delocalized coordinates. (default is .FALSE.) AUTO is a flag to generate all redundant coordinates, automatically. The DLC space will consist of all non-redundant combinations of these which can be found. The list of redundant coordinates will consist of bonds, angles, and torsions only. (default is .FALSE.) NONVDW is an array of atom pairs which are to be joined by a bond, but might be skipped by the routine that automatically includes all distances shorter than the sum of van der Waals radii. Any angles and torsions associated with the new bond(s) are also automatically included. Cases where the AUTO generation of DLC coordinates fails to find the full set of 3N-6 coordinates typically fall 6 short of 3N-6. These cases are invariably due to the system being divided into pieces too far apart to have bonds detected (for example, system A might be H-bonded to system B, finding 3N-12 coordinates only). Adding NONVDW input for that H-bond will tie A and B together, and result in a correct AUTO generation of all 3N-6 coordinates. Falling short by an integer multiple of 6 indicates more than two pieces, requiring several NONVDW pairs. Falling short by 3 coordinates indicates one of the separate systems A or B is likely a single atom, which has no rotational degrees of freedom, again it should be attached by NONVDW. DLC coordinate generation can be fine tuned by IXZMAT, IRZMAT, IFZMAT whose format is the same as IZMAT: IXZMAT is an extra array of simple internal coordinates which you want to have added to the list generated by AUTO. Unlike NONVDW, IXZMAT will add only the coordinate(s) you specify. IRZMAT is an array of simple internal coordinates which you would like to remove from the AUTO list of redundant coordinates. It is sometimes necessary to remove a torsion if other torsions around a bond are being frozen, to obtain a nonsingular G matrix. IFZMAT is an array of simple internal coordinates which you would like to freeze. See also FVALUE below, which is --required-- input when IFZMAT is given. IFZMAT/FVALUE work with ordinary coordinate input using IZMAT, as well as with DLC, but in the former case be careful that IFZMAT specifies coordinates that were already given in IZMAT. In addition, IFZMAT works only for IZMAT=1,2,3 type coordinates. See IFREEZ in $STATPT you wish to freeze regular or natural internal coordinates. FVALUE is an array of values to which the internal coordinates should be constrained. It is not necessary to input $DATA such that the initial values match these desired final values, but it is helpful if the initial values are not too far away. **** SIJ,IJS define natural internal coordinates **** SIJ is a transformation matrix of dimension NZVAR x M, used to transform the NZVAR internal coordinates in IZMAT into M new internal coordinates. SIJ is a sparse matrix, so only the non-zero elements are given, by using the IJS array described below. The columns of SIJ will be normalized by GAMESS. (Default: SIJ = I, unit matrix) IJS is an array of pairs of indices, giving the row and column index of the entries in SIJ. example - if the above triatomic is water, using IJS(1) = 1,1, 3,1, 1,2, 3,2, 2,3 SIJ(1) = 1.0, 1.0, 1.0,-1.0, 1.0 gives the matrix S= 1.0 1.0 0.0 0.0 0.0 1.0 1.0 -1.0 0.0 which defines the symmetric stretch, asymmetric stretch, and bend of water. references for natural internal coordinates: P.Pulay, G.Fogarasi, F.Pang, J.E.Boggs J.Am.Chem.Soc. 101, 2550-2560(1979) G.Fogarasi, X.Zhou, P.W.Taylor, P.Pulay J.Am.Chem.Soc. 114, 8191-8201(1992) reference for delocalized coordinates: J.Baker, A. Kessi, B.Delley J.Chem.Phys. 105, 192-212(1996) ========================================================== ========================================================== $LIBE group (required if linear bends are used in $ZMAT) A degenerate linear bend occurs in two orthogonal planes, which are specified with the help of a point A. The first bend occurs in a plane containing the atoms I,J,K and the user input point A. The second bend is in the plane perpendicular to this, and containing I,J,K. One such point must be given for each pair of bends used. APTS(1)= x1,y1,z1,x2,y2,z2,... for linear bends 1,2,... Note that each linear bend serves as two coordinates, so that if you enter 2 linear bends (HCCH, for example), the correct value of NZVAR is M-2, where M=3N-6 or 3N-5, as appropriate. ========================================================== ========================================================== $SCF group relevant if SCFTYP = RHF, UHF, or ROHF, required if SCFTYP = GVB) This group of parameters provides additional control over the RHF, UHF, ROHF, or GVB SCF steps. It must be given to define GVB open shell or perfect pairing wavefunctions. See $MCSCF for multireference inputs. DIRSCF = a flag to activate a direct SCF calculation, which is implemented for all the Hartree-Fock type wavefunctions: RHF, ROHF, UHF, and GVB. This keyword also selects direct MP2 computation. The default of .FALSE. stores integrals on disk storage for a conventional SCF calculation. FDIFF = a flag to compute only the change in the Fock matrices since the previous iteration, rather than recomputing all two electron contributions. This saves much CPU time in the later iterations. This pertains only to direct SCF, and has a default of .TRUE. This option is implemented only for the RHF, ROHF, UHF cases. Cases with many diffuse functions in the basis set, or large molecules, may sometimes be "mushy" at the end, rather than converging. Increasing ICUT in $CONTRL by one may help this, or consider turning this parameter off. ---- The next flags affect convergence rates. DIIS = selects Pulay's DIIS interpolation. SOSCF = selects second order SCF orbital optimization. Only one of DIIS or SOSCF may be .TRUE. in any run. Which is chosen by default depends on the run: for RHF, GVB, UHF, or ROHF (if Abelian): SOSCF.TRUE. for any DFT, or for non-Abelian groups: DIIS=.TRUE. NOCONV = .TRUE. means neither SOSCF nor DIIS will be used, specify NOCONV=.TRUE. DIIS=.FALSE. SOSCF=.FALSE. in the rare case you don't wish to use either one. NOCONV's default is .FALSE., meaning the program will obey your choice of DIIS/SOSCF, or else pick its default for them. NOSCAL = a flag to avoid scaling rotations in SOSCF. Once either DIIS or SOSCF are initiated, the following less important accelerators are placed in abeyance: EXTRAP = selects Pople extrapolation of the Fock matrix. DAMP = selects Davidson damping of the Fock matrix. SHIFT = selects level shifting of the Fock matrix. RSTRCT = selects restriction of orbital interchanges. DEM = selects direct energy minimization, which is implemented only for RHF. (default=.FALSE.) defaults for EXTRAP DAMP SHIFT RSTRCT DIIS SOSCF ab initio: T F F F F/T T/F semiempirical: T F F F F F The above parameters are implemented for all SCF wavefunction types, except that DIIS will work for GVB only for those cases with NPAIR=0 or NPAIR=1. FSHIFT = a level static shift parameter of the Fock matrix in Hartree (default = 0.0). CUHF = flag requesting Constrained UHF, which causes the occupied beta orbitals of a UHF run to lie entirely within the occupied alpha orbital space. This produces results identical to high spin ROHF! Obviously, this keyword pertains only when using SCFTYP=UHF. The default is .FALSE., meaning a spin-contaminated ordinary UHF solution is sought. Applicable to UHF, UDFT or MRSF-TDDFT energy and gradients, or to UMP2 energy calculations. ---- These parameters fine tune the various convergers. CONV = SCF density convergence criteria. Convergence is reached when the density change between two consecutive SCF cycles is less than this in absolute value. One more cycle will be executed after reaching convergence. Less accuracy in CONV gives questionable gradients. The default is 1.0d-05, except runs involving CI, MP2, CC, or TDDFT use 1.0d-06 to obtain more crisply converged virtual orbitals. SOGTOL = second order gradient tolerance. SOSCF will be initiated when the orbital gradient falls below this threshold. (default=0.25 au) ETHRSH = energy error threshold for initiating DIIS. The DIIS error is the largest element of e=FDS-SDF. Increasing ETHRSH forces DIIS on sooner. (default = 0.5 Hartree) DIRTHR = energy threshold to force a recomputation of the Fock matrix if energy does not decrease with DIRSCF. (deafult = 0 hartree, disabling the option) MAXDII = Maximum size of the DIIS linear equations, so that at most MAXDII-1 Fock matrices are used in the interpolation. (default=10) SWDIIS = density matrix convergence at which to switch from DIIS to SOSCF. A value of zero means to keep using DIIS at all geometries, which is the default. However, it may be useful to have DIIS work only at the first geometry, in the initial iterations, for example transition metal ECP runs which has a less good Huckel guess, and then use SOSCF for the final SCF iterations at the first geometry, and ever afterwards. A suggested usage might be DIIS=.TRUE. ETHRSH=2.0 SWDIIS=0.005. This option is not programmed for GVB. LOCOPT = When set to .TRUE., SCF options are locked and do not change during the following: SOSCF and DIIS switch (SWDIIS), DFT grid switch, RHF -> DFT switch (SWOFF). If .FALSE., any of these switches resets some SCF options, such as SHIFT or DAMP. (Default: .FALSE.) RESET = In UHF, reset SOSCF or DIIS if energy rises (Default: .TRUE. for UHF, .FALSE. otherwise). DEMCUT = Direct energy minimization will not be done once the density matrix change falls below this threshold. (Default=0.5) DMPCUT = Damping factor lower bound cutoff. The damping damping factor will not be allowed to drop below this value. (default=0.0) note: The damping factor need not be zero to achieve valid convergence (see Hsu, Davidson, and Pitzer, J.Chem.Phys., 65, 609 (1976), see the section on convergence control), but it should not be astronomical either. * * * * * * * * * * * * * * * * * * * * * For more info on the convergence methods, see the 'Further Information' section. * * * * * * * * * * * * * * * * * * * * * ---- orbital modification options ---- The four options UHFNOS, VVOS, MVOQ, and ACAVO are mutually exclusive. The latter 3 require RUNTYP=ENERGY, and should not be used with any correlation treatment. UHFNOS = flag controlling generation of the natural orbitals of a UHF function. (default=.FALSE.) VVOS = flag controlling generation of Valence Virtual Orbitals. See J.Chem.Phys. 120, 2629-2637(2004). The default is .FALSE. VVOs are a quantitative realization of the concept of the "lowest unoccupied molecular orbital". The implementation allows any elements H-Xe (non-relativistic) and H-Rn (relativistic), for RHF, ROHF, and GVB wavefunctions, as well as DFT runs (see also VVOS in $MCSCF). IVVTYP in $LOCAL allows for the explicit selection non-relativistic or relativistic orbitals. VVOs are restricted to all electron basis sets and model core potentials; ECPs may not be used. VVOS should be much better MCSCF starting orbitals than either MVOQ or ACAVO type virtuals. MVOQ = 0 Skip MVO generation (default) = n Form modified virtual orbitals, using a cation with n electrons removed. Implemented for RHF, ROHF, and GVB. If necessary to reach a closed shell cation, the program might remove n+1 electrons. Typically, n will be about 6. = -1 The cation used will have each valence orbital half filled, to produce MVOs with valence-like character in all regions of the molecule. Implemented for RHF and ROHF only. ACAVO = Flag to request Approximate Correlation-Adapted Virtual Orbitals. Implemented for RHF, ROHF, and GVB (w/o direct SCF). Default is .FALSE. PACAVO = Parameters used to define the ACAVO generating operator, which is defined as a*T + b*Vne + c*Jcore + d*Jval + e*Kcore + f*Kval The default, PACAVO(1)=0,0,0,0,0,-1, maximizes the exchange interaction with valence MOs (see for example J.L.Whitten, J.Chem.Phys. 56, 5458-5466(1972). The K-orbitals of D.Feller, E.R.Davidson J.Chem.Phys. 74, 3977-3979(1981) are 0.06*F-K(valence), which is PACAVO(1)= 0.06,0.06,0.12,0.12,-0.06,-1.06. Of course, canonical virtuals are PACAVO(1)=1,1,2,2,-1,-1. SFOVVO = Flag to request the substitution of virtual canonical orbitals for VVOs in energy calculations for SF-ORMAS and SF-ORMAS-PDFT. * * * UHFCHK = a flag to perform a RHF --> UHF wavefunction stability test (relaxation of RHF orbitals to unequal alpha and beta orbitals). This option is implemented only for RHF, and involves testing only pairs of orbitals at a single time: this is not foolproof! default is .FALSE. NHOMO = if UHFCHK is selected, the number of orbitals at The top of the occupied orbital space to be checked for instabilities. Default= -1, meaning check HOMO and HOMO-1. NLUMO = if UHFCHK is selected, the number of orbitals at the bottom of the virtual space checked. Basis sets with diffuse functions should check many more orbitals, to get past the diffuse MOs. Default=2, meaning check LUMO, LUMO+1, LUMO+2. * * * MOM = flag enabling the Molecular Overlap Method (MOM). Currently works only with SCFTYP=UHF and ROHF. MOM computes excitations by selecting orbitals at each iteration to resemble earlier iterations. See J.Phys.Chem. A 112, 13164-13171(2008). (default=.FALSE.) note: MOM typically requires an initial MOREAD set of ground-state MO's, reordered by NORDER=1 in conjunction with the IORDER/JORDER arrays. The MOM flag is an algorithm very similar to that enabled by the RSTRCT flag. Note: RSTRCT works with all SCF types. note for MRSF-TDDFT: When there are nearby triplet states the SCF process can choose a wrong triplet state. This can cause a discontinuity in the potential energy surface. MOM=.t. of ROHF SCF tries to keep the same reference ROHF wavefunction of the initial MOREAD set to avoid the complication. However, one can also use this keyword more aggresively for choosing a particular MO set to be optimized, which would produce new respose states that could not have been obtained without the option. This can be especially useful when there are near-degenerate triplet states. KPROJ determines the flavor of the MOM projections p_j obtained from the overlap matrix O(i,j). KPROJ pertains only if the MOM flag is enabled. = 0 p_j = |sum_i O(i,j)| as given in the original MOM article = 1 p_j = sum_i O(i,j) (default) = 2 p_j = sum_i O(i,j)*O(i,j), as implemented in Q-Chem ----- GVB wavefunction input ----- The next parameters define the GVB wavefunction. See also MULT in the $CONTRL input. The GVB wavefunction assumes orbitals are in the order core, open, pairs. NCO = The number of closed shell orbitals. The default almost certainly should be changed! (default=0). NSETO = The number of sets of open shells in the function. Maximum of 10. (default=0) NO = An array giving the degeneracy of each open shell set. Give NSETO values. (default=0,0,0,...). NPAIR = The number of geminal pairs in the -GVB- function. Maximum of 12. The default corresponds to open shell SCF (default=0). CICOEF = An array of ordered pairs of CI coefficients for the -GVB- pairs. (default = 0.90,-0.20,0.90,-0.20,...) For example, a two pair case for water, say, might be CICOEF(1)=0.95,-0.05,0.95,-0.05. If not normalized, as in the default, CICOEF will be. This parameter is useful in restarting a GVB run, with the current CI coefficients. COUPLE = A switch controlling the input of F, ALPHA, and BETA. (Default=.FALSE.) Input for F, ALPHA, BETA will be ignored unless you select this variable as .TRUE. F = An vector of fractional shell occupations. ALPHA = An array of A coupling coefficients, given in lower triangular order. BETA = An array of B coupling coefficients, given in lower triangular order. Note: The default for F, ALPHA, and BETA depends on the state chosen. Defaults for the most commonly occurring cases are internally stored. See "Further Information" for other cases, including degenerate open shells. Note: ALPHA and BETA can be given for -ROHF- orbital canonicalization control, see "Further Information". ----- miscellaneous options ----- NPUNCH = option for output to the PUNCH file = 0 do not punch out the final orbitals = 1 punch out the occupied orbitals = 2 punch out occupied and virtual orbitals The default is NPUNCH = 2. NPREO = energy and orbital printing options, applied after other output options, for example NPRINT=-5 for no orbital output overrules NPREO. NPREO affects only printing, see NPUNCH just above. Default: 1,9999,2,1 (meaning print all orbitals, but no separate list of orbital energies). Orbitals from NPREO(1) to NPREO(2) and orbital energies from NPREO(3) to NPREO(4) are printed. Positive values are explicit orbitals, while negative numbers are relative to the HOMO. Here, HOMO means NE/2, ie RHF-like counting, no matter what the SCFTYP actually is: NPREO(1)=25,35 prints orbitals 25 to 35 NPREO(1)=-3,-4 prints 8 orbitals: from 3 below the HOMO, the HOMO, to 4 above the HOMO. NPREO(3) to NPREO(4) define separate print-out of the orbital energies, by default this is skipped, since the starting value is higher than the end. To print only HOMO and "LUMO" LCAO coefficients, and all orbital energies, enter: NPREO(1)=0,-1,1,9999 ----- options for virial scaling ----- VTSCAL = A flag to request that the virial theorem be satisfied. An analysis of the total energy as an exact sum of orbital kinetic energies is printed. The default is .FALSE. This option is implemented for RHF, UHF, and ROHF, for RUNTYP=ENERGY, OPTIMIZE, or SADPOINT. Related input is: SCALF = initial exponent scale factor when VTSCAL is in use, useful when restarting. The default is 1.0. MAXVT = maximum number of iterations (at a single geometry) to satisfy the energy virial theorem. The default is 20. VTCONV = convergence criterion for the VT, which is satisfied when 2 + + R x dE/dR is less than VTCONV. The default is 1.0D-6 Hartree. For more information on this option, which is most useful during a geometry search, see M.Lehd and F.Jensen, J.Comput.Chem. 12, 1089-1096(1991). * * * * * * * * * * * * * * * * * * * For more discussion of GVB/ROHF input see the 'further information' section * * * * * * * * * * * * * * * * * * * ========================================================== ========================================================== $SCFMI group (optional, relevant if SCFTYP=RHF) The Self Consistent Field for Molecular Interactions (SCF-MI) method is a modification of the usual Roothaan equations that avoids basis set superposition error (BSSE) in intermolecular interaction calculations, by expanding each monomer's orbitals using only its own basis set. Thus, the resulting orbitals are not orthogonal. The presence of a $SCFMI group in the input triggers the use of this option. The implementation is limited to ten monomers, treated at the RHF level. The energy, gradient, and therefore semi-numerical hessian are available. The SCF step may be run in direct SCF mode, and parallel calculation is also enabled. The calculation must use Cartesian Gaussian AOs only, not spherical harmonics. The SCF-MI driver differs from normal RHF calculations, so not all converger methods are available. Finally, this option is not compatible with electron correlation treatments (DFT, MP2, CI, or CC). The first 3 parameters must be given. All atoms of a fragment must appear consecutively in $DATA. NFRAGS = number of distinct fragments present. Both the supermolecule and its constituent monomers must be well described as closed shells by RHF wavefunctions. NF = an array containing the number of doubly occupied MOs for each fragment. MF = an array containing the number of atomic basis functions located on each fragment. ITER = maximum number of SCF-MI cycles, overriding the usual MAXIT value. (default is 50). DTOL = SCF-MI density convergence criteria. (default is 1.0d-10) ALPHA = possible level shift parameter. (default is 0.0, meaning shifting is not used) DIISON = a flag to active the DIIS convergence. (default is .TRUE.) MXDIIS = the maximum number of previous effective Fock and overlap matrices to be used in DIIS (default=10) DIISTL = the density change value at which DIIS starts. (default=0.01) A Huckel guess is localized by the Boys procedure onto each fragment to provide starting orbitals for each: ITLOC = maximum number of iteration in the localization step (Default is 50) CNVLOC = convergence parameter for the localization. (default is .01). IOPT = prints additional debug information. = 0 standard outout (default) = 1 print for each SCF-MI cycle MOs, overlap between the MOs, CPU times. = 2 print some extra informations in secular systems solution. ========================================================== "Modification of Roothan Equations to exclude BSSE from Molecular Interaction calculations" E. Gianinetti, M. Raimondi, E. Tornaghi Int. J. Quantum Chem. 60, 157-166 (1996) "Implementation of Gradient optimization algorithms and Force Constant computations in BSSE-free direct and conventional SCF approaches" A. Famulari, E. Gianinetti, M. Raimondi, M. Sironi Int. J. Quantum Chem. 69, 151-158 (1997) ========================================================== $DFT group (relevant if DFTTYP is chosen) (relevant if SCFTYP=RHF,UHF,ROHF) Note that if DFTTYP=NONE, an ab initio calculation will be performed, rather than density functional theory. This group permits the use of various one electron (usually empirical) operators instead of the true many electron Hamiltonian. Two programs are provided, METHOD= GRID or GRIDFREE. The programs have different functionals available, and so the keyword DFTTYP (which is entered in $CONTRL) and other associated inputs are documented separately below. Every functional that has the same name in both lists is an identical functional, but each METHOD has a few functionals that are missing in the other. The grid free implementation is based on the use of the resolution of the identity to simplify integrals so that they may be analytically evaluated, without using grid quadratures. The grid free DFT computations in their present form have various numerical errors, primarily in the gradient vectors. Please do not use the grid-free DFT program without reading the discussion in the 'Further References' section regarding the gradient accuracy. The grid based DFT uses a typical grid quadrature to compute integrals over the rather complicated functionals, using two possible angular grid types. Achieving a self-consistent field with DFT is rather more difficult than for normal HF, so DIIS is the default converger. Both DFT programs will run in parallel. See the two lists below for possible functionals in the two programs. See also the $TDDFT input group for excited states. METHOD = selects grid based DFT or grid free DFT. = GRID Grid based DFT (default) = GRIDFREE Grid free DFT DFTTYP is given in $CONTRL, not here in $DFT! Possible values for the grid-based program are listed first, ----- options for METHOD=GRID ----- DFTTYP = NONE means ab initio computation (default) Many choices are given below, perhaps the most sensible are local DFT: SVWN pure DFT GGA: BLYP, PW91, B97-D, PBE/PBEsol hybrid DFT GGA: B3LYP, X3LYP, PBE0 pure DFT meta-GGA: revTPSS hybrid DFT meta-GGA: TPSSh, M06 but of course, everyone has their own favorite! pure exchange functionals: = SLATER Slater exchange = BECKE Becke 1988 exchange = GILL Gill 1996 exchange = OPTX Handy-Cohen exchange = PW91X Perdew-Wang 1991 exchange = PBEX Perdew-Burke-Ernzerhof exchange These will be used with no correlation functional at all. pure correlation functionals: = VWN Vosko-Wilk-Nusair correlation, using their electron gas formula 5 (aka VWN5) = VWN3 Vosko-Wilk-Nusair correlation, using their electron gas formula 3 = VWN1RPA Vosko-Wilke-Nusair correlation, using their e- gas formula 1, with RPA params. = PZ81 Perdew-Zener 1981 correlation = P86 Perdew 1986 correlation = LYP Lee-Yang-Parr correlation = PW91C Perdew-Wang 1991 correlation = PBEC Perdew-Burke-Ernzerhof correlation = OP One-parameter Progressive correlation These will be used with 100% HF exchange, if chosen. combinations (partial list): = SVWN SLATER exchange + VWN5 correlation Called LDA/LSDA in physics for RHF/UHF. = SVWN1RPA Slater exchange + VWN1RPA correlation = BLYP BECKE exchange + LYP correlation = BOP BECKE exchange + OP correlation = BP86 BECKE exchange + P86 correlation = GVWN GILL exchange + VWN5 correlation = GPW91 GILL exchange + PW91 correlation = PBEVWN PBE exchange + VWN5 correlation = PBEOP PBE exchange + OP correlation = OLYP OPTX exchange + LYP correlation = PW91 means PW91 exchange + PW91 correlation = PBE means PBE exchange + PBE correlation There's a nearly infinite set of pairings (well, 6*9), so we show only enough to give you the idea. In other words, pairs are formed by abbreviating the exchange functionals SLATER=S, BECKE=B, GILL=G, OPTX=O, PW91X=PW91, PBEX=PBE and matching them with any correlation functional, of which only two are abbreviated when used in combinations, PW91C==>PW91, PBEC==>PBE The pairings shown above only scratch the surface, but clearly, many possibilities, such as PW91PBE, are nonsense! pure DFT GGA functionals: = EDF1 empirical density functional #1, which is a modified BLYP from Adamson/Gill/Pople. = PW91 Perdew/Wang 1991 = PBE Perdew/Burke/Ernzerhof 1996 = revPBE PBE as revised by Zhang/Yang = RPBE PBE as revised by Hammer/Hansen/Norskov = PBEsol PBE as revised by Perdew et al for solids = HCTH Hamprecht/Cohen/Tozer/Handy's 1998 mod to B97, omitting HF exchange (m=4) = HCTH93 Hamprecht/Cohen/Tozer/Handy's 1998 mod to B97, omitting HF exchange, fitting to 93 atoms and molecules = HCTH120 later fit to 120 systems = HCTH147 later fit to 147 systems = HCTH407 later fit to 407 systems (best) = HCTH407P also known as HCTH/407+, designed by Boese/Chandra/Martin/Marx for describing hydrogen bonds = HCTH14 also known as HCTH p=1/4, made by Menconi/Wilson/Tozer = HCTH76 also known as HCTH p=7/6, made by Menconi/Wilson/Tozer = SOGGA PBE revised by Zhao/Truhlar for solids = MOHLYP metal optimized OPTX, half LYP = B97-D Grimme's modified B97, with dispersion correction (this forces DC=.TRUE.) = SOGGA11 optimized with broad applicability for chemistry, by Peverati/Zhao/Truhlar hybrid GGA functionals: = BHHLYP HF and BECKE exchange + LYP correlation = B3PW91 Becke's 3 parameter exchange hybrid, with PW91 correlation functional = B3LYP this is a hybrid method combining five functionals: Becke + Slater + HF exchange (B3), with LYP + VWN5 correlation. B3LYPV5 is a synonym for B3LYP. = B3LYPV1R use VWN1RPA in place of VWN5, matches the e- gas formula chosen by some programs. = B3LYPV3 use VWN3 in place of B3LYP's VWN5 = B3P86 B3-type exchange, P86 correlation, using VWN3 as the LDA part of the correlation. B3P86V3 is a synonym for B3P86. = B3P86V1R use VWN1RPA in place of VWN3 = B3P86V5 use VWN5 in place of VWN3 = B97 Becke's 1997 hybrid functional = B97-1 Hamprecht/Cohen/Tozer/Handy's 1998 reparameterization of B97 = B97-2 Wilson/Bradley/Tozer's 2001 mod to B97 = B97-3 Keal/Tozer's 2005 mod to B97 = B97-K Boese/Martin's 2004 mod for kinetics = B98 Schmider/Becke's 1998 mode to B97, using their best "2c" parameters. = PBE0 a hybrid made from PBE = X3LYP HF+Slater+Becke88+PW91 exchange, and LYP+VWN1RPA correlation. = SOGGA11X a hybrid based on SOGGA11, with 40.15% HF exchange. = APF mixed functional based on PBE0 and B3PW91 Each includes some Hartree-Fock exchange, and also may use a linear combination of many DFT parts. range separated functionals: These are also known as "long-range corrected functionals". LC-BVWN, LC-BOP, LC-BLYP, or LC-BPBE are available by selecting BVWN, BOP, BLYP, or BPBE and also setting the flag LC=.TRUE. (see LC and also MU below). Others are selected by their specific name, without using LC: = CAMB3LYP coulomb attenuated B3LYP = wB97 omega separated form of B97 = wB97X wB97 with short-range HF exchange = wB97X-D dispersion corrected wB97X = CAMQTP00 Verma/Bartlett's reparametrization of CAMB3LYP for vertical ionization potential = CAMQTP01 Jin/Bartlett's reparametrization of CAMB3LYP for vertical ionization potential M11 is also range-separated, but is listed below with the other meta-GGAs. "double hybrid" GGA: = B2PLYP mixes BLYP, HF exchange, and MP2! See related inputs CHF and CMP2 below. "double hybrid" and "range separated": = wB97X-2 intended for use with GBASIS=CCT,CCQ,CC5 = wB97X-2L intended for use with GBASIS=N311 NGAUSS=6 NDFUNC=3 NFFUNC=1 NPFUNC=3 DIFFSP=.T. DIFFS=.T. Note: there are no analytic gradients for "double hybrids". Note: the B2PLYP family uses the conventional MP2 energy and may be used for closed shell or spin-unrestricted open shell cases. The wB97X-2 family uses the SCS-MP2 energy, and thus is limited to closed shell cases at present. meta-GGA functionals: These are not hybridized with HF exchange, unless that is explicitly stated below. = VS98 Voorhis/Scuseria, 1998 = PKZB Perdew/Kurth/Zupan/Blaha, 1999 = tHCTH Boese/Handy's 2002 metaGGA akin to HCTH = tHCTHhyb tHCTH's hybrid with 15% HF exchange = BMK Boese/Martin's 2004 parameterization of tHCTHhyb for kinetics = TPSS Tao/Perdew/Staroverov/Scuseria, 2003 = TPSSh TPSS hybrid with 10% HF exchange = TPSSm TPSS with modified parameter, 2007 = revTPSS revised TPSS, 2009 = dlDF a reparameterized M05-2X, reproducing interaction energies which have had all dispersion removed. This MUST be used with a special -D correction to recover dispersion. See 'Further References'. = M05 Minnesota exchange-correlation, 2005 a hybrid with 28% HF exchange. = M05-2X M05, with doubled HF exchange, to 56% = M06 Minnesota exchange-correlation, 2006 a hybrid with 27% HF exchange. = M06-L M06, with 0% HF exchange (L=local) = M06-2X M06, with doubled HF exchange, to 54% = M06-HF M06 correlation, using 100% HF exchange = M08-HX M08 with 'high HF exchange' = M08-SO M08 with parameters that enforce the correct second order gradient expansion. = M11 M11 range-separated hybrid = M11-L M11 local (0% HF exchange) with dual-range exchange = MN12-L uses a nonseparable functional form aiming to provide balanced peformance for both chemistry and solid-state physics applications = MN12-SX screened-exchanged (SX) hybrid functional with 20% HF exchange for short-range and 0% HF exchange for long-range = MN15 uses a nonseparable functional form with 44% HF exchange. This fuctional is a global hybride (no range-separation). = MN15-L a local functional with 0% HF exchange. = REVM06 revised hybrid M06 = REVM06-L revised local M06-L = REVM11 revised M11 with 40% HF exchange More than one hundred functionals are available by: = USELIBXC See more details in $LIBXC When the M06 family was created, Truhlar recommended M06 for the general situation, but see his "concluding remarks" in the M06 reference about which functional is best for what kind of test data set. The most recent M11 family is probably a better choice, and two functionals fit all the needs of the older M05/M06/M08 families. https://en.wikipedia.org/wiki/Minnesota_functionals An extensive bibliography for all functionals can be found in the 'Further References' section of this manual. Note that only a subset of these functionals can be used for TD-DFT energy or gradients. These subsets are listed in the $TDDFT input group. * * * dispersion corrections * * * Many exchange-correlation functionals fail to compute intra- and inter-molecular dispersion interactions accurately. Two possible correction schemes are provided below, DC and LRDFLG. DC uses empirically chosen C6 and C8 coefficients, while LRDFLG obtains these from the molecular DFT densities. At most, only one of the LRDFLG or DC options below may be chosen. DC = a flag to turn on Grimme's empirical dispersion correction, involving scaled R**(-6) terms. N.B. This empiricism may also be added to plain Hartree-Fock, by choosing DFTTYP=NONE with DC=.T. Four different versions exist, see IDCVER. (default=.FALSE., except if DFTTYP=B97-D, wB97X-D) IDCVER = 1 means 1st 2004 implementation. = 2 means 2nd 2006 implementation DFT-D2, default for B97-D, wB97X-D. = 3 means 3rd 2010 implementation DFT-D3. Default if DC is chosen and IDCVER isn't given. = 4 means modified 3rd implementation DFT-D3(BJ). (-4 is used for DFT-D3(BJ) for HF-3c). Setting IDCVER will force DC=.TRUE. GCP = a flag for the geometric counterpoise scheme correction in HF-3c. SRB = a flag for short-range basis set incompleteness (SRB) correction in HF-3c. DCCHG = a flag to use Chai-Head-Gordon damping function instead of Grimme's 2006 function. Pertinent only for the DFT-D2 method. Forces DC=.TRUE. (default=.FALSE. except for wB97X-D) DCABC = a flag to turn on the computation of the E(3) non- additive energy term. Pertinent only for DFT-D3, it forces DC=.TRUE. (default=.FALSE.) The following parameters govern Grimme's semiempirical dispersion term. They are basis set and functional dependent, so they exist for only a few DFTTYP. Default values are automatically selected and printed out in the output file for many common density functionals. The following keywords are for entering non-standard values. For DFT-D2 values, see also: R.Peverati and K.K.Baldridge J.Chem.Theory Comput. 4, 2030-2048 (2008). For DFT-D3 values, and a detailed explanation of each parameter, see: S. Grimme, J. Antony, S. Ehrlich and H. Krieg, J.Chem.Phys. 132, 154104/1-19(2010) and for DFT-D3(BJ): S. Grimme, S. Ehrlich and L. Goerigk, J.Comput.Chem. 32, 1456-1465 (2011) DCALP = alpha parameter in the DFT-D damping function (same as alpha6 in Grimme's DFT-D3 notation). Note also that alpha8 and alpha10 in DFT-D3 have constrained values of: alpha8 = alpha6 + 2, alpha10 = alpha8 + 2. Default=14.0 for DFT-D3 =20.0 for DFT-D2 =23.0 for DFT-D1 =6.00 for DCCHG=.TRUE. DCSR = sR exponential parameter to scale the van der Waals radii (same as sR,6 in Grimme's DFT-D3 notation). Note also that sR,8 in DFT-D3 have a fixed value of 1.0. Optimized values are automatically selected for some of the more common functionals, otherwise, the default is 1.00 for DFT-D3, 1.10 for DFT-D2, and 1.22 for DFT-D1. DCS6 = s6 linear parameter for scaling the C6 term. Optimized values are automatically selected for some of the more common functionals, otherwise, the default is 1.00. DCS8 = s8 linear parameter for scaling the C8 term of DFT-D3. Pertinent only for DFT-D3. Optimized values are automatically selected for some of the more common functionals, otherwise, the default is 1.00. DCA1 = a1 parameter appearing in the -D3(BJ) dispersion model. Optimized values are automatically selected for a set of known functionals, otherwise the default is 0.50. DCA2 = a2 parameter appearing in the -D3(BJ) dispersion model. Optimized values are automatically selected for a set of known functionals, otherwise the default is 4.00. The old keywords DCPAR and DCEXP were replaced by DCS6 and DCSR in 2010. Similarly, DCOLD has morphed into IDCVER. - - - The Local Response Dispersion (LRD) correction includes atomic pair-wise -C6/R**6, -C8/R**8, and -C10/R**10 terms, whose coefficients are computed from the molecular system's electron density and its nuclear gradient. The nuclear gradient assumes the dispersion coefficients do not vary with geometry, which causes only a very small error in the gradient. Optionally, 3 and 4 center terms may be added, at the 1/R**6 level; in this case, nuclear gradients may not be computed at all. Since the three numerical parameters are presently known only for the long-range exchange corrected BOP functional, calculations may specify simply DFTTYP=LCBOPLRD. The "LCBOPLRD" functional will automatically select the following: DFTTYP=BOP LC=.TRUE. MU=0.47 LRDFLG=.TRUE. LAMBDA=0.232 KAPPA=0.600 RZERO=3.22 leaving only the choice for MLTINT up to you. References for LRD are T.Sato, H.Nakai J.Chem.Phys. 131, 224104/1-12(2009) T.Sato, H.Nakai J.Chem.Phys. 133, 194101/1-9(2010) LRDFLG = flag choosing the Local Response Dispersion (LRD) C6, C8, and C10 corrections. Default=.FALSE. MLTINT = flag to add the 3 and 4 center 6th order terms, the default=.FALSE. Note that nuclear gradients are not available if these multi-center terms are requested. Three numerical parameters may be input. The defaults shown are optimized for the BOP functional with the LC correction for long-range exchange. LAMBDA = parameter adjusting the density gradient correction for the atomic and atomic pair polarizabilities. (default=0.232) KAPPA = parameter in the damping function (default=0.600) RZERO = parameter in the damping function (default=3.22) It may be interesting to see a breakdown of the total dispersion correction, using these keywords: PRPOL = print out atomic effective polarizabilities (default=.FALSE.) PRCOEF = N (default N=0) print out dispersion coefficient to N-th order. PRPAIR = print out atomic pair dispersion energies (default=.FALSE.) * * * range separation * * * LC = flag to turn on the long range correction (LC), which smoothly replaces the DFT exchange by the HF exchange at long inter-electron distances. (default=.FALSE.) This option can be used only with the Becke exchange functional (Becke) and a few correlation functionals: DFTTYP=BVWN, BOP, BLYP, BPBE only. For example, B3LYP has a fixed admixture of HF exchange, so it cannot work with the LC option. See H.Iikura, T.Tsuneda, T.Yanai, and K.Hirao, J.Chem.Phys. 115, 3540 (2001). MU = A parameter for the long range correction scheme. Increasing MU increases the HF exchange used, very small MU produces the DFT limit. (default=0.33) Other range-separated options exist, invoked by naming the functional, such as DFTTYP=CAMB3LYP (see the DFTTYP keyword for a full list). * * * B2x-PLYP double hybrid functionals * * * B2xPLYP Double Hybrid functionals have the general formula: Exc = (1-cHF) * ExGGA + cHF * ExHF + (1-cMP2) * EcGGA + cMP2 * E(2) The next keywords allow the choice of cHF and cMP2. Both values must be between 0 and 1 (0-100%). CHF = amount of HF exchange. (default=0.53) CMP2 = amount of MP2. (default=0.27) Some other common double hybrid functionals are available simply by choosing DFTTYP=B2PLYP, and changing the CHF and CMP2 parameters. Popular parametrizations are: CHF CMP2 ------------------------------------------ B2-PLYP (default) | 0.53 | 0.27 | ------------------------------------------ B2K-PLYP | 0.72 | 0.42 | ------------------------------------------ B2T-PLYP | 0.60 | 0.31 | ------------------------------------------ B2GP-PLYP | 0.65 | 0.36 | ------------------------------------------ * * * Grid Input * * * Only one of the three grid types may be chosen for the run. The default (if no selection is made) is the Lebedev grid. In order to duplicate results obtained prior to April 2008, select the polar coordinate grid NRAD=96 NTHE=12 NPHI=24. Energies can be compared if and only if the identical grid type and density is used, analogous to needing to compare with the identical basis set expansions. See REFS.DOC for more information on grids. See similar inputs in $TDDFT. RADTYP = type of radial quadrature = MHL : Murray, Handy and Laming radial grid (default) = MK3 : Mura and Knowles Log3 radial grid = TA : Treutler and Ahlrichs radial grid Lebedev grid: NRAD = number of radial points in the Euler-MacLaurin quadrature. (default=96) NLEB = number of angular points in the Lebedev grids. (default=302). Possible values are 86, 110, 146, 170, 194, 302, 350, 434, 590, 770, 974, 1202, 1454, 1730, 2030... The same switch can be used to select spherical grids developed by A.S.Popov: 246, 264, 342, 432 : rotaional octahedral symmetry grids 350 and 398 : full octahedral symmetry grids (like Lebedev grids) 212 : rotational icosahedral grid rotational icosahedral grids with 192 and 242 points are currently disabled in the source code Meta-GGA functionals require a tighter grid to achieve the same accuracy. For this reason a tighter default grid of NRAD=99 and NLEB=590 is chosen by default with all meta-GGA functionals. The default for NLEB means that nuclear gradients will be accurate to about the default OPTTOL=0.00010 (see $STATPT), 590 approaches OPTTOL=0.00001, and 1202 is "army grade". The next two specify radial/angular in a single keyword: SG1 = a flag to select the "standard grid 1", which has 24 radial points, and various pruned Lebedev grids, from 194 down to 6. (default=.FALSE. This grid is very fast, but produces gradients whose accuracy reaches only OPTTOL=0.00050. This grid should be VERY USEFUL for the early steps of a geometry optimization. JANS = two unpublished grids due to Curtis Janssen, implemented here differently than in MPQC: = 1 uses 95 radial points for all atoms, and prunes from a Lebedev grid whose largest size is 434, thus using about 15,000 grid points/atom. = 2 uses 155 radial points for all atoms, and prunes from a Lebedev grid whose largest size is 974, thus using about 71,000 grid points/atom. This is a very accurate grid, e.g. "army grade". The information for pruning exists only for H-Ar, so heavier elements will use the large radial/ Lebedev grid without any pruning. polar coordinate grid: NRAD = number of radial points in the Euler-MacLaurin quadrature. (96 is reasonable) NTHE = number of angle theta grids in Gauss-Legendre quadrature (polar coordinates). (12 is reasonable) NPHI = number of angle phi grids in Gauss-Legendre quadrature. NPHI should be double NTHE so points are spherically distributed. (24 is reasonable) The number of angular points will be NTHE*NPHI. The values shown give a gradient accuracy near the default OPTTOL of 0.00010, while NTHE=24 NPHI=48 approaches OPTTOL=0.00001, and "army grade" is NTHE=36 NPHI=72. * * * Grid Switching * * * At the first geometry of the run, pure HF iterations will be performed, since convergence of DFT is greatly improved by starting with the HF density matrix. After DFT engages, most runs (at all geometries, except for PCM or numerical Hessians) will use a coarser grid during the early DFT iterations, before reaching some initial convergence. After that, the full grid will be used. Together, these switchings can save considerable CPU time. SWOFF = turn off DFT, to perform pure SCF iterations, until the density matrix convergence falls below this threshold. This option is independent of SWITCH and can be used with or without it. It is reasonable to pick SWOFF > SWITCH > CONV in $SCF. SWOFF pertains only to the first geometry that the run computes, and is automatically disabled if you choose GUESS=MOREAD to provide initial orbitals. The default is 5.0E-3. SWITCH = when the change in the density matrix between iterations falls below this threshhold, switch to the desired full grid (default=3.0E-4) This keyword is ignored if the SG1 grid is used. Setting to zero disables DFT grid switching. NRAD0 = same as NRAD, but defines initial coarse grid. default = smaller of 24 and NRAD/4 NLEB0 = same as NLEB, but defines initial coarse grid. default = 110 NTHE0 = same as NTHE, but defines initial coarse grid. default = smaller of 8, NTHE/3 NPHI0 = same as NPHI, but defines initial coarse grid. default = smaller of 16, NPHI/3 molecular grid construction parameters: BFCTYP = specify algorithm of molecular grid construction from atomic grids (Becke's fuzzy cell algorithm) = BECKE : use Becke's fuzzy cell method from the original paper. = SSF : use algorithm by Stratmann, Scuseria and Frisch. It usually screens out more points than Becke's algorithm = LEGACY : use legacy DFT code (default) BECKE and SSF switches allow multithreaded DFT run, but does not support DC, LRD and analytical Hessian calculation types, as well as polar spherical grid with non-C1 symmetry. PARTFN = select grid partitioning function relevant only if BFCTYP=SSF or BECKE = BECKE : use polynomial proposed in Becke's original work. Default when BFCTYP=BECKE = SSF : use polynomial proposed by Stratmann, Scuseria and Frisch. Default when BFCTYP=SSF = ERF : use erf-based partitioning function based on unpublished ERF1 algorithm in NWChem = SMSTEP2 : use smoothstep-2 polynomial (very effective, but aggressive partitioning) = SMSTEP3 : use smoothstep-3 polynomial = SMSTEP4 : use smoothstep-4 polynomial = SMSTEP5 : use smoothstep-5 polynomial technical parameters: THRESH = threshold for ignoring small contributions to the Fock matrix. The default is designed to produce no significant energy loss, even when the grid is as good as "army grade". If for some reason you want to turn all threshhold tests off, of course requiring more CPU, enter 1.0e-15. default: 1.0e-4/Natoms/NRAD/NTHE/NPHI GTHRE = threshold applied to gradients, similar to THRESH. < 1 assign this value to all thresholds = 1 use the default thresholds (default). > 1 divide default thresholds by this value. If you wish to increase accuracy, set GTHRE=10. The default introduces an error of roughly 1e-7 (a.u./bohr) in the gradient. WTDER = switch on/off grid weigths derivative contribution to the nuclear gradient. = .TRUE. : compute grid weights derivatives = .FALSE. : do not compute grid weights derivatives Default is .TRUE. if BFCTYP=LEGACY, and .FALSE. otherwise. Note: BFCTYP=BECKE and BFCTYP=SSF currently do not support grid weight derivatives. Grid weights are only important for low-quality grids. The default grid settings are good enough to disable grid weight derivatives calculation. This switch can be turned off to improve convergence of geometry optimization, if the optimized structure is not the one with lowest energy. The keyword $DFTTYP is given in $CONTRL, and may have these values if the grid-free program is chosen: ----- options for METHOD=GRIDFREE ----- DFTTYP = NONE means ab initio computation (default) exchange functionals: = XALPHA X-Alpha exchange (alpha=0.7) = SLATER Slater exchange (alpha=2/3) = BECKE Becke's 1988 exchange = DEPRISTO Depristo/Kress exchange = CAMA Handy et al's mods to Becke exchange = HALF 50-50 mix of Becke and HF exchange correlation functionals: = VWN Vosko/Wilke/Nusair correlation, formula 5 = PWLOC Perdew/Wang local correlation = LYP Lee/Yang/Parr correlation exchange/correlation functionals: = BVWN Becke exchange + VWN5 correlation = BLYP Becke exchange + LYP correlation = BPWLOC Becke exchange + Perdew/Wang correlation = B3LYP hybrid HF/Becke/LYP using VWN formula 5 = CAMB CAMA exchange + Cambridge correlation = XVWN Xalpha exchange + VWN5 correlation = XPWLOC Xalpha exchange + Perdew/Wang correlation = SVWN Slater exchange + VWN5 correlation = SPWLOC Slater exchange + PWLOC correlation = WIGNER Wigner exchange + correlation = WS Wigner scaled exchange + correlation = WIGEXP Wigner exponential exchange + correlation AUXFUN = AUX0 uses no auxiliary basis set for resolution of the identity, limiting accuracy. = AUX3 uses the 3rd generation of RI basis sets, These are available for the elements H to Ar, but have been carefully considered for H-Ne only. (DEFAULT) THREE = a flag to use a resolution of the identity to turn four center overlap integrals into three center integrals. This can be used only if no auxiliary basis is employed. (default=.FALSE.) ========================================================== ========================================================== $TDDFT group (relevant if TDDFT chosen in $CONTRL) This group generates molecular excitation energies by time-dependent density functional theory computations (or time-dependent Hartree-Fock, also known as the Random Phase Approximation). The functional used for the excited states is necessarily the same one that is used for the reference state, specified by DFTTYP in $CONTRL. For conventional TD-DFT (TDDFT=EXCITE in $CONTRL), the orbitals are optimized for RHF or UHF type reference states. Analytic nuclear gradients are available for singlet excited states, while the energy of excited states of other multiplicities can be computed. Two-photon absorption cross-sections may be predicted for singlet excited states. Ground state hyperpolarizabilities may be computed with the TDDFT module. For spin-flip TD-DFT (TDDFT=SPNFLP in $CONTRL), the calculation obtains orbitals for a reference state of either UHF or ROHF type, with MULT in $CONTRL determining the Ms quantum number of the reference. The reference state's Ms is set equal to the S value implied by $CONTRL's MULT=2S+1. The SF-TD-DFT then uses only determinants with Ms=S-1, due to the flip of one alpha spin into a beta spin. This means that target states (which are spin contaminated) will have multiplicities around the range S-1 to S, only. It is quite possible for some of the target states to have a lower energy than the reference!!! Nuclear gradients and properties are available. Like spin-flip TD-DFT, the orbitals for a reference state of Mixed-Reference Spin-Flip (MRSF) TDDFT (TDDFT=MRSF in $CONTRL) are obtained by either UHF or ROHF type, with MULT in $CONTRL determining its Ms quantum number. Only triplet reference (MULT=3 in $CONTRL) state is supported now. While spin-flip TDDFT only utilizes Ms=+1 of triplet state in the generation of response states, MRSF-TDDFT takes both Ms=+1 and -1, which eliminates a bulk of the notorious spin-contamination of SF-TDDFT. Accordingly, MRSF-TDDFT utilizes both alpha to beta as well as beta to alpha spin-filp excitations in the generation of response determinants with Ms=0, 2. Using these, the spin state specific response target states with singlet, triplet as well as quintet spin states can be generated indivisually. Thus, target states have S=0, 1 with Ms=0 and S=2 with Ms=2 (MULT=1, 3, 5 in $TDDFT). Analytic nuclear energy gradients, (RUNTYP=GRADIENT in $CONTRL), 2-state and 3-state conical-intersection search (see $CONICL), ionization potential and electron affinity based on extended Koopmans' theorem, nonadiabatic MD simulation (see $MD), and some properties are available with MRSF-TDDFT. However, gradient is only available for ROHF. The respective papers are Energy and Gradient of MRSF-TDDFT: Lee S, Filatov M, Lee S, Choi CH. J.Chem.Phys.2018 Sep 14;149(10):104101. Lee S, Kim EE, Nakata H, Lee S, Choi CH. J.Chem.Phys.2019 May 14;150(18):184111. 2- and 3 States Conical Intersection: Baek YS, Lee S, Filatov M, Choi CH. J.Phys.Chem.A2021 Mar 2;125(9):1994-2006. Ionization Potential: Pomogaev V, Lee S, Shaik S, Filatov M, Choi CH. J.Phys.Chem.Lett.2021 Oct 7;12:9963-72. Nonadiabatic MD: Park W, Lee S, Huix-Rotllant M, Filatov M, Choi CH. J.Phys.Chem.Lett.2021 Apr 30;12(18):4339-46. =========================================================== See just below for "limitations" below regarding the two different TD-DFT types. TD-DFT is a single excitation theory. All of the caveats listed in the $CIS input group about states with double excitation character, need for Rydberg basis sets, greatly different topology of excited state surfaces, and so on apply here as well. Please read the introduction to the $CIS input group! If you use very large or very small Gaussian exponents, you may need to increase the number of radial grid points (the program prints advice in such cases). TDHF, TDDFT, and CIS are related in the following way: -- Tamm/Dancoff approximation --> | TDHF CIS DFT | V TDDFT TDDFT/TDA Here TDHF means absorption of photons, to produce excited states (TDHF is called RPA in the physics community). This meaning of TDHF should not be confused with the photon scattering processes computed by RUNTYP=TDHF or TDHFX, which generate polarizabilities. Note, in particular, that CITYP=CIS is equivalent to using TDDFT=EXCITE DFTTYP=NONE TAMMD=.TRUE., provided the former is run with no frozen cores. Solvent effects for CIS calculations are therefore available via the TDDFT codes. Excited state properties are calculated using the TDDFT excited state electronic density only during gradient runs, or by setting TDPRP below. The TD-DFT codes excite all electrons, that is, there is no frozen core concept. Please see the 4th chapter of this manual for more information on both types of TD-DFT. "limitations" for TDDFT=EXCITE: Permissible values for DFTTYP are shown below. These include "NONE" which uses TDHF (i.e. the Random Phase Approximation), noting that extra states may need to be solved for in order to be sure of getting the first few states correctly. If nuclear gradients are needed, you may choose any of the following functionals: NONE SVWN, SOP, SLYP, OLYP, BVWN, BOP, BLYP (and their LC=.TRUE. versions) B3LYP, CAMB3LYP, B3LYPV1R, PBE, PBE0 CAMQTP00, CAMQTP01 For evaluation of just the excitation energies, you may use many more functionals, notably including the metaGGAs in the last three lines: NONE SVWN, SVWN1RPA, SPZ81, SP86, SOP, SLYP, BVWN, BVWN1RPA, BPZ81, BP86, BOP, BLYP, OLYP, B3LYP, CAMB3LYP, B3LYPV1R, B3PW91, X3LYP, PW91, PBE, PBE0, APF, VS98, PKZB, M05, M05-2X, M06, M06-HF, M06-L, M06-2X, M08-HX, M08-SO, MN12-L, MN12-SX, MN15, MN15-L, TPSS, TPSSm, TPSSh, and revTPSS The LC flag in $DFT automatically carries over to TDDFT runs. The LC option may be used with the "B" functionals, and (like the similar range-separated CAMB3LYP) is useful in obtaining better descriptions for charge-transfer excitations or Rydberg excitation energies than are the conventional exchange correlation functionals (whether pure or hybrid). The LC flag is also available for excited state gradient computation. Limits specific to the references for TDDFT=EXCITE are: For SCFTYP=RHF, excitation energies can be found for singlet or triplet coupled excited states. For singlet excited states only, analytic gradients and properties can be found, for either full TD-DFT or in the Tamm/Dancoff approximation. For RHF references, solvent effects can be included by EFP1 or PCM (or both together), for both TD-DFT excitation energies and their nuclear gradients. DFTB (possibly combined with PCM) may be chosen as well, and analytic gradients for singlet and triplet are available. For SCFTYP=UHF, excited states with the same spin projection as the ground state are found. MULT in $CONTRL governs the number of alpha and beta electrons, hence Ms=(MULT-1)/2 is the only good quantum number for either the ground or excited states. Since U-TDDFT is a single excitation theory, excited states with values near Ms and near Ms+1 will appear in the calculation. There are no properties other than the excitation energy, nor gradients, nor solvent effects, at present. "limitations" for TDDFT=SPNFLP and MRSF: Spin-flip TDDFT is programmed in the "collinear approximation" which means only the HF exchange term carries a large impact on the excitation energies. Pure DFT functionals may be used, but normally hybrids with large HF exchange fractions are used. Note that spin-flip TD-DFT in the Tamm/Dancoff approximation using DFTTYP=NONE is equivalent to spin-flip CIS. In the case of SFNFLP, MULT below is ignored, since the equation of motion of SFNFLP cannot restrict a particular spin state. Therefore, the spin state of resulted response states should be determined by inspection. On the other hand, MRSF utilizes it to specify target spin state, since MRSF formulates independent equation of motion for each response spin states. Both spin-flip and MRSF codes operate only in the Tamm/Dancoff approximation, so TAMMD below is automatically .TRUE. Nuclear gradients and/or excited state properties are available only in the gas phase. Solvation effects are available for both energy and gradient calculations, for EFP1, C-PCM, or both. --------- NSTATE = Number of states to be found (excluding the reference state). The default is 1 more state. It is important to know that the first response state corresponds to the ground state in the case of TDDFT=SPNFLP or MRSF. IROOT = State used for geometry optimization and property evaluation. (default=1) TDDFT=EXCITE counts the reference as 0, and this should be the lowest state. Hence IROOT=1 means the 1st excited state, just as you might guess. TDDFT=SPNFLP labels the reference state as 0, but this might not be the lowest state overall. The meaning of IROOT=1 is the lowest state, omitting the reference state from consideration. Hence IROOT=1 might specify the ground state! MULT = Multiplicity of the excited states specifying spin state of response states, i.e., 1 (Singlet) or 3 (Triplet) in TDDFT=EXCITE and 1 (Singlet), 3 (Triplet) or 5 (Quintet) in TDDFT=MRSF (default = 1). It should be emphasized that the singlet response states by MRSF include the ground singlet state. However, this parameter is ignored when TDDFT=SPNFLP. Therefore, the spin state of resulted response states should be determined by inspection in the case of TDDFT=SPNFLP. TDPRP = a flag to request property computation for the state IROOT. Properties can only be obtained when the nuclear gradient is computable, see gradient restrictions noted in the introduction to this group. Properties require significant extra computer time, compared to the excitation energy alone, so the default is .FALSE. Properties are always evaluated during nuclear gradient runs, when they are a free by-product. It includes natural orbitals for TDDFT=MRSF. UNRLX = a flag that requests the use of the unrelaxed density in property computations for the state IROOT. UNRLX=.TRUE. also requires TDPRP=.TRUE. ALTRN = a flag that requests property computations on unrelaxed densities of all NSTATE excited states. ALTRN=.TRUE. also requires TDPRP=.TRUE. and UNRLX=.TRUE. TRNSD = a flag that requests the use of the transition density in property computations. The transition density is taken for the transition between the ground state and the state IROOT. All nuclear contributions to properties are set to zero when using the transition density. By default this flag is set to .FALSE. TRNSD=.TRUE. also requires TDPRP=.TRUE. If this flag is combined with ALTRN, properties are calculated for transition densities of the transitions between the ground state and all NSTATE excited states as well as for all state-to-state transitions between the NSTATE excited states. IFEDAT = an array of four integers that define atom ranges of donor/acceptor for FED. The format is: start of donor, end of donor, start of acceptor, end of acceptor. For the donor/acceptor definition to be valid, the ranges must be valid and non overlapping. By default IFEDAT(1)=0,0,0,0 is set. This means that no donor/acceptor definition is given. FED = a flag that requests the calculation of fragment excitation difference (FED), fragment charge difference (FCD), and Multi-FED-FCD couplings. All three calculations will be done by setting FED=.TRUE. . FED=.TRUE also requires UNRLX=.TRUE. TPA = a flag requesting two-photon absorption cross- sections. These are computed for each of the NSTATE excited states, after first evaluating their excitation energies. The TPA calculation is only available for closed shell references, only for singlet excited states (MULT=1), and may not be used with the Tamm/Dancoff approximation. Solvent effects may be treated by EFP. TAMMD is a flag selecting the Tamm/Dancoff approximation be used. This may be used with closed shell excitation energies or gradients, or open shell excitation energies. Default = .FALSE. This parameter is ignored by TDDFT=SPNFLP and MRSF, which is only coded in the Tamm/Dancoff approximation. NONEQ is a flag controlling PCM's solvent behavior: .TRUE. splits the dielectric constant into a bulk value (EPS in $PCM) and a fast component (EPSINF), see Cossi and Barone, 2001. The idea is that NONEQ=.t. is appropriate for vertical excitations, and .f. for adiabatic. (the default is .TRUE.) This keyword is ignored by TDDFT=SPNFLP. MODTD = bit-additive options 1 use fast sort of virt - occ MO energies Default: 0 MREKT = a flag calling Dyson Orbitals (DOs) of Extended Koopmans' Theorem (EKT) using MRSF-DFT Excited State Response for Ionization Energies. (default, .FALSE.) Note, that NPRINT=8 in $CONTRL prints out a final MRSF Density Matrix of a selected electronic state as well as overlaping between MRSF-DOs and initial canonical MOs if there is a $VEC block. MRDEA = a flag calling DOs of EKT using MRSF-DFT Excited State Response for Electron Affinities. (default, .FALSE.) It requares DIRSCF=.FALSE. in $SCF group. * * * ground state polarizability calculation * * * (use TDDFT=HPOL option in $CONTRL) These two frequency dependent polarizability calculations may be requested in the same run (more efficient). These properties are available only for closed shell references, require the default MULT=1 value in this input group, and may not be used with the Tamm/Dancoff approximation. Solvent effects may be treated by EFP. ALPHA = requests the polarizability. Default=.FALSE. If BETA is not chosen, give just one PFREQ. BETA = requests the hyperpolarizability. Default=.FALSE. Two values are required for PFREQ. PFREQ = an array of one or two input frequencies, omega1 and omega2, to yield the polarizability alpha(omega1;omega1) [if BETA=.F.] alpha(omega2;omega2) [if BETA=.T.] alpha(omega3;omega3) [if BETA=.T.] and/or to yield the hyperpolarizability beta(omega3;omega1,omega2). The output photon frequency is determined from omega3=-(omega1+omega2). Useful special cases second harmonic generation beta(-2W;W,W), electro-optic Pockels effect beta(-W;W,0), and optical rectification beta(0;W,-W) are among the possibilities. The default is the static polarizability and/or static hyperpolarizability: PFREQ(1)=0.0,0.0 PFREQ is given in atomic units: PFREQ=45.56/lamda, for wavelength lambda in nm. * * * Grid Selection * * * The grid type and point density used in $TDDFT may be chosen independently of the values in $DFT. Excitation energies accurate to 0.01 eV may be obtained with grids that are much sparser than those needed for the ground state, and this is reflected in the defaults. Prior to April 2008, the default grid was NRAD=24 NTHE=8 NPHI=16. NRAD = number of radial grid points in Euler-Maclaurin quadrature, used in calculations of the second or third derivatives of density functionals. (default=48) NLEB = number of angular points in the Lebedev grid. (default=110) NTHE = number of theta grid points if a polar coordinate grid is used. NPHI = number of phi grid points if a polar coordinate grid is used. NPHI should be twice NTHE. SG1 = flag selecting "standard grid one". (default=.FALSE.) See both $DFT and REFS.DOC for more information on grids. The "army grade" standard for $TDDFT is NRAD=96 combined with either NLEB=302 or NTHE=12/NPHI=24. the remaining parameters are technical in nature: CNVTOL = convergence tolerance in the iterative TD-DFT step. (default=1.0E-7) MAXVEC = the maximum number of expansion vectors used by the solver's iterations, per state (default=50). The total size of the expansion space will be NSTATE*MAXVEC. NTRIAL = the number of initial expansion vectors used. (default is the larger of 5 and NSTATE). ZVTOL = the convergence criteria of Z-Vector routine. (default is 1.0d-10. But it appears to be too tight. So 1.0d-8 is recommended.) IXCORE = an array of occupied orbitals that are included in the excitations during MRSF-TDDFT calculation. This is an option of REW (Restricted Energy Window) technique, which can limit excitations of response theory. MRSF-TDDFT utilizes a ROHF reference, which is composed of doubly occupied (C), two singly occupied (O) and virtual orbitals (V). Only the elements of IXCORE array are allowed in the excitations from C to O and V. For example, IXCORE(1)=1,10 only allows the transitions from 1 and 10th doubly occupied orbitals to O and V. All the transitions that starts from the other orbitals of C are excluded in the response space. ========================================================== ========================================================== $DFTB group (relevant for GBASIS=DFTB) Density-functional tight-binding (DFTB) is turned on by selecting GBASIS=DFTB in $BASIS. $DFTB controls optional parameters for a DFTB calculation. DFTB is formulated in a two-center approximation utilizing implicitly a minimal pseudoatomic orbital basis set with corresponding, pretabulated one- and two-center integrals. Because of this, many properties (for instances, multipoles higher than dipoles) and many options are ignored or not available in the current implementations of DFTB. DFTB also uses an independent SCF driver (SCF in DFTB is also called SCC, see below), so most SCF options are not available for DFTB. Only SCFTYP=RHF and UHF are implemented. SCFTYP=ROHF is available, only when all SPNCST values are zero. DFTB does not explicitly use symmetry (C1 throughout) since integrals are never computed during the calculations. Slater-Koster tables are only defined for spherical functions (5d) so DFTB sets ISPHER=1. Most $GUESS options do not work for DFTB (DFTB does not use initial orbitals in the usual sense). Other than the default (METHOD=HUCKEL, which is ignored), only METHOD=MOREAD works (note that SCC-DFTB can use initial charges on atoms, derived from the orbitals). RUNTYP=OPTIMIZE, HESSIAN and RAMAN are available for full (non-FMO) DFTB and FMO-DFTB. Excited state calculations for full DFTB may be performed through the standard (linear- response) time-dependent formalism (only closed shell). PCM can be used for both ground and excited state calculations, and energy and gradient can be evaluated. In DFTB calculation, the atom type is determined by its name, not its nuclear charge as elsewhere in GAMESS. The nuclear charge (the second column in $DATA) is used only in population analysis, but not in SCF. DFTB uses a notion of "species", which means an atomic type. The species are numbered according to the order in which atoms appear in $DATA. For instances, in water there are two species, O and H. An atomic type of each species needs MAXANG, which for most but not all atoms is set automatically. NDFTB order of the Taylor expansion of the total energy around a reference density in the DFTB model. = 1 NCC-DFTB, also called DFTB1. NCC stands for non-charge-consistent, i.e., no explicity charge-charge interaction term is included in the energy calculation. = 2 SCC-DFTB, also called DFTB2. SCC means a self-charge-consistent approach, and SCC implies that SCF iterations are carried out that converge monopolar charges towards self-consistency. = 3 DFTB3, including 3rd order correction using Hubbard derivatives (HUBDER). In order to reproduce the published DFTB3 approach, it is necessary to also specify DAMPXH=.TRUE. to add other terms. Gaus, M. et al. J. Chem. Theory Comput. 2011, 7, 931-948 is referred to as Gaus2011 below. Default: 2. LCDFTB = a flag to turn on long range corrected DFTB, as derived from LC-DFT. At present, only DFTB2 can be combined with this option. See EMU below. Note that LC-DFTB requires parameters optimized for it. Default: .FALSE. EMU = a long-range separation parameter for LC-DFTB (like MU in $DFT). To use LC-DFTB properly, some non-zero value (for example, 0.3) should be assigned to EMU. The value of EMU can be found in the readme file for the parameters used. Default: 0. DAMPXH = a flag to include the damping function for X-H atomic pair in DFTB3. See also DAMPEX, and eq 21 in Gaus2011. The damping function is used when at least one atom in a pair is "H". "HYDROGEN" and any other name will turn off the damping. Default: .FALSE. DAMPEX = an exponent used in the damping function for X-H atomic pairs. The default value is 4.0 (taken from the 3OB parameter set). SRSCC = a flag to perform shell-resolved SCC calculation. If set to .FALSE., the code uses the Hubbard value for an s orbital for p and d orbitals, ignoring their Hubbard values defined in Slater- Koster tables. Using .TRUE. enables the use of proper Hubbard values for p and d orbitals, implemented only for DFTB1 and DFTB2. Default: .FALSE. ITYPMX Convergence method of SCC calculations. = -1 Use standard GAMESS convergence methods. SOSCF and DIIS are supported, but DEM is not. = 0 Broyden's method. Interpolation is applied for atomic (or shell-resolved when SRSCC=.TRUE.) charges, but not Hamiltonian matrix. = 1 (reserved) = 2 DIIS for charges. Default: 0. ETEMP = electronic temperature in Kelvin. Non-zero values of ETEMP help SCF convergence of nearly-degenerate systems by smearing occupation numbers around the Fermi-level. Only the Fermi-Dirac distribution function is available as a smearing function. The default value is 0 Kelvin, meaning the smearing function is not used. ETEMP is implemented only for SCFTYP=RHF and when FMO is not used. OCCTHR = threshold for identifying doubly occupied, partially occupied or unoccupied molecular orbitals with ETEMP. If an occupation number n is OCCTHR1 (or PARALL=.TRUE.), the default becomes CODE=DDI. However, if FMO is in use, the default for closed shell parallel runs is CODE=IMS. The "SERIAL" code for OSPT=RMP will run with modest scalability when p>1. The many different MP2 programs are written for different hardware situations. Here N is the number of atomic basis functions, and O is the number of correlated orbitals in the run: The original SERIAL programs use N**3 memory, and have larger disk files and generally takes longer than CODE=IMS. The IMS program uses N*O**2 memory, and places most of its data on local disks (so you must have good disk access), and will run in parallel...ideal for small clusters. Using this program on a node where the disks are of poor quality (SATA-type) and with many cores accessing that single disk may be very I/O bound. Adding more memory can make this program run more efficiently. Network traffic is modest when running in parallel. The DDI program uses N**4 memory, but this is distributed across all nodes, and there is essentially no I/O...ideal for large parallel machines where the manufacturer has forgotten to include disk drives. MEMDDI must be given in $SYSTEM for these codes, so large problems may require many nodes to aggregate enough MEMDDI. The network traffic is high, so an Infiniband quality network or better preferred. Scalability is very good, for example, this program has been used up to 4,000 cores on Altix/ICE equipment. All of the programs just mentioned should generate the same numerical results, so select which one best matches your hardware. The RIMP2 program is an approximation to the true MP2 energy, using the "resolution of the identity" to reduce the amount of data stored (in memory and/or on disk), and also the total amount of computation. See the paper on this program for its reduced CPU and memory requirements. Network traffic is modest. The code has options within the $RIMP2 input to govern the use of replicated memory versus shared memory, as well as the use of disk storage versus distributed memory, so you can tune this to your hardware. The OMPRIMP2 program is an OpenMP threaded version of the RIMP2 code. TO use the OMPRIMP2 program you must build GAMESS with GMS_OPENMP set to true in the following files: - $GMS_DIR/install.info - $GMS_DIR/Makefile References for the various programs are given in REFS.DOC. NOSYM = disables the orbital symmetry test completely. This is not recommended, as loss of orbital symmetry is likely to mean a bad calculation. It has the same meaning as the keyword in $CONTRL, but just for the MP2 step. (Default=0) CUTOFF = transformed integral retention threshold, the default is 1.0d-9 (1.0d-12 in FMO runs). The following keyword applies only to RHF references: SCSPT = spin component scaled MP2 energy selection. = NONE - the energy will be the normal MP2 value. This is the default. = SCS - the energy used for the potential surface will be the SCS energy value. Use of SCSPT=SCS causes gradients to be those for the SCS- MP2 potential surface. For CODE=IMS, the nuclear gradient can be evaluated analytically. See NUMGRD in $CONTRL if for some reason you wish to use the other two closed shell codes for SCS-MP2 gradients. The following keywords apply to any CODE=SERIAL MP2 run, or to parallel ROHF+MP2 runs using OSPT=RMP: LMOMP2= a flag to analyze the closed shell MP2 energy in terms of localized orbitals. Any type of localized orbital may be used. This option is implemented only for RHF, and its selection forces use of the METHOD=3 transformation, in serial runs only. The default is .FALSE. CPHFBS = BASISMO solves the response equations during gradient computations in the MO basis. This is programmed only for RHF references without frozen core orbitals, when it is the default. = BASISAO solves the response equations using AO integrals, for frozen core MP2 with a RHF reference, or for ROHF or UHF based MP2. NWORD = controls memory usage. The default uses all available memory. Applies to CODE=SERIAL. (default=0) METHOD= n selects transformation method, 2 being the segmented transformation, and 3 being a more conventional two phase bin sort implementation. 3 requires more disk, but less memory. The default is to attempt method 2 first, and method 3 second. Applies only to CODE=SERIAL. AOINTS= defines AO integral storage during conventional integral transformations, during parallel runs. DUP stores duplicated AO lists on each node, and is the default for parallel computers with slow interprocessor communication, e.g. ethernet. DIST distributes the AO integral file across all nodes, and is the default for parallel computers with high speed communications. Applies only to parallel OSPT=RMP runs. ========================================================== =========================================================== $RIMP2 group (optional, relevant if CODE=RIMP2 in $MP2) This group controls the resolution of the identity MP2 program, which approximately evaluates the MP2 energy. The RI approximation greatly reduces the computer resources required, while suffering only a small error in the energies. Thus, very large atomic basis sets may be used. The input below controls both utilization of the computer resources, and the accuracy of the calculation. See also $AUXBAS, regarding the auxiliary basis set, whose choice also affects the accuracy of the calculation. The program is enabled for parallel calculation, and is tuned to today's SMP nodes. It is limited to energy calculations only, without any solvent effects, for RHF or UHF references. IAUXBF = 0 uses Cartesian Gaussians = 1 uses spherical harmonics for the auxiliary basis set used to expand the MP2 energy expression into products of 3-index matrices. The default is inherited from ISPHER. The next two control computer resources, trading memory for disk storage. GOSMP = flag requesting shared memory use. The default is .TRUE. in multi-core nodes, but .FALSE. in a uniprocessor. This option means only one copy of certain large matrices is stored per node. USEDM = a flag to store two and three center repulsion integrals in distributed memory (.TRUE.), or in disk files (.FALSE., which is the default). Selection of this flag requires MEMDDI in $SYSTEM. The default is .TRUE. The RI approximation reduces CPU time, memory requirements, and total disk storage requirements compared to exact calculation. Experimentation with these two keywords will let you tune the program to your hardware situation. For example, choosing GOSMP=.TRUE. and USEDM=.TRUE. will run without any extra disk files, while setting GOSMP=.TRUE. and USEDM .FALSE. will minimize memory usage (and network usage) at the expense of doing disk I/O. Total memory usage per node can be obtained by running EXETYP=CHECK. Note the largest replicated memory printed during the RIMP2's output, dividing by 1000000 to get the correct input for MWORDS (round up a bit). Note the largest shared memory requirement printed, also dividing by 100000, and rounding up a bit. Note the distributed memory requirement, which is already in megawords, and is the correct input for MEMDDI. Then, assuming you use p total compute process on multiple n-way nodes, the memory per node is GBytes/node= 8(n*MWORDS + shared + n*MEMDDI/p)/1024 Turning off GOSMP reduces the shared memory to 0 but increases MWORDS, which is multiplied by the number of cores per node! Turning off USEDM leads to MEMDDI=0 by using disk storage instead. If additional memory is available, increasing MWORDS can lead to a reduction in the level of the occupied orbital batch, or "LV". Larger MWORDS permits a smaller LV, which will in turn reduce the required computational time, and the required network traffic or disk I/O. The value of LV used is the last line appearing after "CHECKING SIZE OF OCCUPIED ORBITAL BATCH". The next four control numerical accuracy, but see $AUXBAS which is even more influential in regards the accuracy! OTHAUX = flag to orthogonalize the RI basis set by diagonalization of the overlap matrix. If there is reason to suspect linear dependence may exist in the RI basis, select this option to have a more numerically stable result. Larger RI basis sets such as CCT and ACCT, in particular, may benefit from selecting this. (default=.FALSE.) STOL = threshold at which to remove small overlap matrix eigenvectors, ignored if OTHAUX=.FALSE. This keyword is analogous to QMTTOL in $CONTRL for the true AO basis. (default= 1.0d-6) IVMTD = selects the procedure for removing redundancies when inverting the two-center, two-e- matrix. = 0 use Cholesky decomposition (default) = 2 use diagonalization VTOL = threshold at which to remove redundancies. This is ignored unless IVMTD=2 (default= 1.0d-6) Don't forget to see also the $AUXBAS input group! An example of this program follows. The molecule is taxol, with 1032 AOs and MOs in the 6-31G(d) basis, correlating 164 valence orbitals. The RI basis set used is SVP, which matches the true basis set in quality. There are 4175 AOs in the RI basis. The job was run on a single 8-way node (n=8, p=1,2,4,8), using MWORDS=50 (leading to LV=6), MEMDDI=580, and the largest shared memory needed is 95 million words. The total node memory is thus (8 bytes/word)*(8*50 + 95 + 8*580/ 8)/1024 = 8.4 GBytes easily fitting into a modern 16 GByte node. It reduces to (8 bytes/word)*(8*50 + 95 + 8*580/16)/1024 = 6.1 GB/node if two 8-way nodes are used. Scaling is p SCF RI-MP2 job total 1 7391 7919 15366 2 3718 4131 7860 4 1857 2290 4174 8 952 1488 2479 16 486 758 1276 using two 8-way nodes. numerical results are E(RI-MP2)= -2920.607512 versus the exact E(MP2)= -2920.606231 The 0.0013 error should be measured against the total 2nd order correlation energy, which is -8.7855, while noting the time for the 2nd order E is similar to the SCF time. =========================================================== $HPCCHEM group (optional, relevent for libcchem v2.0) This group controls the behavior of libcchem v2.0 NMODE = 0 do not use libcchem v2.0 (default) = 1 enable libcchem v2.0 ETHRSH = energy SCF convergence threshold defult value = 1.0d-06 DTHRSH = starting value for dynamic integrals screening threshold. As calculation proceeds through SCF iterations the value decreses to the minimum value (10e-10). NPRINT = print mode. Default = 1 higher numbers print less detail. See LIBCCHEM 2.0 documentation for details of print modes. SCF = use libcchem's SCF routines (default false) FOCK = use libcchem's fock build routines (default false) =========================================================== $AUXBAS group (required if CODE=RIMP2 in $MP2) This group specifies the auxiliary basis set used to define the resolution of the identity in the RI-MP2 method. The RI methods are formally exact if the RI basis set is complete, so selecting larger bases improves the results. However, this also increases the computational cost of the run! It is reasonable to use smaller RI basis sets when the AO basis is modest, and increase the RI basis when you use very large AO bases. CABNAM specifies built-in basis sets for the RI: = SVP Ahlrich's SVP basis, available H-Kr = TZVP Ahlrich's TZVP basis, available H-Ar = TZVPP Ahlrich's TZVPP basis, available H-Ar = CCD cc-pVDZ basis, available H-Ar = ACCD aug-cc-pVDZ basis, available H-Ar = CCT cc-pVTZ basis, available H-Ar = ACCT aug-cc-pVTZ basis, available H-Ar = XXXXX externally defined: see EXTCAB. CABNAM has no default, this is a required input! In FMO/AP, CABNAM(2) is used to set the auxiliary basis set for the larger basis in AP. Note IAUXBF in $RIMP2 for selecting spherical harmonics versus Cartesian Gaussians. In the RI-G3(MP2,CCSD(T) (also called the PDG method), rimp2 energy and gradient are computed. The gradient level is rimp2/6-31g(d)/aux_grd, the energy level is rimp2/g3l/aux_eng. Here, g3l is the large G3 AO basis set, aux_grd and aux_eng are auxiliary basis sets customized for gradient and energy computations, respectively. For good performance and accuracy aux_grd is usually smaller than aux_eng. To input aux_grd and aux_eng set: cabnam(1)=aux_grd,aux_eng One recommendation in our publication [Phys. Chem. A 2021, 125, 42, 9421-9429] is to use: cabnam(1)=ccd,cct EXTCAB = flag to read the basis from an external file. (default is .FALSE.) This is analogous to EXTBAS in $BASIS: no external files are provided with GAMESS. The value for CABNAM=XXXX must be 8 or fewer letters: avoid the name of any built in auxiliary basis. Your XXXX bases will be read from a file defined by environment variable EXTCAB, in the execution script. Every atom present in your molecule must be defined in the external file by a line giving its chemical symbol and then XXXX. Following this header, give the basis in free format $DATA style, containing only S, P, D, F, and G shells, terminating each atom by the usual blank line. =========================================================== ========================================================== $CCINP group (optional, relevant for any CCTYP) This group controls the coupled-cluster calculation specified by CCTYP in $CONTRL. The reference orbitals may be RHF or high-spin ROHF. If this group is not given, as is often the case, only valence electrons are correlated. This group is also necessary to execute ACP calculations. Several ground-state CCTYP choices obey at least a few of the keywords from $EOMINP, so please see that group too. Excited-state runs, such as CCTYP=EOM-CCSD, CR-EOM, or CR-EOML, and the EA, IP, DEA, and DIP-EOMCC computations read $CCINP to define orbital ranges for the underlying ground-state CCSD steps preceding EOMCC diagonalizations under $EOMINP's control. A number of CCTYP choices have been superseded by more advanced or improved approaches. For example, R-CC and CR-CC were developed prior to their CR-CCL replacement, while CR-EOML supersedes CR-EOM. CR-CCL provides an accurate approximation to the fully iterative CCSDT method, superior to the widely used CCSD(T), especially in bond breaking situations. If it is important to account for the coupling of singly, doubly, and triply excited clusters, which non-iterative corrections to CCSD neglect, or if a highly accurate representation of full CCSDT energetics is desired, it is best to use the CCT3 approach. If one is interested in methods with up to connected triple excitations, the most reasonable choices are: RHF ROHF (high spin) --- ---------------- ground states [properties]: CCSD [CCPRP] CCSD [n/a] CCSD(T) n/a CR-CCL CR-CCL CCSD3A CCSD3A CCT3 CCT3 excited states [properties]: EOM-CCSD [CCPRPE] EOM-CCSD CR-EOML n/a If one is interested in strongly correlated systems (e.g., models of metal-insulator transitions), it is worth considering the ACP options, namely, ACCSD ACC23 ACCSD3A ACCT3 which work for both RHF and ROHF references (see the ACP keyword below for further information). If one is interested in a spin-adapted and accurate description of ground and excited states of an open-shell system obtained by attaching one electron to or removing one electron from the parent closed-shell core, such as a radical, the recommended choices, which should be used with SCFTYP=RHF, are: EA-EOM3A IP-EOM3A If one is interested in a spin-adapted and accurate description of ground and excited states of an open-shell system obtained by attaching two electrons to or removing two electrons from the parent closed-shell core, such as a biradical, the recommended choices, which can work with RHF or ROHF orbitals of the target system or RHF orbitals of the underlying closed-shell core, are: DEAEOM3A DEAEOM3B DIPEOM3A although the less expensive DEAEOM2, DEAEOM2A, and DIPEOM2 options may be useful too, for example, as part of composive schemes aimed at larger systems that mix different DEA- and DIP-EOMCC levels. For a brief explanation of the above acronyms, see CCTYP. See Section 4 of this manual and papers cited in the outputs for more information. When reading CR-EOM and CR-EOML outputs, it is recommended to focus on the DELTA-CR-EOMCCSD(T) excitation energies, their rigorously size-intensive delta-CR-EOMCC(2,3) counterparts, and the delta-CR-EOMCC(2,3) total energies. Parallel computation using DDI is possible for RHF references only and only for CCTYP=CCSD or CCSD(T). Memory use in parallel runs is exotic; use EXETYP=CHECK with PARALL in $SYSTEM set on prints the per node memory requirements. See the "Further Information" section of this manual for more details about coupled-cluster runs. **** The first six keywords pertain to both RHF and ROHF **** NCORE = gives the number of frozen core orbitals to be omitted from the CC calculation. The default is the number of chemical core orbitals. NFZV = the number of frozen virtual orbitals to be omitted from the calculation. (default is 0) MAXCC = defines the maximum number of CCSD (or LCCD, CCD, CCSDt, ACCSD, ACCSDt) iterations. This parameter also applies to ROHF's left CCSD vector solver, but not RHF's left CCSD vector. See MAXCCL for RHF. (default=30 for LCCD, CCD, and CCSD, and 50 for CCSD3A, ACCSD, and ACCSD3A) ICONV = defines the convergence criterion for the cluster amplitudes as 10**(-ICONV). The ROHF reference also uses this for its left eigenstate solver, but see CVGEOM in $EOMINP for RHF references. (default is 7, but it tightens to 8 for FMO-CC.) NACTO = number of active occupied orbitals to be used in CCSDt (CCTYP=CCSD3A), CC(t;3) (CCTYP=CCT3), and selected ACP (CCTYP=ACCSD3A, ACCT3) calculations. This number corresponds to the beta active occupied orbitals, which is the same as the number of active orbitals that are doubly occupied in the RHF or ROHF reference. The alpha active occupied orbitals are automatically calculated using the formula (MULT - 1) + NACTO, which means that the singly occupied alpha orbitals in open-shell systems are automatically treated as active, even if NACTO=0 [CCSDt and CC(t;3) become CCSD and CR-CC(2,3), respectively, by setting NACTO=-1]. NACTU = number of active unoccupied orbitals to be used in CCSDt (CCTYP=CCSD3A), CC(t;3) (CCTYP=CCT3), and selected ACP (CCTYP=ACCSD3A, ACCT3) calculations. This number corresponds to the alpha active unnocupied orbitals, which is the same as the number of active orbitals which are empty in the RHF or ROHF reference. The beta active unoccupied orbitals are automatically calculated according to the formula (MULT - 1) + NACTU. SHIFT = a small positive value (usually <1.0 Hartree) added to the denominators entering the iterative amplitude equations of the CCSDt (CCTYP=CCSD3A), CC(t;3) (CCTYP=CCT3), and ACP (CCTYP=ACCSD, ACC23, ACCSD3A, and ACCT3) systems to help convergence in cases of more severe quasi-degeneracies. (default is 0) ACP = an array of five real numbers specifying the scaling factors multiplying each of the five T2^2 diagrams projected on doubles [see, e.g., Figure 1 of I. Magoulas, J. Shen, and P. Piecuch, Mol. Phys. 120, e2057365 (2022) for the numbering convention of these T2^2 diagrams]. There is no default for the ACP array. When one of the ACP options (CCTYP=ACCSD, ACC23, ACCSD3A, and ACCT3) is selected, an five-number array, excluding (1.0,1.0,1.0,1.0,1.0), must be provided. For the study of strongly correlated systems, the recommended options include ACP(1)=0.0,0.0,0.0,1.0,1.0 ACP(1)=1.0,0.0,0.0,1.0,0.0 and ACP(1)=1.0,0.0,no/(no+nu),nu/(no+nu),0.0 where no and nu are the numbers of correlated occupied and unoccupied orbitals, respectively. This keyword is ignored unless one of the ACP options is selected. **** the next group of keywords pertains to SCFTYP=RHF only **** CCPRP = a flag to select computation of the CCSD-level ground-state density matrix (see, also, CCPRPE in $EOMINP for EOMCCSD-level density matrices involving excited states). The computation of the required left eigenstates of CCSD takes extra time, so the default is .FALSE., except for CCTYP=CR-CCL or CR-EOML, where the work required for properties must be done anyway. This keyword is only available in serial runs. Notes: CCSD is the only level at which properties can be obtained. Therefore, this option can only be chosen for CCTYP=CCSD, CR-CCL, EOM-CCSD, CR-EOM, or CR-EOML. A CCSD run requesting CCPRP=.TRUE. will internally change itself to EOMCCSD to run the left CCSD, but since NSTATE of $EOMINP will still be zero, this remains a ground-state calculation. Note that the convergence criterion for left eigenstates is CVGEOM in $EOMINP, which is set to obtain excitation energies and may need tightening. Use of CCTYP=CR-EOM with CCPRP=.TRUE. or CCTYP=CR-EOML will calculate triples corrections and CCSD-level properties. There is usually little reason to select any of these: MAXCCL = maximum number of iterations for the left CCSD calculations needed to determine CCSD properties or CR-CCL energies in SCFTYP=RHF runs. This is just a synonym for MAXEOM in $EOMINP. If you want to alter the left CCSD state's convergence threshold for RHF, use CVGEOM in $EOMINP. The right CCSD state's convergence is set by MAXCC and ICONV. NWORD = a limit on memory to be used in the CC steps. The default is 0, meaning all memory available will be used. IREST = defines the restart option for CC runs using SCFTYP=RHF. If the value of IREST is greater or equal 3, program will restart from the earlier CC SCFTYP=RHF run. This option requires saving the disk file CCREST from the previous CC run. Values of IREST between 0 and 3 should not be used. The value of IREST also sets the iteration counter in the restarted run. The default is 0, meaning no restart is attempted. IREST does not apply to CCSD3A, CCT3, ACCSD, ACC23, ACCSD3A, and ACCT3 runs using SCFTYP=RHF, which are automatically redirected to SCFTYP=ROHF and MULT=1. The CCSD3A, CCT3, ACCSD, ACC23, ACCSD3A, and ACCT3 runs can be restarted using KREST and, in the case of CCT3, ACC23, and ACCT3, KREST and LREST described below, independent of SCFTYP. MXDIIS = defines the number of cluster amplitude vectors from previous iterations to be included in the DIIS extrapolation during the CCSD (or LCCD, CCD) SCFTYP=RHF run. The default value of MXDIIS is 5 for all but small problems. The DIIS solver can be disengaged by entering MXDIIS = 0. It is not necessary to change the default value of MXDIIS, unless the CC equations do not converge in spite of increasing the value of MAXCC. MXDIIS does not apply to CCSD3A, ACCSD, and ACCSD3A iterations. AMPTSH = defines a threshold for eliminating small cluster amplitudes from the CC calculations. Amplitudes with absolute values smaller than AMPTSH are set to zero. The default is to retain all small amplitudes, meaning fully accurate CC iterations. Default = 0.0. AMPTSH does not apply to CCSD3A, ACCSD, and ACCSD3A iterations. CCERI = defines whether or not a tensor decomposition technique should be applied to the two-electron repulsion integrals (2ERIs) for CC calculations. = STANDARD, no tensor decomposition, the standard 4-center 2ERIs are used (default). = RI, the resolution-of-the-identity approximation is employed and the 4-center 2ERIs are assembled as products of 3-center RI integrals on the fly. The RI approximation largely eliminates memory bottlenecks of CC calculations and offers a better parallel scaling and CPU usage. Don't forget to see also the $RICC input group for additional parameters. **** the next group of keywords pertains to SCFTYP=ROHF only **** There is little reason to select any of these. MULT = spin multiplicity to use in the reference determinant during the CC computation. The value of MULT given in the $CONTRL input determines the spin state for the ROHF orbital optimization, and is the default for the CC. It would be quite unusual to use a different spin in the SCF than in the CC. The MULT keyword in $EOMINP is of greater physical interest. IOPMET = method for the CR-CC(2,3) triples correction. = 0 means try 1 and then try 2 (default) = 1, the high memory option This option uses the most memory, but the least disk storage and the least CPU time. = 2, the high disk option This option uses least memory, by storing a large disk file. Time is slightly more than IOPMET=1, but the disk file is (NO**3 * NU**3)/6 words, where NO = number of correlated occupied spatial orbitals and NU = number of correlated unoccupied spatial orbitals) = 3, the high I/O option This option requires slightly more memory than 2, and slightly more disk than 1, but does much I/O. It is also the slowest of the three choices. Check runs will print memory needed by all three options. KREST = 0 fresh start of the CCSD or CCSDt iterations (default) = 1 restart from AMPROCC ($JOB.F70) file of a previous CCSD run, or from the $JOB.MOE file of a previous ACCSD run when CCTYP=ACCSD or ACC23, a previous CCSDt run when CCTYP=CCSD3A or CCT3, or a previous ACCSDt run when CCTYP=ACCSD3A or ACCT3. KMICRO = n performs DIIS extrapolation of the open shell CCSD, every n iterations (default is 6). Enter 0 to avoid using the DIIS converger. LREST = 0 fresh start of the left CCSD iterations (CCTYP=CCSD with CCPRP=.TRUE. or CR-CCL) or the left CCSD-like iterations using singly and doubly excited cluster amplitudes obtained in CCSDt (CCTYP=CCT3) or ACP (CCTYP=ACC23, ACCT3) calculations. (default) = 1 restart from AMPROCC ($JOB.F70) file of a previous left CCSD (CCTYP=CCSD with CCPRP=.TRUE. or CR-CCL) or left CCSD-like (CCTYP=CCT3, ACC23, and ACCT3) run. LREST=1 must be combined with KREST=1. LMICRO = n performs DIIS extrapolation of the open shell left equations every n iterations (default is 5). Enter 0 to avoid using the DIIS converger. KMICRO and LMICRO are ignored for trivial problem sizes. **** The IOSIZE input parameter applies to all CC **** calculations. IOSIZE = n enables "chunking" of direct access I/O during CC computations. The direct access (DA) record size of several scratch files generated during CC computations can become quite large, but many compilers will not permit DA record lengths in excess of 2GB. If the record size L of a DA file exceeds n, then each record will be written to/read from disk in multiple chunks, each of which is at most n bytes in length. For example, if the logical record length L = 2,400,000,000 and n = 1,000,000,000, then each logical record will be processed in 3 chunks, with respective record lengths of 10^9, 10^9, and 4^8 bytes, thereby avoiding the 2GB record length limit. = 0 disables I/O chunking (default.) An exetyp=check job with iosize=0 will indicate if chunking needs to be enabled and will also list recommended values of n. ========================================================== ========================================================== $RICC group (optional, for SCFTYP=RHF and CCTYP=CCSD/CCSD(T) in $CONTRL and essential if CCERI=RI in $CCINP) This group controls the calculation of ground state energies for closed-shell molecules (SCFTYP=RHF) using the coupled-cluster method within the singles and doubles (CCSD) truncation scheme or CCSD augmented with an approximate triples correction: CCSD(T). These calculations run in parellel using a hybrid MPI/OpenMP parallelization model and employ the resolution-of-the-identity (RI) approximation for two- electron repulsion integrals (2ERIs). The following three paremeters are related to the definition of the auxiliary basis set to be used for the RI approximation to 2ERIs. IAUXBF = 0 uses Cartesian Gaussians = 1 uses spherical harmonics for the auxiliary basis set. The default is inherited from ISPHER. CABNAM specifies built-in auxiliary basis sets: = SVP Ahlrich's SVP basis, available H-Kr = TZVP Ahlrich's TZVP basis, available H-Ar = TZVPP Ahlrich's TZVPP basis, available H-Ar = CCD cc-pVDZ basis, available H-Ar = ACCD aug-cc-pVDZ basis, available H-Ar = CCT cc-pVTZ basis, available H-Ar = ACCT aug-cc-pVTZ basis, available H-Ar = XXXXX externally defined: see EXTCAB. CABNAM has no default, this is a required input! In FMO/AP, CABNAM(2) is used to set the auxiliary basis set for the larger basis in AP. EXTCAB = flag to read the basis from an external file. (default is .FALSE.) This is analogous to EXTBAS in $BASIS: no external files are provided with GAMESS. The value for CABNAM=XXXX must be 8 or fewer letters: avoid the name of any built in auxiliary basis. Your XXXX bases will be read from a file defined by environment variable EXTCAB, in the execution script. Every atom present in your molecule must be defined in the external file by a line giving its chemical symbol and then XXXX. Following this header, give the basis in free format $DATA style, containing only S, P, D, F, and G shells, terminating each atom by the usual blank line. *** Please note that the current RI integral evaluation module in GAMESS can only deal with up to G shells. Therefore, remember to discard basis functions with higher angular momenta if they are present in your external auxiliary basis. *** The auxiliary basis sets optimized for RIMP2 calculations (extension -RI) introduce smaller errors to the RI-CCSD or RI-CCSD(T) energies compared to the RI-JK basis sets. *** The common practice is to use the same class of auxiliary basis set as the AO basis set, e.g., cc-pVDZ (AO) and cc-pVDZ-RI (i.e., GBASIS=CCD in $BASIS and CABNAM=CCD in $RICC). Please note that the chosen auxiliary basis set controls the accuracy of the RI approximation. If a higher accuracy is desired, we recommend using a larger auxiliary basis set, e.g., if GBASIS=CCD in $BASIS, then CABNAM in $RICC may be set to CCT or ACCT. Importantly, increasing the size of the auxiliary basis set only marginally increases the computational cost of CC calculations, as long as the AO basis remains the same. The following five paremeters are related to the evaluation and storage of the three-center RI integrals. OTHAUX = flag to orthogonalize the auxiliary basis set by diagonalization of the overlap matrix. If there is reason to suspect linear dependence may exist in the auxiliary basis, set OTHAUX to .TRUE. to have a more numerically stable result. Larger auxiliary basis sets such as CCT and ACCT, in particular, may benefit from selecting this. Since the CCn/ACCn auxiliary basis sets are very commonly used for coupled-cluster calculations, and they often entail linear dependence, the default is set to .TRUE. STOL = threshold at which to remove small overlap matrix eigenvectors, ignored if OTHAUX=.FALSE. This keyword is analogous to QMTTOL in $CONTRL for the true AO basis. (default= 1.0d-6) IVMTD = selects the procedure for removing redundancies when inverting the two-center, two-e- matrix. = 0 use Cholesky decomposition (default) = 2 use diagonalization VTOL = threshold at which to remove redundancies. This is ignored unless IVMTD=2 (default= 1.0d-6) USEDM3 = flag to control the storage of the three-center RI integrals in distributed memory or in the shared memory of each MPI/DDI process. = .FALSE. if stored in the shared memory of each MPI/ DDI process (default). = .TRUE. if stored in distributed memory. For this choice, the MEMDDI parameter in $SYSTEM is a mandatory input. The three-center RI integrals do not consume a lot of memory and do not present the true memory bottleneck for RI-CC calculations. This makes their storage in the shared memory of each MPI/DDI process feasible. This is also a good choice, as this reduces communication between the MPI/DDI processes. Hence, the default value of USEDM3 is set to .FALSE. For large molecules, however, choosing USEDM3=.TRUE. might be a good option for making calculations feasible. For example, if the shared memory requirement per MPI/DDI processor is nearly equal to the available memory on a single node of your computer cluster, choosing USEDM3=.TRUE. would reduce the shared memory need per MPI/DDI process to some extent. Note that this helps only a little and not a lot, since the real memory bottleneck of RI-CC calculations do not come from the three-center RI integrals. The following paremeter is related to the choice of an algorithm for RI-CCSD iterations. USEDM4 = flag to control whether or not to temporarily store blocks of preassembled four-center 2ERIs in the shared memory of the MPI/DDI proccesses during CCSD iterations. = .TRUE., needs larger memory per MPI/DDI process but the algorithm is faster (default). = .FALSE. a memory-economic algorithm but takes a longer wall time/CCSD iteration. This algorithm completely avoids the storage of four-center 2ERIs and assembles them from the 3-center RI integrals on the fly. Choosing USEDM4=.FALSE. is useful if the computational resources are limited, e.g., for computer clusters with a small number of compute nodes. In this case, selecting USEDM4=.FALSE. would allow you to perform your desired calculation at the expense of some additional wall time. *** We strongly recommend running an EXETYP=CHECK job always prior to submitting a real RI-CC calculation. This provides a detailed information on the memory requirements for the various options described above and a guideline to set the parameters MWORDS and MEMDDI in the $SYSTEM input group. The RI-CC program is OpenMP threaded. Hence, to use this program you must build GAMESS with GMS_OPENMP set to true in the following files: - $GMS_DIR/install.info - $GMS_DIR/Makefile ========================================================== ========================================================== $EOMINP group (optional, for CCTYP=EOM-CCSD, CR-EOM, or CR-EOML) (optional, for CCTYP=EA-EOMx or IP-EOMx) (optional, for CCTYP=DEAEOMx or DIPEOMx) (optional for CCSD properties, or CCTYP=CR-CCL) This group controls the calculation of excited electronic states by the equation-of-motion (EOM) coupled-cluster (CC) with single and double excitations, with optional triples corrections. It also pertains to single and double electron attachment and detachment processes, which may result in the system being left in an excited state. EOM-CCSD, EA-EOMx, IP-EOMx, DEAEOMx, and DIPEOMx can be selected for RHF and ROHF reference states, however the use of EA-EOMx and IP-EOMx for open-shell references is strongly discouraged. The remaining CCTYP options listed above, namely CR-EOM and CR-EOML, can only be used with SCFTYP=RHF. These EOM-type runs consist of an SCF calculation to obtain the appropriate reference state entering a subsequent ground-state CCSD computation (see the $CCINP group to control the ground-state calculation and the correlated orbital range) followed by EOMCCSD, IP-EOMCC, EA-EOMCC, DIP-EOMCC, or DEA-EOMCC calculations to determine the relevant excited, attached, or ionized electronic states of interest (see NSTATE below). In the case of EOMCCSD, triples corrections based on the method of moments approach may follow. This input group permits the user to select how many states are computed, however keep in mind that machine time is linear in the number of requested states. Because the default state selection is rather simplistic (only one state in the totally symmetric representation), it is usually necessary to provide inputs for this group, which includes selecting the NSTATE and IROOT keywords sensibly. Because this input group is used for several CCTYP calculations, not all keywords are used in every case and certain keywords may have slightly different meanings. In particular, a) if the reference type is RHF and CCTYP is EOM-CCSD, CR- EOM, or CR-EOML, keywords MULT and NACT are ignored. b) if the reference type is ROHF and CCTYP is EOM-CCSD, keywords CCPRPE, NACT, MTRIP, MEOM, MCI, CVGCI, MAXCI, and MICCI are ignored. MULT is sometimes ignored. c) if CCTYP is EA-EOMx or IP-EOMx (x=2,3A), keywords CCPRPE, MTRIP, MEOM, MCI, MACT, CVGCI, MAXCI, MICCI are ignored. d) if CCTYP is DEAEOMx (x=2,2A,3,3A,3B) or DIPEOMx (x=2,3,3A), keywords GROUP, IROOT, CCPRPE, MTRIP, MEOM, MCI, MOACT, MICEOM, CVGCI, NAXCI, MICCI are ignored. Additional information on CC and EOMCC methods can be found in the "Further Information" section of this manual. --- spin and spatial symmetry and state selection (relevant to CCTYP=EOM-CCSD, CR-EOM, CR-EOML, EA-EOMx, and IP-EOMx) MULT = Spin multiplicities for the target states. The meaning depends on the particular calculation. The default value for MULT is -1. If any doubly excited states or EA/IP quartets are sought, be sure to select MINIT=1 so that the initial guess routine will include these states. If the calculation is of the IP-EOM or EA-EOM type, the run will pass through an open-shell code even when the reference in $CONTRL is given as RHF. The resulting IP or EA states will be spin-adapted. The IP-EOM and EA-EOM methods can also be used with ROHF references and MULT=1 in $CONTRL and will similarly result in spin-adapted states. Although IP-EOM and EA-EOM calculations will run using genuine open-shell ROHF references (MULT>1 in $CONTRL), the resulting ionized or attached states will not be spin-adapted, so we do not recommend to do this. Therefore, assuming that $CONTRL specifies a closed-shell reference (accomplished by using either SCFTYP=RHF or ROHF with MULT=1), MULT = -1 means determine both doublet and quartet states = 2 means determine only doublet states = 4 means determine only quartet states. In case the run is EOM-CCSD with SCFTYP=ROHF and $CONTRL selects a closed-shell reference using MULT=1, the meaning of MULT in $EOMINP is MULT = -1 determine singlet, triplet, and quintet excited states, = 1 determine only singlet excited states, = 3 consider only triplet excited states, = 5 consider only quintet excited states. For ROHF-based EOMCCSD calculations using an open-shell reference (selected with MULT>1 in $CONTRL), the EOM states will not be perfectly spin adapted. The spin projection of Ms from $CONTRL's (MULT-1)/2 is the only good quantum number. The excited states will have S values near to Ms, Ms+1, or Ms+2 since single and double excitations are treated. Note that target states with spins LOWER than Ms are not generated, even if they exist in nature. The output will not print approximate S or values, but will label spatial symmetry. The MULT input will be be ignored. In case the run is RHF-based EOM-CCSD (possibly followed by the triples corrections obtained using CR-EOM or CR-EOML), the target states are singlets only. In this case, the MULT keyword in $EOMINP will be ignored. GROUP the name of the Abelian group to be used, which may only be one of the groups shown in the table below. The default is taken from $DATA and is reset to C1 if the group is non-Abelian. The purpose is to let the Abelian symmetry be turned off by setting GROUP=C1, if desired. Symmetry is used to help with the initial state selection, for controlling the EOMCC calculations, and for labeling the calculated states in the output. Spatial symmetry is not used not to speed up the calculations. NSTATE an array of up to 8 integers that indicates how many states of each symmetry type should be computed. The default is NSTATE(1)=1,0,0,0,0,0,0,0 meaning 1 totally symmetric state is to be found. The ground state is always computed, and MUST NOT be included in NSTATE's input for excited-state EOMCCSD runs. For EA or IP runs, the NSTATE input MUST include the ground state of the target species and may also include excited states. See also ISELCT below. There is no particular reason the first excited state (or ground state in EA/IP runs) should be totally symmetric, so most calculations require the user to provide a sensible NSTATE input. Up to 10 states can be found in any irrep. Keep in mind that machine time is linear in the number of states to be found, so be realistic! NSTATE uses this order for irreducible representations: irrep 1 2 3 4 5 6 7 8 C1 A C2 A B Cs A' A'' Ci Ag Au C2v A1 A2 B1 B2 C2h Ag Au Bg Bu D2 A B1 B2 B3 D2h Ag Au B1g B1u B2g B2u B3g B3u As an aside, using NSTATE(1)=0,0,0,0,0,0,0,0 for RHF references will calculate the ground state only. This can be used in CR-EOM runs to generate the type I, II, or possibly III CR-CCSD(T) energies, provided that MTRIP (described below) is given appropriately, which are not otherwise available in standard ground-state calculations. IROOT selects the state whose energy is to be saved for further calculations, such as numerical gradients, or whose properties are evaluated (see CCPRPE below). The first integer defines the irrep number, obtained from the same table as NSTATE, and the second integer specifies the state number. For example, if GROUP=C2V, then the input IROOT(1)=3,2 refers to the second B1 state. (default is IROOT(1)=1,0) IROOT's default is moderately sensible for RHF-based excited-state runs, corresponding to the ground state (labeled as state 0), which must belong to the totally symmetric representation. For ROHF-based excitations or EA/IP runs, this is not generally the case and one should select IROOT appropriately! If degenerate EOMCCSD states are detected, only one such state will be considered for subsequent triples corrections, if they are requested using CCTYP=CR-EOM or CR-EOML. The state chosen for possible triples corrections will be the state with the lower irrep number, so make sure IROOT matches this. ISELCT = an array allowing experts to reduce the number of states that are actually solved for. When given, NSTATE determines the number of states generated by the initial guess procedures and ISELECT selects those that will carry into the calculation. For example, NSTATE(1)=2,2,2,2 with ISELCT(1)=1,3,5,7 prepares two guesses in each irrep, but only iterates the EOMCCSD equations for the lowest state in each irrep (the guesses are counted serially). The next keyword addresses triples corrections to EOMCCSD. Note that non-iterative triples corrections to EOMCCSD are only available for SCFTYP=RHF. MTRIP selects the type of non-iterative triples corrections to RHF-based EOMCCSD energies (relevant only to CCTYP=CR-EOM or CR-EOML): 1 = compute the CR-EOMCCSD(T) triples corrections, termed type I and II in the output. This is the default, which skips the iterative CISD calculations needed to construct the CR-EOMCCSD(T) triples corrections of type III. 2 = after performing an additional CISD calculation, evaluate all types of the CR-EOMCCSD(T) triples corrections, including types I, II, and III. This choice of MTRIP uses approximately 50% more memory, but less CPU time than MTRIP=4. 3 = evaluate the CR-EOMCCSD(T) corrections of type III only. As with MTRIP=2, this calculation includes the iterative CISD calculation, which is needed to construct the type III triples corrections, in addition to the EOMCCSD and CR-EOMCCSD(T) calculations. 4 = carry out MTRIP=1 calculations, followed by MTRIP=3 calculations, thus evaluating all types of the CR-EOMCCSD(T) corrections (types I, II, and III in the output). As with MTRIP=2, this calculation includes the CISD iterations, which are needed to construct the type III triples corrections, in addition to the EOMCCSD and CR-EOMCCSD(T) calculations. Note: The solver used in CISD calculations needed by MTRIP=2,3,4 is selected by the keyword MCI described later on. IEOMTP = Type of the EOM-CC energy: 0 default choice (check the output) 1 CR-CC(2,3),D 2 delta-CR-CC(2,3),D The energy chosen by IEOMTP will be used as the final EOM-CC energy, affecting numerical gradients, FMO runs etc. Note that IEOMTP= 1 or 2 is meaningful only when the specified values are available (as determined by CCTYP). Default: 0 NACT pertains only to EA-EOM3A or IP-EOM3A runs (for NACT pertaining to DEAEOMx and DIPEOMx, see the subsection discussing DEA/DIP-EOMCC calculations below): NACT = the number of active MOs used to select the 3p2h or 3h2p excitations in EA-EOMCCSDt (EA-EOM3A) or IP-EOMCCSDt (IP-EOM3A) calculations. For CCTYP=EA-EOM3A, used to describe the (N+1)-electron system, NACT refers to the NACT lowest unoccupied orbitals of the N-electron reference system. For CCTYP=IP-EOM3A, used to describe the (N-1)-electron system, NACT refers to the NACT highest occupied orbitals of the N-electron reference system. The default for NACT is 0, which does not allow any 3p2h or 3h2p operators, thus yielding only EA-EOMCC2 IP-EOMCC2 results. In other words, you should input a nonzero value for NACT! CCPRPE = a flag to select computation of the EOMCCSD-level excited-state density matrices (see, also, CCPRP in $CCINP for the CCSD-level ground-state density matrix). The computation of the required left eigenstates of EOMCCSD takes extra time, so the default is .FALSE., except for CCTYP=CR-EOML, where the work required for properties must be done anyway. Notes: CCPRPE can only be used if SCFTYP=RHF and CCTYP=EOM-CCSD, CR-EOM, or CR-EOML. The property printout includes ground- and excited-state dipole moments, transition moments and oscillator strengths involving pairs of calculated states, occupation numbers for natural orbitals, and selected additional properties for the state indicated by IROOT (cf. $ELMOM). The CVGEOM convergence criterion for the left eigenstates of EOMCCSD, which is set for excitation energies, may need tightening if the calculated properties require higher precision. The CC/EOMCC density matrices are square, but not symmetric, which means that CCSD and EOMCCSD natural orbitals come in left/right pairs. To reduce the amount of output, only the right natural orbitals for the state indicated by IROOT will be found in the output file. The complete set of the calculated CCSD and EOMCCSD density matrices are printed into the PUNCH file for possible use elsewhere. Use of CCTYP=CR-EOM with CCPRPE=.TRUE. or CCTYP=CR-EOML will calculate triples corrections to CCSD and EOMCCSD as well as CCSD/EOMCCSD-level properties. --- iterative solver selection (relevant to CCTYP=EOM-CCSD, CR-EOM, CR-EOML, EA-EOMx, and IP-EOMx) MEOM selects the solver for the EOMCCSD calculations: 0 = one root at a time using a united iterative space for all calculated roots (default) 1 = one root at a time using a separate iterative space for each calculated root 2 = the Hirao-Nakatsuji multi-root solver 3 = one root at a time using a separate iterative space for each computed right/left root. (compare to 1) 4 = one root at a time using a united iterative space for all right/left roots (compare to 0). For ROHF-based EOMCCSD or IP/EA runs, there is only one iterative solver, so MEOM is ignored in these cases. MEOM=0,1,2 obtain all the right eigenvectors first, and then, if properties are being computed, proceed to compute the left eigenvectors. MEOM=3,4 obtain right and left eigenvectors simultaneously, and therefore should only be chosen if you are computing properties (see CCPRP/CCPRPE). MCI selects the solver for the CISD step, which is irrelevant unless MTRIP is greater than 1. 1 = one root at a time using a separate iterative space for each calculated root (default) 2 = the Hirao-Nakatsuji multi-root solver (slower). --- initial guess for EOMCCSD, EA/IP, and possible CISD solvers: The keywords MINIT and MACT below help select the type of initial guess used to begin EOMCCSD and EA/IP-EOMCC iterations. In the context of EOMCCSD runs, "S" and "D" stand for using all singles and doubles, while "s" and "d" mean restricting those excitations to and from a smaller number of active orbitals. Meanwhile, in the context of IP-EOMCCSD or IP-EOMCCSDt and EA-EOMCCSD or EA-EOMCCSDt calculations, "S" and "D" refer to using all 1h and 2h1p excitations in the case of IP-EOMCC and all 1p and 2p1h excitations in the case of EA-EOMCC, while "s" and "d" denote restricting these kinds of excitations to and from a smaller set of active orbitals. Of course, to define the range of orbitals "active" in the initial guess, inputs NOACT and NUACT (and perhaps MOACT) below must be given. The reason that MINIT=1 is preferred is because low-lying states with non-negligible double excitation or significant multi-configurational character are missing in simplistic CIS, 1h, or 1p guesses and, thus, may not appear in the final converged EOMCCSD or EA/IP-EOMCC calculations. MINIT selects the initial guess procedure for EOMCCSD, EA-EOMCC, IP-EOMCC, and possibly CISD iterations, in case MTRIP>1. 1 = Use EOMCCSd to start the EOMCCSD iterations, use IP-EOMCCSd to start the IP-EOMCCSD or IP-EOMCCSDt iterations, use EA-EOMCCSd to start the EA-EOMCCSD or EA-EOMCCSDt iterations, and CISd to start possible CISD iterations. 2 = Use CIS wave functions to start EOMCCSD as well as any possible CISD calculations, use 1h wave functions to start IP-EOMCC iterations, and use 1p wave functions to start EA-EOMCC iterations. MINIT's default is 2, but MINIT=1 is HIGHLY RECOMMENDED! MACT = fine tuning of MINIT's EOMCCSD initial guess. For MINIT=1 MACT=0, use EOMCCSd guess. For MINIT=1 MACT=1, use EOMCCsd guess. For MINIT=2 MACT=0, use CIS guess. For MINIT=2 MACT=1, use CIs guess. The default for MACT is 0. MACT is ignored for EA/IP-EOMCC runs. The next three keywords define active orbitals used in the initial guesses. There are no default values of NOACT and NUACT, so the user MUST provide NOACT and NUACT values if they are needed. NOACT and NUACT are usually small (5 or so), but should be chosen carefully to avoid splitting any degenerate (e.g., pi) orbital shells. NOACT the number of occupied MOs in the active space for little "s" or little "d" initial guesses. NUACT the number of unoccupied MOs in the active space for little "s" or little "d" initial guesses. MOACT an array allowing for explicit selection of the active orbitals used to define the initial guesses. If not provided, the MOACT array is filled such that the NOACT highest occupied and NUACT lowest unoccupied orbitals are selected. If MOACT is given, the number of values provided must be NOACT+NUACT. MOACT is most useful in the virtual space, where the lowest orbitals might be diffuse in nature. An example with 15 occupied orbitals, where the user has searched the virtual space looking for valence-like orbitals, might be MINIT=1 NOACT=3 NUACT=5 MOACT(1)=13,14,15, 19,20,24,25,30 ---- restarts of EOMCCSD, EA-EOMCC, and IP-EOMCC runs: JREST = 0 means no restart is used (default). = 1 means that restart data is read from AMPROCC file. One use for JREST=1 is to request additional states, with the restart taking any converged roots from disk, and doing an initial guess for additional states. You must not change MULT when restarting. --- iteration control: CVGEOM convergence criterion for the EOMCCSD, EA/IP-EOMCC, and DEA/DIP-EOMCC excitation amplitudes, which also applies to the RHF-based left EOMCCSD iterations. (default=1.0d-4). CVGEOM also serves as a convergence criterion for left CCSD calculations using RHF, in which case the default is 10**(-ICONV). MAXEOM maximum number of iterations in the EOMCCSD, EA/IP-EOMCC, and DEA/DIP-EOMCC calculations (default=50). In the case of RHF-basd EOMCCSD runs, this keyword takes on a different meaning depending on the choice of solver selected using MEOM. For MEOM=0 or 1, this is the maximum number of iterations per each calculated state. For MEOM=2, this is the maximum number of iterations for all states of the EOMCCSD multi-root procedure. MICEOM maximum number of microiterations in the EOMCCSD calculations (default=80). Rarely used. Again, for RHF-based EOMCCSD runs, the meaning of this keyword changes based on the value of MEOM. For MEOM=1 (separate iterative space for each root), this is the maximum number of microiterations for each calculated state. For MEOM=0 or 2 (united iterative space for all calculated roots), this is the maximum number of microiterations for all calculated states. It is much better to perform calculations with MICEOM > MAXEOM (i.e., remain within a single iteration cycle). If for some reason the EOMCCSD convergence is very slow and the iterative space becomes too large, it may be worth changing the default MICEOM value to MICEOM < MAXEOM to reduce the disk usage. This is not going to happen too often and normally there is no need to change the default MICEOM value. The next three keywords only apply to RHF-based CR-EOMCCSD(T) triples corrections (CCTYP=CR-EOM) if MTRIP is greater than 1: CVGCI convergence criterion for the CISD expansion coefficients (default=1.0d-4). MAXCI maximum number of iterations in the CISD calculation (default=50). For MCI=1, this is the maximum number of iterations per each calculated CISD state. For MCI=2, this is the maximum number of iterations for all states of the CISD multi-root procedure. MICCI maximum number of microiterations in the CISD calculation (default=80). Rarely used. For MCI=1 (separate iterative space for each root), this is the maximum number of microiterations for each calculated state. For MCI=2 (united iterative space for all calculated roots), this is the maximum number of microiterations for all calculated states. In analogy to MICEOM, it is much better to perform the CISD calculations with MICCI > MAXCI (i.e., in a single iteration cycle). ---- DEA/DIP-EOMCC calculations: DEA-EOMCC and DIP-EOMCC runs use the CVGEOM and MAXEOM keywords described above, in addition to the following keywords described below. The keyword MULT specifies the spin multiplicities for the targeted doubly electron attached or doubly ionized states. If quintet states are sought in the DEA-EOMCC calculations, be sure to set MINIT=1 and MACT=0 or 1. MULT = -1 determine singlet, triplet, and quintet states (default) = 1 determine only singlet states = 3 determine only triplet states = 5 determine only quintet states. NSTATE = the number of states to be computed. The default is NSTATE=1. ISELCT = an array to specify the states to be determined. The default is ISELCT(1)=1,...,NSTATE. For example, MULT=1 NSTATE=2 ISELCT(1)=1,3 will use the first and third roots of the singlet symmetry obtained in the initial guess calculation. NACT = the number of active MOs defining the active-space DEA/DIP-EOMCC methods. In the case of DEA-EOMCC calculations, NACT refers to the number of active unoccupied MOs used to select the active 3p1h excitations in DEA-EOMCC(3p1h){Nu} (CCTYP=DEAEOM2A), the active 3p1h and 4p2h excitations in DEA-EOMCC(3p1h,4p2h){Nu} (CCTYP=DEAEOM3A), or active 4p2h excitations in DEA-EOMCC(4p2h){Nu} (CCTYP=DEAEOM3B). For DIP-EOMCC methods, NACT refers to the number of active occupied MOs used to select the active 4h2p excitations in the DIP-EOMCC(4h2p){No} (CCTYP=DIPEOM3A) calculations. The default value for NACT is 0, which does not allow any 3p1h, 4p2h or 4h2p operators, and thus yields only DEA-EOMCC(2p) (when CCTYP=DEAEOM2A or DEAEOM3A), DEA-EOMCC(3p1h) (when CCTYP=DEAEOM3B), or DIP-EOMCC(3h1p) (when CCTYP=DIPEOM3A) results. MINIT and MACT select the initial guess procedure used to start the DEA-EOMCC iterations. These keywords are not applicable to DIP-EOMCC runs, which are initiated using full 2h and, if NOACT and NUACT are given, active 3h1p guesses. The allowed choices of MINIT and MACT in DEA-EOMCC runs are as follows: For MINIT=1 MACT=0, use full 2p and active 3p1h guess For MINIT=1 MACT=1, use active 2p and active 3p1h guess For MINIT=2 MACT=0, use full 2p guess (default) For MINIT=2 MACT=1, use active 2p guess. MINIT=2 MACT=0 is a default, but MINIT=1 MACT=0 with the value of NOACT and NUACT specified is HIGHLY RECOMMENDED! In order for the DEA-EOMCC calculations to find quintet states using MULT=-1 or 5, MINIT must be set at 1. NOACT = the number of active occupied MOs defining active 3h1p (CCTYP=DIPEOMx) or active 3p1h (CCTYP=DEAEOMx) initial guesses. NUACT = the number of active unoccupied MOs defining active 3h1p (CCTYP=DIPEOMx), active 2p or active 3p1h (CCTYP=DEAEOMx) initial guesses. The default values of NOACT and NUACT are 0, however, values for NOACT and NUACT must be provided unless MINIT=2 and MACT=0. The keyword JREST below controls restarts of DEA-EOMCC and DIP-EOMCC runs. The allowed values of JREST are as follows: JREST = 0 no restart, initial guess procedure will be performed (default) =-1 stop after completing the initial guess procedure = 1 restart data is read from $SCR/$JOB/$JOB.R1 file = 2 initial guess vectors are read from $SCR/$JOB/$JOB.R0 file. Initial guess vectors are written into $SCR/$JOB/$JOB.R0 file when JREST=0 or -1. The reference determinant for the DEA-EOMCC and DIP-EOMCC calculations defined in $CONTRL can be either of the RHF or the ROHF singlet or triplet type. The keyword NREF below determines whether the DEA/DIP-EOMCC calculation is carried out using the MOs of the underlying closed-shell core or the MOs of the doubly electron attached or doubly ionized target species. NREF = 0 singlet RHF MOs (SCFTYP=RHF or SCFTYP=ROHF, MULT=1 in $CONTRL) of the (N-2)-electron (DEA-EOMCC) or (N+2)-electron (DIP-EOMCC) closed-shell reference system will be used. Set ICHARGE for the closed-shell reference system in $CONTRL accordingly. (default) = 1 singlet RHF (SCFTYP=RHF or SCFTYP=ROHF, MULT=1 in $CONTRL) or triplet ROHF (SCFTYP=ROHF, MULT=3 in $CONTRL) MOs of the N-electron target system will be used. In this case, the input in $CONTRL should describe the N-electron target system. This choice of NREF requires DIRSCF=.FALSE. in $SCF group. NSAVE = 0 temporary binary $SCR/$JOB/$JOB.F* files will be deleted after all calculations are completed (default) = 1 temporary binary files are kept; this requires modifying the rungms file to keep all $SCR/$JOB/$JOB.F* files = 2 CCSD amplitudes and H-BAR will be read from existing binary files generated by previous DEA/DIP-EOMCC calculations (with NSAVE=1). All procedures before DEA/DIP-EOMCC, including SCF, CCSD, and the generation of H-BAR, will be skipped. ========================================================== ========================================================== $MOPAC group (relevant if GBASIS=PM3, AM1, or MNDO) This group affects only semi-empirical jobs, which are selected in $BASIS by keyword GBASIS. PEPTID = flag for peptide bond correction. By default a molecular mechanics-style torsion potential term is added for every peptide bond linkage found. The intent is to correct these torsions to be closer to planar than they would otherwise be in the semi-empirical model. Here, the peptide bond means any O H \\ / C----N / \ X One such torsion is added for O-C-N-H and one for O-C-N-X. This term is parameterized as in MOPAC6. Default=.TRUE. ========================================================== ========================================================== $GUESS group (optional, relevant for all SCFTYP's) This group controls the selection of initial molecular orbitals. GUESS = Selects type of initial orbital guess. = HUCKEL Carry out an extended Huckel calculation using a Huzinaga MINI basis set, and project this onto the current basis. This is implemented for atoms up to Rn, and will work for any all electron or core potential basis set. (default for most runs) = HCORE Diagonalize the one electron Hamiltonian to obtain the initial guess orbitals. This method is applicable to any basis set, but does not work as well as the HUCKEL guess. = MOREAD Read in formatted vectors punched by an earlier run. This requires a $VEC deck, and you MUST pay attention to NORB below. = RDMINI Read in a $VEC deck from a converged SCF calculation using GBASIS=MINI, to project the MINI orbitals onto the current basis. The option improves upon the Huckel guess because it involves SCF orbitals, which are typically easily obtained in the small MINI basis. This option doesn't work if the current basis uses core potentials. potentials. The $VEC from the MINI run must contain all virtual orbitals. = MOSAVED (default for restarts) The initial orbitals are read from the DICTNRY file of the earlier run. = SKIP Bypass initial orbital selection. The initial orbitals and density matrix are assumed to be in the DICTNRY file. Mostly used for RUNTYP=HESSIAN when the hessian is being read in from the input. The next options are less general, being for Fragment Molecular Orbital runs, or Divide and Conquer runs: = FMO Read orbitals from the DICTNRY file, from previous FMO run with MODPRP=1. = HUCSUB Perform a Huckel guess in each subsystem of a Divide and Conquer run = DMREAD Read a density matrix from a formatted $DM group, produced by a previous Divide and Conquer run, see NDCPRT in $DANDC. All GUESS types except 'SKIP' permit reordering of the orbitals, carry out an orthonormalization of the orbitals, and generate the correct initial density matrix, for RHF, UHF, ROHF, and GVB, but note that correct computation of the GVB density requires also CICOEF in $SCF. The density matrix cannot be generated from the orbitals alone for MP2, CI, or MCSCF, so property evaluation for these should be RUNTYP=ENERGY rather than RUNTYP=PROP using GUESS=MOREAD. PRTMO = a flag to control printing of the initial guess. (default=.FALSE.) PUNMO = a flag to control punching of the initial guess. (default=.FALSE.) MIX = rotate the alpha and beta HOMO and LUMO orbitals so as to generate inequivalent alpha and beta orbital spaces. This pertains to UHF singlets only. This may require use of NOSYM=1 in $CONTRL depending on your situation. (default=.FALSE.) NORB = The number of orbitals to be read in the $VEC group. This applies only to GUESS=MOREAD. For -RHF-, -UHF-, -ROHF-, and -GVB-, NORB defaults to the number of occupied orbitals. NORB must be given for -CI- and -MCSCF-. For -UHF-, if NORB is not given, only the occupied alpha and beta orbitals should be given, back to back. Otherwise, both alpha and beta orbitals must consist of NORB vectors. NORB may be larger than the number of occupied MOs, if you wish to read in the virtual orbitals. If NORB is less than the number of atomic orbitals, the remaining orbitals are generated as the orthogonal complement to those read. NORDER = Orbital reordering switch. = 0 No reordering (default) = 1 Reorder according to IORDER and JORDER. IORDER = Reordering instructions, giving the new molecular orbital order. This parameter applies to the common orbitals (both alpha and beta) except for UHF, where IORDER only affects the alpha MOs. Examples (let there be 10 occupied orbitals): transposition of HOMO and LUMO: IORDER(10)=11,10 a different transposition: IORDER(10)=15 IORDER(15)=10 a more general permutation: IORDER(8)=11,8,9,10 so the new orbital 10 is the original 9th. The default is IORDER(i)=i. JORDER = Reordering instructions. Same as IORDER, but for the beta MOs of UHF. INSORB = the first INSORB orbitals specified in the $VEC group will be inserted into the Huckel guess, making the guess a hybrid of HUCKEL/MOREAD. This keyword is meaningful only when GUESS=HUCKEL, and it is useful mainly for QM/MM runs where some orbitals (buffer) are frozen and need to be transferred to the initial guess vector set, see $MOFRZ. (default=0) * * * the next are 3 ways to clean up orbitals * * * PURIFY = flag to symmetrize starting orbitals. This is the most soundly based of the possible procedures. However it may fail in complicated groups when the orbitals are very unsymmetric. (default=.FALSE.) TOLZ = level below which MO coefficients will be set to zero. (default=1.0E-7) TOLE = level at which MO coefficients will be equated. This is a relative level, coefficients are set equal if one agrees in magnitude to TOLE times the other. (default=5.0E-5) SYMDEN = project the initial density in order to generate symmetric orbitals. This may be useful if the HUCKEL or HCORE guess types give orbitals of impure symmetry (?'s present). The procedure will generate a fairly high starting energy, and thus its use may not be a good idea for orbitals of the quality of MOREAD. (default=.FALSE.) ========================================================== ========================================================== $VEC group (optional, relevant for all SCFTYP's) (required if GUESS=MOREAD) This group consists of formatted vectors, as written onto file PUNCH in a previous run. It is considered good form to retain the titling comment cards punched before the $VEC card, as a reminder to yourself of the origin of the orbitals. For Morokuma decompositions, the names of this group are $VEC1, $VEC2, ... for each monomer, computed in the identical orientation as the supermolecule. For transition moment or spin-orbit coupling runs, orbitals for states one and possibly two are $VEC1 and $VEC2. ========================================================== $DM group (relevant in Divide and Conquer runs) This group consists of a formatted density matrix, read in exactly the format it was written. See GUESS=DM, and NDCPR in $DANDC. ========================================================== $MOFRZ group (optional, relevant for RHF, ROHF, GVB) This group controls freezing the molecular orbitals of your choice during the SCF procedure. If you choose this option, select DIIS in $SCF since SOSCF will not converge as well. GUESS=MOREAD is required in $GUESS. FRZ = flag which triggers MO freezing. (default=.FALSE.) IFRZ = an array of MOs in the input $VEC set which are to be frozen. There is no default for this. ========================================================== ========================================================== $STATPT group (for RUNTYP=OPTIMIZE or SADPOINT) This group controls the search for stationary points. Note that NZVAR in $CONTRL determines if the geometry search is conducted in Cartesian or internal coordinates. METHOD = optimization algorithm selection. Pick from NR Straight Newton-Raphson iterate. This will attempt to locate the nearest stationary point, which may be of any order. There is no steplength control. RUNTYP can be either OPTIMIZE or SADPOINT RFO Rational Function Optimization. This is one of the augmented Hessian techniques where the shift parameter(s) is(are) chosen by a rational function approximation to the PES. For SADPOINT searches it involves two shift parameters. If the calculated stepsize is larger than DXMAX the step is simply scaled down to size. QA Quadratic Approximation. This is another version of an augmented Hessian technique where the shift parameter is chosen such that the steplength is equal to DXMAX. It is completely equivalent to the TRIM method. (default) SCHLEGEL The quasi-NR optimizer by Schlegel. CONOPT, CONstrained OPTimization. An algorithm which can be used for locating TSs. The starting geometry MUST be a minimum! The algorithm tries to push the geometry uphill along a chosen Hessian mode (IFOLOW) by a series of optimizations on hyperspheres of increasingly larger radii. Note that there currently are no restart capabilitites for this method, not even manually. OPTTOL = gradient convergence tolerance, in Hartree/Bohr. Convergence of a geometry search requires the largest component of the gradient to be less than OPTTOL, and the root mean square gradient less than 1/3 of OPTTOL. (default=0.0001) NSTEP = maximum number of steps to take. Restart data is punched if NSTEP is exceeded. The default is 50 steps for a minimum search, but only 20 for a transition state search, which benefit from relatively frequent Hessian re-evaluations. KDIAGH = eigenvalue solver for the Hessian matrix. Possible choices are listed in $SYSTEM. The default is 2 for FMO or DC, and -1 otherwise (-1 means use KDIAG from $SYSTEM). NPRICO = Print coordinates every NPRICO optimization steps to reduce the output file. (default: 1) --- the next four control the step size --- DXMAX = initial trust radius of the step, in Bohr. For METHOD=RFO, QA, or SCHLEGEL, steps will be scaled down to this value, if necessary. (default=0.3 for OPTIMIZE and 0.2 for SADPOINT) For METHOD=NR, DXMAX is inoperative. For METHOD=CONOPT, DXMAX is the step along the previous two points to increment the hypersphere radius between constrained optimizations. (default=0.1) the next three apply only to METHOD=RFO or QA: TRUPD = a flag to allow the trust radius to change as the geometry search proceeds. (default=.TRUE.) TRMAX = maximum permissible value of the trust radius. (default=0.5 for OPTIMIZE and 0.3 for SADPOINT) TRMIN = minimum permissible value of the trust radius. (default=0.05) --- the next three control mode following --- IFOLOW = Mode selection switch, for RUNTYP=SADPOINT. For METHOD=RFO or QA, the mode along which the energy is maximized, other modes are minimized. Usually referred to as "eigenvector following". For METHOD=SCHLEGEL, the mode whose eigenvalue is (or will be made) negative. All other curvatures will be made positive. For METHOD=CONOPT, the mode along which the geometry is initially perturbed from the minima. (default is 1) In Cartesian coordinates, this variable doesn't count the six translation and rotation degrees. Note that the "modes" aren't from mass-weighting. STPT = flag to indicate whether the initial geometry is considered a stationary point. If .true. the initial geometry will be perturbed by a step along the IFOLOW normal mode with stepsize STSTEP. (default=.false.) The positive direction is taken as the one where the largest component of the Hessian mode is positive. If there are more than one largest component (symmetry), the first is taken as positive. Note that STPT=.TRUE. has little meaning with HESS=GUESS as there will be many degenerate eigenvalues. STSTEP = Stepsize for jumping off a stationary point. Using values of 0.05 or more may work better. (default=0.01) IFREEZ = array of coordinates to freeze. These may be internal or Cartesian coordinates. For example, IFREEZ(1)=1,3 freezes the two bond lengths in the $ZMAT example, which was for a triatomic $CONTRL NZVAR=3 $END $ZMAT IZMAT(1)=1,1,2, 2,1,2,3, 1,2,3 $END while optimizing the angle. If NZVAR=0, so that this value applies to the Cartesian coordinates instead, the input of IFREEZ(1)=4,8 means to freeze the x coordinate of the 2nd and y coordinate of the 3rd atom. See also IFZMAT and FVALUE in $ZMAT, and IFCART below, as IFREEZ does not apply to DLC internals. In a numerical Hessian run, IFREEZ specifies Cartesian displacements to be skipped for a Partial Hessian Analysis. IFREEZ can pertain to EFP particles, but only during RUNTYP=HESSIAN, where the 6 translational and rotational degrees of freedom of each EFP come AFTER the QM atom coordinates. For more information: J.D.Head, Int.J.Quantum Chem. 65, 827, 1997 H.Li, J.H.Jensen Theoret. Chem. Acc. 107, 211-219(2002) IFCART = array of Cartesian coordinates to freeze during a geometry optimization using delocalized internal coordinates. This probably works less well than IFREEZ when it freezes Cartesians. Only one of IFREEZ or IFCART may be chosen in a single run. IACTAT = array of "active atoms", which is a complimentary input to IFREEZ. Any atom *not* included in the list has its Cartesian coordinates frozen. Thus IACTAT(1)=3,-5,107,144,202,-211 allows 15 atoms, namely 3-5, 107, 144, and 202-211 to be optimized, while all other atoms are frozen. NZVAR in $CONTRL must be 0 when this option is chosen. IFREEZ and IACTAT are mutually exclusive. The latter acts by generating a IFREEZ for all atom coordinates not defined as "active", so users can input whichever list is shorter. --- The next two control the hessian matrix quality --- HESS = selects the initial hessian matrix. = GUESS chooses an initial guess for the hessian. (default for RUNTYP=OPTIMIZE) = READ causes the hessian to be read from a $HESS group. (default for RUNTYP=SADPOINT) = RDAB reads only the ab initio part of the hessian, and approximates the effective fragment blocks. = RDALL reads the full hessian, then converts any fragment blocks to 6x6 T+R shape. (this option is seldom used). = CALC compute the hessian, see $FORCE input. IHREP = the number of steps before the hessian is recomputed. If given as 0, the hessian will be computed only at the initial geometry if you choose HESS=CALC, and never again. If nonzero, the hessian is recalculated every IHREP steps, with the update formula used on other steps. (default=0) HSSEND = a flag to control automatic hessian evaluation at the end of a successful geometry search. (default=.FALSE.) --- the next two control the amount of output --- Let 0 mean the initial geometry, L mean the last geometry, and all mean every geometry. Let INTR mean the internuclear distance matrix. Let HESS mean the approximation to the hessian. Note that a directly calculated hessian matrix will always be punched, NPUN refers only to the updated hessians used by the quasi-Newton step. NPRT = 1 Print INTR at all, orbitals at all 0 Print INTR at all, orbitals at 0+L (default) -1 Print INTR at all, orbitals never -2 Print INTR at 0+L, orbitals never NPUN = 3 Punch all orbitals and HESS at all 2 Punch all orbitals at all 1 same as 0, plus punch HESS at all 0 Punch all orbitals at 0+L, otherwise only occupied orbitals (default) -1 Punch occ orbitals at 0+L only -2 Never punch orbitals ---- the next parameters control harmonic constraints --- Harmonic constraints can be added to the current geometry by setting ALL the keywords below. For instance, to harmonically constrain the distance between atom 3 and 12 to a distance of 2.0 Angstrom and a force constant of 500 kcal/mol, the following example can be used: IHMCON(1)=1,3,12 SHMCON(1)=2.0 FHMCON(1)=500.0 The default is all zeros which means do not do this. IHMCON = array of coordinates to constrain. The input is similar to IZMAT in $ZMAT, a code integer, and the atoms involved in the coordinate. The allowed codes are: 1 (stretches), 2 (angles), and 3 (dihedrals) SHMCON = equilibrium constraint values for the distances specified by IHMCON, given in Angstrom (bond stretches) or degrees (for angles and dihedrals) FHMCON = array of force constants for the distances specified by IHMCON, given in kcal/mol. ---- the following parameters are quite specialized ---- PURIFY = a flag to help eliminate the rotational and translational degrees of freedom from the initial hessian (and possibly initial gradient). This is much like the variable of the same name in $FORCE, and will be relevant only if internal coordinates are in use. (default=.FALSE.) PROJCT = a flag to eliminate translation and rotational degrees of freedom from Cartesian optimizations. The default is .TRUE. since this normally will reduce the number of steps, except that this variable is set false when POSITION=FIXED is used during EFP runs. ITBMAT = number of micro-iterations used to compute the step in Cartesians which corresponds to the desired step in internals. The default is 5. UPHESS = SKIP do not update Hessian (not recommended) BFGS default for OPTIMIZE using RFO or QA POWELL default for OPTIMIZE using NR or CONOPT POWELL default for SADPOINT MSP mixed Murtagh-Sargent/Powell update SCHLEGEL only choice for METHOD=SCHLEGEL ---- NNEG, RMIN, RMAX, RLIM apply only to SCHLEGEL ---- NNEG = The number of negative eigenvalues the force constant matrix should have. If necessary the smallest eigenvalues will be reversed. The default is 0 for RUNTYP=OPTIMIZE, and 1 for RUNTYP=SADPOINT. RMIN = Minimum distance threshold. Points whose root mean square distance from the current point is less than RMIN are discarded. (default=0.0015) RMAX = Maximum distance threshold. Points whose root mean square distance from the current point is greater than RMAX are discarded. (default=0.1) RLIM = Linear dependence threshold. Vectors from the current point to the previous points must not be colinear. (default=0.07) ========================================================== * * * * * * * * * * * * * * * * * * * * * See the 'further information' section for some help with OPTIMIZE and SADPOINT runs * * * * * * * * * * * * * * * * * * * * * ========================================================== $TRUDGE group (required for RUNTYP=TRUDGE) This group defines the parameters for a non-gradient optimization of exponents or the geometry. The TRUDGE package is a modified version of the same code from Michel Dupuis' HONDO 7.0 system, origially written by H.F.King. Presently the program allows for the optimization of 10 parameters. Exponent optimization works only for uncontracted primitives, without enforcing any constraints. Two non-symmetry equivalent H atoms would have their p function exponents optimized separately, and so would two symmetry equivalent atoms! A clear case of GIGO. Geometry optimization works only in HINT internal coordinates (see $CONTRL and $DATA inputs). The total energy of all types of SCF wavefunctions can be optimized, although this would be extremely stupid as gradient methods are far more efficient. The main utility is for open shell MP2 or CI geometry optimizations, which may not be done in any other way with GAMESS. If your run requires NOSYM=1 in $CONTRL, you must be sure to use only C1 symmetry in the $DATA input. OPTMIZ = a flag to select optimization of either geometry or exponents of primitive gaussian functions. = BASIS for basis set optimization. = GEOMETRY for geometry optimization (default). This means minima search only, there is no saddle point capability. NPAR = number of parameters to be optimized. IEX = defines the parameters to be optimized. If OPTMIZ=BASIS, IEX declares the serial number of the Gaussian primitives for which the exponents will be optimized. If OPTMIZ=GEOMETRY, IEX define the pointers to the HINT internal coordinates which will be optimized. (Note that not all internal coordinates have to be optimized.) The pointers to the internal coordinates are defined as: (the number of atom on the input list)*10 + (the number of internal coordinate for that atom). For each atom, the HINT internal coordinates are numbered as 1, 2, and 3 for BOND, ALPHA, and BETA, respectively. P = Defines the initial values of the parameters to be optimized. You can use this to reset values given in $DATA. If omitted, the $DATA values are used. If given here, geometric data must be in Angstroms and degrees. A complete example is a TCSCF multireference 6-31G geometry optimization for methylene, $CONTRL SCFTYP=GVB CITYP=GUGA RUNTYP=TRUDGE COORD=HINT $END $BASIS GBASIS=N31 NGAUSS=6 $END $DATA Methylene TCSCF+CISD geometry optimization Cnv 2 C 6. LC 0.00 0.0 0.00 - O K H 1. PCC 1.00 53. 0.00 + O K I $END $SCF NCO=3 NPAIR=1 $END $TRUDGE OPTMIZ=GEOMETRY NPAR=2 IEX(1)=21,22 P(1)=1.08 $END $CIDRT GROUP=C2V SOCI=.TRUE. NFZC=1 NDOC=3 NVAL=1 NEXT=-1 $END using GVB-PP(1), or TCSCF orbitals in the CI. The starting bond length is reset to 1.09, while the initial angle will be 106 (twice 53). Result after 17 steps is R=1.1283056, half-angle=51.83377, with a CI energy of -38.9407538472 Note that you may optimize the geometry for an excited CI state, just specify $GUGDIA NSTATE=5 $END $GUGDM IROOT=3 $END to find the equilibrium geometry of the third state (of five total states) of the symmetry implied by your $CIDRT. ========================================================== ========================================================== $TRURST group (optional, relevant if RUNTYP=TRUDGE) This group specifies restart parameters for TRUDGE runs and accuracy thresholds. KSTART indicates the conjugate gradient direction in which the optimization will proceed. ( default = -1 ) -1 .... indicates that this is a non-restart run. 0 .... corresponds to a restart run. FNOISE accuracy of function values. Variation smaller than FNOISE are not considered to be significant (Def. 0.0005) TOLF accuracy required of the function (Def. 0.001) TOLR accuracy required of conjugate directions (Def. 0.05) For geometry optimization, the values which give better results (closer to the ones obtained with gradient methods) are: TOLF=0.0001, TOLR=0.001, FNOISE=0.00001 ========================================================== ========================================================== $DERIV group (optional for numerical derivatives) This group sets options for numerical gradients. Do not confuse with $FORCE ($FORCE is for second derivatives). Note that numerical gradients for FMO cannot use GDDI (the only option is to use DDI) whereas non-FMO runs can use GDDI). VIBSIZ = Displacement size (in bohr). If negative, do not transform coordinates projecting out translations+rotations; use Cartesian coordiantes and no symmetry for offsets. Otherwise, if positive, project out T+R and use symmetry. For PBC, use a negative VIBSIZ. Default=0.01 NVIB = 1 or 2, the number of offsets per coordinate, similar to NVIB in $FORCE. Default=2 STRESS = a logical flag; if set, then the numerical gradient is computed with respect to cell parameters (defining a stress tensor) in PBC. Currently, only DFTB/PBC uses STRESS. Default=false MODGRN = Options for numerical gradient 0 Use the default general driver. 1 Use built-in FMO numerical gradient driver (STRESS cannot be used). Note that the use of GDDI differs in these two methods; FMO gradients can benefit from the full support of GDDI/2 in MODGRN=1, on par with any other FMO run (no, GDDI/3 is not available). In contrast, the general driver uses GDDI/2 to distribute offsets to groups. Default: 0 ========================================================== ========================================================== $FORCE group (optional, relevant for RUNTYP=HESSIAN,OPTIMIZE,SADPOINT) This group controls the computation of the hessian matrix (the energy second derivative tensor, also known as the force constant matrix), and an optional harmonic vibrational analysis. This can be a very time consuming calculation. However, given the force constant matrix, the vibrational analysis for an isotopically substituted molecule is very cheap. Related input is HESS= in $STATPT, and the $MASS, $HESS, $GRAD, $DIPDR, $VIB inputs. Calculation of the hessian automatically yields the dipole derivative tensor, giving IR frequencies. Raman intensities are obtained by following with RUNTYP=RAMAN. METHOD = chooses the computational method: = ANALYTIC is a fully analytic calculation. This is implemented for SCFTYP=RHF, UHF, ROHF, GVB (for NPAIR=0 or 1, only), and MCSCF (for CISTEP=ALDET or ORMAS, only). R-DFT and U-DFT are also analytic. This is the default for these cases. = SEMINUM does numerical differentiation of analytically computed first derivatives. This is the default for UHF, MCSCF using other CISTEPs, all solvent models, relativistic corrections, and most MP2 or CI runs. = FULLNUM numerically differentiates the energy twice, which can be used by all other cases. It requires many energies (a check run will tell how many) and so it is mainly useful for systems with only very few symmetry unique atoms. The default for METHOD is to pick ANALYTIC over SEMINUM if that is programmed, and SEMINUM otherwise. FULLNUM will never be chosen unless you specifically request it. RDHESS = a flag to read the hessian from a $HESS input, rather than computing it. This variable pertains only to RUNTYP=HESSIAN. See also HESS= in the $STATPT input group. (default is .FALSE.) PURIFY = controls cleanup Given a $ZMAT, the hessian and dipole derivative tensor can be "purified" by transforming from Cartesians to internals and back to Cartesians. This effectively zeros the frequencies of the translation and rotation "modes", along with their IR intensities. The purified quantities are punched out. Purification does change the Hessian slightly, frequencies at a stationary point can change by a wave number or so. The change is bigger at non-stationary points. (default=.FALSE. if $ZMAT is given) PRTIFC = prints the internal coordinate force constants. You MUST have provided $ZMAT input to use this. (Default=.FALSE.) HSSTYP selects the integral package used for integral derivatives in analytic hessians. This is used only for debugging purposes. = BEST use experimental routines, not fully tested. = GENERAL use slower safe code (default). --- the next four apply to numeric differentiation ---- NVIB = The number of displacements in each Cartesian direction for force field computation. This pertains only to METHOD=SEMINUM, as FULLNUM always uses double difference formulae. = 1 Move one VIBSIZ unit in each positive Cartesian direction. This requires 3N+1 evaluations of the wavefunction, energy, and gradient, where N is the number of SYMMETRY UNIQUE atoms given in $DATA. = 2 Move one VIBSIZ unit in the positive direction and one VIBSIZ unit in the negative direction. This requires 6N+1 evaluations of the wavefunction and gradient, and gives a small improvement in accuracy. In particular, the frequencies will change from NVIB=1 results by no more than 10-100 wavenumbers, and usually much less. However, the normal modes will be more nearly symmetry adapted, and the residual rotational and translational "frequencies" will be much closer to zero. (default) VIBSIZ = Displacement size (in Bohrs). This pertains to both SEMINUM and FULLNUM. The second value of VIBSIZ may be specified for Raman with PAVE, with the meaning of EFIELD in $RAMAN. Default=0.01 Let 0 mean the Vib0 geometry, and D mean all the displaced geometries NPRT = 1 Print orbitals at 0 and D = 0 Print orbitals at 0 only (default) NPUN = 2 Punch all orbitals at 0 and D = 1 Punch all orbitals at 0 and occupied orbs at D = 0 Punch all orbitals at 0 only (default) NPRHSS = an option to reduce I/O for RUNTYP=FMOHESS 1 - stop printing the Hessian and dipole derivatives 2 - stop punching the Hessian and dipole derivatives 4 - stop printing/punching normal modes 16 - in PAVE, punch full dimension Hessian default: 0 NZPHA = partition analysis for RUNTYP=FMOHESS -2 do not do it -1 define the number of normals modes to skip automatically, i.e., ignore all complex-valued frequencies. >=0 the number of low normal modes to skip (complex frequencies are always skept. For example, if there are 4 complex values and NZPHA=3 then 4 are skept; if NZPHA=6, then 6 are skept, 4 complex and 2 real). Common values are -2, -1, 0, 3 or 6. (3 is chosen to skip translational modes (usually lowest), 6 is for translational + rotational). Default: -2 JACTAT = atoms in the partition analysis for RUNTYP=FMOHESS The format of JACTAT is the same as IACTAT in STATPT. VDOSPA = partial vibrational density for RUNTYP=FMOHESS. If set to a non-zero value, this density is printed discretized by the VDOSPAR value. ----- the rest control normal coordinate analysis ---- VIBANL = flag to activate vibrational analysis. (the default is .TRUE. for RUNTYP=HESSIAN, and otherwise is .FALSE.) SCLFAC = scale factor for vibrational frequencies, used in calculating the zero point vibrational energy. Some workers correct for the usual overestimate in SCF frequencies by a factor 0.89. ZPE or other methods might employ other factors, see J.P.Merrick, D.Moran, L.Radom J.Phys.Chem.A 111, 11683-11700 (2007). The output always prints unscaled frequencies, so this value is used only during the thermochemical analysis. (Default is 1.0) TEMP = an array of up to ten temperatures at which the thermochemistry should be printed out. The default is a single temperature, 298.15 K. To use absolute zero, input 0.001 degrees. FREQ = an array of vibrational frequencies. If the frequencies are given here, the hessian matrix is not computed or read. You enter any imaginary frequencies as negative numbers, omit the zero frequencies corresponding to translation and rotation, and enter all true vibrational frequencies. Thermodynamic properties will be printed, nothing else is done by the run. PRTSCN = flag to print contribution of each vibrational mode to the entropy. (Default is .FALSE.) DECOMP = activates internal coordinate analysis. Vibrational frequencies will be decomposed into "intrinsic frequencies", by the method of J.A.Boatz and M.S.Gordon, J.Phys.Chem., 93, 1819-1826(1989). If set .TRUE., the $ZMAT input may define more than 3N-6 (3N-5) coordinates. (default=.FALSE.) PROJCT = controls the projection of the hessian matrix. The projection technique is described by W.H.Miller, N.C.Handy, J.E.Adams in J. Chem. Phys. 1980, 72, 99-112. At stationary points, the projection simply eliminates rotational and translational contaminants. At points with non-zero gradients, the projection also ensures that one of the vibrational modes will point along the gradient, so that there are a total of 7 zero frequencies. The other 3N-7 modes are constrained to be orthogonal to the gradient. Because the projection has such a large effect on the hessian, the hessian is punched both before and after projection. For the same reason, the default is .FALSE. to skip the projection, which is mainly of interest in dynamical calculations. ========================================================== There is a program ISOEFF for the calculation of kinetic and equilibrium isotope effects from the group of Piotr Paneth at the Technical University of Lodz. This program will accepts data computed by GAMESS (and other programs), and can be requested from paneth@p.lodz.pl ========================================================== $CPHF group (relevant for analytic RUNTYP=HESSIAN) This group controls the solution of the response equations, also known as coupled perturbed Hartree-Fock (CPHF). POLAR = a flag to request computation of the static polarizability, alpha. Because this property needs 3 additional response vectors, beyond those needed for the hessian, the default is to skip the property. (default = .FALSE.) CPHF = nature of integrals driving the formation of the response equation. The defaults are intelligent, so this is meant mostly for experts/debugging. = AO forms response equations from AO integrals, which usually takes less memory, and is more parallel. This is the default for RHF, UHF, ROHF, R-DFT, or U-DFT. AO-driven mode is not available for any other cases. = MO forms response equations from transformed MO integrals. This is the default for GVB or MCSCF. This is available forRHF or ROHF. = AODDI forms response equations from AO integrals, using distributed memory (see MEMDDI). This does AO integrals about 2x more than AO, but spreads the CPHF memory requirement out across multiple nodes. Coded only for RHF. SOLVER = linear equation solver choice. This is primarily a debugging option. For RHF analytic Hessians, choose from CONJG (default), DIIS, ONDISK, not all of which will work for all CPHF= choices. For imaginary frequency dependent polarizability responses (MAKEFP jobs), choose GMRES (default), biconjugate gradient stabilized BCGST, DODIIS, or an explicit solver GAUSS. Most response equations have only one solver programmed, and thus ignore this keyword. NWORD = controls memory usage for this step. The default uses all available memory. (default=0) NCPHF = maximum number of iteractions in CPHF (currently, applies to restricted methods only). Note that DFTB uses MXCPIT in $DFTB instead of NCPHF. (default=50) TOLCP = tolerance for convergence in CPHF (currently, applies to restricted methods only). Note that DFTB uses CPCONV in $DFTB instead of TOLCP. (default=5.0D-05) NSCZV = maximum numbers of iteractions for SCZV in FMO: NSCZV(1) refers to global iteractions over fragments NSCZV(2) refers to local iteractions in a fragment (default=50,50) TOSCZV = tolerances for convergence for SCZV in FMO TOSCZV(1) refers to global iteractions over fragments TOSCZV(2) refers to local iteractions in a fragment (default=1.0D-06,1.0D-08) ========================================================== ========================================================== $CPMCHF group (relevant for analytic RUNTYP=HESSIAN,NACME,CONICAL) This group controls the solution of the response equations, also called coupled perturbed multiconfiguration Hartree-Fock, for MCSCF wavefunctions. These are needed for analytic hessians (CISTEP=ALDET or ORMAS), or for state-averaged gradients or non-adiabatic coupling matrix element calculations. The full response equations are solved for hessians, while Z-vector equations are used for NACME and SA-gradients. The default converger is a (linear) conjugate gradient (CG) method, but three others may be chosen. Difficult cases might work upwards from the default CG method by: $cpmchf $end $cpmchf ipdir=50 $end $cpmchf gcro=.t. micit=5 kicit=10 $end $cpmchf gcrodr=.t. micit=10 kicit=5 $end $cpmchf gcrodr=.t. micit=30 kicit=10 reclin=.false. $end $cpmchf gcrodr=.t. micit=20 kicit=10 reclin=.false. prcchg=.true. prctol=1.0 $end $cpmchf gcro=.t. micit=50 kicit=100 prcchg=.true. prctol=1.0 $end where the last one is "the sledgehammer". The options shown in the next to last case very often work, and will be considerably faster than the last set, which should always work. GCRO will typically take many fewer iterations than CG, a measure of its robustness, but will use more machine time due to its microiterations. --- the next set apply to any CP-MCHF converger --- MAXIT = maximum iterations. (default=300) CPTOL = accuracy tolerance for cpmchf equations Ax=b. (compared to r/||b||, within orbital, CI, and state averaged components) (default=1.0D-07) PRCCHG = a flag to adjust the linear equation's preconditioning. (Default=.FALSE.) For ORMAS runs in particular, the standard preconditioner might lead to ill-conditioning with very small elements. If selected, any preconditioner elements below PRCTOL will be reset to a value of 1.0 PRCTOL = a tolerance to keep preconditioner elements, if PRCCHG is chosen. Default=1.0D-6 --- for RUNTYP=NACME --- NAPICK = a flag to select running through z-vector setup and choose which linear equations to solve. .TRUE. leads to z-vector linear equations. .FALSE. leads to non-z-vector linear equations. The z-vector equations are advantageous when the degrees of freedom exceeds the no. of electronic states involved in the state-averaging. (The defaults are set to enforce this option.) The z-vector equations are also advantageous when only a few NA couplings out of the total total possible couplings are of interest. If CISTEP=ORMAS, NAPICK=.TRUE. is forced. Selecting NAPICK=.TRUE. requires the choice of NA couplings in the NACST array (see below). (default=varies... see ROUTINE NACMEX for notes) NACST = an array that indicates which NA couplings to Calculate, if NAPICK is chosen. For example, if WSTATE in $DET contains at least four states, the NACME can be limited to state pairs 1<->2 and 3<->4 by NACST(1)=3,4, 1,2. Note that you should pick increasing order within any pair of states! The program always generates the state-specific gradient of every state with a non-zero WSTATE. (default=none) --- the next three choose the other CPMCHF convergers --- If none is selected, conjugate gradient (CG) is used. GCRODR = a flag to select Generalized Conjugate Residual with inner Orthogonalization with Deflated Restarting. (default=.FALSE.) GCRO = a flag to select Generalized Conjugate Residual with inner Orthogonalization. (default=.FALSE.) CGGCRO = a flag to alternate between GCRO and CG solvers. (default=.FALSE.) --- next apply only to GCRODR: RECLIN = a flag to select recycling of Krylov subspaces from the first linear equation. (recycle linear). The first CP-MCHF linear equation is solved and a recyclable subspace is generated. Then, after a projection of approximate solution across the subspace from the first system, the rest of the linear equations are solved. The recycled subspace may or may not give rapid convergence with fewer iterations. See MICIT and KICIT. (default=.TRUE.) --- next apply only to GCRODR,GCRO,CGGCRO: MICIT = total size of the Krylov expansion space, namely the number of micro-iterations within an overall iteration. While the MICIT variable has no limit, fifty or more micro-iterations start to become computationally unmanagable for larger systems. The default often must be increased for large systems or for geometries far from equilibrium. In addition, the GCRODR converger has a slightly modified scheme for the micro-iterations. In the first iteration, MICIT micro-iterations are performed in a GMRES(MICIT) iteration. However, for subsequent iterations, (MICIT-KICIT) micro-iterations are performed. (default=5, often increased to 20 or 30) KICIT = the size of the recyclable Krylov basis saved and modified from iteration to iteration, created from eigenvectors with small eigenvalues. If this space is too small, the run will experience ill-conditioning, but if too large, the search space includes ineffective parts. In case PRCCHG is chosen, make sure the number of vectors reset from small values does not exceed KICIT, if it does, increase KICIT. (default=5, usually small, e.g. 5-10 for GCRODR) (MICIT> KICIT for GCRODR) (KICIT>=MICIT for GCRO and CGGCRO) --- next apply only to (linear) CG: NACFAC = number of iterations before softening the current convergence tolerance by 1/3, for CG only. Use NACFAC=0 (or MAXIT) to prevent raising the from the initial CPTOL. (default=50) IPDIR = number of iterations before resetting the residual from a pseudo-residual to the true residual. This reset also resets the search directions, since keeping the old 'ill-rounded' directions is not very beneficial. If a run almost convergences but struggles in later iterations, IPDIR=50 is recommended. This option appears more useful for NACME rather than Hessian runs. (default=MAXIT) --- next apply only to CGGCRO: ITERA = the number of (linear) CG iterations to apply before alternation to the GCRO converger. (default=5) ITERB = the number of GCRO iterations to apply before alternation to the (linear) CG converger. (default=5) ========================================================== ========================================================== $MASS group (relevant for RUNTYP=HESSIAN, IRC, or DRC) This group permits isotopic substitution during the computation of mass weighted Cartesian coordinates. Of course, the masses affect the frequencies and normal modes of vibration. AMASS = An array giving the atomic masses, in amu. The default is to use the mass of the most abundant isotope. Masses through element 104 are stored. example - $MASS AMASS(3)=2.0140 $END will make the third atom in the molecule a deuterium. ========================================================== ========================================================== $HESS group (relevant for RUNTYP=HESSIAN if RDHESS=.TRUE.) (relevant for RUNTYP=IRC if FREQ,CMODE not given) (relevant for RUNTYP=OPTIMIZE,SADPOINT if HESS=READ) Formatted force constant matrix (FCM), i.e. hessian matrix. This data is punched out by a RUNTYP=HESSIAN job, in the correct format for subsequent runs. The first card in the group must be a title card. $HESS information is always punched in Cartesians. It will be transformed into internal coordinate space if a geometry search uses internals. It will be mass weighted (according to $MASS) for IRC and frequency runs. The initial FCM is updated during the course of a geometry optimization or saddle point search, and will be punched if a run exhausts its time limit. This allows restarts where the job leaves off. You may want to read this FCM back into the program for your restart, or you may prefer to regenerate a new initial hessian. In any case, this updated hessian is absolutely not suitable for frequency prediction! ========================================================== $GRAD group (relevant for RUNTYP=OPTIMIZE or SADPOINT) (relevant for RUNTYP=HESSIAN when RDHESS=.TRUE.) Formatted gradient vector at the $DATA geometry. This data is read in the same format it was punched out. For RUNTYP=HESSIAN, this information is used to determine if you are at a stationary point, and possibly for projection. If omitted, the program pretends the gradient is zero, and otherwise proceeds normally. For geometry searches, this information (if known) can be read into the program so that the first step can be taken instantly. ========================================================== ========================================================== $DIPDR group (relevant for RUNTYP=HESSIAN if RDHESS=.T.) Formatted dipole derivative tensor, punched in a previous RUNTYP=HESSIAN job. If this group is omitted, then a vibrational analysis will be unable to predict the IR intensities, but the run can otherwise proceed. ========================================================== $ALPDR group (relevant for RUNTYP=RAMAN or HESSIAN) Formatted alpha polarizability derivative tensor, punched by a previous RUNTYP=RAMAN job. If both $DIPDR and $ALPDR group are found in the input file, the applied electric field computation will be skipped, to immediately evaluate IR and Raman intensities. This restart may be most relevant for isotopic substitution. If this group is found during RUNTYP=HESSIAN, the Raman intensities will be added to the output. You might want to restart as RUNTYP=HESSIAN instead of RUNTYP=RAMAN in order to have access to PROJCT or the other options available in the $FORCE input group. ========================================================== ========================================================== $VIB group (relevant for RUNTYP=HESSIAN, METHOD=SEMINUM) Formatted restart data, consisting of energies, gradients, and dipole moments. This data is read in the same format by which is was written to the RESTART file. Just add a " $END" card, and place this group into the input file to effect a restart. If the final displacement's gradient was written as zero, delete the entire last data set (energy, gradient, and dipole). In case the numerical hessian was run in groups (see $GDDI), the $VIB entries will be out of order. The unsorted data can be read back in, if and only if the new run is also using processor groups. Note that assembling a complete $VIB could be sent to one core in one group, but use NGROUP=1 to make it look like a "group run". This group can be used to turn a less accurate single differencing run into a more accurate double differencing run (NVIB in $HESS). The mere presence of this group triggers the restart. ========================================================== $VIB2 group (relevant for hessians, METHOD=FULLNUM) (relevant for gradients, with NUMGRD=.TRUE.) Formatted restart information, consisting of energy values, as written to the RESTART file. Just add a " $END" line at the bottom, and place this group into the input file to effect a restart. This group has the same name ($VIB2), but different contents, depending on whether you are restarting a numerical gradient or a fully numerical hessian job. The mere presence of this group triggers the restart. =========================================================== ========================================================== $VSCF group (optional, relevant to RUNTYP=VSCF) This group governs the computation of vibrational frequencies including anharmonic effects. Besides the keywords shown below, the input file must contain a $HESS input (and perhaps a $DIPDR input), to start with previously obtained harmonic vibrational information. The VSCF method requires only energies, so any energy type in GAMESS may be used, perhaps with fully numerical harmonic vibrational information. Energies are sampled along the directions of the harmonic normal modes, and usually along pairs of harmonic normal modes, after which the nuclear vibrational wavefunctions are obtained. The dipole on the grid points may be used to give improved IR intensities. The most accurate calculation computes the potential surface directly, on all grid points, but this involves many energy evaluations. An attractive alternative is the Quartic Force Field approximation of Yagi et al., which computes a fit to the derivatives up to fourth order by computing a specialized set of points, after which this fit is used to generate the full grid of points for the solver. Since there are a great many independent energy evaluations, no matter which type of surface is computed, the VSCF method allows for computations in subgroups (much like the FMO method). Thus any $GDDI input group will be read and acted upon, if found. Vibrational wavefunctions are obtained at an SCF-like level, termed VSCF, using product nuclear wavefunctions, along with an MP2-like correction to the vibrational energy, which is termed correlation corrected (cc-VSCF). In addition, vibrational energy levels based on second order degenerate pertubation theory (see VDPT) or a CI analog (see VCI) may be obtained. Most VSCF applications have been carried out with an electronic structure level of MP2 with triple zeta basis sets. This is thought to give accuracy to 50 wavenumbers for the larger fundamentals. Use of internal coordinates is known to give improved accuracy for lower frequencies, particularly in weakly bound clusters. Restarts involve the $VIBSCF input (which has different formats for each PETYP), and the READV keyword. Restarts are safest on the same machine, where normal mode phases are reproducible. References for the VSCF method, the QFF approximation, and the solvers are given in Section 4 of this manual, along with a number of sample applications. * * * * * The first input variables control the generation of the potential surface on which the nuclear vibrations occur: PETYP = DIRECT computes the full potential energy surface, according to NCOUP/NGRID. The total number of energy/dipole calculations for NCOUP=2 will be M*NGRID + (M*(M-1)/2)*NGRID*NGRID, where M is the number of normal modes. This is the default. = QFF the Quartic Force Field approximation to the potential surface is obtained. This is usually only slightly less accurate, but has a greatly reduced computational burden, namely 6*M + 12*M*(M-1)/2 energy/dipoles. INTCRD = flag setting the coordinate system used for the grids. Any internal coordinates to be used must be defined in $ZMAT, using 3N-6 simple, DLC, or natural internal coordinates. Of course, you must enter NZVAR in $CONTRL as well. The default is to use Cartesians (default .FALSE.) INTTYP = 0 default if INTCRD=.FALSE. (ignore this keyword) = 1 implies that the $ZMAT contains only stretches, bends, and torsions. It also selects an approximate transformation between Cartesian and internal coords. = 2 the other $ZMAT coordinates may be used, and the coordinate transformation will be iterated to convergence. (default if INTCRD=.TRUE.) NCOUP = the order of mode couplings included. = 1 computes 1-D grids along each harmonic mode = 2 adds additionally, 2-D grids along each pair of normal modes. (default=2) = 3 adds additionally, 3-D grids for mode triples, for PETYP=DIRECT only. NGRID = number of grid points to be used in solving for the anharmonic vibrational levels. In the case of PETYP=DIRECT, each of these grid points must be explicitly computed. For PETYP=QFF these grid points are obtained from a fitted quartic force field. Reasonable values are 8 or 16 for DIRECT, with 16 considered significantly more accurate. For PETYP=QFF, the generation of the solver grid is very fast, so use 16 always. (default=16) AMP = step size for PETYP=DIRECT displacements. The maximum distance along each mode is a function of its frequency, amplitude(i)=sqrt(2*(AMP+1/2)/freq(i)) so that AMP resembles a vibrational quantum number. The default goes far enough past the classical turning points of the fundamentals to capture the relevant part of the surface. (default = 7.0) STPSZ = step size for PETYP=QFF displacements. The step along each mode depends on the harmonic frequency, as well as this parameter, whose default is usually satisfactory (default=0.5) In case the user wants to control each normal mode with a separate parameter, arrays of values may be given, using the keywords AMPX(1)=xx,yy,... or STPSZX(1)=xx,yy,zz... IMODE = array of modes for which anharmonic effects will be computed. IMODE(1)=10,19 computes anharmonic energies and wavefunctions for modes 10 and 19, only. In the current implementation, pairs of modes cannot be coupled, so NCOUP is forced to 1 if this option is specified. This approximation is intended for larger molecules, where the whole VSCF calculation is prohibitive. * * * * * The next set of keywords relates to the solver step which finds the vibrational states. The results always include VSCF and cc-VSCF (SCF and non-degenerate MP2-like solutions). Use of the restart option makes comparing the solvers very fast, compared to the time to generate the electronic potential energy surface's points. VDPT = option to use 2nd order degenerate perturbation theory, based on the ground and singly excited vibrational levels. Results for virtual CI within the same singly excited space will also be given. Selection of VDPT turns VCI on, as well. (default=.FALSE.) VCI = option to use the virtual CI solver within a space of the ground and both singly and doubly excited vibrational levels. Selection of VCI turns VDPT off. (default=.FALSE.) The solver always finds the ground vibrational state (v=0) by default, and defaults to finding the fundamentals (v=1 in every mode). It can rapidly find excited levels (such as all v=2) if restarted (see READV) from $VIBSCF, using the following to control the excitation levels: IEXC = 1 obtain fundamental frequencies (default) = 2 instead, obtain first overtones = 3 instead, obtain second overtones IEXC2 = 0 skip combination bands (default) = 1 add one additional quanta in other modes = 2 add two other quanta in one mode at a time. IEXC IEXC2 for H2O, which has only three modes: 0 0 only 000 ground state, no transitions 1 0 000, and 100, 010, 001 (fundamentals) 2 0 000, and 200, 020, 002 (1st overtones) 3 0 000, and 300, 030, 003 (2nd overtones) 1 1 000, and 100, 010, 001, 110, 101, 110 (1st overtones and combinations) 1 2 000, and 100, 010, 001, 210, 201, 021 2 1 000, and 200, 020, 002, 120, 102, 012 between them, 1st and 2nd overtones, and all 2-1-0 combinations. ICAS1, ICAS2 = starting and ending vibrations whose quanta are included. The default is all modes, ICAS1=1 and ICAS2=3N-6 (or 3N-5). SFACT = a numerical cutoff for small contributions in the solver. The default is 1d-4: 5d-3 or 1d-3 may affect accuracy of results, 1d-4 is safer, and 1d-5 might not converge. VCFCT = scaling factor for pair-coupling potential. Sometimes when pair-coupling potential values are larger than the corresponding single mode values, they must be scaled down. It is seldom necessary to select a scaling other than unity. (Default=1.0) * * * * * The next two relate to simplified intensity computation. These simplifications are aimed at speeding up MP2 runs, if one does not care so much about intensities, and would like to eliminate the considerable extra time to compute MP2- level dipoles. DMDR must not be used if overtones are being computed. DMDR = if true, indicates that the harmonic dipole derivative tensor $DIPDR will be read and used, rather than computing dipoles. (default=.FALSE.) MPDIP = If .TRUE. the run will compute MP2 level dipoles for the IR intensity evaluation. Entering .FALSE. uses SCF level dipoles instead. Default=.TRUE. for MP2 runs, except when using the RI-MP2 program, which cannot compute MP2 dipoles, and so chooses .FALSE. here. It is more accurate to use the DMDR flag instead instead of turning off MPDIP, if an MP2 level $DIPDR is available from the MP2 hessian run. * * * * These relate to the initial harmonic mode generation. Normally, a $HESS is provided, from which harmonic modes are obtained. It is possible to give the harmonic data explicitly with the first two: RDFRQ = array of harmonic frequencies, starting from the smallest. CMODE = array of normal mode displacements given in the same order as the frequencies read in RDFRQ. The data should be the x,y,z displacement of the first atom of the first mode, then x,y,z for the second atom, then going on to give each additional mode. PROJCT = controls the projection of the hessian matrix (same meaning as in $FORCE). Default is .TRUE. which removes small mixings between rotations or translations and the harmonic modes. * * * * READV = flag to indicate restart data $VIBSCF should be read in to resume an interrupted calculation, or to obtain overtones in follow-on runs. (default is .FALSE.) GEONLY = option to generate all points on the potential energy surface needed by the VSCF routine, without energy evaluations. The purpose of this is to prepare a set of geometries at which the energy is needed. A possible use for this is to obtain energies from a different program package, which might have an energy unavailable in GAMESS, but which lacks its own VSCF program. (default=.false.) ========================================================== $VIBSCF group (optional, relevant to RUNTYP=VSCF) This is restart data, as written to the disk file RESTART in a complete or partially completed previous run. Append a " $END", and also select READV=.TRUE. to read the data. $VIBSCF's contents are different for PETYP=DIRECT or QFF. The format of this group changed in December 2006, so that old groups can no longer be used. ========================================================== ========================================================== $GAMMA group required if RUNTYP=GAMMA This group governs evaluation of the 3rd derivative of the energy with respect to nuclear coordinates, by finite differentiation of Hessians (see $FORCE options). NFCM = n describes the amount of restart data provided. The default is n=-1, to evaluate everything. A value of n means that n+1 $FCM groups are to be read from the file (hessian #0 means the equilibrium geometry). Restart data is read from a .gamma file, created by an earlier run. DELTA = step size, default=0.01 Bohr PRTALL = flag to print full Hessian and Gamma matrix, the default is .FALSE. PRTSYM = flag to print unsymmetrical Gamma elements, the default is .FALSE. PRTBIG = flag to print large Gamma elements, default = .F. ========================================================== $EQGEOM group required if NFFLVL=2 or 3 in $CONTRL The coordinates of the stationary point, where the hessian and possibly 3rd derivative information was evaluated, in exactly the format it was printed by RUNTYP=GAMMA. ========================================================== $HLOWT group required if NFFLVL=2 or 3 in $CONTRL $GLOWT group required if NFFLVL=3 in $CONTRL These are the lower triangular parts of the hessian and 3rd derivative matrices, read in the same format as printed by an earlier RUNTYP=GAMMA. ========================================================== ========================================================== $IRC group (relevant for RUNTYP=IRC) This group governs the location of the intrinsic reaction coordinate (also called the minimum energy path, MEP), a steepest descent path in mass weighted coordinates, that connects the saddle point to reactants and products. The IRC serves a proof of the mechanism for a reaction, and is a starting point for reaction path dynamics. The IRC may be found for systems with QM atoms, EFP particles, or the combinations of QM and EFP particles, or QM plus the optional SIMOMM plug-in MM atoms. Restart data for RUNTYP=IRC is written into the PUNCH file. Information summarizing the reaction path is written to the TRAJECT file, which should be saved, appending these as various restarts are done. The graphics program MacMolPlt can display a movie of the entire mechanism, if you join the entire forward and entire backwards trajectory files, while changing the path distance parameter in the reverse part to a negative value. ----- there are five integration methods chosen by PACE. PACE = GS2 selects the Gonzalez-Schlegel second order method. This is the default method. Related input is: GCUT cutoff for the norm of the mass-weighted gradient tangent (the default is chosen in the range from 0.00005 to 0.00020, depending on the value for STRIDE chosen below. RCUT cutoff for Cartesian RMS displacement vector. (the default is chosen in the range 0.0005 to 0.0020 Bohr, depending on the value for STRIDE) ACUT maximum angle from end points for linear interpolation (default=5 degrees) MXOPT maximum number of constrained optimization steps for each IRC point (default=20) IHUPD is the hessian update formula. 1 means Powell, 2 means BFGS (default=2) GA is a gradient from the previous IRC point, and is used when restarting. OPTTOL is a gradient cutoff used to determine if the IRC is approaching a minimum. It has the same meaning as the variable in $STATPT. (default=0.0001) PACE = LINEAR selects linear gradient following (Euler's method). Related input is: STABLZ switches on Ishida/Morokuma/Komornicki reaction path stabilization. The default is .TRUE. DELTA initial step size along the unit bisector, if STABLZ is on. Default=0.025 Bohr. ELBOW is the collinearity threshold above which the stabilization is skipped. If the mass weighted gradients at QB and QC are almost collinear, the reaction path is deemed to be curving very little, and stabilization isn't needed. The default is 175.0 degrees. To always perform stabilization, input 180.0. READQB,EB,GBNORM,GB are energy and gradient data already known at the current IRC point. If it happens that a run with STABLZ on decides to skip stabilization because of ELBOW, this data will be punched to speed the restart. PACE = QUAD selects quadratic gradient following. Related input is: SAB distance to previous point on the IRC. GA gradient vector at that historical point. PACE = AMPC4 selects the fourth order Adams-Moulton variable step predictor-corrector. Related input is: GA0,GA1,GA2 which are gradients at previous points. PACE = RK4 selects the 4th order Runge-Kutta variable step method. There is no related input. ----- The next two are used by all PACE choices ----- STRIDE = Determines how far apart points on the reaction path will be. STRIDE is used to calculate the step taken, according to the PACE you choose. The default is good for the GS2 method, which is very robust. Other methods should request much smaller step sizes, such as 0.10 or even 0.05. (default = 0.30 sqrt(amu)-Bohr) NPOINT = The number of IRC points to be located in this run. The default is to find only the next point. (default = 1) ----- constraint ----- Of course, applying a constraint to the saddle point search and the reaction path means that you are not locating the true saddle, nor following the true reaction path. IFREEZ = array of Cartesian coordinates to freeze. The IRC stepper works in mass-weighted Cartesian space, making it impossible to freeze internal coordinates. An input of IFREEZ(1)=4,8 means to freeze the x coordinate of the 2nd atom and the y coordinate of the 3rd atom, that is, we count coordinates x1,y1,z1,x2,y2,z2,x3,y3,z3,... ----- The next two let you choose your output volume ----- Let F mean the first IRC point found in this run, and L mean the final IRC point of this run. Let INTR mean the internuclear distance matrix. NPRT = 1 Print INTR at all, orbitals at all IRC points 0 Print INTR at all, orbitals at F+L (default) -1 Print INTR at all, orbitals never -2 Print INTR at F+L, orbitals never NPUN = 1 Punch all orbitals at all IRC points 0 Punch all orbitals at F+L, only occupied orbitals at IRC points between (default) -1 Punch all orbitals at F+L only -2 Never punch orbitals ----- The next two tally the reaction path results. The defaults are appropriate for starting from a saddle point, restart values are automatically punched out. NEXTPT = The number of the next point to be computed. STOTAL = Total distance along the reaction path to next IRC point, in mass weighted Cartesian space. ----- The following controls jumping off the saddle point. If you give $HESS input, FREQ and CMODE will be generated automatically. SADDLE = A logical variable telling if the coordinates given in the $DATA deck are at a saddle point (.TRUE.) or some other point lying on the IRC (.FALSE.). If SADDLE is true, either a $HESS group or else FREQ and CMODE must be given. (default = .FALSE.) Related input is: TSENGY = A logical variable controlling whether the energy and wavefunction are evaluated at the transition state coordinates given in $DATA. Since you already know the energy from the transition state search and force field runs, the default is .F. FORWRD = A logical variable controlling the direction to proceed away from a saddle point. The forward direction is defined as the direction in which the largest magnitude component of the imaginary normal mode is positive. (default =.TRUE.) EVIB = Desired decrease in energy when following the imaginary normal mode away from a saddle point. (default=0.0005 Hartree) FREQ = The magnitude of the imaginary frequency, given in cm**-1. CMODE = An array of the components of the normal mode whose frequency is imaginary, in Cartesian coordinates. Be careful with the signs! You must give FREQ and CMODE if you don't give a $HESS group, when SADDLE=.TRUE. The option of giving these two variables instead of a $HESS does not apply to the GS2 method, which must have a hessian input, even for restarts. Note also that EVIB is ignored by GS2 runs. * * * * * * * * * * * * * * * * * * For hints about IRC tracking, see the 'further information' section. * * * * * * * * * * * * * * * * * * ========================================================== ========================================================== $DRC group (relevant for RUNTYP=DRC) This group governs "direct dynamics", following the dynamical reaction coordinate, which is a classical trajectory based on quantum chemistry potential energy surfaces. These may be either ab initio or semi-empirical, and are computed "on the fly" as the trajectory proceeds. Because the vibrational period of a normal mode with frequency 500 wavenumbers is 67 fs, a DRC needs to run for many steps in order to sample a representative portion of phase space. Restart data can be found in the job's OUTPUT file, with important results summarized to the TRAJECT file. Almost all DRCs break molecular symmetry, so build your molecule with C1 symmetry in $DATA, or specify NOSYM=1 in $CONTRL. RUNTYP=DRC may not be used with EFP particles. NSTEP = The number of DRC points to be calculated, not including the initial point. (default = 1000) DELTAT = is the time step. (default = 0.1 fs) TOTIME = total duration of the DRC computed in a previous job, in fs. The default is the correct value when initiating a DRC. (default=0.0 fs) * * * In general, a DRC can be initiated anywhere, so $DATA might contain coordinates of the equilibrium geometry, or a nearby transition state, or something else. You must also supply an initial kinetic energy, and the direction of the initial velocity, for which there are a number of options: EKIN = The initial kinetic energy (default = 0.0 kcal/mol) See also ENM, NVEL, and VIBLVL regarding alternate ways to specify the initial value. VEL = an array of velocity components, in Bohr/fs. When NVEL is false, this is simply the direction of the velocity vector. Its magnitude will be automatically adjusted to match the desired initial kinetic energy, and it will be projected so that the translation of the center of mass is removed. Give in the order vx1, vy1, vz1, vx2, vy2, ... NVEL = a flag to compute the initial kinetic energy from the input VEL using the sum of mass*VEL*VEL/2. This flag is usually selected only for restarts. (default=.FALSE.) The next three allow the kinetic energy to be partitioned over all normal modes. The coordinates in $DATA are likely to be from a stationary point! You must also supply $HESS input, which is the nuclear force constant matrix at the starting geometry. VIBLVL = a flag to turn this option on (default=.FALSE.) VIBENG = an array of energies (in units of multiples of the hv of each mode) to be imparted along each normal mode. The default is to assign the zero point energy only, VIBENG(1)=0.5, 0.5, ..., 0.5 when HESS=MIN, and 0.0, 0.5, ..., 0.5 if HESS=TS. If given as a negative number, the initial direction of the velocity vector is along the reverse direction of the mode. "Reverse" means the phase of the normal mode is chosen such that the largest magnitude component is a negative value. An example might be VIBENG(4)=2.5 to add two quanta to mode 4, along with zero point energy in all modes. RCENG = reaction coordinate energy, in kcal/mol. This is the initial kinetic energy given to the imaginary frequency normal mode when HESS=TS. If this is given as a negative value, the direction of the velocity vector will be the "reverse direction", meaning the phase of the normal mode will be chosen so its largest component is negative. * * * The next two pertain to initiating the DRC along a single normal mode of vibration. No kinetic energy is assigned to the other modes. You must also supply $HESS input for the initial geometry. NNM = The number of the normal mode to which the initial kinetic energy is given. The absolute value of NNM must be in the range 1, 2, ..., 3N-6. If NNM is a positive/negative value, the initial velocity will lie in the forward/reverse direction of the mode. "Forward" means the largest normal mode component is a positive value. (default=0) ENM = the initial kinetic energy given to mode NNM, in units of vibrational quanta hv, so the amount depends on mode NNM's vibrational frequency, v. If you prefer to impart an arbitrary initial kinetic energy to mode NNM, specify EKIN instead. (default = 0.0 quanta) To summarize, there are 5 ways to initiate a trajectory: 1. VEL vector with NVEL=.TRUE. This is difficult to specify at your initial point, and so this option is mainly used when restarting your trajectory. The restart information is always in this format. 2. VEL vector and EKIN with NVEL=.FALSE. This will give a desired amount of kinetic energy in the direction of the velocity vector. 3. VIBLVL and VIBENG and possibly RCENG, to give some initial kinetic energy to all normal modes. 4. NNM and ENM to give quanta to a single normal mode. 5. NNM and EKIN to give arbitrary kinetic energy to a single normal mode. * * * The most common use of the next two is to analyze a trajectory with respect to the normal modes of a minimum energy geometry it travels around. NMANAL = a flag to select mapping of the mass-weighted Cartesian DRC coordinates and velocity (conjugate momentum) in terms of normal modes at a nearby reference stationary point (which can be either a minimum or transition state). This reference geometry could in fact be the same as the initial point of the DRC, but does not need to be. If you choose this option, you must supply C0, HESS2, and $HESS2 input corresponding to the reference stationary point. (default=.FALSE.) C0 = an array of the coordinates of the stationary reference point (the coordinates in $DATA might well be some other coordinates). Give in the order x1,y1,z1,x2,y2,... in Angstroms. * * * The next options apply to input choices which may read a $HESS at the initial DRC point, namely NNM or VIBLVL, or to those that read a $HESS2 at some reference geometry (NMANAL). HESS = MIN indicates the hessian supplied for the initial geometry corresponds to a minimum (default). = TS indicates the hessian is for a saddle point. HESS2 = MIN (default) or TS, the same meaning, for the reference geometry. These are used to decide if modes 1-6 (minimum) or modes 2-7 (TS) are to be excluded from the hessian as the translational and rotational contaminants. If the initial and reference geometries are the same, these two hessians will be duplicates of each other. The next variables can cause termination of a run, if molecular fragments get too far apart or close together. NFRGPR = Number of atom pairs whose distance will be checked. (default is 0) IFRGPR = Array of the atom pairs. 2 times NFRGPR values. FRGCUT = Array for a boundary distance (in Bohr) for atom pairs to end DRC calculations. The run will stop if any distance exceeds the tolerance, or if a value is given as a negative number, if the distance becomes shorter than the absolute value. In case the trajectory starts outside the bounds specified, they do not apply until after the trajectory reaches a point where the criteria are satisfied, and then goes outside again. Give NFRGPR values. * * * The final variables control the volume of output. Let F mean the first DRC point found in this run, and L mean the last DRC point of this run. NPRTSM = summarize the DRC results every NPRTSM steps, to the TRAJECT file. (default = 1) NPRT = 1 Print orbitals at all DRC points 0 Print orbitals at F+L (default) -1 Never print orbitals NPUN = 2 Punch all orbitals at all DRC points 1 Punch all orbitals at F+L, and occupied orbitals at DRC points between 0 Punch all orbitals at F+L only (default) -1 Never punch orbitals ========================================================== =========================================================== $MEX group (relevant if RUNTYP=MEX) This group governs a search for the lowest energy on the 3N-7 dimensional "seam" of intersection of two different electronic potential energy surfaces. Such Minimum Energy Crossing Points are important for processes such as spin- orbit coupling that involve transfer from one surface to another, and thus are analogous to transition states on a single surface. The present program requires that the two surfaces differ in spin quantum number, or space symmetry, or both. Analytic gradients are used in the search. In case the two potential surfaces have identical spin and space symmetry, this kind of intersection point is referred to as a Conical Intersection. See $CONICL using RUNTYP=CONICAL instead. SCF1, SCF2 = define the molecular wavefunction types, possibly in conjunction with the usual MPLEVL and DFTTYP keywords. MULT1, MULT2 = give the spin multiplicity of the states. Permissible combinations of wavefunctions are RHF with ROHF/UHF ROHF with ROHF UHF with UHF as well as their MP2 and DFT counterparts, and GVB with ROHF/UHF MCSCF with MCSCF (CISTEP=ALDET or GUGA only) NSTEP = maximum number of search steps (default=50) STPSZ = Step size during the search (default = 0.1D+00) NRDMOS = Initial orbitals can be read in = 0 No initial orbitals (default) = 1 Read in orbitals for first state (in $VEC1) = 2 Read in orbitals for second state (in $VEC2) = 3 Read in orbitals for both ($VEC1 and $VEC2) NMOS1 = Number of orbitals for first state's $VEC1. NMOS2 = Number of orbitals for second state's $VEC2. NPRT = Printing orbitals = 0 No orbital printed out except at the first geometry (default) = 1 Orbitals are printed each geometry. If MCSCF is used, CI expansions are also printed. Finer control of the convergence criterion: TDE = energy difference between two states (default = 1.0D-05) TDXMAX = maximum displacement of coordinates (default = 2.0D-03) TDXRMS = root mean square displacement (default = 1.5D-03) TGMAX = maximum of effective gradient between the two states (default = 5.0D-04) TGRMS = root mean square effective gradient tolerance (default = 3.0D-04) =========================================================== Usage notes: 1. Normally $CONTRL will not give SCFTYP or MULT keywords. SCF1 and SCF2 can be given in any order. The combinations permitted ensure roughly equal sophistication in the treatment of electron correlation. 2. After reading $MEX, SCFTYP and MULT will be set to the more complex of the two choices, which is considered to be RHF < ROHF < UHF < GVB < MCSCF. This permits the $SCF input defining a GVB wavefunction to be read and tested for correctness, in a GVB+ROHF run. Since only one SCFTYP is stored while reading the input, you might need to provide some keywords that are normally set by default for the other (such as ensuring DIIS is selected in $SCF if either of the states is UHF). 3. It is safest by far to prepare and read $VEC1 and $VEC2 groups so that you know what electronic states you start with. It is a good idea to regenerate both states at the end of the MEX search, to be sure that they remain as you began. 4. It is your responsibility to make sure that the states have a different space symmetry, or a different spin symmetry (or both). That is why note 3 is so important. 5. $GRAD1 and/or $GRAD2 groups containing gradients may be given to speed up the first geometry of the MEX search. 6. The search is even trickier than a saddle point search, for it involves the peaks and valleys of BOTH surfaces being generated. Starting geometries may be guessed as lying between the minima of the two surfaces, but the lowest energy on the crossing seam may turn out to be somewhere else. Be prepared to restart! 7. The procedure is a Newton-Raphson search, conducted in Cartesian coordinates, with a Lagrange multiplier imposing the constraint of equal energy upon the two states. The hessian matrices in the search are guessed at, and subjected to BFGS updates. Internal coordinates will be printed (for monitoring purposes) if you define $ZMAT, but the stepper operates in Cartesian coordinates only. No geometry constraints can be applied, apart from the point group in $DATA. A good paper to read about this kind of search is A.Farazdel, M.Dupuis J.Comput.Chem. 12, 276-282(1991) =========================================================== $CONICL group (relevant if RUNTYP=CONICAL) This group governs a search for the lowest energy on the 3N-7 dimensional "seam" of intersection of two (three) electronic potential energy surfaces of the same spin and space symmetry. Such Conical Intersections (CI) are important in photochemistry, where they serve as "funnels" for the transfer from an excited state to a lower state. See RUNTYP=MEX and the $MEX input for the simpler case where the two surfaces differ by either space or spin symmetry. Three search procedures are given, one of which requires the non-adiabatic coupling matrix element (NACME), and two others which do not require NACME information. The conical intersection search is available only for MCSCF (for which NACME are available) or for TD-DFT potential surfaces (where NACME are not available). The TD-DFT must be used in the Tamm/Dancoff approximation (see TAMMD in $TDDFT), but can be either conventional or spin-flip. The search utilizes some of the options of $STATPT, but note that the Schlegel stepper and HESS=CALC are not permitted. It may be reasonable to try the RFO stepper sometimes. The search can only be run in Cartesian coordinates. Restarts are possible only by updating the coordinates in $DATA. At present, the only solvation model that is supported is conventional TD-DFT with EFP1. OPTTYP = search procedure choice, see references below! = GPWNAC Gradient Projection with NACME, so this is only available for MCSCF. = BPUPD branching plane updating method (default) = PENALTY penalty-constrained optimization method = PENALTT penalty-constrained optimization method for three state conical intersection search See J. Phys.Chem.A 125.9 (2021) Note that for MCSCF surfaces, if state-averaging is used, the program executes the code needed to produce NACME vectors, to producing the state-averaged gradients. There is essentially no extra time required to produce also the NACME, hence the GPWNAC stepper might as well be used. IXROOT = array of two states whose CI point is sought. For example, this might be IXROOT(1)=2,3 The roots are counted exactly the same as IROOT in the $DET or $TDDFT input groups. For the latter case, set IXROOT to 0 if you want the ground state to be one of the two surfaces searched on. There is no default for IXROOT! ITROOT = array of three states whose CI point is sought. For example, this might be ITROOT(1)=2,3,4 The roots are counted exactly the same as IROOT in the $DET or $TDDFT input groups. There is no default for ITROOT. Only OPTTYP=PENALTY is supported. SYMOFF = flag to switch off point group symmetry, the default is .TRUE. DEBUG = flag to print debugging info, default is .FALSE. The following are meaningful only for OPTTYP=PENALTY: TOLSTP = energy difference tolerance default=1d-6 Hartree TOLGRD = gradient convergence tolerance default=5d-3 Hartree/Bohr ALPHA = parameter ensuring a singularity free penalty, default=0.02 Hartree SIGMA = Lagrange multiplier for the penalty term. In case the energy gap between the states is not acceptable at the CI point, increase the value. default = 3.5 (unitless) An understanding of the search procedures can be gained by reading the following papers: Gradient Projection with NACME: M.J.Bearpark, M.A.Robb, H.B.Schlegel Chem.Phys.Lett. 223, 269(1994) Branching Plane Updating method: S.Maeda, K.Ohno, K.Morokuma J.Chem.Theor Comput. 6, 1538(2010) Penalty constrained update method: B.G.Levine, C.Ko, J.Quenneville, T.J.Martinez Mol.Phys. 104, 1039(2006) B.G.Levine, J.D.Coe, T.J.Martinez J.Phys.Chem.B 112, 405(2008) A comparative study of the first two procedures is T.W.Keal, A.Koslowski, W.Thiel Theoret.Chem.Acc. 118, 837(2007) =========================================================== =========================================================== $MD group (relevant if RUNTYP=MD) This group controls the molecular dynamics trajectory for a collection of quantum mechanical atoms and/or Effective Fragment Potential particles. A typical MD simulation starts with an equilibration phase, running long enough to produce a randomized structure and velocity distribution. Typically equilibration is done with an NVT ensemble, allowing the system to equilibrate to a desired temperature. A production run restarts with the positions and the velocity and quaternion data from the equilibration run, might use either a NVE or NVT ensemble, and collects radial distribution functions and other properties. Only a few properties are computed from the MD trajectory, apart from correct radial distribution functions. In particular, the pressures, diffusion constants, and heats of vaporization that appear on the printout (presently only for pure EFP runs) are from a preliminary code, which has not yet been verified. If the system contains only EFP particles, it may be placed in a periodic box, according to the minimum image convention. The optional periodic boundary conditions, along with cut-offs, are given in the $EFRAG input. See also the $EWALD input group for long-range electrostatic treatment if PBC is used. The first keywords relate to the steps: MDINT = MD integrator selection. = FROG (leapfrog). This is less accurate, and lacks the special ensemble stepper option NVTNH. = VVERLET (velocity Verlet) - default. DT = MD time step size, in seconds, default=1.0d-15, which is a femtosecond. NVTNH selects a integrator step appropriate to the desired ensemble. This is only implemented for velocity Verlet. = 0 means use NVE Verlet stepping = 1 means use NVT Verlet stepping = 2 means use Nose/Hoover chain NVT Verlet stepping The default is 2 if either NVT option RSTEMP or RSRAND is chosen, but is 0 otherwise. NSTEPS = number of MD time steps to be found in this run, default=10000. TTOTAL = total time elapsed in the previous part of a MD trajectory which is being restarted (READ=.TRUE.). The default means this trajectory is a new one, or perhaps the start of a production phase of the MD. (default=0.0 seconds) * * * BATHT = bath temperature, in Kelvin (default=300.0) This value is used during NVT runs, or if the MD is initialized to a Maxwell-Boltzmann velocity distribution. In REMD, BATHT is an array with the size given by the number of replicas. * * * Two options exist to create NVT runs, to bring the system to a desired bath temperature. If neither is chosen, the ensemble is NVE: RSTEMP = flag to rescale the temperature. default=.FALSE. DTEMP = temperature range for the RSTEMP option. The velocities are rescaled to the bath temperature if T < (BATHT-DTEMP) or T > (BATHT+DTEMP). The default is DTEMP=100.0 degrees. RSRAND = flag to reset to Maxwell-Boltzmann distribution, using random numbers (same algorithm as MBT and MBR) to choose individual velocity magnitudes and directions. default=.FALSE. NRAND = number of steps for the RSRAND option. Reassign velocities (translational and rotational) every NRAND time steps. Default=1000. NVTOFF = step number at which to turn off either NVT thermostat, and switch to NVE. At this point, the NVTNH parameter will be reset to 0, and the PROD flag will be turned on, so that the production run will start (gathering and printing the RDF information to .log file). This keyword is also useful in NVE runs to postpone the accumulation of production information. The default means no switch to NVE (default=0). JEVERY = report simulation quantities (write info such as energies, temps, etc. to .log file) and collect RDF info each JEVERY time step. Default=10 KEVERY = write coordinates (to log and TRAJECT files), velocity/quaternion restart info (to the TRAJECT file and RDFs (to log file) at each KEVERY step. default=100 PROD = production run, at present this means only that information for radial distribution functions is collected, and printed. default=.FALSE. DELR = spacing for radial bins in RDF calculations, default=0.02 Angstroms. NPROP = step number at which to begin collecting data for the other properties, such as pressure and diffusion constants. This should be a value between 1 and NSTEPS, as it counts off the current run's steps. Default=0. PBCOUT = print PBC coordinates in the end of simulation (i.e. all molecules will be contained in one box) Default=.FALSE. * * * SSBP = flag to add spherical solvent boundary condition using the harmonic restraint potential(V) V=0.5*SFORCE*(R-DROFF)**2. (Default=.FALSE) SFORCE = the force constant for SSBP in kcal/mol-A**2 (default: 0.0). CCMS = flag to add a harmonic potential to constrain the center of mass of the QM subsystem. This will keep the QM subsystem in QM/MM-MD with spherical boundary conditions near the center of sphere. (Default=.FALSE) CFORCE = the force constant for CCMS in kcal/mol-A**2 (default: 0.0). DROFF = is an array of distances such that V=0 if R gives the number of RDFs which should be computed. Line 2. gives a string for the printout (a good choice involves both atoms, such as ClCl), the name of the $FRAGNAME containing the first atom of the pair, the name of the $FRAGNAME input group with the second atom of the pair, and how many such pairs exist. Line 3.