howto:gw
Differences
This shows you the differences between two versions of the page.
Next revision | Previous revisionNext revisionBoth sides next revision | ||
howto:gw [2021/04/21 19:51] – created jwilhelm | howto:gw [2023/10/18 18:00] – [5. GW for 2D materials: Example monolayer MoS2] jwilhelm | ||
---|---|---|---|
Line 1: | Line 1: | ||
====== GW method for computing electronic levels ====== | ====== GW method for computing electronic levels ====== | ||
- | The purpose of this section is to explain how to compute the energy of a molecular | + | The purpose of this section is to explain how to compute the energy of molecular |
- | The GW implementation in CP2K is based on the developments described in [[doi> | + | The GW implementation in CP2K is based on the developments described in [[doi> |
- | + | ||
- | In this tutorial, GW values from the GW100 benchmark set [[doi> | + | |
Since the calculations are rather small, please use a single MPI rank for the calculation: | Since the calculations are rather small, please use a single MPI rank for the calculation: | ||
Line 14: | Line 12: | ||
===== 1. Reproducing values from the GW100 set ===== | ===== 1. Reproducing values from the GW100 set ===== | ||
- | 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 " | + | 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 " |
- | 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. [[doi> | + | For checking the basis set convergence, |
<code - H2O_GW100.inp> | <code - H2O_GW100.inp> | ||
& | & | ||
Line 45: | Line 43: | ||
! GW is part of the WF_CORRELATION section | ! GW is part of the WF_CORRELATION section | ||
& | & | ||
- | ! 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 | &RI_RPA | ||
- | ! use 100 quadrature | + | ! use 100 points to perform the frequency integration in GW |
- | ! frequency integration in GW | + | |
- | | + | ! SIZE_FREQ_INTEG_GROUP is a group size for parallelization and |
- | ! SIZE_FREQ_INTEG_GROUP is a group size for parallelization and | + | |
! should be increased for large calculations to prevent out of memory. | ! should be increased for large calculations to prevent out of memory. | ||
! maximum for SIZE_FREQ_INTEG_GROUP is the number of MPI tasks | ! maximum for SIZE_FREQ_INTEG_GROUP is the number of MPI tasks | ||
- | | + | &GW |
- | GW | + | ! compute the G0W0@PBE energy of HOMO-9, |
- | | + | ! HOMO-8, ... , HOMO-1, HOMO |
- | ! compute the G0W0@PBE energy of HOMO-9, | + | |
- | ! HOMO-8, ... , HOMO-1, HOMO | + | |
| | ||
- | ! compute the G0W0@PBE energy of LUMO, | + | ! compute the G0W0@PBE energy of LUMO, |
! LUMO+1, ... , LUMO+20 | ! LUMO+1, ... , LUMO+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 | ||
- | | ||
- | ! for solving the quasiparticle equation, the Newton method | ||
- | ! is used as in the GW100 benchmark | ||
- | | ||
! use the RI approximation for the exchange part of the self-energy | ! use the RI approximation for the exchange part of the self-energy | ||
| | ||
- | & | + | & |
&END RI_RPA | &END RI_RPA | ||
- | ! NUMBER_PROC is a group size for parallelization and should | ||
- | ! be increased for large calculations | ||
- | NUMBER_PROC 1 | ||
&END | &END | ||
&END XC | &END XC | ||
Line 99: | Line 79: | ||
&KIND H | &KIND H | ||
! def2-QZVP is the basis which has been used in the GW100 paper | ! def2-QZVP is the basis which has been used in the GW100 paper | ||
- | BASIS_SET def2-QZVP | + | BASIS_SET |
- | ! just use a very large RI basis to ensure excellent | + | ! just use a very large RI basis to ensure excellent |
! convergence with respect to the RI basis | ! convergence with respect to the RI basis | ||
- | | + | |
POTENTIAL ALL | POTENTIAL ALL | ||
&END KIND | &END KIND | ||
&KIND O | &KIND O | ||
- | BASIS_SET def2-QZVP | + | BASIS_SET |
- | | + | |
POTENTIAL ALL | POTENTIAL ALL | ||
&END KIND | &END KIND | ||
Line 117: | Line 97: | ||
PRINT_LEVEL | PRINT_LEVEL | ||
&END GLOBAL | &END GLOBAL | ||
+ | |||
</ | </ | ||
Line 130: | Line 111: | ||
< | < | ||
&KIND H | &KIND H | ||
- | BASIS_SET cc-DZVP-all | + | BASIS_SET |
- | | + | |
POTENTIAL ALL | POTENTIAL ALL | ||
&END KIND | &END KIND | ||
&KIND O | &KIND O | ||
- | BASIS_SET cc-DZVP-all | + | BASIS_SET |
- | | + | |
POTENTIAL ALL | POTENTIAL ALL | ||
&END KIND | &END KIND | ||
Line 159: | Line 140: | ||
</ | </ | ||
- | For the extrapolation, | + | For the extrapolation, |
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. | 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. | 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 [[doi> | + | The basis set extrapolated values from the table above deviate from the values reported in the GW100 paper [[doi> |
- | 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, | + | 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, |
===== 3. Input for large-scale calculations ===== | ===== 3. Input for large-scale calculations ===== | ||
An exemplary input for a parallel calculation can be found in the supporting information of [[doi> | An exemplary input for a parallel calculation can be found in the supporting information of [[doi> | ||
- | ===== 4. Periodic | + | ===== 4. Self-consistent |
- | For periodic GW calculations, | + | |
- | The basis can be found in {{exercises: | + | The G0W0@PBE HOMO value of the H2O molecule (~ -12.0 eV) is not in good agreement with the experimental ionization potential (12.62 eV). Benchmarks on molecules and solids indicate that self-consistency |
- | < | + | You can run GW0 calculations in CP2K by putting |
+ | <code> | ||
+ | &GW | ||
+ | SC_GW0_ITER | ||
+ | CORR_OCC | ||
+ | CORR_VIRT | ||
+ | RI_SIGMA_X | ||
+ | &END GW | ||
+ | </ | ||
+ | " | ||
+ | |||
+ | ===== 5. GW for 2D materials: Example monolayer MoS2 ===== | ||
+ | There is also a periodic GW implementation [[doi> | ||
+ | |||
+ | |||
+ | For computing the G0W0@LDA quasiparticle energy levels of monolayer MoS2, please use the input file | ||
+ | <code> | ||
+ | & | ||
+ | PROJECT | ||
+ | RUN_TYPE ENERGY | ||
+ | &END GLOBAL | ||
& | & | ||
METHOD Quickstep | METHOD Quickstep | ||
&DFT | &DFT | ||
- | BASIS_SET_FILE_NAME | + | BASIS_SET_FILE_NAME |
- | POTENTIAL_FILE_NAME | + | POTENTIAL_FILE_NAME |
+ | SORT_BASIS EXP | ||
&MGRID | &MGRID | ||
- | CUTOFF | + | CUTOFF |
- | REL_CUTOFF | + | REL_CUTOFF |
&END MGRID | &END MGRID | ||
&QS | &QS | ||
METHOD GPW | METHOD GPW | ||
- | EPS_DEFAULT 1.0E-15 | + | EPS_DEFAULT 1.0E-12 |
- | EPS_PGF_ORB 1.0E-20 | + | EPS_PGF_ORB 1.0E-12 |
- | EPS_FILTER_MATRIX 0.0e0 | + | |
&END QS | &END QS | ||
&SCF | &SCF | ||
- | EPS_SCF 1.0E-6 | + | |
+ | | ||
MAX_SCF 100 | MAX_SCF 100 | ||
+ | &MIXING | ||
+ | METHOD BROYDEN_MIXING | ||
+ | ALPHA 0.1 | ||
+ | BETA 1.5 | ||
+ | NBROYDEN 8 | ||
+ | &END | ||
&END SCF | &END SCF | ||
&XC | &XC | ||
- | & | + | & |
&END XC_FUNCTIONAL | &END XC_FUNCTIONAL | ||
- | & | ||
- | METHOD | ||
- | &RI_RPA | ||
- | RPA_NUM_QUAD_POINTS | ||
- | GW | ||
- | & | ||
- | | ||
- | | ||
- | ! activate the periodic correction | ||
- | | ||
- | | ||
- | | ||
- | &END RI_G0W0 | ||
- | ! HF calculation for the exchange part of the self-energy | ||
- | ! Here, the truncation of the Coulomb operator works | ||
- | &HF | ||
- | & | ||
- | ! for other materials, a smaller EPS_SCHWARZ might be necessary | ||
- | EPS_SCHWARZ 1.0E-6 | ||
- | SCREEN_ON_INITIAL_P TRUE | ||
- | &END | ||
- | & | ||
- | POTENTIAL_TYPE TRUNCATED | ||
- | ! the truncation radius is half the cell size | ||
- | CUTOFF_RADIUS | ||
- | T_C_G_DATA t_c_g.dat | ||
- | &END | ||
- | &MEMORY | ||
- | MAX_MEMORY | ||
- | &END | ||
- | &END | ||
- | &END RI_RPA | ||
- | NUMBER_PROC | ||
- | &END | ||
&END XC | &END XC | ||
&END DFT | &END DFT | ||
+ | & | ||
+ | & | ||
+ | &DOS | ||
+ | ! k-point mesh for the self-energy | ||
+ | KPOINTS 2 2 1 | ||
+ | &END | ||
+ | &GW | ||
+ | ! for details on parameters, please consult | ||
+ | ! manual.cp2k.org/ | ||
+ | NUM_TIME_FREQ_POINTS | ||
+ | MEMORY_PER_PROC | ||
+ | EPS_FILTER | ||
+ | &END | ||
+ | &SOC | ||
+ | &END | ||
+ | &END | ||
+ | &END PROPERTIES | ||
&SUBSYS | &SUBSYS | ||
&CELL | &CELL | ||
- | ABC 4.084 4.084 4.084 | + | ABC 3.15 3.15 15.0 |
+ | ALPHA_BETA_GAMMA | ||
+ | PERIODIC XY | ||
+ | ! the calculation is on a 9x9 supercell with 243 atoms | ||
+ | MULTIPLE_UNIT_CELL 9 9 1 | ||
&END CELL | &END CELL | ||
- | &COORD | + | &TOPOLOGY |
- | | + | |
- | Li 2.042 2.042 0 | + | & |
- | Li 2.042 0 2.042 | + | |
- | Li 0 2.042 2.042 | + | & |
- | H 0 2.042 0 | + | BASIS_SET |
- | H 0 0 2.042 | + | |
- | H 2.042 0 0 | + | POTENTIAL |
- | H 2.042 2.042 2.042 | + | |
- | & | + | |
- | & | + | |
- | BASIS_SET | + | |
- | | + | |
- | POTENTIAL GTH-PBE-q1 | + | |
&END KIND | &END KIND | ||
- | | + | |
- | BASIS_SET | + | |
- | | + | BASIS_SET |
- | POTENTIAL GTH-PBE-q3 | + | |
+ | POTENTIAL | ||
&END KIND | &END KIND | ||
- | &END SUBSYS | ||
- | &END FORCE_EVAL | ||
- | &GLOBAL | ||
- | PROJECT | ||
- | PRINT_LEVEL MEDIUM | ||
- | RUN_TYPE ENERGY | ||
- | &END GLOBAL | ||
- | </ | ||
+ | &KIND Mo | ||
+ | BASIS_SET ORB TZVP-MOLOPT-GTH_upscaled | ||
+ | BASIS_SET RI_AUX RI | ||
+ | POTENTIAL | ||
+ | &END KIND | ||
- | ===== 5. Cubic-scaling GW calculations ===== | + | &KIND W |
- | 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. | + | |
+ | BASIS_SET RI_AUX RI | ||
+ | POTENTIAL | ||
+ | &END KIND | ||
- | 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). | ||
- | |||
- | <code - H2O_GW100_cubic_scaling.inp> | ||
- | & | ||
- | 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 | ||
- | & | ||
- | PERIODIC NONE | ||
- | PSOLVER MT | ||
- | &END | ||
- | &SCF | ||
- | EPS_SCF 1.0E-6 | ||
- | SCF_GUESS ATOMIC | ||
- | MAX_SCF 200 | ||
- | &END SCF | ||
- | &XC | ||
- | & | ||
- | &END XC_FUNCTIONAL | ||
- | & | ||
- | METHOD RI_RPA_GPW | ||
- | ERI_METHOD OS | ||
- | ! cubic-scaling GW only works with overlap metric RI | ||
- | RI OVERLAP | ||
- | & | ||
- | ! 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 | ||
- | &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/ | ||
- | ! (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 | ||
- | ! imaginary time flag enables cubic-scaling RPA or | ||
- | ! GW calculations | ||
- | IM_TIME | ||
- | & | ||
- | ! EPS_FILTER_IM_TIME should be tuned for the specific | ||
- | ! application: | ||
- | ! 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 | ||
- | GW | ||
- | &END IM_TIME | ||
- | & | ||
- | | ||
- | | ||
- | | ||
- | | ||
- | | ||
- | | ||
- | &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 | &COORD | ||
- | O | + | Mo 0.00000 1.81865 3.07500 |
- | | + | |
- | | + | |
&END COORD | &END COORD | ||
- | & | ||
- | & | ||
- | &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 SUBSYS | ||
&END FORCE_EVAL | &END FORCE_EVAL | ||
- | &GLOBAL | ||
- | RUN_TYPE | ||
- | PROJECT | ||
- | PRINT_LEVEL | ||
- | &END GLOBAL | ||
</ | </ | ||
+ | Running the input file requires access to a large computer (the calculation took 2.5 hours on 32 nodes on Noctua2 cluster in Paderborn). You find the input and output files here: | ||
+ | |||
+ | https:// | ||
+ | |||
+ | The quasiparticle levels are printed to the files SCF_and_G0W0_band_structure_for_kpoint_xyz. | ||
+ | |||
+ | Some remarks: | ||
+ | |||
+ | * You can find the G0W0 bandgap in the cp2k output file in the line | ||
+ | < | ||
+ | G0W0 indirect band gap (eV): 2.470 | ||
+ | </ | ||
+ | * For adjusting the keywords NUM_TIME_FREQ_POINTS, | ||
+ | |||
+ | * The computational parameters from the input file reach numerical convergence of the band gap within ~ 50 meV (TZVP basis set, 10 time/ | ||
+ | |||
+ | * The code also outputs SOC splittings of the levels based on the SOC parameters from Hartwigsen-Goedecker-Hutter pseudopotentials [[doi> | ||
+ | |||
+ | * The code prints restart files with ending .matrix that can be used to restart a crashed calculation. | ||
+ | |||
+ | In case anything does not work, please feel free to contact jan.wilhelm (at) ur.de. |
howto/gw.txt · Last modified: 2024/01/14 12:15 by oschuett