7.7. SCF Convergence

SCF convergence is a pressing problem in any electronic structure package because the total execution times increases linearly with the number of iterations. Thus, it remains true that the best way to enhance the performance of an SCF program is to make it converge better. In some cases, especially for open-shell transition metal complexes, convergence may be very difficult. ORCA makes a dedicated effort to achieve reasonable SCF convergence for these cases without compromising efficiency.

Another issue is whether the solution found by ORCA is stable, i.e. a minimum on the surface of orbital rotations. Especially for open-shell singlets it can be hard to achieve a broken-symmetry solution. The SCF stability analysis (section SCF Stability Analysis) may be able to help in such situations. Please also note that if ! TRAH is used the solution must be a true local minimum though not necessarily a global.

7.7.1. Convergence Tolerances

Before discussing how to converge a SCF calculation it should be defined what is meant by “converged”. ORCA has a variety of options to control the target precision of the energy and the wavefunction that can be selected in the % scf block, or with a simple input line keyword that merges the criterion label with “SCF”, e.g. ! StrongSCF or ! VeryTightSCF:

%scf
  Convergence             # The default convergence is between medium and strong
               Sloppy     # very weak convergence
               Loose      # still weak convergence
               Medium     # intermediate accuracy
               Strong     # stronger
               Tight      # still stronger
               VeryTight  # even stronger
               Extreme    # close to numerical zero of the computer
                          # in double precision arithmetic
 end

Like other keys, Convergence is a compound key that assigns default values to a variety of other variables given in the box below. In table Threshold choices for compound convergence keys we present the chosen values for each compound key. If the corresponding simple inputs are given (StrongSCF, VeryTightSCF, …etc), then in addition the values of table Additional threshold choices set by the simple input keys (strongSCF, …etc.) are also set. The default convergence criteria are reasonable and should be sufficient for most purposes. For a cursory look at populations weaker convergence may be sufficient, whereas other cases may require stronger than default convergence. Note that Convergence does not only affect the target convergence tolerances but also the integral accuracy as discussed in the section about direct SCF and alike. This is very important: if the error in the integrals is larger than the convergence criterion, a direct SCF calculation cannot possibly converge.

The convergence criteria are always printed in the output. Given below is a list of the convergence criteria for ! TightSCF, which is often used for calculations on transition metal complexes.

%scf
  TolE        1e-8  # energy change between two cycles
  TolRMSP     5e-9  # RMS density change
  TolMaxP     1e-7  # maximum density change
  TolErr      5e-7  # DIIS error convergence
  TolG        1e-5  # orbital gradient convergence
  TolX        1e-5  # orbital rotation angle convergence
  ConvCheckMode  2  # = 0: check all convergence criteria
                    # = 1: stop if one of criterion is met, this is sloppy!
                    # = 2: check change in total energy and in one-electron energy
                    #      Converged if delta(Etot)<TolE and delta(E1)<1e3*TolE
  ConvForced        # = 0: convergence not mandatory for next calculation step
                    # = 1: break, if you did not meet the convergence criteria
end
Table 7.10 Threshold choices for compound convergence keys

Sloppy

TolE

3e-5

TolMAXP

1e-4

TolRMSP

1e-5

TolErr

1e-4

Thresh

1e-9

TCut

1e-10

DFTGrid.BFCut

1e-10

TolG

3e-4

TolX

3e-4

Z_Tol

5e-3

Loose

TolE

1e-5

TolMAXP

1e-3

TolRMSP

1e-4

TolErr

5e-4

Thresh

1e-9

TCut

1e-10

DFTGrid.BFCut

1e-10

TolG

1e-4

TolX

1e-4

Z_Tol

3e-3

Medium

SCFConvMode

_CONVCHECK_ENERGY

TolE

1e-6

TolMAXP

1e-5

TolRMSP

1e-6

TolErr

1e-5

Thresh

1e-10

TCut

1e-11

DFTGrid.BFCut

1e-10

TolG

5e-5

TolX

5e-5

Z_Tol

1e-3

Strong

SCFConvMode

_CONVCHECK_ENERGY

