HF: Hartree–Fock Theory¶

the Hartree–Fock energy is, by Slater’s rules,

\[E_<\mathrm> = \langle \Psi_0 | \hat H | \Psi_0 \rangle = \sum_ \langle i | \hat h | i \rangle + \frac 1 2 \sum_ [ii|jj] - [ij|ji]\]

Here \(H\) is the AO-basis one-electron potential, encapsulating both electron-nuclear attraction and kinetic energy,

\(D\) is the AO-basis density matrix, build from the occupied orbital coefficients,

and \(F\) is the Fock matrix, which is the effective one-body potential at the current value of the density,

Here the tensor \((\mu\nu|\lambda\sigma)\) is an AO Electron-Repulsion Integral (ERI) in chemists’ notation,

\[(\mu\nu|\lambda\sigma) = \iint_<\mathbb^6> \phi_ <\mu>(\vec r_1) \phi_ (\vec r_1) \frac> \phi_ <\lambda>(\vec r_2) \phi_ (\vec r_2) \ \mathrm^3 r_1 \ \mathrm^3 r_2.\]

The MO coefficients are found as the generalized eigenvectors of the Fock Matrix,

The eigenvalues \(\epsilon\) are the orbital energies, and the metric matrix \(S\) is the AO-basis overlap matrix

Note that the Fock Matrix depends on the density (both alpha and beta), and therefore the orbitals. Because of this, SCF is a nonlinear procedure, which terminates when the generating orbitals are self-consistent with the Fock matrix they generate.

The formation of the Coulomb matrix \(J\) and the exchange matrix \(K^\) dominate the computational effort of the SCF procedure. For very large systems, diagonalization of the Fock matrix can also present a significant hurdle.

Minimal Input¶

Minimal input for a Hartree–Fock computation is a molecule block, basis set option, and a call to energy('scf') :

molecule  He > set basis sto-3g energy('scf') 

This will run a Restricted Hartree–Fock (RHF) on neutral singlet Helium in \(D_\) spatial symmetry with a minimal STO-3G basis, 1.0E-6 energy and density convergence criteria (since single-point, see SCF Convergence & Algorithm ), a DF ERI algorithm, symmetric orthogonalization, DIIS, and a core Hamiltonian guess. For more information on any of these options, see the relevant section below.

Spin/Symmetry Treatment¶

P SI 4 implements the most popular spin specializations of Hartree–Fock theory, including:

Restricted Hartree–Fock (RHF) [Default]

Appropriate only for closed-shell singlet systems, but twice as efficient as the other flavors, as the alpha and beta densities are constrained to be identical.

Unrestricted Hartree–Fock (UHF)

Appropriate for most open-shell systems and fairly easy to converge. The spatial parts of the alpha and beta orbitals are fully independent of each other, which allows a considerable amount of flexibility in the wavefunction. However, this flexibility comes at the cost of spin symmetry; UHF wavefunctions need not be eigenfunctions of the \(\hat S^2\) operator. The deviation of this operator from its expectation value is printed on the output file. If the deviation is greater than a few hundredths, it is advisable to switch to a ROHF to avoid this “spin-contamination” problem.

Restricted Open-Shell Hartree–Fock (ROHF)

Appropriate for open-shell systems where spin-contamination is problem. Sometimes more difficult to converge, and assumes uniformly positive spin polarization (the alpha and beta doubly-occupied orbitals are identical).

Constrained Unrestricted Hartree–Fock (CUHF)

A variant of ROHF that starts from a UHF ansatz and is therefore often easier to converge.

These can be invoked by the REFERENCE keyword, which defaults to RHF . The charge and multiplicity may either be specified in the molecule definition:

molecule h  0 2 # Neutral doublet H > 

or, dynamically, by setting the relevant attributes in the Python molecule object:

h.set_molecular_charge(0) h.set_multiplicity(2) 

Abelian spatial symmetry is fully supported in P SI 4 and can be used to obtain physical interpretation of the molecular orbitals, to assist in difficult convergence cases, and, in some methods, to obtain significant performance gains. The point group of the molecule is inferred when reading the molecule section, and may be overridden by the symmetry flag, as in:

molecule h  0 2 H symmetry c1 > 

or by the reset_point_group Python molecule attribute:

h.reset_point_group('c2v') 

During the SCF procedure, the occupation of orbitals is typically determined by the Aufbau principal across all spatial symmetries. This may result in the occupation shifting between iterations. If the occupations are known a priori, they may be clamped throughout the procedure by using the DOCC and SOCC options. For instance, all good quantum chemists know that \(C_\) water is actually,:

molecule h2o  0 1 O H 1 1.0 H 1 1.0 2 104.5 > set  docc [3, 0, 1, 1] # 1A1 2A1 1B1 3A1 1B2 basis cc-pvdz > energy('scf') 

Broken Symmetry¶

For certain problems, such diradicals, allowing the spin-up and spin-down orbitals to differ in closed-shell computations can be advantageous; this is known as symmetry breaking. The resulting unrestricted wavefunction will often provide superior energetics, due to the increased flexibility, but it will suffer non-physical spin contamination from higher multiplicity states. A convenient approach to break symmetry is to perform a UHF or UKS calculation with the guess HOMO and LUMO orbitals mixed. Mixing of the guess orbitals can be requested by setting the GUESS_MIX keyword to true:

set reference uhf set guess_mix true energy('scf') 

Orthogonalization¶

One of the first steps in the SCF procedure is the determination of an orthogonal basis (known as the OSO basis) from the atomic orbital basis (known as the AO basis). The Molecular Orbital basis (MO basis) is then built as a particular unitary transformation of the OSO basis. In P SI 4 , the determination of the OSO basis is accomplished via either symmetric, canonical, or partial Cholesky orthogonalization.

Symmetric orthogonalization uses the symmetric inverse square root of the overlap matrix for the orthogonalization matrix. Use of symmetric orthogonalization always yields the same number of OSO functions (and thereby MOs) as AO functions. However, this may lead to numerical problems if the overlap matrix has small eigenvalues, which may occur for large systems or for systems where diffuse basis sets are used.

This problem may be avoided by using canonical orthogonalization, in which an asymmetric inverse square root of the overlap matrix is formed, with numerical stability enhanced by the elimination of eigenvectors corresponding to very small eigenvalues. As a few combinations of AO basis functions may be discarded, the number of canonical-orthogonalized OSOs and MOs may be slightly smaller than the number of AOs.

When the basis set is too overcomplete, the eigendecomposition of the overlap matrix is no longer numerically stable. In this case the partial Cholesky decomposition can be used to pick a subset of basis functions that span a sufficiently complete set, see [Lehtola:2019:241102] and [Lehtola:2020:032504] . This subset can then be orthonormalized as usual; the rest of the basis functions are hidden from the calculation. The Cholesky approach allows reaching accurate energies even in the presence of significant linear dependencies [Lehtola:2020:134108] .

In P SI 4 , symmetric orthogonalization is used by default, unless the smallest overlap eigenvalue falls below the user-supplied double option S_TOLERANCE , which defaults to 1E-7. If the smallest eigenvalue is below this cutoff, canonical orthogonalization is forced, and all eigenvectors corresponding to eigenvalues below the cutoff are eliminated.

If the eigendecomposition is detected to be numerically unstable - the reciprocal condition number of the overlap matrix to be smaller than the machine epsilon - the partial Cholesky decomposition is undertaken until S_CHOLESKY_TOLERANCE , which defaults to 1E-8.

Use of symmetric, canonical, and partial Cholesky orthogonalization can be forced by setting the S_ORTHOGONALIZATION option to SYMMETRIC , CANONICAL , or PARTIALCHOLESKY , respectively.

Note that in practice, the MOs and OSOs are built separately within each irrep from the symmetry-adapted combinations of AOs known as Unique Symmetry Orbitals (USOs). For canonical orthogonalization, this implies that the number of MOs and OSOs per irrep may be slightly smaller than the number of USOs per irrep.

A contrived example demonstrating OSOs/MOs vs. AOs with symmetry is shown below:

molecule h2o  0 1 O H 1 1.0 H 1 1.0 2 104.5 symmetry c2 # Two irreps is easier to comprehend > set  s_tolerance 0.0001 # Set an unreasonably tight # tolerance to force canonical basis aug-cc-pv5z # This diffuse basis will have # small-ish eigenvalues for even H2O print 3 > energy('scf') 

==> Pre-Iterations  ------------------------------------------------------- Irrep Nso Nmo Nalpha Nbeta Ndocc Nsocc ------------------------------------------------------- A 145 145 0 0 0 0 B 142 142 0 0 0 0 ------------------------------------------------------- Total 287 287 5 5 5 0 ------------------------------------------------------- . Minimum eigenvalue in the overlap matrix is 1.6888063568E-05. Using Canonical Orthogonalization with cutoff of 1.0000000000E-04. Irrep 0, 1 of 145 possible MOs eliminated. Irrep 1, 2 of 142 possible MOs eliminated. Overall, 3 of 287 possible MOs eliminated. 

