Hartree-Fock exchange

The purpose of this section is to explain how to compute hybrid functionals (or Hartree-Fock exchange, HFX) with CP2K in condensed phase systems. It is based on the developments described in 10.1021/ct900494g and 10.1063/1.2931945, and its efficient extension (ADMM) described in 10.1021/ct1002225.

Hartree-Fock exchange in CP2K is based on four center integrals, these are computed with an external library (libint). To do these exercises, CP2K must be linked to this library.

This approach has a computational cost that depends strongly on the nature of the basis, unless combined with ADMM (see below), do not use MOLOPT basis sets with HFX. We use basis sets from HFX_BASIS, which are suitable.

Truncated Coulomb operator

To enable HFX in the condensed phase (described at the Gamma point only), CP2K employs a truncated Coulomb operator for the exchange part. The physical picture is that we do not want to have 'self-exchange interactions' of an electron with its image in neighboring unit cells. As a rule of thumb, the maximum range (truncation radius) is L/2 where L is the smallest edge of the unit cell. The convergence of the exchange energy is exponential wrt. this radius. Typically, 5-6A provides good results, but this depends on the nature of the system, i.e. the band gap or the range of the maximally localized Wannier orbitals.

1st task : GGA restart wfn

Using the water input from the previous exercise, we will perform a single point GGA calculation to generate an initial wavefunction (wfn) restart. HFX calculations benefit from this.

Change the input to:

RUN_TYPE ENERGY
IOLEVEL MEDIUM
RESTART ON
comment section &EXT_RESTART

Run the input and rename the generated wfn file (WATER-RESTART.wfn) to WATER-RESTART-GGA.wfn. Also make a note of the HOMO - LUMO gap [eV]

2nd task: PBE0-D3 water

To do a hybrid calculation, we just change the &XC section. Various examples can be found in the regtests, but here we employ a section equivalent to PBE0-D3.

Change the input to:

SCF_GUESS RESTART
WFN_RESTART_FILE_NAME WATER-RESTART-GGA.wfn

And employ the following &XC section:

    ! specify the exchange and correlation treatment
    &XC
      ! use a PBE0 functional 
      &XC_FUNCTIONAL
       &PBE
         ! 75% GGA exchange
         SCALE_X 0.75
         ! 100% GGA correlation
         SCALE_C 1.0
       &END PBE
      &END XC_FUNCTIONAL
      &HF
        ! 25 % HFX exchange
        FRACTION 0.25
        &SCREENING
          ! important parameter to get stable HFX calcs
          EPS_SCHWARZ 1.0E-6
          ! needs a good (GGA) initial guess 
          SCREEN_ON_INITIAL_P TRUE
        &END
        &INTERACTION_POTENTIAL
          ! for condensed phase systems
          POTENTIAL_TYPE TRUNCATED
          ! should be less than halve the cell
          CUTOFF_RADIUS 6.0
          ! data file needed with the truncated operator
          T_C_G_DATA ./t_c_g.dat
        &END
        &MEMORY
          ! In MB per MPI rank.. use as much as need to get in-core operation
          MAX_MEMORY 4000
          EPS_STORAGE_SCALING 0.1
        &END
      &END
      ! adding Grimme's D3 correction (by default without C9 terms) 
      &VDW_POTENTIAL
         POTENTIAL_TYPE PAIR_POTENTIAL
         &PAIR_POTENTIAL
            PARAMETER_FILE_NAME dftd3.dat
            TYPE DFTD3
            REFERENCE_FUNCTIONAL PBE0
            R_CUTOFF [angstrom] 16
         &END
      &END VDW_POTENTIAL
    &END XC

Topics:

EPS_SCHWARZ, EPS_PGF_ORB, SCREEN_ON_INITIAL_P : parameters to guarantee stable SCF.
Fraction of exchange (SCALE_X, FRACTION).

