# 21.4. pbc.df — PBC denisty fitting¶

## 21.4.1. Introduction¶

The `pbc.df`

module provides the fundamental functions to handle the
density fitting (DF) integral tensors required by the gamma-point and k-point
PBC calculations. There are four types of DF methods available for PBC
systems. They are FFTDF (plane-wave density fitting with fast Fourier
transformation), AFTDF (plane-wave density fitting with analytical Fourier
transformation), GDF (Gaussian density fitting) and MDF (mixed density fitting).
The Coulomb integrals and nuclear attraction integrals in the PBC calculations
are all computed with DF technique. The default scheme is FFTDF.

The characters of these PBC DF methods are summarized in the following table

Subject | FFTDF | AFTDF | GDF | MDF |

Initialization | No | No | Slow | Slow |

HF Coulomb matrix (J) | Fast | Slow | Fast | Moderate |

HF exchange matrix (K) | Slow | Slow | Fast | Moderate |

Building ERIs | Slow | Slow | Fast | Moderate |

All-electron calculation | Huge error | Large error | Accurate | Most accurate |

Low-dimension system | N/A | 0D,1D,2D | 0D,1D,2D | 0D,1D,2D |

### 21.4.1.1. FFTDF — FFT-based density fitting¶

FFTDF represents the method to compute electron repulsion integrals in reciprocal space with the Fourier transformed Coulomb kernel

\(\mathbf{G}\) is the plane wave vector. \(\rho_{ij}(\mathbf{G})\) is the Fourier transformed orbital pair

Here are some examples to initialize FFTDF object:

```
>>> import numpy as np
>>> from pyscf.pbc import gto, df, scf
>>> cell = gto.M(atom='He 1 1 1', a=np.eye(3)*2, basis='3-21g')
>>> fftdf = df.FFTDF(cell)
>>> print(fftdf)
<pyscf.pbc.df.fft.FFTDF object at 0x7f599dbd6450>
>>> mf = scf.RHF(cell)
>>> print(mf.with_df)
<pyscf.pbc.df.fft.FFTDF object at 0x7f59a1a10c50>
```

As the default integral scheme of PBC calculations, FFTDF is created when
initializing the PBC mean-field object and held in the attribute `with_df`

.

#### 21.4.1.1.1. Nuclear type integrals¶

PBC nuclear-electron interaction and pseudo-potential (PP) integrals can be
computed with the FFTDF methods `FFTDF.get_nuc()`

and `FFTDF.get_pp()`

.
`FFTDF.get_nuc()`

function only evaluates the integral of the point charge.
If PP was specified in the cell object, `FFTDF.get_nuc()`

produces the
integrals of the point nuclei with the effective charges. If PP was not
defined in the cell object, `FFTDF.get_pp()`

and `FFTDF.get_nuc()`

produce the same integrals. Depending on the input k-point(s),
the two functions can produce the nuclear-type integrals for a single k-point or
a list of nuclear-type integrals for the k-points. By default, they compute the
nuclear-type integrals of Gamma point:

```
>>> vnuc = fftdf.get_pp()
>>> print(vnuc.shape)
(2, 2)
>>> kpts = cell.make_kpts([2,2,2])
>>> vnuc = fftdf.get_pp(kpts)
>>> print(vnuc.shape)
(8, 2, 2)
>>> vnuc = fftdf.get_pp(kpts)
>>> print(vnuc.shape)
(2, 2)
```

#### 21.4.1.1.2. Hartree-Fock Coulomb and exchange¶

`FFTDF`

class provides a method `FFTDF.get_jk()`

to compute
Hartree-Fock Coulomb matrix (J) and exchange matrix (K). This method can take
one density matrix or a list of density matrices as input and return the J and K
matrices for each density matrix:

```
>>> dm = numpy.random.random((2,2))
>>> j, k = fftdf.get_jk(dm)
>>> print(j.shape)
(2, 2)
>>> dm = numpy.random.random((3,2,2))
>>> j, k = fftdf.get_jk(dm)
>>> print(j.shape)
(3, 2, 2)
```