TolE

3e-7

TolMAXP

3e-6

TolRMSP

1e-7

TolErr

3e-6

Thresh

1e-10

TCut

3e-11

DFTGrid.BFCut

3e-11

TolG

2e-5

TolX

2e-5

Z_Tol

7e-4

Tight

SCFConvMode

_CONVCHECK_ENERGY

TolE

1e-8

TolMAXP

1e-7

TolRMSP

5e-9

TolErr

5e-7

Thresh

2.5e-11

TCut

2.5e-12

DFTGrid.BFCut

1e-11

TolG

1e-5

TolX

1e-5

Z_Tol

1e-4

VeryTight

SCFConvMode

_CONVCHECK_ENERGY

TolE

1e-9

TolMAXP

1e-8

TolRMSP

1e-9

TolErr

1e-8

Thresh

1e-12

TCut

1e-14

DFTGrid.BFCut

1e-12

TolG

2e-6

TolX

2e-6

Z_Tol

3e-5

Extreme

SCFConvMode

_CONVCHECK_ALL

TolE

1e-14

TolMAXP

1e-14

TolRMSP

1e-14

TolErr

1e-14

Thresh

3e-16

TCut

3e-16

TolG

1e-09

TolX

1e-09

Z_Tol

3e-06

DFTGrid.BFCut

3e-16

Table 7.11 Additional threshold choices set by the simple input keys (strongSCF, …etc.)

SloppySCF

DCAS.TolG

5.0e-3

DCAS.TolE

1.0e-6

DMDCI.STol

1.0e-4

DMRCI.ETol

1.0e-5

DMRCI.RTol

1.0e-5

DCIS.ETol

1.0e-5

DCIS.RTol

1.0e-5

LooseSCF

IN.DCAS.TolG

5.0e-3

DCAS.TolE

1.0e-6

DMDCI.STol

1.0e-4

DMRCI.ETol

1.0e-5

DMRCI.RTol

1.0e-5

DCIS.ETol

1.0e-5

DCIS.RTol

1.0e-5

NORMALSCF

DCAS.TolG

1.0e-3

DCAS.TolE

1.0e-7

DMDCI.STol

2.5e-5

DMRCI.ETol

1.0e-6

DMRCI.RTol

1.0e-6

DCIS.ETol

1.0e-6

DCIS.RTol

1.0e-6

STRONGSCF

DCAS.TolG

5.00e-4

DCAS.TolE

6.66e-8

DMDCI.STol

7.50e-6

DMRCI.ETol

6.66e-7

DMRCI.RTol

6.66e-7

DCIS.ETol

6.66e-7

DCIS.RTol

6.66e-7

TIGHTSCF

DCAS.TolG

2.5e-4

DCAS.TolE

2.5e-8

DMDCI.STol

1.0e-5

DMRCI.ETol

2.5e-7

DMRCI.RTol

2.5e-7

DCIS.ETol

2.5e-7

DCIS.RTol

2.5e-7

VERYTIGHTSCF

DCAS.TolG

1.0e-5

DCAS.TolE

1.0e-8

DMDCI.STol

1.0e-6

DMRCI.ETol

1.0e-7

DMRCI.RTol

1.0e-7

DCIS.ETol

1.0e-7

DCIS.RTol

1.0e-7

EXTREMESCF

DCAS.TolG

1.0e-9

DCAS.TolE

1.0e-12

DMDCI.STol

1.0e-9

DMDCI.TCutInt

0.0

DMRCI.ETol

1.0e-12

DMRCI.RTol

1.0e-12

DCIS.ETol

1.0e-12

DCIS.RTol

1.0e-12

If ConvCheckMode=0, all convergence criteria have to be satisfied for the program to accept the calculation as converged, which is a quite rigorous criterion. In this mode, the program also has mechanisms to decide that a calculation is converged even if one convergence criterion is not fulfilled but the others are overachieved. ConvCheckMode=1 means that one criterion is enough. This is quite dangerous, so ensure that none of the criteria are too weak, otherwise the result will be unreliable. The default ConvCheckMode=2 is a check of medium rigor — the program checks for the change in total energy and for the change in the one-electron energy. If the ratio of total energy and one-electron energy is constant, the self-consistent field does not fluctuate anymore and the calculation can be considered converged. If you have small eigenvalues of the overlap matrix, the density may not be converged to the number of significant figures requested by TolMaxP and TolRMSP.