In this example, there are 287 AO basis functions after spherical harmonics are applied. These are used to produce 287 symmetry adapted USOs, 145 of which are assigned to irrep A, and 142 of which are assigned to irrep B. Within irrep A, 144 OSOs fall above the eigenvalue cutoff, and within irrep B 140 OSOs fall above the eigenvalue cutoff. In total, 284 molecular orbitals are chosen from 287 AOs/USOs.

Initial Guess¶

In each step of the SCF procedure, a new Fock or Kohn–Sham potential is built according to the previous density, following which the potential is diagonalized to produce new molecular orbitals, from which a new density is computed. This procedure is continued until either convergence is reached or a preset maximum number of iterations is exceeded. Convergence is determined by both change in energy and root-mean-square change in density matrix values, which must be below the user-specified E_CONVERGENCE and D_CONVERGENCE , respectively. The maximum number of iterations is specified by the MAXITER option. It should be noted that SCF is a chaotic process, and, as such, often requires careful selection of initial orbitals and damping during iterations to ensure convergence. This is particularly likely for large systems, metallic systems, multireference systems, open-shell systems, anions, and systems with diffuse basis sets.

For initial orbital selection, several options are available. These include:

Diagonalization of the core Hamiltonian, removing even mean-field electron repulsion. Simple, but often too far from the final solution for larger systems. This is the default for single atoms.

Superposition of Atomic Densities. Builds the initial density as the spin-averaged sum of atomic UHF computations in the current basis. If an open-shell system, uniform scaling of the spin-averaged density matrices is performed. If orbitals are needed (e.g., in density fitting), a partial Cholesky factorization of the density matrices is used. Often extremely accurate, particularly for closed-shell systems. This is the default for systems of more than one atom.

Natural orbitals from Superposition of Atomic Densities. Similar to the above, but it forms natural orbitals from the SAD density matrix to get proper orbitals which are used to start the calculation, see [Lehtola:2019:1593] .

A generalized Wolfsberg-Helmholtz modification of the core Hamiltonian matrix. Usually less accurate than the core guess: the latter is exact for one-electron systems, GWH is not; see [Lehtola:2019:1593] ).

An extended Hückel guess based on on-the-fly atomic UHF calculations alike SAD, see [Lehtola:2019:1593] .

Like HUCKEL, an extended Hückel guess based on on-the-fly atomic UHF calculations alike SAD, see [Lehtola:2019:1593] . This variant employs an updated rule for the generalized Wolfsberg-Helmholz formula from [Ammeter:1978:3686] .

Read the previous orbitals from a wfn file, casting from one basis to another if needed. Useful for starting anion computations from neutral orbitals, or after small geometry changes. At present, casting from a different molecular point group is not supported. This becomes the default for the second and later iterations of geometry optimizations.

Superposition of Atomic Potentials. This is essentially a modification of the core Hamiltonian, which includes screening effects by using a radially screened effective atomic charge. The screening effects have been calculated at the complete basis set limit with finite-element calculations, see [Lehtola:2019:25945] and [Lehtola:2020:012516] . The guess and its implementation have been described in [Lehtola:2019:1593] . The guess is evaluated on a DFT quadrature grid, so the guess energy depends slightly on the used DFT quadrature. The current implementation is based on exchange-only local density calculations that are but nanohartree away from the complete basis set limit [Lehtola:2020:012516] .

Superposition of Atomic Potentials, but using error function based fits to the atomic radial potentials as discussed in [Lehtola:2020:144105] . The main difference to the SAP guess discussed above [Lehtola:2019:25945] is that the SAPGAU scheme is analytic, and can be efficiently formed in terms of three-center two-electron integrals [Lehtola:2020:144105] . The potential in the SAPGAU scheme is passed with the SAPGAU_BASIS keyword. The default potential is given by the large fit to the HelFEM potential, sap_helfem_large, described in [Lehtola:2020:144105] . Note that this guess is known in the DIRAC program as .SCRPOT and in the ERKALE program as SAPFIT.

