# 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 |