HF: Hartree–Fock Theory¶
Code author: Justin M. Turney, Robert M. Parrish, and Andrew C. Simmonett
Section author: Robert M. Parrish and Jerome F. Gonthier
Module: Keywords, PSI Variables, LIBSCF_SOLVER, LIBMINTS, LIBFOCK, LIBDIIS
Introduction¶
SelfConsistentField (SCF) theory forms the cornerstone of ab initio quantum chemistry. Here SCF refers both to conventional Hartree–Fock (HF) molecular orbital theory and also to generalized Kohn–Sham Density Functional Theory (KSDFT). PSI4 contains a wholly rewritten SCF code, including many of the most popular spin specializations, several efficient numerical methods for treating Fock Matrix construction, and a brand new KSDFT code featuring many of the most popular DFT functional technologies.
An illustrative example of using the SCF module is as follows:
molecule {
0 3
O
O 1 1.21
}
set {
basis ccpvdz
guess sad
reference uhf
scf_type direct
}
energy('scf')
This will run a UHF computation for triplet molecular oxygen (the ground state) using a Direct algorithm for the Electron Repulsion Integrals (ERI) and starting from a Superposition of Atomic Densities (SAD) guess. DF integrals are automatically used to converge the DFSCF solution before the Direct algorithm is activated. After printing all manner of titles, geometries, sizings, and algorithm choices, the SCF finally reaches the iterations:
Total Energy Delta E RMS [F,P]
@DFUHF iter 0: 149.80032977420572 1.49800e+02 1.48808e01
@DFUHF iter 1: 149.59496320631871 2.05367e01 2.58009e02
@DFUHF iter 2: 149.62349901753706 2.85358e02 6.68980e03 DIIS
@DFUHF iter 3: 149.62639942687878 2.90041e03 2.19285e03 DIIS
@DFUHF iter 4: 149.62689561367233 4.96187e04 5.99497e04 DIIS
@DFUHF iter 5: 149.62694151275420 4.58991e05 1.27338e04 DIIS
@DFUHF iter 6: 149.62694337910040 1.86635e06 1.65616e05 DIIS
@DFUHF iter 7: 149.62694340915198 3.00516e08 2.68990e06 DIIS
@DFUHF iter 8: 149.62694340999315 8.41169e10 2.61249e07 DIIS
DF guess converged.
...
@UHF iter 9: 149.62730705472407 3.63645e04 8.63697e05 DIIS
@UHF iter 10: 149.62730737348096 3.18757e07 1.50223e05 DIIS
@UHF iter 11: 149.62730738537113 1.18902e08 3.80466e06 DIIS
@UHF iter 12: 149.62730738624032 8.69193e10 7.06634e07 DIIS
The first set of iterations are from the DF portion of the computation, the second set uses the exact (but much slower) Direct algorithm. Within the DF portion of the computation, the zerothiteration uses a nonidempotent density matrix obtained from the SAD guess, so the energy is unphysically low. However, the first true iteration is quite close to the final DF energy, highlighting the efficiency of the SAD guess. Pulay’s DIIS procedure is then used to accelerate SCF convergence, with the DF phase reaching convergence in eight true iterations. When used together, SAD and DIIS are usually sufficient to converge the SCF for all but the most difficult systems. Additional convergence techniques are available for more difficult cases, and are detailed below. At this point, the code switches on the requested Direct integrals technology, which requires only four full iterations to reach convergence, starting from the DF guess. This hybrid DF/Direct procedure can significantly accelerate SCF computations requiring exact integrals.
After the iterations are completed, a number of oneelectron properties are printed, and some bookkeeping is performed to set up possible correlated computations. Additional oneelectron properties are available by increasing the PRINT option. Also printed are the occupied and virtual orbital energies, which are useful in elucidating the stability and reactivity of the system.
Theory¶
The objective of Hartree–Fock (HF) Theory is to produce the optimized Molecular Orbitals (MOs) \(\{\psi_i\}\),
Here, \(\{\phi_{\mu}\}\) are the basis functions, which, in PSI4 are contracted Cartesian Gaussian functions often referred to as Atomic Orbitals (AOs). The matrix \(C_{\mu i}\) contains the MO coefficients, which are the constrained variational parameters in Hartree–Fock. The molecular orbitals are used to build the simplest possible antisymmetric wavefunction, a single Slater determinant,
This form for the Hartree–Fock wavefunction is actually entirely equivalent to treating the electron correlation as a mean field repulsion in \(\mathbb{R}^6\) instead of a more complicated effect in \(\mathbb{R}^N\).
Considering the electronic Hamiltonian,
the Hartree–Fock energy is, by Slater’s rules,
Here \(H\) is the AObasis oneelectron potential, encapsulating both electronnuclear attraction and kinetic energy,
\(D\) is the AObasis density matrix, build from the occupied orbital coefficients,
and \(F\) is the Fock matrix, which is the effective onebody potential at the current value of the density,
Here the tensor \((\mu\nu\lambda\sigma)\) is an AO ElectronRepulsion Integral (ERI) in chemists’ notation,
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 AObasis 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 selfconsistent with the Fock matrix they generate.
The formation of the Coulomb matrix \(J\) and the exchange matrix \(K^{\alpha}\) 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 sto3g
energy('scf')
This will run a Restricted Hartree–Fock (RHF) on neutral singlet Helium in
\(D_{2h}\) spatial symmetry with a minimal STO3G
basis, 1.0E6
energy and density convergence criteria (since singlepoint, 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¶
PSI4 implements the most popular spin specializations of Hartree–Fock theory, including:
 Restricted Hartree–Fock (RHF) [Default]
Appropriate only for closedshell 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 openshell 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 “spincontamination” problem.
 Restricted OpenShell Hartree–Fock (ROHF)
Appropriate for openshell systems where spincontamination is problem. Sometimes more difficult to converge, and assumes uniformly positive spin polarization (the alpha and beta doublyoccupied 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 PSI4 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_{2v}\) 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 ccpvdz
}
energy('scf')
Broken Symmetry¶
For certain problems, such diradicals, allowing the spinup and spindown orbitals to differ in closedshell 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 nonphysical 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 PSI4, 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 canonicalorthogonalized 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 PSI4, symmetric orthogonalization is used by default, unless the smallest overlap eigenvalue falls below the usersupplied double option S_TOLERANCE, which defaults to 1E7. 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 1E8.
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 symmetryadapted 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 augccpv5z # This diffuse basis will have
# smallish eigenvalues for even H2O
print 3
}
energy('scf')
Output:
==> PreIterations <==

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.6888063568E05.
Using Canonical Orthogonalization with cutoff of 1.0000000000E04.
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 rootmeansquare change in density matrix values, which must be below the userspecified 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, openshell systems, anions, and systems with diffuse basis sets.
For initial orbital selection, several options are available. These include:
 CORE
