6. scf — Mean-field methods

6.1. Stability analysis

6.2. Addons

Special treatments may be required to the SCF methods in some situations. These special treatments cannot be universally applied for all SCF models. They were defined in the scf.addons module. For example, in an UHF calculation, we may want the \(S_z\) value to be changed (the numbers of alpha and beta electrons not conserved) during SCF iteration while conserving the total number of electrons. scf.addons.dynamic_sz_() can provide this functionality:

from pyscf import gto, scf
mol = gto.M(atom='O 0 0 0; O 0 0 1')
mf = scf.UHF(mol)
mf.verbose=4
mf = scf.addons.dynamic_sz_(mf)
mf.kernel()
print('S^2 = %s, 2S+1 = %s' % mf.spin_square())

This function automatically converges the ground sate of oxygen molecule to triplet state although we didn’t specify spin state in the mol object.

Note

Function scf.addons.dynamic_sz_() has side effects. It changes the underlying mean-field object.

The addons mechanism increases the flexibility of PySCf program. You can define various addons to customize the default behaviour of pyscf program. For example, if you’d like to track the changes of the density (the diagonal term of density matrix) of certain basis during the SCF iteration, you can write the following addon to output the required density:

def output_density(mf, basis_label):
    ao_labels = mf.mol.ao_labels()
    old_make_rdm1 = mf.make_rdm1
    def make_rdm1(mo_coeff, mo_occ):
        dm = old_make_rdm1(mo_coeff, mo_occ)
        print('AO         alpha             beta')
        for i,s in enumerate(ao_labels):
            if basis_label in s:
                print(s, dm[0][i,i], dm[1][i,i])
        return dm
    mf.make_rdm1 = make_rdm1
    return mf
from pyscf import gto, scf
mol = gto.M(atom='O 0 0 0; O 0 0 1')
mf = scf.UHF(mol)
mf.verbose=4
mf = scf.addons.dynamic_sz_(mf)
mf = output_density(mf, 'O 2p')
mf.kernel()

6.3. Caching two-electron integrals

When memory is enough (specified by the max_memory of SCF object), the SCF object generates all two-electron integrals in memory and cache them in _eri array. The default max_memory (defined in lib.parameters.MAX_MEMORY, see Maximum memory) is 4 GB. It roughly corresponds to two-electron real integrals for 250 orbitals. For small systems, the cached integrals usually provide the best performance. If you have enough main memory in your computer, you can increase the max_memory of SCF object to cache the integrals in memory.

The cached integrals _eri are treated as a dense tensor. When system becomes larger and the two-electron integral tensor becomes sparse, caching integrals may lose performance advantage. This is mainly due to the fact that the implementation of J/K build for the cached integrals did not utilize the sparsity of the integral tensor. Also, the data locality was not considered in the implementation which sometimes leads to bad OpenMP multi-threading speed up. For large system, the AO-driven direct SCF method is more favorable.

6.4. Customizing Hamiltonian

This integral object _eri is not merely used by the mean-field calculation. Along with the get_hcore() method, this two-electron integral object will be treated as the Hamiltonian in the post-SCF code whenever possible. This mechanism provides a way to model arbitrary fermion system in PySCF. You can customize a system by changing the 1-electron Hamiltonian and the mean-field _eri attribute. For example, the following code solves a model system:

import numpy
from pyscf import gto, scf, ao2mo, ccsd
mol = gto.M()
n = 10
mol.nelectron = 10
mf = scf.RHF(mol)

t = numpy.zeros((n,n))
for i in range(n-1):
    t[i,i+1] = t[i+1,i] = -1.0
t[n-1,0] = t[0,n-1] = 1.0  # anti-PBC
eri = numpy.zeros((n,n,n,n))
for i in range(n):
    eri[i,i,i,i] = 4.0

mf.get_hcore = lambda *args: t
mf.get_ovlp = lambda *args: numpy.eye(n)
# ao2mo.restore(8, eri, n) to get 8-fold symmetry of the integrals
# ._eri only supports the 2-electron integrals in 4-fold or 8-fold symmetry.
mf._eri = ao2mo.restore(8, eri, n)

mf.kernel()

mycc = ccsd.RCCSD(mf).run()
e,v = mycc.ipccsd(nroots=3)
print('IP = ', e)
e,v = mycc.eaccsd(nroots=3)
print('EA = ', e)

Some post-SCF methods require the 4-index MO integrals. Depending the available memory (affected by the value of max_memory in each class), these methods may not use the “AO integrals” cached in _eri. To ensure the post mean-field methods to use the _eri integrals no matter whether the actual memory usage is over the max_memory limite, you can set the flag incore_anyway in Mole class to True before calling the kernel() function of the post-SCF methods. In the following example, without setting incore_anyway=True, the CCSD calculations may crash:

import numpy
from pyscf import gto, scf, ao2mo, ccsd
mol = gto.M()
n = 10
mol.nelectron = n
mol.max_memory = 0
mf = scf.RHF(mol)

t = numpy.zeros((n,n))
for i in range(n-1):
    t[i,i+1] = t[i+1,i] = 1.0
t[n-1,0] = t[0,n-1] = -1.0
eri = numpy.zeros((n,n,n,n))
for i in range(n):
    eri[i,i,i,i] = 4.0

mf.get_hcore = lambda *args: t
mf.get_ovlp = lambda *args: numpy.eye(n)
mf._eri = ao2mo.restore(8, eri, n)
mf.kernel()

mol.incore_anyway = True
mycc = ccsd.RCCSD(mf).run()
e,v = mycc.ipccsd(nroots=3)
print('IP = ', e)
e,v = mycc.eaccsd(nroots=3)
print('EA = ', e)

Holding the entire two-particle interactions matrix elements in memory often leads to high memory usage. In the SCF calculation, the memory usage can be optimized if _eri is sparse. The SCF iterations requires only the Fock matrix which in turn calls the J/K build function SCF.get_jk() to compute the Coulomb and HF-exchange matrix. Overwriting the SCF.get_jk() function can reduce the memory footprint of the SCF part in the above example:

import numpy
from pyscf import gto, scf, ao2mo, ccsd
mol = gto.M()
n = 10
mol.nelectron = n
mol.max_memory = 0
mf = scf.RHF(mol)

t = numpy.zeros((n,n))
for i in range(n-1):
    t[i,i+1] = t[i+1,i] = 1.0
t[n-1,0] = t[0,n-1] = -1.0

mf.get_hcore = lambda *args: t
mf.get_ovlp = lambda *args: numpy.eye(n)
def get_jk(mol, dm, *args):
    j = numpy.diag(dm.diagonal()) * 4.
    k = numpy.diag(dm.diagonal()) * 4.
    return j, k
mf.get_jk = get_jk
mf.kernel()

Another way to handle the two-particle interactions of large model system is to use the density fitting/Cholesky decomposed integrals. See also Saving/Loading DF integral tensor.

6.5. Program reference

6.5.1. Hartree-Fock

Simple usage:

>>> from pyscf import gto, scf
>>> mol = gto.M(atom='H 0 0 0; H 0 0 1')
>>> mf = scf.RHF(mol).run()

scf.RHF() returns an instance of SCF class. There are some parameters to control the SCF method.

verbose
: int
Print level. Default value equals to Mole.verbose
max_memory
: float or int
Allowed memory in MB. Default value equals to Mole.max_memory
chkfile
: str
checkpoint file to save MOs, orbital energies etc.
conv_tol
: float
converge threshold. Default is 1e-10
max_cycle
: int
max number of iterations. Default is 50
init_guess
: str
initial guess method. It can be one of ‘minao’, ‘atom’, ‘1e’, ‘chkfile’. Default is ‘minao’
DIIS
: class listed in scf.diis
Default is diis.SCF_DIIS. Set it to None/False to turn off DIIS.
diis
: bool
whether to do DIIS. Default is True.
diis_space
: int
DIIS space size. By default, 8 Fock matrices and errors vector are stored.
diis_start_cycle
: int
The step to start DIIS. Default is 0.
level_shift_factor
: float or int
Level shift (in AU) for virtual space. Default is 0.
direct_scf
: bool
Direct SCF is used by default.
direct_scf_tol
: float
Direct SCF cutoff threshold. Default is 1e-13.
callback
: function
callback function takes one dict as the argument which is generated by the builtin function locals(), so that the callback function can access all local variables in the current envrionment.
conv_check
: bool
An extra cycle to check convergence after SCF iterations.
nelec
: (int,int), for UHF/ROHF class
freeze the number of (alpha,beta) electrons.
irrep_nelec
: dict, for symmetry- RHF/ROHF/UHF class only
to indicate the number of electrons for each irreps. In RHF, give {‘ir_name’:int, ...} ; In ROHF/UHF, give {‘ir_name’:(int,int), ...} . It is effective when Mole.symmetry is set True.
auxbasis
: str, for density fitting SCF only

Auxiliary basis for density fitting.

>>> mf = scf.density_fit(scf.UHF(mol))
>>> mf.scf()

Density fitting can be applied to all non-relativistic HF class.

with_ssss
: bool, for Dirac-Hartree-Fock only
If False, ignore small component integrals (SS|SS). Default is True.
with_gaunt
: bool, for Dirac-Hartree-Fock only
If False, ignore Gaunt interaction. Default is False.

Saved results

converged
: bool
SCF converged or not
e_tot
: float
Total HF energy (electronic energy plus nuclear repulsion)
mo_energy
:
Orbital energies
mo_occ
Orbital occupancy
mo_coeff
Orbital coefficients

6.5.2. Non-relativistic Hartree-Fock

class pyscf.scf.hf.SCF(mol)[source]

SCF base class. non-relativistic RHF.

Attributes:
verbose
: int
Print level. Default value equals to Mole.verbose
max_memory
: float or int
Allowed memory in MB. Default equals to Mole.max_memory
chkfile
: str
checkpoint file to save MOs, orbital energies etc. Writing to chkfile can be disabled if this attribute is set to None or False.
conv_tol
: float
converge threshold. Default is 1e-9
conv_tol_grad
: float
gradients converge threshold. Default is sqrt(conv_tol)
max_cycle
: int
max number of iterations. Default is 50
init_guess
: str
initial guess method. It can be one of ‘minao’, ‘atom’, ‘hcore’, ‘1e’, ‘chkfile’. Default is ‘minao’
DIIS
: DIIS class
The class to generate diis object. It can be one of diis.SCF_DIIS, diis.ADIIS, diis.EDIIS.
diis
: boolean or object of DIIS class defined in scf.diis.
Default is the object associated to the attribute self.DIIS. Set it to None/False to turn off DIIS. Note if this attribute is inialized as a DIIS object, the SCF driver will use this object in the iteration. The DIIS informations (vector basis and error vector) will be held inside this object. When kernel function is called again, the old states (vector basis and error vector) will be reused.
diis_space
: int
DIIS space size. By default, 8 Fock matrices and errors vector are stored.
diis_start_cycle
: int
The step to start DIIS. Default is 1.
diis_file: ‘str’
File to store DIIS vectors and error vectors.
level_shift
: float or int
Level shift (in AU) for virtual space. Default is 0.
direct_scf
: bool
Direct SCF is used by default.
direct_scf_tol
: float
Direct SCF cutoff threshold. Default is 1e-13.
callback
: function(envs_dict) => None
callback function takes one dict as the argument which is generated by the builtin function locals(), so that the callback function can access all local variables in the current envrionment.
conv_check
: bool
An extra cycle to check convergence after SCF iterations.
check_convergence
: function(envs) => bool
A hook for overloading convergence criteria in SCF iterations.

Saved results:

converged
: bool
SCF converged or not
e_tot
: float
Total HF energy (electronic energy plus nuclear repulsion)
mo_energy :
Orbital energies
mo_occ
Orbital occupancy
mo_coeff
Orbital coefficients

Examples:

>>> mol = gto.M(atom='H 0 0 0; H 0 0 1.1', basis='cc-pvdz')
>>> mf = scf.hf.SCF(mol)
>>> mf.verbose = 0
>>> mf.level_shift = .4
>>> mf.scf()
-1.0811707843775884
analyze(verbose=None, with_meta_lowdin=True, **kwargs)[source]

Analyze the given SCF object: print orbital energies, occupancies; print orbital coefficients; Mulliken population analysis; Diople moment.

as_scanner(mf)

Generating a scanner/solver for HF PES.

The returned solver is a function. This function requires one argument “mol” as input and returns total HF energy.

The solver will automatically use the results of last calculation as the initial guess of the new calculation. All parameters assigned in the SCF object (DIIS, conv_tol, max_memory etc) are automatically applied in the solver.

Note scanner has side effects. It may change many underlying objects (_scf, with_df, with_x2c, ...) during calculation.

Examples:

>>> from pyscf import gto, scf
>>> hf_scanner = scf.RHF(gto.Mole().set(verbose=0)).as_scanner()
>>> hf_scanner(gto.M(atom='H 0 0 0; F 0 0 1.1'))
-98.552190448277955
>>> hf_scanner(gto.M(atom='H 0 0 0; F 0 0 1.5'))
-98.414750424294368
canonicalize(mf, mo_coeff, mo_occ, fock=None)

Canonicalization diagonalizes the Fock matrix within occupied, open, virtual subspaces separatedly (without change occupancy).

dip_moment(mol=None, dm=None, unit='Debye', verbose=3, **kwargs)[source]

Dipole moment calculation

\[\begin{split}\mu_x = -\sum_{\mu}\sum_{\nu} P_{\mu\nu}(\nu|x|\mu) + \sum_A Q_A X_A\\ \mu_y = -\sum_{\mu}\sum_{\nu} P_{\mu\nu}(\nu|y|\mu) + \sum_A Q_A Y_A\\ \mu_z = -\sum_{\mu}\sum_{\nu} P_{\mu\nu}(\nu|z|\mu) + \sum_A Q_A Z_A\end{split}\]

where \(\mu_x, \mu_y, \mu_z\) are the x, y and z components of dipole moment

Args:
mol: an instance of Mole dm : a 2D ndarrays density matrices
Return:
A list: the dipole moment on x, y and z component
eig(h, s)[source]

Solver for generalized eigenvalue problem

\[HC = SCE\]
energy_elec(mf, dm=None, h1e=None, vhf=None)

Electronic part of Hartree-Fock energy, for given core hamiltonian and HF potential

... math:

E = \sum_{ij}h_{ij} \gamma_{ji}
  + \frac{1}{2}\sum_{ijkl} \gamma_{ji}\gamma_{lk} \langle ik||jl\rangle
Args:
mf : an instance of SCF class
Kwargs:
dm
: 2D ndarray
one-partical density matrix
h1e
: 2D ndarray
Core hamiltonian
vhf
: 2D ndarray
HF potential
Returns:
Hartree-Fock electronic energy and the Coulomb energy

Examples:

>>> from pyscf import gto, scf
>>> mol = gto.M(atom='H 0 0 0; H 0 0 1.1')
>>> mf = scf.RHF(mol)
>>> mf.scf()
>>> dm = mf.make_rdm1()
>>> scf.hf.energy_elec(mf, dm)
(-1.5176090667746334, 0.60917167853723675)
>>> mf.energy_elec(dm)
(-1.5176090667746334, 0.60917167853723675)
energy_tot(mf, dm=None, h1e=None, vhf=None)

Total Hartree-Fock energy, electronic part plus nuclear repulstion See scf.hf.energy_elec() for the electron part

from_chk(chkfile=None, project=None)[source]

Read the HF results from checkpoint file, then project it to the basis defined by mol

Returns:
Density matrix, 2D ndarray
get_fock(mf, h1e=None, s1e=None, vhf=None, dm=None, cycle=-1, diis=None, diis_start_cycle=None, level_shift_factor=None, damp_factor=None)

F = h^{core} + V^{HF}

Special treatment (damping, DIIS, or level shift) will be applied to the Fock matrix if diis and cycle is specified (The two parameters are passed to get_fock function during the SCF iteration)

Kwargs:
h1e
: 2D ndarray
Core hamiltonian
s1e
: 2D ndarray
Overlap matrix, for DIIS
vhf
: 2D ndarray
HF potential matrix
dm
: 2D ndarray
Density matrix, for DIIS
cycle
: int
Then present SCF iteration step, for DIIS
diis
: an object of SCF.DIIS class
DIIS object to hold intermediate Fock and error vectors
diis_start_cycle
: int
The step to start DIIS. Default is 0.
level_shift_factor
: float or int
Level shift (in AU) for virtual space. Default is 0.
get_grad(mo_coeff, mo_occ, fock=None)[source]

RHF orbital gradients

Args:
mo_coeff
: 2D ndarray
Obital coefficients
mo_occ
: 1D ndarray
Orbital occupancy
fock_ao
: 2D ndarray
Fock matrix in AO representation
Returns:
Gradients in MO representation. It’s a num_occ*num_vir vector.
get_j(mol=None, dm=None, hermi=1)[source]

Compute J matrix for the given density matrix.

get_jk(mol=None, dm=None, hermi=1)[source]

Compute J, K matrices for the given density matrix

Args:

mol : an instance of Mole

dm
: ndarray or list of ndarrays
A density matrix or a list of density matrices
Kwargs:
hermi
: int

Whether J, K matrix is hermitian

0 : not hermitian and not symmetric
1 : hermitian or symmetric
2 : anti-hermitian
vhfopt :
A class which holds precomputed quantities to optimize the computation of J, K matrices
Returns:
Depending on the given dm, the function returns one J and one K matrix, or a list of J matrices and a list of K matrices, corresponding to the input density matrices.

Examples:

>>> from pyscf import gto, scf
>>> from pyscf.scf import _vhf
>>> mol = gto.M(atom='H 0 0 0; H 0 0 1.1')
>>> dms = numpy.random.random((3,mol.nao_nr(),mol.nao_nr()))
>>> j, k = scf.hf.get_jk(mol, dms, hermi=0)
>>> print(j.shape)
(3, 2, 2)
get_k(mol=None, dm=None, hermi=1)[source]

