You should always start an MCSCF run with orbitals from some other run, by means of GUESS=MOREAD. Do not expect to be able to use HCORE, MINGUESS, or EXTGUESS! Example 6 is a poor example, converging only because methylene has so much symmetry, and the basis is so small. If you are just beginning your MCSCF problem, use orbitals from some appropriate converged SCF run. Thus, a realistic example of beginning a MCSCF calculation is examples 8 and 9. Once you get an MCSCF to converge, you can and should use these MCSCF MOs (which will be Schmidt orthogonalized) at other nearby geometries.

Starting from SCF orbitals can take a little bit of care. Most of the time (but not always) the orbitals you want to correlate will be the highest occupied orbitals in the SCF. Fairly often, however, the correlating orbitals you wish to use will not be the lowest unoccupied virtuals of the SCF. You will soon become familiar with NORDER=1 in $GUESS, as well as the nearly mandatory GUESS=MOREAD. You should also explore SYMLOC in $LOCAL, as the canonical SCF MOs are often spread out over regions of the molecule other than "where the action is".

Convergence of MCSCF is by no means guaranteed. Poor convergence can invariably be traced back to a poor initial selection of orbitals. My advice is, before you even start: "Look at the orbitals. Then look at the orbitals again". Later, if you have any trouble: "Look at the orbitals some more". Even if you don't have any trouble, look at the orbitals to see if they converged to what you expected, and have reasonable occupation numbers. MCSCF is by no means the sort of "black box" that RHF is these days, so LOOK VERY CAREFULLY AT YOUR RESULTS. Since orbitals are very hard to look at on the screen, this means print them out on paper, and inspect the individual MO expansion coefficients! Better yet, plot them and look at their nodal characteristics.

The iterative Process

Each iteration contains MICIT Newton-Raphson orbital improvements. After the first such NR step, a micro- iteration consists solely of an integral transformation over the new MOs, followed immediately by a NR orbital improvment, reusing the old second order density matrix. This avoids the CI diagonalization step, so may be of some use in MCSCF calculations with large CSF expansions. If the number of CSFs is small, using additional micro- iterations is a stupid idea

.

Miscellaneous Hints

Note that the partial transformation is coded only for small in-memory jobs, and even then may be slower than the full transformation. This is particularly true on vector machines, as the full transformation is vectorized by matrix multiply calls. Thus the default is the full integral transformation.

A very common MCSCF wavefunction, the simplest possible function describing a singlet diradical, has 2 electrons in 2 active MOs. While this wavefunction can be obtained in a MCSCF run (using NDOC=1 NVAL=1), it can be obtained much faster by use of the GVB code, with one GVB pair. This GVB-PP(1) wavefunction is also known in the literature as two configuration SCF (TCSCF). The two configuration form of this GVB is equivalent to the three configurations obtained in this MCSCF, as optimization in natural form (20 and 02) causes the coefficient of the 11 CSF to vanish.

One way to improve upon the SCF orbitals as starting MOs is to generate modified virtual orbitals (MVOs). MVOs are obtained by diagonalizing the Fock operator of a very positive ion, within the virtual orbital space only. As implemented in GAMESS, MVOs can be obtained at the end of any RHF, ROHF, or GVB run by setting MVOQ in $SCF nonzero, at the cost of a single SCF cycle. (Generating MVOs in a run feeding in converged orbitals is more expensive, as the 2 electron integrals must be recomputed). Generating MVOs will never change any of the occupied orbitals, but will generate more valence like LUMOs.

It is easy to generate large numbers of CSFs with the $DRT input (especially SOCI), which produces a bottleneck in the CI diagonalization step. If the number of CSFs is small (say under 10,000), then the most time consuming step in a CI run is the integral transformation, and in a MCSCF run is both the integral transformation and the Newton-Raphson orbital improvement. However, larger numbers of CSFs will require a great deal of time in the diagonalization step, as well as external disk storage for the pieces of the Hamiltonian matrix. The largest CI calculation performed at NDSU had about 85,000 CSFs, although we have heard of a CISD with over 200,000 CSFs that required 8 IBM 3380 disks to hold the Hamiltonian matrix. The largest MCSCF performed at NDSU had about 20,000 CSFs, although we have heard of one with about 35,000 CSFs. MCSCF typically is smaller, simply because this kind of run must perform a CI diagonalization once every iteration.

MCSCF Implementation

There are three different methods for obtaining the two electron contributions to the Lagrangian matrix and the orbital hessian in GAMESS (collectively called the NR matrices). These methods compute exactly the same NR matrices, so the overall number of iterations is the same for all methods.

The computational problem in forming the NR matrices is the combination of the second order density matrix, and the transformed two electron integrals. Both of these are 4 index matrices, but the DM2 is much smaller since its indices run only over occupied orbitals. The TEIs are much more numerous, as two of the indices run over all orbitals, including unoccupied ones.

METHOD=TEI is driven by the two electron integrals. After reading all of the DM2 into memory, this method begins to read the TEI file, and combines these integrals with the available DM2 elements. This method uses the least memory, reads the DM2 and TEI files only once each, and is usually very much slower than the other two methods.