When k-points are specified, the input density matrices should have the correct shape that matches the number of k-points:

```
>>> kpts = cell.make_kpts([1,1,3])
>>> dm = numpy.random.random((3,2,2))
>>> j, k = fftdf.get_jk(dm, kpts=kpts)
>>> print(j.shape)
(3, 2, 2)
>>> dm = numpy.random.random((5,3,2,2))
>>> j, k = fftdf.get_jk(dm, kpts=kpts)
>>> print(j.shape)
(5, 3, 2, 2)
```

#### 21.4.1.1.3. 4-index ERI tensor and integral transformation¶

4-index electron repulsion integrals can be computed with `FFTDF.get_eri()`

and `FFTDF.ao2mo()`

methods. Given 4 k-points(s) (corresponding to the 4
AO indices), `FFTDF.get_eri()`

method produce the regular 4-index ERIs
\((ij|kl)\) in AO basis. The 4 k-points should follow the law of momentum
conservation

By default, four \(\Gamma\)-points are assigned to the four AO indices. As the format of molecular ERI tensor, the PBC ERI tensor is reshaped to a 2D array:

```
>>> eri = fftdf.get_eri()
>>> print(eri.shape)
(4, 4)
>>> eri = fftdf.get_eri([kpts[0],kpts[0],kpts[1],kpts[1]])
>>> print(eri.shape)
(4, 4)
```

`FFTDF.ao2mo()`

function applies integral transformation for the given four
sets of orbital coefficients, four input k-points. The four k-points need to
follow the momentum conservation law. Similar to `FFTDF.get_eri()`

, the
returned integral tensor is shaped to a 2D array:

```
>>> orbs = numpy.random.random((4,2,2))
>>> eri_mo = fftdf.get_eri(orbs, [kpts[0],kpts[0],kpts[1],kpts[1]])
>>> print(eri_mo.shape)
(4, 4)
```

#### 21.4.1.1.4. Kinetic energy cutoff¶

The accuracy of FFTDF integrals are affected by the kinetic energy cutoff. The
default kinetic energy cutoff is a conservative estimation based on the basis
set and the lattice parameter. You can adjust the attribute `FFTDF.gs`

(the numbers of grid points in each positive direction) to change the kinetic
energy cutoff. If any values in `FFTDF.gs`

is too small to reach the
required accuracy `cell.precision`

, `FFTDF`

may output a warning
message, eg:

```
WARN: ke_cutoff/gs (12.437 / [3, 4, 4]) is not enough for FFTDF to get integral accuracy 1e-08.
Coulomb integral error is ~ 2.6 Eh.
Recomended ke_cutoff/gs are 538.542 / [20 20 20].
```

In this warning message, `Coulomb integral error`

is a rough estimation for
the largest error of the matrix elements of the two-electron Coulomb integrals.
The overall computational error may be varied by 1 - 2 orders of magnitude.

### 21.4.1.2. AFTDF — AFT-based density fitting¶

AFTDF mans that the Fourier transform of the orbital pair is computed analytically

To enable AFTDF in the calculation, `AFTDF`

object can be initialized
and assigned to `with_df`

object of mean-field object:

```
>>> import numpy as np
>>> from pyscf.pbc import gto, df, scf
>>> cell = gto.M(atom='He 1 1 1', a=np.eye(3)*2, basis='3-21g')
>>> aft = df.AFTDF(cell)
>>> print(aft)
<pyscf.pbc.df.aft.AFTDF object at 0x7ff8b1893d90>
>>> mf = scf.RHF(cell)
>>> mf.with_df = aft
```

Generally, AFTDF is slower than FFTDF method.

`AFTDF`

class offers the same methods as the `FFTDF`

class.
Nuclear and PP integrals, Hartree-Fock J and K matrices, electron repulsion
integrals and integral transformation can be computed with functions
`AFTDF.get_nuc()`

, `AFTDF.get_pp()`

, `AFTDF.get_jk()`

,
`AFTDF.get_eri()`

and `AFTDF.ao2mo()`