Compute K matrix for the given density matrix.

get_occ(mf, mo_energy=None, mo_coeff=None)

Label the occupancies for each orbital

Kwargs:
mo_energy
: 1D ndarray
Obital energies
mo_coeff
: 2D ndarray
Obital coefficients

Examples:

>>> from pyscf import gto, scf
>>> mol = gto.M(atom='H 0 0 0; F 0 0 1.1')
>>> mf = scf.hf.SCF(mol)
>>> energy = numpy.array([-10., -1., 1, -2., 0, -3])
>>> mf.get_occ(energy)
array([2, 2, 0, 2, 2, 2])
get_veff(mol=None, dm=None, dm_last=0, vhf_last=0, hermi=1)[source]

Hartree-Fock potential matrix for the given density matrix

Args:

mol : an instance of Mole

dm
: ndarray or list of ndarrays
A density matrix or a list of density matrices
Kwargs:
dm_last
: ndarray or a list of ndarrays or 0
The density matrix baseline. If not 0, this function computes the increment of HF potential w.r.t. the reference HF potential matrix.
vhf_last
: ndarray or a list of ndarrays or 0
The reference HF potential matrix.
hermi
: int

Whether J, K matrix is hermitian

0 : no hermitian or symmetric
1 : hermitian
2 : anti-hermitian
vhfopt :
A class which holds precomputed quantities to optimize the computation of J, K matrices
Returns:
matrix Vhf = 2*J - K. Vhf can be a list matrices, corresponding to the input density matrices.

Examples:

>>> import numpy
>>> from pyscf import gto, scf
>>> from pyscf.scf import _vhf
>>> mol = gto.M(atom='H 0 0 0; H 0 0 1.1')
>>> dm0 = numpy.random.random((mol.nao_nr(),mol.nao_nr()))
>>> vhf0 = scf.hf.get_veff(mol, dm0, hermi=0)
>>> dm1 = numpy.random.random((mol.nao_nr(),mol.nao_nr()))
>>> vhf1 = scf.hf.get_veff(mol, dm1, hermi=0)
>>> vhf2 = scf.hf.get_veff(mol, dm1, dm_last=dm0, vhf_last=vhf0, hermi=0)
>>> numpy.allclose(vhf1, vhf2)
True
init_guess_by_1e(mol=None)[source]

Generate initial guess density matrix from core hamiltonian

Returns:
Density matrix, 2D ndarray
init_guess_by_atom(mol=None)[source]

Generate initial guess density matrix from superposition of atomic HF density matrix. The atomic HF is occupancy averaged RHF

Returns:
Density matrix, 2D ndarray
init_guess_by_chkfile(chkfile=None, project=None)[source]

Read the HF results from checkpoint file, then project it to the basis defined by mol

Returns:
Density matrix, 2D ndarray
init_guess_by_minao(mol=None)[source]

Generate initial guess density matrix based on ANO basis, then project the density matrix to the basis set defined by mol

Returns:
Density matrix, 2D ndarray

Examples:

>>> from pyscf import gto, scf
>>> mol = gto.M(atom='H 0 0 0; H 0 0 1.1')
>>> scf.hf.init_guess_by_minao(mol)
array([[ 0.94758917,  0.09227308],
       [ 0.09227308,  0.94758917]])
kernel(*args, **kwargs)

SCF main driver

Kwargs:
dm0
: ndarray
If given, it will be used as the initial guess density matrix

Examples:

>>> import numpy
>>> from pyscf import gto, scf
>>> mol = gto.M(atom='H 0 0 0; F 0 0 1.1')
>>> mf = scf.hf.SCF(mol)
>>> dm_guess = numpy.eye(mol.nao_nr())
>>> mf.kernel(dm_guess)
converged SCF energy = -98.5521904482821
-98.552190448282104
make_rdm1(mo_coeff=None, mo_occ=None, **kwargs)[source]

One-particle density matrix in AO representation

Args:
mo_coeff
: 2D ndarray
Orbital coefficients. Each column is one orbital.
mo_occ
: 1D ndarray
Occupancy
mulliken_meta(mol=None, dm=None, verbose=5, pre_orth_method='ANO', s=None)[source]

Mulliken population analysis, based on meta-Lowdin AOs. In the meta-lowdin, the AOs are grouped in three sets: core, valence and Rydberg, the orthogonalization are carreid out within each subsets.

Args:

mol : an instance of Mole

dm
: ndarray or 2-item list of ndarray
Density matrix. ROHF dm is a 2-item list of 2D array
Kwargs:

verbose : int or instance of lib.logger.Logger

pre_orth_method
: str

Pre-orthogonalization, which localized GTOs for each atom. To obtain the occupied and unoccupied atomic shells, there are three methods

‘ano’ : Project GTOs to ANO basis
‘minao’ : Project GTOs to MINAO basis
‘scf’ : Fraction-averaged RHF
mulliken_pop(mol=None, dm=None, s=None, verbose=5)[source]

Mulliken population analysis

\[M_{ij} = D_{ij} S_{ji}\]

Mulliken charges

\[\delta_i = \sum_j M_{ij}\]
mulliken_pop_meta_lowdin_ao(*args, **kwargs)

Mulliken population analysis, based on meta-Lowdin AOs. In the meta-lowdin, the AOs are grouped in three sets: core, valence and Rydberg, the orthogonalization are carreid out within each subsets.

Args:

mol : an instance of Mole

dm
: ndarray or 2-item list of ndarray
Density matrix. ROHF dm is a 2-item list of 2D array
Kwargs:

verbose : int or instance of lib.logger.Logger

pre_orth_method
: str

Pre-orthogonalization, which localized GTOs for each atom. To obtain the occupied and unoccupied atomic shells, there are three methods

‘ano’ : Project GTOs to ANO basis
‘minao’ : Project GTOs to MINAO basis
‘scf’ : Fraction-averaged RHF
nuc_grad_method()[source]

Hook to create object for analytical nuclear gradients.

pop(*args, **kwargs)[source]

Mulliken population analysis, based on meta-Lowdin AOs. In the meta-lowdin, the AOs are grouped in three sets: core, valence and Rydberg, the orthogonalization are carreid out within each subsets.

Args:

mol : an instance of Mole

dm
: ndarray or 2-item list of ndarray
Density matrix. ROHF dm is a 2-item list of 2D array
Kwargs:

verbose : int or instance of lib.logger.Logger

pre_orth_method
: str

Pre-orthogonalization, which localized GTOs for each atom. To obtain the occupied and unoccupied atomic shells, there are three methods

‘ano’ : Project GTOs to ANO basis
‘minao’ : Project GTOs to MINAO basis
‘scf’ : Fraction-averaged RHF
scf(dm0=None, **kwargs)[source]

SCF main driver

Kwargs:
dm0
: ndarray
If given, it will be used as the initial guess density matrix

Examples:

>>> import numpy
>>> from pyscf import gto, scf
>>> mol = gto.M(atom='H 0 0 0; F 0 0 1.1')
>>> mf = scf.hf.SCF(mol)
>>> dm_guess = numpy.eye(mol.nao_nr())
>>> mf.kernel(dm_guess)
converged SCF energy = -98.5521904482821
-98.552190448282104
update(chkfile=None)

Read attributes from the chkfile then replace the attributes of current object. It’s an alias of function update_from_chk_.

update_(chkfile=None)[source]

Read attributes from the chkfile then replace the attributes of current object. It’s an alias of function update_from_chk_.

update_from_chk(chkfile=None)

Read attributes from the chkfile then replace the attributes of current object. It’s an alias of function update_from_chk_.

update_from_chk_(chkfile=None)

Read attributes from the chkfile then replace the attributes of current object. It’s an alias of function update_from_chk_.


class pyscf.scf.hf.RHF(mol)[source]

SCF base class. non-relativistic RHF.

Attributes:
verbose
: int
Print level. Default value equals to Mole.verbose
max_memory
: float or int
Allowed memory in MB. Default equals to Mole.max_memory
chkfile
: str
checkpoint file to save MOs, orbital energies etc. Writing to chkfile can be disabled if this attribute is set to None or False.
conv_tol
: float
converge threshold. Default is 1e-9
conv_tol_grad
: float
gradients converge threshold. Default is sqrt(conv_tol)
max_cycle
: int
max number of iterations. Default is 50
init_guess
: str
initial guess method. It can be one of ‘minao’, ‘atom’, ‘hcore’, ‘1e’, ‘chkfile’. Default is ‘minao’
DIIS
: DIIS class
The class to generate diis object. It can be one of diis.SCF_DIIS, diis.ADIIS, diis.EDIIS.
diis
: boolean or object of DIIS class defined in scf.diis.
Default is the object associated to the attribute self.DIIS. Set it to None/False to turn off DIIS. Note if this attribute is inialized as a DIIS object, the SCF driver will use this object in the iteration. The DIIS informations (vector basis and error vector) will be held inside this object. When kernel function is called again, the old states (vector basis and error vector) will be reused.
diis_space
: int
DIIS space size. By default, 8 Fock matrices and errors vector are stored.
diis_start_cycle
: int
The step to start DIIS. Default is 1.
diis_file: ‘str’
File to store DIIS vectors and error vectors.
level_shift
: float or int
Level shift (in AU) for virtual space. Default is 0.
direct_scf
: bool
Direct SCF is used by default.
direct_scf_tol
: float
Direct SCF cutoff threshold. Default is 1e-13.
callback
: function(envs_dict) => None
callback function takes one dict as the argument which is generated by the builtin function locals(), so that the callback function can access all local variables in the current envrionment.
conv_check
: bool
An extra cycle to check convergence after SCF iterations.
check_convergence
: function(envs) => bool
A hook for overloading convergence criteria in SCF iterations.

Saved results:

converged
: bool
SCF converged or not
e_tot
: float
Total HF energy (electronic energy plus nuclear repulsion)
mo_energy :
Orbital energies
mo_occ
Orbital occupancy
mo_coeff
Orbital coefficients

Examples:

>>> mol = gto.M(atom='H 0 0 0; H 0 0 1.1', basis='cc-pvdz')
>>> mf = scf.hf.SCF(mol)
>>> mf.verbose = 0
>>> mf.level_shift = .4
>>> mf.scf()
-1.0811707843775884
convert_from_(mf)[source]

Convert given mean-field object to RHF/ROHF

get_jk(mol=None, dm=None, hermi=1)[source]

Compute J, K matrices for the given density matrix

Args:

mol : an instance of Mole

dm
: ndarray or list of ndarrays
A density matrix or a list of density matrices
Kwargs:
hermi
: int

Whether J, K matrix is hermitian

0 : not hermitian and not symmetric
1 : hermitian or symmetric
2 : anti-hermitian
vhfopt :
A class which holds precomputed quantities to optimize the computation of J, K matrices
Returns:
Depending on the given dm, the function returns one J and one K matrix, or a list of J matrices and a list of K matrices, corresponding to the input density matrices.

Examples:

>>> from pyscf import gto, scf
>>> from pyscf.scf import _vhf
>>> mol = gto.M(atom='H 0 0 0; H 0 0 1.1')
>>> dms = numpy.random.random((3,mol.nao_nr(),mol.nao_nr()))
>>> j, k = scf.hf.get_jk(mol, dms, hermi=0)
>>> print(j.shape)
(3, 2, 2)
spin_square(mo_coeff=None, s=None)[source]

Spin square and multiplicity of RHF determinant

stability(internal=True, external=False, verbose=None)[source]

RHF/RKS stability analysis.

See also pyscf.scf.stability.rhf_stability function.

Kwargs:
internal
: bool
Internal stability, within the RHF optimization space.
external
: bool
External stability. Including the RHF -> UHF and real -> complex stability analysis.
Returns:
New orbitals that are more close to the stable condition. The return value includes two set of orbitals. The first corresponds to the internal stablity and the second corresponds to the external stability.

class pyscf.scf.rohf.ROHF(mol)[source]

SCF base class. non-relativistic RHF.

Attributes:
verbose
: int
Print level. Default value equals to Mole.verbose
max_memory
: float or int
Allowed memory in MB. Default equals to Mole.max_memory
chkfile
: str
checkpoint file to save MOs, orbital energies etc. Writing to chkfile can be disabled if this attribute is set to None or False.
conv_tol
: float
converge threshold. Default is 1e-9
conv_tol_grad
: float
gradients converge threshold. Default is sqrt(conv_tol)
max_cycle
: int
max number of iterations. Default is 50
init_guess
: str
initial guess method. It can be one of ‘minao’, ‘atom’, ‘hcore’, ‘1e’, ‘chkfile’. Default is ‘minao’
DIIS
: DIIS class
The class to generate diis object. It can be one of diis.SCF_DIIS, diis.ADIIS, diis.EDIIS.
diis
: boolean or object of DIIS class defined in scf.diis.
Default is the object associated to the attribute self.DIIS. Set it to None/False to turn off DIIS. Note if this attribute is inialized as a DIIS object, the SCF driver will use this object in the iteration. The DIIS informations (vector basis and error vector) will be held inside this object. When kernel function is called again, the old states (vector basis and error vector) will be reused.
diis_space
: int
DIIS space size. By default, 8 Fock matrices and errors vector are stored.
diis_start_cycle
: int
The step to start DIIS. Default is 1.
diis_file: ‘str’
File to store DIIS vectors and error vectors.
level_shift
: float or int
Level shift (in AU) for virtual space. Default is 0.
direct_scf
: bool
Direct SCF is used by default.
direct_scf_tol
: float
Direct SCF cutoff threshold. Default is 1e-13.
callback
: function(envs_dict) => None
callback function takes one dict as the argument which is generated by the builtin function locals(), so that the callback function can access all local variables in the current envrionment.
conv_check
: bool
An extra cycle to check convergence after SCF iterations.
check_convergence
: function(envs) => bool
A hook for overloading convergence criteria in SCF iterations.

Saved results:

converged
: bool
SCF converged or not
e_tot
: float
Total HF energy (electronic energy plus nuclear repulsion)
mo_energy :
Orbital energies
mo_occ
Orbital occupancy
mo_coeff
Orbital coefficients

Examples:

>>> mol = gto.M(atom='H 0 0 0; H 0 0 1.1', basis='cc-pvdz')
>>> mf = scf.hf.SCF(mol)
>>> mf.verbose = 0
>>> mf.level_shift = .4
>>> mf.scf()
-1.0811707843775884
analyze(verbose=None, with_meta_lowdin=True, **kwargs)[source]

Analyze the given SCF object: print orbital energies, occupancies; print orbital coefficients; Mulliken population analysis

canonicalize(mf, mo_coeff, mo_occ, fock=None)

Canonicalization diagonalizes the Fock matrix within occupied, open, virtual subspaces separatedly (without change occupancy).

check_sanity()

Check input of class/object attributes, check whether a class method is overwritten. It does not check the attributes which are prefixed with “_”. The return value of method set is the object itself. This allows a series of functions/methods to be executed in pipe.

eig(fock, s)[source]

Solver for generalized eigenvalue problem

\[HC = SCE\]
get_fock(mf, h1e=None, s1e=None, vhf=None, dm=None, cycle=-1, diis=None, diis_start_cycle=None, level_shift_factor=None, damp_factor=None)

Build fock matrix based on Roothaan’s effective fock. See also get_roothaan_fock()

get_grad(mo_coeff, mo_occ, fock=None)[source]

ROHF gradients is the off-diagonal block [co + cv + ov], where [ cc co cv ] [ oc oo ov ] [ vc vo vv ]

get_occ(mf, mo_energy=None, mo_coeff=None)

Label the occupancies for each orbital. NOTE the occupancies are not assigned based on the orbital energy ordering. The first N orbitals are assigned to be occupied orbitals.

Examples:

>>> mol = gto.M(atom='H 0 0 0; O 0 0 1.1', spin=1)
>>> mf = scf.hf.SCF(mol)
>>> energy = numpy.array([-10., -1., 1, -2., 0, -3])
>>> mf.get_occ(energy)
array([2, 2, 2, 2, 1, 0])
get_veff(mol=None, dm=None, dm_last=0, vhf_last=0, hermi=1)[source]

Unrestricted Hartree-Fock potential matrix of alpha and beta spins, for the given density matrix

\[\begin{split}V_{ij}^\alpha &= \sum_{kl} (ij|kl)(\gamma_{lk}^\alpha+\gamma_{lk}^\beta) - \sum_{kl} (il|kj)\gamma_{lk}^\alpha \\ V_{ij}^\beta &= \sum_{kl} (ij|kl)(\gamma_{lk}^\alpha+\gamma_{lk}^\beta) - \sum_{kl} (il|kj)\gamma_{lk}^\beta\end{split}\]
Args:

mol : an instance of Mole

dm
: a list of ndarrays
A list of density matrices, stored as (alpha,alpha,...,beta,beta,...)
Kwargs:
dm_last
: ndarray or a list of ndarrays or 0
The density matrix baseline. When it is not 0, this function computes the increment of HF potential w.r.t. the reference HF potential matrix.
vhf_last
: ndarray or a list of ndarrays or 0
The reference HF potential matrix.
hermi
: int

Whether J, K matrix is hermitian

0 : no hermitian or symmetric
1 : hermitian
2 : anti-hermitian
vhfopt :
A class which holds precomputed quantities to optimize the computation of J, K matrices
Returns:
\(V_{hf} = (V^\alpha, V^\beta)\). \(V^\alpha\) (and \(V^\beta\)) can be a list matrices, corresponding to the input density matrices.

Examples:

>>> import numpy
>>> from pyscf import gto, scf
>>> mol = gto.M(atom='H 0 0 0; H 0 0 1.1')
>>> dmsa = numpy.random.random((3,mol.nao_nr(),mol.nao_nr()))
>>> dmsb = numpy.random.random((3,mol.nao_nr(),mol.nao_nr()))
>>> dms = numpy.vstack((dmsa,dmsb))
>>> dms.shape
(6, 2, 2)
>>> vhfa, vhfb = scf.uhf.get_veff(mol, dms, hermi=0)
>>> vhfa.shape
(3, 2, 2)
>>> vhfb.shape
(3, 2, 2)
make_rdm1(mo_coeff=None, mo_occ=None, **kwargs)[source]