METHOD=DM2 is the reverse strategy, with all the integrals in core, and is driven by reading DM2 elements from disk. Because the number of TEIs is large, the implementation actually divides the TEIs into four sets, (ij/kl), (aj/kl), (ab/kl), and (aj/bl), where i,j,k,l run over the internal orbitals, and a and b are external orbitals. The program must have enough memory to hold all the TEI's belonging to a particular set in memory. Thus, construction of the NR matrices requires 4 reads of the TEI file, and 4 reads of the DM2. METHOD=DM2 requires more memory than METHOD=TEI, and has 4 times the I/O, but is much faster.

METHOD=FORMULA tries to combine the merits of both methods. This method constructs a Newton formula tape on the first MCSCF iteration at the first geometry computed. This formula tape contains integer pointers telling how TEIs should be combined with DM2 elements to form the NR matrices. If all of the TEIs will not fit into memory, the formula tape must be sorted after it is generated. The formula tape can be saved for use in later runs on the same molecule, see FMLFIL in $MCSCF. This can be a very large disk file, which may diminish your enthusiasm for saving it.

Once the formula tape is constructed, the "passes" (say M in number) begin. A pass consists of reading as large a block of the TEI file as will fit into the available memory, then reading the related portion of the formula tape and the entire DM2 file for information about how to efficiently combine the TEI's and the DM2. When all of this TEI block is processed, the next block of TEI's is read to begin the next "pass". An iteration thus performs one complete read through the TEI and formula files, and M reads of the DM2. METHOD=FORMULA uses as much memory as the program has available (above a certain minimum), by adjusting the chunk size used for the TEI file), has an I/O requirement somewhat larger than METHOD=TEI, and is perhaps a bit faster than METHOD=DM2. It does require additional disk space for the Newton formula tape.

Unlike METHOD=DM2 and METHOD=TEI, METHOD=FORMULA requires a canonical order file of transformed integrals. The other two methods use Steve Elbert's reverse canonical ordering, which allows the integral transformation to run 14% faster. Thus, METHOD=FORMULA may lose in the transformation whatever it gains over METHOD=DM2 in the Newton-Raphson step.

As all three methods have advantages and disadvantages, all three are available. The amount of memory required for each method will be included in your output (including EXETYP=CHECK jobs) to assist you in selecting a method.

References

1. "The Multiconfiguration (MC) SCF Method"
   B.O.Roos, in "Methods in Computational Molecular Physics", edited by G.H.F.Diercksen and S.Wilson D.Reidel Publishing, Dordrecht, Netherlands, 1983, pp 161-187.

2. "Optimization and Characterization of a MCSCF State"
   Olsen, Yeager, Jorgensen in "Advances in Chemical Physics", vol.54, edited by I.Prigogine and S.A.Rice, Wiley Interscience, New York, 1983, pp 1-176.

3. "Matrix Formulated Direct MCSCF and Multiconfiguration Reference CI Methods"
   H.-J.Werner, in "Advances in Chemical Physics", vol.69, edited by K.P.Lawley, Wiley Interscience, New York, 1987, pp 1-62.

4. "The MCSCF Method" R.Shepard, ibid., pp 63-200.

Some papers particularly germane to the implementation in GAMESS are

5. "General second order MCSCF theory: A Density Matrix Directed Algorithm"
   B.H.Lengsfield, III, J.Chem.Phys. 73,382-390(1980).

6. "The use of the Augmented Matrix in MCSCF Theory"
   D.R.Yarkony, Chem.Phys.Lett. 77,634-635(1981).

The next set of references discuss the choice of orbitals for the active space. This is a matter of some importance, and needs to be understood well by the GAMESS user.

7. "The CASSCF Method and its Application in Electronic Structure Calculations"
   B.O.Roos, in "Advances in Chemical Physics", vol.69, edited by K.P.Lawley, Wiley Interscience, New York, 1987, pp 339-445.

8. "Are Atoms Intrinsic to Molecular Electronic Wavefunctions?"
   K.Ruedenberg, M.W.Schmidt, M.M.Dombek, S.T.Elbert Chem.Phys. 71, 41-49, 51-64, 65-78 (1982).

The final references are simply some examples of FORS MCSCF applications, the latter done with GAMESS.

9. D.F.Feller, M.W.Schmidt, K.Ruedenberg,
   J.Am.Chem.Soc. 104, 960-967(1982).

10. M.W.Schmidt, P.N.Truong, M.S.Gordon,
    J.Am.Chem.Soc. 109, 5217-5227(1987).

GAMESS locates minima by geometry optimization, RUNTYP=OPTIMIZE, and transition states by saddle point searches, RUNTYP=SADPOINT. In many ways these are similar, and in fact nearly identical FORTRAN code is used for both. The term "geometry search" is used here to describe features which are common to both procedures. The input to control both RUNTYPs is found in the $STATPT group.

As will be noted in the symmetry section below, an OPTIMIZE run does not always find a minimum, and a SADPOINT run may not find a transtion state, even though the gradient is brought to zero. You can prove you have located a minimum or saddle point only by examining the local curvatures of the potential energy surface. This can be done by following the geometry search with a RUNTYP=HESSIAN job, which should be a matter of routine.

Quasi-Newton Searches