These are all set by the GUESS keyword. Also, an automatic Python procedure has been developed for converging the SCF in a small basis, and then casting up to the true basis. This can be done by adding BASIS_GUESS = SMALL_BASIS to the options list. We recommend the 3-21G or pcseg-0 basis for the small basis due to its efficient mix of flexibility and compactness. An example of performing an RHF solution of water by SAD guessing in a 3-21G basis and then casting up to cc-pVTZ is shown below:

molecule h2o  0 1 O H 1 1.0 H 1 1.0 2 104.5 > set  basis cc-pvtz basis_guess 3-21G guess sad > energy('scf') 

Restarting the SCF¶

Reading orbital data from a previous calculations is done via the restart_file option, where the actual file is a serialized wfn object (see saving the wfn ) By default, the orbital data file of the converged SCF( psi.PID.name.180.npy ) is deleted after P SI 4 exits or the clean() function is called. The orbital guess is automatically set to READ when restart_file is given a wfn file. To write the orbitals after every iteration and keep the orbitals from the last iteration, the write_orbitals options is available:

energy('scf', write_orbitals='my_mos'), 

which writes a Wavefunction object converted (serialized) to a numpy file called my_mos.npy . The restart can then be done as follows:

energy('scf', restart_file='my_mos') 

Specifying the .npy suffix when writing and reading restart files is optional.

Alternatively, the restart can also be done from any previously saved wfn object.

energy, scf_wfn = energy('scf',return_wfn=True) scf_wfn.to_file('my_wfn') energy('scf', restart_file='my_wfn') 

For advanced users manipulating or writing custom wavefunction files, note that P SI 4 expects the numpy file on disk to have the .npy extension, not, e.g., .npz .

Convergence Stabilization¶

A summary of Psi’s supported convergence stabilization techniques is presented below:

DIIS [On by Default]

DIIS uses previous iterates of the Fock matrix together with an error criterion based on the orbital gradient to produce an informed estimate of the next Fock Matrix. DIIS is almost always necessary to converge the SCF procedure and is therefore turned on by default. In rare cases, the DIIS algorithm may need to be modified or turned off altogether, which may be accomplished via options .

ADIIS [On by Default]

ADIIS uses previous iterates of the Fock and density matrices to produce an informed estimate of the next Fock matrix. ADIIS estimates are based on minimizing an energy estimate rather than zeroing the residual, so this performs best in the early iterations. By default, Psi will start using ADIIS before blending the ADIIS step with the DIIS step, eventually using the pure DIIS step. The closely-related EDIIS procedure may be used instead by setting SCF_INITIAL_ACCELERATOR . This is formally identical to ADIIS for HF, but the methods will differ for more general DFT.

MOM [Off by Default]

MOM was developed to combat a particular class of convergence failure: occupation flipping. In some cases, midway though the SCF procedure, a partially converged orbital which should be occupied in the fully-optimized SCF solution has a slightly higher orbital eigenvalue than some other orbital which should be destined to be a virtual orbital. This results in the virtual orbital being spuriously occupied for one or more iterations. Sometimes this resolves itself without help, other times the occupation flips back and forth between two, four, or more orbitals. This is typically visible in the output as a non-converging SCF which eventually settles down to steady oscillation between two (or more) different total energies. This behavior can be ameliorated by choosing occupied orbitals by “shape” instead of by orbital eigenvalue, i.e., by choosing the set of new orbitals which looks most like some previously known “good” set. The “good” set is typically the occupied orbitals from one of the oscillating iterations with the lowest total energy. For an oscillating system where the lowest total energy occurs on iterations \(N,N+2,\ldots\) , invoking MOM_START N can often rescue the convergence of the SCF. MOM can be used in concert with DIIS, though care should be taken to not turn MOM on until the oscillatory behavior begins.

Damping [Off by Default]

In some cases, a static mixing of Fock Matrices from adjacent iterations can quench oscillations. This mixing, known as “damping” can be activated by setting the DAMPING_PERCENTAGE keyword to a nonzero percent. Damping is turned off when the DIIS error is smaller than DAMPING_CONVERGENCE .

Level shifting [Off by default]

A commonly used alternative to damping is to use level shifting, which decreases the mixing of occupied and unoccupied orbitals in the SCF update by moving the unoccupied orbitals up in energy. It can be shown that the SCF procedure always converges with a suitably large level shift; however, the larger the shift is, the slower the convergence becomes, and the calculation may end up converging onto a higher lying SCF solution. Because of this, in practice level shifting is most useful in the initial phase of the calculation to reduce the orbital error enough for DIIS to work well. The level shift is controlled by the parameter LEVEL_SHIFT , and it is turned off when the DIIS error is smaller than LEVEL_SHIFT_CUTOFF . Reasonable values for the shift and convergence threshold are 5.0 and 1e-2, respectively.