One-particle densit matrix. mo_occ is a 1D array, with occupancy 1 or 2.

spin_square(mo_coeff=None, s=None)[source]

Spin square and multiplicity of RHF determinant

stability(internal=True, external=False, verbose=None)[source]

ROHF/ROKS stability analysis.

See also pyscf.scf.stability.rohf_stability function.

Kwargs:
internal
: bool
Internal stability, within the RHF optimization space.
external
: bool
External stability. It is not available in current version.
Returns:
The return value includes two set of orbitals which are more close to the required stable condition.

class pyscf.scf.uhf.UHF(mol)[source]

SCF base class. non-relativistic RHF.

Attributes:
verbose
: int
Print level. Default value equals to Mole.verbose
max_memory
: float or int
Allowed memory in MB. Default equals to Mole.max_memory
chkfile
: str
checkpoint file to save MOs, orbital energies etc. Writing to chkfile can be disabled if this attribute is set to None or False.
conv_tol
: float
converge threshold. Default is 1e-9
conv_tol_grad
: float
gradients converge threshold. Default is sqrt(conv_tol)
max_cycle
: int
max number of iterations. Default is 50
init_guess
: str
initial guess method. It can be one of ‘minao’, ‘atom’, ‘hcore’, ‘1e’, ‘chkfile’. Default is ‘minao’
DIIS
: DIIS class
The class to generate diis object. It can be one of diis.SCF_DIIS, diis.ADIIS, diis.EDIIS.
diis
: boolean or object of DIIS class defined in scf.diis.
Default is the object associated to the attribute self.DIIS. Set it to None/False to turn off DIIS. Note if this attribute is inialized as a DIIS object, the SCF driver will use this object in the iteration. The DIIS informations (vector basis and error vector) will be held inside this object. When kernel function is called again, the old states (vector basis and error vector) will be reused.
diis_space
: int
DIIS space size. By default, 8 Fock matrices and errors vector are stored.
diis_start_cycle
: int
The step to start DIIS. Default is 1.
diis_file: ‘str’
File to store DIIS vectors and error vectors.
level_shift
: float or int
Level shift (in AU) for virtual space. Default is 0.
direct_scf
: bool
Direct SCF is used by default.
direct_scf_tol
: float
Direct SCF cutoff threshold. Default is 1e-13.
callback
: function(envs_dict) => None
callback function takes one dict as the argument which is generated by the builtin function locals(), so that the callback function can access all local variables in the current envrionment.
conv_check
: bool
An extra cycle to check convergence after SCF iterations.
check_convergence
: function(envs) => bool
A hook for overloading convergence criteria in SCF iterations.

Saved results:

converged
: bool
SCF converged or not
e_tot
: float
Total HF energy (electronic energy plus nuclear repulsion)
mo_energy :
Orbital energies
mo_occ
Orbital occupancy
mo_coeff
Orbital coefficients

Examples:

>>> mol = gto.M(atom='H 0 0 0; H 0 0 1.1', basis='cc-pvdz')
>>> mf = scf.hf.SCF(mol)
>>> mf.verbose = 0
>>> mf.level_shift = .4
>>> mf.scf()
-1.0811707843775884
Attributes for UHF:
nelec
: (int, int)
If given, freeze the number of (alpha,beta) electrons to the given value.
level_shift
: number or two-element list
level shift (in Eh) for alpha and beta Fock if two-element list is given.

Examples:

>>> mol = gto.M(atom='O 0 0 0; H 0 0 1; H 0 1 0', basis='ccpvdz', charge=1, spin=1, verbose=0)
>>> mf = scf.UHF(mol)
>>> mf.kernel()
-75.623975516256706
>>> print('S^2 = %.7f, 2S+1 = %.7f' % mf.spin_square())
S^2 = 0.7570150, 2S+1 = 2.0070027
canonicalize(mf, mo_coeff, mo_occ, fock=None)

Canonicalization diagonalizes the UHF Fock matrix within occupied, virtual subspaces separatedly (without change occupancy).

convert_from_(mf)[source]

Create UHF object based on the RHF/ROHF object

det_ovlp(mo1, mo2, occ1, occ2, ovlp=None)[source]

Calculate the overlap between two different determinants. It is the product of single values of molecular orbital overlap matrix.

\[S_{12} = \langle \Psi_A | \Psi_B \rangle = (\mathrm{det}\mathbf{U}) (\mathrm{det}\mathbf{V^\dagger})\prod\limits_{i=1}\limits^{2N} \lambda_{ii}\]

where \(\mathbf{U}, \mathbf{V}, \lambda\) are unitary matrices and single values generated by single value decomposition(SVD) of the overlap matrix \(\mathbf{O}\) which is the overlap matrix of two sets of molecular orbitals:

\[\mathbf{U}^\dagger \mathbf{O} \mathbf{V} = \mathbf{\Lambda}\]
Args:
mo1, mo2
: 2D ndarrays
Molecualr orbital coefficients
occ1, occ2: 2D ndarrays
occupation numbers
Return:
A list: the product of single values: float
x_a, x_b: 1D ndarrays \(\mathbf{U} \mathbf{\Lambda}^{-1} \mathbf{V}^\dagger\) They are used to calculate asymmetric density matrix
energy_elec(mf, dm=None, h1e=None, vhf=None)

Electronic energy of Unrestricted Hartree-Fock

Returns:
Hartree-Fock electronic energy and the 2-electron part contribution
get_jk(mol=None, dm=None, hermi=1)[source]

Coulomb (J) and exchange (K)

Args:
dm
: a list of 2D arrays or a list of 3D arrays
(alpha_dm, beta_dm) or (alpha_dms, beta_dms)
get_veff(mol=None, dm=None, dm_last=0, vhf_last=0, hermi=1)[source]

Unrestricted Hartree-Fock potential matrix of alpha and beta spins, for the given density matrix

\[\begin{split}V_{ij}^\alpha &= \sum_{kl} (ij|kl)(\gamma_{lk}^\alpha+\gamma_{lk}^\beta) - \sum_{kl} (il|kj)\gamma_{lk}^\alpha \\ V_{ij}^\beta &= \sum_{kl} (ij|kl)(\gamma_{lk}^\alpha+\gamma_{lk}^\beta) - \sum_{kl} (il|kj)\gamma_{lk}^\beta\end{split}\]
Args:

mol : an instance of Mole

dm
: a list of ndarrays
A list of density matrices, stored as (alpha,alpha,...,beta,beta,...)
Kwargs:
dm_last
: ndarray or a list of ndarrays or 0
The density matrix baseline. When it is not 0, this function computes the increment of HF potential w.r.t. the reference HF potential matrix.
vhf_last
: ndarray or a list of ndarrays or 0
The reference HF potential matrix.
hermi
: int

Whether J, K matrix is hermitian

0 : no hermitian or symmetric
1 : hermitian
2 : anti-hermitian
vhfopt :
A class which holds precomputed quantities to optimize the computation of J, K matrices
Returns:
\(V_{hf} = (V^\alpha, V^\beta)\). \(V^\alpha\) (and \(V^\beta\)) can be a list matrices, corresponding to the input density matrices.

Examples:

>>> import numpy
>>> from pyscf import gto, scf
>>> mol = gto.M(atom='H 0 0 0; H 0 0 1.1')
>>> dmsa = numpy.random.random((3,mol.nao_nr(),mol.nao_nr()))
>>> dmsb = numpy.random.random((3,mol.nao_nr(),mol.nao_nr()))
>>> dms = numpy.vstack((dmsa,dmsb))
>>> dms.shape
(6, 2, 2)
>>> vhfa, vhfb = scf.uhf.get_veff(mol, dms, hermi=0)
>>> vhfa.shape
(3, 2, 2)
>>> vhfb.shape
(3, 2, 2)
init_guess_by_minao(mol=None, breaksym=True)[source]

Initial guess in terms of the overlap to minimal basis.

make_asym_dm(mo1, mo2, occ1, occ2, x)[source]

One-particle asymmetric density matrix

Args:
mo1, mo2
: 2D ndarrays
Molecualr orbital coefficients
occ1, occ2: 2D ndarrays
Occupation numbers
x: 2D ndarrays
\(\mathbf{U} \mathbf{\Lambda}^{-1} \mathbf{V}^\dagger\). See also det_ovlp()
Return:
A list of 2D ndarrays for alpha and beta spin

Examples:

>>> mf1 = scf.UHF(gto.M(atom='H 0 0 0; F 0 0 1.3', basis='ccpvdz')).run()
>>> mf2 = scf.UHF(gto.M(atom='H 0 0 0; F 0 0 1.4', basis='ccpvdz')).run()
>>> s = gto.intor_cross('int1e_ovlp_sph', mf1.mol, mf2.mol)
>>> det, x = det_ovlp(mf1.mo_coeff, mf1.mo_occ, mf2.mo_coeff, mf2.mo_occ, s)
>>> adm = make_asym_dm(mf1.mo_coeff, mf1.mo_occ, mf2.mo_coeff, mf2.mo_occ, x)
>>> adm.shape
(2, 19, 19)
make_rdm1(mo_coeff=None, mo_occ=None, **kwargs)[source]

One-particle density matrix

Returns:
A list of 2D ndarrays for alpha and beta spins
spin_square(mo_coeff=None, s=None)[source]

Spin square and multiplicity of UHF determinant

\[S^2 = \frac{1}{2}(S_+ S_- + S_- S_+) + S_z^2\]

where \(S_+ = \sum_i S_{i+}\) is effective for all beta occupied orbitals; \(S_- = \sum_i S_{i-}\) is effective for all alpha occupied orbitals.

  1. There are two possibilities for \(S_+ S_-\)
    1. same electron \(S_+ S_- = \sum_i s_{i+} s_{i-}\),
    \[\sum_i \langle UHF|s_{i+} s_{i-}|UHF\rangle = \sum_{pq}\langle p|s_+s_-|q\rangle \gamma_{qp} = n_\alpha\]

    2) different electrons \(S_+ S_- = \sum s_{i+} s_{j-}, (i\neq j)\). There are in total \(n(n-1)\) terms. As a two-particle operator,

    \[\langle S_+ S_- \rangle = \langle ij|s_+ s_-|ij\rangle - \langle ij|s_+ s_-|ji\rangle = -\langle i^\alpha|j^\beta\rangle \langle j^\beta|i^\alpha\rangle\]
  2. Similarly, for \(S_- S_+\)
    1. same electron
    \[\sum_i \langle s_{i-} s_{i+}\rangle = n_\beta\]
    1. different electrons
    \[\langle S_- S_+ \rangle = -\langle i^\beta|j^\alpha\rangle \langle j^\alpha|i^\beta\rangle\]
  3. For \(S_z^2\)
    1. same electron
    \[\langle s_z^2\rangle = \frac{1}{4}(n_\alpha + n_\beta)\]
    1. different electrons
    \[\begin{split}&\frac{1}{2}\sum_{ij}(\langle ij|2s_{z1}s_{z2}|ij\rangle -\langle ij|2s_{z1}s_{z2}|ji\rangle) \\ &=\frac{1}{4}(\langle i^\alpha|i^\alpha\rangle \langle j^\alpha|j^\alpha\rangle - \langle i^\alpha|i^\alpha\rangle \langle j^\beta|j^\beta\rangle - \langle i^\beta|i^\beta\rangle \langle j^\alpha|j^\alpha\rangle + \langle i^\beta|i^\beta\rangle \langle j^\beta|j^\beta\rangle) \\ &-\frac{1}{4}(\langle i^\alpha|j^\alpha\rangle \langle j^\alpha|i^\alpha\rangle + \langle i^\beta|j^\beta\rangle\langle j^\beta|i^\beta\rangle) \\ &=\frac{1}{4}(n_\alpha^2 - n_\alpha n_\beta - n_\beta n_\alpha + n_\beta^2) -\frac{1}{4}(n_\alpha + n_\beta) \\ &=\frac{1}{4}((n_\alpha-n_\beta)^2 - (n_\alpha+n_\beta))\end{split}\]

In total

\[\begin{split}\langle S^2\rangle &= \frac{1}{2} (n_\alpha-\sum_{ij}\langle i^\alpha|j^\beta\rangle \langle j^\beta|i^\alpha\rangle +n_\beta -\sum_{ij}\langle i^\beta|j^\alpha\rangle\langle j^\alpha|i^\beta\rangle) + \frac{1}{4}(n_\alpha-n_\beta)^2 \\\end{split}\]
Args:
mo
: a list of 2 ndarrays
Occupied alpha and occupied beta orbitals
Kwargs:
s
: ndarray
AO overlap
Returns:
A list of two floats. The first is the expectation value of S^2. The second is the corresponding 2S+1

Examples:

>>> mol = gto.M(atom='O 0 0 0; H 0 0 1; H 0 1 0', basis='ccpvdz', charge=1, spin=1, verbose=0)
>>> mf = scf.UHF(mol)
>>> mf.kernel()
-75.623975516256706
>>> mo = (mf.mo_coeff[0][:,mf.mo_occ[0]>0], mf.mo_coeff[1][:,mf.mo_occ[1]>0])
>>> print('S^2 = %.7f, 2S+1 = %.7f' % spin_square(mo, mol.intor('int1e_ovlp_sph')))
S^2 = 0.7570150, 2S+1 = 2.0070027
stability(internal=True, external=False, verbose=None)[source]

Stability analysis for RHF/RKS method.

See also pyscf.scf.stability.uhf_stability function.

Args:
mf : UHF or UKS object
Kwargs:
internal
: bool
Internal stability, within the UHF space.
external
: bool
External stability. Including the UHF -> GHF and real -> complex stability analysis.
Returns:
New orbitals that are more close to the stable condition. The return value includes two set of orbitals. The first corresponds to the internal stablity and the second corresponds to the external stability.

Hartree-Fock

class pyscf.scf.hf.RHF(mol)[source]

SCF base class. non-relativistic RHF.

Attributes:
verbose
: int
Print level. Default value equals to Mole.verbose
max_memory
: float or int
Allowed memory in MB. Default equals to Mole.max_memory
chkfile
: str
checkpoint file to save MOs, orbital energies etc. Writing to chkfile can be disabled if this attribute is set to None or False.
conv_tol
: float
converge threshold. Default is 1e-9
conv_tol_grad
: float
gradients converge threshold. Default is sqrt(conv_tol)
max_cycle
: int
max number of iterations. Default is 50
init_guess
: str
initial guess method. It can be one of ‘minao’, ‘atom’, ‘hcore’, ‘1e’, ‘chkfile’. Default is ‘minao’
DIIS
: DIIS class
The class to generate diis object. It can be one of diis.SCF_DIIS, diis.ADIIS, diis.EDIIS.
diis
: boolean or object of DIIS class defined in scf.diis.
Default is the object associated to the attribute self.DIIS. Set it to None/False to turn off DIIS. Note if this attribute is inialized as a DIIS object, the SCF driver will use this object in the iteration. The DIIS informations (vector basis and error vector) will be held inside this object. When kernel function is called again, the old states (vector basis and error vector) will be reused.
diis_space
: int
DIIS space size. By default, 8 Fock matrices and errors vector are stored.
diis_start_cycle
: int
The step to start DIIS. Default is 1.
diis_file: ‘str’
File to store DIIS vectors and error vectors.
level_shift
: float or int
Level shift (in AU) for virtual space. Default is 0.
direct_scf
: bool
Direct SCF is used by default.
direct_scf_tol
: float
Direct SCF cutoff threshold. Default is 1e-13.
callback
: function(envs_dict) => None
callback function takes one dict as the argument which is generated by the builtin function locals(), so that the callback function can access all local variables in the current envrionment.
conv_check
: bool
An extra cycle to check convergence after SCF iterations.
check_convergence
: function(envs) => bool
A hook for overloading convergence criteria in SCF iterations.

Saved results:

converged
: bool
SCF converged or not
e_tot
: float
Total HF energy (electronic energy plus nuclear repulsion)
mo_energy :
Orbital energies
mo_occ
Orbital occupancy
mo_coeff
Orbital coefficients

Examples:

>>> mol = gto.M(atom='H 0 0 0; H 0 0 1.1', basis='cc-pvdz')
>>> mf = scf.hf.SCF(mol)
>>> mf.verbose = 0
>>> mf.level_shift = .4
>>> mf.scf()
-1.0811707843775884
convert_from_(mf)[source]

Convert given mean-field object to RHF/ROHF

get_jk(mol=None, dm=None, hermi=1)[source]

Compute J, K matrices for the given density matrix

Args:

mol : an instance of Mole

dm
: ndarray or list of ndarrays
A density matrix or a list of density matrices
Kwargs:
hermi
: int

Whether J, K matrix is hermitian

0 : not hermitian and not symmetric
1 : hermitian or symmetric
2 : anti-hermitian
vhfopt :
A class which holds precomputed quantities to optimize the computation of J, K matrices
Returns:
Depending on the given dm, the function returns one J and one K matrix, or a list of J matrices and a list of K matrices, corresponding to the input density matrices.

Examples:

>>> from pyscf import gto, scf
>>> from pyscf.scf import _vhf
>>> mol = gto.M(atom='H 0 0 0; H 0 0 1.1')
>>> dms = numpy.random.random((3,mol.nao_nr(),mol.nao_nr()))
>>> j, k = scf.hf.get_jk(mol, dms, hermi=0)
>>> print(j.shape)
(3, 2, 2)
spin_square(mo_coeff=None, s=None)[source]

Spin square and multiplicity of RHF determinant

stability(internal=True, external=False, verbose=None)[source]

RHF/RKS stability analysis.