GAMESS contains two different implementations of the quasi-Newton procedure for finding stationary points, namely METHOD=SCHLEGEL and METHOD=BAKER. For minima searches, the two procedures vary mainly in the algorithm that is used to approximately update the hessian as the search proceeds. For saddle point runs, Baker's method treats the unique direction of negative curvature in a more sophisticated way than the Schlegel algorithm. For both kinds of searches, Baker's method seems to be the more robust, boasting the ability (in some cases) to locate transition states even when started in a region of all positive local curvatures.

You should read the papers by Baker and Schlegel to understand how these methods are designed to work. See section 2 of the paper by Dunning's group for the Argonne approach to transition searches. A good summary of various search procedures is given by Bell and Crighton, and in the review by Schlegel. These papers are:

1. J.Baker J.Comput.Chem. 7,385-395(1986)

2. H.B.Schlegel J.Comput.Chem. 3,214(1982)

3. T.H.Dunning, L.B.Harding, R.A.Blair, R.A.Eades, R.L.Shepard J.Phys.Chem. 90, 344-356, (1986).

4. S.Bell, J.S.Crighton J.Chem.Phys. 80, 2464-2475, (1984).

5. H.B.Schlegel Advances in Chemical Physics (Ab Initio Methods in Quantum Chemistry, Part I), volume 67, K.P.Lawley, Ed. Wiley, New York, 1987, pp 249-286.

There is a procedure contained within GAMESS for guessing a diagonal, positive definite hessian matrix, HESS=GUESS. This guess is more sophisticated when internal coordinates are defined, as Badger's rules can be used to estimate stretching force constants. This option often works well for minima, but cannot possibly help you find transition states (because it is positive definite). Therefore, GUESS cannot be selected for SADPOINT runs.

Two options for providing a more accurate hessian are HESS=READ and CALC. In the latter, the true hessian is obtained by direct calculation at the initial geometry, and then the geometry search begins, all in one run. The READ option allows you to feed in the hessian in a $HESS group, as obtained by a RUNTYP=HESSIAN job. The second procedure may actually be preferable, as you get a chance to see the frequencies. Then, if the local curvatures look good, you can commit to the geometry search. Be sure to include the $GRAD group (if the exact gradient is available) in the HESS=READ job so that GAMESS can take its first step immediately.

Note also that you can compute the hessian at a lower basis set and/or wavefunction level, and read this into a higher level geometry search. This trick works because the hessian is 3Nx3N for N atoms, no matter what atomic basis is used. The gradient from the lower level is of course worthless, as the geometry search must work with the exact gradient.

You get what you pay for here. HESS=GUESS is free, but may lead to significantly more steps in the geometry search. The other two options are more expensive at the beginning, but may pay back by rapid convergence to the stationary point.

The procedure for updating the hessian as the geometry search progresses frequently improves the hessian for a few steps (especially for HESS=GUESS), and then breaks down. The symptoms of this are a nice lowering of the energy or RMS gradient for maybe 10 steps, followed by crazy steps. You can help by putting the best coordinates into $DATA, and resubmitting the job, which will then make a fresh determination of the hessian.

Internal Coordinates

Optimizations in Cartesian space are relatively uncoupled, whereas internal coordinates are frequently strongly coupled. Because of this, Jerry Boatz calls them "infernal coordinates"! An very common example to illustrate this might be a bond length in a ring, and an angle on the opposite side of the ring. Clearly, changing one changes the other simultaneously. A more mathematical definition of "coupled" might be to say that there is a large off diagonal element in the hessian. In this case convergence might well be unsatisfactory, especially with the diagonal GUESS hessian, but possibly even with an accurately calculated one.

Particularly poorly chosen coordinates may not correspond to a quadratic surface, thereby ending all hope that a quasi-Newton method will converge.

The best procedure is to pick a set of internal coordinates which are as independent as possible. For a straight chain molecule of N atoms, one might do this according to the usual "model builder" Z-matrix input, using N-1 bond lengths, N-2 bond angles, and N-3 torsions. However, you should note that GAMESS does not require this particular mixture of the common types. The only requirement is that you provide a total of 3N-6 coordinates, chosen from these three basic types or several more exotic possibilities. (Of course, we mean 3N-5 throughout for linear molecules).

For example, the most effective choice of coordinates for the atoms in a four membered ring is to define all four sides, any one of the internal angles, and a dihedral defining the ring pucker. For a six membered ring, the best coordinates seem to be 6 sides, 3 angles, and 3 torsions. The angles should be every other internal angle, so that the molecule can "breathe" freely. The torsions should be arranged so that the central bond of each is placed on alternating sides of the ring, as if they were pi bonds in Kekule benzene. For a five membered ring, we suggest all 5 sides, 2 internal angles, alternating every other one, and 2 dihedrals to fill in. The internal angles of necessity skip two atoms where the ring closes. Larger rings should generalize on the idea of using all sides and only alternating angles.

Rings and more especially fused rings can be tricky. In fact, these systems may work better in Cartesian coordinates, and also may benefit from decreasing DXMAX to avoid large steps. Rings nearly always benefit from a computed hessian, rather than the diagonal GUESS one, even when looking for a minimum.

GAMESS also has several other types of internal coordinates, such as linear bends for 3 collinear atoms, out of plane bends, and so on. If you experience difficulty with one set of internals, consider using some of these more exotic types. A little thought about a set of coordinates which are natural for your problem can be the difference between good convergence of the geometry search, and poor or even no convergence.