SOSCF [Off by Default]

ERI Algorithms¶

The key difficulty in the SCF procedure is treatment of the four-index ERI contributions to the Fock Matrix. A number of algorithms are available in P SI 4 for these terms. The algorithm is selected by the SCF_TYPE keyword. Most consist of a single algorithm applied to the construction of both the Coulomb and Exchange parts of the Fock Matrix:

An out-of-core, presorted algorithm using exact ERIs. Quite fast for a zero-error algorithm if enough memory is available. Integrals are generated only once, and symmetry is utilized to reduce number of integrals.

An out-of-core, unsorted algorithm using exact ERIs. Overcomes the memory bottleneck of the current PK algorithm. Integrals are generated only once, and symmetry is utilized to reduce number of integrals.

A threaded, sieved, integral-direct algorithm, with full permutational symmetry. This algorithm is brand new, but seems to be reasonably fast up to 1500 basis functions, uses zero disk (if DF pre-iterations are turned off), and can obtain significant speedups with negligible error loss if INTS_TOLERANCE is set to 1.0E-8 or so.

A density-fitted algorithm designed for computations with thousands of basis functions. This algorithm is highly optimized, and is threaded with a mixture of parallel BLAS and OpenMP. Note that this algorithm should use the -JKFIT series of auxiliary bases, not the -RI or -MP2FIT bases. The default guess for auxiliary basis set should work for most bases, otherwise the DF_BASIS_SCF keyword can be used to manually specify the auxiliary basis. This algorithm is preferred unless either absolute accuracy is required [ \(\gtrsim\) CCSD(T)] or a -JKFIT auxiliary basis is unavailable for the orbital basis/atoms involved.

A threaded algorithm using approximate ERIs obtained by Cholesky decomposition of the ERI tensor. The accuracy of the Cholesky decomposition is controlled by the keyword CHOLESKY_TOLERANCE . This algorithm is similar to the DF algorithm, but it is not suitable for gradient computations. The algorithm to obtain the Cholesky vectors is not designed for computations with thousands of basis functions.

P SI 4 also features the capability to use “composite” Fock matrix build algorithms - arbitrary combinations of specialized algorithms that construct either the Coulomb or the Exchange matrix separately. In general, since separate Coulomb and Exchange matrix build algorithms exploit properties specific to their respective matrix, composite algorithms display lower scaling factors than their combined Fock build counterparts. However, composite algorithms also introduce redundant ERI computations into the calculation. Therefore, composite Fock build algorithms tend to perform better for larger systems, but worse for smaller systems. Arbitrary composite algorithms can be accessed by setting SCF_TYPE to J_alg+K_alg , where J_alg and K_alg are the names of the separate Coulomb and Exchange construction algorithms to use, respectively. Alternatively, if one is using DFT with non-hybrid functionals, a composite Coulomb construction algorithm can be specified solo by setting SCF_TYPE to J_alg , without the need to set an associated K_alg .

Specialized algorithms available to construct the Coulomb term within a composite framework are as follows:

An integral-direct algorithm constructing the Coulomb term based on [Weigend:2002:4285] The DFDIRJ algorithm combines the benefits of integral-direct SCF approaches with that of density-fitting. Specifically, DFJ utilizes no I/O and displays strong performance with large system size through a combination of effective parallelization and utilization of density-fitting to minimize ERI computational cost. See the Integral-Direct Density-Fitted Coulomb Construction section for more information.

Specialized algorithms available to construct the Exchange term within a composite framework are as follows:

An algorithm based on the semi-numerical “chain of spheres exchange” (COSX) approach described in [Neese:2009:98] . The coulomb term is computed with a direct density-fitting algorithm. The COSX algorithm uses no I/O, scales well with system size, and requires minimal memory, making it ideal for large systems and multi-core CPUs. See COSX Exchange for more information.

An implementation of the linear-scaling “Linear Exchange” (LinK) algorithm described in [Ochsenfeld:1998:1663] . The LINK algorithm provides many of the benefits of integral-direct SCF algorithms, including no disk I/O, low memory usage, and effective parallelization. Additionally, the LINK implementation scales well with system size while simultaneously providing a formally-exact computation of the Exchange term. See Linear Exchange for more information.