See also pyscf.scf.stability.rhf_stability function.

Kwargs:
internal
: bool
Internal stability, within the RHF optimization space.
external
: bool
External stability. Including the RHF -> UHF and real -> complex stability analysis.
Returns:
New orbitals that are more close to the stable condition. The return value includes two set of orbitals. The first corresponds to the internal stablity and the second corresponds to the external stability.
class pyscf.scf.hf.SCF(mol)[source]

SCF base class. non-relativistic RHF.

Attributes:
verbose
: int
Print level. Default value equals to Mole.verbose
max_memory
: float or int
Allowed memory in MB. Default equals to Mole.max_memory
chkfile
: str
checkpoint file to save MOs, orbital energies etc. Writing to chkfile can be disabled if this attribute is set to None or False.
conv_tol
: float
converge threshold. Default is 1e-9
conv_tol_grad
: float
gradients converge threshold. Default is sqrt(conv_tol)
max_cycle
: int
max number of iterations. Default is 50
init_guess
: str
initial guess method. It can be one of ‘minao’, ‘atom’, ‘hcore’, ‘1e’, ‘chkfile’. Default is ‘minao’
DIIS
: DIIS class
The class to generate diis object. It can be one of diis.SCF_DIIS, diis.ADIIS, diis.EDIIS.
diis
: boolean or object of DIIS class defined in scf.diis.
Default is the object associated to the attribute self.DIIS. Set it to None/False to turn off DIIS. Note if this attribute is inialized as a DIIS object, the SCF driver will use this object in the iteration. The DIIS informations (vector basis and error vector) will be held inside this object. When kernel function is called again, the old states (vector basis and error vector) will be reused.
diis_space
: int
DIIS space size. By default, 8 Fock matrices and errors vector are stored.
diis_start_cycle
: int
The step to start DIIS. Default is 1.
diis_file: ‘str’
File to store DIIS vectors and error vectors.
level_shift
: float or int
Level shift (in AU) for virtual space. Default is 0.
direct_scf
: bool
Direct SCF is used by default.
direct_scf_tol
: float
Direct SCF cutoff threshold. Default is 1e-13.
callback
: function(envs_dict) => None
callback function takes one dict as the argument which is generated by the builtin function locals(), so that the callback function can access all local variables in the current envrionment.
conv_check
: bool
An extra cycle to check convergence after SCF iterations.
check_convergence
: function(envs) => bool
A hook for overloading convergence criteria in SCF iterations.

Saved results:

converged
: bool
SCF converged or not
e_tot
: float
Total HF energy (electronic energy plus nuclear repulsion)
mo_energy :
Orbital energies
mo_occ
Orbital occupancy
mo_coeff
Orbital coefficients

Examples:

>>> mol = gto.M(atom='H 0 0 0; H 0 0 1.1', basis='cc-pvdz')
>>> mf = scf.hf.SCF(mol)
>>> mf.verbose = 0
>>> mf.level_shift = .4
>>> mf.scf()
-1.0811707843775884
analyze(verbose=None, with_meta_lowdin=True, **kwargs)[source]

Analyze the given SCF object: print orbital energies, occupancies; print orbital coefficients; Mulliken population analysis; Diople moment.

as_scanner(mf)

Generating a scanner/solver for HF PES.

The returned solver is a function. This function requires one argument “mol” as input and returns total HF energy.

The solver will automatically use the results of last calculation as the initial guess of the new calculation. All parameters assigned in the SCF object (DIIS, conv_tol, max_memory etc) are automatically applied in the solver.

Note scanner has side effects. It may change many underlying objects (_scf, with_df, with_x2c, ...) during calculation.

Examples:

>>> from pyscf import gto, scf
>>> hf_scanner = scf.RHF(gto.Mole().set(verbose=0)).as_scanner()
>>> hf_scanner(gto.M(atom='H 0 0 0; F 0 0 1.1'))
-98.552190448277955
>>> hf_scanner(gto.M(atom='H 0 0 0; F 0 0 1.5'))
-98.414750424294368
canonicalize(mf, mo_coeff, mo_occ, fock=None)

Canonicalization diagonalizes the Fock matrix within occupied, open, virtual subspaces separatedly (without change occupancy).

dip_moment(mol=None, dm=None, unit='Debye', verbose=3, **kwargs)[source]

Dipole moment calculation

\[\begin{split}\mu_x = -\sum_{\mu}\sum_{\nu} P_{\mu\nu}(\nu|x|\mu) + \sum_A Q_A X_A\\ \mu_y = -\sum_{\mu}\sum_{\nu} P_{\mu\nu}(\nu|y|\mu) + \sum_A Q_A Y_A\\ \mu_z = -\sum_{\mu}\sum_{\nu} P_{\mu\nu}(\nu|z|\mu) + \sum_A Q_A Z_A\end{split}\]

where \(\mu_x, \mu_y, \mu_z\) are the x, y and z components of dipole moment

Args:
mol: an instance of Mole dm : a 2D ndarrays density matrices
Return:
A list: the dipole moment on x, y and z component
eig(h, s)[source]

Solver for generalized eigenvalue problem

\[HC = SCE\]
energy_elec(mf, dm=None, h1e=None, vhf=None)

Electronic part of Hartree-Fock energy, for given core hamiltonian and HF potential

... math:

E = \sum_{ij}h_{ij} \gamma_{ji}
  + \frac{1}{2}\sum_{ijkl} \gamma_{ji}\gamma_{lk} \langle ik||jl\rangle
Args:
mf : an instance of SCF class
Kwargs:
dm
: 2D ndarray
one-partical density matrix
h1e
: 2D ndarray
Core hamiltonian
vhf
: 2D ndarray
HF potential
Returns:
Hartree-Fock electronic energy and the Coulomb energy

Examples:

>>> from pyscf import gto, scf
>>> mol = gto.M(atom='H 0 0 0; H 0 0 1.1')
>>> mf = scf.RHF(mol)
>>> mf.scf()
>>> dm = mf.make_rdm1()
>>> scf.hf.energy_elec(mf, dm)
(-1.5176090667746334, 0.60917167853723675)
>>> mf.energy_elec(dm)
(-1.5176090667746334, 0.60917167853723675)
energy_tot(mf, dm=None, h1e=None, vhf=None)

Total Hartree-Fock energy, electronic part plus nuclear repulstion See scf.hf.energy_elec() for the electron part

from_chk(chkfile=None, project=None)[source]

Read the HF results from checkpoint file, then project it to the basis defined by mol

Returns:
Density matrix, 2D ndarray
get_fock(mf, h1e=None, s1e=None, vhf=None, dm=None, cycle=-1, diis=None, diis_start_cycle=None, level_shift_factor=None, damp_factor=None)

F = h^{core} + V^{HF}

Special treatment (damping, DIIS, or level shift) will be applied to the Fock matrix if diis and cycle is specified (The two parameters are passed to get_fock function during the SCF iteration)

Kwargs:
h1e
: 2D ndarray
Core hamiltonian
s1e
: 2D ndarray
Overlap matrix, for DIIS
vhf
: 2D ndarray
HF potential matrix
dm
: 2D ndarray
Density matrix, for DIIS
cycle
: int
Then present SCF iteration step, for DIIS
diis
: an object of SCF.DIIS class
DIIS object to hold intermediate Fock and error vectors
diis_start_cycle
: int
The step to start DIIS. Default is 0.
level_shift_factor
: float or int
Level shift (in AU) for virtual space. Default is 0.
get_grad(mo_coeff, mo_occ, fock=None)[source]

RHF orbital gradients

Args:
mo_coeff
: 2D ndarray
Obital coefficients
mo_occ
: 1D ndarray
Orbital occupancy
fock_ao
: 2D ndarray
Fock matrix in AO representation
Returns:
Gradients in MO representation. It’s a num_occ*num_vir vector.
get_j(mol=None, dm=None, hermi=1)[source]

Compute J matrix for the given density matrix.

get_jk(mol=None, dm=None, hermi=1)[source]

Compute J, K matrices for the given density matrix

Args:

mol : an instance of Mole

dm
: ndarray or list of ndarrays
A density matrix or a list of density matrices
Kwargs:
hermi
: int

Whether J, K matrix is hermitian

0 : not hermitian and not symmetric
1 : hermitian or symmetric
2 : anti-hermitian
vhfopt :
A class which holds precomputed quantities to optimize the computation of J, K matrices
Returns:
Depending on the given dm, the function returns one J and one K matrix, or a list of J matrices and a list of K matrices, corresponding to the input density matrices.

Examples:

>>> from pyscf import gto, scf
>>> from pyscf.scf import _vhf
>>> mol = gto.M(atom='H 0 0 0; H 0 0 1.1')
>>> dms = numpy.random.random((3,mol.nao_nr(),mol.nao_nr()))
>>> j, k = scf.hf.get_jk(mol, dms, hermi=0)
>>> print(j.shape)
(3, 2, 2)
get_k(mol=None, dm=None, hermi=1)[source]

Compute K matrix for the given density matrix.

get_occ(mf, mo_energy=None, mo_coeff=None)

Label the occupancies for each orbital

Kwargs:
mo_energy
: 1D ndarray
Obital energies
mo_coeff
: 2D ndarray
Obital coefficients

Examples:

>>> from pyscf import gto, scf
>>> mol = gto.M(atom='H 0 0 0; F 0 0 1.1')
>>> mf = scf.hf.SCF(mol)
>>> energy = numpy.array([-10., -1., 1, -2., 0, -3])
>>> mf.get_occ(energy)
array([2, 2, 0, 2, 2, 2])
get_veff(mol=None, dm=None, dm_last=0, vhf_last=0, hermi=1)[source]

Hartree-Fock potential matrix for the given density matrix

Args:

mol : an instance of Mole

dm
: ndarray or list of ndarrays
A density matrix or a list of density matrices
Kwargs:
dm_last
: ndarray or a list of ndarrays or 0
The density matrix baseline. If not 0, this function computes the increment of HF potential w.r.t. the reference HF potential matrix.
vhf_last
: ndarray or a list of ndarrays or 0
The reference HF potential matrix.
hermi
: int

Whether J, K matrix is hermitian

0 : no hermitian or symmetric
1 : hermitian
2 : anti-hermitian
vhfopt :
A class which holds precomputed quantities to optimize the computation of J, K matrices
Returns:
matrix Vhf = 2*J - K. Vhf can be a list matrices, corresponding to the input density matrices.

Examples:

>>> import numpy
>>> from pyscf import gto, scf
>>> from pyscf.scf import _vhf
>>> mol = gto.M(atom='H 0 0 0; H 0 0 1.1')
>>> dm0 = numpy.random.random((mol.nao_nr(),mol.nao_nr()))
>>> vhf0 = scf.hf.get_veff(mol, dm0, hermi=0)
>>> dm1 = numpy.random.random((mol.nao_nr(),mol.nao_nr()))
>>> vhf1 = scf.hf.get_veff(mol, dm1, hermi=0)
>>> vhf2 = scf.hf.get_veff(mol, dm1, dm_last=dm0, vhf_last=vhf0, hermi=0)
>>> numpy.allclose(vhf1, vhf2)
True
init_guess_by_1e(mol=None)[source]

Generate initial guess density matrix from core hamiltonian

Returns:
Density matrix, 2D ndarray
init_guess_by_atom(mol=None)[source]

Generate initial guess density matrix from superposition of atomic HF density matrix. The atomic HF is occupancy averaged RHF

Returns:
Density matrix, 2D ndarray
init_guess_by_chkfile(chkfile=None, project=None)[source]

Read the HF results from checkpoint file, then project it to the basis defined by mol

Returns:
Density matrix, 2D ndarray
init_guess_by_minao(mol=None)[source]

Generate initial guess density matrix based on ANO basis, then project the density matrix to the basis set defined by mol

Returns:
Density matrix, 2D ndarray

Examples:

>>> from pyscf import gto, scf
>>> mol = gto.M(atom='H 0 0 0; H 0 0 1.1')
>>> scf.hf.init_guess_by_minao(mol)
array([[ 0.94758917,  0.09227308],
       [ 0.09227308,  0.94758917]])
kernel(*args, **kwargs)

SCF main driver

Kwargs:
dm0
: ndarray
If given, it will be used as the initial guess density matrix

Examples:

>>> import numpy
>>> from pyscf import gto, scf
>>> mol = gto.M(atom='H 0 0 0; F 0 0 1.1')
>>> mf = scf.hf.SCF(mol)
>>> dm_guess = numpy.eye(mol.nao_nr())
>>> mf.kernel(dm_guess)
converged SCF energy = -98.5521904482821
-98.552190448282104
make_rdm1(mo_coeff=None, mo_occ=None, **kwargs)[source]

One-particle density matrix in AO representation

Args:
mo_coeff
: 2D ndarray
Orbital coefficients. Each column is one orbital.
mo_occ
: 1D ndarray
Occupancy
mulliken_meta(mol=None, dm=None, verbose=5, pre_orth_method='ANO', s=None)[source]

Mulliken population analysis, based on meta-Lowdin AOs. In the meta-lowdin, the AOs are grouped in three sets: core, valence and Rydberg, the orthogonalization are carreid out within each subsets.

Args:

mol : an instance of Mole

dm
: ndarray or 2-item list of ndarray
Density matrix. ROHF dm is a 2-item list of 2D array
Kwargs:

verbose : int or instance of lib.logger.Logger

pre_orth_method
: str

Pre-orthogonalization, which localized GTOs for each atom. To obtain the occupied and unoccupied atomic shells, there are three methods

‘ano’ : Project GTOs to ANO basis
‘minao’ : Project GTOs to MINAO basis
‘scf’ : Fraction-averaged RHF
mulliken_pop(mol=None, dm=None, s=None, verbose=5)[source]

Mulliken population analysis

\[M_{ij} = D_{ij} S_{ji}\]

Mulliken charges

\[\delta_i = \sum_j M_{ij}\]
mulliken_pop_meta_lowdin_ao(*args, **kwargs)

Mulliken population analysis, based on meta-Lowdin AOs. In the meta-lowdin, the AOs are grouped in three sets: core, valence and Rydberg, the orthogonalization are carreid out within each subsets.

Args:

mol : an instance of Mole

dm
: ndarray or 2-item list of ndarray
Density matrix. ROHF dm is a 2-item list of 2D array
Kwargs:

verbose : int or instance of lib.logger.Logger

pre_orth_method
: str

Pre-orthogonalization, which localized GTOs for each atom. To obtain the occupied and unoccupied atomic shells, there are three methods

‘ano’ : Project GTOs to ANO basis
‘minao’ : Project GTOs to MINAO basis
‘scf’ : Fraction-averaged RHF
nuc_grad_method()[source]

Hook to create object for analytical nuclear gradients.

pop(*args, **kwargs)[source]

Mulliken population analysis, based on meta-Lowdin AOs. In the meta-lowdin, the AOs are grouped in three sets: core, valence and Rydberg, the orthogonalization are carreid out within each subsets.

Args:

mol : an instance of Mole

dm
: ndarray or 2-item list of ndarray
Density matrix. ROHF dm is a 2-item list of 2D array
Kwargs:

verbose : int or instance of lib.logger.Logger

pre_orth_method
: str

Pre-orthogonalization, which localized GTOs for each atom. To obtain the occupied and unoccupied atomic shells, there are three methods

‘ano’ : Project GTOs to ANO basis
‘minao’ : Project GTOs to MINAO basis
‘scf’ : Fraction-averaged RHF
scf(dm0=None, **kwargs)[source]

SCF main driver

Kwargs:
dm0
: ndarray
If given, it will be used as the initial guess density matrix

Examples:

>>> import numpy
>>> from pyscf import gto, scf
>>> mol = gto.M(atom='H 0 0 0; F 0 0 1.1')
>>> mf = scf.hf.SCF(mol)
>>> dm_guess = numpy.eye(mol.nao_nr())
>>> mf.kernel(dm_guess)
converged SCF energy = -98.5521904482821
-98.552190448282104
update(chkfile=None)

Read attributes from the chkfile then replace the attributes of current object. It’s an alias of function update_from_chk_.

update_(chkfile=None)[source]

Read attributes from the chkfile then replace the attributes of current object. It’s an alias of function update_from_chk_.

update_from_chk(chkfile=None)

Read attributes from the chkfile then replace the attributes of current object. It’s an alias of function update_from_chk_.

update_from_chk_(chkfile=None)

Read attributes from the chkfile then replace the attributes of current object. It’s an alias of function update_from_chk_.

pyscf.scf.hf.analyze(mf, verbose=5, with_meta_lowdin=True, **kwargs)[source]

Analyze the given SCF object: print orbital energies, occupancies; print orbital coefficients; Mulliken population analysis; Diople moment.

pyscf.scf.hf.as_scanner(mf)[source]

Generating a scanner/solver for HF PES.

The returned solver is a function. This function requires one argument “mol” as input and returns total HF energy.

The solver will automatically use the results of last calculation as the initial guess of the new calculation. All parameters assigned in the SCF object (DIIS, conv_tol, max_memory etc) are automatically applied in the solver.

Note scanner has side effects. It may change many underlying objects (_scf, with_df, with_x2c, ...) during calculation.

Examples:

>>> from pyscf import gto, scf
>>> hf_scanner = scf.RHF(gto.Mole().set(verbose=0)).as_scanner()
>>> hf_scanner(gto.M(atom='H 0 0 0; F 0 0 1.1'))
-98.552190448277955
>>> hf_scanner(gto.M(atom='H 0 0 0; F 0 0 1.5'))
-98.414750424294368
pyscf.scf.hf.canonicalize(mf, mo_coeff, mo_occ, fock=None)[source]

Canonicalization diagonalizes the Fock matrix within occupied, open, virtual subspaces separatedly (without change occupancy).

pyscf.scf.hf.dip_moment(mol, dm, unit='Debye', verbose=3, **kwargs)[source]

Dipole moment calculation

