comp-chem/cp2k
1
0
Fork 0
mirror of https://github.com/cp2k/cp2k.git synced 2026-07-28 06:05:29 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

Docs: Extend documentation of post-Hartree-Fock methods

Browse source
This commit is contained in:
Frederick Stein 2024-06-13 12:04:15 +02:00 • committed by Frederick Stein
parent 9d45d454af
commit 59ce2abea8
5 changed files with 248 additions and 13 deletions
Download patch file Download diff file Expand all files Collapse all files

3
.gitignore vendored
View file

@ -246,6 +246,9 @@ docs/_build/
docs/CP2K_INPUT/
docs/CP2K_INPUT.md
docs/bibliography.md
docs/cp2k_input.xml
docs/references.html
docs/units.html
# PyBuilder
.pybuilder/

1
docs/methods/post_hartree_fock/index.md
View file

@ -5,6 +5,7 @@
titlesonly:
maxdepth: 2
---
preliminaries
mp2
rpa
low-scaling

76
docs/methods/post_hartree_fock/mp2.md
View file

@ -6,20 +6,78 @@ $$
E^{(2)} = - \sum_{ijab} \frac{(ia|jb)[2(ia|jb)-(ib|ja)]}{\epsilon_{a}+\epsilon_{b}-\epsilon_{i}-\epsilon_{j}}
$$
with the GPW method as described in [](#DelBen2013) and [](#DelBen2012).
with the GPW method as described in [](#DelBen2012) and [](#DelBen2013). CP2K also implements double-hybrid functionals.
MP2 is controlled via the [WF_CORRELATION](#CP2K_INPUT.ATOM.METHOD.XC.WF_CORRELATION) section.
CP2K provides three different implementations: A canonical implementation, a GPW-based implementation and a RI-based implementation. Different functionals require rescaling of different contributions (singulet, triplet). This is achieved using the keywords [SCALE_S](#CP2K_INPUT.FORCE_EVAL.DFT.XC.WF_CORRELATION.SCALE_S) and [SCALE_T](#CP2K_INPUT.FORCE_EVAL.DFT.XC.WF_CORRELATION.SCALE_T) for all implementations.
## Canonical Implementation
The canonical implementation is available for GPW and GAPW reference calculations and thus allows for core-corrections. It is the most costly implementation. Gradients (Nuclear forces, stress tensors) are not analytically available. A simple WF_CORRELATION section looks like
```none
&WF_CORRELATION
METHOD RI_MP2_GPW
&WFC_GPW
CUTOFF 300
REL_CUTOFF 50
EPS_FILTER 1.0E-12
EPS_GRID 1.0E-8
&END
# See prelimaries
MEMORY 1200
NUMBER_PROC 1
# Only if not the Coulomb potential is required
&INTEGRALS
&END
&MP2
METHOD DIRECT_CANONICAL
&END
&END
```
## GPW-based Implementation
The MP2-based implementation is available only for GPW calculations.[](#DelBen2012) It is the cheaper than the canonical implementation. Gradients (Nuclear forces, stress tensors) are not analytically available. A simple WF_CORRELATION section looks like
```none
&WF_CORRELATION
# See prelimaries
MEMORY 1200
NUMBER_PROC 1
# Only if not the Coulomb potential is required
&INTEGRALS
&WFC_GPW
# To be optimized
CUTOFF 200
REL_CUTOFF 50
&END
&END
&MP2
METHOD MP2_GPW
&END
&END
```
## RI-based Implementation
The RI-based implementation is the most affordable MP2-implementation.[](#DelBen2013) Analytical gradients are available for GPW-based integrals (not for GAPW reference calculations).[](#DelBen2015b) MME-integration does not support stress tensor calculations, Obara-Saika integration does not implement any gradients. A simple WF_CORRELATION section looks like
```none
&WF_CORRELATION
# See prelimaries
MEMORY 1200
NUMBER_PROC 1
# Only if not the Coulomb potential is required
&INTEGRALS
&WFC_GPW
# To be optimized
CUTOFF 200
REL_CUTOFF 50
&END
&END
&RI_MP2
# Larger block sizes require more memory but reduce communication, -1 (default) let CP2K choose it
BLOCK_SIZE 2
# This keyword determines how many copies of the rank-three tensor are kept in the memory of all CPUs. Large number of groups increasethe memory demands but reduce communication (default: -1, automatic determination)
NUMBER_INTEGRATION_GROUPS 2
&END
&END
```
## Performance Considerations
MP2-calculations are generally very expensive for large systems due to their quintical scaling with respect to the number of atoms. Thus, for larger systems, one should switch to the cheaper RPA or SOS-MP2 methods.
CP2K employs local matrix-matrix multiplications for contractions. These can be accelerated with the SpLA library if CP2K was linked against it. In that way, GPU acceleration is possible if no GPU-accelerated DGEMM implementation is available and SpLA was configured accordingly.

101
docs/methods/post_hartree_fock/preliminaries.md Normal file
View file

@ -0,0 +1,101 @@
# Preliminaries
Here, we explain the relevant preparation steps and things to consider with all post-Hartree-Fock methods regarding choice of basis sets and pseudopotentials. Because of the different building blocks, we have to distinguish at least two kinds of basis sets: the primary basis set (PBS) to represent the orbital functions and the RI basis set to expand products of PBS functions into a set of auxiliary functions such that the four-center electron repulsion integrals (in Mulliken-notation) are approximated by
$$
(ia|jb)=\sum_{PQ}(ia|P) (P|Q)^{-1} (Q|jb)
$$
On top of the PBS and the RI basis set, the ADMM approximation of HF may introduce a third kind of basis set (see there for further information).
## Step 1: Choice of Primary Basis Sets
Because of the slow convergence of the calculated properties from post-Hartree-Fock calculations with respect to the size of the PBS, post-Hartree-Fock methods demand for larger PBSs than ordinary hybrid DFT or HF calculations and for very accurate results a basis set convergence scheme employing consistent PBS for accurate results. For this purpose, post-Hartree-Fock calculations are commonly performed with correlation-consistent (cc)-basis sets. Without an extrapolation scheme, we recommend basis sets of at least triple-zeta (cc-TZ) or augmented double-zeta (aug-cc-DZ) quality. With an extrapolation scheme, one needs a second basis set for each atom kind of different size. The usual extrapolation formula for correlation energies is given by
$$
E(X) = A + \frac{B}{X^3}
$$
with $X$ being 2, 3 or 4 in case of DZ, TZ or QZ zeta basis sets, respectively. If possible, one should refrain from non-augmented DZ basis sets as their results are usually not converged enough for extrapolation. Suitable basis sets can be found in the `BASIS_RI_cc-TZ` and `BASIS_ccGRB` basis set files in the data directory.
## Step 2: Choice of RI Basis Sets (Only for RI-based Methods)
The RI approximation asks for an additional auxiliary basis set. They are usually specified as [BASIS_SET RI_AUX](#CP2K_INPUT.FORCE_EVAL.SUBSYS.KIND.BASIS_SET) in the input file. CP2K provides several options to find suitable RI basis sets:
- For a few main group elements, the `BASIS_RI_cc-TZ` file contains RI basis sets for the given PBSs. The literature may help with missing elements.
- Optimize RI basis sets yourself. This is recommended if a larger number of calculations is required to reduce the computational costs significantly. Check the [OPT_RI_BASIS](#CP2K_INPUT.FORCE_EVAL.DFT.XC.WF_CORRELATION.RI.OPT_RI_BASIS) and the regtest suite (QS/regtest-ri-opt) fur further information.
- In all other cases, CP2K can generate RI basis sets automatically (see the [AUTO_BASIS RI_AUX](#CP2K_INPUT.FORCE_EVAL.DFT.AUTO_BASIS). These basis sets provide a decent accuracy but are roughly twice as large as optimized basis sets.
## Step 3: Pseudopotentials
Pseudopotentials should reflect the underlying hybrid or HF calculation used to determine orbitals and orbital energies (see the respective sections).
## Step 4: Preoptimize Reference Orbitals
It is recommended to restart the orbital calculation step from preoptimized HF or hybrid orbitals. This reduces computational costs because DFT and HF calculations do not perform well in case of the large number of CPU cores required for post-Hartree-Fock methods.
## Step 5: Setup of Integral Calculation Step
The integration step is setup with the [INTEGRALS](#CP2K_INPUT.FORCE_EVAL.DFT.XC.WF_CORRELATION.INTEGRALS) section. Because most densities occurring in the formalisms of the post-Hartree-Fock methods have a zero net-charge, such that the Coulomb operator does not need to be truncated it is required by HF. CP2K has its own integration routines for post-Hartree-Fock methods exploiting this circumstance. The available integration methods (compare the [ERI_METHOD](#CP2K_INPUT.FORCE_EVAL.DFT.XC.WF_CORRELATION.INTEGRALS.ERI_METHOD) keyword) are Gaussian-Plane-Wave integration (GPW), Minimax-Ewald integration (MME) and Obara Saika integration (OS).
GPW is fully supported, MME does not implement stress tensors and does not support short-ranged operators (truncated Coulomb, erfc-Coulomb, ...), OS implements gradients only in case of low scaling methods and does not support automatically generated RI basis sets. The GPW and MME integration methods are further configured in the respective [WFC_GPW](#CP2K_INPUT.FORCE_EVAL.DFT.XC.WF_CORRELATION.INTEGRALS.WFC_GPW) and [ERI_MME](#CP2K_INPUT.FORCE_EVAL.DFT.XC.WF_CORRELATION.INTEGRALS.ERI_MME) sections. In case of the GPW integration method, the respective CUTOFF parameters should be tuned as known from ordinary G(A)PW integrations in CP2K although the primary cutoff parameter can be chosen much smaller (150-300 Ry) then usual.
The RI metric can be set by the user with the [RI](#CP2K_INPUT.FORCE_EVAL.DFT.XC.WF_CORRELATION.RI) section. In case of low scaling methods (with activated [LOW_SCALING](#CP2K_INPUT.FORCE_EVAL.DFT.XC.WF_CORRELATION.LOW_SCALING) section), the metric is by default the overlap metric, else it will be the potential operator will be used. If the default is not requested, the RI section should be set.
## General Structure of Post-Hartree-Fock Methods
The general structure is very similar. Each post-Hartree-Fock method has its own section which has to be activated to start the respective calculation. Each section has subsections to further configure the respective kind of calculation. Most post-Hartree-Fock calculations require the following structure
```none
&WF_CORRELATION
# Determine the available memory
MEMORY 1000
# Determines the size of sub groups used in most methods
# The lower the number, the less communication but the more memory (only relevant for very large systems)
GROUP_SIZE 1
&INTEGRALS
# Here, we show the setup with the GPW integration
ERI_METHOD GPW
# Tune the accuracy and performance of the integral calculation step (default is usually sufficient)
SIZE_LATTICE_SUM 5
&WFC_GPW
# Set it as large as necessary and as small as possible (check energies)
CUTOFF 200
REL_CUTOFF 50
&END
# Only necessary for special functionals
&INTERACTION_POTENTIAL
&END
# Relevant only if low-scaling methods or a potential operator without short-range contributions is requested
&RI
# Set up the RI metric
&RI_METRIC
&END
&END
&END
# Set the section corresponding to the requested method (mixtures of different methods are not possible)
&<NAME_OF_METHOD>
<OPTIONS_TO_CONFIGURE_METHOD>
&END
&END
```
## Gradient Calculations
Gradient calculations are available for RI-MP2, RI-RPA and RI-SOS-MP2 calculations. In case of low-scaling implementations, consult the [CPHF](#CP2K_INPUT.FORCE_EVAL.DFT.XC.WF_CORRELATION.LOW_SCALING.CPHF) section in the LOW_SCALING section, otherwise the [CANONICAL_GRADIENTS](#CP2K_INPUT.FORCE_EVAL.DFT.XC.WF_CORRELATION.CANONICAL_GRADIENTS) section for further information. The setup of the CANONICAL_GRADIENTS section islooks like
```none
&CANONICAL_GRADIENTS
# This threshold switches to explicit contractions for almost degenerate orbital pairs (excluding diagonal elements)
EPS_CANONICAL 1.0E-6
# Try to turn it on, but it may crash in the execution step
FREE_HFX_BUFFER .FALSE.
&CPHF
# This threshold is to be optimized (lower values increase accuracy but require more time)
EPS_CONV 1.0E-6
MAX_ITER 20
&END
&END
```
## Basis Set Superposition Error (BSSE)
Energy calculations using post-Hartree-Fock methods suffer from severe basis set superposition error (BSSE). In CP2K, this is enhanced by the many different basis sets in use. Check the [BSSE](#CP2K_INPUT.FORCE_EVAL.BSSE) section for further information on the automatic setup of BSSE calculations and the [GHOST](#CP2K_INPUT.FORCE_EVAL.SUBSYS.KIND.GHOST) for the manual setup of BSSE calculations.

80
docs/methods/post_hartree_fock/rpa.md
View file

@ -1,7 +1,79 @@
# Random Phase Approximation
# Random-Phase Approximation and Laplace-Transformed Scaled-Opposite-Spin-MP2
Unfortunately, nobody has gotten around to writing this page yet :-(
The direct RI-RPA method computes
In the meantime, the following links might be helpful:
$$
E^{RI-dRPA} &= - \frac{1}{4\pi}\int_{-\infty}^{\infty}d\omega\text{Tr}\left(\text{log}(1+Q(\omega))-Q(\omega)\right)\\
Q_{RS}(\omega) &= 2\sum_{ia}B_{iaR}\frac{\epsilon_a-\epsilon_i}{\left(\epsilon_a-\epsilon_i\right)^2+\omega^2}B_{iaS}
$$
- <https://www.cp2k.org/exercises:2017_uzh_cp2k-tutorial:wfc>
whereas the Laplace-transformed Scaled-Opposite-Spin-MP2 (LT-RI-SOS-MP2) method computes
$$
E^{LT-RI-SOS-MP2} &= - \int_{0}^{\infty}d\tau\text{Tr}\left(\overline{Q}(\tau)^2\right)\\
\overline{Q}_{RS}(\tau) &= \sum_{ia}B_{iaR}e^{-\left(\epsilon_a-\epsilon_i\right)\tau}B_{iaS}
$$
as implemented in [](#DelBen2013). The integration is performed numerically. Double-hybrid functionals are available in CP2K.
CP2K provides two different implementations: a quartically-scaling and a low-scaling cubically-scaling implementation. Different functionals require rescaling which is achieved using the keywords [SCALE_RPA](#CP2K_INPUT.FORCE_EVAL.DFT.XC.WF_CORRELATION.RI_RPA.SCALE_RPA) in case of RPA and [SCALE_S](#CP2K_INPUT.FORCE_EVAL.DFT.XC.WF_CORRELATION.SCALE_S) in case of LT-RI-SOS-MP2.
Here, we will focus on the quartically-scaling implementation only.
## RI-dRPA
CP2K implements two quadrature schemes: Clenshaw-Curtis and Minimax. The first requires 30-40 quadrature points, the latter 6-8 quadrature points. Minimax quadrature rules have to be preoptimized such that not all possible numbers of quadrature points are available.
RPA correlation energies are usually combined with exact exchange energies. In CP2K, this is available by activating the [HF](#CP2K_INPUT.FORCE_EVAL.DFT.XC.WF_CORRELATION.RI_RPA.HF) section which is setup like an ordinary HF section.
Several Beyond-RPA schemes are available: The Renormalized Screened Exchange (RSE) correction, the Approximate Exchange Kernel (AXK) and the Second-Order Screened Exchange corrections (SOSEX). These are enabled with the [RSE](#CP2K_INPUT.FORCE_EVAL.DFT.XC.WF_CORRELATION.RI_RPA.RSE) keyword and the [EXCHANGE_CORRECTION](#CP2K_INPUT.FORCE_EVAL.DFT.XC.WF_CORRELATION.RI_RPA.EXCHANGE_CORRECTION) section.
This results into the following possible RPA section
```none
&RI_RPA
# Choose it as large as necessary and as small as possible
QUADRATURE_POINTS 6
# Choose it as large as possible (must be a divisor of the number of quadrature points and the number of processes)
# -1 is default and let CP2K decide on that value
# Larger values increase the memory demands but reduce communication
NUM_INTEG_GROUPS -1
# The RSE correction is only relevant for a non-HF reference
RSE .TRUE.
# Exchange corrections may be quite costly
&EXCHANGE_CORRECTION [NONE|AXK|SOSEX]
# The Hartree-Fock-based implementation scales better for larger systems but introduces more noise
USE_HFX_IMPLEMENTATION F
# This parameter is ignored if USE_HFX_IMPLEMENTATION is set to T
# Larger values improve performance but increase the memory demands
BLOCK_SIZE 16
&END
&END
```
## RI-dRPA Gradient Calculations
Analytical gradients are only available for RPA calculations using a minimax grid, but not for the beyond-RPA methods (RSE, AXK, SOSEX). The general setup is similar to RI-MP2 gradients. In addition, there are two more relevant keywords: [DOT_PRODUCT_BLKSIZE](#CP2K_INPUT.FORCE_EVAL.DFT.XC.WF_CORRELATION.CANONICAL_GRADIENTS.DOT_PRODUCT_BLKSIZE) and [MAX_PARALLEL_COMM](#CP2K_INPUT.FORCE_EVAL.DFT.XC.WF_CORRELATION.CANONICAL_GRADIENTS.MAX_PARALLEL_COMM). The first splits the contraction along the auxiliary index to improve numerical stability. By default, this feature is turned off. The second keyword determines the number of parallel communication channels for non-blocking communication. Larger numbers allow more overlap but increase the memory requirements. Larger values than 3 are commonly not necessary.
## LT-RI-SOS-MP2
CP2K implements two quadrature schemes only the Minimax scheme requiring usually 6-8 quadrature points. Minimax quadrature rules have to be preoptimized such that not all possible numbers of quadrature points are available. Scaling of the energy contribution is possible with the [SCALE_S](#CP2K_INPUT.FORCE_EVAL.DFT.XC.WF_CORRELATION.SCALE_S) keyword which is also used by MP2 calculations.
The input is simpler than in the RPA-case
```none
&RI_SOS_MP2
# Works similar than in RPA
NUM_INTEG_GROUPS -1
# Larger values improve the accuracy but increase the costs
QUADRATURE_POINTS 6
&END
```
## LT-RI-SOS-MP2 Gradient Calculations
Analytical gradients are available and work similarly to RI-MP2 calculations and RI-RPA calculations. As in case of RI-MP2, it makes use of the [EPS_CANONICAL](#CP2K_INPUT.FORCE_EVAL.DFT.XC.WF_CORRELATION.CANONICAL_GRADIENTS.EPS_CANONICAL) keyword. Larger values improve the numerical accuracy but increase the computational costs significantly because CP2K assumes that the number of relevant pairs is negligible.
## Performance Considerations
Without employing the low-scaling implementation, RI-RPA and RI-SOS-MP2 calculations scale quartically with respect to the number of atoms. The implementation relies on parallel matrix-matrix multiplications. Their costs can be reduced with the COSMA library and accelerated on GPUs if COSMA was configured accordingly.
HF calculations may be accelerated using ADMM which is recommended if very diffuse basis functions are employed. The RSE correction employs the same HF section as for the calculation of the HF exchange energy. The costs of AXK and SOSEX depend on the implementation in use. The HF-based implementation is cheaper for larger systems. Its performance depends on the available memory per process group (see [GROUP_SIZE](#CP2K_INPUT.FORCE_EVAL.DFT.XC.WF_CORRELATION.GROUP_SIZE) and [MAX_MEMORY](#CP2K_INPUT.FORCE_EVAL.DFT.XC.WF_CORRELATION.RI_RPA.HF.MEMORY.MAX_MEMORY)). The costs becomes neglegible for large systems. Without the HF implementation, CP2K relies on SpLA-accelerated local matrix-matrix multiplications. These become the bottleneck for larger systems but introduce less numerical noise than the HF-based implementation.