In some cases the above algorithms have multiple implementations that return the same result, but are optimal under different molecules sizes and hardware configurations. Psi4 will automatically detect the correct algorithm to run and only expert users should manually select the below implementations. The DF algorithm has the following two implementations

A DF algorithm optimized around memory layout and is optimal as long as there is sufficient memory to hold the three-index DF tensors in memory. This algorithm may be faster for builds that require disk if SSDs are used.

A DF algorithm (the default DF algorithm before Psi4 1.2) optimized to minimize Disk IO by sacrificing some performance due to memory layout.

Note that these algorithms have both in-memory and on-disk options, but performance penalties up to a factor of 2.5 can be found if the incorrect algorithm is chosen. It is therefore highly recommended that the keyword “DF” be selected in all cases so that the correct implementation can be selected by P SI 4 ‘s internal routines. Expert users can manually switch between MEM_DF and DISK_DF; however, they may find documented exceptions during use as several post SCF algorithms require a specific implementation. Additionally, expert users can manually switch between the in-memory and on-disk options within MEM_DF or DISK_DF using the SCF_SUBTYPE option. Using SCF_SUBTYPE = AUTO , where P SI 4 automatically selects the in-memory or on-disk option for MEM_DF/DISK_DF based on memory and molecule, is the default and recommended option. However, the in-memory or on-disk algorithms for MEM_DF and DISK_DF can be forced by using SCF_SUBTYPE = INCORE or SCF_SUBTYPE = OUT_OF_CORE , respectively. Note that an exception will be thrown if SCF_SUBTYPE = INCORE is used without allocating sufficient memory to P SI 4 .

For some of these algorithms, Schwarz and/or density sieving can be used to identify negligible integral contributions in extended systems. To activate sieving, set the INTS_TOLERANCE keyword to your desired cutoff (1.0E-12 is recommended for most applications). To choose the type of sieving, set the SCREENING keyword to your desired option. For Schwarz screening, set it to SCHWARZ , for CSAM, CSAM , and for density matrix-based screening, DENSITY .

Uses the Cauchy-Schwarz inequality to calculate an upper bounded value of a shell quartet,

An extension of the Schwarz estimate that also screens over the long range 1/r operator, described in [Thompson:2017:144101] .

An extension of the Schwarz estimate that also screens over elements of the density matrix. For the RHF case, described in [Haser:1989:104]

When using density-matrix based integral screening, it is useful to build the J and K matrices incrementally, also described in [Haser:1989:104] , using the difference in the density matrix between iterations, rather than the full density matrix. To turn on this option, set INCFOCK to true .

We have added the automatic capability to use the extremely fast DF code for intermediate convergence of the orbitals, for SCF_TYPE DIRECT . At the moment, the code defaults to cc-pVDZ-JKFIT as the auxiliary basis, unless the user specifies DF_BASIS_SCF manually. For some atoms, cc-pVDZ-JKFIT is not defined, so a very large fitting basis of last resort will be used. To avoid this, either set DF_BASIS_SCF to an auxiliary basis set defined for all atoms in the system, or set DF_SCF_GUESS to false, which disables this acceleration entirely.

Integral-Direct Density-Fitted Coulomb Construction¶

The Resolution of the Identity (RI) can be used to decompose the normally 4-center ERI tensor into a combination of 3-center and 2-center components. By reducing the dimensionality of the ERI tensor, application of the RI (often referred to as density-fitting, or DF) can be used to greatly speed up SCF calculations. The reduction in ERI tensor rank also makes DF an appealing option for conventional SCF calculations, where the ERIs are stored in core or on disk. However, even when using DF, I/O becomes a significant bottleneck for systems of a sufficient size when performing conventional SCF calculations. In principle, though, DF approaches can be utilized in an integral-direct context, gaining the benefits of DF methods without suffering the I/O bottlenecks that conventional DF methods will eventually run into. One such approach, outlined by Weigend in [Weigend:2002:4285] , is available for use in Psi4 for the separate construction of the Coulomb contribution to the Fock matrix. This implementation can be used alongside Psi4’s separate Exchange construction algorithms for composite Fock matrix construction by using the keyword DFDIRJ as the Coulomb construction algorithm when specifying SCF_TYPE to use a composite algorithm combination ( DFDIRJ+K_alg in general, or DFDIRJ for DFT with non-hybrid functionals).