\[\begin{split}\mu_x = -\sum_{\mu}\sum_{\nu} P_{\mu\nu}(\nu|x|\mu) + \sum_A Q_A X_A\\ \mu_y = -\sum_{\mu}\sum_{\nu} P_{\mu\nu}(\nu|y|\mu) + \sum_A Q_A Y_A\\ \mu_z = -\sum_{\mu}\sum_{\nu} P_{\mu\nu}(\nu|z|\mu) + \sum_A Q_A Z_A\end{split}\]

where \(\mu_x, \mu_y, \mu_z\) are the x, y and z components of dipole moment

Args:
mol: an instance of Mole dm : a 2D ndarrays density matrices
Return:
A list: the dipole moment on x, y and z component
pyscf.scf.hf.dot_eri_dm(eri, dm, hermi=0)[source]

Compute J, K matrices in terms of the given 2-electron integrals and density matrix:

J ~ numpy.einsum(‘pqrs,qp->rs’, eri, dm) K ~ numpy.einsum(‘pqrs,qr->ps’, eri, dm)

Args:
eri
: ndarray
8-fold or 4-fold ERIs or complex integral array with N^4 elements (N is the number of orbitals)
dm
: ndarray or list of ndarrays
A density matrix or a list of density matrices
Kwargs:
hermi
: int

Whether J, K matrix is hermitian

0 : no hermitian or symmetric
1 : hermitian
2 : anti-hermitian
Returns:
Depending on the given dm, the function returns one J and one K matrix, or a list of J matrices and a list of K matrices, corresponding to the input density matrices.

Examples:

>>> from pyscf import gto, scf
>>> from pyscf.scf import _vhf
>>> mol = gto.M(atom='H 0 0 0; H 0 0 1.1')
>>> eri = _vhf.int2e_sph(mol._atm, mol._bas, mol._env)
>>> dms = numpy.random.random((3,mol.nao_nr(),mol.nao_nr()))
>>> j, k = scf.hf.dot_eri_dm(eri, dms, hermi=0)
>>> print(j.shape)
(3, 2, 2)
pyscf.scf.hf.eig(h, s)[source]

Solver for generalized eigenvalue problem

\[HC = SCE\]
pyscf.scf.hf.energy_elec(mf, dm=None, h1e=None, vhf=None)[source]

Electronic part of Hartree-Fock energy, for given core hamiltonian and HF potential

... math:

E = \sum_{ij}h_{ij} \gamma_{ji}
  + \frac{1}{2}\sum_{ijkl} \gamma_{ji}\gamma_{lk} \langle ik||jl\rangle
Args:
mf : an instance of SCF class
Kwargs:
dm
: 2D ndarray
one-partical density matrix
h1e
: 2D ndarray
Core hamiltonian
vhf
: 2D ndarray
HF potential
Returns:
Hartree-Fock electronic energy and the Coulomb energy

Examples:

>>> from pyscf import gto, scf
>>> mol = gto.M(atom='H 0 0 0; H 0 0 1.1')
>>> mf = scf.RHF(mol)
>>> mf.scf()
>>> dm = mf.make_rdm1()
>>> scf.hf.energy_elec(mf, dm)
(-1.5176090667746334, 0.60917167853723675)
>>> mf.energy_elec(dm)
(-1.5176090667746334, 0.60917167853723675)
pyscf.scf.hf.energy_tot(mf, dm=None, h1e=None, vhf=None)[source]

Total Hartree-Fock energy, electronic part plus nuclear repulstion See scf.hf.energy_elec() for the electron part

pyscf.scf.hf.get_fock(mf, h1e=None, s1e=None, vhf=None, dm=None, cycle=-1, diis=None, diis_start_cycle=None, level_shift_factor=None, damp_factor=None)[source]

F = h^{core} + V^{HF}

Special treatment (damping, DIIS, or level shift) will be applied to the Fock matrix if diis and cycle is specified (The two parameters are passed to get_fock function during the SCF iteration)

Kwargs:
h1e
: 2D ndarray
Core hamiltonian
s1e
: 2D ndarray
Overlap matrix, for DIIS
vhf
: 2D ndarray
HF potential matrix
dm
: 2D ndarray
Density matrix, for DIIS
cycle
: int
Then present SCF iteration step, for DIIS
diis
: an object of SCF.DIIS class
DIIS object to hold intermediate Fock and error vectors
diis_start_cycle
: int
The step to start DIIS. Default is 0.
level_shift_factor
: float or int
Level shift (in AU) for virtual space. Default is 0.
pyscf.scf.hf.get_grad(mo_coeff, mo_occ, fock_ao)[source]

RHF orbital gradients

Args:
mo_coeff
: 2D ndarray
Obital coefficients
mo_occ
: 1D ndarray
Orbital occupancy
fock_ao
: 2D ndarray
Fock matrix in AO representation
Returns:
Gradients in MO representation. It’s a num_occ*num_vir vector.
pyscf.scf.hf.get_hcore(mol)[source]

Core Hamiltonian

Examples:

>>> from pyscf import gto, scf
>>> mol = gto.M(atom='H 0 0 0; H 0 0 1.1')
>>> scf.hf.get_hcore(mol)
array([[-0.93767904, -0.59316327],
       [-0.59316327, -0.93767904]])
pyscf.scf.hf.get_init_guess(mol, key='minao')[source]

Generate density matrix for initial guess

Kwargs:
key
: str
One of ‘minao’, ‘atom’, ‘hcore’, ‘1e’, ‘chkfile’.
pyscf.scf.hf.get_jk(mol, dm, hermi=1, vhfopt=None)[source]

Compute J, K matrices for the given density matrix

Args:

mol : an instance of Mole

dm
: ndarray or list of ndarrays
A density matrix or a list of density matrices
Kwargs:
hermi
: int

Whether J, K matrix is hermitian

0 : not hermitian and not symmetric
1 : hermitian or symmetric
2 : anti-hermitian
vhfopt :
A class which holds precomputed quantities to optimize the computation of J, K matrices
Returns:
Depending on the given dm, the function returns one J and one K matrix, or a list of J matrices and a list of K matrices, corresponding to the input density matrices.

Examples:

>>> from pyscf import gto, scf
>>> from pyscf.scf import _vhf
>>> mol = gto.M(atom='H 0 0 0; H 0 0 1.1')
>>> dms = numpy.random.random((3,mol.nao_nr(),mol.nao_nr()))
>>> j, k = scf.hf.get_jk(mol, dms, hermi=0)
>>> print(j.shape)
(3, 2, 2)
pyscf.scf.hf.get_occ(mf, mo_energy=None, mo_coeff=None)[source]

Label the occupancies for each orbital

Kwargs:
mo_energy
: 1D ndarray
Obital energies
mo_coeff
: 2D ndarray
Obital coefficients

Examples:

>>> from pyscf import gto, scf
>>> mol = gto.M(atom='H 0 0 0; F 0 0 1.1')
>>> mf = scf.hf.SCF(mol)
>>> energy = numpy.array([-10., -1., 1, -2., 0, -3])
>>> mf.get_occ(energy)
array([2, 2, 0, 2, 2, 2])
pyscf.scf.hf.get_ovlp(mol)[source]

Overlap matrix

pyscf.scf.hf.get_veff(mol, dm, dm_last=None, vhf_last=None, hermi=1, vhfopt=None)[source]

Hartree-Fock potential matrix for the given density matrix

Args:

mol : an instance of Mole

dm
: ndarray or list of ndarrays
A density matrix or a list of density matrices
Kwargs:
dm_last
: ndarray or a list of ndarrays or 0
The density matrix baseline. If not 0, this function computes the increment of HF potential w.r.t. the reference HF potential matrix.
vhf_last
: ndarray or a list of ndarrays or 0
The reference HF potential matrix.
hermi
: int

Whether J, K matrix is hermitian

0 : no hermitian or symmetric
1 : hermitian
2 : anti-hermitian
vhfopt :
A class which holds precomputed quantities to optimize the computation of J, K matrices
Returns:
matrix Vhf = 2*J - K. Vhf can be a list matrices, corresponding to the input density matrices.

Examples:

>>> import numpy
>>> from pyscf import gto, scf
>>> from pyscf.scf import _vhf
>>> mol = gto.M(atom='H 0 0 0; H 0 0 1.1')
>>> dm0 = numpy.random.random((mol.nao_nr(),mol.nao_nr()))
>>> vhf0 = scf.hf.get_veff(mol, dm0, hermi=0)
>>> dm1 = numpy.random.random((mol.nao_nr(),mol.nao_nr()))
>>> vhf1 = scf.hf.get_veff(mol, dm1, hermi=0)
>>> vhf2 = scf.hf.get_veff(mol, dm1, dm_last=dm0, vhf_last=vhf0, hermi=0)
>>> numpy.allclose(vhf1, vhf2)
True
pyscf.scf.hf.init_guess_by_1e(mol)[source]

Generate initial guess density matrix from core hamiltonian

Returns:
Density matrix, 2D ndarray
pyscf.scf.hf.init_guess_by_atom(mol)[source]

Generate initial guess density matrix from superposition of atomic HF density matrix. The atomic HF is occupancy averaged RHF

Returns:
Density matrix, 2D ndarray
pyscf.scf.hf.init_guess_by_chkfile(mol, chkfile_name, project=None)[source]

Read the HF results from checkpoint file, then project it to the basis defined by mol

Returns:
Density matrix, 2D ndarray
pyscf.scf.hf.init_guess_by_minao(mol)[source]

Generate initial guess density matrix based on ANO basis, then project the density matrix to the basis set defined by mol

Returns:
Density matrix, 2D ndarray

Examples:

>>> from pyscf import gto, scf
>>> mol = gto.M(atom='H 0 0 0; H 0 0 1.1')
>>> scf.hf.init_guess_by_minao(mol)
array([[ 0.94758917,  0.09227308],
       [ 0.09227308,  0.94758917]])
pyscf.scf.hf.kernel(mf, conv_tol=1e-10, conv_tol_grad=None, dump_chk=True, dm0=None, callback=None, conv_check=True, **kwargs)[source]

kernel: the SCF driver.

Args:
mf
: an instance of SCF class

mf object holds all parameters to control SCF. One can modify its member functions to change the behavior of SCF. The member functions which are called in kernel are

mf.get_init_guess
mf.get_hcore
mf.get_ovlp
mf.get_veff
mf.get_fock
mf.get_grad
mf.eig
mf.get_occ
mf.make_rdm1
mf.energy_tot
mf.dump_chk
Kwargs:
conv_tol
: float
converge threshold.
conv_tol_grad
: float
gradients converge threshold.
dump_chk
: bool
Whether to save SCF intermediate results in the checkpoint file
dm0
: ndarray
Initial guess density matrix. If not given (the default), the kernel takes the density matrix generated by mf.get_init_guess.
callback
: function(envs_dict) => None
callback function takes one dict as the argument which is generated by the builtin function locals(), so that the callback function can access all local variables in the current envrionment.
Returns:

A list : scf_conv, e_tot, mo_energy, mo_coeff, mo_occ

scf_conv
: bool
True means SCF converged
e_tot
: float
Hartree-Fock energy of last iteration
mo_energy
: 1D float array
Orbital energies. Depending the eig function provided by mf object, the orbital energies may NOT be sorted.
mo_coeff
: 2D array
Orbital coefficients.
mo_occ
: 1D array
Orbital occupancies. The occupancies may NOT be sorted from large to small.

Examples:

>>> from pyscf import gto, scf
>>> mol = gto.M(atom='H 0 0 0; H 0 0 1.1', basis='cc-pvdz')
>>> conv, e, mo_e, mo, mo_occ = scf.hf.kernel(scf.hf.SCF(mol), dm0=numpy.eye(mol.nao_nr()))
>>> print('conv = %s, E(HF) = %.12f' % (conv, e))
conv = True, E(HF) = -1.081170784378
pyscf.scf.hf.level_shift(s, d, f, factor)[source]

Apply level shift \(\Delta\) to virtual orbitals

\begin{align} FC &= SCE \\ F &= F + SC \Lambda C^\dagger S \\ \Lambda_{ij} &= \begin{cases} \delta_{ij}\Delta & i \in \text{virtual} \\ 0 & \text{otherwise} \end{cases} \end{align}
Returns:
New Fock matrix, 2D ndarray
pyscf.scf.hf.make_rdm1(mo_coeff, mo_occ, **kwargs)[source]

One-particle density matrix in AO representation

Args:
mo_coeff
: 2D ndarray
Orbital coefficients. Each column is one orbital.
mo_occ
: 1D ndarray
Occupancy
pyscf.scf.hf.mulliken_meta(mol, dm, verbose=5, pre_orth_method='ANO', s=None)[source]

Mulliken population analysis, based on meta-Lowdin AOs. In the meta-lowdin, the AOs are grouped in three sets: core, valence and Rydberg, the orthogonalization are carreid out within each subsets.

Args:

mol : an instance of Mole

dm
: ndarray or 2-item list of ndarray
Density matrix. ROHF dm is a 2-item list of 2D array
Kwargs:

verbose : int or instance of lib.logger.Logger

pre_orth_method
: str

Pre-orthogonalization, which localized GTOs for each atom. To obtain the occupied and unoccupied atomic shells, there are three methods

‘ano’ : Project GTOs to ANO basis
‘minao’ : Project GTOs to MINAO basis
‘scf’ : Fraction-averaged RHF
pyscf.scf.hf.mulliken_pop(mol, dm, s=None, verbose=5)[source]

Mulliken population analysis

\[M_{ij} = D_{ij} S_{ji}\]

Mulliken charges

\[\delta_i = \sum_j M_{ij}\]
pyscf.scf.hf.mulliken_pop_meta_lowdin_ao(mol, dm, verbose=5, pre_orth_method='ANO', s=None)

Mulliken population analysis, based on meta-Lowdin AOs. In the meta-lowdin, the AOs are grouped in three sets: core, valence and Rydberg, the orthogonalization are carreid out within each subsets.

Args:

mol : an instance of Mole

dm
: ndarray or 2-item list of ndarray
Density matrix. ROHF dm is a 2-item list of 2D array
Kwargs:

verbose : int or instance of lib.logger.Logger

pre_orth_method
: str

Pre-orthogonalization, which localized GTOs for each atom. To obtain the occupied and unoccupied atomic shells, there are three methods

‘ano’ : Project GTOs to ANO basis
‘minao’ : Project GTOs to MINAO basis
‘scf’ : Fraction-averaged RHF
pyscf.scf.hf.pack_uniq_var(x, mo_occ)[source]

Extract the unique variables from the full orbital-gradients (or orbital-rotation) matrix

pyscf.scf.hf.uniq_var_indices(mo_occ)[source]

Indicies of the unique variables for the orbital-gradients (or orbital-rotation) matrix.

pyscf.scf.hf.unpack_uniq_var(dx, mo_occ)[source]

Fill the full orbital-gradients (or orbital-rotation) matrix with the unique variables.


class pyscf.scf.uhf.UHF(mol)[source]

SCF base class. non-relativistic RHF.

Attributes:
verbose
: int
Print level. Default value equals to Mole.verbose
max_memory
: float or int
Allowed memory in MB. Default equals to Mole.max_memory
chkfile
: str
checkpoint file to save MOs, orbital energies etc. Writing to chkfile can be disabled if this attribute is set to None or False.
conv_tol
: float
converge threshold. Default is 1e-9
conv_tol_grad
: float
gradients converge threshold. Default is sqrt(conv_tol)
max_cycle
: int
max number of iterations. Default is 50
init_guess
: str
initial guess method. It can be one of ‘minao’, ‘atom’, ‘hcore’, ‘1e’, ‘chkfile’. Default is ‘minao’
DIIS
: DIIS class
The class to generate diis object. It can be one of diis.SCF_DIIS, diis.ADIIS, diis.EDIIS.
diis
: boolean or object of DIIS class defined in scf.diis.
Default is the object associated to the attribute self.DIIS. Set it to None/False to turn off DIIS. Note if this attribute is inialized as a DIIS object, the SCF driver will use this object in the iteration. The DIIS informations (vector basis and error vector) will be held inside this object. When kernel function is called again, the old states (vector basis and error vector) will be reused.
diis_space
: int
DIIS space size. By default, 8 Fock matrices and errors vector are stored.
diis_start_cycle
: int
The step to start DIIS. Default is 1.
diis_file: ‘str’
File to store DIIS vectors and error vectors.
level_shift
: float or int
Level shift (in AU) for virtual space. Default is 0.
direct_scf
: bool
Direct SCF is used by default.
direct_scf_tol
: float
Direct SCF cutoff threshold. Default is 1e-13.
callback
: function(envs_dict) => None
callback function takes one dict as the argument which is generated by the builtin function locals(), so that the callback function can access all local variables in the current envrionment.
conv_check
: bool
An extra cycle to check convergence after SCF iterations.
check_convergence
: function(envs) => bool
A hook for overloading convergence criteria in SCF iterations.

Saved results:

converged
: bool
SCF converged or not
e_tot
: float
Total HF energy (electronic energy plus nuclear repulsion)
mo_energy :
Orbital energies
mo_occ
Orbital occupancy
mo_coeff
Orbital coefficients

Examples:

>>> mol = gto.M(atom='H 0 0 0; H 0 0 1.1', basis='cc-pvdz')
>>> mf = scf.hf.SCF(mol)
>>> mf.verbose = 0
>>> mf.level_shift = .4
>>> mf.scf()
-1.0811707843775884
Attributes for UHF:
nelec
: (int, int)
If given, freeze the number of (alpha,beta) electrons to the given value.
level_shift
: number or two-element list
level shift (in Eh) for alpha and beta Fock if two-element list is given.

Examples:

>>> mol = gto.M(atom='O 0 0 0; H 0 0 1; H 0 1 0', basis='ccpvdz', charge=1, spin=1, verbose=0)
>>> mf = scf.UHF(mol)
>>> mf.kernel()
-75.623975516256706
>>> print('S^2 = %.7f, 2S+1 = %.7f' % mf.spin_square())
S^2 = 0.7570150, 2S+1 = 2.0070027
canonicalize(mf, mo_coeff, mo_occ, fock=None)