Diagonalization of the core Hamiltonian, removing even meanfield electron repulsion. Simple, but often too far from the final solution for larger systems. This is the default for single atoms.
 SAD [Default]
Superposition of Atomic Densities. Builds the initial density as the spinaveraged sum of atomic UHF computations in the current basis. If an openshell system, uniform scaling of the spinaveraged 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 closedshell systems. This is the default for systems of more than one atom.
 SADNO
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].
 GWH
A generalized WolfsbergHelmholtz modification of the core Hamiltonian matrix. Usually less accurate than the core guess: the latter is exact for oneelectron systems, GWH is not; see [Lehtola:2019:1593]).
 HUCKEL
An extended Hückel guess based on onthefly atomic UHF calculations alike SAD, see [Lehtola:2019:1593].
 MODHUCKEL
Like HUCKEL, an extended Hückel guess based on onthefly atomic UHF calculations alike SAD, see [Lehtola:2019:1593]. This variant employs an updated rule for the generalized WolfsbergHelmholz formula from [Ammeter:1978:3686].
 READ
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. SAP
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 finiteelement 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 exchangeonly local density calculations that are but nanohartree away from the complete basis set limit [Lehtola:2020:012516].
 SAPGAU
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 threecenter twoelectron 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 321G or pcseg0 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 321G basis and then casting up to ccpVTZ is shown below:
molecule h2o {
0 1
O
H 1 1.0
H 1 1.0 2 104.5
}
set {
basis ccpvtz
basis_guess 321G
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 PSI4 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 PSI4 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 closelyrelated 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 fullyoptimized 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 nonconverging 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 1e2, respectively.
 SOSCF [Off by Default]
ERI Algorithms¶
The key difficulty in the SCF procedure is treatment of the fourindex ERI contributions to the Fock Matrix. A number of algorithms are available in PSI4 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:
 PK [Default]
An outofcore, presorted algorithm using exact ERIs. Quite fast for a zeroerror algorithm if enough memory is available. Integrals are generated only once, and symmetry is utilized to reduce number of integrals.
 OUT_OF_CORE
An outofcore, 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.
 DIRECT
A threaded, sieved, integraldirect 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 preiterations are turned off), and can obtain significant speedups with negligible error loss if INTS_TOLERANCE is set to 1.0E8 or so.
 DF [Default]
