Modules: pyscf.scf, pyscf.pbc.scf, pyscf.soscf
Introduction¶
Self-consistent field (SCF) methods include both Hartree-Fock (HF) theoryand Kohn-Sham (KS) density functional theory (DFT). Self-consistentfield theories only depend on the electronic density matrices, and arethe simplest level of quantum chemical models. Details that arespecific to DFT can be found in Density functional theory (DFT).
In both HF and KS-DFT, the ground-state wavefunction is expressed as asingle Slater determinant \(\Phi_0\) of molecular orbitals (MOs)\(\psi\), \(\Phi_0 = \mathcal{A}|\psi_1(1)\psi_2(2) \ldots\psi_N(N)|\). The total electronic energy\(E=\langle\Psi_0|\hat{H}|\Psi_0\rangle\) is then minimized,subject to orbital orthogonality; this is equivalent to thedescription of the electrons as independent particles that onlyinteract via each others’ mean field.
It can be shown that the minimization of the total energy within agiven basis set (see e.g. [1] or any standardtextbook on quantum chemistry like [2]) leads to theequation
\[\mathbf{F} \mathbf{C} = \mathbf{S} \mathbf{C} \mathbf{E}\]
where \(\mathbf{C}\) is the matrix of molecular orbitalcoefficients, \(\mathbf{E}\) is a diagonal matrix of thecorresponding eigenenergies, and \(\mathbf{S}\) is the atomicorbital overlap matrix. The Fock matrix \(\mathbf{F}\) is definedas
\[\mathbf{F} = \mathbf{T} + \mathbf{V} + \mathbf{J} + \mathbf{K}\;\]
where \(\mathbf{T}\) is the kinetic energy matrix,\(\mathbf{V}\) is the external potential, \(\mathbf{J}\) isthe Coulomb matrix, and \(\mathbf{K}\) is the exchange matrix.
Initial guess¶
The Coulomb and exchange matrices depend on the occupied orbitals,meaning that the SCF equation \(\mathbf{F C}=\mathbf{S C E}\)needs to be solved self-consistently by some sort of iterativeprocedure, which begins from some initial guess. The accuracy ofseveral initial guesses for SCF calculations has recently beenassessed in [3], to which we refer fordetailed discussion on initial guesses.
Several of initial guess have been implemented in PySCF; the usedvariant is controlled by the init_guess attribute of the SCFsolver. The following values are possible
'minao'(default)A superposition of atomic densities[4, 5] technique, inwhich the guess is obtained by projecting the minimal basis of thefirst contracted functions in the cc-pVTZ or cc-pVTZ-PP basis setonto the orbital basis set, and then forming the densitymatrix. The guess orbitals are obtained by diagonalizing the Fockmatrix that arises from the spin-restricted guess density.
'1e'The one-electron guess, also known as the core guess, obtains theguess orbitals from the diagonalization of the core Hamiltonian\(\mathbf{H}_0 = \mathbf{T} + \mathbf{V}\), thereby ignoringall interelectronic interactions and the screening of nuclearcharge. The 1e guess should only be used as a last resort, becauseit is so bad for molecular systems; see[3].
'atom'Superposition of atomic densities[4, 5]. Employsspin-restricted atomic HF calculations that employ sphericallyaveraged fractional occupations with ground states determined withfully numerical calculations at the complete basis set limit in[6].
'huckel'This is the parameter-free Hückel guess described in[3], which is based on on-the-fly atomicHF calculations that are performed analogously to
'atom'. Thespherically averaged atomic spin-restricted Hartree-Fockcalculations yield a minimal basis of atomic orbitals and orbitalenergies, which are used to build a Hückel type matrix that isdiagonalized to obtain guess orbitals.'vsap'Superposition of atomic potentials as described in[3]. A sum of pretabulated, fullynumerical atomic potentials determined with the approach of[6] is used to build a guess potentialon a DFT quadrature grid; this potential is then used to obtainthe orbitals. Note this option is only available for DFTcalculations in PySCF.
'chk'Read in the orbitals from the checkpoint file and use them as theinitial guess (see below for more details).
Alternatively, the user can also override the initial guess densitymatrix for an SCF calculation through the dm0 argument. Forexample, the following script first computes the HF density matrix forthe \(\rm Cr^{6+}\) cation, and then uses it as an initial guessfor a HF calculation of the \(\rm Cr\) atom.
# First calculate the Cr6+ cationmol = gto.Mole()mol.build( symmetry = 'D2h', atom = [['Cr',(0, 0, 0)], ], basis = 'cc-pvdz', charge = 6, spin = 0,)mf = scf.RHF(mol)mf.kernel()dm1 = mf.make_rdm1()# Now switch to the neutral atom in the septet statemol.charge = 0mol.spin = 6mol.build(False,False)mf = scf.RHF(mol)mf.kernel(dm0=dm1)
More examples can be found inexamples/scf/15-initial_guess.py, andexamples/scf/31-cr_atom_rohf_tune_init_guess.py.
Restart from an old calculation¶
Although alike many other quantum chemistry codes, there is norestart mechanism available in PySCF package, calculations can stillbe “restarted” by reading in an earlier wave function as the initialguess for the wave function. The initial guess can be prepared inmany ways. One is to read the chkpoint file which is generated inthe previous or other calculations:
>>> from pyscf import scf>>> mf = scf.RHF(mol)>>> mf.chkfile = '/path/to/chkfile'>>> mf.init_guess = 'chkfile'>>> mf.kernel()
/path/to/chkfile can be found in the output in the calculation (ifmol.verbose >= 4, the filename of the chkfile will be dumped inthe output). If the results of the calculation are needed at a laterstage (e.g. for an eventual restart or use as an initial guess for alarger calculation), the chkfile attribute should be setexplicitly as the chkfile might otherwise be deleted upon successfulcompletion of the calculation, see comments inexamples/scf/14-restart.py.By setting chkfile and init_guess, theSCF module can read the molecular orbitals from the givenchkfile and rotate them to representation of the requiredbasis.
The initial guess can also be fed in directly to the calculation. Forexample, we can read in the density matrix from a checkpoint file, andpass it directly to the SCF solver with:
>>> from pyscf import scf>>> mf = scf.RHF(mol)>>> dm = scf.hf.from_chk(mol, '/path/to/chkfile')>>> mf.kernel(dm)
This approach leads to the same result as setting init_guessto chkfile.
N.B. The chkfile initial guess is not limited to calculations onthe same molecule or the same basis set. One can first do a cheaperSCF calculation with smaller basis sets, or run an SCF calculation ona model system (e.g. drop a few atoms or run the same system in aneasier charge/spin state), then use scf.hf.from_chk() to projectthe results to the target basis sets.
Converging SCF iterations¶
Even with a very good initial guess, making the SCF procedure convergeis sometimes challenging. PySCF implements two kinds of approaches forSCF, namely, direct inversion in the iterative subspace (DIIS) andsecond-order SCF (SOSCF).
DIIS (default)
With DIIS, the Fock matrix at each iteration is extrapolated usingFock matrices from the previous iterations, by minimizing the normof the commutator \([\mathbf{F},\mathbf{PS}]\) where\(\mathbf{P}\) is the density matrix[7, 8]. Two variants of DIISare implemented in PySCF, namely, EDIIS [9]and ADIIS [10]. Examples of selectingdifferent DIIS schemes can be found inexamples/scf/24-tune_diis.py.
SOSCF
To achieve quadratic convergence in the orbital optimization,PySCF implements a general second-order solver called theco-iterative augmented hessian (CIAH) method[11, 12]. This method can be invoked bydecorating the SCF objects with the newton() method:
mf = scf.RHF(mol).newton()
More examples can be found inexamples/scf/22-newton.py.
Damping
The Fock matrix can be damped before DIIS acceleration kicks in.This is achieved by setting the attributes damp anddiis_start_cycle. For example,
mf.damp = 0.5mf.diis_start_cycle = 2
means that DIIS will start at the second cycle, and that the Fockmatrix is damped by 50% in the first cycle.
Level shifting
A level shift increases the gap between the occupied and virtualorbitals, thereby slowing down and stabilizing the orbital update.A level shift can help to converge SCF in the case of systems withsmall hom*o-LUMO gaps. Level shifting is invoked by setting theattribute level_shift. See examples inexamples/scf/03-level_shift.py, andexamples/scf/52-dynamically_control_level_shift.py.
See AlsoTrade finance – what is it?Fractional occupations
Fractional occupations can also be invoked to help the SCFconverge for small gap systems. See the example inexamples/scf/54-fractional_occupancy.py.
Smearing
Smearing sets fractional occupancies according to a temperaturefunction. See the example examples/pbc/23-smearing.py.
Stability analysis¶
Even when the SCF converges, the wave function that is found may notcorrespond to a local minimum; calculations can sometimes alsoconverge onto saddle points. Since saddle points are also extrema ofthe energy functional, the orbital gradient vanishes and the SCFequation \(\mathbf{F C}=\mathbf{S C E}\) is satisfied[1]. However, in such cases the energy can bedecreased by perturbing the orbitals away from the saddle point, whichmeans that the wave function is unstable.
Instabilities in the wave function are conventionally classified aseither internal or external [13]. Externalinstabilities mean that the energy can be decreased by loosening someconstraints on the wave function, such as allowing restrictedHartree-Fock orbitals to transform into unrestricted Hartree-Fock,whereas internal instabilities mean that the SCF has converged onto anexcited state instead of the ground state. PySCF allows detecting bothinternal and external instabilities for a given SCF calculation; seethe examples in examples/scf/17-stability.py.
Property calculations¶
Various properties can be computed by calling the correspondingfunctions, for example,
dipole moment:
mf.dip_moment()
Mulliken population:
mf.mulliken_pop()
nuclear gradients:
g = mf.Gradients()g.kernel()
Also several response properties are available in PySCFwith the properties extension, see theexamples there.
Spin-restricted, spin-unrestricted, restricted open-shell, and generalized calculations¶
The general spin-orbital used in the HF or KS-DFT wave function can bewritten as
\[\psi_i(1) = \phi_{i\alpha}(r)\alpha + \phi_{i\beta}(r)\beta \;,\]
Four variants of the ansatz \(\psi(1)\) are commonly used inquantum chemistry; they are also all available in PySCF.
Restricted (RHF/RKS)
The spin-orbitals are either alpha (spin-up) or beta (spin-down),\(\psi_i =\phi_i(r)\alpha\) or \(\psi_i = \phi_i(r)\beta\),and the alpha and beta orbitals share the same spatial orbital\(\phi_i(r)\). The closed-shell determinant is thus\(\Phi=\mathcal{A}|\phi_1(r_1)\alpha \phi_1(r_2)\beta \ldots\phi_{N/2}(r_{N-1})\alpha \phi_{N/2}(r_N)\beta|\) and \(S=0\).
Unrestricted (UHF/UKS)
The orbitals can have either alpha or beta spin, but the alpha andbeta orbitals may have different spatial components. The determinantis thus \(\Phi=\mathcal{A}|\phi_1(r_1)\sigma_1\phi_2(r_2)\sigma_2 \ldots \phi_{N}(r_N)\sigma_N|\) where\(\sigma \in \{\alpha,\beta\}\). Spin contamination isintroduced for states that don’t have maximal \(S_z\).
Restricted open-shell (ROHF/ROKS)
Equivalent to RHF/RKS for \(N_\alpha = N_\beta\). For\(N_\alpha > N_\beta\), the first \(N_\beta\) orbitals havethe same spatial components for both \(\alpha\) and\(\beta\) spin. The remaining \(N_\alpha - N_\beta\)orbitals are of \(\alpha\) spin. \(\Phi=\mathcal{A}|\phi_1\alpha \phi_1\beta \ldots \phi_{N_\beta} \alpha \phi_{N_\beta}\beta\phi_{N_\beta+1}\alpha \ldots \phi_{N}\alpha|\) The finalwavefunction is an eigenfunction of the \(\hat{S}^2\) operatorwith \(S_z=S\).
Generalized (GHF/GKS)
The general form of the spin-orbital \(\psi\) is used. GHF/GKSis useful when none of the previous methods provide stable solutions(see examples/scf/17-stability.py), or when theHamiltonian does not commute with \(\hat{S}_z\) (e.g. in thepresence of spin-orbit coupling, seeexamples/scf/44-soc_ecp.py).
Calculations with these methods can be invoked by creating an instanceof the corresponding class:
mf = scf.RHF(mol).run()mf = scf.UHF(mol).run()mf = scf.ROHF(mol).run()mf = scf.GHF(mol).run()mf = scf.RKS(mol).run()mf = scf.UKS(mol).run()mf = scf.ROKS(mol).run()mf = scf.GKS(mol).run()
More examples can be found inexamples/scf/00-simple_hf.py,examples/scf/01-h2o.py,examples/scf/02-rohf_uhf.py, andexamples/scf/02-ghf.py.
Linear dependencies¶
Most quantum chemistry programs solve the self-consistent fieldequations
\[\mathbf{F} \mathbf{C} = \mathbf{S} \mathbf{C} \mathbf{E}\]
in an orthonormal basis, which is formally obtained as
\[\mathbf{C} = \mathbf{X} \tilde{\mathbf{C}}\]
where the orthogonalizing matrix \(\mathbf{X}\) is typicallychosen as \(\mathbf{X}=\mathbf{S}^{-1/2}\). By expressing theorbitals in terms of the half-inverse overlap matrix, the generalizedeigenproblem \(\mathbf{F} \mathbf{C} = \mathbf{S} \mathbf{C}\mathbf{E}\) can be rewritten as a regular eigenproblem\(\tilde{\mathbf{F}} \tilde{\mathbf{C}} = \tilde{\mathbf{C}}\mathbf{E}\) [1].
Moreover, as the half-inverse overlap matrix \(\mathbf{S}^{-1/2}\)is typically formed by the canonical orthonormalization procedure[14] in which eigenvectors of the overlap matrixwith small eigenvalues are thrown out, this procedure typicallyresults in a better-conditioned basis since linearly dependentcombinations of the basis functions are excluded by the procedure.
At variance, PySCF relies on SciPy’s generalized eigenvalue solver bydefault, which may fail for poorly conditioned basis sets. One can,however, switch to the use of canonical orthonormalization by togglinge.g.:
mf = scf.RHF(mol).apply(scf.addons.remove_linear_dep_)
In the presence of truly pathological linear dependencies, such asthose that occur in molecular calculations with multiply augmentedbasis sets, and at extreme molecular geometries where two nuclei areclose to each other, also canonical orthonormalization fails.However, the addons module also implements the partial Choleskyorthonormalization technique[15, 16], which has beenshown to work reliably even in the presence of such truly pathologicallinear dependencies.
Scalar relativistic correction¶
Scalar relativistic effects can be applied on the one-body operatorsthrough spin-free eXact-2-component (SFX2C) Hamiltonian[17]. The SFX2C Hamiltonian can be invoked bydecorating the SCF objects with the x2c() method, three otherequivalent function names are also listed below:
mf = scf.RHF(mol).x2c()mf = scf.RHF(mol).x2c1e()mf = scf.RHF(mol).sfx2c()mf = scf.RHF(mol).sfx2c1e()
Note that the SFX2C Hamiltonian only changes the one-body operators,and it only accounts for the mass-velocity effect, while picturechange effect and spin-orbit coupling are not included. Once the SCFobject is decorated by x2c() method, the corresponding post-SCFobjects will also automatically have the SFX2C Hamiltonian applied.To turn it off explicitly, one can do:
mf.with_x2c = False
More examples can be found inexamples/scf/21-x2c.py.