In many cases the "model builder" set of coordinates (N-1 bonds, N-2 angles, and N-3 torsions) is actually an appropriate one. If you have entered your molecule in $DATA using one of the Z-matrix options, GAMESS can automatically build a $ZMAT using the "model builder" coordinates given in $DATA. This is done by setting NZVAR to 3N-6, while not entering a $ZMAT. Note that this option will not work if you have used any dummy atoms to enter the molecule. If you want to override the "model builder" choice with some other set, simply leave NZVAR at 3N-6, and provide the appropriate $ZMAT.

A poor choice of internal coordinates will sometimes cause GAMESS to stop with a unfriendly message about "Cartesian displacement exceeds 3*DXMAX". This is often caused by three angles around an atom summing to nearly 360 degrees (an almost planar center), but can be due to other causes as well. The best fix for this is to look for a set of internals which is better behaved. Sometimes providing a better hessian than GUESS will help, but a willingness to think about alternate coordinates is even better. A little brain time can save a lot of CPU time.

If the molecule has symmetry, be sure to define coordinates which are symmetrically related.

Sometimes it is handy to constrain the geometry search by freezing one or more coordinates, via the IFREEZ array. For example, such constrained optimizations might be useful while trying to determine what area of the potential surface contains a saddle point. GAMESS will not actually hold such frozen coordinates truly constant, instead they sometimes drift a small amount. This does not trouble us in the least, as a structure obtained under such a constraint is neither a minimum or a saddle point anyway. If you are trying to freeze coordinates with an automatically generated $ZMAT, you need to know that the order of the coordinates is

      y
      y  x r1
      y  x r2  x a3
      y  x r4  x a5  x w6
      y  x r7  x a8  x w9

The Role of Symmetry

This apparent mystery is due to the fact that the gradient vector transforms under the totally symmetric representation of the molecular point group. As a direct consequence, a geometry search is point group conserving. (For a proof of these statements, see J.W.McIver and A.Komornicki, Chem.Phys.Lett., 10,303-306(1971)). In simpler terms, the molecule will remain in whatever point group you select in $DATA, even if the true minimum is in some lower point group. Since a geometry search only explores totally symmetric degrees of freedom, the only way to learn about the curvatures for all degrees of freedom is RUNTYP=HESSIAN.

As an example, consider disilene, the silicon analog of ethene. It is natural to assume that this molecule is planar like ethene, and an OPTIMIZE run in D2h symmetry will readily locate a stationary point. However, as a calculation of the hessian will readily show, this structure is a transition state (one imaginary frequency), and the molecule is really trans-bent (C2h). A careful worker will always characterize a stationary point as either a minimum, a transition state, or some higher order stationary point (which is not of great interest!) by performing a RUNTYP=HESSIAN.

The point group conserving properties of a geometry search can be annoying, as in the preceeding example, or advantageous. For example, assume you wish to locate the transition state for rotation about the double bond in ethene. A little thought will soon reveal that ethene is D2h, the 90 degrees twisted structure is D2d, and structures in between are D2. Since the saddle point is actually higher symmetry than the rest of the rotational surface, you can locate it by RUNTYP=OPTIMIZE within D2d symmetry. You can readily find this stationary point with the diagonal guess hessian! In fact, if you attempt to do a RUNTYP=SADPOINT within D2d symmetry, there will be no totally symmetric modes with negative curvatures, and it is unlikely that the geometry search will be very well behaved.

Although a geometry search cannot lower the symmetry, the gain of symmetry is quite possible. For example, if you initiate a water molecule optimization with a trial structure which has unequal bond lengths, the geometry search will come to a structure that is indeed C2v (to within OPTTOL, anyway). However, GAMESS leaves it up to you to realize that a gain of symmetry has occurred.

In general, Mother Nature usually chooses more symmetrical structures over less symmetrical structures. Therefore you are probably better served to assume the higher symmetry, perform the geometry search, and then check the stationary point's curvatures. The alternative is to start with artificially lower symmetry and see if your system regains higher symmetry. The problem with this approach is that you don't necessarily know which subgroup is appropriate, and you lose the great speedups GAMESS can obtain from proper use of symmetry. It is good to note here that "lower symmetry" does not mean simply changing the name of the point group and entering more atoms in $DATA, instead the nuclear coordinates themselves must actually be of lower symmetry.

Practical Matters

If a geometry search should happen to run out of time, or to exceed NSTEP before meeting the OPTTOL criterion, it can easily be restarted. For RUNTYP=OPTIMIZE, the restart should use the coordinates with the lowest total energy (do a string search on "FINAL"). For RUNTYP=SADPOINT, the restart should use the coordinates which have the smallest gradient (do a string search on "RMS", which means root mean square). Note that these are not necessarily the last geometry considered by the program!

The "restart" should actually be a normal run, that is you should not try to use the restart options in $CONTRL (which anyway may not work). A geometry search can be restarted by extracting the desired coordinates for $DATA from the printout, and by extracting the corresponding $GRAD group from the PUNCH file. If the $GRAD group is supplied, the program is able to save the time it would ordinarily take to compute the wavefunction and gradient at the initial point, which can be a substantial savings. There is no input to trigger reading of a $GRAD group: if found, it is read and used. Be careful that your $GRAD group actually corresponds to the coordinates in $DATA, as GAMESS has no check for this.