DFDIRJ supports multiple capabilities to improve performance. Specifically, DFDIRJ allows for a combination of density-matrix based ERI screening (set SCREENING to DENSITY ) and incremental Fock matrix construction (set INCFOCK to TRUE ). These two, when combined, enable more aggressive screening of ERI contributions to the Coulomb matrix and thus greatly improve performance.

COSX Exchange¶

The semi-numerical COSX algorithm described in [Neese:2009:98] evaluates two-electron ERIs analytically over one electron coordinate and numerically over the other electron coordinate, and belongs to the family of pseudospectral methods originally suggested by Friesner. In COSX, numerical integration is performed on standard DFT quadrature grids, which are described in DFT: Density Functional Theory . Both the accuracy of the COSX algorithm and also the computational cost are directly determined by the size of the integration grid, so selection of the grid is important. This COSX implementation uses two separate grids. By default, the SCF algorithm is first converged on a smaller grid, followed by a number of SCF iterations up to a maximum value (controlled by the COSX_MAXITER_FINAL keyword) on a larger grid. By default, COSX_MAXITER_FINAL is set to 1, a single SCF iteration, which results in numerical errors comparable to performing the entire SCF on the expensive larger grid at a computational cost much closer to the smaller grid. Setting COSX_MAXITER_FINAL to 0 disables the larger grid entirely. Setting COSX_MAXITER_FINAL to -1 allows for the SCF to fully converge on the larger grid, useful for the study of wavefunction properties such as gradients. The size of the initial grid is controlled by the keywords COSX_RADIAL_POINTS_INITIAL and COSX_SPHERICAL_POINTS_INITIAL . The final grid is controlled by COSX_RADIAL_POINTS_FINAL and COSX_SPHERICAL_POINTS_FINAL . Currently, the default grids are very crude, allowing for high performance at the cost of accuracy. If high-accuracy calculations are desired with COSX, the grid sizes should be increased.

Screening thresholds over integrals, densities, and basis extents are set with the COSX_INTS_TOLERANCE , COSX_DENSITY_TOLERANCE , and COSX_BASIS_TOLERANCE keywords, respectively. COSX_INTS_TOLERANCE is the most consequential of the three thresholds in both cost and accuracy. This keyword determines screening of negligible one-electron integrals. COSX_DENSITY_TOLERANCE controls the threshold for significant shell pairs in the density matrix. Lastly, COSX_BASIS_TOLERANCE is a cutoff for the value of basis functions at grid points. This keyword is used to determine the radial extent of the each basis shell, and it is the COSX analogue to DFT_BASIS_TOLERANCE .

The INCFOCK keyword (defaults to false ) increases performance by constructing the Fock matrix from differences in the density matrix, which are more amenable to screening. This option is disabled by default because of potential SCF convergence issues, particularly when using diffuse basis functions. The COSX_OVERLAP_FITTING keyword (defaults to true ) reduces numerical integration errors using the method described in [Izsak:2011:144105] and is always recommended.

Linear Exchange¶

Large SCF calculations can benefit from specialized screening procedures that further reduce the scaling of the ERI contribution to the Fock matrix. LinK, the linear-scaling exchange method described in [Ochsenfeld:1998:1663] , is available in Psi4 in conjunction with composite algorithms that build J ( SCF_TYPE set to J_alg+LINK ). LinK achieves linear-scaling by exploiting shell pair sparsity in the density matrix and overlap sparsity between shell pairs. Specifically, LinK exploits the fact that the Exchange term requires only a linear-scaling number of significant elements through reformulating the shell quartet screening process to scale linearly with system size. LinK is most competitive when used with non-diffuse orbital basis sets, since orbital and density overlaps decay slower with diffuse functions. LinK is especially powerful when combined with density-matrix based ERI screening (set SCREENING to DENSITY ) and incremental Fock builds (set INCFOCK to TRUE ), which decrease the number of significant two-electron integrals to calculate.

To control the LinK algorithm, here are the list of options provided.

LINK_INTS_TOLERANCE : The integral screening tolerance used for sparsity-prep in the LinK algorithm. Defaults to the INTS_TOLERANCE option.

Second-order Convergence¶

Second-order convergence takes into account both the gradient and Hessian to take a full Newton step with respect to the orbital parameters. This results in quadratic convergence with respect to density for SCF methods. For cases where normal acceleration methods either fail or take many iterations to converge, second-order can reduce the total time to solution.