A densityfitted 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.
 CD
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.
PSI4 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 nonhybrid 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:
 DFDIRJ
An integraldirect algorithm constructing the Coulomb term based on [Weigend:2002:4285] The DFDIRJ algorithm combines the benefits of integraldirect SCF approaches with that of densityfitting. Specifically, DFJ utilizes no I/O and displays strong performance with large system size through a combination of effective parallelization and utilization of densityfitting to minimize ERI computational cost. See the IntegralDirect DensityFitted Coulomb Construction section for more information.
Specialized algorithms available to construct the Exchange term within a composite framework are as follows:
 COSX
An algorithm based on the seminumerical “chain of spheres exchange” (COSX) approach described in [Neese:2009:98]. The coulomb term is computed with a direct densityfitting algorithm. The COSX algorithm uses no I/O, scales well with system size, and requires minimal memory, making it ideal for large systems and multicore CPUs. See COSX Exchange for more information.
 LINK
An implementation of the linearscaling “Linear Exchange” (LinK) algorithm described in [Ochsenfeld:1998:1663]. The LINK algorithm provides many of the benefits of integraldirect 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 formallyexact 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
 MEM_DF
A DF algorithm optimized around memory layout and is optimal as long as there is sufficient memory to hold the threeindex DF tensors in memory. This algorithm may be faster for builds that require disk if SSDs are used.
 DISK_DF
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 inmemory and ondisk 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
PSI4‘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 inmemory and ondisk options within MEM_DF or DISK_DF using
the SCF_SUBTYPE option. Using SCF_SUBTYPE = AUTO
, where PSI4
automatically selects the inmemory or ondisk option for MEM_DF/DISK_DF based on memory and molecule, is the default
and recommended option. However, the inmemory or ondisk 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
PSI4.
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.0E12 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 matrixbased screening, DENSITY
.
 SCHWARZ
Uses the CauchySchwarz inequality to calculate an upper bounded value of a shell quartet,
 CSAM