ConvForced is a flag to prevent time consuming calculations on non-converged wave functions. It will default to ConvForced=1 for Post-HF methods, Excited States runs and Broken Symmetry calculations. You can overwrite this default behavior by setting ConvForced=0.

Irrespective of the ConvForced value that has been chosen, properties or numerical calculations (NumGrad, NumFreq) will not be performed on non-converged wavefunctions!

7.7.2. Dynamic and Static Damping

Damping is the oldest and simplest convergence aid. It was already invented by Douglas Hartree when he did his famous atomic calculations. Damping consists of mixing the old density with the new density as:

(7.45)\[P_{\text{new, damped} } =\left({ 1-\alpha } \right)P_{\text{new} } +\alpha P_{\text{old} } \]

where \(\alpha\) is the damping factor, which must have a value of less than 1. Thus the permissible range (not checked by the program) is 0 …0.999999. For \(\alpha\) values larger than 1, the calculation cannot proceed since no new density is admixed. Damping is important in the early stages of a calculation where \(P_{\text{old} }\) and \(P_{\text{new} }\) are very different from each other and the energy is strongly fluctuating. Many schemes have been suggested that vary the damping factor dynamically to give strong damping in the beginning and no damping in the end of an SCF. The scheme implemented in ORCA is that by Hehenberger and Zerner [911] and is invoked with CNVZerner=true. Static damping is invoked with CNVDamp=true. These convergers are mutually exclusive. They can be used in the beginning of a calculation when it is not within the convergence radius of DIIS or SOSCF. Damping works reasonably well, but most other convergers in ORCA are more powerful.

If damping used in conjunction with DIIS or SOSCF, the value of DampErr is important: once the DIIS error falls below DampErr, the damping is turned off. In case SOSCF is used, DampErr refers to the orbital gradient value at which the damping is turned off. The default value is 0.1 Eh. In difficult cases, however, it is a good idea to choose DampErr much smaller, e.g. 0.001. This is — to some extent — chosen automatically together with the keyword ! SlowConv.

%scf
  # control of the Damping procedure
  CNVDamp   true    # default: true
  CNVZerner false   # default: false
  DampFac   0.98    # default: 0.7
  DampErr   0.05    # default: 0.1
  DampMin   0.1     # default: 0.0
  DampMax   0.99    # default: 0.98
  # more convenient:
  Damp fac 0.98 ErrOff 0.05 Min 0.1 Max 0.99 end
end

7.7.3. Level Shifting

Level shifting is a frequently used technique. The basic idea is to shift the energies of the virtual orbitals such that after diagonalization the occupied and virtual orbitals mix less strongly and the calculation converges more smoothly towards the desired state. Also, level shifting should prevent flipping of electronic states in near-degenerate cases. In a special context it has been shown by Saunders and Hillier [335, 751] to be equivalent to damping.

Similar to DampErr described in the previous section, ShiftErr refers to the DIIS error at which the level shifting is turned off.

%scf
  # control of the level shift procedure
  CNVShift   true  # default: true
  LShift     0.1   # default: 0.25, energy unit is Eh.
  ShiftErr   0.1   # default: 0.0
  # more convenient:
  Shift Shift 0.1 ErrOff 0.1 end
end

7.7.4. Direct Inversion in Iterative Subspace (DIIS)

