====== Hartree-Fock exchange ======
The purpose of this section is to explain how to perform hybrid functional calculations (or Hartree-Fock exchange, HFX) with CP2K in condensed phase systems. It is based on the developments described in [[doi>10.1021/ct900494g]] and [[doi>10.1063/1.2931945]], and its efficient extension (ADMM) described in [[doi>10.1021/ct1002225]].
Hartree-Fock exchange in CP2K is based on four center electron repulsion integrals (ERI), these are computed with an external library ([[http://sourceforge.net/projects/libint/|libint]]).
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.
This tutorial is also based on the ''GPW'' method and the calculations will be run using the GTH pseudopotentials.
Tutorial re-adapted from the [[https://www.cp2k.org/exercises:2016_summer_school:hfx| Tutorial on HFX and ADMM 2016 at KCL]]
For more info see also the slides from Joost VandeVondele [[pdf>https://www.cecam.org/upload/talk/presentation_5766.pdf| presentation 2011]], Sanliang Ling [[pdf>https://www.cp2k.org/_media/events:2016_user_meeting:cp2k-uk-2016-ling.pdf | presentation 2016]] and
Matt Watkins [[https://mattatlincoln.github.io/cp2k_workshop_2017 | presentation 2017]]
===== Truncated Coulomb operator =====
To enable HFX in the condensed phase 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 ====
We will perform a single point PBE-D3 calculation to generate an initial wavefunction (wfn) restart. HFX calculations benefit from this.
Below there is an input file for a single point (energy) calculation with CP2K. Copy this input file to a directory where you will run the calculation, and in the same directory extract the file {{ :exercises:2017_uzh_cp2k-tutorial:water.tar.gz |}}, which contains the structure of a box containing 64 water molecules.
&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
! amount of information printed to output
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
! 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-7
! 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
&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 ATOMIC ! can be used to RESTART an interrupted calculation
MAX_SCF 30
EPS_SCF 1.0E-6 ! accuracy of the SCF procedure, for OT typically 1.0E-6 - 1.0E-7, for diagonalization may have to be smaller
&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 PBE functional
&XC_FUNCTIONAL
&PBE
&END
&END XC_FUNCTIONAL
! 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 PBE
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
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]'' (the experimental water band-gap is of about 6.9 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. One can also find more functionals like meta-GGA etc [[http://octopus-code.org/wiki/Libxc|library of exchange-correlation functionals]] ([[https://manual.cp2k.org/cp2k-2_6-branch/CP2K_INPUT/ATOM/METHOD/XC/XC_FUNCTIONAL/LIBXC.html#list_SECTION_PARAMETERS|&LIBXC]]), and mix with different ratio of HF contribution to formulate various hybrid functionals.
Change the input to:
* ''SCF_GUESS RESTART'' in the ''&SCF'' section
* ''WFN_RESTART_FILE_NAME WATER-RESTART-GGA.wfn'' in the ''&DFT'' section
And replace the ''&XC...'' section with the following one:
! 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
! Important to improve scaling from O(N^4) to O(N)
&SCREENING
! important parameter to get stable HFX calcs (contributions to hfx smaller than EPS_SCHWARZ are not considered)
EPS_SCHWARZ 1.0E-6
! needs a good (GGA) initial guess
! screening on the product between maximum of density matrix elements and ERI
SCREEN_ON_INITIAL_P TRUE
&END
&INTERACTION_POTENTIAL
! for condensed phase systems
POTENTIAL_TYPE TRUNCATED
! should be less than half 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_PGF_ORB'' controls the sparse pattern of the overlap matrix, any contribution of the density matrix smaller than ''EPS_PGF_ORB'' is treated as zero;
* ''EPS_SCHWARZ'' parameter for Schwarz screening, where contributions to HFX smaller than ''EPS_SCHWARZ'' treated as zero. Typical values range between $10^{-6} and 10^{-9}$;
* ''SCREEN_ON_INITIAL_P'' uses the density matrix initially available to SCF and starts to use screening on that;
* ''EPS_PGF_ORB'', ''EPS_SCHWARZ'', ''SCREEN_ON_INITIAL_P'' and ''EPS_FILTER_MATRIX'' are parameters to guarantee stable SCF.
* Fraction of exchange (''SCALE_X'', ''FRACTION'').
Have a look at the output, in the section where CP2K is performing the SCF loop, using the OT method.
The first iteration should look something like this:
----------------------------------- OT ---------------------------------------
Step Update method Time Convergence Total energy Change
------------------------------------------------------------------------------
Trace(PS): 512.0000000000
Electronic density on regular grids: -511.9999881229 0.0000118771
Core density on regular grids: 511.9999857029 -0.0000142971
Total charge density on r-space grids: -0.0000024200
Total charge density g-space grids: -0.0000024200
HFX_MEM_INFO| Est. max. program size before HFX [MiB]: 388
*** WARNING in hfx_energy_potential.F:600 :: The Kohn Sham matrix is not ***
*** 100% occupied. This may result in incorrect Hartree-Fock results. Try ***
*** to decrease EPS_PGF_ORB and EPS_FILTER_MATRIX in the QS section. For ***
*** more information see FAQ: https://www.cp2k.org/faq:hfx_eps_warning ***
HFX_MEM_INFO| Number of cart. primitive ERI's calculated: 11551170770
HFX_MEM_INFO| Number of sph. ERI's calculated: 2681710894
HFX_MEM_INFO| Number of sph. ERI's stored in-core: 1115419808
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 [MiB]: 1075
HFX_MEM_INFO| Whereof max-vals [MiB]: 44
HFX_MEM_INFO| Total compression factor ERI's RAM: 7.92
HFX_MEM_INFO| Total memory consumption ERI's disk [MiB]: 0
HFX_MEM_INFO| Total compression factor ERI's disk: 0.00
HFX_MEM_INFO| Size of density/Fock matrix [MiB]: 10
HFX_MEM_INFO| Size of buffers [MiB]: 2
HFX_MEM_INFO| Number of periodic image cells considered: 27
HFX_MEM_INFO| Est. max. program size after HFX [MiB]: 664
1 OT DIIS 0.80E-01 121.7 0.00096329 -1102.1377186463 -1.10E+03
You can see from the output the appearance of a ''WARNING'' that is related to the
tags ''EPS_PGF_ORB'' and ''EPS_FILTER_MATRIX''. It is recommended to check and understand every warning in
the output and if it is OK to discard it or not. This particular warning is related to the stability of the SCF cycle. For more information see [[https://www.cp2k.org/faq:hfx_eps_warning]].
**Questions:**
* Look at the output where the HOMO-LUMO gap has been printed out. 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 ?
* The most expensive part of the whole SCF cycle is represented by the first step, while the other steps are much faster. Why is that?
* __Optional__ You can check if the SCF cycle is stable by decreasing the values of ''EPS_PGF_ORB'', ''EPS_FILTER_MATRIX'' as well as ''EPS_SCHWARZ''. Re-run the calculation and see how this affects the ''ENERGY''.
* __Optional__ CP2K tries to store the ERI in-core and avoid to calculate them at each SCF step. Especially for large systems that can be run on large HCP machines it is important to run in-core operation and fit the calculations of the ERI into memory. To see the effect of not having enough memory on the time for the ''SCF'' cycle, modify the tag of ''MAX_MEMORY'' to 40 and rerun the calculation. How do the timings compare with those where ''MAX_MEMORY'' was larger?
It is recommended to run HFX calculations in-core, either by setting larger ''MAX_MEMORY'' or using more MPI processes.
===== Truncated Coulomb operator with long range correction =====
In HSE and other screened hybrid functionals, the 1/r exchange interaction potential can be separated into a long range and a short range part, where only the short range part is computed using exact exchange and the long-range part is computed using GGA exchange. This is what is done also when using a truncated Coulomb operator. 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 ''CUTOFF_RADIUS'' 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? Remember to use the same cutoff radius under the ''&PBE_HOLE_T_C_LR'' and ''&INTERACTION_POTENTIAL'' sections.
===== 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 to the input:
* insert an additional line ''BASIS_SET_FILE_NAME BASIS_ADMM''.
* 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
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 [[doi>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:
! Spin polarization, 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 [[doi>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 [[doi>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.
&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 1
MULTIPLICITY 2
&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