How to do RHF, ROHF, UHF, and GVB calculations general considerations These four SCF wavefunctions are all based on Fock operator techniques, even though some GVB runs use more than one determinant. Thus all of these have an intrinsic N**4 time dependence, because they are all driven by integrals in the AO basis. This similarity makes it convenient to discuss them all together. In this section we will use the term HF to refer generically to any of these four wavefunctions, including the multi-determinate GVB-PP functions. $SCF is the main input group for all these HF wavefunctions. As will be discussed below, in GAMESS the term ROHF refers to high spin open shell SCF only, but other open shell coupling cases are possible using the GVB code. Analytic gradients are implemented for every possible HF type calculation possible in GAMESS, and therefore numerical hessians are available for each. Analytic hessian calculation is implemented for RHF, ROHF, and any GVB case with NPAIR=0 or NPAIR=1. Analytic hessians are more accurate, and much more quickly computed than numerical hessians, but require additional disk storage to perform an integral transformation, and also more physical memory. The second order Moller-Plesset energy correction (MP2) is implemented for RHF, UHF, ROHF, and MCSCF wavefunctions. Analytic gradients may be obtained for MP2 with RHF, UHF, or ROHF reference wavefunctions, and MP2 level properties are therefore available for these, see MP2PRP in $MP2. All other cases give properties for the SCF function. Direct SCF is implemented for every possible HF type calculation. The direct SCF method may not be used with DEM convergence. Direct SCF may be used during energy, gradient, numerical or analytic hessian, CI or MP2 energy correction, or localized orbital computations. direct SCF Normally, HF calculations proceed by evaluating a large number of two electron repulsion integrals, and storing these on a disk. This integral file is read in once during each HF iteration to form the appropriate Fock operators. In a direct HF, the integrals are not stored on disk, but are instead reevaluated during each HF iteration. Since the direct approach *always* requires more CPU time, the default for DIRSCF in $SCF is .FALSE. Even though direct SCF is slower, there are at least two reasons why you may want to consider using it. The first is that it may not be possible to store all of the integrals on the disk drives attached to your computer. Second, what you are really interested in is reducing the wall clock time to obtain your answer, not the CPU time. Workstations, particularly nodes with multiple CPUs and only one disk subsystem, may have modest hardware I/O capabilities. Other environments such as a mainframe shared by many users may also have very poor CPU/wall clock performance for I/O bound jobs such as conventional HF. You can estimate the disk storage requirements for conventional HF using a P or PK file by the following formulae: nint = 1/sigma * 1/8 * N**4 Mbytes = nint * x / 1024**2 Here N is the total number of basis functions in your run, which you can learn from an EXETYP=CHECK run. The 1/8 accounts for permutational symmetry within the integrals. Sigma accounts for the point group symmetry, and is difficult to estimate accurately. Sigma cannot be smaller than 1, in no symmetry (C1) calculations. For benzene, sigma would be almost six, since you generate 6 C's and 6 H's by entering only 1 of each in $DATA. For water sigma is not much larger than one, since most of the basis set is on the unique oxygen, and the C2v symmetry applies only to the H atoms. The factor x is 12 bytes per integral for basis sets smaller than 255, and 16 otherwise. Finally, since integrals that are very close to zero need not be stored on disk, the actual power dependence is not as bad as N**4, and in fact in the limit of very large molecules can be as low as N**2. Thus plugging in sigma=1 should give you an upper bound to the actual disk space needed. If the estimate exceeds your available disk storage, your only recourse is direct HF. What are the economics of direct HF? Naively, if we assume the run takes 10 iterations to converge, we must spend 10 times more CPU time computing the integrals on each iteration. However, we do not have to waste any CPU time reading blocks of integrals from disk, or in unpacking their indices. We also do not have to waste any wall clock time waiting for a relatively slow mechanical device such as a disk to give us our data. There are some less obvious savings too, as first noted by Almlof. First, since the density matrix is known while we are computing integrals, we can use the Schwarz inequality to avoid doing some of the integrals. In a conventional SCF this inequality is used to avoid doing small integrals. In a direct SCF it can be used to avoid doing integrals whose contribution to the Fock matrix is small (density times integral=small). Secondly, we can form the Fock matrix by calculating only its change since the previous iteration. The contributions to the change in the Fock matrix are equal to the change in the density times the integrals. Since the change in the density goes to zero as the run converges, we can use the Schwarz screening to avoid more and more integrals as the calculation progresses. The input option FDIFF in $SCF selects formation of the Fock operator by computing only its change from iteration to iteration. The FDIFF option is not implemented for GVB since there are too many density matrices from the previous iteration to store, but is the default for direct RHF, ROHF, and UHF. So, in our hypothetical 10 iteration case, we do not spend as much as 10 times more time in integral evaluation. Additionally, the run as a whole will not slow down by whatever factor the integral time is increased. A direct run spends no additional time summing integrals into the Fock operators, and no additional time in the Fock diagonalizations. So, generally speaking, a RHF run with 10-15 iterations will slow down by a factor of 2-4 times when run in direct mode. The energy gradient time is unchanged by direct HF, and this is a large time compared to HF energy, so geometry optimizations will be slowed down even less. This is really the converse of Amdahl's law: if you slow down only one portion of a program by a large amount, the entire program slows down by a much smaller factor. To make this concrete, here are some times for GAMESS for a job which is a RHF energy for a SbC4O2NH4. These timings were obtained an extremely long time ago, on a DECstation 3100 under Ultrix 3.1, which was running only these tests, so that the wall clock times are meaningful. This system is typical of Unix workstations in that it uses SCSI disks, and the operating system is not terribly good at disk I/O. By default GAMESS stores the integrals on disk in the form of a P supermatrix, because this will save time later in the SCF cycles. By choosing NOPK=1 in $INTGRL, an ordinary integral file can be used, which typically contains many fewer integrals, but takes more CPU time in the SCF. Because the DECstation is not terribly good at I/O, the wall clock time for the ordinary integral file is actually less than when the supermatrix is used, even though the CPU time is longer. The run takes 13 iterations to converge, the times are in seconds. P supermatrix ordinary file # nonzero integrals 8,244,129 6,125,653 # blocks skipped 55,841 55,841 CPU time (ints) 709 636 CPU time (SCF) 1289 1472 CPU time (job total) 2123 2233 wall time (job total) 3468 3200 When the same calculation is run in direct mode (integrals are processed like an ordinary integral disk file when running direct), iteration 1: FDIFF=.TRUE. FDIFF=.FALSE. # nonzero integrals 6,117,416 6,117,416 # blocks skipped 60,208 60,208 iteration 13: # nonzero integrals 3,709,733 6,122,912 # blocks skipped 105,278 59,415 CPU time (job total) 6719 7851 wall time (job total) 6764 7886 Note that elimination of the disk I/O dramatically increases the CPU/wall efficiency. Here's the bottom line on direct HF: best direct CPU / best disk CPU = 6719/2123 = 3.2 best direct wall/ best disk wall= 6764/3200 = 2.1 Direct SCF is slower than conventional disk SCF, but not outrageously so! From the data in the tables, we can see that the best direct method spends about 6719-1472 = 5247 seconds doing integrals. This is an increase of about 5247/636 = 8.2 in the time spent doing integrals, in a run that does 13 iterations (13 times evaluating integrals). 8.2 is less than 13 because the run avoids all CPU charges related to I/O, and makes efficient use of the Schwarz inequality to avoid doing many of the integrals in its final iterations. convergence accelerators Generally speaking, the simpler the HF function, the better its convergence. In our experience, the majority of RHF and ROHF runs converge readily from GUESS=HUCKEL. UHF often takes considerably more iterations than either of these, due to the extremely common occurrence of heavy spin contamination. GVB runs typically require GUESS=MOREAD, although the Huckel guess usually works for NPAIR=0. GVB cases with NPAIR greater than one are particularly difficult. Unfortunately, not all HF runs converge readily. The best way to improve your convergence is to provide better starting orbitals! In many cases, this means to MOREAD orbitals from some simpler HF case. For example, if you want to do a doublet ROHF, and the HUCKEL guess does not seem to converge, do this: Do an RHF on the +1 cation. RHF is typically more stable than ROHF, UHF, or GVB, and cations are usually readily convergent. Then MOREAD the cation's orbitals into the neutral calculation which you wanted to do at first. GUESS=HUCKEL does not always start with the correct electronic configuration. It may be useful to use PRTMO in $GUESS during a CHECK run to examine the starting orbitals, and then reorder them with NORDER if that seems appropriate. Of course, by default GAMESS uses the convergence procedures which are usually most effective. Still, there are cases which are difficult, so the $SCF group permits you to select several alternative methods for improving convergence. Briefly, these are EXTRAP. This extrapolates the three previous Fock matrices, in an attempt to jump ahead a bit faster. This is the most powerful of the old-fashioned accelerators, and normally should be used at the beginning of any SCF run. When an extrapolation occurs, the counter at the left of the SCF printout is set to zero. DAMP. This damps the oscillations between several successive Fock matrices. It may help when the energy is seen to oscillate wildly. Thinking about which orbitals should be occupied initially may be an even better way to avoid oscillatory behaviour. SHIFT. Level shifting moves the diagonal elements of the virtual part of the Fock matrix up, in an attempt to uncouple the unoccupied orbitals from the occupied ones. At convergence, this has no effect on the orbitals, just their orbital energies, but will produce different (and hopefully better) orbitals during the iterations. RSTRCT. This limits mixing of the occupied orbitals with the empty ones, especially the flipping of the HOMO and LUMO to produce undesired electronic configurations or states. This should be used with caution, as it makes it very easy to converge on incorrectly occupied electronic configurations, especially if DIIS is also used. If you use this, be sure to check your final orbital energies to see if they are sensible. A lower energy for an unoccupied orbital than for one of the occupied ones is a sure sign of problems. DIIS. Direct Inversion in the Iterative Subspace is a modern method, due to Pulay, using stored error and Fock matrices from a large number of previous iterations to interpolate an improved Fock matrix. This method was developed to improve the convergence at the final stages of the SCF process, but turns out to be quite powerful at forcing convergence in the initial stages of SCF as well. By giving ETHRSH as 10.0 in $SCF, you can practically guarantee that DIIS will be in effect from the first iteration. The default is set up to do a few iterations with conventional methods (extrapolation) before engaging DIIS. This is because DIIS can sometimes converge to solutions of the SCF equations that do not have the lowest possible energy. For example, the 3-A-2 small angle state of SiLi2 (see M.S.Gordon and M.W.Schmidt, Chem.Phys.Lett., 132, 294-8(1986)) will readily converge with DIIS to a solution with a reasonable S**2, and an energy about 25 milliHartree above the correct answer. A SURE SIGN OF TROUBLE WITH DIIS IS WHEN THE ENERGY RISES TO ITS FINAL VALUE. However, if you obtain orbitals at one point on a PES without DIIS, the subsequent use of DIIS with MOREAD will probably not introduce any problems. Because DIIS is quite powerful, EXTRAP, DAMP, and SHIFT are all turned off once DIIS begins to work. DEM and RSTRCT will still be in use, however. SOSCF. Approximate second-order (quasi-Newton) SCF orbital optimization. SOSCF will converge about as well as DIIS at the initial geometry, and slightly better at subsequent geometries. There's a bit less work solving the SCF equations, too. The method kicks in after the orbital gradient falls below SOGTOL. Some systems, particularly transition metals with ECP basis sets, may have Huckel orbitals for which the gradient is much larger than SOGTOL. In this case it is probably better to use DIIS instead, with a large ETHRSH, rather than increasing SOGTOL, since you may well be outside the quadratic convergence region. SOSCF does not exhibit true second order convergence since it uses an approximation to the inverse hessian. SOSCF will work for MOPAC runs, but is slower in this case. SOSCF will work for UHF, but its convergence may be better than DIIS. SOSCF will work for non-Abelian cases, but may encounter problems if the open shell is degenerate. It should be clear that SOSCF and DIIS are the two work-horse convergers, with DAMP (and possibly SHIFT) useful in cases where the initial guess is such that these two are not engaged immediately. If you compute many different types of molecules, you will find cases where SOSCF works but DIIS does not, but also cases where DIIS works but SOSCF does not (although often both will work). If you do not obtain convergence with one of these, try the other one! If you still have problems, attempt to get better starting orbitals. DEM. Direct energy minimization should be your last recourse. It explores the "line" between the current orbitals and those generated by a conventional change in the orbitals, looking for the minimum energy on that line. DEM should always lower the energy on every iteration, but is very time consuming, since each of the points considered on the line search requires evaluation of a Fock operator. DEM will be skipped once the density change falls below DEMCUT, as the other methods should then be able to affect final convergence. While DEM is working, RSTRCT is held to be true, regardless of the input choice for RSTRCT. Because of this, it behooves you to be sure that the initial guess is occupying the desired orbitals. DEM is available only for RHF. The implementation in GAMESS resembles that of R.Seeger and J.A.Pople, J.Chem.Phys. 65, 265-271(1976). Simultaneous use of DEM and DIIS resembles the ADEM-DIOS method of H.Sellers, Chem.Phys.Lett. 180, 461-465(1991). DEM does not work with direct SCF. high spin open shell SCF (ROHF) Open shell SCF calculations are performed in GAMESS by both the ROHF code and the GVB code. Note that when the GVB code is executed with no pairs, the run is NOT a true GVB run, and should be referred to in publications and discussion as a ROHF calculation. Low spin couplings are possible with the GVB program. The ROHF module in GAMESS can handle any number of open shell electrons, provided these have a high spin coupling. For example: one open shell, doublet: $CONTRL SCFTYP=ROHF MULT=2 $END two open shells, triplet: $CONTRL SCFTYP=ROHF MULT=3 $END m open shells, high spin: $CONTRL SCFTYP=ROHF MULT=m+1 $END John Montgomery (who was then at United Technologies) is responsible for the ROHF implementation in GAMESS. The following discussion is due to him, dating from 1988 when his method of forming a combined Fock operator was included in GAMESS. Other choices (Euler and two "canonical" sets) were added to the table in 2009/2010. The Fock matrix in the MO basis has the form closed open virtual closed F2 | Fb | (Fa+Fb)/2 ----------------------------------- open Fb | F1 | Fa ----------------------------------- virtual (Fa+Fb)/2 | Fa | F0 where Fa and Fb are the usual alpha and beta Fock matrices any UHF program produces. All ROHF methods agree on these, as they are the variational conditions that separate the doubly occupied, alpha occupied, and empty orbital spaces. The diagonal blocks can be written F2 = Acc*Fa + Bcc*Fb F1 = Aoo*Fa + Boo*Fb F0 = Avv*Fa + Bvv*Fb Some choices for the canonicalization coefficients to define the diagonal blocks are Acc Bcc Aoo Boo Avv Bvv Guest and Saunders 1/2 1/2 1/2 1/2 1/2 1/2 Roothaan single matrix -1/2 3/2 1/2 1/2 3/2 -1/2 Davidson/1988 1/2 1/2 1 0 1 0 Binkley, Pople, Dobosh 1/2 1/2 1 0 0 1 McWeeny and Diercksen 1/3 2/3 1/3 1/3 2/3 1/3 Faegri and Manne 1/2 1/2 1 0 1/2 1/2 GVB program/Euler 1/2 1/2 1/2 0 1/2 1/2 "canonical 1" 0 1 1 0 1 0 "canonical 2" (2S+1)/2S -1/2S 0 1 -1/2S (2S+1)/2S See below for how these last two rows connect to ionization events. The 1988 GAMESS ROHF program using a now deleted Davidson-type ROHF produced final orbitals matching the line "1988" above. This differs from the choices made in Davidson's own MELD program. The MELD program itself has always done a "cleanup" of the virtual space, after convergence, using Avv=Bvv=1/2, producing orbitals which are the same as Faegri/Manne. If MELD's MP2 option is chosen, the occupied space is also altered after convergence, using Aoo=Boo=1/2, which is the Guest/Saunders line above. Thus the term "Davidson orbitals" is used here to refer to the behavior of the now-deleted 1988 ROHF code in GAMESS, which didn't have either type of final orbital cleanup. The choice of the diagonal blocks is arbitrary, as ROHF is converged when the off diagonal blocks go to zero. The exact choice for these blocks can however have an effect on the convergence rate. This choice also affects the MO coefficients, and orbital energies, as the different choices produce different canonical orbitals within the three subspaces. All methods, however, will give identical total wavefunctions, and hence identical properties such as nuclear gradients and hessians. Some of the perturbation theories for open shell cases are defined in terms of a particular canonicalization, if so, GAMESS automatically canonicalizes after convergence so the desired orbitals and energies are given to the perturbation theory codes. The default coupling case in GAMESS is the Roothaan single matrix set. Note that pre-1988 versions of GAMESS produced "Davidson/1988" orbitals. If you would like to fool around with any of these other canonicalizations, the Acc, Aoo, Avv and Bcc, Boo, Bvv parameters can be input as the first three elements of ALPHA and BETA in $SCF. For example, the McWeeny/Diercksen canonicalization, in double precision, is obtained by $scf couple=.true. alpha(1)=0.333333333333333333, 0.333333333333333333, 0.666666666666666667 beta(1)=0.666666666666666667, 0.333333333333333333, 0.333333333333333333 $end Here is some idea of the range of eigenvalues that result from using the various canonicalization schemes in the table above. The system is 6-31G nitrogen atom, all runs give E= -54.3820511123 (matching all 10 decimals): 1s 2s 2p "3p" "3s" Roothaan -15.5514 -0.5306 -0.1774 +0.7666 +0.8704 McWeeny,Diercksen -15.6214 -0.8745 -0.1183 +0.8984 +0.9684 Davidson -15.6355 -0.9432 -0.5657 +0.8457 +0.9292 Guest,Saunders -15.6355 -0.9432 -0.1774 +0.9248 +0.9880 Binkley,Pople,Dob. -15.6355 -0.9432 -0.5657 +1.0039 +1.0469 Faegri,Manne -15.6355 -0.9432 -0.5657 +0.9248 +0.9880 GVB (aka Euler) -15.6355 -0.9432 -0.2828 +0.9248 +0.9880 "canonical 1" -15.5933 -0.7370 -0.5657 +0.8457 +0.9292 "canonical 2" -15.7065 -1.2863 +0.2109 +1.0567 +1.0861 Since the ionization potential (IP) for a 2p electron in Nitrogen is 0.53 Hartree, it is clear that most of the orbital energies above do not approximately predict this IP. Recent work by Boris Plakhutin and co-workers (see the 3 papers in the ROHF references above) leads to two sets of orbitals and eigenvalues, for the prediction of both IP and electron affinities (EA), for various ionization events, starting from a high spin ROHF state: canonical 1 (to produce high spin final states) A1 = remove beta e- from filled space, B1 = remove alpha e- from open shell, C1 = attach alpha e- to virtual space. canonical 2 (to produce low spin final states) A2 = remove alpha e- from filled space, B2 = attach beta e- to filled space, C2 = attach beta e- to virtual space. Correct handling of spin requires the value of the spin S of the initial state in the 2nd canonicalization set. Two different ROHF runs are necessary to get all six EA and IP processes. Additional discussion about ROHF orbital energies may be found in K.R.Glaesemann, M.W.Schmidt J.Phys.Chem.A 114, 8772-8777(2010) [available on-line] along with a demonstration of the non-uniqueness of ROHF- based perturbation theories. High-spin ROHF results are automatically obtained if a UHF calculation is constrained so that the occupied beta orbitals lie entirely within the occupied alpha orbital space. Such a constraint is called CUHF, and while it will appear to have different alpha/beta orbitals, and alpha/beta eigenvalues, its spin expectation value will be
=2S+1, and its energy will be exactly that of ROHF. The CUHF method is given by G.E.Scuseria, T.Tsuchimochi J.Chem.Phys. 134, 064101/1-14(2011) CUHF runs as a UHF calculation, solving two distinct SCF equations in the alpha and beta space, so its convergence behavior and its final eigenvalues will differ from that for any of the ROHF canonicalizations. See the section on perturbation theory in this chapter regarding CUHF's perturbation theory, CUMP2. other open shell SCF cases (GVB) Genuine GVB-PP "perfect pairing" runs will be discussed later in this section. First, we will consider how to do open shell SCF with the GVB part of the program, with NPAIR=0. This is an extremely powerful open shell SCF program, capable of handling low spin couplings and/or degenerate open shells, with the energy formula E = 2 sum F-i*h-ii + sum sum ALPHA-ij*J-ij + BETA-ij*Kij i i j Here F-i (meaning subscript i) is a fractional occupancy, and is 1.0 for the filled shell that is usually present below the open shell(s). A d**2 shell has F-i=0.2. The h,J,K are the usual one electron, coulomb, and exchange operators. For a few common cases, the F, ALPHA, BETA are internally stored in the program, one open shell, doublet: $CONTRL SCFTYP=GVB MULT=2 $END $SCF NCO=xx NSETO=1 NO(1)=1 $END two open shells, triplet: $CONTRL SCFTYP=GVB MULT=3 $END $SCF NCO=xx NSETO=2 NO(1)=1,1 $END two open shells, singlet: $CONTRL SCFTYP=GVB MULT=1 $END $SCF NCO=xx NSETO=2 NO(1)=1,1 $END Note that the first two cases duplicate runs which the high spin ROHF module can readily do. The last case is also an ROHF calculation, albeit for a low-spin coupling. One should note that the open shell singlet case is a variational calculation IF AND ONLY IF the two singly degenerate orbitals have different space symmetry. A publication should refer to any run that has NPAIR=0 as being ROHF type, rather than GVB, since there are no perfect pairs in use. Note, however, that the GVB program can, if you wish, have both open shells and genuine valence bond pairs at the same time! See the following section about the use of GVB's perfect pairs. The use of open shells and pairs is illustrated in one of the GAMESS standard example inputs. If you would like to do any cases other than those shown above, including many cases with orbital degeneracy, you must derive the coupling coefficients ALPHA and BETA, and input them with the occupancies F in the $SCF group. Fortunately, many cases can be looked up! Mariusz Klobukowski of the University of Alberta has shown how to obtain coupling coefficients for the GVB open shell program for many such open shell states. These can be obtained from Appendix A of the book "A General SCF Theory" by Ramon Carbo and Josep M. Riera, Springer-Verlag (1978). The rule is (1) F(i) = 1/2 * omega(i) (2) ALPHA(i) = alpha(i) (3) BETA(i) = - beta(i), where omega, alpha, beta are symbols used in these Tables. The keyword NSETO gives the number of open shells, and keyword NO gives the degeneracy of each open shell. Thus the 5-S state of carbon (from configuration s1,p3) would enter NSETO=2 and NO(1)=1,3. NCO is an easy keyword: that is the number of filled orbitals. The total number of electrons in an open shell GVB run is NE = 2*NCO + sum 2*F(i)*NO(i) i and this value will be checked against your ICHARG input. - - - - Specific input for all terms found in the atomic p**N configurations follow. Be sure to enter 14 digits, as these values are part of a double precision energy formula! Values for the excited terms in the p**N configurations were extracted from C.F.Jackels and E.R.Davidson Int. J. Quantum Chem. 8, 707-714(1974), which explains the concept of averaging over equivalent determinants to enforce a symmetric density matrix, which preserves radial symmetry in the atomic orbitals. The ALPHA and BETA values can also be extracted from Roothaan's 1960 A and B coefficients by the prescription detailed below. ! p**1 2-P state $CONTRL SCFTYP=GVB MULT=2 $END $SCF NCO=xx NSETO=1 NO=3 COUPLE=.TRUE. F(1)= 1.0 0.16666666666667 ALPHA(1)= 2.0 0.33333333333333 0.00000000000000 BETA(1)= -1.0 -0.16666666666667 -0.00000000000000 $END ! p**2 3-P state $CONTRL SCFTYP=GVB MULT=3 $END $SCF NCO=xx NSETO=1 NO=3 COUPLE=.TRUE. F(1)= 1.0 0.33333333333333 ALPHA(1)= 2.0 0.66666666666667 0.16666666666667 BETA(1)= -1.0 -0.33333333333333 -0.16666666666667 $END For the 1-D excited state, change the open-shell parameters to ALPHA(3)=0.1 and BETA(3)= +0.03333333333333 For the 1-S excited state, change the open-shell parameters to ALPHA(3)=0.0 and BETA(3)= +0.33333333333333 ! p**3 4-S state $CONTRL SCFTYP=ROHF MULT=4 $END which is equivalent to $CONTRL SCFTYP=GVB MULT=4 $END $SCF NCO=xx NSETO=1 NO=3 COUPLE=.TRUE. F(1)= 1.0 0.50000000000000 ALPHA(1)= 2.0 1.00000000000000 0.50000000000000 BETA(1)= -1.0 -0.50000000000000 -0.50000000000000 $END For 2-D, use ALPHA(3)= 0.4 and BETA(3)= -0.2 for 2-P, use ALPHA(3)= 0.33333333333333 and BETA(3)= 0.0 ! p**4 3-P state $CONTRL SCFTYP=GVB MULT=3 $END $SCF NCO=xx NSETO=1 NO=3 COUPLE=.TRUE. F(1)= 1.0 0.66666666666667 ALPHA(1)= 2.0 1.33333333333333 0.83333333333333 BETA(1)= -1.0 -0.66666666666667 -0.50000000000000 $END For 1-D, use ALPHA(3)= 0.76666666666667 and BETA(3)= -0.3 for 1-S, use ALPHA(3)= 0.66666666666667 and BETA(3)= 0.0 ! p**5 2-P state $CONTRL SCFTYP=GVB MULT=2 $END $SCF NCO=xx NSETO=1 NO=3 COUPLE=.TRUE. F(1)= 1.0 0.83333333333333 ALPHA(1)= 2.0 1.66666666666667 1.33333333333333 BETA(1)= -1.0 -0.83333333333333 -0.66666666666667 $END - - - - Coupling constants for the highest spin state(s) in d**N configurations are taken from "Handbook of Gaussian Basis Sets", R.Poirier, R.Kari, I.G.Csizmadia, Elsevier, Amsterdam, 1985. ! d**1 2-D state $CONTRL SCFTYP=GVB MULT=2 $END $SCF NCO=xx NSETO=1 NO=5 COUPLE=.TRUE. F(1)=1.0,0.1 ALPHA(1)= 2.0, 0.20, 0.00 BETA(1)=-1.0,-0.10, 0.00 $END ! d**2 average of 3-F and 3-P states $CONTRL SCFTYP=GVB MULT=3 $END $SCF NCO=xx NSETO=1 NO=5 COUPLE=.TRUE. F(1)=1.0,0.2 ALPHA(1)= 2.0, 0.40, 0.05 BETA(1)=-1.0,-0.20,-0.05 $END Note: "average" means a degeneracy weighted combination of determinants with Sz=S: here, 2 electrons in 5 orbitals coupled as a triplet is ten determinants with Sz=1. The energy from these parameters is E=[7xE(3-F)+3xE(3-P)]/10. ! d**3 average of 4-F and 4-P states $CONTRL SCFTYP=GVB MULT=4 $END $SCF NCO=xx NSETO=1 NO=5 COUPLE=.TRUE. F(1)=1.0,0.3 ALPHA(1)= 2.0, 0.60, 0.15 BETA(1)=-1.0,-0.30,-0.15 $END ! d**4 5-D state $CONTRL SCFTYP=GVB MULT=5 $END $SCF NCO=xx NSETO=1 NO=5 COUPLE=.TRUE. F(1)=1.0,0.4 ALPHA(1)= 2.0, 0.80, 0.30 BETA(1)=-1.0,-0.40,-0.30 $END ! d**5 6-S state $CONTRL SCFTYP=ROHF MULT=6 $END ! d**6 5-D state $CONTRL SCFTYP=GVB MULT=5 $END $SCF NCO=xx NSETO=1 NO=5 COUPLE=.TRUE. F(1)=1.0,0.6 ALPHA(1)= 2.0, 1.20, 0.70 BETA(1)=-1.0,-0.60,-0.50 $END ! d**7 average of 4-F and 4-P states $CONTRL SCFTYP=GVB MULT=4 $END $SCF NCO=xx NSETO=1 NO=5 COUPLE=.TRUE. F(1)=1.0,0.7 ALPHA(1)= 2.0, 1.40, 0.95 BETA(1)=-1.0,-0.70,-0.55 $END ! d**8 average of 3-F and 3-P states $CONTRL SCFTYP=GVB MULT=3 $END $SCF NCO=xx NSETO=1 NO=5 COUPLE=.TRUE. F(1)=1.0,0.8 ALPHA(1)= 2.0, 1.60, 1.25 beta(1)=-1.0,-0.80,-0.65 $end ! d**9 2-D state $CONTRL SCFTYP=GVB MULT=2 $END $SCF NCO=xx NSETO=1 NO=5 COUPLE=.TRUE. F(1)=1.0,0.9 ALPHA(1)= 2.0, 1.80, 1.60 BETA(1)=-1.0,-0.90,-0.80 $END Note that GAMESS can do a proper calculation on the ground term for the d**2, d**3, d**7, d**8 configurations only by means of state averaged MCSCF. For d**8, use $contrl scftyp=mcscf mult=3 $end $det group=c1 ncore=xx nact=5 nels=8 nstate=10 wstate(1)=1,1,1,1,1,1,1,0,0,0 $end to correctly average the 7 lowest roots (3-F) with no weight given to the highest three roots (3-P). Although this is done with the SCFTYP=MCSCF program, this is still a SCF calculation since only the orbitals, but not any CI coefficients, are optimized. Dirk Andrae in Berlin has provided a great many examples of using the MCSCF program to do terms found in the p**n, d**n, and even f**n configurations: http://userpage.fu-berlin.de/~dandrae/ openshell/openls/openls.html Type that web reference on just one line, of course. Open shell cases such as s**1,d**n are probably most easily tackled with the state-averaged MCSCF program. The ORMAS CI code may be convenient in fixing the number of electrons found in each open shell. If you are a true afficionado of atomic calculations, including open shell f configurations, look for Russ Pitzer's version of the famous Roothaan/Bagus ATOM-SCF program R.M.Pitzer, Comput.Phys.Commun. 170, 239-264(2005) - - - - The literature contains an alternate energy formula, involving the so-called Roothaan A and B parameters: E = 2 sum f-i*h-ii + sum sum f-i*f-j[2Aij*Jij - Bij*Kij] i i j Comparing this to the GVB program's energy formula above, we immediately see that its fractional occupancy f is the same as GVB's F, and that ALPHA-ij = +2 * f-i * f-j * A-ij BETA-ij = -1 * f-i * f-j * B-ij Let c=the closed shell orbitals, and o=a single open shell. Tables such as Roothaan's 1960 paper give only the A-oo and B-oo values, since A-cc = B-cc = A-co = B-co = 1.0. It is easy to convert such A,B values to all ALPHA-cc,-co,-oo and BETA-cc,-co,-oo values. Degenerate open shells in linear molecules follow immediately from Roothaan's 1960 table: ! 2-Pi from configuration pi**1: $contrl scftyp=gvb mult=2 $end $scf nco=xx nseto=1 no=2 couple=.true. f(1)=1.0,0.25 alpha(1)= 2.0, 0.50, 0.0 beta(1)=-1.0,-0.25, 0.0 $end ! 3-Sigma-minus from configuration pi**2: $contrl scftyp=gvb mult=3 $end $scf nco=xx nseto=1 no=2 couple=.true. f(1)=1.0,0.50 alpha(1)= 2.00, 1.00, 0.50 beta(1)=-1.00,-0.50,-0.50 $end ! 1-Delta from configuration pi**2: $contrl scftyp=gvb mult=1 $end $scf nco=xx nseto=1 no=2 couple=.true. f(1)=1.0,0.50 alpha(1)= 2.00, 1.00, 0.00 beta(1)=-1.00,-0.50,+0.50 $end ! 1-Sigma-plus from configuration pi**2: $contrl scftyp=gvb mult=1 $end $scf nco=xx nseto=1 no=2 couple=.true. f(1)=1.0,0.50 alpha(1)= 2.00, 1.00, 0.25 beta(1)=-1.00,-0.50, 0.00 $end ! 2-Pi from configuration pi**3: $contrl scftyp=gvb mult=2 $end $scf nco=xx nseto=1 no=2 couple=.true. f(1)=1.0,0.75 alpha(1)= 2.00, 1.50, 1.00 beta(1)=-1.00,-0.75,-0.50 $end The same inputs apply to delta**N configurations: only the names of the energy terms (states) change! It is possible to do every term arising from the s**1,p**N configurations. Roothaan-style A,B parameters are as follows: Acc = Bcc = 1.0 (as always) Aco = Bco = 1.0, for o=s and for o=p (as always) Ass = Bss = 0.0, since there is no other s e- to repel Asp = 1.0 Bsp = -3K101 + 1 App and Bpp: take from the parent term in config p**N and then apply the rule turning A,B into ALPHA,BETA, for Fc=1, Fs=0.5, Fp=N/6. K101 should be taken from Roothaan and Bagus in Methods Comput. Phys. 2, 49(1963). - - - - It may be useful to obtain a "configuration average", especially whenever a specific term energy cannot be computed. Boris Plakhutin in Novosibirsk has kindly provided a recipe for the Roothaan A,B values that give the average energy of configuration O**N, meaning N electrons in some open shell O. Let D be the degeneracy of the open shell O: atoms: O=p, d, f... has D= 3, 5, 7... non-linear molecules: O=e, t, g, h has D= 2, 3, 4, 5. The fractional occupation number F-o = N/(2*D). The open shell Roothaan couplings are A-oo = B-oo = (N-1)/(N-F). For example, the case d**2 has F-o=1/5 so A-oo=B-oo=5/9. Using the formulae to generate all ALPHA and BETA causes the GVB program to optimize the following spin-and-space weighted average of all terms appearing in d**2: E-avg = [21xE(3-F) + 9xE(3-P) + 9xE(1-G) + 5xE(1-D) + E(1-S)]/45 by this input $contrl scftyp=gvb mult=3 $end $scf nco=xx nseto=1 no=5 couple=.true. f(1)= 1.0, 0.2 alpha(1)= 2.00, 0.40, 0.044444444444444444 beta(1)=-1.00,-0.20,-0.022222222222222222 $end Note that here, as in all other runs with the GVB module, the spin in $CONTRL is irrelevant: the energy that is optimized is determined by F, ALPHA, and BETA. In this last case, these data implicitly include both triplet and singlet determinants in the energy averaging. Boris Plakhutin also provided a formula in case you want to average over a particular spin S within an open shell O**N with some degeneracy D implying a fractional occupancy of F. The result is N**2(N-1) + FxN(N-4) + 4FxS(S+1) A-oo = -------------------------------- N(N**2-4F**2) and 4F(N-1) + N(N-4) + 4S(S+1) B-oo = -------------------------- (N**2-4F**2) For the d**2 configuration, the space-and-spin degeneracy weighted average of the 3-F and 3-P terms has Aoo=5/8 and Boo=10/8, while the average of the 1-G, 1-D, and 1-S terms has Aoo=5/12 and Boo=-10/12. These optimize orbitals for the average E-triplet = [21xE(3-F) + 9xE(3-P)]/30 which is not quite the same as the space degeneracy average shown above, and for the average E-singlet = [9xE(1-G) + 5xE(1-D) + E(1-S)]/15 The literature reference for the overall and the spin- specific configurational averages is B.N.Plakhutin, J. Mathematical Chem. 22, 203-233(1997) true GVB perfect pairing runs True GVB runs are obtained by choosing NPAIR nonzero. If you wish to have some open shell electrons in addition to the geminal pairs, you may add the pairs to the end of any of the GVB coupling cases shown above. The GVB module assumes that you have reordered your MOs into the order: NCO double occupied orbitals, NSETO sets of open shell orbitals, and NPAIR sets of geminals (with NORDER=1 in the $GUESS group). Each geminal consists of two orbitals and contains two singlet coupled electrons (perfect pairing). The first MO of a geminal is probably heavily occupied (such as a bonding MO u), and the second is probably weakly occupied (such as an antibonding, correlating orbital v). If you have more than one pair, you must be careful that the initial MOs are ordered u1, v1, u2, v2..., which is -NOT- the same order that RHF starting orbitals will be found in. Use NORDER=1 to get the correct order. These pair wavefunctions are actually a limited form of MCSCF. GVB runs are much faster than MCSCF runs, because the natural orbital u,v form of the wavefunction permits a Fock operator based optimization. However, convergence of the GVB run is by no means assured. The same care in selecting the correlating orbitals that you would apply to an MCSCF run must also be used for GVB runs. In particular, look at the orbital expansions when choosing the starting orbitals, and check them again after the run converges. GVB runs will be carried out entirely in orthonormal natural u,v form, with strong orthogonality enforced on the geminals. Orthogonal orbitals will pervade your thinking in both initial orbital selection, and the entire orbital optimization phase (the CICOEF values give the weights of the u,v orbitals in each geminal). However, once the calculation is converged, the program will generate and print the nonorthogonal, generalized valence bond orbitals. These GVB orbitals are an entirely equivalent way of presenting the wavefunction, but are generated only after the fact. Convergence of true GVB runs is by no means as certain as convergence of RHF, UHF, ROHF, or GVB with NPAIR=0. You can assist convergence by doing a preliminary RHF or ROHF calculation, and use these orbitals for GUESS=MOREAD. Few, if any, GVB runs with NPAIR non-zero will converge without using GUESS=MOREAD. Generation of MVOs during the preliminary SCF can also be advantageous. In fact, all the advice outlined for MCSCF computations below is germane, for GVB-PP is a type of MCSCF computation. The total number of electrons in the GVB wavefunction is given by the following formula: NE = 2*NCO + sum 2*F(i)*NO(i) + 2*NPAIR i The charge is obtained by subtracting the total number of protons given in $DATA. The multiplicity is implicit in the choice of alpha and beta constants. Note that ICHARG and MULT must be given correctly in $CONTRL anyway, as the number of electrons from this formula is double checked against the ICHARG value. the special case of TCSCF The wavefunction with NSETO=0 and NPAIR=1 is called GVB-PP(1) by Goddard, two configuration SCF (TCSCF) by Schaefer or Davidson, and CAS-SCF with two electrons in two orbitals by others. Note that this is just semantics, as these are identical. This is a very important type of wavefunction, as TCSCF is the minimum acceptable treatment for singlet biradicals. The TCSCF wavefunction can be obtained with SCFTYP=MCSCF, but it is usually much faster to use the Fock based SCFTYP=GVB. Because of its importance, the TCSCF function (together with open shells, if desired) permits analytic hessian computation. a caution about symmetry Caution! Some exotic calculations with the GVB program do not permit the use of symmetry. The symmetry algorithm in GAMESS was "derived assuming that the electronic charge density transforms according to the completely symmetric representation of the point group", Dupuis/King, JCP, 68, 3998(1978). This may not be true for certain open shell cases, and in fact during GVB runs, it may not be true for closed shell singlet cases! First, consider the following correct input for the singlet-delta state of NH: $CONTRL SCFTYP=GVB NOSYM=1 $END $SCF NCO=3 NSETO=2 NO(1)=1,1 $END for the x**1y**1 state, or for the x**2-y**2 state, $CONTRL SCFTYP=GVB NOSYM=1 $END $SCF NCO=3 NPAIR=1 CICOEF(1)=0.707,-0.707 $END Neither gives correct results, unless you enter NOSYM=1. The electronic term symbol is degenerate, a good tip off that symmetry cannot be used. However, some degenerate states can still use symmetry, because they use coupling constants averaged over all degenerate states within a single term, as is done in EXAM15 and EXAM16. Here the "state averaged SCF" leads to a charge density which is symmetric, and these runs can exploit symmetry. Secondly, since GVB runs exploit symmetry for each of the "shells", or type of orbitals, some calculations on totally symmetric states may not be able to use symmetry. An example is CO or N2, using a three pair GVB to treat the sigma and pi bonds. Individual configurations such as (sigma)**2,(pi-x)**2,(pi-y*)**2 do not have symmetric charge densities since neither the pi nor pi* level is completely filled. Correct answers for the sigma-plus ground states result only if you input NOSYM=1. Problems of the type mentioned should not arise if the point group is Abelian, but will be fairly common in linear molecules. Since GAMESS cannot detect that the GVB electronic state is not totally symmetric (or averaged to at least have a totally symmetric density), it is left up to you to decide when to input NOSYM=1. If you have any question about the use of symmetry, try it both ways. If you get the same energy, both ways, it remains valid to use symmetry to speed up your run. And beware! Brain dead computations, such as RHF on singlet O2, which actually is a half filled degenerate shell, violate the symmetry assumptions, and also violate nature. Use of partially filled degenerate shells always leads to very wild oscillations in the RHF energy, which is how the program tries to tell you to think first, and compute second. Configurations such as pi**2, e**1, or f2u**4 can be treated, but require GVB wavefunctions and F, ALPHA, BETA values from the sources mentioned.