The direct inversion in iterative subspace (DIIS) is a technique that was invented by Pulay [702, 703]. It has become the de facto standard in most modern electronic structure programs, because DIIS is robust, efficient and easy to implement. Basically DIIS uses a criterion to judge how far a given trial density is from self-consistency. The commutator of the Fock and density matrices [F,P] is a convenient measure for this error. With this information, an extrapolated Fock matrix from the present and previous Fock matrices is constructed, which should be much closer to self-consistency. In practice this is usually true, and better than linear convergence has been observed with DIIS. In some rare (open-shell) cases however, DIIS convergence is slow or absent after some initial progress. As self-consistency is approached, the set of linear equations to be solved for DIIS approaches linear dependency and it is useful to bias DIIS in favor of the SCF cycle that had the lowest energy using the factor DIISBfac. This is achieved by multiplying all diagonal elements of the DIIS matrix with this factor unless it is the Fock matrix/density which leads to the lowest energy. The default value for DIISBfac is 1.05.

The value of DIISMaxEq is the maximum number of old Fock matrices to remember. Values of 5-7 have been recommended, while other users store 10-15 Fock matrices. Should the standard DIIS not achieve convergence, some experimentation with this parameter can be worthwhile. In cases where DIIS causes problems in the beginning of the SCF, it may have to be invoked at a later stage. The start of the DIIS procedure is controlled by DIISStart. It has a default value of 0.2 Eh, which usually starts DIIS after 0-3 cycles. A different way of controlling the DIIS start is adjusting the value DIISMaxIt, which sets the maximum number of cycles after which DIIS will be started irrespective of the error value.

%scf
  # control of the DIIS procedure
  CNVDIIS     true  # default: true
  DIISStart   0.1   # default: 0.2
  DIISMaxIt   5     # default: 12
  DIISMaxEq   7     # default: 5
  DIISBFac    1.2   # default: 1.05
  DIISMaxC    15.0  # default: 10.0
  # more convenient:
  DIIS Start 0.1  MaxIt 5  MaxEq 7  BFac 1.2  MaxC 15.0  end
end

Note that for troublesome or lacking SCF convergence the TRAH algorithm should be used (see Sec. Trust-Region Augmented Hessian (TRAH) SCF). If not turned off explicitly, TRAH is switched on automatically whenever convergence problems are present by means of the AutoTRAH feature (see Sec. Trust-Region Augmented Hessian (TRAH) SCF).

7.7.5. An alternative DIIS algorithm: KDIIS

An alternative algorithm that makes use of the DIIS concept is called KDIIS (Kolmar’s DIIS[452, 453]) in ORCA. The KDIIS algorithm is designed to bring the orbital gradient of any energy expression to zero using a combination of DIIS extrapolation and first order perturbation theory. Thus, the method is diagonalization-free. In our hands it is superior to the standard DIIS algorithm in many cases, but not always. The algorithm is invoked with the keyword ! KDIIS and is available for RHF, UHF and CASSCF.

7.7.6. Approximate Second Order SCF (SOSCF)

SOSCF is an approximately quadratically convergent variant of the SCF procedure [264, 608]. The theory is relatively involved and will not be described here. In short – SOSCF computes an initial guess to the inverse orbital Hessian and then uses the BFGS formula in a recursive way to update orbital rotation angles. As information from a few iterations accumulates, the guess to the inverse orbital Hessian becomes better and better and the calculation reaches a regime where it converges superlinearly. As implemented, the procedure converges as well or slightly better than DIIS and takes a somewhat less time. However, it is also a lot less robust, so that DIIS is the method of choice for many problems (see also the description of the full second-order trust-region augmented Hessian (TRAH) procedure in the next section). On the other hand, SOSCF is useful when DIIS gets stuck at some error around \(\sim\) 0.001 or 0.0001. Such cases were the primary motive for the implementation of SOSCF into ORCA.

The drawback of SOSCF is the following: in the beginning of the SCF, the orbital gradient (the derivative of the total energy with respect to rotations that describe the mixing of occupied and virtual MOs) is large, so that one is far from the quadratic regime. In such cases, the procedure is not successful and may even wildly diverge. Therefore it is recommended to only invoke the SOSCF procedure in the very end of the SCF where DIIS may lead to “trailing” convergence. SOSCF is controlled by the variables SOSCFStart and SOSCFMaxIt. SOSCFStart is a threshold for the orbital gradient. When the orbital gradient, or equivalently the DIIS Error, fall below SOSCFStart, the SOSCF procedure is initiated. SOSCFMaxIt is the latest iteration to start the SOSCF even if the orbital gradient is still above SOSCFStart.

