SCF and DFT methods#

Modules: pyscf.pbc.scf, pyscf.pbc.dft

Introduction#

PySCF supports periodic Hartree-Fock and density functional theory calculations with Brillouin zone sampling. The results of these calculations serve as input to a variety of periodic post-HF and post-DFT calculations. A minimal example of a periodic HF calculation on diamond with a 2x2x2 sampling of the Brillouin zone is shown below. Note that the kpts keyword argument must be in physical units (inverse bohr), which can be easily achieved using the Cell.make_kpts() method:

from pyscf.pbc import gto, scf
cell = gto.M(
    atom = '''C 0.0000 0.0000 0.0000
              C 0.8917 0.8917 0.8917''',
    a = '''0.0000 1.7834 1.7834
           1.7834 0.0000 1.7834
           1.7834 1.7834 0.0000''',
    pseudo = 'gth-pade',
    basis = 'gth-szv'
    )

kmf = scf.KRHF(cell, kpts=cell.make_kpts([2,2,2])).run()
# converged SCF energy = -10.9308552994574

See Crystal structure and Periodic boundary conditions for more details about the Cell object and Brillouin zone sampling.

Various integral schemes have been developed for periodic calculations in PySCF. A regular periodic calculation can be started with some type of density fitting. For large systems, multi-grid algorithm (see Multigrid) or range-separation (RS) integral algorithm (see Range-separation integration) can be considered to reduce the cost of Fock build. See How to choose integral scheme for the comparison of different integral schemes in periodic calculations.

Density fitting#

The default behavior is to use plane-wave density fitting (FFTDF). The number of plane-waves used as an auxiliary basis is controlled by kinetic energy cutoff, which is specified by the Cell.ke_cutoff parameter. The default value of this parameter is chosen to provide many digits of precision in the ERIs and the subsequent energies. If reduced precision is tolerable, this parameter can be set manually, resulting in significant speedups and memory savings. For the diamond example above, lowering the kinetic energy cutoff to 100 hartree changes the SCF energy by about 5 microhartree:

cell.ke_cutoff = 100.0 # in Hartree
cell.build()
kmf = scf.KRHF(...).run()
# converged SCF energy = -10.9308602914696

For many systems, Gaussian density fitting (GDF) is more economical, although it incurs larger errors than FFTDF. Periodic GDF can be activated in the same way as for molecules:

kmf = scf.KRHF().density_fit().run()

If a corresponding auxiliary basis is found for the chosen atomic orbital basis, it will be used. Otherwise, an even-tempered Gaussian basis will be used.

See Density fitting for crystalline calculations for more details about periodic density fitting.

Finite-size effects#

The long-ranged nature of the Coulomb interaction is responsible for a number of divergent contributions to the energy. For charge-neutral unit cells, the divergence of the nuclear repulsion energy, the electron-nuclear attraction energy, and the electronic Hartree energy cancel one another.

The nonlocal exact exchange energy exhibits an integrable divergence at \(G=0\) that is responsible for the leading-order finite-size error of HF and hybrid DFT calculations [60, 61, 62]. In PySCF, this exchange divergence can be addressed in a number of ways using the exxdiv keyword argument to the mean-field class, with the following possible values.

'ewald' (default)

The \(G=0\) value of the Coulomb potential is the supercell Madelung constant [60, 61, 62], which is evaluated by Ewald summation. The finite-size error of the exchange energy decays as \(N_k^{-1}\), where \(N_k\) is the number of k-points sampled in the Brillouin zone.
None

The \(G=0\) value of the Coulomb potential is set to zero. The finite-size error of the exchange energy decays slowly as \(N_k^{-1/3}\).
'vcut_sph'

The Coulomb potential is spherically truncated in real space at a radius equal to half of the supercell side length [63]. The finite-size error of the exchange energy decays as \(\exp(-aN_k^{1/3})\). Only supported with plane-wave density fitting (FFTDF).
'vcut_ws'

The Coulomb potential is truncated outside of the Wigner-Seitz supercell [62], which is more appropriate than spherical truncation for anisotropic cells. The finite-size error of the exchange energy decays as \(\exp(-aN_k^{1/3})\). Only supported with plane-wave density fitting (FFTDF).

An example calculation with exchange treated with the spherically truncated Coulomb potential is shown here:

kmf = scf.KRHF(cell, kpts=kpts, exxdiv='vcut_sph').run()

Band structure calculations#

After an SCF calculation has been performed, the band structure can be calculated non-self-consistently along a k-point path using the SCF.get_bands(kpts) function, where kpts is a list of k-points along which the band structure is desired.

Warning

For Hartree-Fock or hybrid DFT, the discontinuity of the exchange potential at \(G=0\) is problematic for band structure calculations. Using exxdiv='vcut_sph' with FFTDF is recommended instead. Alternatively, the SCF procedure can be repeated at each k-point, which is much more expensive but allows the use of any exxdiv or density fitting.

See examples/pbc/09-band_ase.py for an example DFT band structure calculation.

