7.8. Choice of Wavefunction and Integral Handling

7.8.1. Choice of Wavefunction Type

The basic variable that controls the type of wavefunction to be computed is the variable HFTyp in the %scf block. If nothing is specified for HFTyp, the program will check the multiplicity given in the input: for closed-shell molecules with multiplicity 1, RHF/RKS is assumed; for open shell molecules with multiplicity larger than 1, UHF/UKS is invoked. RHF will lead to a spin restricted closed-shell type computation [839]. For DFT calculations, RKS, UKS and ROKS can be used as synonyms for RHF, UHF and ROHF. The restricted open-shell DFT method (ROKS) is only operative for high-spin states that have \(n\) unpaired electrons and \(S=n/2\). UKS wavefunctions will not be spin-purified.

%scf
  HFTyp  RHF     # closed-shell (RKS for DFT)
         UHF     # unrestricted open-shell (UKS for DFT)
         ROHF    # restricted open-shell (ROKS for DFT)
         CASSCF  # complete active space SCF
end

In certain cases you may want to run open-shell molecules with RHF/RKS to get a “half-electron” type wavefunction [205]. The total energy is not corrected! Sometimes these half-electron computations lead to acceptable convergence, and the resulting orbitals may be used as input for ROHF, UHF or MRCI calculations. Especially for transition metal complexes the orbitals are quite different from ROHF or UHF orbitals, so that it is not recommended to over-interpret the wavefunctions from such calculations. The calculation is set up in the following way:

%method AllowRHF true end
# or simply: ! AllowRHF

7.8.2. ROHF Options

For ROHF calculations[102, 112, 125, 141, 243, 452, 453, 574, 597] the program will try to figure out what type of open-shell situation is present on the basis of the initial guess orbitals and their energies. Most “simple” cases are well recognized, but sometimes a little help from the user is needed.

The simplest ROHF case is the HIGHSPIN case, where all unpaired electrons in the open-shell are coupled parallel to each other, resulting in the highest multiplicity possible. The user can request this case as follow:

%scf
  HFTyp       ROHF
  ROHF_case   HIGHSPIN
  ROHF_NEl[1] 4          # number of electrons in the open-shell
end

The ROHF code also has a very powerful feature that goes back to insights of Mike Zerner [817, 906]. It can average over either all states of a given configuration (CAHF) or all states of a given spin for a given configuration (SAHF). Especially the SAHF feature gives you easy access to most degenerate high symmetry situations and the orbitals resulting from such calculations will be very convenient as input for CI calculations.

%scf
  HFTyp         ROHF
  ROHF_case     CAHF    # configuration averaged HF
                SAHF    # spin averaged HF
  ROHF_NumOp    3       # number of operators (3, 2 or 1)
  ROHF_NOrb[1]  2,1     # number of orbitals in each open-shell
  ROHF_NEl[1]   1,1     # number of electrons in each open-shell
end

The hypothetical example below could represent an excited state of an octahedral d\(^{3}\) transition metal complex. In this case there are five open-shell orbitals. The first three open-shell orbitals contain two electrons and the last two one electron. The input for a SAHF calculation is identical, just replace CAHF with SAHF.

%scf
  HFTyp         ROHF
  ROHF_case     CAHF   # configuration averaged HF
  ROHF_NumOp    3      # 3 operators in this case: closed, open1, open2
  ROHF_NOrb[1]  3,2    # 3 orbitals in first open shell, 2 in the second
  ROHF_NEl[1]   2,1    # 2 electrons in first open shell, 1 in the second
end

Another feature of the ROHF code is the ability to converge the SCF to a given Configuration State Function (CSF-ROHF) [512]. In this way one can approach results from MCSCF calculations. This can be requested in two ways.

The user can give a specific coupling situation.

%scf
  HFTyp      ROHF
  ROHF_CASE  USER_CSF          # User defined CSF
  ROHF_REF   {1 1 -1 -1} end   # CSF to be converged to
end

Or the user can give how many orbitals per shell. Where each open-shell will couple with antiparallel spin with the previous one.

%scf
  HFTyp       ROHF
  ROHF_CASE   AF_CSF          # User defined CSF
  ROHF_AFORBS 2,2             # Coupling Situation
end

As an example, one can think of a Fe(III) dimer, where each center is locally high spin, but they couple antiferromagnetically to each other. In order to get the ROHF solution for this system, first one need a set of guess orbitals. The guess orbitals can be obtained either from the QROs of an UHF calculation, or a high spin ROHF calculation, or even a SAHF or CAHF. Independently of the method used, the orbitals need to be localized and ordered in a way that the 5 3d orbitals of each iron are grouped together in sequence. From this, one can run a CSF-ROHF calculation for the antiferromagnetic CSF as shown bellow:

%scf
  HFTyp      ROHF
  ROHF_CASE  USER_CSF                         # User defined CSF
  ROHF_REF   {1 1 1 1 1 -1 -1 -1 -1 -1} end   # 2 open shells coupling antiparallel
end

or

%scf
  HFTyp       ROHF
  ROHF_CASE   AF_CSF                          # User defined CSF
  ROHF_AFORBS 5,5                             # 2 open shells coupling antiparallel
end

The CSF-ROHF procedure can recognize doubly occupied and virtual orbitals in the definition of the CSF when the USER_CSF case is invoked. When detected, these orbitals will be rotated out of the open-shells defined in the ROHF method and the calculation will run normally:

%scf
  HFTyp      ROHF
  ROHF_CASE  USER_CSF                   # User defined CSF
  ROHF_REF   {1 1 2 1  -1 -1 0 -1} end  # the DOMO will be rotated to the closed-shell and
                                        # the VMO will be rotated to the virtual space. 
end

The user can also directly input the ROHF variables by means of the ROHFOPT Case User keyword. For example for the high spin case with three electrons in three orbitals gives two operators with vector coupling coefficients \(a=1\) and \(b=2\) (Zerner convention).

%scf
  HFTyp      ROHF
  ROHFOP  Case    User  # manual input of ROHF variables
    Nop     2           # number of operators
    Norb[1] 3           # number of open-shell orbitals
    Nel[1]  3           # number of open-shell electrons
    A[1,1]  1           # Coulomb vector in the open shell
    B[1,1]  2           # Exchange vector in the open shell
  end
end

One awkward feature of the ROHF theory is that the Fock operator is somewhat arbitrarily defined. Different choices lead to the same wavefunction, but have different convergence properties that may vary from system to system. ORCA thus lets the user choose the desired variant. Playing around with these choices may turn a divergent or slowly converging ROHF calculation into a successful calculation!

The ROHF_Restrict feature is another feature that may be useful. If you suspect that the ROHF calculation does not converge because an open-shell and a closed-shell orbital are flipping back and forth, you can try to avoid this behavior by choosing ROHF_Restrict true. Of course there is no guarantee that it will work, and no guarantee that the system stays in the desired state. However, it decreases the chances of large, uncontrolled steps.

%scf
  ROHF_Mode 0   # construct F according to Pulay (default)
            1   # construct F as in the Gamess program
            2   # construct F according to Kollmar

  ROHF_Restrict false # restrict orbital interchanges and off-diagonal elements
                      # (default=false)
                      
  # a complete list of ROHF variables
  ROHFOP 
    Case      User   # manual input of ROHF variables
    Nop       2      # number of operators
    Norb[1]   3      # number of open-shell orbitals
    Nel[1]    3      # number of open-shell electrons
    A[1,1]    1      # Coulomb vector in the open shell
    B[1,1]    2      # Exchange vector in the open shell
    Mode      2      # use the Kollmar operator
    Restrict  false  # do not restrict
  end
end

7.8.3. UHF Natural Orbitals

The program can produce the UHF natural orbitals (UNOs). With these, the open-shell wavefunction can be pictured conveniently. The syntax is simple:

%scf
  UHFNO true
end
# or simply: ! UNO

There are various printing options for UNOs described in the output section Population Analyses and Control of Output. The UNOs can also be plotted as described in the plots section Orbital and Density Plots. In general the program stores a file BaseName.uno, where BaseName is by default the name of you input file with .inp stripped off. Accordingly, the gbw file is named BaseName.gbw. The .uno file is a normal gbw file that contains the geometry, basis set and the UNO orbitals. It could be used, for example, to start a ROHF calculation.

7.8.4. Integral Handling (Conventional and Direct)

As the number of nonzero integrals grows very rapidly and reaches easily hundreds of millions even with medium sized basis sets in medium sized molecules, storage of all integrals is not generally feasible. This desperate situation prevented SCF calculations on larger molecules for quite some time, so that Almlöf [19, 20, 21] made the insightful suggestion to repeat the integral calculation, which was already the dominant step, in every SCF cycle to solve the storage problem. Naively, one would think that this raises the effort for the calculation to \(n_{\text{iter} }t_{\text{integrals} }\) (where \(n_{\text{iter} }\) is the number of iterations and \(t_{\text{integrals} }\) is the time needed to generate the nonzero integrals). However, this is not the case because only the change in the Fock matrix is required from one iteration to the next, but not the Fock matrix itself. As the calculations starts to converge, more and more integrals can be skipped. The integral calculation time will still dominate the calculation quite strongly, so that ways to reduce this burden are clearly called for. As integrals are calculated in batches[1] the cost of evaluating the given batch of shells \(p, q, r, s\) may be estimated as:

(7.46)\[\text{cost}\approx n_{p} n_{q} n_{r} n_{s} \left({ 2l_{p} +1} \right)\left({ 2l_{q} +1} \right)\left({ 2l_{r} +1} \right)\left({ 2l_{s} +1} \right)\]

Here, \(n_{p}\) is the number of primitives involved in shell \(p\), and \(l_{p}\) is the angular momentum for this shell. Large integrals are also good candidates for storage, because small changes in the density that multiply large integrals are likely to give a nonzero contribution to the changes in the Fock matrix.

ORCA thus features two possibilities for integral handling, which are controlled by the variable SCFMode. In the mode Conventional, all integrals above a given threshold are stored on disk (in a compressed format that saves much disk space). In the mode Direct, all two-electron integrals are recomputed in each iteration.

Two further variables are of importance: In the Conventional mode the program may write enormous amounts of data to disk. To ensure this stays within bounds, the program first performs a so-called “statistics run” that gives a pessimistic estimate of how large the integral files will be. Oftentimes, the program will overestimate the amount of disk space required by a factor of two or more. The maximum amount of disk space that is allowed for the integral files is given by MaxDisk (in Megabytes).

On the other hand, if the integral files in Conventional run are small enough to fit into the central memory, it is faster to do this since it avoids I/O bottlenecks. The maximum amount of memory allocated for integrals in this way is specified by MaxIntMem (in Megabytes). If the integral files are larger than MaxIntMem, no integrals will be read into memory.

%scf
  MaxIter   100    # Max. no. of SCF iterations
  SCFMode Direct   # default, other choice: Conventional
  Thresh  1e-8     # Threshold for neglecting integrals / Fock matrix contributions
                   # Depends on the chosen convergence tolerance (in Eh).
  TCut    1e-10    # Threshold for neglecting primitive batches. If the prefactor
                   # in the integral is smaller than TCut, the contribution of the
                   # primitive batch to the total batch is neglected.
  UseCheapInts false   # default: false
  DirectResetFreq 20   # default: 15
  MaxDisk  2500        # Max. amount of disk for 2 el. ints. (MB)
  MaxIntMem 400        # Max. amount of RAM for 2 el. ints. (MB)
end

The flag UseCheapInts has the following meaning: In a Direct SCF calculation, the oscillations in the total energy and density are initially quite large. High accuracy in the integrals is therefore not crucial. If UseCheapInts is switched on, the program loosens the threshold for the integrals and thus saves a lot of computational time. After having obtained a reasonable initial convergence, the thresholds are tightened to the target accuracy. One pitfall with this method is that the number of cycles required to reach convergence may be larger relative to a calculation with full integral accuracy throughout.[2] When restarting calculations that are close to convergence, it is recommended to switch UseCheapInts off. UseCheapInts has no meaning in a conventional SCF.

The value of DirectResetFreq sets the number of incremental Fock matrix builds after which the program should perform a full Fock matrix build in a Direct SCF calculation. To prevent numerical instabilities that arise from accumulated errors in the recursively build Fock matrix, the value should not be too large, since this will adversely affect the SCF convergence. If the value is too small, the program will update more frequently, but the calculation will take considerably longer, since a full Fock matrix build is more expensive than a recursive one.

The thresholds TCut and Thresh also deserve a closer explanation. Thresh is a threshold that determines when to neglect two-electron integrals. If a given integral is smaller than Thresh Eh, it will not be stored or used in Fock matrix construction. Additionally, contributions to the Fock matrix that are smaller than Thresh Eh will not be calculated in a Direct SCF.

Clearly, it would be wasteful to calculate an integral, then find out it is good for nothing and thus discard it. A useful feature would be an efficient way to estimate the size of the integral before it is even calculated, or even have an estimate that is a rigorous upper bound on the value of the integral. Häser and Ahlrichs [346] were the first to recognize that such an upper bound is actually rather easy to calculate. They showed that:

(7.47)\[\left|{ \left\langle { ij\left|{ kl} \right.} \right\rangle} \right|\leqslant \sqrt{ \left\langle { ij\left|{ ij} \right.} \right\rangle} \sqrt{ \left\langle { kl\left|{ kl} \right.} \right\rangle} \]

where:

(7.48)\[\left\langle { ij\left|{ kl} \right.} \right\rangle=\int{ \int{ \phi_{i} \left({ \vec{{r} }_{1} } \right)\phi_{j} \left({ \vec{{r} }_{1} } \right)r_{12}^{-1} \phi_{k} \left({ \vec{{r} }_{2} } \right)\phi_{l} \left({\vec{{r} }_{2} } \right)d\vec{{r} }_{1} d\vec{{r} }_{2} } } \]

Thus, in order to compute an upper bound for the integral only the right hand side of this equation must be known. This involves only two index quantities, namely the matrix of two center exchange integrals \(\left\langle { ij\left|{ ij} \right.} \right\rangle\). These integrals are easy and quick to calculate and they are all \(\geqslant\)0 so that there is no trouble with the square root. Thus, one has a powerful device to avoid computation of small integrals. In an actual calculation, the Schwartz prescreening is not used on the level of individual basis functions but on the level of shell batches because integrals are always calculated in batches. To realize this, the largest exchange integral of a given exchange integral block is looked for and its square root is stored in the so called pre-screening matrix K (that is stored on disk in ORCA). In a Direct SCF this matrix is not recalculated in every cycle, but simply read from disk whenever it is needed. The matrix of exchange integrals on the level of individual basis function is used in Conventional calculations to estimate the disk requirements (the “statistics” run).

Once it has been determined that a given integral batch survives it may be calculated as:

(7.49)\[\left\langle { ij\left|{ kl} \right.} \right\rangle=\sum\limits_p { d_{pi} \sum\limits_q { d_{qj} \sum\limits_r { d_{kr} \sum\limits_s { d_{sl} \left\langle { i_{p} j_{q} \left|{ k_{r} l_{s} } \right.} \right\rangle} } } } \]

where the sums \(p,q,r,s\) run over the primitive Gaussians in each basis function \(i,j,k,l\) and the \(d\)’s are the contraction coefficients. There are more powerful algorithms than this one and they are also used in ORCA. However, if many terms in the sum can be skipped and the total angular momentum is low, it is still worthwhile to compute contracted integrals in this straightforward way. In equation (7.49), each primitive integral batch \(\left\langle { i_{p} j_{q} \left|{ k_{r} l_{s} } \right.} \right\rangle\) contains a prefactor \(I_{Gaussians}\) that depends on the position of the four Gaussians and their orbital exponents. Since a contracted Gaussian usually has orbital exponents over a rather wide range, it is clear that many of these primitive integral batches will contribute negligibly to the final integral values. In order to reduce the overhead, the parameter TCut is introduced. If the common prefactor \(I_{pqrs}\) is smaller than TCut, the primitive integral batch is skipped. However, \(I_{pqrs}\) is not a rigorous upper bound to the true value of the primitive integral. Thus, one has to be more conservative with TCut than with Thresh. In practice it appears that choosing TCut=0.01*Thresh provides sufficient accuracy, but the user is encouraged to determine the influence of TCut if it is suspected that the accuracy reached in the integrals is not sufficient.

Hint

• If the direct SCF calculation is close to convergence but fails to finally converge, this maybe related to a numerical problem with the Fock matrix update procedure – the accumulated numerical noise from the update procedure prevents sharp convergence. In this case, set Thresh and TCut lower and/or let the calculation more frequently reset the Fock matrix (DirectResetFreq).

Note

• For a Direct calculation, there is no way to have Thresh larger than TolE. If the errors in the Fock matrix are larger than the requested convergence of the energy, the change in energy can never reach TolE. The program checks for that.

• The actual disk space used for all temporary files may easily be larger than MaxDisk. MaxDisk only pertains to the two-electron integral files. Other disk requirements are not currently checked by the program and appear to be uncritical.

---

[1]

A batch is a set of integrals that arises from all components of the shells involved in the integral. For example a \(\left\langle pp|pp \right\rangle\) batch gives rise to \(3 \times 3\times 3 \times 3 = 81\) integrals due to all possible combinations of \(p_{x}\), \(p_{y}\) and \(p_{z}\) functions in the four shells. Computations based on batches lead to great computational advantages because the 81 integrals involved in the \(\left\langle pp|pp \right\rangle\) batch share many common intermediate quantities.

[2]

This might be an undesirable feature of the current implementation.