%scf
  # control of the SOSCF procedure
  CNVSOSCF    true  # default: false
  SOSCFStart  0.1   # default: 0.01
  SOSCFMaxIt  5     # default: 1000
  # more convenient:
  SOSCF Start 0.1  MaxIt 5  end
end

For many calculations on transition metal complexes, it is a good idea to be conservative in the startup criterion for SOSCF, it may diverge otherwise. A choice of 0.01 or lower is recommended.

7.7.7. Trust-Region Augmented Hessian (TRAH) SCF

The trust-region augmented Hessian (! TRAH) approach[64, 349, 381, 735] should be used if the standard SCF solver[264, 608, 703, 704] DIIS+SOSCF fails or is expected to fail. TRAH-SCF should converge for any system. Convergence to the electronic ground state is also guaranteed because information of the electronic Hessian is exploited.

The TRAH approach constructs a quadratic model of the SCF energy as a function of the orbital rotation parameters \(\mathbf{x}\),

\[\begin{aligned} E(\mathbf{x}) = E_0 + \mathbf{g}^T \mathbf{x} + \frac{1}{2}\mathbf{x}^{T} \mathbf{H} \mathbf{x}\end{aligned}\]

and minimizes \(E(\mathbf{x})\) w.r.t \(\mathbf{x}\) with the constraint that orbital rotations should lie within a trust region \(h\)

\[\begin{aligned} ||\mathbf{x}|| \le h \text{ .}\end{aligned}\]

Such a constraint minimization leads to the level-shifted Newton equations

\[\begin{aligned} \left( \mathbf{H} - \mu \mathbf{I} \right) \mathbf{x} = - \mathbf{g} \text{,}\end{aligned}\]

which have the two unknowns \((\mu, \mathbf{x})\). Instead of solving the level-shifted Newtons equations, the eigenvalues and eigenvectors of the scaled augmented Hessian are solved,

\[\begin{split}\begin{aligned} & \begin{pmatrix} 0 & \alpha \mathbf{g}^T \\ \alpha \mathbf{g} & \mathbf{H} \end{pmatrix} \begin{pmatrix} 1 \\ \tilde{\mathbf{x} } \end{pmatrix} = \mu \begin{pmatrix} 1 \\ \tilde{\mathbf{x} } \end{pmatrix} \text{,}\end{aligned}\end{split}\]

The TRAH eigenvalue equations are solved iteratively with the Davidson algorithm until the residual norm for \(\tilde{\mathbf{x} }\) is below a scaled (by TolFacMicro) norm of the electronic gradient. The scaling parameter \(\alpha\) is adjusted in every Davidson iteration (micro iteration) such that

\[\begin{aligned} ||\mathbf{x}|| \approx || \frac{1}{\alpha} \tilde{\mathbf{x} } || \le h\end{aligned}\]

by using a bisection search within \([\alpha_0,\alpha_1]\) and ensures that the orbital rotation (update) vectors are within the trust region. Once \(\mathbf{x}\) is found, the orbitals are updated and a new macro iteration starts with the SCF energy and electronic gradient computation \(\mathbf{g}\). TRAH terminates if the gradient norm \(||\mathbf{g}||\) is below a user-given threshold TolG.

The most time consuming steps of the algorithm are the computation of the electronic gradient \(\mathbf{g}\) and the linear transformations of the electronic Hessian \(\mathbf{H}\) with some trial vectors during the Davidson micro iterations. Both intermediates can be computed efficiently in the atomic-orbital basis using AO-Fock matrices as done for TD-DFT or CP-SCF. Hence, TRAH-SCF can be also used for very large molecules. However, in contrast to a standard DIIS approach, difference density matrices cannot be used which makes the shell pair screening based on Schwarz estimates and density matrices less effective for TRAH. To accelerate the various sigma vector computations, we choose for those steps in the micro iterations as for CP-SCF smaller grids for evaluation of XC functionals and semi-numerical exchange COSX, which is controlled via