Sometimes when you are fairly close to the minimum, an OPTIMIZE run will take a first step which raises the energy, with subsequent steps improving the energy and perhaps finding the minimum. The erratic first step is caused by the GUESS hessian, and the subsequent convergence by the updating procedures correcting this guess. It can be helpful to limit the size of this wrong first step, by reducing the radius of the step, DXMAX. Conversely, if you are far from the minimum, sometimes you can decrease the total number of steps by increasing DXMAX.

Saddle Points

In contrast, finding saddle points is a black art. The diagonal guess hessian will never work, so you must provide a computed one. The hessian should be computed at your best guess as to what the transition state should be. It is safest to do this in two steps as outlined above, rather than HESS=CALC. This lets you verify you have guessed a structure with one and only one negative curvature. Guessing a good trial structure is the hardest part of a RUNTYP=SADPOINT!

The RUNTYP=HESSIAN's normal coordinate analysis is rigorously valid only at stationary points on the surface. This means the frequencies from the hessian at your trial geometry are untrustworthy, in particular the six "zero" frequencies corresponding to translational and rotational degrees of freedom will usually be 300-500 cm**-1, and possibly imaginary. The Sayvetz conditions on the printout will help you distinguish the T&R "contaminants" from the real vibrational modes. If you have defined a $ZMAT, the PURIFY option in $FORCE will help zap out these T&R contaminants), but this is not really necessary.

If the hessian at your assumed geometry does not have one and only one imaginary frequency (taking into account that the "zero" frequencies can sometimes be 300i!), then it will probably be difficult to find the saddle point. Instead you need to compute a hessian at a better guess for the initial geometry, or else read about mode following below.

If you need to restart your run, do so with the coordinates which have the smallest RMS gradient. Note that the energy does not necessarily have to decrease in a SADPOINT run, in contrast to an OPTIMIZE run. It is often necessary to do several restarts, involving recomputation of the hessian, before actually locating the saddle point.

Assuming you do find the T.S., it is always a good idea to recompute the hessian at this structure. As described in the discussion of symmetry, only totally symmetric vibrational modes are probed in a geometry search. Thus it is fairly common to find that at your "T.S." there is a second imaginary frequency, which corresponds to a non-totally symmetric vibration. This means you haven't found the correct T.S., and are back to the drawing board. The proper procedure is to lower the point group symmetry by distorting along the symmetry breaking "extra" imaginary mode, by a reasonable amount. Don't be overly timid in the amount of distortion, or the next run will come back to the invalid structure.

The real trick here is to find a good guess for the transition structure. The closer you are, the better. It is often difficult to guess these structures. One way around this is to compute a linear least motion (LLM) path. This connects the reactant structure to the product structure by linearly varying each coordinate. If you generate about ten structures intermediate to reactants and products, and compute the energy at each point, you will in general find that the energy first goes up, and then down. The maximum energy structure is a "good" guess for the true T.S. structure. Actually, the success of this method depends on how curved the reaction path is.

Mode Following

Mode following is really not a substitute for the ability to intuit regions of the PES with a single local negative curvature. Baker's case D works well, but case E is much trickier. When you start near a minimum, it matters a great deal which side of the minima you start from, as the direction of the search depends on the sign of the gradient. Be sure to use your brain when you use IFOLOW!

MOPAC Calculations in GAMESS

1. the one- and two-electron integrals,

2. the two-electron part of the Fock matrix,

3. the cartesian energy derivatives, and

4. the ZDO atomic charges and molecular dipole.

Everything else is actually GAMESS: coordinate input (including point group symmetry), the SCF convergence procedures, the matrix diagonalizer, the geometry searcher, the numerical hessian driver, and so on. Most of the output will look like an ab initio output.

It is extremely simple to perform one of these calculations. All you need to do is specify GBASIS=MNDO, AM1, or PM3 in the $BASIS group. Note that this not only picks a particular Slater orbital basis, it also selects a particular "hamiltonian", namely a particular parameter set.

MNDO, AM1, and PM3 will not work with every option in GAMESS. Currently the semiempirical wavefunctions support SCFTYP=RHF, UHF, and ROHF in any combination with RUNTYP=ENERGY, GRADIENT, OPTIMIZE, SADPOINT, HESSIAN, and IRC. Note that all hessian runs are numerical finite differencing. The MOPAC CI and half electron methods are not supported.

Because the majority of the implementation is GAMESS rather than MOPAC you will notice a few improvments. Dynamic memory allocation is used, so that GAMESS uses far less memory for a given size of molecule. The starting orbitals for SCF calculations are generated by a Huckel initial guess routine. Spin restricted (high spin) ROHF can be performed. Converged SCF orbitals will be labeled by their symmetry type. Numerical hessians will make use of point group symmetry, so that only the symmetry unique atoms need to be displaced. Infrared intensities will be calculated at the end of hessian runs. We have not at present used the block diagonalizer during intermediate SCF iterations, so that the run time for a single geometry

point in GAMESS is usually longer than in MOPAC. However, the geometry optimizer in GAMESS can frequently optimize the structure in fewer steps than the procedure in MOPAC. Orbitals and hessians are punched out for convenient reuse in subsequent calculations. Your molecular orbitals can be drawn with the PLTORB graphics program.