Add-ons#

All molecular SCF add-ons are also available for periodic SCF but must be accessed through the molecular pyscf.scf.addons module. Here we highlight a few of the most useful add-ons.

Linear dependencies#

The dense packing of atoms in solids combined with the use of diffuse atom-centered basis functions is responsible for frequent linear dependencies. The linear dependency can be eliminated by Cholesky orthogonalization:

from pyscf import scf as mol_scf
kmf = scf.KRHF(cell, kpts=kpts)
kmf = mol_scf.addons.remove_linear_dep_(kmf).run()

Smearing#

For metals or small band gap semiconductors, it can be useful to smear the orbital occupation numbers away from integer values. This can improve SCF convergence and can expedite convergence to the thermodynamic limit with k-point sampling. Because this approach assumes a finite electronic temperature, it yields an entropy and free energy:

kmf = scf.KRHF(cell, kpts=kpts)
kmf = scf.addons.smearing_(kmf, sigma=0.01, method='fermi').run()
print('Entropy = %s' % kmf.entropy)
print('Free energy = %s' % kmf.e_free)
print('Zero temperature energy = %s' % ((kmf.e_tot+kmf.e_free)/2))

Fermi-Dirac smearing (method='fermi') and Gaussian smearing (method='gauss') are supported.

Warning

Because most functions in PySCF assume integer occupations, they may fail if combined with a mean-field calculation that was performed with smearing.

Stability#

Periodic SCF solutions can be checked with stability analysis:

kmf = scf.KRHF(cell).run()
kmf.stability()

Multigrid#

For pure DFT calculations, multi-grid algorithm becomes efficient if the size of unit cell is relatively large (e.g. more than 10 atoms) or the kinetic energy cutoff is relatively high (e.g. more than 100k plane waves). This algorithm is implemented in the pyscf.pbc.dft.multigrid module:

from pyscf.pbc.dft import multigrid
mf = multigrid.multigrid(dft.KRKS(cell))
mf.run()

More examples can be found

#!/usr/bin/env python

'''
Use multi-grid to accelerate DFT numerical integration.
'''

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

cell = gto.M(
    verbose = 4,
    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 = 'sto3g',
    #basis = 'ccpvdz',
    #basis = 'gth-dzvp',
    #pseudo = 'gth-pade'
)

mf = dft.UKS(cell)
mf.xc = 'lda,vwn'

#
# There are two ways to enable multigrid numerical integration
#
# Method 1: use multigrid.multigrid function to update SCF object
#
mf = multigrid.multigrid(mf)
mf.kernel()

#
# Method 2: MultiGridFFTDF is a DF object.  It can be enabled by overwriting
# the default with_df object.
#
kpts = cell.make_kpts([4,4,4])
mf = dft.KRKS(cell, kpts)
mf.xc = 'lda,vwn'
mf.with_df = multigrid.MultiGridFFTDF(cell, kpts)
mf.kernel()

#
# MultiGridFFTDF can be used with second order SCF solver.
#
mf = mf.newton()
mf.kernel()

The multi-grid algorithm does not support exact exchange.

Range-separation integration#

This algorithm computes most of the four-center integrals in real space. For small systems, it is less efficient than most density fitting algorithms. When you need to handle a huge unit cell, or to compute exact exchange for many k points, or to use all-electron basis sets, you can consider to invoke this algorithm. For periodic SCF models, you can call the method jk_method:

kmf = scf.KRHF(cell, kpts).jk_method('RS')
kmf.kernel()

This algorithm has much smaller memory footprint than density fitting algorithm. It is also a good choice if your computer does not have enough memory to call density fitting methods.

Note

The implementation in current release does not support band structure calculations

How to choose integral scheme#

In the tables below, we provide a very rough estimation of the capability and characters for each integral algorithm in typical periodic mean-field calculations.

Number of basis functions for Gamma point calculations

Scheme	Pure DFT	hybrid DFT
FFTDF	< 5000	< 300
GDF	< 5000	< 5000
MDF	< 5000	< 2000
Multi-grid	100 - 100000	N/A
RS	100 - 10000	100 - 10000

Assuming 100 AOs per unit cell, the number of k-points

Scheme	Pure DFT	hybrid DFT
FFTDF	< 2000	< 10
GDF	< 2000	< 100
MDF	< 2000	< 50
Multi-grid	5000	N/A
RS	1000	1000

Relative performance for systems of 100 AOs per unit cell, 10 k-points

Scheme	Pure DFT	hybrid DFT
FFTDF	1x	1000x
GDF	10x	100x
MDF	30x	300x
Multi-grid	1x	N/A
RS	100x	100x

Other metrics

Scheme	All-electron	Warm-up cost	Accuracy	Memory cost	IO cost
FFTDF	Extremely slow	Low	High	Extremely high	None
GDF	Well supported	High	Low	Moderate	Huge
MDF	Well supported	High	Medium	Moderate	Huge
Multi-grid	Limited cases	Low	High	High	None
RS	Well supported	Low	High	Low	None