%method
  Z_GridXC   1      // Lebedev Grid
  Z_IntAccXC 3.467  // eps parameter of radial grid
  Z_GridX    1      // Lebedev Grid
  Z_IntAccX  3.067  // eps parameter of radial grid
end

TRAH-SCF is currently implemented for restricted closed-shell (RHF and RKS) and unrestricted open-shell determinants (UHF and UKS) and can be accelerated with RIJ, RIJONX, RIJK, or RIJCOSX. Solvation effects can also be accounted for with the C-PCM model. The implementation is also parallelized with MPI. A restricted open-shell implementation (ROHF and ROKS) is not yet available.

The default preconditioner (diag) for the Davidson algorithm uses approximate matrix element of the diagonal Hessian. We have also added an improved preconditioner (full) that uses the exact orbital Hessian for a subset of the most important occupied-virtual MO pairs (with smallest orbital-energy difference). This number PreconMaxRed (default 250) cannot be too large.

%scf
  trah
    Precond         Diag   # DIAG, FULL, or NONE
    PreconMaxRed    250    # maximum dimension of FULL Hessian for preconditioning
  end
end

Otherwise, additional computational bottlenecks would be introduced when transforming the two-electron integrals in the MO basis or when diagonalizing this reduced-space Hessian. Note that for the integral transformation the RI approximation is used and an auxiliary /C basis must be provided. Please also note that, so far, we have only implemented an XC Hessian contribution for LDA functionals. From our experience, the “full” precondition is very advantageous for RHF and UHF calculations of small molecules but does not provide any advantage for other (TRAH-)SCF calculations.

In cases for which the conventional SCF procedures (DIIS/KDIIS/SOSCF) struggle, we invoke TRAH-SCF automatically (AutoTRAH). For this purpose, we perform a linear interpolation of the norm of the electronic gradient on the \({\log}_{10}\) scale after a minimum number of SCF iterations AutoTRAHIter (default 20). The number of iterations for interpolations is controlled by AutoTRAHNInter (default 10). If the slope \(s\) of the interpolated gradient norm \(10^{-s}\) becomes smaller than AutoTRAHTol (default 1.125), conventional SCF is shut down and TRAH-SCF start from the current set of MOs. Those parameters were optimized for a benchmark set with the purpose to minimize the calculation times even for SCF calculations that are hard to converge. There is not really a need to modify the AutoTRAH parameters except for turning TRAH off entirely.

! NoTRAH

The accuracy of the SCF calculation is controlled by via the simple keywords NORMALSCF, LOOSESCF, etc. The accuracy threshold that checks the gradient norm for TRAH-SCF calculations is is read from the SCF input block. Note that checking the Frobensius norm of \(||\mathbf{g}||\) is the only convergence check in TRAH.

%scf
  TolG 1.e-6  # gradient norm threshold (converging macro iterations)
end

Below a complete list of input parameters is given for TRAH. Please note that all parameters influence the convergence and should not be changed carelessly!

%scf
 # AutoTRAH parameter
 AutoTRAH       true
 AutoTRAHTol    1.125
 AutoTRAHIter   20
 AutoTRAHNInter 10

 # TRAH parameter
   trah
     MaxRed          16     # maximum number of Davidson micro iterations
     NStart          2      # number of start vectors for Davidson (at least 2)
     TolFacMicro     0.1    # Scaling factor for Davidson convergence
                            #  threshold =  TolFacMicro * || G ||    
     MinTolMicro     1.e-2  # minimum accuracy of micro iterations
     QuadRegionStart 1.e-3  # start Newton-Raphson if || G || < QuadRegionStart 
     tradius         0.4    # initial trust radius
     AlphaMin        0.1    # lower bound of gradient scaling parameter
     AlphaMax        1000.  # upper bound of gradient scaling parameter
     Randomize       true   # add white noise to Davidson start vectors
     PseudoRand      false  # use pseudo random numbers for comparibility
     MaxNoise        1.e-2  # maximum random number magnitude
     OrbUpdate       Taylor # orbital update algorithm (TAYLOR or CAYLEY)
     InactiveMOs     Canonical # CANONICAL or NOTSET
     Precond         Diag   # DIAG, FULL, or NONE
     PreconMaxRed    250    # maximum dimension of FULL Hessian for preconditioning
   end
 end