To reduce CPU time, only the EXTRAP convergence accelerator is used by the SCF procdures. For difficult cases, the DIIS, RSTRCT, and/or SHIFT options will work, but may add significantly to the run time. With the Huckel guess, most calculations will converge acceptably without these special options.

MOPAC parameters exist for the following elements. The quote means that these elements are treated as "sparkles" rather than as atoms with genuine basis functions. For MNDO:

 H
Li  *          B  C  N  O  F
Na' *         Al Si  P  S Cl
 K' * ...  Zn  * Ge  *  * Br
Rb' * ...   *  * Sn  *  *  I
*   * ...  Hg  * Pb  *

         For AM1:                         For PM3:
 H                               H
 *  *          B  C  N  O  F     *  Be         *  C  N  O  F
Na' *         Al Si  P  S Cl    Na' Mg        Al Si  P  S Cl
 K' * ...  Zn  * Ge  *  * Br     K' * ...  Zn Ga Ge As Se Br
Rb' * ...   *  * Sn  *  *  I    Rb' * ...  Cd In Sn Sb Te  I
*   * ...  Hg  *  *  *          *   * ...  Hg Tl Pb Bi

It is important to exercise caution as semiempirical methods can be dead wrong! The reasons for this are bad parameters (in certain chemical situations), and the underlying minimal basis set. A good question to ask before using MNDO, AM1 or PM3 is "how well is my system modeled with an ab initio minimal basis sets, such as STO-3G?" If the answer is "not very well" there is a good chance that a semiempirical description is equally poor.

Molecular Properties and Conversion Factors

All units are derived from the atomic units for distance and the monopole electric charge, as given below.

distance               - 1 au = 5.291771E-09 cm

monopole               - 1 au = 4.803242E-10 esu
                        1 esu = sqrt(g-cm**3)/sec

dipole                 - 1 au = 2.541766E-18 esu-cm
                      1 Debye = 1.0E-18 esu-cm

quadrupole             - 1 au = 1.345044E-26 esu-cm**2
                 1 Buckingham = 1.0E-26 esu-cm**2

octupole               - 1 au = 7.117668E-35 esu-cm**3

electric potential     - 1 au = 9.076814E-02 esu/cm

electric field         - 1 au = 1.715270E+07 esu/cm**2
                  1 esu/cm**2 = 1 dyne/esu

electric field gradient- 1 au = 3.241390E+15 esu/cm**3

The atomic unit for spin density is excess alpha spins per unit volume, h/4*pi*bohr**3. Only the expectation value is computed, with no constants premultiplying it.

IR intensities are printed in Debye**2/amu-Angstrom**2. These can be converted into intensities as defined by Wilson, Decius, and Cross's equation 7.9.25, in km/mole, by multiplying by 42.255. If you prefer 1/atm-cm**2, use a conversion factor of 171.65 instead. A good reference for deciphering these units is A.Komornicki, R.L.Jaffe J.Chem.Phys. 1979, 71, 2150-2155. A reference showing how IR intensities change with basis improvements at the HF level is Y.Yamaguchi, M.Frisch, J.Gaw, H.F.Schaefer, J.S.Binkley, J.Chem.Phys. 1986, 84, 2262-2278.

Localization Tips

The method due to Edmiston and Ruedenberg works by maximizing the sum of the orbitals' two electron self repulsion integrals. Most people who think about the different localization criteria end up concluding that this one seems superior. The method requires the two electron integrals, transformed into the molecular orbital basis. Because only the integrals involving the orbitals to be localized are needed, the integral transformation is actually not very time consuming.

The Boys method maximizes the sum of the distances between the orbital centroids, that is the difference in the orbital dipole moments.

The population method due to Pipek and Mezey maximizes a certain sum of gross atomic Mulliken populations. This procedure will not mix sigma and pi bonds to yield localized banana bonds. Hence it is rather easy to find cases where this method give rather different results than the Ruedenberg or Boys approach.

GAMESS will localize orbitals for any kind of RHF, UHF, ROHF, or MCSCF wavefunctions. The localizations will automatically restrict any rotation that would cause the energy of the wavefunction to be changed (the total wavefunction is left invariant). As discussed below, localizations for GVB or CI functions are not permitted.

The default is to freeze core orbitals. The localized valence orbitals are scarcely changed if the core orbitals are included, and it is usually convenient to leave them out. Therefore, the default localizations are: RHF functions localize all doubly occupied valence orbitals. UHF functions localize all valence alpha, and then all valence beta orbitals. ROHF functions localize all doubly occupied orbitals, and all singly occupied orbitals, but do not mix these two orbital spaces. MCSCF functions localize all valence MCC type orbitals, and localize all active orbitals, but do not mix these two orbital spaces. To recover invariant an MCSCF function, you must be using a FORS=.TRUE. wavefunction, and you must set GROUP=C1 in $DRT as the localized orbitals possess no symmetry.

In general, GVB functions are invariant only to localizations of the NCO doubly occupied orbitals. Any pairs must be written in natural form, so pair orbitals cannot be localized. The open shells may be degenerate, so in general these should not be mixed. If for some reason you feel you must localize the doubly occupied space, do a RUNTYP=PROP job. Feed in the GVB orbitals, but tell the program it is SCFTYP=RHF, and enter a negative ICHARG so that GAMESS thinks all orbitals occupied in the GVB are occupied in this fictitous RHF. Use NINA or NOUTA to localize the desired doubly occupied orbitals. Orbital localization is not permitted for CI, because we cannot imagine why you would want to do that anyway.