Run this input, it should lead to an output like:

  HFX_MEM_INFO| Number of cart. primitive ERI's calculated:          20780449251
  HFX_MEM_INFO| Number of sph. ERI's calculated:                      4626861713
  HFX_MEM_INFO| Number of sph. ERI's stored in-core:                  1440639962
  HFX_MEM_INFO| Number of sph. ERI's stored on disk:                           0
  HFX_MEM_INFO| Number of sph. ERI's calculated on the fly:                    0
  HFX_MEM_INFO| Total memory consumption ERI's RAM [MB's]:                  1368
  HFX_MEM_INFO| Whereof max-vals [MB's]:                                      78
  HFX_MEM_INFO| Total compression factor ERI's RAM:                         8.03
  HFX_MEM_INFO| Total memory consumption ERI's disk [MB's]:                    0
  HFX_MEM_INFO| Total compression factor ERI's disk:                        0.00
  HFX_MEM_INFO| Size of density/Fock matrix [MB's]:                           14
  HFX_MEM_INFO| Size of buffers [MB's]:                                        2
  HFX_MEM_INFO| Number of periodic image cells considered:                    27
  HFX_MEM_INFO| Est. max. program size after HFX  [MB's]:                    243

     1 OT DIIS     0.80E-01   68.5     0.00096341     -1102.1375916328 -1.10E+03

  Trace(PS):                                  512.0000000000
  Electronic density on regular grids:       -511.9999999886        0.0000000114
  Core density on regular grids:              511.9999999836       -0.0000000164
  Total charge density on r-space grids:       -0.0000000051
  Total charge density g-space grids:          -0.0000000051

     2 OT DIIS     0.80E-01    4.1     0.00064929     -1102.1607534171 -2.32E-02

Topics:

in-core operation.
How to detect instabilities due to too aggressive screening.

Question: What is the HOMO-LUMO gap for this configuration ? How does this compare to the GGA result ? Adjust the fraction of exchange (modify the input in two places!) to 20% and/or 30%, how does this influence the gap ?

Truncated Coulomb operator with long range correction

Like in the HSE functional, the difference between the operator used for exchange and 1/r, can be accounted for by a special GGA exchange functional. Also for the truncated coulomb operator this is possible, and allows for xc functionals that embed very short range exchange operators only. This can be used to speedup the calculation, while retaining the benefits of HFX. The functional employed in this way smoothly goes from PBE to PBE0 as the range goes from 0 to Infinity.

3rd task

Add to the &XC_FUNCTIONAL section (i.e. in addition to &PBE) the following section:

       &PBE_HOLE_T_C_LR
         CUTOFF_RADIUS 2.5
         SCALE_X 0.25
       &END

and employ the same CUTOFF_RADIUS for the INTERACTION_POTENTIAL.

Rerun the single point energy calculation and note the band gap.

Is such a short range sufficient to have a sizable effect on the band gap ?
is HFX_MEM_INFO| Number of cart. primitive ERI's calculated very different for calculations with 2.5 and 6.0A truncation radius ?

Auxiliary Density Matrix Methods (ADMM)

ADMM is an approach to mitigate the cost of HFX for large basis sets. In particular, if MOLOPT basis sets are used, standard HFX becomes too expensive (CP2K can not deal efficiently with highly contracted AOs). In ADMM, an AUX_FIT_BASIS_SET is introduced, which is used to create an auxiliary density matrix (ADM) by projection. HFX is evaluated for this ADM, while the error introduced by using an ADM is corrected for with a GGA exchange functional.

4rd task : introduce ADMM

Make the following changes:

insert and additional line BASIS_SET_FILE_NAME BASIS_ADMM (and copy that file from cp2k/data as needed).
insert for each &KIND a line AUX_FIT_BASIS_SET cFIT3
insert a secion &AUXILIARY_DENSITY_MATRIX_METHOD

    ! use ADMM
    &AUXILIARY_DENSITY_MATRIX_METHOD
      ! recommended, i.e. use a smaller basis for HFX
      ! each kind will need an AUX_FIT_BASIS_SET.
      METHOD BASIS_PROJECTION
      ! recommended, this method is stable and allows for MD. 
      ! can be expensive for large systems
      ADMM_PURIFICATION_METHOD MO_DIAG
    &END

In this tutorial we combine ADMM using a very small basis set (cFIT3) with a small primary basis (DZVP-GTH), so gains are small at best (and the results not very accurate). ADMM is most useful with good quality primary basis sets, such as e.g. MOLOPTs

For CP2K 2.7, a new library of ADMM basis sets is available.

Run the input, what's the HOMO - LUMO gap ?

Chasing charge localization in liquid water

The combination of truncated exchange and ADMM results in the most effective way to run AIMD with hybrid functionals. In some systems the difference between GGA DFT and hybrids is very large. One such systems is liquid water after ionization (i.e. charge +1), where only with hybrids the expected species (OH radicals) are formed. See 10.1063/1.3664746.

5th task : ionized water

This task is optional, and can be done near the end if time is available.

adapt the admm input for water to reflect the ionized state:

    ! Charge and multiplicity
    LSD
    CHARGE 1
    MULTIPLICITY 2

because the system is electronically very difficult initially, we'll reduce the convergence threshold EPS_SCF 1.0E-5 (twice).

Note that the WFN_RESTART_FILE_NAME must point to a GGA calculation of the same charge and multiplicity (do this calculations first).

Run single point energy calculations varying the fraction of exchange from 0.25 to 0.50, does the mulliken spin population reproduce Fig. 2 in 10.1021/ct1002225 ?

For the fraction 0.5, run AIMD for about 50-100fs (if time permits), what happens with the water molecule on which the spin was localized ? Do you results agree with 10.1063/1.3664746 ?

Required files

No new files are required for this exercise. If you're stuck, you can use the following worked out examples.

water_pbe0_cheating.inp

&GLOBAL
  ! the project name is made part of most output files... useful to keep order 
  PROJECT WATER
  ! various runtypes (energy, geo_opt, etc.) available.
  RUN_TYPE ENERGY
  ! limit the runs to 30min
  WALLTIME 1800
  ! reduce the amount of IO
  IOLEVEL  MEDIUM 
&END GLOBAL

&FORCE_EVAL
  ! the electronic structure part of CP2K is named Quickstep
  METHOD Quickstep
  &DFT
    ! basis sets and pseudopotential files can be found in cp2k/data
    BASIS_SET_FILE_NAME HFX_BASIS
    POTENTIAL_FILE_NAME GTH_POTENTIALS            
    ! GGA restart to provide a good initial density matrix
    WFN_RESTART_FILE_NAME WATER-RESTART-GGA.wfn 

    ! Charge and multiplicity
    CHARGE 0
    MULTIPLICITY 1

    &MGRID
       ! PW cutoff ... depends on the element (basis) too small cutoffs lead to the eggbox effect.
       ! certain calculations (e.g. geometry optimization, vibrational frequencies,
       ! NPT and cell optimizations, need higher cutoffs)
       CUTOFF [Ry] 400 
    &END

    &QS
       ! use the GPW method (i.e. pseudopotential based calculations with the Gaussian and Plane Waves scheme).
       METHOD GPW 
       ! default threshold for numerics ~ roughly numerical accuracy of the total energy per electron,
       ! sets reasonable values for all other thresholds.
       EPS_DEFAULT 1.0E-10 
       ! used for MD, the method used to generate the initial guess.
       EXTRAPOLATION ASPC 
    &END

    &POISSON
       PERIODIC XYZ ! the default, gas phase systems should have 'NONE' and a wavelet solver
    &END

    &PRINT
       ! at the end of the SCF procedure generate cube files of the density
       &E_DENSITY_CUBE OFF
       &END E_DENSITY_CUBE
       ! compute eigenvalues and homo-lumo gap each 10nd MD step
       &MO_CUBES
          NLUMO 4
          NHOMO 4
          WRITE_CUBE .FALSE.
          &EACH
            MD 10
          &END
       &END
    &END

    ! use the OT METHOD for robust and efficient SCF, suitable for all non-metallic systems.
    &SCF                              
      SCF_GUESS RESTART ! can be used to RESTART an interrupted calculation
      MAX_SCF 30
      EPS_SCF 1.0E-6 ! accuracy of the SCF procedure typically 1.0E-6 - 1.0E-7
      &OT
        ! an accurate preconditioner suitable also for larger systems
        PRECONDITIONER FULL_SINGLE_INVERSE
        ! the most robust choice (DIIS might sometimes be faster, but not as stable).
        MINIMIZER DIIS
      &END OT
      &OUTER_SCF ! repeat the inner SCF cycle 10 times
        MAX_SCF 10
        EPS_SCF 1.0E-6 ! must match the above
      &END
      ! do not store the wfn during MD
      &PRINT
        &RESTART ON
        &END
      &END
    &END SCF

    ! specify the exchange and correlation treatment
    &XC
      ! use a PBE0 functional 
      &XC_FUNCTIONAL 
       &PBE
         ! 75% GGA exchange
         SCALE_X 0.75
         ! 100% GGA correlation
         SCALE_C 1.0
       &END PBE
      &END XC_FUNCTIONAL
      &HF
        ! 25 % HFX exchange
        FRACTION 0.25
        &SCREENING
          ! important parameter to get stable HFX calcs
          EPS_SCHWARZ 1.0E-6
          ! needs a good (GGA) initial guess 
          SCREEN_ON_INITIAL_P TRUE
        &END
        &INTERACTION_POTENTIAL
          ! for condensed phase systems
          POTENTIAL_TYPE TRUNCATED
          ! should be less than halve the cell
          CUTOFF_RADIUS 6.0
          ! data file needed with the truncated operator
          T_C_G_DATA ./t_c_g.dat
        &END
        &MEMORY
          ! In MB per MPI rank.. use as much as need to get in-core operation
          MAX_MEMORY 4000
          ! additional accuracy for storing compressed results
          EPS_STORAGE_SCALING 0.1
        &END
      &END
      ! adding Grimme's D3 correction (by default without C9 terms) 
      &VDW_POTENTIAL
         POTENTIAL_TYPE PAIR_POTENTIAL 
         &PAIR_POTENTIAL
            PARAMETER_FILE_NAME dftd3.dat
            TYPE DFTD3
            REFERENCE_FUNCTIONAL PBE0
            R_CUTOFF [angstrom] 16
         &END
      &END VDW_POTENTIAL
    &END XC
  &END DFT
 
  ! description of the system
  &SUBSYS
    &CELL 
      ! unit cells that are orthorhombic are more efficient with CP2K
      ABC [angstrom] 12.42 12.42 12.42
    &END CELL

    ! atom coordinates can be in the &COORD section,
    ! or provided as an external file.
    &TOPOLOGY
      COORD_FILE_NAME water.xyz
      COORD_FILE_FORMAT XYZ
    &END

    ! MOLOPT basis sets are fairly costly,
    ! but in the 'DZVP-MOLOPT-SR-GTH' available for all elements
    ! their contracted nature makes them suitable
    ! for condensed and gas phase systems alike.
    &KIND H                              
      BASIS_SET DZVP-GTH
      POTENTIAL GTH-PBE-q1             
    &END KIND
    &KIND O
      BASIS_SET DZVP-GTH
      POTENTIAL GTH-PBE-q6
    &END KIND
  &END SUBSYS
&END FORCE_EVAL

! how to propagate the system, selection via RUN_TYPE in the &GLOBAL section
&MOTION
 &GEO_OPT
   OPTIMIZER BFGS ! Good choice for 'small' systems (use LBFGS for large systems)
   MAX_ITER  100
   MAX_DR    [bohr] 0.003 ! adjust target as needed
   &BFGS
   &END
 &END
 &MD
   ENSEMBLE NVT  ! sampling the canonical ensemble, accurate properties might need NVE
   TEMPERATURE [K] 300
   TIMESTEP [fs] 0.5
   STEPS 1000
   # GLE thermostat as generated at http://epfl-cosmo.github.io/gle4md 
   # GLE provides an effective NVT sampling.
   &THERMOSTAT
     REGION MASSIVE
     TYPE GLE
     &GLE
       NDIM 5
       A_SCALE [ps^-1] 1.00
       A_LIST    1.859575861256e+2   2.726385349840e-1   1.152610045461e+1  -3.641457826260e+1   2.317337581602e+2
       A_LIST   -2.780952471206e-1   8.595159180871e-5   7.218904801765e-1  -1.984453934386e-1   4.240925758342e-1
       A_LIST   -1.482580813121e+1  -7.218904801765e-1   1.359090212128e+0   5.149889628035e+0  -9.994926845099e+0
       A_LIST   -1.037218912688e+1   1.984453934386e-1  -5.149889628035e+0   2.666191089117e+1   1.150771549531e+1
       A_LIST    2.180134636042e+2  -4.240925758342e-1   9.994926845099e+0  -1.150771549531e+1   3.095839456559e+2
     &END GLE
   &END THERMOSTAT
 &END
  &PRINT
   &TRAJECTORY
     &EACH
       MD 1
     &END EACH
   &END TRAJECTORY
   &VELOCITIES OFF
   &END VELOCITIES
   &FORCES OFF
   &END FORCES
   &RESTART_HISTORY
     &EACH
       MD 500
     &END EACH
   &END RESTART_HISTORY
   &RESTART
     BACKUP_COPIES 3
     &EACH
       MD 1
     &END EACH
   &END RESTART
  &END PRINT
&END

CP2K Open Source Molecular Dynamics

Table of Contents