using the same calling APIs as the
analogy functions in FFTDF — FFT-based density fitting.

#### 21.4.1.2.1. Kinetic energy cutoff¶

`AFTDF`

also makes estimation on the kinetic energy cutoff. When the
any values of `AFTDF.gs`

are too small for required accuracy
`cell.precision`

, this class also outputs the
`Coulomb integral error`

warning message as the `FFTDF`

class.

### 21.4.1.3. GDF — Gaussian density fitting¶

GDF is an analogy of the conventional density fitting method with periodic
boundary condition. The auxiliary fitting basis in PBC GDF is periodic Gaussian
function (To ensure the long range Coulomb integrals converging in the real
space lattice summation, the multipoles are removed from the auxiliary basis).
`GDF`

object can be initialized and enabled in the SCF calculation in two
ways:

```
>>> import numpy as np
>>> from pyscf.pbc import gto, df, scf
>>> cell = gto.M(atom='He 1 1 1', a=np.eye(3)*2, basis='3-21g')
>>> gdf = df.GDF(cell)
>>> mf = scf.RHF(cell)
>>> mf.with_df = gdf
>>> mf.run()
>>> # Using SCF.density_fit method
>>> mf = scf.RHF(cell).density_fit().run()
>>> print(mf.with_df)
<pyscf.pbc.df.df.GDF object at 0x7fec7722aa10>
```

Similar to the molecular code, `SCF.density_fit()`

method returns a
mean-field object with `GDF`

as the integral engine.

In the `GDF`

method, the DF-integral tensor is precomputed and stored
on disk. `GDF`

method supports both the \(\Gamma\)-point ERIs and
the ERIs of different k-points. `GDF.kpts`

should be specified before
initializing `GDF`

object. `GDF`

class provides the same APIs as
the `FFTDF`

class to compute nuclear integrals and electron Coulomb
repulsion integrals:

```
>>> import numpy as np
>>> from pyscf.pbc import gto, df, scf
>>> cell = gto.M(atom='He 1 1 1', a=np.eye(3)*2, basis='3-21g')
>>> gdf = df.GDF(cell)
>>> gdf.kpts = cell.make_kpts([2,2,2])
>>> gdf.get_eri([kpts[0],kpts[0],kpts[1],kpts[1]])
```

In the mean-field calculation, assigning `kpts`

attribute to mean-field
object updates the `kpts`

attribute of the underlying DF method:

```
>>> import numpy as np
>>> from pyscf.pbc import gto, df, scf
>>> cell = gto.M(atom='He 1 1 1', a=np.eye(3)*2, basis='3-21g')
>>> mf = scf.KRHF(cell).density_fit()
>>> kpts = cell.make_kpts([2,2,2])
>>> mf.kpts = kpts
>>> mf.with_df.get_eri([kpts[0],kpts[0],kpts[1],kpts[1]])
```

Once the GDF integral tensor was initialized, the `GDF`

can be only used
with certain k-points calculations. An incorrect `kpts`

argument can lead
to a runtime error:

```
>>> import numpy as np
>>> from pyscf.pbc import gto, df, scf
>>> cell = gto.M(atom='He 1 1 1', a=np.eye(3)*2, basis='3-21g')
>>> gdf = df.GDF(cell, kpts=cell.make_kpts([2,2,2]))
>>> kpt = np.random.random(3)
>>> gdf.get_eri([kpt,kpt,kpt,kpt])
RuntimeError: j3c for kpts [[ 0.53135523 0.06389596 0.19441766]
[ 0.53135523 0.06389596 0.19441766]] is not initialized.
You need to update the attribute .kpts then call .build() to initialize j3c.
```

The GDF initialization is very expensive. To reduce the initialization cost in a series of calculations, it would be useful to cache the GDF integral tensor in a file then load them into the calculation when needed. The GDF integral tensor can be saved and loaded the same way as we did for the molecular DF method (see Saving/Loading DF integral tensor):