Boys localization of the core orbitals in molecules having elements from the third or higher row almost never succeeds. Boys localization including the core for second row atoms will often work, since there is only one inner shell on these. The Ruedenberg method should work for any element, although including core orbitals in the integral transformation is more expensive.

The easiest way to do localization is in the run which determines the wavefunction, by selecting LOCAL=xxx in the $CONTRL group. However, localization may be conveniently done at any time after determination of the wavefunction, by executing a RUNTYP=PROP job. This will require only $CONTRL, $BASIS/$DATA, $GUESS (pick MOREAD), the converged $VEC, possibly $SCF or $DRT to define your wavefunction, and optionally some $LOCAL input.

There is an option to restrict all rotations that would mix orbitals of different symmetries. SYMLOC=.TRUE. yields only partially localized orbitals, but these still possess symmetry. They are therefore very useful as starting orbitals for MCSCF or GVB-PP calculations. Because they still have symmetry, these partially localized orbitals run as efficiently as the canonical orbitals. Because it is much easier for a user to pick out the bonds which are to be correlated, a significant number of iterations can be saved, and convergence to false solutions is less likely.

Intrinsic Reaction Coordinate (IRC) Methods

The IRC is becoming an important part of polyatomic dynamics research, as it is hoped that only knowledge of the PES in the vicinity of the IRC is needed for prediction of reaction rates, at least at threshhold energies. The IRC has a number of uses for electronic structure purposes as well. These include the proof that a certain transition structure does indeed connect a particular set of reactants and products, as the structure and imaginary frequency normal mode at the saddle point do not always unambiguously identify the reactants and products. The study of the electronic and geometric structure along the IRC is also of interest. For example, one can obtain localized orbitals along the path to determine when bonds break or form.

The accuracy to which the IRC is determined is dictated by the use one intends for it. Dynamical calculations require a very accurate determination of the path, as derivative information (second derivatives of the PES at various IRC points, and path curvature) is required later. Thus, a sophisticated integration method (such as AMPC4 or RK4), and small step sizes (STRIDE=0.05, 0.01, or even smaller) will be needed. In addition to this, care should be taken to locate the transition state carefully (perhaps decreasing OPTTOL by a factor of 10), and in the initiation of the IRC run. The latter might require a hessian matrix obtained by double differencing, certainly the hessian should be PURIFY'd. Note also that EVIB must be chosen carefully, as decribed below.

On the other hand, identification of reactants and products allows for much larger step sizes, and cruder integration methods. In this type of IRC one might want to be careful in leaving the saddle point (perhaps STRIDE should be reduced to 0.10 or 0.05 for the first few steps away from the transition state), but once a few points have been taken, larger step sizes can be employed. In general, the defaults in the $IRC group are set up for this latter, cruder quality IRC. The STRIDE value for the GS2 method can usually be safely larger than for other methods, no matter what your interest in accuracy is.

The simplest method of determining an IRC is linear gradient following, PACE=LINEAR. This method is also known as Euler's method. If you are employing PACE=LINEAR, you can select "stabilization" of the reaction path by the Ishida, Morokuma, Komornicki method. This type of corrector has no apparent mathematical basis, but works rather well since the bisector usually intersects the reaction path at right angles (for small step sizes). The ELBOW variable allows for a method intermediate to LINEAR and stabilized LINEAR, in that the stabilization will be skipped if the gradients at the original IRC point, and at the result of a linear prediction step form an angle greater than ELBOW. Set ELBOW=180 to always perform the stabilization.

A closely related method is PACE=QUAD, which fits a quadratic polynomial to the gradient at the current and immediately previous IRC point to predict the next point. This pace has the same computational requirement as LINEAR, and is slightly more accurate due to the reuse of the old gradient. However, stabilization is not possible for this pace, thus a stabilized LINEAR path is usually more accurate than QUAD.

Two rather more sophisticated methods for integrating the IRC equations are the fourth order Adams-Moulton predictor-corrector (PACE=AMPC4) and fourth order Runge- Kutta (PACE=RK4). AMPC4 takes a step towards the next IRC point (prediction), and based on the gradient found at this point (in the near vincinity of the next IRC point) obtains a modified step to the desired IRC point (correction). AMPC4 uses variable step sizes, based on the input STRIDE. RK4 takes several steps part way toward the next IRC point, and uses the gradient at these points to predict the next IRC point. RK4 is the most accurate integration method implemented in GAMESS, and is also the most time consuming.

The Gonzalez-Schlegel 2nd order method finds the next IRC point by a constrained optimization on the surface of a hypersphere, centered at 1/2 STRIDE along the gradient vector leading from the previous IRC point. By construction, the reaction path between two successive IRC points is thus a circle tangent to the two gradient vectors. The algorithm is much more robust for large steps than the other methods, so you may wish to try STRIDE=0.4. The number of energy and gradients need to find the next point depends on the difficulty of the constrained optimization, but is normally not very many points. The defaults for GCUT and RCUT are meant for use with STRIDE=0.4, so if you decide to use a smaller stride you should reduce these parameters accordingly. Be sure to provide the updated hessian from the previous run when restarting PACE=GS2.