An extension of the Schwarz estimate that also screens over the long range 1/r operator, described in [Thompson:2017:144101].
 DENSITY
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 densitymatrix 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 ccpVDZJKFIT as the
auxiliary basis, unless the user specifies DF_BASIS_SCF manually. For
some atoms, ccpVDZJKFIT 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.
IntegralDirect DensityFitted Coulomb Construction¶
The Resolution of the Identity (RI) can be used to decompose the normally 4center ERI tensor into a combination of 3center and 2center components.
By reducing the dimensionality of the ERI tensor, application of the RI (often referred to as densityfitting, 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 integraldirect 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 nonhybrid functionals).
DFDIRJ supports multiple capabilities to improve performance. Specifically, DFDIRJ allows for a combination of densitymatrix 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 seminumerical COSX algorithm described in [Neese:2009:98] evaluates twoelectron 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 highaccuracy 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 oneelectron 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 linearscaling 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 linearscaling 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 linearscaling number of significant elements through reformulating the
shell quartet screening process to scale linearly with system size.
LinK is most competitive when used with nondiffuse orbital basis sets, since orbital and density overlaps decay slower with diffuse functions.
LinK is especially powerful when combined with densitymatrix based ERI screening (set SCREENING to DENSITY
) and incremental Fock builds (set INCFOCK to TRUE
), which decrease the number of significant twoelectron integrals to calculate.
To control the LinK algorithm, here are the list of options provided.
LINK_INTS_TOLERANCE: The integral screening tolerance used for sparsityprep in the LinK algorithm. Defaults to the INTS_TOLERANCE option.
Secondorder Convergence¶
Secondorder 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, secondorder can reduce the total time to solution.
Solving secondorder (SO) methods exactly would require an inversion of the orbital Hessian (an expensive \(\mathbb{N}^6\) operation); however, these equations are normally solved iteratively where each iteration costs the same as a normal Fock build (\(\mathbb{N}^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, preSOSCF SCF iterations use normal gradientbased 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 gradientbased convergence acceleration, the cost of solving the microiterations typically results in the overall cost being greater for SOSCF than for gradientbased 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 lowerenergy solution, there is an external instability. In PSI4, the only external instability that can be checked at present is the RHF \(\rightarrow\) UHF one.
Currently, two algorithms exist in PSI4 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}{2}\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.
Note
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 rerun the calculation with REFERENCE set to UHF
.
The main algorithm available in PSI4 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 Hessianvector 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:
Algorithm 
Stability checked 


Direct Inversion 
RHF 
Internal, External (\(\rightarrow\) UHF) 

PK only 
ROHF 
Internal 

PK only 

Davidson 
UHF 
Internal 

Anything 
The best algorithm is automatically selected, i.e. Davidson for UHF \(\rightarrow\) UHF and Direct Inversion otherwise.
In addition to the options available for Direct Inversion, the Davidson algorithm can automatically adapt FOLLOW_STEP_SCALE to find a new SCF minimum. If MAX_ATTEMPTS > 1, additional attempts will automatically increment FOLLOW_STEP_SCALE by 0.2 every time the SCF falls back to the previously found unstable minimum. The increment can be adjusted by setting FOLLOW_STEP_INCREMENT. The default value is 0.2; adjust if needed to try different values of FOLLOW_STEP_SCALE in a single computation.
The Davidson solver for the eigenvalues is controlled through several keywords.
SOLVER_MAXITER: maximum number of iterations
SOLVER_CONVERGENCE: eigenvector convergence threshold
SOLVER_N_ROOT: Solve for N eigenvectors in each irreducible representation
SOLVER_ROOTS_PER_IRREP: The number of eigenvectors to solve in each irreducible representation. An array of as many integers as there are irreducible representations.
SOLVER_N_GUESS: Use N guess vectors, this needs to be larger than the number of roots so that the lowest ones can be captured reliably. Defaults to 4 guess vectors per root.
Warning
Prior to Dec 2022, v1.7, Psi4 had a different set of keywords controlling instability analysis, and those were included in the CPHF module rather than the SCF module.
Warning
Extending Davidson instability analysis to KohnSham references is under development. As of 1.7, only LDA functionals are currently supported.
In case convergence problems are encountered during the Davidson procedure, file a bug report. If the solver seems to converge on the wrong eigenvalue, try increasing SOLVER_N_GUESS. Otherwise, if the solver is almost converged but reaches the maximum number of iterations, try increasing SOLVER_MAXITER.
Effective core potentials (ECPs)¶
PSI4 supports the use of effective core potentials to describe the innermost electrons in heavy elements. ECPs are only available if PSI4 is compiled with the LibECPInt library. If a basis set is designed to use an effective core potential, the ECP definition should be simply placed alongside the orbital basis set definition, c.f. Inputting effective core potentials (ECPs). All information related to the definition and number of core electrons will automatically be detected and no further input is required to use the ECPcontaining basis set. See scfecp and dfmp2ecp for examples of computations with ECPcontaining basis sets.
Warning
Prior to May 2022, v1.6, Psi4 used a builtin ECP code. Analytic derivatives of ECPs were not available. The HF and DFT derivatives were implemented in a seminumerical scheme, where numerical ECP gradients were added to analytic SCF gradients. For postSCF methods, the entire gradient computation needed to be run as finite difference of energies.
Warning
As of May 2022, v1.6, Psi4 uses the LibECPInt library, and analytic derivatives and Hessians of ECPs are available. Analytic derivatives of molecular systems including ECPs should be available whenever the method has analytic derivatives, but these have so far only been verified for HF and DFT.
Warning
ECPs have not been tested with projected basis set guesses or with FISAPT calculations. If you require this functionality, please contact the developers on GitHub and/or the forum.
External potentials and QM/MM¶
In addition to the implementation of EFP for accurate QM/MM computations, PSI4 can perform more rudimentary QM/MM procedures via the EXTERN keyword. The following snippet, extracted from the extern1 test case, demonstrates its use for a TIP3P external potential:
import numpy as np
external_potentials = [
[0.834, np.array([1.649232019048,0.0,2.356023604706]) / psi_bohr2angstroms],
[ 0.417, np.array([0.544757019107,0.0,3.799961446760]) / psi_bohr2angstroms],
[ 0.417, np.array([0.544757019107,0.0,0.912085762652]) / psi_bohr2angstroms]]
gradient('scf', external_potentials=external_potentials)
The external_potentials
array has three rows for three separate
particles, and it is passed to the SCF code on the last line. The
rows are composed of the atomic charge, x coordinate, y coordinate,
and z coordinate in that order. The atomic charge and coordinates are
specified in atomic units, [e] and [a0]. Add as many particle rows as
needed to describe the full MM region.
Caution
In PSI4 previous to Spring 2022 and v1.6, setting an external potential like the above looked like
Chrgfield = QMMM()
Chrgfield.extern.addCharge(0.834, 1.649232019048, 0.0, 2.356023604706)
Chrgfield.extern.addCharge( 0.417, 0.544757019107, 0.0, 3.799961446760)
Chrgfield.extern.addCharge( 0.417, 0.544757019107, 0.0, 0.912085762652)
psi4.set_global_option_python('EXTERN', Chrgfield.extern)
gradient('scf')
The main differences are that (1) the specification of
charge locations in the old way used the units of the active
molecule, whereas the new way always uses Bohr and (2) the
specification of the charge and locations in the old way used the
psi4.driver.QMMM
class directly and added one charge
per command, whereas the new way consolidates all into an array and
passes it by keyword argument to the calculation.
The successor to the psi4.driver.QMMM
class,
psi4.driver.QMMMbohr
, is operable, but it is discouraged
from being used directly.
To run a computation in a constant dipole field, the PERTURB_H, PERTURB_WITH and PERTURB_DIPOLE keywords can be used. As an example, to add a dipole field of magnitude 0.05 a.u. in the y direction and 0.1 a.u. in the z direction, we can use the following keywords:
set perturb_h true
set perturb_with dipole
set perturb_dipole [ 0, 0.05, 0.1 ]
Note that if any specified fields do not fall along a symmetry axis, the symmetry of the calculation should be reduced accordingly; if in doubt run the calculation in C1 symmetry. For examples of SCF and MP2 calculations in an external field, see scf7 and dfmp2grad5.
Convergence and Algorithm Defaults¶
Ab Initio Method 
Calculation Type 


SCF of HF or DFT 
energy 
6 
6 
DF 
optimization 
8 
8 

frequency [7] 
8 
8 

SCF of postHF 
energy 
8 
8 
PK [3] 
optimization 
10 
10 

frequency [7] 
10 
10 

CC property [2] 
10 
10 
Ab Initio Method 
Calculation Type 
E_CONVERGENCE [5] 
R_CONVERGENCE [6] 

postHF of postHF 
energy 
6 

optimization 
8 

frequency [7] 
8 

CC property [2] 
8 
Footnotes
Recommendations¶
The SCF code is quite flexible and powerful. We have tried as much as possible to keep the number of options to a minimum and to allow all options to be used in the presence of all other options. Below are some rough words of advice about using the SCF code for practical calculations:
For GUESS, the
SAD
guess is usually your friend, even for openshell systems (at the very least, it gets the right number of electrons, unlike some other programs). For instance, we have found that a simple SAD guess is often as good as doing a full SCF in a 321G basis and then performing a castup, at a fraction of the cost. However, SAD and DOCC/SOCC arrays do not play very well together at the moment.For wall time,
DF
may be a factor of ten or more faster than the exact integral technologies available in PSI4. UseDF
unless you need absolute accuracy or do not have a JKFIT auxiliary set for your orbital basis/atom type. Then useDIRECT
.Don’t mess with the DIIS convergence options unless convergence is a problem. We have optimized the parameters for efficiency over a wide array of system types.
Buy a developer a beer!
The “bestpractice” input file for HF is:
memory 1 GB # As much as you've got, the DF algorithm can use
molecule {
O
H 1 1.0
H 1 1.0 2 104.5
}
set {
basis ccpvdz
scf_type df
guess sad
ints_tolerance 1.0E10 # Even this is epically tight, 1.0E8 is OK
}
energy('scf')