OBS.: The maximum number of macro iterations is defined by MAXITER under the %SCF block.

7.7.8. Finite Temperature HF/KS-DFT

A finite temperature can be used to apply a Fermi-like occupation number smearing over the orbitals of the system, which may sometimes help to get convergence of the SCF equations in near hopeless cases. Through the smearing, the electrons are distributed according to Fermi statistics among the available orbitals. The “chemical potential” is found through the condition that the total number of electrons remains correct. Gradients can be computed in the presence of occupation number smearing.

%scf SmearTemp 5000  # ``temperature'' in Kelvin
 end

Note

7.7.8.1. Fractional Occupation Numbers

Only a very basic implementation of fractional occupation numbers is presently provided. It is meant to deal with orbitally degenerate states in the UHF/UKS method. Mainly it was implemented to avoid symmetry breaking in DFT calculations on orbitally degenerate molecules and atoms. The program checks the orbital energies of the initial guess orbitals, finds degenerate sets and averages the occupation numbers among them. Currently the criterion for degenerate orbitals is \(10^{-3}\) Eh. The fractional occupation number option is invoked by:

%scf
  FracOcc true
end

Clearly, the power of fractional occupation numbers goes far beyond what is presently implemented in the program and future releases will likely make more use of them. The program prints a warning whenever it uses fractional occupation numbers. The fractionally occupied orbitals should be checked to ensure they are actually the intended ones.

Note

  • Using GuessMode = CMatrix will cause problems because there are no orbital energies for the initial guess orbitals. The program will then average over all orbitals — which makes no sense at all.

7.7.8.2. Fractional Occupation Number Weighted Electron Density (FOD)

Many approximate QC methods do not yield reliable results for systems with significant static electron correlation (SEC) but it is often difficult to predict if the system in question suffers from SEC or not. Existing scalar SEC diagnostics (e.g., the \(T_1\) diagnostic) do not provide any information where the SEC is located in the molecule. Furthermore, often quite expensive calculations have to be performed first (e.g., CCSD) in order to judge the reliability of the results based on a single number. Molecular systems with strong SEC (e.g. covalent bond-breaking, biradicals, open-shell transition metal complexes) are usually characterized by small energy gaps between frontier orbitals, and hence, the appearance of many equally important determinants in their electronic wavefunction. This finding is used in the FOD analysis[327] which is based on finite temperature KS-DFT where the fractional occupation numbers are determined from the Fermi distribution (“Fermi smearing”)

\[\begin{gathered} f_i=\frac{1}{e^{(\varepsilon_i-E_F)/kT_{el} }+1}\end{gathered}\]

The central quantity of the FOD analysis is the fractional occupation number weighted electron density (\(\rho^{FOD}\)), a real-space function of the position vector \(r\):

\[\begin{gathered} \rho^{FOD}(r)=\sum_i^N (\delta_1-\delta_2f_i) |\varphi_i(r)|^2\end{gathered}\]

(\(\delta_1\) and \(\delta_2\) are unity if the level is lower than \(E_F\) while they are 0 and \(-1\), respectively, for levels higher than \(E_F\)). The \(f_i\) represent the fractional occupation numbers (0\(\leq f_i\leq\)1; sum over all electronic single-particle levels obtained by solving self-consistently the KS-SCF equations minimizing the free-electronic energy).

\(\rho^{FOD}(r)\) can be plotted using a pre-defined contour surface value (see  FOD plots). FOD plots only show the contribution of the ‘hot’ (strongly correlated) electrons and can thus be used to choose a suitable QC method for the system in question based on some rules of thumb (see FOD plots). Mulliken reduced orbital charges based on \(\rho^{FOD}(r)\) (see Mulliken Population Analysis) offer a fast alternative to get the information of the FOD plot.