The number of wavefunction evaluations, and energy gradients needed to jump from one point on the IRC to the next point are summarized in the following table:

     PACE      # energies   # gradients
     ----      ----------   -----------
    LINEAR        1             1
stabilized
    LINEAR        3             2
    QUAD          1             1  (+ reuse of historical
                                            gradient)
    AMPC4         2             2  (see note)
    RK4           4             4
    GS2          2-4           2-4 (equal numbers)

* * *

The only method which proves that a properly converged IRC has been obtained is to regenerate the IRC with a smaller step size, and check that the IRC is unchanged. Again, note that the care with which an IRC must be obtained is highly dependent on what use it is intended for.

A small program which converts the IRC results punched to file IRCDATA into a format close to that required by the POLYRATE VTST dynamics program written in Don Truhlar's group is available. For a copy of this IRCED program, contact Mike Schmidt. The POLYRATE program must be obtained from the Truhlar group.

    Some key IRC references are:
K.Ishida, K.Morokuma, A.Komornicki
      J.Chem.Phys.  66, 2153-2156 (1977)
K.Muller
      Angew.Chem., Int.Ed.Engl.19, 1-13 (1980)
M.W.Schmidt, M.S.Gordon, M.Dupuis
      J.Am.Chem.Soc.  107, 2585-2589 (1985)
B.C.Garrett, M.J.Redmon, R.Steckler, D.G.Truhlar,
K.K.Baldridge, D.Bartol, M.W.Schmidt, M.S.Gordon
      J.Phys.Chem.  92, 1476-1488(1988)
K.K.Baldridge, M.S.Gordon, R.Steckler, D.G.Truhlar
      J.Phys.Chem.  93, 5107-5119(1989)
C.Gonzales, H.B.Schlegel
      J.Phys.Chem.  94, 5523-5527(1990)
      J.Chem.Phys.  95, 5853-5860(1991)

Transition Moments and Spin-Orbit Coupling

The orbitals for the CI can be one common set of orbitals used by both CI states. If one set of orbitals is used, the transition moment or spin-orbit coupling can be found for any type of CI wavefunction GAMESS can compute.

Alternatively, two sets of orbitals (obtained from separate MCSCF orbital optimizations) can be used. Two separate CIs will be carried out (in 1 run). The two MO sets must share a common set of frozen core orbitals, and the CI -must- be of the complete active space type. These restrictions are needed to leave the CI wavefunctions invariant under the necessary transformation to corresponding orbitals. The non-orthogonal procedure implemented in GAMESS is a GUGA driven equivalent to the method of Lengsfield, et al. Note that the FOCI and SOCI methods described by these workers are not available in GAMESS.

If you would like to use separate orbitals for the states, use the FCORE option in $MCSCF in a SCFTYP=MCSCF optimization of the orbitals. Typically you would optimize the ground state completely, and use these MCSCF orbitals in an optimization of the excited state, under the constraint of FCORE=.TRUE. The core orbitals of such a MCSCF optimization should be declared MCC, rather than FZC, even though they are frozen.

In the case of transition moments either one or two CI calculations are performed, necessarily on states of the same multiplicity. Thus, only a $DRT1 is read. A spin-orbit coupling run almost always does two CI calculations, as the states are usually of different multiplicities. So, spin-orbit runs might read only $DRT1, but usually read both $DRT1 and $DRT2. The first CI calculation, defined by $DRT1, must be for the state of lower multiplicity, with $DRT2 defining the state of higher multiplicity.

You will probably have to lower the symmetry in $DRT1 and $DRT2 to C1. You may use full spatial symmetry in the DRT groups only if the two states happen to have the same total spatial symmetry.

The transition moment and spin orbit coupling driver is a rather restricted path through GAMESS, in that

1. GUESS=SKIP is forced, regardless of what you put in a $GUESS group. The rest of this group is ignored. The NOCC orbitals are read (as GUESS=MOREAD) in a $VEC1 group, and possibly a $VEC2 group. It is not possible to reorder the orbitals.

2. $CIINP is not read. The CI is hardwired to consist of DRT generation, integral transformation and sorting, Hamiltonian generation, and diagonalization. This means $DRT1 (and maybe $DRT2), $TRNSFM, $CISORT, $GUGEM, and $GUGDIA input is read, and acted upon.

3. The density matrices are not generated, and so no properties (other than the transition moment or the spin-orbit coupling) are computed.

4. There is no restart capability provided.

5. DRT input is given in $DRT1 and maybe $DRT2.

6. IROOTS will determine the number of eigenvectors found in each CI calculation, except NSTATE in $GUGDIA will override IROOTS if it is larger.

References:

F.Weinhold, J.Chem.Phys. 54,1874-1881(1970)

B.H.Lengsfield, III, J.A.Jafri, D.H.Phillips, C.W.Bauschlicher, Jr. J.Chem.Phys. 74,6849-6856(1981)

"Intramediate Quantum Mechanics, 3rd Ed." Hans A. Bethe, Roman Jackiw Benjamin/Cummings Publishing, Menlo Park, CA (1986), chapters 10 and 11.

For an application, S.Koseki, M.S.Gordon J.Mol.Spectrosc. 123, 392-404(1987)

* * *