The purpose of this section is to explain how to compute the energy of a molecular orbital from GW for molecules and condensed phase systems with CP2K. In DFT, the energy of a molecular orbital corresponds to an eigenvalue of the Kohn-Sham matrix. In GW, the procedure for getting the level energies is to first perform a DFT calculation (commonly with the PBE or PBE0 functional) to get the molecular orbital wavefunctions and then compute a new GW energy for the molecular orbitals of interest. For an introduction into the concept of GW, please read Sec. II and the introduction to Sec. III in 10.1103/PhysRevB.87.235132 [1].
The GW implementation in CP2K is based on the developments described in 10.1021/acs.jctc.6b00380 [2] which is very similar to the GW implementation in FHI-aims 10.1088/1367-2630/14/5/053020 [3], Fiesta 10.1103/PhysRevB.83.115103 [4] and molgw 10.1016/j.cpc.2016.06.019 [5]. The computational cost of GW is comparable to RPA and MP2 total energy calculations and therefore high. The computational cost of a canonical GW implementation grows as $N^4$ with the system size $N$, while the memory scales as $N^3$ with the system size. The basis set convergence of GW is slow and therefore has to be carefully examined.
In this tutorial, GW values from the GW100 benchmark set 10.1021/acs.jctc.5b00453 [6] are reproduced (Section 1), the basis set extrapolation is examined for the water molecule (Section 2), an input for large-scale, parallel calculations is given (Section 3), periodic GW calculations are presented 10.1103/PhysRevB.95.235123 [7] (Section 4), and cubic-scaling GW calculations are shown which can be more efficient for systems with hundrets of atoms (Section 5).
Since the calculations are rather small, please use a single MPI rank for the calculation:
mpirun -n 1 cp2k.popt H2O_GW100.inp | tee cp2k.out
See below the input for a G0W0@PBE calculation of the water molecule in a def2-QZVP basis: A PBE calculation is used for computing the molecular orbitals which can be seen from the keyword “XC_FUNCTIONAL PBE”. The input parameters for G0W0 are commented below. While the calculation is running, you can look up the G0W0@PBE value for the highest occupied molecular orbital (HOMO) and the lowest unoccupied molecular orbital (LUMO) in 10.1021/acs.jctc.5b00453 [6] [Table 2 and 3 (column AIMS-P16/TM-no RI, molecule index 76) which should be -11.97 eV and 2.37 eV for HOMO and LUMO, respectively]. CP2K should be able to exactly reproduce these values. In the output of CP2K, the G0W0@PBE results are listed after the SCF after the headline GW quasiparticle energies.
The G0W0@PBE HOMO value is not in good agreement with the experimental ionization potential of water (12.62 eV). A possible explanation is that PBE may not be a good starting point for G0W0 calculations for molecules in the gas phase, see e.g. 10.1021/ct300835h [8]. For checking the basis set convergence, we refer to a detailed analysis in 10.1021/acs.jctc.6b00380 [2] in Fig. 2 for benzene. An extensive table of basis set extrapolated GW levels can be found in 10.1021/acs.jctc.5b00453 [6] (column EXTRA).
&FORCE_EVAL
METHOD Quickstep
&DFT
BASIS_SET_FILE_NAME BASIS_def2_QZVP_RI_ALL
POTENTIAL_FILE_NAME POTENTIAL
&MGRID
CUTOFF 400
REL_CUTOFF 50
&END MGRID
&QS
! all electron calculation since GW100 is all-electron test
METHOD GAPW
&END QS
&POISSON
PERIODIC NONE
PSOLVER MT
&END
&SCF
EPS_SCF 1.0E-6
SCF_GUESS ATOMIC
MAX_SCF 200
&END SCF
&XC
&XC_FUNCTIONAL PBE
&END XC_FUNCTIONAL
! GW is part of the WF_CORRELATION section
&WF_CORRELATION
! RPA is used to compute the density response function
METHOD RI_RPA_GPW
! Use Obara-Saika integrals instead of GPW integrals
! since OS is much faster
ERI_METHOD OS
&RI_RPA
! use 100 quadrature points to perform the
! frequency integration in GW
RPA_NUM_QUAD_POINTS 100
! SIZE_FREQ_INTEG_GROUP is a group size for parallelization and
! should be increased for large calculations to prevent out of memory.
! maximum for SIZE_FREQ_INTEG_GROUP is the number of MPI tasks
SIZE_FREQ_INTEG_GROUP 1
GW
&RI_G0W0
! compute the G0W0@PBE energy of HOMO-9,
! HOMO-8, ... , HOMO-1, HOMO
CORR_OCC 10
! compute the G0W0@PBE energy of LUMO,
! LUMO+1, ... , LUMO+20
CORR_VIRT 20
! fit a Pade approximant to the correlation self-energy
! as function of imaginary frequency. this has been done
! in the GW100 benchmark set and turned out to be reliable
ANALYTIC_CONTINUATION PADE
! for solving the quasiparticle equation, the Newton method
! is used as in the GW100 benchmark
CROSSING_SEARCH NEWTON
! use the RI approximation for the exchange part of the self-energy
RI_SIGMA_X
&END RI_G0W0
&END RI_RPA
! NUMBER_PROC is a group size for parallelization and should
! be increased for large calculations
NUMBER_PROC 1
&END
&END XC
&END DFT
&SUBSYS
&CELL
ABC 10.0 10.0 10.0
PERIODIC NONE
&END CELL
&COORD
O 0.0000 0.0000 0.0000
H 0.7571 0.0000 0.5861
H -0.7571 0.0000 0.5861
&END COORD
&TOPOLOGY
&CENTER_COORDINATES
&END
&END TOPOLOGY
&KIND H
! def2-QZVP is the basis which has been used in the GW100 paper
BASIS_SET def2-QZVP
! just use a very large RI basis to ensure excellent
! convergence with respect to the RI basis
RI_AUX_BASIS RI-5Z
POTENTIAL ALL
&END KIND
&KIND O
BASIS_SET def2-QZVP
RI_AUX_BASIS RI-5Z
POTENTIAL ALL
&END KIND
&END SUBSYS
&END FORCE_EVAL
&GLOBAL
RUN_TYPE ENERGY
PROJECT ALL_ELEC
PRINT_LEVEL MEDIUM
&END GLOBAL
In this section, the slow basis set convergence of GW calculations is examined. We compute the G0W0@PBE HOMO and LUMO level of the water molecule with Dunning's cc-DZVP, cc-TZVP, cc-QZVP and cc-5ZVP all-electron basis sets and extrapolate these values to the complete basis set limit. To do so, download the cc basis sets cc_basis_h2o.tar which has been taken from the EMSL basis set database. Run the input from Sec. 1 using the cc-DZVP to cc-5ZVP basis set (in total four calculations) by replacing the basis sets:
BASIS_SET_FILE_NAME BASIS_def2_QZVP_RI_ALL BASIS_SET_FILE_NAME ./BASIS_H2O
&KIND H BASIS_SET cc-DZVP-all RI_AUX_BASIS RI-5Z POTENTIAL ALL &END KIND &KIND O BASIS_SET cc-DZVP-all RI_AUX_BASIS RI-5Z POTENTIAL ALL &END KIND
Employ the RI-5Z basis set as RI-basis which ensures excellent convergence for the RI basis. In practice, smaller RI basis sets can be used from the EMSL database (just check the convergence with respect to the RI basis by using smaller and larger RI basis sets).
The results for the G0W0@PBE HOMO and LUMO from CP2K should be as follows:
Basis set G0W0@PBE HOMO (eV) G0W0@PBE LUMO (eV) N_basis N_card cc-DZVP -12.480 4.770 23 2 cc-TZVP -12.417 3.424 57 3 cc-QZVP -12.180 2.773 114 4 cc-5ZVP -12.108 2.088 200 5 Extrapolation using cc-TZVP to cc-5ZVP with 1/N_card^3 -12.02 +/- 0.01 1.90 +/- 0.29 with 1/N_basis -11.97 +/- 0.02 1.71 +/- 0.29 GW100 -12.05 2.01
For the extrapolation, two schemes have been used as described in the GW100 paper and its supporting information 10.1021/acs.jctc.5b00453 [6]. The first scheme employs a linear fit on the HOMO or LUMO values when they are plotted against the inverse cardinal number $N_\text{card}$ of the basis set while the second scheme extrapolates versus the inverse number of basis functions $N_\text{basis}$ which can be computed as sum of the number of occupied orbitals and the number of virtual orbitals as printed in RI_INFO in the output. You can check the extrapolation from the table above with your tool of choice.
The basis set extrapolated values from the table above deviate from the values reported in the GW100 paper 10.1021/acs.jctc.5b00453 [6], probably because only two basis sets (def2-TZVP, def2-QZVP) have been used in 10.1021/acs.jctc.5b00453 [6] for the extrapolation. The extrapolation for the LUMO is not working well because one would need much more diffuse functions to represent unbound electronic levels (with positive energy).
Often, the HOMO-LUMO gap is of interest. In this case, augmented basis sets (e.g. from the EMSL database) can offer an alternative for very fast basis set convergence, see also Fig. 2b in 10.1021/acs.jctc.6b00380 [2].
An exemplary input for a parallel calculation can be found in the supporting information of 10.1021/acs.jctc.6b00380 [2]. The emphasis is on the parameters SIZE_FREQ_INTEG_GROUP and NUMBER_PROC which should be increased for larger calculations. In case of a too small number, the code will crash due to out of memory while a too large number results in slow speed. Typically, one starts for large-scale calculations from a small molecule. When increasing the system size, the parameters SIZE_FREQ_INTEG_GROUP and NUMBER_PROC should be both increased to avoid a crash due to out of memory. The maximum number for both parameters is the number of MPI tasks. Also, the number of nodes should be increased with $N^3_\text{atoms}$ to account for the scaling of the memory of GW.
For periodic GW calculations, a special correction scheme is necessary. A similar problem is appearing in Hartree-Fock calculations. In HF, an easy way out is given by the truncation of the Coulomb operator which works due to the convenient form of the HF equations. GW does not exhibit this convenient form and therefore, this truncation does not work for GW calculations. The theory why a correction is necessary and the correction scheme is described in 10.1103/PhysRevB.95.235123 [7].
The basis can be found in lih_basis_ri.tar. Then run the calculation as listed below. The computed gap of LiH with the input from below should be 6.05 eV which is still significantly off from the converged result of 4.7 eV from 10.1103/PhysRevB.95.235123 [7]. Here, a basis set extrapolation (as shown above) and a larger supercell are necessary to get closer to the result of 4.7 eV. Please redo the calculation without the flag PERIODIC in the GW section and see that the resulting gap of 11.25 eV is much more off than the gap with correction.
&FORCE_EVAL
METHOD Quickstep
&DFT
BASIS_SET_FILE_NAME ./LiH_BASIS_RI
POTENTIAL_FILE_NAME POTENTIAL
&MGRID
CUTOFF 600
REL_CUTOFF 60
&END MGRID
&QS
METHOD GPW
EPS_DEFAULT 1.0E-15
EPS_PGF_ORB 1.0E-20
EPS_FILTER_MATRIX 0.0e0
&END QS
&SCF
EPS_SCF 1.0E-6
MAX_SCF 100
&END SCF
&XC
&XC_FUNCTIONAL PBE
&END XC_FUNCTIONAL
&WF_CORRELATION
METHOD RI_RPA_GPW
&RI_RPA
RPA_NUM_QUAD_POINTS 100
GW
&RI_G0W0
CORR_OCC 5
CORR_VIRT 5
! activate the periodic correction
PERIODIC
ANALYTIC_CONTINUATION PADE
CROSSING_SEARCH NEWTON
&END RI_G0W0
! HF calculation for the exchange part of the self-energy
! Here, the truncation of the Coulomb operator works
&HF
&SCREENING
! for other materials, a smaller EPS_SCHWARZ might be necessary
EPS_SCHWARZ 1.0E-6
SCREEN_ON_INITIAL_P TRUE
&END
&INTERACTION_POTENTIAL
POTENTIAL_TYPE TRUNCATED
! the truncation radius is half the cell size
CUTOFF_RADIUS 2.00
T_C_G_DATA t_c_g.dat
&END
&MEMORY
MAX_MEMORY 0
&END
&END
&END RI_RPA
NUMBER_PROC 1
&END
&END XC
&END DFT
&SUBSYS
&CELL
ABC 4.084 4.084 4.084
&END CELL
&COORD
Li 0 0 0
Li 2.042 2.042 0
Li 2.042 0 2.042
Li 0 2.042 2.042
H 0 2.042 0
H 0 0 2.042
H 2.042 0 0
H 2.042 2.042 2.042
&END COORD
&KIND H
BASIS_SET cc-DZVP-GTH
RI_AUX_BASIS_SET RI_DZ_opt_basis
POTENTIAL GTH-PBE-q1
&END KIND
&KIND Li
BASIS_SET cc-DZVP-GTH
RI_AUX_BASIS_SET RI_DZ_opt_basis
POTENTIAL GTH-PBE-q3
&END KIND
&END SUBSYS
&END FORCE_EVAL
&GLOBAL
PROJECT LiH_bulk_2x2x2_DZVP
PRINT_LEVEL MEDIUM
RUN_TYPE ENERGY
&END GLOBAL
Cubic-scaling GW calculations could be a more efficient alternative for large systems. See below an exemplary input for one water molecule. Compare the results to the ones from Sec. 1. In general, small deviations (< 0.05 eV) for GW levels can be expected from cubic-scaling GW calculations compared to canonical GW calculations due to additional approximations in cubic-scaling GW.
Please observe that the input below is much slower than the input for canonical GW. Therefore, it can be beneficial to run it with more MPI tasks. The beneficial scaling of cubic-scaling GW only pays off for large systems where it is more efficient as canonical GW calculations (rule of thumb: cubic-scaling GW can be more efficient for systems with more than 100 atoms if the filter parameters are well set).
&FORCE_EVAL
METHOD Quickstep
&DFT
! retrieve basis set from the CP2K trunk version
BASIS_SET_FILE_NAME BASIS_def2_QZVP_RI_ALL
POTENTIAL_FILE_NAME POTENTIAL
&MGRID
CUTOFF 400
REL_CUTOFF 50
&END MGRID
&QS
! all electron calculation since GW100 is all-electron test
METHOD GAPW
&END QS
&POISSON
PERIODIC NONE
PSOLVER MT
&END
&SCF
EPS_SCF 1.0E-6
SCF_GUESS ATOMIC
MAX_SCF 200
&END SCF
&XC
&XC_FUNCTIONAL PBE
&END XC_FUNCTIONAL
&WF_CORRELATION
METHOD RI_RPA_GPW
ERI_METHOD OS
! cubic-scaling GW only works with overlap metric RI
RI OVERLAP
&WFC_GPW
! EPS_FILTER should be tuned for the specific application:
! the computational cost strongly depends on EPS_FILTER
EPS_FILTER 1.0E-12
! EPS_GRID may be tuned since memory is weakly
! dependent on it
EPS_GRID 1.0E-12
&END WFC_GPW
&RI_RPA
! cubic-scaling GW only works with the minimax grid
! in imag. time and frequency
MINIMAX
! If the HOMO-LUMO gap of the system is small, 20
! points for the time/frequency grid should be used
! (flag RPA_NUM_QUAD_POINTS). The time and frequency grid
! are equally large. The maximum grid size is 20.
! For large-gap systems (as the water molecule), 12 points
! should be sufficient
RPA_NUM_QUAD_POINTS 12
! imaginary time flag enables cubic-scaling RPA or
! GW calculations
IM_TIME
&IM_TIME
! EPS_FILTER_IM_TIME should be tuned for the specific
! application: the computational cost strongly
! depends on EPS_FILTER
EPS_FILTER_IM_TIME 1.0E-12
! for large systems, increase GROUP_SIZE_3C
! to prevent out of memory (OOM)
GROUP_SIZE_3C 1
! for extremely large systems, increase GROUP_SIZE_P
! to prevent OOM
! for very large systems, it is also recommended
! to use OMP threads to prevent OOM
GROUP_SIZE_P 1
! for larger systems, MEMORY_CUT must be increased
! to prevent out of memory (OOM)
MEMORY_CUT 1
GW
&END IM_TIME
&RI_G0W0
CORR_OCC 20
CORR_VIRT 20
ANALYTIC_CONTINUATION PADE
NPARAM_PADE 16
CROSSING_SEARCH NEWTON
RI_SIGMA_X
&END RI_G0W0
&END RI_RPA
&END
&END XC
&END DFT
&SUBSYS
&CELL
ABC 10.0 10.0 10.0
PERIODIC NONE
&END CELL
&COORD
O 0.0000 0.0000 0.0000
H 0.7571 0.0000 0.5861
H -0.7571 0.0000 0.5861
&END COORD
&TOPOLOGY
&CENTER_COORDINATES
&END
&END TOPOLOGY
&KIND H
BASIS_SET def2-QZVP
RI_AUX_BASIS RI-5Z
POTENTIAL ALL
&END KIND
&KIND O
BASIS_SET def2-QZVP
RI_AUX_BASIS RI-5Z
POTENTIAL ALL
&END KIND
&END SUBSYS
&END FORCE_EVAL
&GLOBAL
RUN_TYPE ENERGY
PROJECT ALL_ELEC
PRINT_LEVEL MEDIUM
&END GLOBAL