Density fitting (DF)

Modules: df, pbc.df

Introduction

Density fitting (DF), sometimes also called the resolution of identity (RI) approximation, is a method to approximate the four-index electron repulsion integrals (ERIs) by two- and three-index tensors. [1] In DF, the atomic orbital (AO) product space is expanded in terms of an auxiliary basis set. PySCF has extensive support for using DF in both molecular and crystalline calculations. This section covers the basic usage of DF in molecular applications. See Density fitting for crystalline calculations for the use of DF in periodic calculations.

Using DF in SCF calculations

DF is not used by default when initializing a SCF or MCSCF object, but can be invoked via the density_fit() method as follows:

from pyscf import gto, scf, mcscf
mol = gto.M(atom='N 0 0 0; N 0 0 1.2', basis='def2-tzvp')
mf = scf.RHF(mol).density_fit().run()    # output: -108.943773290737
mc = mcscf.CASSCF(mf, 8, 10).density_fit().run()    # output: -109.077600181849

Alternatively, one can initialize a DF object and overwrite the with_df attribute of the SCF object:

from pyscf import gto, scf, df
mol = gto.M(atom='N 0 0 0; N 0 0 1.2', basis='def2-tzvp')
mydf = df.DF(mol).build()
mf = scf.RHF(mol).density_fit()
mf.with_df = mydf
mf.kernel()    # output: -108.943773290737

Setting with_df to None turns DF off.

Using DF in post-SCF calculations

A post-SCF calculation will use DF automatically if the preceding SCF calculation is performed using DF. The example below performs a TDDFT/TDA calculation using DF:

from pyscf import gto, dft, tddft
mol = gto.M(atom='N 0 0 0; N 0 0 1.2', basis='def2-tzvp')
mf = dft.RKS(mol, xc="b3lyp").density_fit().run()
td = tddft.TDA(mf).run()
print(td.e)    # output: [0.27294247 0.28837975 0.2883798 ]

Alternatively, one can use the density_fit() method provided by most post-SCF classes in the same way as in a SCF calculation described above. This is useful when one wants to use DF to accelerate only the post-SCF calculations or to use a different auxiliary basis for the post-SCF calculations than that used by SCF. The example below performs a DF-MP2 calculation using the cc-pvtz-ri basis [2] on top of a SCF calculation that does not use DF:

from pyscf import gto, scf, mp
mol = gto.M(atom='N 0 0 0; N 0 0 1.2', basis='cc-pvtz')
mf = scf.RHF(mol).run()
mmp = mp.MP2(mf).density_fit(auxbasis="cc-pvtz-ri")
e_corr = mmp.kernel()[0]
print(e_corr)    # output: -0.4288734824009443

Choice of auxiliary basis

General considerations

The choice of auxiliary basis depends on the AO basis. By default, PySCF uses a pre-defined auxiliary basis set optimized for the used AO basis set, if one exists. This includes many commonly used AO basis sets in electronic structure calculations, e.g., the Ahlrichs def2 family, [3][4] and the Dunning cc family, [5].

When a pre-defined auxiliary basis set is not available, an even-tempered basis (ETB) set is generated by the following rule

\[\varphi_i = r^l \exp(-\alpha_l \beta^i r^2), \quad i = 0, 1, \cdots, n_l\]

where both \(\alpha_l\) and \(n_l\) are determined automatically from the AO basis. Specifically, \(\alpha_l\) is set to the smallest exponent of the AO products with angular momentum \(l\), and \(n_l\) is chosen so that the largest exponent of the AO products with angular momentum \(l\) lies in between \(\alpha_l \beta^{n_l-1}\) and \(\alpha_l \beta^{n_l}\). The parameter \(\beta\) controls the “density” of the fitting functions: smaller \(\beta\) gives an ETB of larger size and vice versa. The default is \(\beta = 2.0\).

The user can overwrite the default choice of the auxiliary basis by setting auxbasis upon initialization or at a later stage:

from pyscf import gto, scf, df
mol = gto.M(atom='N 0 0 0; N 0 0 1.2', basis='cc-pvdz')
mf = scf.RHF(mol).density_fit(auxbasis="weigend")
mf.kernel() # -108.910953335055
mf.with_df.auxbasis = "cc-pvdz-jkfit"   # this is the default
mf.kernel() # -108.913710743723
mf.with_df.auxbasis = df.aug_etb(mol, beta=1.7) # ETB with beta = 1.7
mf.kernel() # -108.914059329528

The user can print the used auxiliary basis set at any time by:

print(mf.with_df.auxmol.basis)

More examples on using customized auxiliary basis can be found in examples/df/01-auxbasis.py.

Special considerations for DFT

For DFT calculations with pure exchange-correlation functionals (i.e., LDA and GGA), the default auxiliary basis, which is designed for fitting both the Coulomb and the Hartree-Fock exchange integrals, may be unnecessarily large. We recommend using the def2-universal-jfit basis [6] for a more cost-effective choice as shown in the following example:

from pyscf import gto, dft
mol = gto.M(atom='N 0 0 0; N 0 0 1.2', basis='def2-tzvpp')
mf = dft.RKS(mol, xc="pbe").density_fit().run() # -109.432329411505
print(mf.with_df.auxmol.basis)                  # default: def2-tzvpp-jkfit
print(mf.with_df.auxmol.nao_nr())               # 154 aux basis functions
mf.with_df.auxbasis = "def2-universal-jfit"     # same as "weigend"
mf.kernel()                                     # -109.432334646585
print(mf.with_df.auxmol.nao_nr())               # 98 aux basis functions

Saving and reusing DF tensors

A key result of the DF class is the Cholesky decomposition of the ERIs (CDERIs) in terms of the 3-index tensor \(d_{L\mu\nu}\),

\[(\mu\nu|\lambda\sigma) = \sum_{L} d^*_{L\mu\nu} d_{L\lambda\sigma}\]

The build() method of DF computes the CDERIs. By default, the 3-index tensor defining the CDERIs is discarded once the calculation finishes. Sometimes it is useful to save it to disk for re-use in later calculations. This can be achieved by specifying a HDF5 file by setting _cderi_to_save either at the SCF stage:

mf = scf.RHF(mol).density_fit()
mf.with_df._cderi_to_save = 'saved_cderi.h5'
mf.kernel()

or by initializing a DF object separately:

mydf = df.DF(mol)
mydf.auxbasis = df.make_auxbasis(mol)
mydf._cderi_to_save = 'saved_cderi.h5'
mydf.build()

The saved DF tensor can be used later by setting _cderi to the HDF5 file:

mf = scf.RHF(mol).density_fit()
mf.with_df._cderi = 'saved_cderi.h5'
mf.kernel()

More examples on saving and using DF integrals and tensors can be found in examples/df/10-access_df_integrals.py, examples/df/11-access_df_tensor.py, and examples/df/40-precompute_df_integrals.py.

Advanced examples

More advanced examples of using the df module include

References

1

Brett I. Dunlap. Robust and variational fitting. Phys. Chem. Chem. Phys., 2:2113–2116, 2000. doi:10.1039/B000027M.

2

Florian Weigend, Andreas Köhn, and Christof Hättig. Efficient use of the correlation consistent basis sets in resolution of the identity mp2 calculations. J. Chem. Phys., 116:3175–3183, 2002. doi:10.1063/1.1445115.

3

Arnim Hellweg, Christof Hättig, Sebastian Höfener, and Wim Klopper. Optimized accurate auxiliary basis sets for ri-mp2 and ri-cc2 calculations for the atoms rb to rn. Theor. Chem. Acc., 117:587–597, 2007. doi:10.1007/s00214-007-0250-5.

4

Florian Weigend, Marco Häser, Holger Patzelt, and Reinhart Ahlrichs. Ri-mp2: optimized auxiliary basis sets and demonstration of efficiency. Chem. Phys. Lett., 294:143–152, 1998. doi:10.1016/S0009-2614(98)00862-8.

5

Florian Weigend. A fully direct ri-hf algorithm: implementation, optimised auxiliary basis sets, demonstration of accuracy and efficiency. Phys. Chem. Chem. Phys., 4:4285–4291, 2002. doi:10.1039/B204199P.

6

Florian Weigend. Accurate coulomb-fitting basis sets for h to rn. Phys. Chem. Chem. Phys., 8:1057–1065, 2006. doi:10.1039/B515623H.