SCF and DFT methods

Modules: scf, dft, pbc.scf, 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 user_pbc_multigrid) or range-separation (RS) integral algorithm (see user_pbc_rsjk) can be considered to reduce the cost of Fock build. See user_pbc_df_comparison 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 [1][2][3]. 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 [1][2][3], 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 [4]. 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 [3], 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

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.

  1. 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

  1. 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

  1. 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

  1. 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

References

1(1,2)

Joachim Paier, Robin Hirschl, Martijn Marsman, and Georg Kresse. The perdew–burke–ernzerhof exchange-correlation functional applied to the g2-1 test set using a plane-wave basis set. J. Chem. Phys., 122(23):234102, 2005. doi:10.1063/1.1926272.

2(1,2)

Peter Broqvist, Audrius Alkauskas, and Alfredo Pasquarello. Hybrid-functional calculations with plane-wave basis sets: effect of singularity correction on total energies, energy eigenvalues, and defect energy levels. Phys. Rev. B, 80:085114, Aug 2009. doi:10.1103/PhysRevB.80.085114.

3(1,2,3)

Ravishankar Sundararaman and T. A. Arias. Regularization of the coulomb singularity in exact exchange by wigner-seitz truncated interactions: towards chemical accuracy in nontrivial systems. Phys. Rev. B, 87:165122, Apr 2013. doi:10.1103/PhysRevB.87.165122.

4

James Spencer and Ali Alavi. Efficient calculation of the exact exchange energy in periodic systems using a truncated coulomb potential. Phys. Rev. B, 77:193110, May 2008. doi:10.1103/PhysRevB.77.193110.