```
import numpy as np
from pyscf.pbc import gto, df, scf
cell = gto.M(atom='He 1 1 1', a=np.eye(3)*2, basis='3-21g')
gdf = df.GDF(cell, kpts=cell.make_kpts([2,2,2]))
gdf._cderi_to_save = 'df_ints.h5' # To save the GDF integrals
gdf.build()
mf = scf.KRHF(cell, kpts=cell.make_kpts([2,2,2])).density_fit()
mf.with_df._cderi = 'df_ints.h5' # To load the GDF integrals
mf.run()
```

#### 21.4.1.3.1. Auxiliary Gaussian basis¶

GDF method requires a set of Gaussian functions as the density fitting auxiliary basis. See also DF auxiliary basis and Even-tempered auxiliary Gaussian basis for the choices of DF auxiliary basis in PySCF GDF code. There are not many optimized auxiliary basis sets available for PBC AO basis. You can use the even-tempered Gaussian functions as the auxiliary basis in the PBC GDF method:

```
import numpy as np
from pyscf.pbc import gto, df, scf
cell = gto.M(atom='He 1 1 1', a=np.eye(3)*2, basis='3-21g')
gdf = df.GDF(cell, kpts=cell.make_kpts([2,2,2]))
gdf.auxbasis = df.aug_etb(cell, beta=2.0)
gdf.build()
```

#### 21.4.1.3.2. Kinetic energy cutoff¶

GDF method does not require the specification of kinetic energy cutoff.
`cell.ke_cutoff`

and `cell.gs`

are ignored in the `GDF`

class. Internally, a small set of planewaves is used in the GDF method to
accelerate the convergence of GDF integrals in the real space lattice summation.
The estimated energy cutoff is generated in the `GDF`

class and stored in
the attribute `GDF.gs`

. It is not recommended to change this parameter.

### 21.4.1.4. MDF — mixed density fitting¶

MDF method combines the AFTDF and GDF in the same framework. The MDF auxiliary
basis is Gaussian and plane-wave mixed basis. `MDF`

object can be
created in two ways:

```
>>> import numpy as np
>>> from pyscf.pbc import gto, df, scf
>>> cell = gto.M(atom='He 1 1 1', a=np.eye(3)*2, basis='3-21g', ke_cutoff=10)
>>> mdf = df.MDF(cell)
>>> print(mdf)
<pyscf.pbc.df.mdf.MDF object at 0x7f4025120a10>
>>> mf = scf.RHF(cell).mix_density_fit().run()
>>> print(mf.with_df)
<pyscf.pbc.df.mdf.MDF object at 0x7f7963390a10>
```

The kinetic energy cutoff is specified in this example to constrain the number of
planewaves. The number of planewaves can also be controlled by through
attribute `MDF.gs`

.

In principle, the accuracy of MDF method can be increased by adding
more plane waves in the auxiliary basis. In practice, the linear dependency
between plane waves and Gaussians may lead to numerical stability issue.
The optimal accuracy (with reasonable computational cost) requires a reasonable
size of plan wave basis with a reasonable linear dependency threshold. A
threshold too large would remove many auxiliary functions while a threshold too
small would cause numerical instability.
.. In our preliminary test, `ke_cutoff=10`

is able to produce 0.1 mEh accuracy in
.. total energy.
The default linear dependency threshold is 1e-10. The threshold can be adjusted
through the attribute `MDF.linear_dep_threshold`

.

Like the GDF method, it is also very demanding to initialize the 3-center
Gaussian integrals in the MDF method. The 3-center Gaussian integral tensor can
be cached in a file and loaded to `MDF`

object at the runtime:

```
import numpy as np
from pyscf.pbc import gto, df, scf
cell = gto.M(atom='He 1 1 1', a=np.eye(3)*2, basis='3-21g')
mdf = df.MDF(cell, kpts=cell.make_kpts([2,2,2]))
mdf._cderi_to_save = 'df_ints.h5' # To save the GDF integrals
mdf.build()
mf = scf.KRHF(cell, kpts=cell.make_kpts([2,2,2])).mix_density_fit()
mf.with_df._cderi = 'df_ints.h5' # To load the GDF integrals
mf.run()
```

### 21.4.1.5. All-electron calculation¶

All-electron calculations with FFTDF or AFTDF methods requires high energy cutoff for most elements. It is recommended to use GDF or MDF methods in the all-electron calculations. In fact, GDF and MDF can also be used in PP calculations to reduce the number of planewave basis if steep functions are existed in the AO basis.

### 21.4.1.6. Low-dimension system¶

`AFTDF`

supports the systems with 0D (molecule), 1D and 2D periodic
boundary conditions. When computing the integrals of low-dimension systems, an
infinite vacuum is placed on the free boundary. You can set the
`cell.dimension`

, to enable the integral algorithms for
low-dimension systems in `AFTDF`

class:

```
import numpy as np
from pyscf.pbc import gto, df, scf
cell = gto.M(atom='He 1 1 1', a=np.eye(3)*2, basis='3-21g', dimension=1)
aft = df.AFTDF(cell)
aft.get_eri()
```

`GDF`

and `MDF`

all support the integrals of low-dimension system.
Similar to the usage of AFTDF method, you need to set `cell.dimension`

for
the low-dimension systems:

```
import numpy as np
from pyscf.pbc import gto, df, scf
cell = gto.M(atom='He 1 1 1', a=np.eye(3)*2, basis='3-21g', dimension=1)
gdf = df.GDF(cell)
gdf.get_eri()
```

See more examples in `examples/pbc/31-low_dimensional_pbc.py`

## 21.4.2. Interface to molecular DF-post-HF methods¶

PBC DF object is compatible to the molecular DF object. The \(\Gamma\)-point PBC SCF object can be directly passed to molecular DF post-HF methods for an electron correlation calculations in PBC:

```
import numpy as np
from pyscf.pbc import gto, df, scf
from pyscf import cc as mol_cc
cell = gto.M(atom='He 1 1 1', a=np.eye(3)*2, basis='3-21g', dimension=1)
mf = scf.RHF(cell).density_fit()
mol_cc.RCCSD(mf).run()
```

## 21.4.3. Examples¶

DF relevant examples can be found in the PySCF examples directory:

```
examples/pbc/10-gamma_point_scf.py
examples/pbc/11-gamma_point_all_electron_scf.py
examples/pbc/12-gamma_point_post_hf.py
examples/pbc/20-k_points_scf.py
examples/pbc/21-k_points_all_electron_scf.py
examples/pbc/30-ao_integrals.py
examples/pbc/30-ao_value_on_grid.py
examples/pbc/30-mo_integrals.py
examples/pbc/31-low_dimensional_pbc.py
```

## 21.4.4. Program reference¶

### 21.4.4.1. FFTDF class¶

### 21.4.4.2. FFTDF helper functions¶

JK with discrete Fourier transformation

Integral transformation with FFT

- (ij|kl) = int dr1 dr2 i*(r1) j(r1) v(r12) k*(r2) l(r2)
- = (ij|G) v(G) (G|kl)
- i*(r) j(r) = 1/N sum_G e^{iGr} (G|ij)
- = 1/N sum_G e^{-iGr} (ij|G)
- “forward” FFT:
- (G|ij) = sum_r e^{-iGr} i*(r) j(r) = fft[ i*(r) j(r) ]
- “inverse” FFT:
- (ij|G) = sum_r e^{iGr} i*(r) j(r) = N * ifft[ i*(r) j(r) ]
- = conj[ sum_r e^{-iGr} j*(r) i(r) ]

### 21.4.4.3. AFTDF class¶

### 21.4.4.4. AFTDF helper functions¶

JK with analytic Fourier transformation

Integral transformation with analytic Fourier transformation

### 21.4.4.5. GDF class¶

### 21.4.4.6. GDF helper functions¶

Density fitting with Gaussian basis Ref: J. Chem. Phys. 147, 164119 (2017)

### 21.4.4.7. MDF class¶

### 21.4.4.8. MDF helper functions¶

Exact density fitting with Gaussian and planewaves Ref: J. Chem. Phys. 147, 164119 (2017)