# Møller–Plesset Perturbation Theory

The MP2 method computes

\[ 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 DelBen2012 and DelBen2013. CP2K also implements double-hybrid functionals.

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

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

    &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

    &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 increase the memory demands but reduce communication (default: -1, automatic determination)
        NUMBER_INTEGRATION_GROUPS 2
      &END
    &END

## RI-MP2 Gradient Calculations

Analytical gradients (nuclear gradients, stress tensors) are available for RI-MP2 and related double-hybrid functionals.DelBen2015b,Stein2022 HF calculations may also be accelerated by the ADMM approximation which is recommended if very diffuse basis functions are included. The computational costs of the MP2 part is 2-3 times as high as for energy-only calculations. A typical input file of a RI-MP2 calculation using ADMM-accelerated HF is

    &GLOBAL
      PRINT_LEVEL LOW
      PROJECT example
      RUN_TYPE FORCE_EVAL
    &END GLOBAL

    &FORCE_EVAL
      METHOD Quickstep
      &DFT
        # This file contains basis sets for common elements (H, C, N, O; for F, Al, Cl, Si, P, S only TZ basis sets)
        BASIS_SET_FILE_NAME BASIS_RI_cc-TZ
        BASIS_SET_FILE_NAME BASIS_ADMM
        POTENTIAL_FILE_NAME POTENTIAL
        # Plug in your preconverged wfn file here
        WFN_RESTART_FILE_NAME ./example_HF.wfn
        # Use if no RI basis set is available (possible values: SMALL, MEDIUM, LARGE, HUGE)
        AUTO_BASIS RI_AUX LARGE

        &AUXILIARY_DENSITY_MATRIX_METHOD
          # Purification is not implemented
          ADMM_PURIFICATION_METHOD NONE
          # Try different options (check with reference values or use the default)
          EXCH_CORRECTION_FUNC PBEX
          # other methods are not implemented
          METHOD BASIS_PROJECTION
        &END AUXILIARY_DENSITY_MATRIX_METHOD
        &MGRID
          # Adjust as usual
          CUTOFF 600
          REL_CUTOFF 50
        &END MGRID
        &POISSON
          PERIODIC XYZ
          POISSON_SOLVER WAVELET
        &END POISSON
        &QS
          # Choose as tight as necessary
          EPS_DEFAULT 1.0E-10
          METHOD GPW
        &END QS
        &SCF
          # Choose as tight as necessary
          EPS_SCF 1.0E-8
          MAX_SCF 100
          SCF_GUESS RESTART
        &END SCF
        &XC
          &HF
            FRACTION 1.0000000
            &INTERACTION_POTENTIAL
              # Adjust the cutoff radius according to your cell
              CUTOFF_RADIUS 1.5
              POTENTIAL_TYPE TRUNCATED
              T_C_G_DATA t_c_g.dat
            &END INTERACTION_POTENTIAL
            &SCREENING
              # Tune these parameters
              EPS_SCHWARZ 1.0E-10
              EPS_SCHWARZ_FORCES 1.0E-5
              SCREEN_ON_INITIAL_P .FALSE.
            &END SCREENING
          &END HF
          &WF_CORRELATION
            MEMORY 500
            NUMBER_PROC 1
            &CANONICAL_GRADIENTS
              # The default is usually good enough
              EPS_CANONICAL 1E-6
              # This is the option which should always work
              # Try to set it to .TRUE. (may segfault)
              FREE_HFX_BUFFER .FALSE.
              &CPHF
                # Choose as tight as necessary
                EPS_CONV 1.0E-6
                # Smaller values of EPS_CONV may require more iterations
                MAX_ITER 10
              &END CPHF
            &END CANONICAL_GRADIENTS
            &INTEGRALS
              &WFC_GPW
                # Adjust these parameters to your needs of accuracy
                CUTOFF 200
                EPS_FILTER 1.0E-12
                EPS_GRID 1.0E-8
                REL_CUTOFF 50
              &END WFC_GPW
            &END INTEGRALS
            &RI_MP2
              # Determine an automatic block size, if memory is low, set it to 1
              BLOCK_SIZE -1
            &END RI_MP2
          &END WF_CORRELATION
          &XC_FUNCTIONAL NONE
          &END XC_FUNCTIONAL
        &END XC
      &END DFT
      &PRINT
        &FORCES
        &END FORCES
      &END PRINT
      &SUBSYS
        # Adjust cell information and coordinations as necessary
        &CELL
          ABC [angstrom] 5.0 5.0 5.0
          PERIODIC XYZ
        &END CELL
        &COORD
          O       0.000000    0.000000    -0.211000
          H       0.000000   -0.844000     0.495000
          H       0.000000    0.744000     0.495000
        &END COORD
        &KIND H
          # orbital and RI basis set should match
          BASIS_SET cc-TZ
          BASIS_SET RI_AUX RI_TZ
          # Use the largest affordable one
          BASIS_SET AUX_FIT cpFIT3
          POTENTIAL GTH-HF-q1
        &END KIND
        &KIND O
          BASIS_SET cc-TZ
          BASIS_SET RI_AUX RI_TZ
          BASIS_SET AUX_FIT cpFIT3
          POTENTIAL GTH-HF-q6
        &END KIND
        &TOPOLOGY
          &CENTER_COORDINATES
          &END CENTER_COORDINATES
        &END TOPOLOGY
      &END SUBSYS
    &END FORCE_EVAL

A few more words are required regarding systems with many degenerate occupied orbital pairs as they are found in structures with many symmetry-equivalent atoms. In that case, non-diagonal elements of the density matrix of almost degenerate occupied-orbital pairs have to be calculated explicitly for numerical reason. This is significantly more expensive than for non-degenerate orbital pairs. This behavior can be tuned using the EPS_CANONICAL keyword.

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