Canonicalization diagonalizes the UHF Fock matrix within occupied, virtual subspaces separatedly (without change occupancy).

convert_from_(mf)[source]

Create UHF object based on the RHF/ROHF object

det_ovlp(mo1, mo2, occ1, occ2, ovlp=None)[source]

Calculate the overlap between two different determinants. It is the product of single values of molecular orbital overlap matrix.

\[S_{12} = \langle \Psi_A | \Psi_B \rangle = (\mathrm{det}\mathbf{U}) (\mathrm{det}\mathbf{V^\dagger})\prod\limits_{i=1}\limits^{2N} \lambda_{ii}\]

where \(\mathbf{U}, \mathbf{V}, \lambda\) are unitary matrices and single values generated by single value decomposition(SVD) of the overlap matrix \(\mathbf{O}\) which is the overlap matrix of two sets of molecular orbitals:

\[\mathbf{U}^\dagger \mathbf{O} \mathbf{V} = \mathbf{\Lambda}\]
Args:
mo1, mo2
: 2D ndarrays
Molecualr orbital coefficients
occ1, occ2: 2D ndarrays
occupation numbers
Return:
A list: the product of single values: float
x_a, x_b: 1D ndarrays \(\mathbf{U} \mathbf{\Lambda}^{-1} \mathbf{V}^\dagger\) They are used to calculate asymmetric density matrix
energy_elec(mf, dm=None, h1e=None, vhf=None)

Electronic energy of Unrestricted Hartree-Fock

Returns:
Hartree-Fock electronic energy and the 2-electron part contribution
get_jk(mol=None, dm=None, hermi=1)[source]

Coulomb (J) and exchange (K)

Args:
dm
: a list of 2D arrays or a list of 3D arrays
(alpha_dm, beta_dm) or (alpha_dms, beta_dms)
get_veff(mol=None, dm=None, dm_last=0, vhf_last=0, hermi=1)[source]

Unrestricted Hartree-Fock potential matrix of alpha and beta spins, for the given density matrix

\[\begin{split}V_{ij}^\alpha &= \sum_{kl} (ij|kl)(\gamma_{lk}^\alpha+\gamma_{lk}^\beta) - \sum_{kl} (il|kj)\gamma_{lk}^\alpha \\ V_{ij}^\beta &= \sum_{kl} (ij|kl)(\gamma_{lk}^\alpha+\gamma_{lk}^\beta) - \sum_{kl} (il|kj)\gamma_{lk}^\beta\end{split}\]
Args:

mol : an instance of Mole

dm
: a list of ndarrays
A list of density matrices, stored as (alpha,alpha,...,beta,beta,...)
Kwargs:
dm_last
: ndarray or a list of ndarrays or 0
The density matrix baseline. When it is not 0, this function computes the increment of HF potential w.r.t. the reference HF potential matrix.
vhf_last
: ndarray or a list of ndarrays or 0
The reference HF potential matrix.
hermi
: int

Whether J, K matrix is hermitian

0 : no hermitian or symmetric
1 : hermitian
2 : anti-hermitian
vhfopt :
A class which holds precomputed quantities to optimize the computation of J, K matrices
Returns:
\(V_{hf} = (V^\alpha, V^\beta)\). \(V^\alpha\) (and \(V^\beta\)) can be a list matrices, corresponding to the input density matrices.

Examples:

>>> import numpy
>>> from pyscf import gto, scf
>>> mol = gto.M(atom='H 0 0 0; H 0 0 1.1')
>>> dmsa = numpy.random.random((3,mol.nao_nr(),mol.nao_nr()))
>>> dmsb = numpy.random.random((3,mol.nao_nr(),mol.nao_nr()))
>>> dms = numpy.vstack((dmsa,dmsb))
>>> dms.shape
(6, 2, 2)
>>> vhfa, vhfb = scf.uhf.get_veff(mol, dms, hermi=0)
>>> vhfa.shape
(3, 2, 2)
>>> vhfb.shape
(3, 2, 2)
init_guess_by_minao(mol=None, breaksym=True)[source]

Initial guess in terms of the overlap to minimal basis.

make_asym_dm(mo1, mo2, occ1, occ2, x)[source]

One-particle asymmetric density matrix

Args:
mo1, mo2
: 2D ndarrays
Molecualr orbital coefficients
occ1, occ2: 2D ndarrays
Occupation numbers
x: 2D ndarrays
\(\mathbf{U} \mathbf{\Lambda}^{-1} \mathbf{V}^\dagger\). See also det_ovlp()
Return:
A list of 2D ndarrays for alpha and beta spin

Examples:

>>> mf1 = scf.UHF(gto.M(atom='H 0 0 0; F 0 0 1.3', basis='ccpvdz')).run()
>>> mf2 = scf.UHF(gto.M(atom='H 0 0 0; F 0 0 1.4', basis='ccpvdz')).run()
>>> s = gto.intor_cross('int1e_ovlp_sph', mf1.mol, mf2.mol)
>>> det, x = det_ovlp(mf1.mo_coeff, mf1.mo_occ, mf2.mo_coeff, mf2.mo_occ, s)
>>> adm = make_asym_dm(mf1.mo_coeff, mf1.mo_occ, mf2.mo_coeff, mf2.mo_occ, x)
>>> adm.shape
(2, 19, 19)
make_rdm1(mo_coeff=None, mo_occ=None, **kwargs)[source]

One-particle density matrix

Returns:
A list of 2D ndarrays for alpha and beta spins
spin_square(mo_coeff=None, s=None)[source]

Spin square and multiplicity of UHF determinant

\[S^2 = \frac{1}{2}(S_+ S_- + S_- S_+) + S_z^2\]

where \(S_+ = \sum_i S_{i+}\) is effective for all beta occupied orbitals; \(S_- = \sum_i S_{i-}\) is effective for all alpha occupied orbitals.

  1. There are two possibilities for \(S_+ S_-\)
    1. same electron \(S_+ S_- = \sum_i s_{i+} s_{i-}\),
    \[\sum_i \langle UHF|s_{i+} s_{i-}|UHF\rangle = \sum_{pq}\langle p|s_+s_-|q\rangle \gamma_{qp} = n_\alpha\]

    2) different electrons \(S_+ S_- = \sum s_{i+} s_{j-}, (i\neq j)\). There are in total \(n(n-1)\) terms. As a two-particle operator,

    \[\langle S_+ S_- \rangle = \langle ij|s_+ s_-|ij\rangle - \langle ij|s_+ s_-|ji\rangle = -\langle i^\alpha|j^\beta\rangle \langle j^\beta|i^\alpha\rangle\]
  2. Similarly, for \(S_- S_+\)
    1. same electron
    \[\sum_i \langle s_{i-} s_{i+}\rangle = n_\beta\]
    1. different electrons
    \[\langle S_- S_+ \rangle = -\langle i^\beta|j^\alpha\rangle \langle j^\alpha|i^\beta\rangle\]
  3. For \(S_z^2\)
    1. same electron
    \[\langle s_z^2\rangle = \frac{1}{4}(n_\alpha + n_\beta)\]
    1. different electrons
    \[\begin{split}&\frac{1}{2}\sum_{ij}(\langle ij|2s_{z1}s_{z2}|ij\rangle -\langle ij|2s_{z1}s_{z2}|ji\rangle) \\ &=\frac{1}{4}(\langle i^\alpha|i^\alpha\rangle \langle j^\alpha|j^\alpha\rangle - \langle i^\alpha|i^\alpha\rangle \langle j^\beta|j^\beta\rangle - \langle i^\beta|i^\beta\rangle \langle j^\alpha|j^\alpha\rangle + \langle i^\beta|i^\beta\rangle \langle j^\beta|j^\beta\rangle) \\ &-\frac{1}{4}(\langle i^\alpha|j^\alpha\rangle \langle j^\alpha|i^\alpha\rangle + \langle i^\beta|j^\beta\rangle\langle j^\beta|i^\beta\rangle) \\ &=\frac{1}{4}(n_\alpha^2 - n_\alpha n_\beta - n_\beta n_\alpha + n_\beta^2) -\frac{1}{4}(n_\alpha + n_\beta) \\ &=\frac{1}{4}((n_\alpha-n_\beta)^2 - (n_\alpha+n_\beta))\end{split}\]

In total

\[\begin{split}\langle S^2\rangle &= \frac{1}{2} (n_\alpha-\sum_{ij}\langle i^\alpha|j^\beta\rangle \langle j^\beta|i^\alpha\rangle +n_\beta -\sum_{ij}\langle i^\beta|j^\alpha\rangle\langle j^\alpha|i^\beta\rangle) + \frac{1}{4}(n_\alpha-n_\beta)^2 \\\end{split}\]
Args:
mo
: a list of 2 ndarrays
Occupied alpha and occupied beta orbitals
Kwargs:
s
: ndarray
AO overlap
Returns:
A list of two floats. The first is the expectation value of S^2. The second is the corresponding 2S+1

Examples:

>>> mol = gto.M(atom='O 0 0 0; H 0 0 1; H 0 1 0', basis='ccpvdz', charge=1, spin=1, verbose=0)
>>> mf = scf.UHF(mol)
>>> mf.kernel()
-75.623975516256706
>>> mo = (mf.mo_coeff[0][:,mf.mo_occ[0]>0], mf.mo_coeff[1][:,mf.mo_occ[1]>0])
>>> print('S^2 = %.7f, 2S+1 = %.7f' % spin_square(mo, mol.intor('int1e_ovlp_sph')))
S^2 = 0.7570150, 2S+1 = 2.0070027
stability(internal=True, external=False, verbose=None)[source]

Stability analysis for RHF/RKS method.

See also pyscf.scf.stability.uhf_stability function.

Args:
mf : UHF or UKS object
Kwargs:
internal
: bool
Internal stability, within the UHF space.
external
: bool
External stability. Including the UHF -> GHF and real -> complex stability analysis.
Returns:
New orbitals that are more close to the stable condition. The return value includes two set of orbitals. The first corresponds to the internal stablity and the second corresponds to the external stability.
pyscf.scf.uhf.analyze(mf, verbose=5, with_meta_lowdin=True, **kwargs)[source]

Analyze the given SCF object: print orbital energies, occupancies; print orbital coefficients; Mulliken population analysis; Dipole moment

pyscf.scf.uhf.canonicalize(mf, mo_coeff, mo_occ, fock=None)[source]

Canonicalization diagonalizes the UHF Fock matrix within occupied, virtual subspaces separatedly (without change occupancy).

pyscf.scf.uhf.det_ovlp(mo1, mo2, occ1, occ2, ovlp)[source]

Calculate the overlap between two different determinants. It is the product of single values of molecular orbital overlap matrix.

\[S_{12} = \langle \Psi_A | \Psi_B \rangle = (\mathrm{det}\mathbf{U}) (\mathrm{det}\mathbf{V^\dagger})\prod\limits_{i=1}\limits^{2N} \lambda_{ii}\]

where \(\mathbf{U}, \mathbf{V}, \lambda\) are unitary matrices and single values generated by single value decomposition(SVD) of the overlap matrix \(\mathbf{O}\) which is the overlap matrix of two sets of molecular orbitals:

\[\mathbf{U}^\dagger \mathbf{O} \mathbf{V} = \mathbf{\Lambda}\]
Args:
mo1, mo2
: 2D ndarrays
Molecualr orbital coefficients
occ1, occ2: 2D ndarrays
occupation numbers
Return:
A list: the product of single values: float
x_a, x_b: 1D ndarrays \(\mathbf{U} \mathbf{\Lambda}^{-1} \mathbf{V}^\dagger\) They are used to calculate asymmetric density matrix
pyscf.scf.uhf.energy_elec(mf, dm=None, h1e=None, vhf=None)[source]

Electronic energy of Unrestricted Hartree-Fock

Returns:
Hartree-Fock electronic energy and the 2-electron part contribution
pyscf.scf.uhf.get_grad(mo_coeff, mo_occ, fock_ao)[source]

UHF Gradients

pyscf.scf.uhf.get_veff(mol, dm, dm_last=0, vhf_last=0, hermi=1, vhfopt=None)[source]

Unrestricted Hartree-Fock potential matrix of alpha and beta spins, for the given density matrix

\[\begin{split}V_{ij}^\alpha &= \sum_{kl} (ij|kl)(\gamma_{lk}^\alpha+\gamma_{lk}^\beta) - \sum_{kl} (il|kj)\gamma_{lk}^\alpha \\ V_{ij}^\beta &= \sum_{kl} (ij|kl)(\gamma_{lk}^\alpha+\gamma_{lk}^\beta) - \sum_{kl} (il|kj)\gamma_{lk}^\beta\end{split}\]
Args:

mol : an instance of Mole

dm
: a list of ndarrays
A list of density matrices, stored as (alpha,alpha,...,beta,beta,...)
Kwargs:
dm_last
: ndarray or a list of ndarrays or 0
The density matrix baseline. When it is not 0, this function computes the increment of HF potential w.r.t. the reference HF potential matrix.
vhf_last
: ndarray or a list of ndarrays or 0
The reference HF potential matrix.
hermi
: int

Whether J, K matrix is hermitian

0 : no hermitian or symmetric
1 : hermitian
2 : anti-hermitian
vhfopt :
A class which holds precomputed quantities to optimize the computation of J, K matrices
Returns:
\(V_{hf} = (V^\alpha, V^\beta)\). \(V^\alpha\) (and \(V^\beta\)) can be a list matrices, corresponding to the input density matrices.

Examples:

>>> import numpy
>>> from pyscf import gto, scf
>>> mol = gto.M(atom='H 0 0 0; H 0 0 1.1')
>>> dmsa = numpy.random.random((3,mol.nao_nr(),mol.nao_nr()))
>>> dmsb = numpy.random.random((3,mol.nao_nr(),mol.nao_nr()))
>>> dms = numpy.vstack((dmsa,dmsb))
>>> dms.shape
(6, 2, 2)
>>> vhfa, vhfb = scf.uhf.get_veff(mol, dms, hermi=0)
>>> vhfa.shape
(3, 2, 2)
>>> vhfb.shape
(3, 2, 2)
pyscf.scf.uhf.init_guess_by_chkfile(mol, chkfile_name, project=None)[source]

Read SCF chkfile and make the density matrix for UHF initial guess.

Kwargs:
project
: None or bool

Whether to project chkfile’s orbitals to the new basis. Note when the geometry of the chkfile and the given molecule are very different, this projection can produce very poor initial guess. In PES scanning, it is recommended to swith off project.

If project is set to None, the projection is only applied when the basis sets of the chkfile’s molecule are different to the basis sets of the given molecule (regardless whether the geometry of the two molecules are different). Note the basis sets are considered to be different if the two molecules are derived from the same molecule with different ordering of atoms.

pyscf.scf.uhf.init_guess_by_minao(mol, breaksym=True)[source]

Generate initial guess density matrix based on ANO basis, then project the density matrix to the basis set defined by mol

Returns:
Density matrices, a list of 2D ndarrays for alpha and beta spins
pyscf.scf.uhf.make_asym_dm(mo1, mo2, occ1, occ2, x)[source]

One-particle asymmetric density matrix

Args:
mo1, mo2
: 2D ndarrays
Molecualr orbital coefficients
occ1, occ2: 2D ndarrays
Occupation numbers
x: 2D ndarrays
\(\mathbf{U} \mathbf{\Lambda}^{-1} \mathbf{V}^\dagger\). See also det_ovlp()
Return:
A list of 2D ndarrays for alpha and beta spin

Examples:

>>> mf1 = scf.UHF(gto.M(atom='H 0 0 0; F 0 0 1.3', basis='ccpvdz')).run()
>>> mf2 = scf.UHF(gto.M(atom='H 0 0 0; F 0 0 1.4', basis='ccpvdz')).run()
>>> s = gto.intor_cross('int1e_ovlp_sph', mf1.mol, mf2.mol)
>>> det, x = det_ovlp(mf1.mo_coeff, mf1.mo_occ, mf2.mo_coeff, mf2.mo_occ, s)
>>> adm = make_asym_dm(mf1.mo_coeff, mf1.mo_occ, mf2.mo_coeff, mf2.mo_occ, x)
>>> adm.shape
(2, 19, 19)
pyscf.scf.uhf.make_rdm1(mo_coeff, mo_occ, **kwargs)[source]

One-particle density matrix

Returns:
A list of 2D ndarrays for alpha and beta spins
pyscf.scf.uhf.mulliken_meta(mol, dm_ao, verbose=5, pre_orth_method='ANO', s=None)[source]

Mulliken population analysis, based on meta-Lowdin AOs.

pyscf.scf.uhf.mulliken_pop(mol, dm, s=None, verbose=5)[source]

Mulliken population analysis

pyscf.scf.uhf.mulliken_pop_meta_lowdin_ao(mol, dm_ao, verbose=5, pre_orth_method='ANO', s=None)

Mulliken population analysis, based on meta-Lowdin AOs.

pyscf.scf.uhf.spin_square(mo, s=1)[source]

Spin square and multiplicity of UHF determinant

\[S^2 = \frac{1}{2}(S_+ S_- + S_- S_+) + S_z^2\]

