21.2. pbc.scf — Mean-field with periodic boundary condition

This module is an analogy to molecular pyscf.scf module to handle mean-filed calculation with periodic boundary condition.

21.2.1. Gamma point and single k-point calculation

The usage of gamma point Hartree-Fock program is very close to that of the molecular program. In the PBC gamma point calculation, one needs intialize Cell object and the corresponding pyscf.pbc.scf.hf.RHF class:

from pyscf.pbc import gto, scf
cell = gto.M(
    atom = '''H     0.      0.      0.
              H     0.8917  0.8917  0.8917''',
    basis = 'sto3g',
    h = '''
    0       1.7834  1.7834
    1.7834  0       1.7834
    1.7834  1.7834  0     ''',
    gs = [10]*3,
    verbose = 4,
)
mf = scf.RHF(cell).run()

Comparing to the pyscf.scf.hf.RHF object for molecular calculation, the PBC-HF calculation with pyscf.pbc.scf.hf.RHF or pyscf.pbc.scf.uhf.UHF has three differences

  • psycf.pbc.scf.hf.RHF is the single k-point PBC HF class. By default, it creates the gamma point calculation. You can change to other k-point by setting the kpt attribute:

    mf = scf.RHF(cell)
    mf.kpt = cell.get_abs_kpts([.25,.25,.25])  # convert from scaled kpts
    mf.kernel()
    
  • The exchange integrals of the PBC Hartree-Fock method has slow convergence with respect to the number of k-points. Proper treatments for the divergent part of exchange integrals can effectively improve the convergence. Attribute exxdive is used to control the method to handle exchange divergent term. The default exxdiv='ewald' is favored in most scenario. However, if the molecular post-HF methods was mixed with the gamma point HF method (see Mixing with molecular program for post-HF methods, you might need set exxdiv=None to get consistent total energy (see Exchange divergence treatment).

  • In the finite-size system, one can obtain right answer without considering the model to evaluate 2-electron integrals. But the integral scheme might need to be updated in the PBC calculations. The default integral scheme is accurate for pseudo-potential. In the all-electron calculation, you may need change the with_df attribute to mixed density fitting (MDF) method for better accuracy (see with_df for density fitting). Here is an example to update with_df

#!/usr/bin/env python

'''
Gamma point Hartree-Fock/DFT for all-electron calculation

The default FFT-based 2-electron integrals may not be accurate enough for
all-electron calculation.  It's recommended to use MDF (mixed density fitting)
technique to improve the accuracy.

See also
examples/df/00-with_df.py
examples/df/01-auxbasis.py
examples/df/40-precomupte_df_ints.py
'''

import numpy
from pyscf.pbc import gto, scf, dft

cell = gto.M(
    a = numpy.eye(3)*3.5668,
    atom = '''C     0.      0.      0.    
              C     0.8917  0.8917  0.8917
              C     1.7834  1.7834  0.    
              C     2.6751  2.6751  0.8917
              C     1.7834  0.      1.7834
              C     2.6751  0.8917  2.6751
              C     0.      1.7834  1.7834
              C     0.8917  2.6751  2.6751''',
    basis = '6-31g',
    verbose = 4,
)

mf = scf.RHF(cell).density_fit()
mf.kernel()

# Mixed density fitting is another option for all-electron calculations
mf = scf.RHF(cell).mix_density_fit()
mf.with_df.mesh = [10]*3  # Tune #PWs in MDF for performance/accuracy balance
mf.kernel()

# Or use even-tempered Gaussian basis as auxiliary fitting functions.
# The following auxbasis is generated based on the expression
#    alpha = a * 1.7^i   i = 0..N
# where a and N are determined by the smallest and largest exponents of AO basis.
import pyscf.df
auxbasis = pyscf.df.aug_etb(cell, beta=1.7)
mf = scf.RHF(cell).density_fit(auxbasis=auxbasis)
mf.kernel()

#
# Second order SCF solver can be used in the PBC SCF code the same way in the
# molecular calculation
#
mf = dft.RKS(cell).density_fit(auxbasis='weigend')
mf.xc = 'bp86'
# You should first set mf.xc then apply newton method (see also
# examples/scf/22-newton.py)
mf = mf.newton()
mf.kernel()

#
# The computational costs to initialize PBC DF object is high.  The density
# fitting integral tensor created in the initialization can be cached for
# future use.  See also examples/df/40-precomupte_df_ints.py
#
mf = dft.RKS(cell).density_fit(auxbasis='weigend')
mf.with_df._cderi_to_save = 'df_ints.h5'
mf.kernel()
#
# The DF integral tensor can be preloaded in an independent calculation.
#
mf = dft.RKS(cell).density_fit()
mf.with_df._cderi = 'df_ints.h5'
mf.kernel()

21.2.1.1. Mixing with molecular program for post-HF methods

The gamma point HF code adopts the same code structure, the function and method names and the arguments’ convention as the molecule SCF code. This desgin allows one mixing PBC HF object with the existed molecular post-HF code for PBC electron correlation problems. A typical molecular post-HF calculation starts from the finite-size HF method with the Mole object:

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

from pyscf import cc
cc.CCSD(mf).run()

The PBC gamma point post-HF calculation requires the Cell object and PBC HF object:

from pyscf.pbc import gto, scf
cell = gto.M(atom='H 0 0 0; H 0 0 1', basis='ccpvdz',
             h=numpy.eye(3)*2, gs=[10,10,10])
mf = scf.RHF(cell).run()

from pyscf import cc
cc.CCSD(mf).run()

The differences are the the mol or cell object to create and the scf module to import. With the system-specific mean-field object, one can carray out various post-HF methods (MP2, Coupled cluster, CISD, TDHF, TDDFT, ...) using the same code for finite-size and extended systems. See Mixing PBC and molecular modules for more details of the interface between PBC and molecular modules.

21.2.2. k-point sampling

21.2.2.1. Newton solver

21.2.2.2. Smearing

21.2.3. Exchange divergence treatment

The PBC Hartree-Fock has slow convergence of exchange integral with respect to the number of k-points. In the single k-point calculation, Generally, exxdiv leads to a shift in the total energy and the spectrum of orbital energy. It should not affect the following correlation energy in the post-HF calculation. In practice, when the gamma-point calculation is mixed with molecular program eg the FCI solver, the exxdiv attribute may leads to inconsistency in the total energy.

21.2.4. with_df for density fitting

Placing the with_df attribute in SCF object to get the compatibility to molecule DF-SCF methods.

21.2.5. Stability analysis

21.2.6. Program reference

Hartree-Fock for periodic systems at a single k-point

See Also:
pyscf.pbc.scf.khf.py : Hartree-Fock for periodic systems with k-point sampling
class pyscf.pbc.scf.hf.SCF(cell, kpt=array([ 0., 0., 0.]), exxdiv='ewald')[source]

SCF base class adapted for PBCs.

Attributes:
kpt
: (3,) ndarray
The AO k-point in Cartesian coordinates, in units of 1/Bohr.
exxdiv
: str

Exchange divergence treatment, can be one of

None : ignore G=0 contribution in exchange
‘ewald’ : Ewald probe charge correction (JCP, 122, 234102)
with_df
: density fitting object
Default is the FFT based DF model. For all-electron calculation, MDF model is favored for better accuracy. See also pyscf.pbc.df.
direct_scf
: bool
When this flag is set to true, the J/K matrices will be computed directly through the underlying with_df methods. Otherwise, depending the available memory, the 4-index integrals may be cached and J/K matrices are computed based on the 4-index integrals.
dip_moment(cell=None, dm=None, unit='Debye', verbose=3, **kwargs)[source]

Dipole moment in the unit cell.

Args:

cell : an instance of Cell

dm (ndarray) : density matrix

Return:
A list: the dipole moment on x, y and z components
get_bands(mf, kpts_band, cell=None, dm=None, kpt=None)

Get energy bands at the given (arbitrary) ‘band’ k-points.

Returns:
mo_energy
: (nmo,) ndarray or a list of (nmo,) ndarray
Bands energies E_n(k)
mo_coeff
: (nao, nmo) ndarray or a list of (nao,nmo) ndarray
Band orbitals psi_n(k)
get_j(cell=None, dm=None, hermi=1, kpt=None, kpts_band=None)[source]

Compute J matrix for the given density matrix and k-point (kpt). When kpts_band is given, the J matrices on kpts_band are evaluated.

J_{pq} = sum_{rs} (pq|rs) dm[s,r]

where r,s are orbitals on kpt. p and q are orbitals on kpts_band if kpts_band is given otherwise p and q are orbitals on kpt.

get_jk(cell=None, dm=None, hermi=1, kpt=None, kpts_band=None)[source]

Get Coulomb (J) and exchange (K) following scf.hf.RHF.get_jk_(). for particular k-point (kpt).

When kpts_band is given, the J, K matrices on kpts_band are evaluated.

J_{pq} = sum_{rs} (pq|rs) dm[s,r] K_{pq} = sum_{rs} (pr|sq) dm[r,s]

where r,s are orbitals on kpt. p and q are orbitals on kpts_band if kpts_band is given otherwise p and q are orbitals on kpt.

get_jk_incore(cell=None, dm=None, hermi=1, verbose=5, kpt=None)[source]

Get Coulomb (J) and exchange (K) following scf.hf.RHF.get_jk_().

Incore version of Coulomb and exchange build only. Currently RHF always uses PBC AO integrals (unlike RKS), since exchange is currently computed by building PBC AO integrals.

get_k(cell=None, dm=None, hermi=1, kpt=None, kpts_band=None)[source]

Compute K matrix for the given density matrix.

get_rho(mf, dm=None, grids=None, kpt=None)

Compute density in real space

get_veff(cell=None, dm=None, dm_last=0, vhf_last=0, hermi=1, kpt=None, kpts_band=None)[source]

Hartree-Fock potential matrix for the given density matrix. See scf.hf.get_veff() and scf.hf.RHF.get_veff()

pyscf.pbc.scf.hf.dip_moment(cell, dm, unit='Debye', verbose=3, grids=None, rho=None, kpt=array([ 0., 0., 0.]), origin=None)[source]

Dipole moment in the unit cell.

Args:

cell : an instance of Cell

dm (ndarray) : density matrix

Return:
A list: the dipole moment on x, y and z components
pyscf.pbc.scf.hf.get_bands(mf, kpts_band, cell=None, dm=None, kpt=None)[source]

Get energy bands at the given (arbitrary) ‘band’ k-points.

Returns:
mo_energy
: (nmo,) ndarray or a list of (nmo,) ndarray
Bands energies E_n(k)
mo_coeff
: (nao, nmo) ndarray or a list of (nao,nmo) ndarray
Band orbitals psi_n(k)
pyscf.pbc.scf.hf.get_hcore(cell, kpt=array([ 0., 0., 0.]))[source]

Get the core Hamiltonian AO matrix.

pyscf.pbc.scf.hf.get_j(cell, dm, hermi=1, vhfopt=None, kpt=array([ 0., 0., 0.]), kpts_band=None)[source]

Get the Coulomb (J) AO matrix for the given density matrix.

Args:
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
vhfopt :
A class which holds precomputed quantities to optimize the computation of J, K matrices
kpt
: (3,) ndarray
The “inner” dummy k-point at which the DM was evaluated (or sampled).
kpts_band
: (3,) ndarray or (*,3) ndarray
An arbitrary “band” k-point at which J is evaluated.
Returns:
The function returns one J matrix, corresponding to the input density matrix (both order and shape).
pyscf.pbc.scf.hf.get_jk(mf, cell, dm, hermi=1, vhfopt=None, kpt=array([ 0., 0., 0.]), kpts_band=None)[source]

Get the Coulomb (J) and exchange (K) AO matrices for the given density matrix.

Args:
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
vhfopt :
A class which holds precomputed quantities to optimize the computation of J, K matrices
kpt
: (3,) ndarray
The “inner” dummy k-point at which the DM was evaluated (or sampled).
kpts_band
: (3,) ndarray or (*,3) ndarray
An arbitrary “band” k-point at which J and K are evaluated.
Returns:
The function returns one J and one K matrix, corresponding to the input density matrix (both order and shape).
pyscf.pbc.scf.hf.get_nuc(cell, kpt=array([ 0., 0., 0.]))[source]

Get the bare periodic nuc-el AO matrix, with G=0 removed.

See Martin (12.16)-(12.21).

pyscf.pbc.scf.hf.get_ovlp(cell, kpt=array([ 0., 0., 0.]))[source]

Get the overlap AO matrix.

pyscf.pbc.scf.hf.get_rho(mf, dm=None, grids=None, kpt=None)[source]

Compute density in real space

pyscf.pbc.scf.hf.get_t(cell, kpt=array([ 0., 0., 0.]))[source]

Get the kinetic energy AO matrix.

pyscf.pbc.scf.hf.init_guess_by_chkfile(cell, chkfile_name, project=None, kpt=None)[source]

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

Returns:
Density matrix, (nao,nao) ndarray
pyscf.pbc.scf.hf.makov_payne_correction(mf)[source]

Makov-Payne correction (Phys. Rev. B, 51, 4014)

Unrestricted Hartree-Fock for periodic systems at a single k-point

See Also:
pyscf/pbc/scf/khf.py : Hartree-Fock for periodic systems with k-point sampling
class pyscf.pbc.scf.uhf.UHF(cell, kpt=array([ 0., 0., 0.]), exxdiv='ewald')[source]

UHF class for PBCs.

convert_from_(mf)[source]

Convert given mean-field object to RHF/ROHF

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

Dipole moment in the unit cell.

Args:

cell : an instance of Cell

dm_kpts (a list of ndarrays) : density matrices of k-points

Return:
A list: the dipole moment on x, y and z components
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

get_bands(kpts_band, cell=None, dm=None, kpt=None)[source]

Get energy bands at the given (arbitrary) ‘band’ k-points.

Returns:
mo_energy
: (nmo,) ndarray or a list of (nmo,) ndarray
Bands energies E_n(k)
mo_coeff
: (nao, nmo) ndarray or a list of (nao,nmo) ndarray
Band orbitals psi_n(k)
get_j(cell=None, dm=None, hermi=1, kpt=None, kpts_band=None)

Compute J matrix for the given density matrix and k-point (kpt). When kpts_band is given, the J matrices on kpts_band are evaluated.

J_{pq} = sum_{rs} (pq|rs) dm[s,r]

where r,s are orbitals on kpt. p and q are orbitals on kpts_band if kpts_band is given otherwise p and q are orbitals on kpt.

get_jk(cell=None, dm=None, hermi=1, kpt=None, kpts_band=None)

Get Coulomb (J) and exchange (K) following scf.hf.RHF.get_jk_(). for particular k-point (kpt).

When kpts_band is given, the J, K matrices on kpts_band are evaluated.

J_{pq} = sum_{rs} (pq|rs) dm[s,r] K_{pq} = sum_{rs} (pr|sq) dm[r,s]

where r,s are orbitals on kpt. p and q are orbitals on kpts_band if kpts_band is given otherwise p and q are orbitals on kpt.

get_jk_incore(cell=None, dm=None, hermi=1, verbose=5, kpt=None)

Get Coulomb (J) and exchange (K) following scf.hf.RHF.get_jk_().

Incore version of Coulomb and exchange build only. Currently RHF always uses PBC AO integrals (unlike RKS), since exchange is currently computed by building PBC AO integrals.

get_k(cell=None, dm=None, hermi=1, kpt=None, kpts_band=None)

Compute K matrix for the given density matrix.

get_rho(mf, dm=None, grids=None, kpt=None)

Compute density in real space

pyscf.pbc.scf.uhf.dip_moment(cell, dm, unit='Debye', verbose=3, grids=None, rho=None, kpt=array([ 0., 0., 0.]))[source]

Dipole moment in the unit cell.

Args:

cell : an instance of Cell

dm_kpts (a list of ndarrays) : density matrices of k-points

Return:
A list: the dipole moment on x, y and z components
pyscf.pbc.scf.uhf.get_rho(mf, dm=None, grids=None, kpt=None)[source]

Compute density in real space

pyscf.pbc.scf.uhf.init_guess_by_chkfile(cell, chkfile_name, project=None, kpt=None)[source]

Read the HF results from checkpoint file and make the density matrix for UHF initial guess.

Returns:
Density matrix, (nao,nao) ndarray

Hartree-Fock for periodic systems with k-point sampling

See Also:
hf.py : Hartree-Fock for periodic systems at a single k-point
class pyscf.pbc.scf.khf.KSCF(cell, kpts=array([[ 0., 0., 0.]]), exxdiv='ewald')[source]

SCF base class with k-point sampling.

Compared to molecular SCF, some members such as mo_coeff, mo_occ now have an additional first dimension for the k-points, e.g. mo_coeff is (nkpts, nao, nao) ndarray

Attributes:
kpts
: (nks,3) ndarray
The sampling k-points in Cartesian coordinates, in units of 1/Bohr.
dip_moment(cell=None, dm=None, unit='Debye', verbose=3, **kwargs)[source]

Dipole moment in the unit cell.

Args:

cell : an instance of Cell

dm_kpts (a list of ndarrays) : density matrices of k-points

Return:
A list: the dipole moment on x, y and z components
energy_elec(mf, dm_kpts=None, h1e_kpts=None, vhf_kpts=None)

Following pyscf.scf.hf.energy_elec()

get_bands(kpts_band, cell=None, dm_kpts=None, kpts=None)[source]

Get energy bands at the given (arbitrary) ‘band’ k-points.

Returns:
mo_energy
: (nmo,) ndarray or a list of (nmo,) ndarray
Bands energies E_n(k)
mo_coeff
: (nao, nmo) ndarray or a list of (nao,nmo) ndarray
Band orbitals psi_n(k)
get_fermi(mf, mo_energy_kpts=None, mo_occ_kpts=None)

Fermi level

get_grad(mo_coeff_kpts, mo_occ_kpts, fock=None)[source]

returns 1D array of gradients, like non K-pt version note that occ and virt indices of different k pts now occur in sequential patches of the 1D array

get_hcore(mf, cell=None, kpts=None)

Get the core Hamiltonian AO matrices at sampled k-points.

Args:
kpts : (nkpts, 3) ndarray
Returns:
hcore : (nkpts, nao, nao) ndarray
get_occ(mf, mo_energy_kpts=None, mo_coeff_kpts=None)

Label the occupancies for each orbital for sampled k-points.

This is a k-point version of scf.hf.SCF.get_occ

get_ovlp(mf, cell=None, kpts=None)

Get the overlap AO matrices at sampled k-points.

Args:
kpts : (nkpts, 3) ndarray
Returns:
ovlp_kpts : (nkpts, nao, nao) ndarray
get_rho(mf, dm=None, grids=None, kpts=None)

Compute density in real space

get_veff(cell=None, dm_kpts=None, dm_last=0, vhf_last=0, hermi=1, kpts=None, kpts_band=None)[source]

Hartree-Fock potential matrix for the given density matrix. See scf.hf.get_veff() and scf.hf.RHF.get_veff()

pyscf.pbc.scf.khf.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.pbc.scf.khf.dip_moment(cell, dm_kpts, unit='Debye', verbose=3, grids=None, rho=None, kpts=array([[ 0., 0., 0.]]))[source]

Dipole moment in the unit cell.

Args:

cell : an instance of Cell

dm_kpts (a list of ndarrays) : density matrices of k-points

Return:
A list: the dipole moment on x, y and z components
pyscf.pbc.scf.khf.energy_elec(mf, dm_kpts=None, h1e_kpts=None, vhf_kpts=None)[source]

Following pyscf.scf.hf.energy_elec()

pyscf.pbc.scf.khf.get_fermi(mf, mo_energy_kpts=None, mo_occ_kpts=None)[source]

Fermi level

pyscf.pbc.scf.khf.get_grad(mo_coeff_kpts, mo_occ_kpts, fock)[source]

returns 1D array of gradients, like non K-pt version note that occ and virt indices of different k pts now occur in sequential patches of the 1D array

pyscf.pbc.scf.khf.get_hcore(mf, cell=None, kpts=None)[source]

Get the core Hamiltonian AO matrices at sampled k-points.

Args:
kpts : (nkpts, 3) ndarray
Returns:
hcore : (nkpts, nao, nao) ndarray
pyscf.pbc.scf.khf.get_j(mf, cell, dm_kpts, kpts, kpts_band=None)[source]

Get the Coulomb (J) AO matrix at sampled k-points.

Args:
dm_kpts
: (nkpts, nao, nao) ndarray or a list of (nkpts,nao,nao) ndarray
Density matrix at each k-point. If a list of k-point DMs, eg, UHF alpha and beta DM, the alpha and beta DMs are contracted separately. It needs to be Hermitian.
Kwargs:
kpts_band
: (k,3) ndarray
A list of arbitrary “band” k-points at which to evalute the matrix.
Returns:
vj : (nkpts, nao, nao) ndarray or list of vj if the input dm_kpts is a list of DMs
pyscf.pbc.scf.khf.get_jk(mf, cell, dm_kpts, kpts, kpts_band=None)[source]

Get the Coulomb (J) and exchange (K) AO matrices at sampled k-points.

Args:
dm_kpts
: (nkpts, nao, nao) ndarray
Density matrix at each k-point. It needs to be Hermitian.
Kwargs:
kpts_band
: (3,) ndarray
A list of arbitrary “band” k-point at which to evalute the matrix.
Returns:
vj : (nkpts, nao, nao) ndarray vk : (nkpts, nao, nao) ndarray or list of vj and vk if the input dm_kpts is a list of DMs
pyscf.pbc.scf.khf.get_occ(mf, mo_energy_kpts=None, mo_coeff_kpts=None)[source]

Label the occupancies for each orbital for sampled k-points.

This is a k-point version of scf.hf.SCF.get_occ

pyscf.pbc.scf.khf.get_ovlp(mf, cell=None, kpts=None)[source]

Get the overlap AO matrices at sampled k-points.

Args:
kpts : (nkpts, 3) ndarray
Returns:
ovlp_kpts : (nkpts, nao, nao) ndarray
pyscf.pbc.scf.khf.get_rho(mf, dm=None, grids=None, kpts=None)[source]

Compute density in real space

pyscf.pbc.scf.khf.init_guess_by_chkfile(cell, chkfile_name, project=None, kpts=None)[source]

Read the KHF results from checkpoint file, then project it to the basis defined by cell

Returns:
Density matrix, 3D ndarray
pyscf.pbc.scf.khf.make_rdm1(mo_coeff_kpts, mo_occ_kpts, **kwargs)[source]

One particle density matrices for all k-points.

Returns:
dm_kpts : (nkpts, nao, nao) ndarray
pyscf.pbc.scf.khf.mulliken_meta(cell, dm_ao_kpts, verbose=5, pre_orth_method='ANO', s=None)[source]

Mulliken population analysis, based on meta-Lowdin AOs.

Note this function only computes the Mulliken population for the gamma point density matrix.

Hartree-Fock for periodic systems with k-point sampling

See Also:
hf.py : Hartree-Fock for periodic systems at a single k-point
class pyscf.pbc.scf.kuhf.KUHF(cell, kpts=array([[ 0., 0., 0.]]), exxdiv='ewald')[source]

UHF class with k-point sampling.

canonicalize(mf, mo_coeff_kpts, mo_occ_kpts, fock=None)

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

convert_from_(mf)[source]

Convert given mean-field object to KUHF

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

Dipole moment in the unit cell.

Args:

cell : an instance of Cell

dm_kpts (two lists of ndarrays) : KUHF density matrices of k-points

Return:
A list: the dipole moment on x, y and z components
energy_elec(mf, dm_kpts=None, h1e_kpts=None, vhf_kpts=None)

Following pyscf.scf.hf.energy_elec()

get_bands(kpts_band, cell=None, dm_kpts=None, kpts=None)[source]

Get energy bands at the given (arbitrary) ‘band’ k-points.

Returns:
mo_energy
: (nmo,) ndarray or a list of (nmo,) ndarray
Bands energies E_n(k)
mo_coeff
: (nao, nmo) ndarray or a list of (nao,nmo) ndarray
Band orbitals psi_n(k)
get_fermi(mf, mo_energy_kpts=None, mo_occ_kpts=None)

A pair of Fermi level for spin-up and spin-down orbitals

get_hcore(mf, cell=None, kpts=None)

Get the core Hamiltonian AO matrices at sampled k-points.

Args:
kpts : (nkpts, 3) ndarray
Returns:
hcore : (nkpts, nao, nao) ndarray
get_occ(mf, mo_energy_kpts=None, mo_coeff_kpts=None)

Label the occupancies for each orbital for sampled k-points.

This is a k-point version of scf.hf.SCF.get_occ

get_ovlp(mf, cell=None, kpts=None)

Get the overlap AO matrices at sampled k-points.

Args:
kpts : (nkpts, 3) ndarray
Returns:
ovlp_kpts : (nkpts, nao, nao) ndarray
get_rho(mf, dm=None, grids=None, kpts=None)

Compute density in real space

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
pyscf.pbc.scf.kuhf.canonicalize(mf, mo_coeff_kpts, mo_occ_kpts, fock=None)[source]

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

pyscf.pbc.scf.kuhf.dip_moment(cell, dm_kpts, unit='Debye', verbose=3, grids=None, rho=None, kpts=array([[ 0., 0., 0.]]))[source]

Dipole moment in the unit cell.

Args:

cell : an instance of Cell

dm_kpts (two lists of ndarrays) : KUHF density matrices of k-points

Return:
A list: the dipole moment on x, y and z components
pyscf.pbc.scf.kuhf.energy_elec(mf, dm_kpts=None, h1e_kpts=None, vhf_kpts=None)[source]

Following pyscf.scf.hf.energy_elec()

pyscf.pbc.scf.kuhf.get_fermi(mf, mo_energy_kpts=None, mo_occ_kpts=None)[source]

A pair of Fermi level for spin-up and spin-down orbitals

pyscf.pbc.scf.kuhf.get_occ(mf, mo_energy_kpts=None, mo_coeff_kpts=None)[source]

Label the occupancies for each orbital for sampled k-points.

This is a k-point version of scf.hf.SCF.get_occ

pyscf.pbc.scf.kuhf.get_rho(mf, dm=None, grids=None, kpts=None)[source]

Compute density in real space

pyscf.pbc.scf.kuhf.init_guess_by_chkfile(cell, chkfile_name, project=None, kpts=None)[source]

Read the KHF results from checkpoint file, then project it to the basis defined by cell

Returns:
Density matrix, 3D ndarray
pyscf.pbc.scf.kuhf.make_rdm1(mo_coeff_kpts, mo_occ_kpts, **kwargs)[source]

Alpha and beta spin one particle density matrices for all k-points.

Returns:
dm_kpts : (2, nkpts, nao, nao) ndarray
pyscf.pbc.scf.kuhf.mulliken_meta(cell, dm_ao_kpts, verbose=5, pre_orth_method='ANO', s=None)[source]

Mulliken population analysis, based on meta-Lowdin AOs.

Note this function only computes the Mulliken population for the gamma point density matrix.