Solving second-order (SO) methods exactly would require an inversion of the orbital Hessian (an expensive \(\mathbb^6\) operation); however, these equations are normally solved iteratively where each iteration costs the same as a normal Fock build ( \(\mathbb^4\) ). The overall SOSCF operation is thus broken down into micro- and macroiterations where the microiterations refer to solving the SOSCF equations and macroiterations are the construction of a new Fock matrix based on the orbitals from a SOSCF step.

SOSCF requires that all elements of the gradient to be less than one before the method is valid. To this end, pre-SOSCF SCF iterations use normal gradient-based extrapolation procedures (e.g., DIIS) until the gradient conditions are met. Note that while the total number of macroiterations will be less for SOSCF than gradient-based convergence acceleration, the cost of solving the microiterations typically results in the overall cost being greater for SOSCF than for gradient-based methods. Therefore, SOSCF should only be used if it is difficult to locate a stable minimum.

SOSCF is available for all HF and DFT references with the exception of meta- GGA functionals. To enable, set the option SOSCF to true . Additional options to modify the number of microiterations taken are as follows:

SOSCF_START_CONVERGENCE : when to start SOSCF based on the current density RMS

SOSCF_MAX_ITER : the maximum number of SOSCF microiterations per macroiteration

SOSCF_CONV : the relative convergence tolerance of the SOSCF microiterations

SOSCF_PRINT : option to print the microiterations or not

Stability Analysis¶

SCF algorithms attempt to minimize the gradient of the energy with respect to orbital variation parameters. At convergence, the gradient should be approximately zero given a convergence criterion. Although this is enough to make sure the SCF converged to a stationary point, this is not a sufficient condition for a minimal SCF solution. It may be a saddle point or a maximum.

To ensure that a minimum has been found, the electronic Hessian, i.e. the matrix of second derivatives of the energy with respect to orbital variation parameters, must be computed. If one or more eigenvalues of the electronic Hessian are negative, the SCF solution is not a minimum. In that case, orbital parameters can be varied along the lowest Hessian eigenvector to lower the energy.

Orbital variation parameters are usually constrained. For example, in RHF the spatial parts of the \(\alpha\) and \(\beta\) orbitals are the same. In UHF, the orbital coefficients are usually constrained to be real. A stability analysis can check whether a lower SCF solution exists while respecting the constraints of the original solution; this is an internal instability. If one or more constraints have to be relaxed to reach a lower-energy solution, there is an external instability. In P SI 4 , the only external instability that can be checked at present is the RHF \(\rightarrow\) UHF one.

Currently, two algorithms exist in P SI 4 for stability analysis: the original Direct Inversion and the newly implemented Davidson algorithms. We will first describe options common to both algorithms. To request a stability analysis at the end of the SCF, set the keyword STABILITY_ANALYSIS . Value CHECK only computes the electronic Hessian eigenvalue and checks if an actual SCF minimum has been found, while value FOLLOW rotates the converged orbitals along the lowest eigenvector, then invokes the SCF procedure again to lower the energy. In case the minimization does not succeed or ends up on the same unstable solution, you can tune the scale factor for the orbital rotation through the keyword FOLLOW_STEP_SCALE . The rotation angle is \(\frac<\pi>\mbox < >\cdot\) ( FOLLOW_STEP_SCALE ). The default value of 0.5 usually provides a good guess, and modification is only recommended in difficult cases. The default behavior for the stability code is to stop after trying to reoptimize the orbitals once if the instability still exists. For more attempts, set MAX_ATTEMPTS ; the default value of 1 is recommended. In case the SCF ends up in the same minimum, modification of FOLLOW_STEP_SCALE is recommended over increasing MAX_ATTEMPTS .

Setting the option STABILITY_ANALYSIS to FOLLOW is only avalible for UHF. When using RHF and ROHF instabilities can be checked, but not followed. If you want to attempt to find a lower energy solution you should re-run the calculation with REFERENCE set to UHF .

The main algorithm available in P SI 4 is the Direct Inversion algorithm. It can only work with SCF_TYPE PK , and it explicitly builds the full electronic Hessian matrix before explicitly inverting it. As such, this algorithm is very slow and it should be avoided whenever possible. Direct Inversion is automatically invoked if the newer algorithm is not available.

The Davidson algorithm for stability analysis was implemented recently. Only the lowest eigenvalues of the electronic Hessian are computed, and Hessian-vector products are computed instead of the full Hessian. This algorithm is thus much more efficient than the Direct Inversion, but at present, it is only available for UHF \(\rightarrow\) UHF stability analysis. The capabilities of both algorithms are summarized below: