# Low-scaling post Hartree-Fock

## Brief theory recap

The low-scaling post-HF methods implemented in CP2K are SOS-MP2 and (direct) RPA. To achieve the desired low-scaling (sub-cubic), all equations are formulated in the atomic orbital (AO) basis, using the resolution-of-the-identity technique with a short range metric, and a numerical Laplace transform. A central quantity of post-HF is the collection of 4-center 2-electron electron repulsion integrals (ERIs). In the AO basis and within the RI approximation:

\[ (\mu\sigma|\nu\lambda) = (\mu\sigma\lfloor P)\ (P\lfloor Q)^{-1}\ (Q|R)\ (R\lfloor S)^{-1}\ (S\lfloor \nu\lambda) \]

where \(\lfloor\) denotes the short range RI metric, and \(P,Q,R,S\) are RI basis elements. A short range metric increases the sparsity of the 3-center RI integrals \((\mu\sigma\lfloor P)\), resulting in more efficient sparse tensor contractions. The shortest possible ranged RI metric is the overlap. Using the truncated Coulomb with a short truncation radius (1.5-2.0 Angstroms) is recommended in sparser systems (e.g. water), while the overlap is typically sufficient for dense solids. The 3- and 2-center RI ERIs are calculated analytically with the libint library. The 2-center potential ERIs, \((Q|R)\), are calculated numerically with the GPW scheme, because it uses the (very) long range 1/r Coulomb interaction.

During post-HF calculations, the ERIs discussed above are contracted with various other quantities, with a specific order of operations for maximum performance. The details can be found in Wilhelm2016b for the initial RPA energy implementation. Bussy2023 covers the detailed implementation of SOS-MP2 and RPA nuclear gradients. Calculations for the energy, forces, MD, cell and gemoetry optimization are possible. While these low-scaling methods have more favorable scaling than their standard implementations (sub-cubic vs quartic), they suffer from a heavy prefactor. Therefore, they are only interesting for large to very large systems (hundreds of atoms), and will still consume a lot of resources. For smaller system, the standard MP2 and RPA implementations are recomended (they also have gradients implemented).

## Simple examples

We cover two simple examples of low-scaling post-HF calculations: SOS-MP2 water MD, and titania RPA@PBE0 cell optimization. These examples are meant to be representative of the available options, and the parameters to pay attention to. They are also relatively afforable (although the titania input might require a few nodes on an HPC). To maintain a resonable cost, the systems are small and with low quality basis sets. In practice, such small systems should be simulated with the standard quartic scaling implementation. Moreover, it is never recommended to run post-HF calculations with double-zeta quality basis sets. All files required to run these examples can be found here.

### SOS-MP2 liquid water

This example is a short 3 steps MD liquid water calculation (32 molecules). The RI-MP2 optimized cc-TZ and RI_TZ basis sets are used. Whenever available, one should use pre-optimized RI basis sets for better performance. By default, the RI basis is generated on the fly. While it yields accurate results, the automatically generated RI basis sets are larger, and therefore less efficient. Using a short range RI metric improves efficiency by introducing sparsity. However, in somewhat sparse systems like water, the overlap metric lacks accuracy (Bussy2023). Using the truncated Coulomb operator with a cutoff radius of 1.5-2.0 Angstroms as a RI metric is safer.

All post-HF calculations need an accurate and well converged SCF as a starting point. In this case, a pure Hartree-Fock calculation is performed. Note the use of the truncated Coulomb operator as HFX potential, with cutoff radius < L/2, as always required in periodic HF. Here, the value of 4.5 Angstroms is somewhat small, and simulating a larger water box would be advised.

A numerical Laplace transform is used to reduce the scaling of the method. The integral is replaced by a weighted sum over a few grid points, placed according to the MINIMAX algortihm. Typically, using 6-8 points is enough. CP2K can go up to 20 points, although more points do not necessarily increase accuracy (and sometime bring numerical instability). The MINIMAX quadrature depends on the band gap of the system and the total spread of the SCF eigenvalues. If the seleceted number of grid point is innapropriate, a warning is issued.

    &GLOBAL
      PROJECT water32
      RUN_TYPE MD
      PRINT_LEVEL MEDIUM
    &END GLOBAL
    &MOTION
      &MD
        STEPS 3
      &END MD
    &END MOTION
    &FORCE_EVAL
      &DFT
        !cc-TZ: RI-MP2 optimized basis sets
        BASIS_SET_FILE_NAME BASIS_RI_cc-TZ
        POTENTIAL_FILE_NAME POTENTIAL
        SORT_BASIS EXP

        &MGRID
          CUTOFF 600
          REL_CUTOFF 50
          NGRIDS 5
        &END MGRID

        &SCF
          SCF_GUESS RESTART
          EPS_SCF 1.0E-6
          MAX_SCF 40
        &END SCF

        &XC
          &XC_FUNCTIONAL NONE
          &END XC_FUNCTIONAL
          &HF
            FRACTION 1.0
            &INTERACTION_POTENTIAL
              !TC potential with cutoff < L/2 for periodic HFX
              POTENTIAL_TYPE TRUNCATED
              CUTOFF_RADIUS 4.5
            &END INTERACTION_POTENTIAL
            &MEMORY
              !maximum memory allocated to HFX ERI storage, per MPI rank
              !the optimal number depends on the specifics of the computer
              MAX_MEMORY 4000
            &END MEMORY
          &END HF
          &WF_CORRELATION
            !explicit opposite-spin scaling
            SCALE_S 1.3
            &RI_SOS_MP2
              !MINIMAX quadrature by default
              QUADRATURE_POINTS 6
            &END RI_SOS_MP2
            !Enabling low-scaling
            &LOW_SCALING
              MEMORY_CUT 3
            &END LOW_SCALING
            &RI
              &RI_METRIC
                !Short range RI metric for SOS-MP2
                POTENTIAL_TYPE TRUNCATED
                CUTOFF_RADIUS 1.5
              &END RI_METRIC
            &END RI
            &INTEGRALS
              ERI_METHOD GPW
              &WFC_GPW
                !Safe yet faster than default values
                CUTOFF 200
                REL_CUTOFF 40
              &END  WFC_GPW
            &END INTEGRALS
          &END WF_CORRELATION
        &END XC
      &END DFT
      &SUBSYS
        &CELL
          ABC 9.8528 9.8528 9.8528
        &END CELL
        &TOPOLOGY
          COORD_FILE_FORMAT XYZ
          COORD_FILE_NAME ./H2O-32.xyz
        &END TOPOLOGY
        &KIND H
          BASIS_SET cc-DZ
          BASIS_SET RI_AUX RI_DZ
          POTENTIAL GTH-HF
        &END KIND
        &KIND O
          BASIS_SET cc-DZ
          BASIS_SET RI_AUX RI_DZ
          POTENTIAL GTH-HF
        &END KIND
      &END SUBSYS
    &END FORCE_EVAL

### RPA bulk titania

This second example is 2 steps of a cell optimization of titania at the RPA@PBE0 level (72 atoms). All low-scaling post-HF methods implemented in CP2K can be used with the ADMM (ADMM2 flavor) approximation to speed up the Hartree-Fock calculations for the SCF, the response forces, and EXX (in case of RPA). Here, the admm-dzp basis is used as an auxiliary ADMM basis set, and ccGRB-D as the primary basis. The RI basis is generated automatically.

The overlap RI metric is used in this example. For dense systems, this choice of metric offers a good tradeoff between accuracy and efficiency. Note that the `&HF` input section is repeated; it appears once for the SCF specification, and once for the RPA. In RPA@PBE0, the SCF is first converged at the PBE0 level of theory. Then, the RPA correlation energy is computed using the PBE0 orbitals and eigenvalues, and finally, the PBE0 exchange-correlation energy is replaced by the exact exchange energy (EXX) calculated with the same orbitals. The `&HF` section in RPA refers to that last step. If taken to be the same as the SCF (except for `FRACTION` and `MEMORY`), the stored ERIs can be reused and rescaled instead of recomputed, thus saving precious computational time. The `ADMM` keyword in RPA specifies that the EXX should be calculated within the ADMM approximation. Note the truncated Coulomb potential used for HFX, and its cutoff radius of 4.5 Angstroms (< L/2), necessary for periodic HFX calculations. This cutoff radius is on the shorter size, and using a larger simulation cell would be advised.

The `MINIMAX_QUADRATURE` is explicitely requested in the input. This is a necessary step for efficiency, as the default grid follows a Clenshaw-Curtis scheme, which requires many more integration points for the same accuracy.

    &GLOBAL
      PROJECT TiO2
      RUN_TYPE CELL_OPT
      PRINT_LEVEL MEDIUM
    &END GLOBAL
    &MOTION
      &CELL_OPT
        MAX_ITER 2
      &END CELL_OPT
    &END MOTION
    &FORCE_EVAL
      STRESS_TENSOR ANALYTICAL
      &DFT
        BASIS_SET_FILE_NAME BASIS_ccGRB_UZH
        BASIS_SET_FILE_NAME BASIS_ADMM_UZH
        POTENTIAL_FILE_NAME POTENTIAL_UZH
        SORT_BASIS EXP
        !automatically generated RI basis set
        AUTO_BASIS RI_AUX SMALL

        !enabling the ADMM2 approximation
        &AUXILIARY_DENSITY_MATRIX_METHOD
          METHOD BASIS_PROJECTION
          ADMM_PURIFICATION_METHOD NONE
          EXCH_CORRECTION_FUNC PBEX
        &END AUXILIARY_DENSITY_MATRIX_METHOD

        &MGRID
          CUTOFF 600
          REL_CUTOFF 50
          NGRIDS 5
        &END MGRID

        &SCF
          EPS_SCF 1.0E-6
          MAX_SCF 20
          &OT
            PRECONDITIONER FULL_ALL
            MINIMIZER DIIS
          &END OT
          &OUTER_SCF
            MAX_SCF 5
            EPS_SCF 1.0E-6
          &END OUTER_SCF
        &END SCF

        &XC
          &XC_FUNCTIONAL
            &PBE
              SCALE_X 0.75
            &END PBE
          &END XC_FUNCTIONAL
          &HF
            FRACTION 0.25
            !always use a short range potential < L/2 in periodic HFX
            &INTERACTION_POTENTIAL
              POTENTIAL_TYPE TRUNCATED
              CUTOFF_RADIUS 4.5
            &END INTERACTION_POTENTIAL
            &MEMORY
              !maximum memory allocated to HFX ERI storage, per MPI rank
              !the optimal number depends on the specifics of the computer
              MAX_MEMORY 4000
            &END MEMORY
          &END HF
          &WF_CORRELATION
            &RI_RPA
              MINIMAX_QUADRATURE
              QUADRATURE_POINTS 6
              !calculate EXX with ADMM, using the same HF section as in SCF
              !so that the integrals can be resued without recomputing
              ADMM
              &HF
                FRACTION 1.0
                &INTERACTION_POTENTIAL
                  POTENTIAL_TYPE TRUNCATED
                  CUTOFF_RADIUS 4.5
                &END INTERACTION_POTENTIAL
              &END HF
            &END RI_RPA
            !enabling low-scaling
            &LOW_SCALING
              MEMORY_CUT 3
            &END LOW_SCALING
            &RI
              !overlap RI metric is appropriate for dense systems
              &RI_METRIC
                POTENTIAL_TYPE IDENTITY
              &END RI_METRIC
            &END RI
            &INTEGRALS
              ERI_METHOD GPW
              !Safe yet faster than default values
              &WFC_GPW
                CUTOFF 200
                REL_CUTOFF 40
              &END WFC_GPW
            &END INTEGRALS
          &END WF_CORRELATION
        &END XC
      &END DFT
      &SUBSYS
        &CELL
          ABC 9.330 9.330 9.107
        &END CELL
        &TOPOLOGY
          COORD_FILE_FORMAT XYZ
          COORD_FILE_NAME ./TiO2.xyz
        &END TOPOLOGY
        &KIND Ti
          BASIS_SET ccGRB-D-q12
          BASIS_SET AUX_FIT admm-dz-q12
          POTENTIAL GTH-PBE0-q12
        &END KIND
        &KIND O
          BASIS_SET ccGRB-D-q6
          BASIS_SET AUX_FIT admm-dz-q6
          POTENTIAL GTH-PBE0-q6
        &END KIND
      &END SUBSYS
    &END FORCE_EVAL

## Important input parameters

There are a few importnat input parameters for low-scaling post-HF calculations:

- `SORT_BASIS`, in the DFT input section, should be set to EXP. This way, basis elements are sorted according to their exponents, leading to increased sparsity and performance.

- `EPS_FILTER`, in the LOW_SCALING input section. The bottleneck of the low-scaling post-HF calculations consists of sparse tensor contractions. This parameter is the threshold for block sparsity during calculations, with a safe default of \(1.0\times10^{-9}\). A looser threshold might significantly speedup calculations (not recommended to go higher than \(1.0\times10^{-8}\)).

- `MEMORY_CUT`, in the LOW_SCALING input section. This keyword influences the batching strategy for large tensor contractions. In post-HF, some tensor contractions are done in multiple steps, with large intermediate results. Storing these intermediates can lead to memory shortage. The value of `MEMORY_CUT` (default of 5) indicates how large tensors are split into batches, such that smaller intermediates can be stored. Note that MEMORY_CUT 3 does not mean that the total memory consumption of the program is divided by three (initial and final tensors are not affected, as well as the data concerning the rest of the calculation). A higher value reduces the memory footprint, but leads to performance overheads.

- `RI METRIC`: the choice of RI metric is crucial for performance in periodic calculations. A
  shorter ranged RI metric will generally be more efficient, while a longer ranged metric (e.g. TC
  with a cutoff radius of 1.5 Angstrom) can be more accurate. See Bussy2023 for a full discussion.