How to do MCSCF (and CI) calculations Multi-configuration self consistent field (MCSCF) wavefunctions are the most general SCF possible. MCSCF allows for a natural description of chemical processes involving the separation of electrons (bond breaking, electronic excitation, etc), which are often not well represented using the single configuration SCF methods. MCSCF wavefunctions, as the name implies, contain more than one configuration, each of which is multiplied by a configuration interaction (CI) coefficient, optimized to determine its weight. In addition, the shapes of the orbitals used to form each of the configurations are optimized, just as in a simpler SCF, to self consistency. Typically every different chemical problem requires that an MCSCF wavefunction be designed to treat it, on a case by case basis, by choosing an "active space". For example, one may be interested in describing the reactivity of a particular functional group, instead of elsewhere in the molecule. So, the active electrons and active orbitals will be those that are "active" on that functional group. Orbitals elsewhere in the molecule just remain doubly occupied, as for RHF. This means some attention must be paid to orbitals in order to obtain the desired results. Procedures for the selection of configurations (which amounts to choosing the number of active electrons and active orbitals), for the two mathematical optimizations just mentioned, ways to interpret the resulting MCSCF wavefunction, and treatments for dynamical electron correlation of MCSCF wavefunctions are the focus of a review article: "The Construction and Interpretation of MCSCF wavefunctions" M.W.Schmidt and M.S.Gordon, Annu.Rev.Phys.Chem. 49,233-266(1998) One section of this article is devoted to the problem of designing the correct active space to treat your problem. Additional reading is listed at the end of this section. A much newer and very valuable review article is "Multiconfiguration self-consistent field and multireference configuration interaction methods and applications" P.G.Szalay, T.Muller, G.Gidofalvi, H.Lischka, R.Shepard Chem.Rev. 112, 108-181(2012) These pages describe a powerful and mature MCSCF program, allowing computation of the MCSCF energy, nuclear gradient, and nuclear hessian for pure states. Practical procedures for generating starting orbitals are available. Localized orbital analysis of the final active orbitals is provided. State-averaged energies and their analytic gradients can be obtained. Minimum energy crossing points and conical intersections between surfaces may be found, see RUNTYP=MEX and CONINT. Non-adiabatic coupling matrix elements (NACME) can be computed, see RUNTYP=NACME. Diabatic potential energy surfaces can be obtained, see DIABAT in $MCSCF. If desired, spin-orbit couplings or transition dipole moments can be found, see elsewhere in this chapter. Efficient perturbative treatments of the dynamical correlation energy for all electrons, whether in active and filled orbitals, are provided. Of course, parallel computation has been enabled. The most efficient technique implemented in GAMESS for finding the dynamic correlation energy of MCSCF is second order perturbation theory, in the variant known as MCQDPT (known as MRMP for one state). MCQDPT is discussed in a different section of this chapter. The use of CI, probably in the form of second order CI, will be described below, en passant, during discussion of the input defining the configurations for MCSCF. Selection of a CI following any type of SCF (except UHF) is made with CITYP in the $CONTRL group, and masterminded by $CIINP. MCSCF implementation With the exception of the QUAD converger, the MCSCF program is of the type termed "unfolded two-step" by Roos. This means the orbital and CI coefficient optimizations are separated. The latter are obtained in a conventional CI diagonalization, while the former are optimized by a separate orbital improvement step. Each MCSCF iteration (except for the JACOBI and QUAD convergers) consists of the following steps: 1) transformation of AO integrals to the current MO basis, 2) generation of the Hamiltonian matrix and optimization of the CI coefficients by a Davidson diagonalization, 3) generation of the first and second order density matrix, 4) improvement of the molecular orbitals. During the first iteration at the first geometry, you will receive verbose output from each of these steps, but each subsequent iteration produce only a single summary line. The CI problem in steps two and three has four options for the many electron basis, namely ALDET, ORMAS, or GENCI using determinants, or GUGA using CSFs. This choice is made with the keyword CISTEP in $MCSCF. Much more will be said below about the differences between determinants and CSFs. The word "configuration" will be used throughout this section to refer to either determinants or CSFs, when a generic term is needed for the many-electron functions. Most people use CSF and configuration interchangeably, so please note the distinction made here. The orbital update in step four has five options, namely FOCAS, SOSCF, FULLNR, JACOBI, and QUAD, listed here in roughly the order of their increasing mathematical sophistication, convergence characteristics, and of course, their computer resource requirements. Again, these are chosen by keywords in the $MCSCF group. More will be said just below about the relative merits of these. Depending on the converger chosen, the program will select the appropriate kind of integral transformation at step one. There's seldom need to try to fine tune this, but note that the $TRANS group allows you to choose an AO integral direct transformation, with the DIRTRF flag. The type of CI and the type of orbital converger are to some extent "mix and match". This is particularly true for the two full CI programs, GUGA or ALDET, where either produces exactly the same CI density matrices. Here is a chart of the ways to combine CI and orbital optimizers: parallel run's orbital transformation CI computation via CISTEP= converger memory GUGA ALDET GENCI ORMAS --------- -------------- ---- ----- ----- ----- FOCAS replicated ok ok silly silly SOSCF replicated ok ok ok ok FULLNR distributed ok ok ok ok QUAD serial ok xx xx xx JACOBI serial ok ok ok ok "xx" means QUAD converger is coded only for CISTEP=GUGA. "silly" means that this converger ignores active-active rotations, so these runs are likely to be divergent, or perhaps converge to a false solution. "serial" means this can only run sequentially at present. The next two sections provide more information on the two mathematical optimizations, first how the orbital shape is refined, and then the determinantion of CI coefficients. orbital updates There are presently five orbital improvement options, namely FOCAS, SOSCF, FULLNR, JACOBI, and QUAD. All but the JACOBI update run in parallel. Each converger is discussed briefly below, in order of increasing robustness. The most commonly used convergers are SOSCF and FULLNR. The input to control the orbital update step is the $MCSCF group, where you can pick the convergence procedure. Most of the input in this group is rather specialized, but note in particular MAXIT and ACURCY, which control the convergence behavior. FOCAS is a first order, complete active space MCSCF optimization procedure. It is based on a novel approach due to Meier and Staemmler, using very fast but numerous microiterations to improve the convergence of what is intrinsically a first order method. Since FOCAS requires only one virtual orbital in the integral transformation to compute the Lagrangian (whose asymmetry is the orbital gradient, and must fall below ACURCY at convergence), the total MCSCF job may take less time than a second order method, even though it may require many more iterations to converge. The use of microiterations is crucial to FOCAS' ability to converge. It is important to take a great deal of care choosing the starting orbitals. SOSCF is a method built upon the FOCAS code, which seeks to combine the speed of FOCAS with approximate second order convergence properties. Thus SOSCF is an approximate Newton-Raphson, based on a diagonal guess at the orbital hessian, and in fact has much in common with the SOSCF option in $SCF. Its time requirements per iteration are like FOCAS, with a convergence rate better than FOCAS but not as good as true second order. Storage of only the diagonal of the orbital hessian allows the SOSCF method to be used with much larger basis sets than exact second order methods. Because SOSCF usually requires the least CPU time, disk space, and memory needs, it is the default. Good convergence by the SOSCF method requires that you prepare starting orbitals carefully, and read in all MOs in $VEC, as providing canonicalized virtual orbitals increases the diagonal dominance of the orbital hessian. Parallel computations are possible with SOSCF, but only to a modest number of nodes. FULLNR means a full Newton-Raphson orbital improvement step is taken, using the exact orbital hessian. FULLNR is a robust convergence method, and normally takes the fewest iterations to converge. Computing the exact orbital hessian requires two virtual orbital indices be included in the integral transformation, making this step quite time consuming, and of course memory for storage of the orbital hessian must be available. Because both the transformation and orbital improvement steps of FULLNR are time consuming, FULLNR is not the default. You may want to try FULLNR when convergence is difficult, assuming you have already tried preparing good starting orbitals by the hints below. The serial FULLNR code uses the augmented hessian matrix approach to solve the Newton-Raphson equations. There are two suboptions for computation of the orbital hessian: DM2 is faster, but takes more memory than TEI. The parallel implementation of FULLNR avoids explicit storage of the orbital hessian, by recomputing the product of the hessian times orbital rotation vector during the subiterations solving the Newton-Raphson problem. The partial integral transformation used to set up the FULLNR converger has been changed to use distributed memory, and will scale like the MP2 energy/gradient programs, to many nodes. Parallel FULLNR requires large MEMORY only for the CI step (if the active space is big), but always requires a large MEMDDI. The parallel FULLNR program is essentially diskless, apart from storage of converged CI vectors. The JACOBI method uses a series of 2 by 2 orbital rotations by an angle predicted to lower the energy. This should essentially ensure convergence after sweeping through all possible orbital pairs enough times. The procedure was created to converge selected (general) determinant MCSCF functions, but of course it can be used will full lists as well in difficult cases. The JACOBI calculation will consist of a full four index transformation over all MOs before the iterations begin. Each iteration consists of 1. a small 4 index transformation over active orbitals 2. optimization of the CI vector 3. generation of the 1e- and 2e- density matrices 4. sweeps over Jacobi rotations, using MO integrals in memory to generate each rotation, with a subsequent update after each pair is rotated. 5. when sufficient energy lowering has been achieved, begin a new iteration. This procedure never generates the orbital Lagrangian! Unfortunately this means that at present it is not possible to compute nuclear gradients. Due to lack of a Lagrangian, ACURCY is of course irrelevant, so the convergence test is on ENGTOL. QUAD uses a fully quadratic, or second order approach and is thus the most powerful MCSCF converger. The QUAD code is programmed only for CISTEP=GUGA. QUAD runs begin with normal unfolded FULLNR iterations, until the orbitals approach convergence sufficiently. QUAD then begins the simultaneous optimization of CI coefficients and orbitals, and convergence should be obtained in 3-4 additional MCSCF iterations. The QUAD method requires building the full electronic hessian, including orbital/orbital, orbital/CI, and CI/CI blocks, which is a rather big matrix. In principle, this is the most robust method available, but it is limited to perhaps 50-100 CSFs only, because it is a memory hog. QUAD may be helpful in converging excited electronic states, but note that you may not use state averaging with QUAD. In practice, QUAD has not received very much use compared to the unfolded convergers. CI coefficient optimization Determinants or configuration state functions (CSFs) may be used to form the many electron basis set. It is necessary to explain these in a bit of detail so that you can understand the advantages of each. A determinant is a simple object: a product of spin orbitals with a given Sz quantum number, that is, the number of alpha spins and number of beta spins are a constant. Matrix elements involving determinants are correspondingly simple, but unfortunately determinants are not necessarily eigenfunctions of the S**2 operator. To expand on this point, consider the four familiar 2e- functions which satisfy the Pauli principle. Here u, v are space orbitals, and a, b are the alpha and beta spin functions. As you know, the singlet and triplets are: S1 = (uv + vu)/sqrt(2) * (ab - ba)/sqrt(2) T1 = (uv - vu)/sqrt(2) * aa T2 = (uv - vu)/sqrt(2) * (ab + ba)/sqrt(2) T3 = (uv - vu)/sqrt(2) * bb It is a simple matter to multiply out S1 and T2, and to expand the two determinants which have Sz=0, D1 = |ua vb| D2 = |va ub| This reveals that S1 = (D1+D2)/sqrt(2) or D1 = (S1 + T2)/sqrt(2) T2 = (D1-D2)/sqrt(2) D2 = (S1 - T2)/sqrt(2) Thus, one must take a linear combination of determinants in order to have a wavefunction with the desired total spin. There are two important points to note: a) A two by two Hamiltonian matrix over D1 and D2 has eigenfunctions with -different- spins, S=0 and S=1. b) use of all determinants with Sz=0 does allow for the construction of spin adapted states. D1+D2, or D1-D2, are -not- spin contaminated. By itself, a determinant such as D1 is said to be "spin contaminated", being a fifty-fifty admixture of singlet and triplet. (It is curious that calculations with just one such determinant are often called "singlet UHF", when this is half triplet!). Of course, some determinants are spin adapted all by themselves, for example the spin adapted functions T1 and T3 above are single determinants, as are the closed shells S2 = (uu) * (ab - ba)/sqrt(2). S3 = (vv) * (ab - ba)/sqrt(2). It is possible to perform a triplet calculation, with no singlet states present, by choosing determinants with Sz=1 such as T1, since then no state with Sz=0 exists in the determinant basis set (as is required when S=0). To summarize, the eigenfunctions of a Hamiltonian formed by determinants with any particular Sz will be spin states with S=Sz, S=Sz+1, S=Sz+2, ... but will not contain any S values smaller than Sz. CSFs are an antisymmetrized combination of a space orbital product, and a spin adapted linear combination of simple alpha-beta products. Namely, the following CSF C1 = A (uv) * (ab-ba)/sqrt(2) which has a singlet spin function is identical to S1 above if you write out what the antisymmetrizer A does, and the CSFs C2 = A (uv) * aa C3 = A (uv - vu)/sqrt(2) * (ab + ba)/sqrt(2) C4 = A (uv) * bb equal T1-T3. Since the three triplet CSFs have the same energy, GAMESS works with the simpler form C2. Singlet and triplet computations using CSFs are done in separate runs, because when spin-orbit coupling is not considered, the Hamiltonian is block diagonal in a CSF basis. Technical information about the CSFs is that they use Yamanouchi- Kotani spin couplings, and matrix elements are obtained using a GUGA, or graphical unitary group approach. Determinant and CSF are both primarily used for MCSCF wavefunctions, but can be used in CI (see CITYP in $CONTRL). Other comparisons between the determinant and CSF implementations, as they exist in GAMESS today, are determinants CSFs parallel execution yes yes direct CI yes no use Abelian group symmetry yes yes state average mixed spins yes no first order density yes yes state averaged densities yes yes analytic nuclear hessian yes no can form CI Lagrangian no yes In nearly every circumstance the determinant CI will run faster than GUGA, so it is the default. Here are timings for N electrons in N orbitals, no symmetry used: N in N ALDET GENCI --- GUGA --- 8 1 1 1 0 10 8 38 19 33 12 228 3122 534 2209 14 7985 -- 15377 130855 The reason there are two numbers under GUGA is that the first is for writing the loops (Hamiltonian data) to disk, and the second is for the actual diagonalization. Note that the GUGA Hamiltonian time alone is greater than the entire ALDET computation! The ALDET program does not store anything on disk, and so runs at CPU/wall ratios of 1. The quality of the initial guess of the CI eigenvector in the various determinant codes is much better than in the GUGA CSF code, so the chances of converging to an incorrect excited state root is much less. Finally, determinant input is generally easier. No wonder determinants are now the default configurations! Two of the determinant CI programs, namely ALDET or ORMAS (but not GENCI) have been changed to use replicated memory parallelism, with modest scaling. The GUGA program is parallel in its solving, but not its Hamiltonian generation, and as already noted, its H formation takes more time than the entire determinant CI (translation: use CISTEP=ALDET, not CISTEP=GUGA). The GENCI program will run as a serial bottleneck (no speedup) in parallel runs. The next two sections describe in detail the input for specification of the configurations, either determinants or CSFs. determinant CI Three determinant CI codes are provided for MCSCF, one for full CI spaces (ALDET), another named the Occupation Restricted Multiple Active Spaces (ORMAS), and finally there is a program for arbitrary (selected) determinant lists (GENCI). For straight CI, but not MCSCF, there is a fourth program, the full second order CI (CITYP=FSOCI), whose purpose is MR-CISD. ALDET is a full CI within the chosen active space. It is possible to go up to 16 electrons in 16 orbitals, if your computer has a lot of memory. ALDET is the only CISTEP for which analytic nuclear hessian is possible, and it is also the most scalable CI code (using replicated memory to store CI vectors). A sample input for ALDET is $DET STSYM=B1 NSTATE=3 NCORE=xx NELS=8 NACT=6 $END Keywords in this group actually relate to all determinant programs, and are described below. The $DET input group is basic to all determinant CI codes. Keywords GROUP and STSYM specify the desired spatial symmetry of the determinants. Most runs need give only the orbital and electron counts: NCORE, NACT, and NELS. The number of electrons is 2*NCORE+NELS, and will be checked against the charge implied by ICHARG. The MULT given in $CONTRL is used to determine the desired Sz value, by extracting S from MULT=2S+1, then by default Sz=S. If you wish to include lower spin multiplicities, which will increase the CPU time of the run, but will let you know what the energies of such states are, just input a smaller value for SZ. The states whose orbitals will be MCSCF optimized will be those having the requested MULT value, unless you choose otherwise with the PURES flag. The remaining parameters in the $DET group give extra control over the diagonalization process. Most are not given in normal circumstances, except NSTATE, which you may need to adjust to produce enough roots of the desired MULT value. The only important keyword which has not been discussed is the WSTATE array, giving the weights for each state in forming the first and second order density matrix elements, which drive the orbital update methods. Note that analytic gradients are available only when the WSTATE array is a unit vector, corresponding to a pure state, such as WSTATE(1)=0,1,0 which permits gradients of the first excited state to be computed. When used for state averaged MCSCF, WSTATE is normalized to a unit sum, thus input of WSTATE(1)=1,1,1 really means a weight of 0.33333... for the each of the states being averaged. ORMAS (Occupation Restricted Multiple Active Space) is a program designed to limit the size of the full CI problem, and may be useful when the number of active orbitals is 10 or higher. By dividing your total active space into multiple subspaces, and specifying a range of electrons to occupy each subspace, most of the full CI's effect can be included. ORMAS generates a full CI within each orbital subspace, taking the product of each small full CI to generate the determinant list. Here are some ideas on how to use ORMAS, which is a very flexible CI program: a) single reference, arbitrary excition level CI-X, from a closed shell reference: $det ncore=y nact=z nels=10 (y+z=entire basis) $ormas nspace=2 mstart(1)=y+1,y+6 mine(1)=10-x,0 maxe(1)=10,x This excites the 5 doubly occupied orbitals, to the desired excitation level of X. An open shell example of CI-SD from ...22111 might be $contrl mult=4 $det ncore=y nact=z nels=7 (y+z=entire basis) $ormas nspace=3 mstart(1)=y+1,y+3,y+6 mine(1)=2,1,0 maxe(1)=4,5,2 No more than 2e- are allowed to be promoted from the doubly occupied or singly occupied spaces, and no more than 2 are allowed to enter the singly occupied or empty spaces. b) simple product of active spaces For example, consider furan, with two active subspaces. Keeping the 5 true core and the 4 CH bonds in the core space, the sigma subspace might contains 5 ring sigma, one oxygen lone pair, and 5 ring sigma antibonds, with a total of 12 e-. The pi active space contains 5 pi orbitals and 6 e-: $det ncore=9 nact=16 nels=18 $ormas nspace=2 mstart(1)=10,21 mine(1)=12,6 maxe(1)=12,6 Having the minimum and maximum electron counts the same is what makes this the simple product of two separate active spaces. In other words, this is similar to the QCAS procedure of Nakano and Hirao, but ORMAS limits only the total electron counts, not separately the numbers of alpha and beta e-, in other words all spin couplings are used. c) flexible occupancy between active subspaces Imagine that you are interested in excited states of formaldehyde, some of which will have Rydberg character, dominated by single excitations into diffuse orbitals. H2CO's valence states arise from 3 orbitals, the CO pi and pi* and one oxygen lone pair. Placing the O sp lone pair and three sigma bonds into the filled space, and centering diffuse s,p,d shells on the carbon: $det ncore=6 nact=12 nels=4 $ormas nspace=2 mstart(1)=7,10 mine(1)=3,0 maxe(1)=4,1 This is a 4e-, 3 orbital n,pi,pi* space to describe valence states, and excites one electron into the 9 diffuse orbitals to describe Rydberg states. It is many fewer determinants than a 4e- in 12 orbital FCI. d) RAS-like CI The previous example is reminiscent of Roos' RAS-SCF. In fact ORMAS can do RAS-SCF, which is three spaces: the lowest space is allowed to excite only a few electrons, a middle space that is the rest, and a top space into which only a few electrons can be excited. Suppose there are 10 e-, 10 orbitals, that the bottom and top spaces involve 3 orbitals, and that a "few" means specifically 2 e-: $det ncore=20 nact=10 nels=10 $end $ormas nspace=3 mstart(1)=21,24,28 mine(1)=4,2,0 maxe(1)=6,6,2 However, ORMAS can use more than 3 orbital subspaces. e) first or second order CI. Consider C2H4, with a 4 orbital active space of CC sigma, pi, pi*, and sigma*. In order to correlate the four valence CH orbitals by double excitations, an MCSCF based on $DET, followed by SOCI based on $CIDET and $ORMAS, is: $contrl scftyp=mcscf cityp=ormas $mcscf cistep=aldet $det ncore=6 nact=4 nels=4 $cidet ncore=2 nact=y nels=12 (y=rest of basis) $ormas nspace=3 mstart(1)=3,7,11 mine(1)=6,2,0 maxe(1)=8,6,2 which permits singles and doubles out of the CH and CC spaces, into the CC and external spaces. ORMAS is a full CI (or several full CI's) within each orbital subspace. However, ORMAS does not generate all excitation levels between spaces (just those implied by the minimum and maximum electron counts you give). This means ORMAS MCSCF runs must optimize active-active rotations between the subspaces, and therefore you should expect better convergence from FULLNR than SOSCF. ORMAS is sure to require orbital reordering. For the furan example just mentioned, there is no reason to expect that the RHF occupied orbitals will not have the filled sigma and pi orbitals intermingled. You must use the NORDER and IORDER keywords in $GUESS to carefully partition starting orbitals into sigma and pi subspaces. The selected (general) determinant list is used if CISTEP=GENCI, and the list is controlled by two input groups. The first is $GEN, which is identical to $DET except for inclusion of an additional keyword GLIST=INPUT. This reads the determinants (as space products) from an additional input group $GCILST. Completely arbitrary choices for the space products may be made, but peculiar lists may lead to poor MCSCF convergence. The FOCAS converger should not be used, as it assumes full CI spaces. If you are doing straight CI calculations, the required input for each determinant CITYP is: ALDET needs $CIDET ORMAS needs $CIDET and $ORMAS GENCI needs $CIDET and $CIGEN and probably $GCILST FSOCI needs $CIDET and $SODET In other words, $CIDET replaces $DET, and $CIGEN replaces $GEN, but the keywords in the groups mean the same thing. The reason for different names is to allow CI calculations to follow MCSCF in the same run, without clashing input group names. CSF CI The GUGA-based CSF package was originally a set of different programs, so its input is spread over several input groups. The CSFs are specified by a $CIDRT group in the case of CITYP=GUGA, and by a $DRT group for MCSCF wavefunctions. Thus it is possible to perform an MCSCF defined by a $DRT input (or perhaps using $DET during the MCSCF), and follow this with a second order CI defined by a $CIDRT group, in the same run. The remaining input groups used by the GUGA CSFs are $CISORT, $GUGEM, $GUGDIA, and $GUGDM2 for MCSCF runs, with the latter two being the most important, and in the case of CI computations, $GUGDM and possibly $LAGRAN groups are relevant. Perhaps the most interesting variables outside the $DRT/$CIDRT group are NSTATE in $GUGDIA to include excited states in the CI computation, IROOT in $GUGDM to select the CI state for properties, and WSTATE in $GUGDM2 to control which state's orbitals are optimized, and possible state-averaging. The $DRT and $CIDRT groups are almost the same, with the only difference being orbitals restricted to double occupancy are called MCC in $DRT, and FZC in $CIDET. Therefore the rest of this section refers only to "$DRT". The CSFs are specified by giving a reference CSF, together with a maximum degree of electron excitation from that single CSF. The MOs in the reference CSF are filled in the order MCC or FZC first, followed by DOC, AOS, BOS, ALP, VAL, and EXT (the Aufbau principle). AOS, BOS, and ALP are singly occupied MOs. ALP means a high spin alpha coupling, while AOS/BOS are an alpha/beta coupling to an open shell singlet. This requires the value NAOS=NBOS, and their MOs alternate. An example is NFZC=1 NDOC=2 NAOS=2 NBOS=2 NALP=1 NVAL=3 which gives the reference CSF FZC,DOC,DOC,AOS,BOS,AOS,BOS,ALP,VAL,VAL,VAL This is a doublet state with five unpaired electrons. VAL orbitals are unoccupied only in the reference CSF, they will become occupied as the other CSFs are generated. This is done by giving an excitation level, either explicitly by the IEXCIT variable, or implicitly by the FORS, FOCI, or SOCI flags. One of these four keywords must be chosen, and during MCSCF runs, this is usually FORS. Consider another simpler example, for an MCSCF run, NMCC=3 NDOC=3 NVAL=2 which gives the reference CSF MCC,MCC,MCC,DOC,DOC,DOC,VAL,VAL having six electrons in five active orbitals. MCSCF calculations are usually of the Full Optimized Reaction Space (FORS) type. Some workers refer to FORS as CASSCF, complete active space SCF. These are the same, but the keyword is spelled FORS in GAMESS. In the present instance, choosing FORS=.TRUE. gives an excitation level of 4, as the 6 valence electrons have only 4 holes available for excitation. MCSCF runs typically have only a small number of VAL orbitals. It is common to summarize this example as "six electrons in five orbitals". The next example is a first or second order multi- reference CI wavefunction, where NFZC=3 NDOC=3 NVAL=2 NEXT=-1 leads to the reference CSF FZC,FZC,FZC,DOC,DOC,DOC,VAL,VAL,EXT,EXT,... FOCI or SOCI is chosen by selecting the appropriate flag, the correct excitation level is automatically generated. Note that the -1 for NEXT causes all remaining MOs to be included in the external orbital space. One way of viewing FOCI and SOCI wavefunctions is as all singles, or all singles and doubles, from the entire MCSCF wavefunction as a reference. An equivalent way of saying this is that all CSFs with N electrons (in this case N=6) distributed in the valence orbitals in all ways (that is the FORS MCSCF wavefunction) make up the reference wavefunction. To this, FOCI adds all CSFs with N-1 electrons in active and 1 electron in external orbitals. SOCI adds all CSFs with N-2 electrons in active orbitals and 2 in external orbitals. SOCI is often prohibitively large, but is also a very accurate wavefunction. SOCI can also be performed with determinants, as CITYP=FSOCI, or CITYP=ORMAS. The latter may be the most efficient way to generate SOCI energies. For larger molecules, where SOCI is impractical, the most effective way to recover dynamic correlation energy is the multireference perturbation method. Sometimes people use the CI package for ordinary single reference CI calculations, such as NFZC=3 NDOC=5 NVAL=34 which means the reference RHF wavefunction is FZC FZC FZC DOC DOC DOC VAL VAL ... VAL and in this case NVAL is a large number conveying the total number of -virtual- orbitals into which electrons are excited. The excitation level would be given as IEXCIT=2, perhaps, to perform a SD-CI. All excitations smaller than the value of IEXCIT are automatically included in the CI. Note that NVAL's spelling was chosen to make the most sense for MCSCF calculations, and so it is a bit of a misnomer here. Before going on, there is a quirk related to single reference CI that should be mentioned. Whenever the single reference contains unpaired electrons, such as NFZC=3 NDOC=4 NALP=2 NVAL=33 some "extra" CSFs will be generated. The reference here can be abbreviated 2222 11 000 000 000 000 000 000 000 000 000 000 000 Supposing IEXCIT=2, the following CSF 2200 22 000 011 000 000 000 000 000 000 000 000 000 will be generated and used in the CI. Most people would prefer to think of this as a quadruple excitation from the reference, but acting solely on the reasoning that no more than two electrons went into previously vacant NVAL orbitals, the GUGA CSF package decides it is a double. So, an open shell SD-CI calculation with GAMESS will not give the same result as other programs, although the result for any such calculation with these "extras" is correctly computed. Note that if you also select the INTACT option, the extra space products are eliminated, but that some of the spin couplings for the truly IEXCIT'd space products are also eliminated. Note that this kind of problem does not arise if you use ORMAS! As was discussed above, the CSFs are automatically spin-symmetry adapted, with S implicit in the reference CSF. The spin quantum number you appear to be requesting in $DRT (basically, S = NALP/2) will be checked against the value of MULT in $CONTRL. The total number of electrons, 2*NMCC(or NFZC) + 2*NDOC + NAOS + NBOS + NALP will be checked against the input given for ICHARG. The CSF package is also able to exploit spatial symmetry, which like the spin and charge, is implicitly determined by the choice of the reference CSF. The keyword GROUP in $DRT governs the use of spatial symmetry. The CSF program works with Abelian point groups, which are D2h and any of its subgroups. However, $DRT allows the input of some (but not all) higher point groups. For non- Abelian groups, the program automatically assigns the orbitals to an irrep in the highest possible Abelian subgroup. For the other non-Abelian groups, you must at present select GROUP=C1. Note that when you are computing a Hessian matrix, many of the displaced geometries are asymmetric, hence you must choose C1 in $DRT (however, be sure to use the highest symmetry possible in $DATA!). The symmetry of the reference CSF given in your $DRT is one way to determine the symmetry of the CSFs which are generated. As an example, consider a molecule with Cs symmetry, and these two reference CSFs ...MCC...DOC DOC VAL VAL ...MCC...DOC AOS BOS VAL Suppose that the 2nd and 3rd active MOs have symmetries a' and a". Both of these generate singlet wavefunctions, with 4 electrons in 4 active orbitals, but the former constructs 1-A' CSFs, while the latter generates 1-A" CSFs. However, if the 2nd and 3rd orbitals have the same symmetry type, an identical list of CSFs is generated. The alternative is to enter the spatial symmetry with the STSYM keyword. In cases with high point group symmetry, it may be possible to generate correct state degeneracies only by using no symmetry (GROUP=C1) when generating CSFs. As an example, consider the 2-pi ground state of NO. If you use GROUP=C4V, which will be mapped into its highest Abelian subgroup C2v, the two components of the pi state will be seen as belonging to different irreps, B1 and B2. The only way to ensure that both sets of CSFs are generated is to enforce no symmetry at all, so that CSFs for both components of the pi level are generated. This permits state averaging (WSTATE(1)=0.5,0.5) to preserve cylindrical symmetry. It is however perfectly feasible to use C4v or D4h symmetry in $DRT when treating sigma states. The use of spatial symmetry decreases the number of CSFs, and thus the size of the Hamiltonian that must be computed. In molecules with high symmetry, this may lead to faster run times with the GUGA CSF code, compared to the determinant code. starting orbitals The first step is to partition the orbital space into core, active, and external sets, in a manner which is sensible for your chemical problem. This is a bit of an art, and the user is referred to the references quoted at the end of this section. Having decided what MCSCF to perform, you now must consider the more pedantic problem of what orbitals to begin the MCSCF calculation with. 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 HUCKEL! At the start of a MCSCF problem, use orbitals from some appropriate converged SCF run. A realistic example of an MCSCF calculation is GAMESS examples 8 and 9. Once you get an MCSCF to converge, you can and should use these MCSCF MOs at other nearby geometries (MOREAD will apply an appropriate Schmidt orthogonalization). 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 reordering is needed in 50% or more cases. The occupied and especially the virtual canonical SCF MOs are often spread out over regions of the molecule other than "where the action is". Orbitals which remedy this can generated by two additional options at almost no CPU cost. The best way to improve upon the SCF canonical virtual orbitals as starting MOs is to generate valence virtual orbitals (VVOs), after any RHF, ROHF, or GVB calculation. These are constructed by projection of internally stored atomic core and valence orbitals onto the SCF external orbitals, so by construction, the resulting VVOs are valence in character. See VVOS in $SCF. An alternative choice, usually not as good as VVOs, are the modified virtual orbitals. 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. Typically, we use MVOQ=+6. Generating MVOs does not change any of the occupied SCF orbitals of the original neutral, but gives more valence-like LUMOs. Another way to improve SCF starting orbitals is by a partial localization of the occupied orbitals. Typically MCSCF active orbitals are concentrated in the part of the molecule where bonds are breaking, etc. Canonical SCF MOs are normally more spread out. By choosing LOCAL=BOYS along with SYMLOC=.TRUE. in $LOCAL, you can get orbitals which are localized, but still retain orbital symmetry to help speed the MCSCF along. In groups with an inversion center, a SYMLOC Boys localization does not change the orbitals, but you can instead use LOCAL=POP. LOCAL=RUEDNBRG may also be used, but requires more machine resources, and the other two localizations are normally good enough for starting orbital purposes. Localization tends to order the orbitals fairly randomly, so be prepared to inspect them, and then reorder them appropriately. In case VVOS are generated in the same run as the localization, the localization is also applied within the valence virtual space. The effect is to localize occupied orbitals into lone pairs and bond pairs, and in the VVOs space, to find localized antibond pairs. When you choose the bonds and antibonds of chemical interest as the starting orbitals for the active space, convergence to the desired MCSCF wavefunction should follow. If you take the time to design your active space sensibly, select appropriate starting orbitals from the occupied and VVO unoccupied spaces, possibly by localizing these two subspaces, and carefully inspect your converged results, you will be able to carry out MCSCF computations correctly. Convergence of MCSCF is by no means guaranteed. Poor convergence can invariably be traced back to either a poor initial selection of orbitals, or poor design of the active space. The best 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". Few people are able to see the orbital shapes in the LCAO matrix in a log file, and so need a visualization program. In particular, you should download a copy of MacMolPlt from http://code.google.com/p/wxmacmolplt This runs on all popular desktop operating systems (Apple, Linux, and Windows), making it easy to see your final MCSCF orbital shapes. 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. It is particularly useful to check the oriented localized MCSCF orbitals (see the discussion of this in the section on localized orbitals in this section for more information). MCSCF is by no means the sort of "black box" that RHF is these days, so please look very carefully at your final results. miscellaneous hints It is very helpful to execute a EXETYP=CHECK run before doing any MCSCF or CI run. The CHECK run will tell you the total number of configurations and check the charge and multiplicity and electronic state symmetry, based on your input. The CHECK run also lets the program feel out the memory that will be required to actually do the run. Thus the CHECK run can potentially prevent costly mistakes, or tell you when a calculation is prohibitively large. A very common MCSCF wavefunction has 2 electrons in 2 active MOs. This is the simplest possible wavefunction describing a singlet diradical. While this function can be obtained in an MCSCF run (using NACT=2 NELS=2 or 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, or TCSCF. The two configurations of this GVB are equivalent to the three configurations used in this MCSCF, as orbital optimization in natural form (configurations 20 and 02) causes the coefficient of the 11 configuration to vanish. If you are using a large active space (say, 12 or more orbitals), the main bottleneck in the MCSCF calculation is the formation and diagonalization of the Hamiltonian, not the integral transformation and orbital updates. Of course, since determinants are much faster than CSFs, and do not use large disk files, you should use determinants for large active spaces. In this case, you would be wise to switch to FULLNR, which will minimize the total number of iterations, and thus the number of CI calculations. Note that by selecting ITERMX=5 in $DET or $GEN, you can avoid fully converging the CI during each MCSCF iteration, saving a bit of time. Since each iteration's CI calculation starts with the previous iteration's result, the CI vectors will become fully converged during the MCSCF cycles. The total run time may decrease, although a few additional MCSCF iterations may be required. For small active spaces, where the CI step takes trivial time, you should use a bigger ITERMX to ensure fully converged CI states are generated every iteration. If you choose to use ORMAS, a general determinant CI, or if you select an CSF excitation level IEXCIT smaller than that needed to generate the FORS space, you must use the SOSCF, JACOBI, or FULLNR method as these can optimize active-active rotations. Be sure to set FORS=.FALSE. in $MCSCF when for non-full CI cases, or else very poor convergence will result. Actually, the convergence for incomplete active spaces is likely to be poorer than for full active spaces, anyway. A good way to check the active space is to localize the orbitals, to see if they resemble the atomic orbitals which you imagined formed the bonds, antibonds, and lone pairs in the active space. The ORIENT keyword in $LOCAL will print a density matrix analysis, showing active electron bonding and antibonding patterns (see reference 18/19 below). - - - - - The MCSCF technology in GAMESS is the result of some considerable programming effort: The FOCAS, serial FULLNR, and QUAD convergers were adapted from Michel Dupuis' HONDO program. The SOSCF converger was written by Galina Chaban, the parallel FULLNR converger is due to Graham Fletcher, and the JACOBI converger is due to Joe Ivanic. The GUGA CI programs were written by Bernie Brooks and others, while all determinant CI codes (ALDET, GENCI, ORMAS, and FSOCI) stem from Joe Ivanic. Analytic nuclear Hessians were programmed by Tim Dudley. The CSF-based multireference pertubation program was written by Haruyuki Nakano, with a determinant implementation provided by Joe Ivanic. Shiro Koseki and Dmitri Fedorov are responsible for the spin- orbit coupling and transition moment codes. The expertise of Klaus Ruedenberg in MCSCF wavefunctions has been the inspiration for many of these developments! MCSCF references There are several review articles about MCSCF listed below. Of these, the first two are a nice overview of the subject, the final 3 are more technical. 1. "The Construction and Interpretation of MCSCF wavefunctions" M.W.Schmidt and M.S.Gordon, Ann.Rev.Phys.Chem. 49,233-266(1998) 2a. "The Multiconfiguration 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. 2b. "The Multiconfiguration SCF Method" B.O.Roos, in "Lecture Notes in Quantum Chemistry", edited by B.O.Roos, Lecture Notes in Chemistry v58, Springer-Verlag, Berlin, 1994, pp 177-254. 3. "Optimization and Characterization of a MCSCF State" J.Olsen, D.L.Yeager, P.Jorgensen Adv.Chem.Phys. 54, 1-176(1983). 4. "Matrix Formulated Direct MCSCF and Multiconfiguration Reference CI Methods" H.-J.Werner, Adv.Chem.Phys. 69, 1-62(1987). 5. "The MCSCF Method" R.Shepard, Adv.Chem.Phys. 69, 63-200(1987). There is an entire section on the choice of active spaces in Reference 1. As this is a matter of great importance, here are two alternate presentations of the design of active spaces: 6. "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. 7. "Are Atoms Intrinsic to Molecular Electronic Wavefunctions?" K.Ruedenberg, M.W.Schmidt, M.M.Gilbert, S.T.Elbert Chem.Phys. 71, 41-49, 51-64, 65-78 (1982). Two papers germane to the FOCAS implementation are 8. "An Efficient first-order CASSCF method based on the renormalized Fock-operator technique." U.Meier, V.Staemmler Theor.Chim.Acta 76, 95-111(1989) 9. "Modern tools for including electron correlation in electronic structure studies" M.Dupuis, S.Chen, A.Marquez, in "Relativistic and Electron Correlation Effects in Molecules and Solids", edited by G.L.Malli, Plenum, NY 1994 The paper germane to the the SOSCF converger is 10. "Approximate second order method for orbital optimization of SCF and MCSCF wavefunctions" G.Chaban, M.W.Schmidt, M.S.Gordon Theor.Chem.Acc. 97: 88-95(1997) Two papers germane to the FULLNR converger, and two discussing implementation details are 11. "General second order MCSCF theory: A Density Matrix Directed Algorithm" B.H.Lengsfield, III, J.Chem.Phys. 73,382-390(1980). 12. "The use of the Augmented Matrix in MCSCF Theory" D.R.Yarkony, Chem.Phys.Lett. 77,634-635(1981). 13. M.Dupuis, P.Mougenot, J.D.Watts, in "Modern Techniques in Theoretical Chemistry", E.Clementi, editor, ESCOM, Leiden, 1989, chapter 7. 14. "A parallel multi-configuration self-consistent field algorithm" G.D.Fletcher, Mol.Phys. 105, 2971-2976(2007) The paper describing the JACOBI converger is 15. "A MCSCF method for ground and excited states based on full optimizatons of successive Jacobi rotations" J.Ivanic, K.Ruedenberg J.Comput.Chem. 24, 1250-1262(2003) For determinant CI codes, see 16. "Identification of deadwood in configuration spaces through general direct configuration interaction" J.Ivanic, K.Ruedenberg Theoret.Chem.Acc. 106, 339-351(2001) 17. "Direct configuration interaction and multi- configurational self-consistent-field method for multiple active spaces with variable occupancies. Part I.Method Part II.Applications" J.Ivanic J.Chem.Phys. 119, 9364-9376 and 9377-9385(2003) For CSFs, see 18. "GUGA approach to the electron correlation problem" B.R.Brooks, H.F.Schaefer J.Chem.Phys. 70, 5092-5106(1979) Orientation of localized MCSCF active orbitals for bonding analysis: 19. J.Ivanic, G.M.Atchity, K.Ruedenberg Theoret.Chem.Acc. 120, 281-294(2008) 20. J.Ivanic, K.Ruedenberg Theoret.Chem.Acc. 120, 295-305(2008)