The integration of \(\rho^{FOD}\) over all space yields as additional information a single size-extensive number termed \(N_{FOD}\) which correlates well with other scalar SEC diagnostics and can be used to globally quantify SEC effects in the molecule.

\(\rho^{FOD}\) (and \(N_{FOD}\)) strongly depend on the orbital energy gap which itself depends almost linearly on the amount of the non-local Fock exchange admixture \(a_x\). The following (empirical) function of the optimal electronic temperature \(T_{el}\) on \(a_x\)

\[\begin{gathered} T_{el}=20000\, { \textup K} \times a_x +5000\, { \textup K}\end{gathered}\]

is used to ensure that similar results of the FOD analysis are obtained with various functionals. For example, the SmearTemp has to be 5000 K for TPSS (\(a_x\) = 0), 9000 K for B3LYP (\(a_x\) = 20%), 10000 K for PBE0 (\(a_x\) = 25%), and 15800 K for M06-2x (\(a_x\) = 54%). The result of the FOD analysis is not strongly dependent on the employed basis set (see supplementary information of the original publication[327]). TPSS/def2-TZVP/TightSCF was chosen as the default since it is fast and robust. The FOD analysis is a very efficient and practicable tool to get information about the amount and localization of SEC in the system of question. It is called by a simple keyword:

# ground state of p-benzyne
! FOD 

* xyz 0 1
C 0.0000000  1.2077612  0.7161013
C 0.0000000  0.0000000  1.3596219
C 0.0000000 -1.2077612  0.7161013
C 0.0000000 -1.2077612 -0.7161013
C 0.0000000  0.0000000 -1.3596219
C 0.0000000  1.2077612 -0.7161013
H 0.0000000  2.1606260  1.2276695
H 0.0000000 -2.1606260  1.2276695
H 0.0000000 -2.1606260 -1.2276695
H 0.0000000  2.1606260 -1.2276695
*

The respective output reads:

-------------------------------------------------------------------------------------------
                                      ORCA LEAN-SCF
                              memory conserving SCF solver
-------------------------------------------------------------------------------------------

----------------------------------------D-I-I-S--------------------------------------------
Iteration    Energy (Eh)           Delta-E    RMSDP     MaxDP     DIISErr   Damp  Time(sec)
-------------------------------------------------------------------------------------------
               ***  Starting incremental Fock matrix formation  ***
    1    -230.8982516003082139     0.00e+00  5.01e-03  1.01e-01  1.12e-01  0.700   1.7
Warning: op=0 Small HOMO/LUMO gap (  -0.021) - skipping pre-diagonalization
         Will do a full diagonalization
    2    -230.9463607195993120    -4.81e-02  1.15e-03  2.59e-02  4.06e-02  0.700   1.6
                               ***Turning on AO-DIIS***
  ... etc.
   12    -231.0033984839932089    -5.02e-09  3.33e-07  7.37e-06  7.95e-06  0.000   1.1
                 **** Energy Check signals convergence ****

FOD:
 Fermi smearing:E(HOMO(Eh)) = -0.201252 MUE = -0.179318 gap=   1.119 eV

 N_FOD =  0.920364

The default functional and basis set are TPSS and def2-TZVP respectively. If the FOD analysis should be done employing a different functional, one has to explicitly specify the functional and basis set in the simple keyword line and adjust the SmearTemp accordingly.

# ozone molecule
! B3LYP def2-TZVP TightSCF

%scf
 SmearTemp 9000
end

* xyz 0 1
  O   0.00000000017911      0.00000000000000      0.43702029765776
  O  -1.09512651993192      0.00000000000000     -0.21851064888325
  O   1.09512651975281      0.00000000000000     -0.21851064877451
*

The FOD analysis may also be useful for finding a suitable active space for e.g. CASSCF calculations.

Note

  • The FOD analysis will be always printed (including Mulliken reduced orbital charges based on \(\rho^{FOD}\)) if SmearTemp \(>\) 0 K \(\rho^{FOD}\) is stored on disk in the file Basename.scfp_fod which is included in the general Basename.densities container).

  • Since the \(\hat{S}^2\) expectation value is not defined for fractional occupation numbers, its printout is omitted.