where \(S_+ = \sum_i S_{i+}\) is effective for all beta occupied orbitals; \(S_- = \sum_i S_{i-}\) is effective for all alpha occupied orbitals.

  1. There are two possibilities for \(S_+ S_-\)
    1. same electron \(S_+ S_- = \sum_i s_{i+} s_{i-}\),
    \[\sum_i \langle UHF|s_{i+} s_{i-}|UHF\rangle = \sum_{pq}\langle p|s_+s_-|q\rangle \gamma_{qp} = n_\alpha\]

    2) different electrons \(S_+ S_- = \sum s_{i+} s_{j-}, (i\neq j)\). There are in total \(n(n-1)\) terms. As a two-particle operator,

    \[\langle S_+ S_- \rangle = \langle ij|s_+ s_-|ij\rangle - \langle ij|s_+ s_-|ji\rangle = -\langle i^\alpha|j^\beta\rangle \langle j^\beta|i^\alpha\rangle\]
  2. Similarly, for \(S_- S_+\)
    1. same electron
    \[\sum_i \langle s_{i-} s_{i+}\rangle = n_\beta\]
    1. different electrons
    \[\langle S_- S_+ \rangle = -\langle i^\beta|j^\alpha\rangle \langle j^\alpha|i^\beta\rangle\]
  3. For \(S_z^2\)
    1. same electron
    \[\langle s_z^2\rangle = \frac{1}{4}(n_\alpha + n_\beta)\]
    1. different electrons
    \[\begin{split}&\frac{1}{2}\sum_{ij}(\langle ij|2s_{z1}s_{z2}|ij\rangle -\langle ij|2s_{z1}s_{z2}|ji\rangle) \\ &=\frac{1}{4}(\langle i^\alpha|i^\alpha\rangle \langle j^\alpha|j^\alpha\rangle - \langle i^\alpha|i^\alpha\rangle \langle j^\beta|j^\beta\rangle - \langle i^\beta|i^\beta\rangle \langle j^\alpha|j^\alpha\rangle + \langle i^\beta|i^\beta\rangle \langle j^\beta|j^\beta\rangle) \\ &-\frac{1}{4}(\langle i^\alpha|j^\alpha\rangle \langle j^\alpha|i^\alpha\rangle + \langle i^\beta|j^\beta\rangle\langle j^\beta|i^\beta\rangle) \\ &=\frac{1}{4}(n_\alpha^2 - n_\alpha n_\beta - n_\beta n_\alpha + n_\beta^2) -\frac{1}{4}(n_\alpha + n_\beta) \\ &=\frac{1}{4}((n_\alpha-n_\beta)^2 - (n_\alpha+n_\beta))\end{split}\]

In total

\[\begin{split}\langle S^2\rangle &= \frac{1}{2} (n_\alpha-\sum_{ij}\langle i^\alpha|j^\beta\rangle \langle j^\beta|i^\alpha\rangle +n_\beta -\sum_{ij}\langle i^\beta|j^\alpha\rangle\langle j^\alpha|i^\beta\rangle) + \frac{1}{4}(n_\alpha-n_\beta)^2 \\\end{split}\]
Args:
mo
: a list of 2 ndarrays
Occupied alpha and occupied beta orbitals
Kwargs:
s
: ndarray
AO overlap
Returns:
A list of two floats. The first is the expectation value of S^2. The second is the corresponding 2S+1

Examples:

>>> mol = gto.M(atom='O 0 0 0; H 0 0 1; H 0 1 0', basis='ccpvdz', charge=1, spin=1, verbose=0)
>>> mf = scf.UHF(mol)
>>> mf.kernel()
-75.623975516256706
>>> mo = (mf.mo_coeff[0][:,mf.mo_occ[0]>0], mf.mo_coeff[1][:,mf.mo_occ[1]>0])
>>> print('S^2 = %.7f, 2S+1 = %.7f' % spin_square(mo, mol.intor('int1e_ovlp_sph')))
S^2 = 0.7570150, 2S+1 = 2.0070027

Non-relativistic restricted Hartree-Fock with point group symmetry.

The symmetry are not handled in a separate data structure. Note that during the SCF iteration, the orbitals are grouped in terms of symmetry irreps. But the orbitals in the result are sorted based on the orbital energies. Function symm.label_orb_symm can be used to detect the symmetry of the molecular orbitals.

pyscf.scf.hf_symm.RHF

alias of SymAdaptedRHF

pyscf.scf.hf_symm.ROHF

alias of SymAdaptedROHF

class pyscf.scf.hf_symm.SymAdaptedRHF(mol)[source]

SCF base class. non-relativistic RHF.

Attributes:
verbose
: int
Print level. Default value equals to Mole.verbose
max_memory
: float or int
Allowed memory in MB. Default equals to Mole.max_memory
chkfile
: str
checkpoint file to save MOs, orbital energies etc. Writing to chkfile can be disabled if this attribute is set to None or False.
conv_tol
: float
converge threshold. Default is 1e-9
conv_tol_grad
: float
gradients converge threshold. Default is sqrt(conv_tol)
max_cycle
: int
max number of iterations. Default is 50
init_guess
: str
initial guess method. It can be one of ‘minao’, ‘atom’, ‘hcore’, ‘1e’, ‘chkfile’. Default is ‘minao’
DIIS
: DIIS class
The class to generate diis object. It can be one of diis.SCF_DIIS, diis.ADIIS, diis.EDIIS.
diis
: boolean or object of DIIS class defined in scf.diis.
Default is the object associated to the attribute self.DIIS. Set it to None/False to turn off DIIS. Note if this attribute is inialized as a DIIS object, the SCF driver will use this object in the iteration. The DIIS informations (vector basis and error vector) will be held inside this object. When kernel function is called again, the old states (vector basis and error vector) will be reused.
diis_space
: int
DIIS space size. By default, 8 Fock matrices and errors vector are stored.
diis_start_cycle
: int
The step to start DIIS. Default is 1.
diis_file: ‘str’
File to store DIIS vectors and error vectors.
level_shift
: float or int
Level shift (in AU) for virtual space. Default is 0.
direct_scf
: bool
Direct SCF is used by default.
direct_scf_tol
: float
Direct SCF cutoff threshold. Default is 1e-13.
callback
: function(envs_dict) => None
callback function takes one dict as the argument which is generated by the builtin function locals(), so that the callback function can access all local variables in the current envrionment.
conv_check
: bool
An extra cycle to check convergence after SCF iterations.
check_convergence
: function(envs) => bool
A hook for overloading convergence criteria in SCF iterations.

Saved results:

converged
: bool
SCF converged or not
e_tot
: float
Total HF energy (electronic energy plus nuclear repulsion)
mo_energy :
Orbital energies
mo_occ
Orbital occupancy
mo_coeff
Orbital coefficients

Examples:

>>> mol = gto.M(atom='H 0 0 0; H 0 0 1.1', basis='cc-pvdz')
>>> mf = scf.hf.SCF(mol)
>>> mf.verbose = 0
>>> mf.level_shift = .4
>>> mf.scf()
-1.0811707843775884
Attributes for symmetry allowed RHF:
irrep_nelec
: dict
Specify the number of electrons for particular irrep {‘ir_name’:int,...}. For the irreps not listed in this dict, the program will choose the occupancy based on the orbital energies.

Examples:

>>> mol = gto.M(atom='O 0 0 0; H 0 0 1; H 0 1 0', basis='ccpvdz', symmetry=True, verbose=0)
>>> mf = scf.RHF(mol)
>>> mf.scf()
-76.016789472074251
>>> mf.get_irrep_nelec()
{'A1': 6, 'A2': 0, 'B1': 2, 'B2': 2}
>>> mf.irrep_nelec = {'A2': 2}
>>> mf.scf()
-72.768201804695622
>>> mf.get_irrep_nelec()
{'A1': 6, 'A2': 2, 'B1': 2, 'B2': 0}
analyze(verbose=None, with_meta_lowdin=True, **kwargs)[source]

Analyze the given SCF object: print orbital energies, occupancies; print orbital coefficients; Occupancy for each irreps; Mulliken population analysis

canonicalize(mf, mo_coeff, mo_occ, fock=None)

Canonicalization diagonalizes the Fock matrix in occupied, open, virtual subspaces separatedly (without change occupancy).

eig(mf, h, s)

Solve generalized eigenvalue problem, for each irrep. The eigenvalues and eigenvectors are not sorted to ascending order. Instead, they are grouped based on irreps.

get_irrep_nelec(mol=None, mo_coeff=None, mo_occ=None, s=None)[source]

Electron numbers for each irreducible representation.

Args:
mol
: an instance of Mole
To provide irrep_id, and spin-adapted basis
mo_coeff
: 2D ndarray
Regular orbital coefficients, without grouping for irreps
mo_occ
: 1D ndarray
Regular occupancy, without grouping for irreps
Returns:
irrep_nelec
: dict
The number of electrons for each irrep {‘ir_name’:int,...}.

Examples:

>>> mol = gto.M(atom='O 0 0 0; H 0 0 1; H 0 1 0', basis='ccpvdz', symmetry=True, verbose=0)
>>> mf = scf.RHF(mol)
>>> mf.scf()
-76.016789472074251
>>> scf.hf_symm.get_irrep_nelec(mol, mf.mo_coeff, mf.mo_occ)
{'A1': 6, 'A2': 0, 'B1': 2, 'B2': 2}
get_occ(mo_energy=None, mo_coeff=None)[source]

We assumed mo_energy are grouped by symmetry irreps, (see function self.eig). The orbitals are sorted after SCF.

class pyscf.scf.hf_symm.SymAdaptedROHF(mol)[source]

SCF base class. non-relativistic RHF.

Attributes:
verbose
: int
Print level. Default value equals to Mole.verbose
max_memory
: float or int
Allowed memory in MB. Default equals to Mole.max_memory
chkfile
: str
checkpoint file to save MOs, orbital energies etc. Writing to chkfile can be disabled if this attribute is set to None or False.
conv_tol
: float
converge threshold. Default is 1e-9
conv_tol_grad
: float
gradients converge threshold. Default is sqrt(conv_tol)
max_cycle
: int
max number of iterations. Default is 50
init_guess
: str
initial guess method. It can be one of ‘minao’, ‘atom’, ‘hcore’, ‘1e’, ‘chkfile’. Default is ‘minao’
DIIS
: DIIS class
The class to generate diis object. It can be one of diis.SCF_DIIS, diis.ADIIS, diis.EDIIS.
diis
: boolean or object of DIIS class defined in scf.diis.
Default is the object associated to the attribute self.DIIS. Set it to None/False to turn off DIIS. Note if this attribute is inialized as a DIIS object, the SCF driver will use this object in the iteration. The DIIS informations (vector basis and error vector) will be held inside this object. When kernel function is called again, the old states (vector basis and error vector) will be reused.
diis_space
: int
DIIS space size. By default, 8 Fock matrices and errors vector are stored.
diis_start_cycle
: int
The step to start DIIS. Default is 1.
diis_file: ‘str’
File to store DIIS vectors and error vectors.
level_shift
: float or int
Level shift (in AU) for virtual space. Default is 0.
direct_scf
: bool
Direct SCF is used by default.
direct_scf_tol
: float
Direct SCF cutoff threshold. Default is 1e-13.
callback
: function(envs_dict) => None
callback function takes one dict as the argument which is generated by the builtin function locals(), so that the callback function can access all local variables in the current envrionment.
conv_check
: bool
An extra cycle to check convergence after SCF iterations.
check_convergence
: function(envs) => bool
A hook for overloading convergence criteria in SCF iterations.

Saved results:

converged
: bool
SCF converged or not
e_tot
: float
Total HF energy (electronic energy plus nuclear repulsion)
mo_energy :
Orbital energies
mo_occ
Orbital occupancy
mo_coeff
Orbital coefficients

Examples:

>>> mol = gto.M(atom='H 0 0 0; H 0 0 1.1', basis='cc-pvdz')
>>> mf = scf.hf.SCF(mol)
>>> mf.verbose = 0
>>> mf.level_shift = .4
>>> mf.scf()
-1.0811707843775884
Attributes for symmetry allowed ROHF:
irrep_nelec
: dict
Specify the number of alpha/beta electrons for particular irrep {‘ir_name’:(int,int), ...}. For the irreps not listed in these dicts, the program will choose the occupancy based on the orbital energies.

Examples:

>>> mol = gto.M(atom='O 0 0 0; H 0 0 1; H 0 1 0', basis='ccpvdz', symmetry=True, charge=1, spin=1, verbose=0)
>>> mf = scf.RHF(mol)
>>> mf.scf()
-75.619358861084052
>>> mf.get_irrep_nelec()
{'A1': (3, 3), 'A2': (0, 0), 'B1': (1, 1), 'B2': (1, 0)}
>>> mf.irrep_nelec = {'B1': (1, 0)}
>>> mf.scf()
-75.425669486776457
>>> mf.get_irrep_nelec()
{'A1': (3, 3), 'A2': (0, 0), 'B1': (1, 0), 'B2': (1, 1)}
canonicalize(mo_coeff, mo_occ, fock=None)[source]

Canonicalization diagonalizes the Fock matrix in occupied, open, virtual subspaces separatedly (without change occupancy).

eig(fock, s)[source]

Solve generalized eigenvalue problem, for each irrep. The eigenvalues and eigenvectors are not sorted to ascending order. Instead, they are grouped based on irreps.

pyscf.scf.hf_symm.analyze(mf, verbose=5, with_meta_lowdin=True, **kwargs)[source]

Analyze the given SCF object: print orbital energies, occupancies; print orbital coefficients; Occupancy for each irreps; Mulliken population analysis

pyscf.scf.hf_symm.canonicalize(mf, mo_coeff, mo_occ, fock=None)[source]

Canonicalization diagonalizes the Fock matrix in occupied, open, virtual subspaces separatedly (without change occupancy).

pyscf.scf.hf_symm.eig(mf, h, s)[source]

Solve generalized eigenvalue problem, for each irrep. The eigenvalues and eigenvectors are not sorted to ascending order. Instead, they are grouped based on irreps.

pyscf.scf.hf_symm.get_irrep_nelec(mol, mo_coeff, mo_occ, s=None)[source]

Electron numbers for each irreducible representation.

Args:
mol
: an instance of Mole
To provide irrep_id, and spin-adapted basis
mo_coeff
: 2D ndarray
Regular orbital coefficients, without grouping for irreps
mo_occ
: 1D ndarray
Regular occupancy, without grouping for irreps
Returns:
irrep_nelec
: dict
The number of electrons for each irrep {‘ir_name’:int,...}.

Examples:

>>> mol = gto.M(atom='O 0 0 0; H 0 0 1; H 0 1 0', basis='ccpvdz', symmetry=True, verbose=0)
>>> mf = scf.RHF(mol)
>>> mf.scf()
-76.016789472074251
>>> scf.hf_symm.get_irrep_nelec(mol, mf.mo_coeff, mf.mo_occ)
{'A1': 6, 'A2': 0, 'B1': 2, 'B2': 2}
pyscf.scf.hf_symm.so2ao_mo_coeff(so, irrep_mo_coeff)[source]

Transfer the basis of MO coefficients, from spin-adapted basis to AO basis


Non-relativistic unrestricted Hartree-Fock with point group symmetry.

class pyscf.scf.uhf_symm.SymAdaptedUHF(mol)[source]

SCF base class. non-relativistic RHF.

Attributes:
verbose
: int
Print level. Default value equals to Mole.verbose
max_memory
: float or int
Allowed memory in MB. Default equals to Mole.max_memory
chkfile
: str
checkpoint file to save MOs, orbital energies etc. Writing to chkfile can be disabled if this attribute is set to None or False.
conv_tol
: float
converge threshold. Default is 1e-9
conv_tol_grad
: float
gradients converge threshold. Default is sqrt(conv_tol)
max_cycle
: int
max number of iterations. Default is 50
init_guess
: str
initial guess method. It can be one of ‘minao’, ‘atom’, ‘hcore’, ‘1e’, ‘chkfile’. Default is ‘minao’
DIIS
: DIIS class
The class to generate diis object. It can be one of diis.SCF_DIIS, diis.ADIIS, diis.EDIIS.
diis
: boolean or object of DIIS class defined in scf.diis.
Default is the object associated to the attribute self.DIIS. Set it to None/False to turn off DIIS. Note if this attribute is inialized as a DIIS object, the SCF driver will use this object in the iteration. The DIIS informations (vector basis and error vector) will be held inside this object. When kernel function is called again, the old states (vector basis and error vector) will be reused.
diis_space
: int
DIIS space size. By default, 8 Fock matrices and errors vector are stored.
diis_start_cycle
: int
The step to start DIIS. Default is 1.
diis_file: ‘str’
File to store DIIS vectors and error vectors.
level_shift
: float or int
Level shift (in AU) for virtual space. Default is 0.
direct_scf
: bool
Direct SCF is used by default.
direct_scf_tol
: float
Direct SCF cutoff threshold. Default is 1e-13.
callback
: function(envs_dict) => None
callback function takes one dict as the argument which is generated by the builtin function locals(), so that the callback function can access all local variables in the current envrionment.
conv_check
: bool
An extra cycle to check convergence after SCF iterations.
check_convergence
: function(envs) => bool
A hook for overloading convergence criteria in SCF iterations.

Saved results:

converged
: bool
SCF converged or not
e_tot
: float
Total HF energy (electronic energy plus nuclear repulsion)
mo_energy :
Orbital energies
mo_occ
Orbital occupancy
mo_coeff
Orbital coefficients

Examples:

>>> mol = gto.M(atom='H 0 0 0; H 0 0 1.1', basis='cc-pvdz')
>>> mf = scf.hf.SCF(mol)
>>> mf.verbose = 0
>>> mf.level_shift = .4
>>> mf.scf()
-1.0811707843775884
Attributes for UHF:
nelec
: (int, int)
If given, freeze the number of (alpha,beta) electrons to the given value.
level_shift
: number or two-element list
level shift (in Eh) for alpha and beta Fock if two-element list is given.

Examples:

>>> mol = gto.M(atom='O 0 0 0; H 0 0 1; H 0 1 0', basis='ccpvdz', charge=1, spin=1, verbose=0)
>>> mf = scf.UHF(mol)
>>> mf.kernel()
-75.623975516256706
>>> print('S^2 = %.7f, 2S+1 = %.7f' % mf.spin_square())
S^2 = 0.7570150, 2S+1 = 2.0070027
Attributes for symmetry allowed UHF:
irrep_nelec
: dict
Specify the number of alpha/beta electrons for particular irrep {‘ir_name’:(int,int), ...}. For the irreps not listed in these dicts, the program will choose the occupancy based on the orbital energies.

Examples:

>>> mol = gto.M(atom='O 0 0 0; H 0 0 1; H 0 1 0', basis='ccpvdz', symmetry=True, charge=1, spin=1, verbose=0)
>>> mf = scf.RHF(mol)
>>> mf.scf()
-75.623975516256692
>>> mf.get_irrep_nelec()
{'A1': (3, 3), 'A2': (0, 0), 'B1': (1, 1), 'B2': (1, 0)}
>>> mf.irrep_nelec = {'B1': (1, 0)}
>>> mf.scf()
-75.429189192031131
>>> mf.get_irrep_nelec()
{'A1': (3, 3), 'A2': (0, 0), 'B1': (1, 0), 'B2': (1, 1)}
canonicalize(mf, mo_coeff, mo_occ, fock=None)

Canonicalization diagonalizes the UHF Fock matrix in occupied, virtual subspaces separatedly (without change occupancy).

get_irrep_nelec(mol=None, mo_coeff=None, mo_occ=None, s=None)[source]

Alpha/beta electron numbers for each irreducible representation.

Args:
mol
: an instance of Mole
To provide irrep_id, and spin-adapted basis
mo_occ
: a list of 1D ndarray
Regular occupancy, without grouping for irreps
mo_coeff
: a list of 2D ndarray
Regular orbital coefficients, without grouping for irreps
Returns:
irrep_nelec
: dict
The number of alpha/beta electrons for each irrep {‘ir_name’:(int,int), ...}.

Examples:

>>> mol = gto.M(atom='O 0 0 0; H 0 0 1; H 0 1 0', basis='ccpvdz', symmetry=True, charge=1, spin=1, verbose=0)
>>> mf = scf.UHF(mol)
>>> mf.scf()
-75.623975516256721
>>> scf.uhf_symm.get_irrep_nelec(mol, mf.mo_coeff, mf.mo_occ)
{'A1': (3, 3), 'A2': (0, 0), 'B1': (1, 1), 'B2': (1, 0)}
get_occ(mo_energy=None, mo_coeff=None)[source]

We assumed mo_energy are grouped by symmetry irreps, (see function self.eig). The orbitals are sorted after SCF.

pyscf.scf.uhf_symm.UHF

alias of SymAdaptedUHF

pyscf.scf.uhf_symm.canonicalize(mf, mo_coeff, mo_occ, fock=None)[source]

Canonicalization diagonalizes the UHF Fock matrix in occupied, virtual subspaces separatedly (without change occupancy).

pyscf.scf.uhf_symm.get_irrep_nelec(mol, mo_coeff, mo_occ, s=None)[source]

Alpha/beta electron numbers for each irreducible representation.

Args:
mol
: an instance of Mole
To provide irrep_id, and spin-adapted basis
mo_occ
: a list of 1D ndarray
Regular occupancy, without grouping for irreps
mo_coeff
: a list of 2D ndarray
Regular orbital coefficients, without grouping for irreps
Returns:
irrep_nelec
: dict
The number of alpha/beta electrons for each irrep {‘ir_name’:(int,int), ...}.

Examples:

>>> mol = gto.M(atom='O 0 0 0; H 0 0 1; H 0 1 0', basis='ccpvdz', symmetry=True, charge=1, spin=1, verbose=0)
>>> mf = scf.UHF(mol)
>>> mf.scf()
-75.623975516256721
>>> scf.uhf_symm.get_irrep_nelec(mol, mf.mo_coeff, mf.mo_occ)
{'A1': (3, 3), 'A2': (0, 0), 'B1': (1, 1), 'B2': (1, 0)}

6.5.3. Relativistic Hartree-Fock

Dirac Hartree-Fock

pyscf.scf.dhf.DHF

alias of UHF

class pyscf.scf.dhf.RHF(mol)[source]

Dirac-RHF

make_rdm1(mo_coeff=None, mo_occ=None, **kwargs)[source]

D/2 = psi_i^dagpsi_i = psi_{Ti}^dagpsi_{Ti} D(UHF) = psi_i^dagpsi_i + psi_{Ti}^dagpsi_{Ti} RHF average the density of spin up and spin down: D(RHF) = (D(UHF) + T[D(UHF)])/2

class pyscf.scf.dhf.UHF(mol)[source]

SCF base class. non-relativistic RHF.

Attributes:
verbose
: int
Print level. Default value equals to Mole.verbose
max_memory
: float or int
Allowed memory in MB. Default equals to Mole.max_memory
chkfile
: str
checkpoint file to save MOs, orbital energies etc. Writing to chkfile can be disabled if this attribute is set to None or False.
conv_tol
: float
converge threshold. Default is 1e-9
conv_tol_grad
: float
gradients converge threshold. Default is sqrt(conv_tol)
max_cycle
: int
max number of iterations. Default is 50
init_guess
: str
initial guess method. It can be one of ‘minao’, ‘atom’, ‘hcore’, ‘1e’, ‘chkfile’. Default is ‘minao’
DIIS
: DIIS class
The class to generate diis object. It can be one of diis.SCF_DIIS, diis.ADIIS, diis.EDIIS.
diis
: boolean or object of DIIS class defined in scf.diis.
Default is the object associated to the attribute self.DIIS. Set it to None/False to turn off DIIS. Note if this attribute is inialized as a DIIS object, the SCF driver will use this object in the iteration. The DIIS informations (vector basis and error vector) will be held inside this object. When kernel function is called again, the old states (vector basis and error vector) will be reused.
diis_space
: int
DIIS space size. By default, 8 Fock matrices and errors vector are stored.
diis_start_cycle
: int
The step to start DIIS. Default is 1.
diis_file: ‘str’
File to store DIIS vectors and error vectors.
level_shift
: float or int
Level shift (in AU) for virtual space. Default is 0.
direct_scf
: bool
Direct SCF is used by default.
direct_scf_tol
: float
Direct SCF cutoff threshold. Default is 1e-13.
callback
: function(envs_dict) => None
callback function takes one dict as the argument which is generated by the builtin function locals(), so that the callback function can access all local variables in the current envrionment.
conv_check
: bool
An extra cycle to check convergence after SCF iterations.
check_convergence
: function(envs) => bool
A hook for overloading convergence criteria in SCF iterations.

Saved results:

converged
: bool
SCF converged or not
e_tot
: float
Total HF energy (electronic energy plus nuclear repulsion)
mo_energy :
Orbital energies
mo_occ
Orbital occupancy
mo_coeff
Orbital coefficients

Examples:

>>> mol = gto.M(atom='H 0 0 0; H 0 0 1.1', basis='cc-pvdz')
>>> mf = scf.hf.SCF(mol)
>>> mf.verbose = 0
>>> mf.level_shift = .4
>>> mf.scf()
-1.0811707843775884
Attributes for Dirac-Hartree-Fock
with_ssss
: bool, for Dirac-Hartree-Fock only
If False, ignore small component integrals (SS|SS). Default is True.
with_gaunt
: bool, for Dirac-Hartree-Fock only
Default is False.
with_breit
: bool, for Dirac-Hartree-Fock only
Gaunt + gauge term. Default is False.

Examples:

>>> mol = gto.M(atom='H 0 0 0; H 0 0 1', basis='ccpvdz', verbose=0)
>>> mf = scf.RHF(mol)
>>> e0 = mf.scf()
>>> mf = scf.DHF(mol)
>>> e1 = mf.scf()
>>> print('Relativistic effects = %.12f' % (e1-e0))
Relativistic effects = -0.000008854205
dip_moment(mol=None, dm=None, unit='Debye', verbose=3, **kwargs)[source]

Dipole moment calculation

\[\begin{split}\mu_x = -\sum_{\mu}\sum_{\nu} P_{\mu\nu}(\nu|x|\mu) + \sum_A Q_A X_A\\ \mu_y = -\sum_{\mu}\sum_{\nu} P_{\mu\nu}(\nu|y|\mu) + \sum_A Q_A Y_A\\ \mu_z = -\sum_{\mu}\sum_{\nu} P_{\mu\nu}(\nu|z|\mu) + \sum_A Q_A Z_A\end{split}\]

where \(\mu_x, \mu_y, \mu_z\) are the x, y and z components of dipole moment

Args:
mol: an instance of Mole dm : a 2D ndarrays density matrices
Return:
A list: the dipole moment on x, y and z component
get_veff(mol=None, dm=None, dm_last=0, vhf_last=0, hermi=1)[source]

Dirac-Coulomb

init_guess_by_minao(mol=None)[source]

Initial guess in terms of the overlap to minimal basis.

mulliken_pop(mol=None, dm=None, s=None, verbose=5)[source]

Mulliken population analysis

\[M_{ij} = D_{ij} S_{ji}\]

Mulliken charges

\[\delta_i = \sum_j M_{ij}\]
pyscf.scf.dhf.dip_moment(mol, dm, unit='Debye', verbose=3, **kwargs)[source]

Dipole moment calculation

\[\begin{split}\mu_x = -\sum_{\mu}\sum_{\nu} P_{\mu\nu}(\nu|x|\mu) + \sum_A Q_A X_A\\ \mu_y = -\sum_{\mu}\sum_{\nu} P_{\mu\nu}(\nu|y|\mu) + \sum_A Q_A Y_A\\ \mu_z = -\sum_{\mu}\sum_{\nu} P_{\mu\nu}(\nu|z|\mu) + \sum_A Q_A Z_A\end{split}\]

where \(\mu_x, \mu_y, \mu_z\) are the x, y and z components of dipole moment

Args:
mol: an instance of Mole dm : a 2D ndarrays density matrices
Return:
A list: the dipole moment on x, y and z component
pyscf.scf.dhf.get_grad(mo_coeff, mo_occ, fock_ao)[source]

DHF Gradients

pyscf.scf.dhf.get_init_guess(mol, key='minao')[source]

Generate density matrix for initial guess

Kwargs:
key
: str
One of ‘minao’, ‘atom’, ‘hcore’, ‘1e’, ‘chkfile’.
pyscf.scf.dhf.init_guess_by_1e(mol)[source]

Initial guess from one electron system.

pyscf.scf.dhf.init_guess_by_atom(mol)[source]

Initial guess from atom calculation.

pyscf.scf.dhf.init_guess_by_chkfile(mol, chkfile_name, project=None)[source]

Read SCF chkfile and make the density matrix for 4C-DHF initial guess.

Kwargs:
project
: None or bool

Whether to project chkfile’s orbitals to the new basis. Note when the geometry of the chkfile and the given molecule are very different, this projection can produce very poor initial guess. In PES scanning, it is recommended to swith off project.

If project is set to None, the projection is only applied when the basis sets of the chkfile’s molecule are different to the basis sets of the given molecule (regardless whether the geometry of the two molecules are different). Note the basis sets are considered to be different if the two molecules are derived from the same molecule with different ordering of atoms.

pyscf.scf.dhf.init_guess_by_minao(mol)[source]

Initial guess in terms of the overlap to minimal basis.

pyscf.scf.dhf.kernel(mf, conv_tol=1e-09, conv_tol_grad=None, dump_chk=True, dm0=None, callback=None, conv_check=True)[source]

the modified SCF kernel for Dirac-Hartree-Fock. In this kernel, the SCF is carried out in three steps. First the 2-electron part is approximated by large component integrals (LL|LL); Next, (SS|LL) the interaction between large and small components are added; Finally, converge the SCF with the small component contributions (SS|SS)

pyscf.scf.dhf.mulliken_pop(mol, dm, s=None, verbose=5)[source]

Mulliken population analysis

\[M_{ij} = D_{ij} S_{ji}\]

Mulliken charges

\[\delta_i = \sum_j M_{ij}\]
pyscf.scf.dhf.time_reversal_matrix(mol, mat)[source]

T(A_ij) = A[T(i),T(j)]^*

6.5.4. addons

pyscf.scf.addons.convert_to_ghf(mf, out=None, remove_df=False)[source]

Convert the given mean-field object to the generalized HF/KS object

Note if mf is an second order SCF object, the second order object will not be converted (in other words, only the underlying SCF object will be converted)

Args:
mf : SCF object
Kwargs
remove_df
: bool
Whether to convert the DF-SCF object to the normal SCF object. This conversion is not applied by default.
Returns:
An generalized SCF object
pyscf.scf.addons.convert_to_rhf(mf, out=None, remove_df=False)[source]

Convert the given mean-field object to the restricted HF/KS object

Note if mf is an second order SCF object, the second order object will not be converted (in other words, only the underlying SCF object will be converted)

Args:
mf : SCF object
Kwargs
remove_df
: bool
Whether to convert the DF-SCF object to the normal SCF object. This conversion is not applied by default.
Returns:
An unrestricted SCF object
pyscf.scf.addons.convert_to_uhf(mf, out=None, remove_df=False)[source]

Convert the given mean-field object to the unrestricted HF/KS object

Note if mf is an second order SCF object, the second order object will not be converted (in other words, only the underlying SCF object will be converted)

Args:
mf : SCF object
Kwargs
remove_df
: bool
Whether to convert the DF-SCF object to the normal SCF object. This conversion is not applied by default.
Returns:
An unrestricted SCF object
pyscf.scf.addons.dynamic_level_shift(mf, factor=1.0)

Dynamically change the level shift in each SCF cycle. The level shift value is set to (HF energy change * factor)

pyscf.scf.addons.dynamic_level_shift_(mf, factor=1.0)[source]

Dynamically change the level shift in each SCF cycle. The level shift value is set to (HF energy change * factor)

pyscf.scf.addons.dynamic_sz_(mf)

For UHF, allowing the Sz value being changed during SCF iteration. Determine occupation of alpha and beta electrons based on energy spectrum

pyscf.scf.addons.float_occ(mf)

For UHF, allowing the Sz value being changed during SCF iteration. Determine occupation of alpha and beta electrons based on energy spectrum

pyscf.scf.addons.float_occ_(mf)[source]

For UHF, allowing the Sz value being changed during SCF iteration. Determine occupation of alpha and beta electrons based on energy spectrum

pyscf.scf.addons.get_ghf_orbspin(mo_energy, mo_occ, is_rhf=None)[source]

Spin of each GHF orbital when the GHF orbitals are converted from RHF/UHF orbitals

For RHF orbitals, the orbspin corresponds to first occupied orbitals then unoccupied orbitals. In the occupied orbital space, if degenerated, first alpha then beta, last the (open-shell) singly occupied (alpha) orbitals. In the unoccupied orbital space, first the (open-shell) unoccupied (beta) orbitals if applicable, then alpha and beta orbitals

For UHF orbitals, the orbspin corresponds to first occupied orbitals then unoccupied orbitals.

pyscf.scf.addons.mom_occ(mf, occorb, setocc)

Use maximum overlap method to determine occupation number for each orbital in every iteration. It can be applied to unrestricted HF/KS and restricted open-shell HF/KS.

pyscf.scf.addons.mom_occ_(mf, occorb, setocc)[source]

Use maximum overlap method to determine occupation number for each orbital in every iteration. It can be applied to unrestricted HF/KS and restricted open-shell HF/KS.

pyscf.scf.addons.project_dm_nr2nr(mol1, dm1, mol2)[source]

Project density matrix representation from basis set 1 (mol1) to basis set 2 (mol2).

\[ \begin{align}\begin{aligned}|AO2\rangle DM_AO2 \langle AO2|\\= |AO2\rangle P DM_AO1 P \langle AO2|\\DM_AO2 = P DM_AO1 P\\P = S_{AO2}^{-1}\langle AO2|AO1\rangle\end{aligned}\end{align} \]

There are three relevant functions: project_dm_nr2nr() is the projection for non-relativistic (scalar) basis. project_dm_nr2r() projects from non-relativistic to relativistic basis. project_dm_r2r() is the projection between relativistic (spinor) basis.

pyscf.scf.addons.project_mo_nr2nr(mol1, mo1, mol2)[source]

Project orbital coefficients from basis set 1 (C1 for mol1) to basis set 2 (C2 for mol2).

\[ \begin{align}\begin{aligned}|\psi1\rangle = |AO1\rangle C1\\|\psi2\rangle = P |\psi1\rangle = |AO2\rangle S^{-1}\langle AO2| AO1\rangle> C1 = |AO2\rangle> C2\\C2 = S^{-1}\langle AO2|AO1\rangle C1\end{aligned}\end{align} \]

There are three relevant functions: project_mo_nr2nr() is the projection for non-relativistic (scalar) basis. project_mo_nr2r() projects from non-relativistic to relativistic basis. project_mo_r2r() is the projection between relativistic (spinor) basis.

pyscf.scf.addons.remove_linear_dep(mf, threshold=1e-08, lindep=1e-10)
Args:
threshold
: float
The threshold under which the eigenvalues of the overlap matrix are discarded to avoid numerical instability.
lindep
: float
The threshold that triggers the special treatment of the linear dependence issue.
pyscf.scf.addons.remove_linear_dep_(mf, threshold=1e-08, lindep=1e-10)[source]
Args:
threshold
: float
The threshold under which the eigenvalues of the overlap matrix are discarded to avoid numerical instability.
lindep
: float
The threshold that triggers the special treatment of the linear dependence issue.

pyscf.scf.chkfile.dump_scf(mol, chkfile, e_tot, mo_energy, mo_coeff, mo_occ, overwrite_mol=True)[source]

save temporary results


DIIS

class pyscf.scf.diis.ADIIS(dev=None, filename=None, incore=False)[source]

Ref: JCP, 132, 054109

class pyscf.scf.diis.EDIIS(dev=None, filename=None, incore=False)[source]

SCF-EDIIS Ref: JCP 116, 8255

pyscf.scf.diis.get_err_vec(s, d, f)[source]

error vector = SDF - FDS