simplified GW self-consistency input

src/core_ppnl.F

@@ -233,7 +233,7 @@ SUBROUTINE build_core_ppnl(matrix_h, matrix_p, force, virial, calculate_forces,
             IF (ASSOCIATED(gth_potential)) THEN
                gpotential(ikind)%gth_potential => gth_potential
                IF (do_soc .AND. (.NOT. gth_potential%soc)) THEN
-                  CPWARN("Spin-orbit coupling selected, but GTH potential without SOC parameters provided")
+                  CPABORT("Spin-orbit coupling selected, but GTH potential without SOC parameters provided")
                END IF
             ELSE IF (ASSOCIATED(sgp_potential)) THEN
                spotential(ikind)%sgp_potential => sgp_potential

src/gw_utils.F

@@ -262,6 +262,7 @@ SUBROUTINE read_gw_input_parameters(bs_env, bs_sec)
       CALL section_vals_val_get(gw_sec, "NUM_TIME_FREQ_POINTS", i_val=bs_env%num_time_freq_points)
       CALL section_vals_val_get(gw_sec, "EPS_FILTER", r_val=bs_env%eps_filter)
       CALL section_vals_val_get(gw_sec, "REGULARIZATION_RI", r_val=bs_env%input_regularization_RI)
+      CALL section_vals_val_get(gw_sec, "CUTOFF_RADIUS_RI", r_val=bs_env%ri_metric%cutoff_radius)
       CALL section_vals_val_get(gw_sec, "MEMORY_PER_PROC", r_val=bs_env%input_memory_per_proc_GB)
       CALL section_vals_val_get(gw_sec, "APPROX_KP_EXTRAPOL", l_val=bs_env%approx_kp_extrapol)
       CALL section_vals_val_get(gw_sec, "HEDIN_SHIFT", l_val=bs_env%do_hedin_shift)
@@ -1484,8 +1485,7 @@ SUBROUTINE set_heuristic_parameters(bs_env, qs_env)
       ! resolution of the identity with the truncated Coulomb metric, cutoff radius 3 Angström
       bs_env%ri_metric%potential_type = do_potential_truncated
       bs_env%ri_metric%omega = 0.0_dp
-      ! cutoff radius = 3 Angström
-      bs_env%ri_metric%cutoff_radius = 3.0_dp/angstrom
+      ! cutoff radius is specified in the input
       bs_env%ri_metric%filename = "t_c_g.dat"
 
       bs_env%eps_eigval_mat_RI = 0.0_dp
@@ -1525,6 +1525,8 @@ SUBROUTINE set_heuristic_parameters(bs_env, qs_env)
       IF (u > 0) THEN
          WRITE (u, FMT="(T2,2A,F21.1,A)") "Cutoff radius for the truncated Coulomb ", &
             "operator in Σ^x:", bs_env%trunc_coulomb%cutoff_radius*angstrom, " Å"
+         WRITE (u, FMT="(T2,2A,F15.1,A)") "Cutoff radius for the truncated Coulomb ", &
+            "operator in RI metric:", bs_env%ri_metric%cutoff_radius*angstrom, " Å"
          WRITE (u, FMT="(T2,A,ES48.1)") "Regularization parameter of RI ", bs_env%regularization_RI
       END IF
 
@@ -1564,6 +1566,7 @@ SUBROUTINE print_header_and_input_parameters(bs_env)
          WRITE (u, '(T2,A,I45)') 'Input: Number of time/freq. points', bs_env%num_time_freq_points
          WRITE (u, '(T2,A,ES27.1)') 'Input: Filter threshold for sparse tensor operations', &
             bs_env%eps_filter
+         WRITE (bs_env%unit_nr, FMT="(T2,A,L62)") "Apply Hedin shift", bs_env%do_hedin_shift
       END IF
 
       CALL timestop(handle)
@@ -2880,11 +2883,6 @@ SUBROUTINE analyt_conti_and_print(bs_env, Sigma_c_ikp_n_freq, Sigma_x_ikp_n, V_x
       is_bandstruc_kpoint = (ikp > bs_env%nkp_only_DOS)
       print_ikp = print_DOS_kpoints .OR. is_bandstruc_kpoint
 
-      IF (bs_env%unit_nr > 0) THEN
-         WRITE (bs_env%unit_nr, FMT="(T2,A)") " "
-         WRITE (bs_env%unit_nr, FMT="(T2,A,L80)") "Apply Hedin shift", bs_env%do_hedin_shift
-      END IF
-
       IF (bs_env%para_env%is_source() .AND. print_ikp) THEN
 
          IF (print_DOS_kpoints) THEN

src/input_constants.F

@@ -1088,7 +1088,10 @@ MODULE input_constants
                                                int_ldos_z = 21, &
                                                int_ldos_none = 22, &
                                                small_cell_full_kp = 31, &
-                                               large_cell_Gamma = 32
+                                               large_cell_Gamma = 32, &
+                                               G0W0 = 35, &
+                                               evGW0 = 36, &
+                                               evGW = 37
 
    ! periodic RESP parameters
    INTEGER, PARAMETER, PUBLIC               :: do_resp_x_dir = 0, &

src/input_cp2k_mp2.F

@@ -37,7 +37,7 @@ MODULE input_cp2k_mp2
       ot_precond_s_inverse, ri_default, ri_rpa_g0w0_crossing_bisection, &
       ri_rpa_g0w0_crossing_newton, ri_rpa_g0w0_crossing_z_shot, soc_lda, soc_none, soc_pbe, &
       wfc_mm_style_gemm, wfc_mm_style_syrk, z_solver_cg, z_solver_pople, z_solver_richardson, &
-      z_solver_sd, rpa_exchange_none, rpa_exchange_axk, rpa_exchange_sosex
+      z_solver_sd, rpa_exchange_none, rpa_exchange_axk, rpa_exchange_sosex, G0W0, evGW0, evGW
    USE input_cp2k_hfx, ONLY: create_hfx_section
    USE input_cp2k_kpoints, ONLY: create_kpoint_set_section
    USE input_keyword_types, ONLY: keyword_create, &
@@ -694,9 +694,25 @@ SUBROUTINE create_ri_g0w0(section)
       CALL section_add_keyword(section, keyword)
       CALL keyword_release(keyword)
 
+      CALL keyword_create(keyword, __LOCATION__, name="SELF_CONSISTENCY", &
+                          description="Decide the level of self-consistency of eigenvalues "// &
+                          "(= quasiparticle energies = single-electron energies) in GW. "// &
+                          "Updates of Kohn-Sham orbitals (for example qsGW) are not implemented. "// &
+                          "For details which type of eigenvalue self-consistency might be good, "// &
+                          "please consult Golze, Dvorak, Rinke, Front. Chem. 2019.", &
+                          usage="GW_SELF_CONSISTENCY evGW0", &
+                          enum_c_vals=s2a("G0W0", "evGW0", "evGW"), &
+                          enum_i_vals=(/G0W0, evGW0, evGW/), &
+                          enum_desc=s2a("Use DFT eigenvalues; not update.", &
+                                        "Update DFT eigenvalues in G, not in W.", &
+                                        "Update DFT eigenvalues in G and W."), &
+                          default_i_val=G0W0)
+      CALL section_add_keyword(section, keyword)
+      CALL keyword_release(keyword)
+
       CALL keyword_create(keyword, __LOCATION__, name="CORR_MOS_OCC", &
                           variants=(/"CORR_OCC"/), &
-                          description="Number of occupied MOs whose energies are corrected by RI-G0W0. "// &
+                          description="Number of occupied MOs whose energies are corrected in GW. "// &
                           "Counting beginning from HOMO, e.g. 3 corrected occ. MOs correspond "// &
                           "to correction of HOMO, HOMO-1 and HOMO-2. Numerical effort and "// &
                           "storage of RI-G0W0 increase linearly with this number. In case you "// &
@@ -709,7 +725,7 @@ SUBROUTINE create_ri_g0w0(section)
 
       CALL keyword_create(keyword, __LOCATION__, name="CORR_MOS_VIRT", &
                           variants=(/"CORR_VIRT"/), &
-                          description="Number of virtual MOs whose energies are corrected by RI-G0W0. "// &
+                          description="Number of virtual MOs whose energies are corrected by GW. "// &
                           "Counting beginning from LUMO, e.g. 3 corrected occ. MOs correspond "// &
                           "to correction of LUMO, LUMO+1 and LUMO+2. Numerical effort and "// &
                           "storage of RI-G0W0 increase linearly with this number. In case you "// &
@@ -773,16 +789,21 @@ SUBROUTINE create_ri_g0w0(section)
       CALL keyword_release(keyword)
 
       CALL keyword_create(keyword, __LOCATION__, name="EV_GW_ITER", &
-                          description="Maximum number of iterations for eigenvalue self-consistency cycle. The "// &
-                          "computational effort of GW scales linearly with this number. ", &
+                          description="Maximum number of iterations for eigenvalue "// &
+                          "self-consistency cycle. The computational effort of GW scales "// &
+                          "linearly with this number. In case of choosing "// &
+                          "GW_SELF_CONSISTENCY EVGW, the code sets EV_GW_ITER 10.", &
                           usage="EV_GW_ITER 3", &
                           default_i_val=1)
       CALL section_add_keyword(section, keyword)
       CALL keyword_release(keyword)
 
       CALL keyword_create(keyword, __LOCATION__, name="SC_GW0_ITER", &
-                          description="Maximum number of iterations for GW0 self-consistency cycle. The "// &
-                          "computational effort of GW is not much affected by the number of scGW0 cycles. ", &
+                          description="Maximum number of iterations for GW0 "// &
+                          "self-consistency cycle. The computational effort "// &
+                          "of GW is not much affected by the number of scGW0 cycles. "// &
+                          "In case of choosing "// &
+                          "GW_SELF_CONSISTENCY EVGW0, the code sets SC_GW0_ITER 10.", &
                           usage="SC_GW0_ITER 3", &
                           default_i_val=1)
       CALL section_add_keyword(section, keyword)

src/input_cp2k_properties_dft.F

@@ -2073,6 +2073,18 @@ SUBROUTINE create_gw_section(section)
       CALL section_add_keyword(section, keyword)
       CALL keyword_release(keyword)
 
+      CALL keyword_create(keyword, __LOCATION__, name="CUTOFF_RADIUS_RI", &
+                          description="The cutoff radius (in Angstrom) for the truncated  "// &
+                          "Coulomb operator. The larger the cutoff radius, the faster "// &
+                          "converges the resolution of the identity (RI) with respect to the "// &
+                          "RI basis set size. Larger cutoff radius means higher computational "// &
+                          "cost.", &
+                          usage="CUTOFF_RADIUS_RI 3.0", &
+                          default_r_val=cp_unit_to_cp2k(value=3.0_dp, unit_str="angstrom"), &
+                          type_of_var=real_t, unit_str="angstrom")
+      CALL section_add_keyword(section, keyword)
+      CALL keyword_release(keyword)
+
       CALL keyword_create(keyword, __LOCATION__, name="MEMORY_PER_PROC", &
                           description="Specify the available memory per MPI process. Set "// &
                           "`MEMORY_PER_PROC` as accurately as possible for good performance. If "// &
@@ -2112,7 +2124,7 @@ SUBROUTINE create_gw_section(section)
                           "similar GW eigenvalues as evGW0; at a lower "// &
                           "computational cost.", &
                           usage="HEDIN_SHIFT", &
-                          default_l_val=.TRUE., &
+                          default_l_val=.FALSE., &
                           lone_keyword_l_val=.TRUE.)
       CALL section_add_keyword(section, keyword)
       CALL keyword_release(keyword)

src/mp2_setup.F

@@ -21,8 +21,9 @@ MODULE mp2_setup
                                               cp_print_key_unit_nr
    USE cp_parser_methods,               ONLY: read_float_object
    USE input_constants,                 ONLY: &
-        do_eri_mme, do_potential_short, mp2_method_direct, mp2_method_gpw, mp2_method_none, &
-        mp2_ri_optimize_basis, ri_mp2_laplace, ri_mp2_method_gpw, ri_rpa_method_gpw
+        do_eri_mme, do_potential_short, evGW, evGW0, mp2_method_direct, mp2_method_gpw, &
+        mp2_method_none, mp2_ri_optimize_basis, ri_mp2_laplace, ri_mp2_method_gpw, &
+        ri_rpa_method_gpw
    USE input_section_types,             ONLY: section_vals_get_subs_vals,&
                                               section_vals_type,&
                                               section_vals_val_get
@@ -65,7 +66,8 @@ SUBROUTINE read_mp2_section(input, mp2_env)
       CHARACTER(LEN=default_string_length), &
          DIMENSION(:), POINTER                           :: string_pointer
       CHARACTER(LEN=max_line_length)                     :: error_message
-      INTEGER                                            :: handle, i, i_special_kp, ival, unit_nr
+      INTEGER                                            :: gw_sc, handle, i, i_special_kp, ival, &
+                                                            unit_nr
       INTEGER, DIMENSION(:), POINTER                     :: tmplist
       LOGICAL                                            :: do_mp2, do_opt_ri_basis, do_ri_mp2, &
                                                             do_ri_sos_mp2, do_rpa
@@ -124,6 +126,8 @@ SUBROUTINE read_mp2_section(input, mp2_env)
       gw_section => section_vals_get_subs_vals(mp2_section, "RI_RPA%GW")
       CALL section_vals_val_get(gw_section, "_SECTION_PARAMETERS_", &
                                 l_val=mp2_env%ri_rpa%do_ri_g0w0)
+      CALL section_vals_val_get(mp2_section, "RI_RPA%GW%SELF_CONSISTENCY", &
+                                i_val=gw_sc)
       CALL section_vals_val_get(mp2_section, "RI_RPA%GW%CORR_MOS_OCC", &
                                 i_val=mp2_env%ri_g0w0%corr_mos_occ)
       CALL section_vals_val_get(mp2_section, "RI_RPA%GW%CORR_MOS_VIRT", &
@@ -144,8 +148,10 @@ SUBROUTINE read_mp2_section(input, mp2_env)
                                 l_val=mp2_env%ri_g0w0%do_hedin_shift)
       CALL section_vals_val_get(mp2_section, "RI_RPA%GW%EV_GW_ITER", &
                                 i_val=mp2_env%ri_g0w0%iter_evGW)
+      IF (gw_sc == evGW) mp2_env%ri_g0w0%iter_evGW = 10
       CALL section_vals_val_get(mp2_section, "RI_RPA%GW%SC_GW0_ITER", &
                                 i_val=mp2_env%ri_g0w0%iter_sc_GW0)
+      IF (gw_sc == evGW0) mp2_env%ri_g0w0%iter_sc_GW0 = 10
       CALL section_vals_val_get(mp2_section, "RI_RPA%GW%EPS_ITER", &
                                 r_val=mp2_env%ri_g0w0%eps_iter)
       CALL section_vals_val_get(mp2_section, "RI_RPA%GW%PRINT_EXX", &
@@ -154,6 +160,7 @@ SUBROUTINE read_mp2_section(input, mp2_env)
                                 l_val=mp2_env%ri_g0w0%print_self_energy)
       CALL section_vals_val_get(mp2_section, "RI_RPA%GW%RI_SIGMA_X", &
                                 l_val=mp2_env%ri_g0w0%do_ri_Sigma_x)
+
       NULLIFY (r_vals)
       CALL section_vals_val_get(mp2_section, "RI_RPA%GW%IC_CORR_LIST", &
                                 r_vals=r_vals)

src/rpa_gw.F

@@ -5095,15 +5095,25 @@ SUBROUTINE print_and_update_for_ev_sc(vec_gw_energ, &
                WRITE (unit_nr, '(T3,A,F42.4)') 'GW direct gap at current kpoint (eV)', E_GAP_GW*evolt
             ELSE IF (do_alpha) THEN
                WRITE (unit_nr, '(T3,A)') ' '
-               WRITE (unit_nr, '(T3,A,F36.4)') 'Alpha GW direct gap at current kpoint (eV)', E_GAP_GW*evolt
+               WRITE (unit_nr, '(T3,A,F36.4)') 'Alpha GW direct gap at current kpoint (eV)', &
+                  E_GAP_GW*evolt
             ELSE IF (do_beta) THEN
                WRITE (unit_nr, '(T3,A)') ' '
-               WRITE (unit_nr, '(T3,A,F37.4)') 'Beta GW direct gap at current kpoint (eV)', E_GAP_GW*evolt
+               WRITE (unit_nr, '(T3,A,F37.4)') 'Beta GW direct gap at current kpoint (eV)', &
+                  E_GAP_GW*evolt
             END IF
          ELSE
             IF (do_closed_shell) THEN
                WRITE (unit_nr, '(T3,A)') ' '
-               WRITE (unit_nr, '(T3,A,F57.4)') 'GW HOMO-LUMO gap (eV)', E_GAP_GW*evolt
+               IF (count_ev_sc_GW > 1) THEN
+                  WRITE (unit_nr, '(T3,A,I3,A,F39.4)') 'HOMO-LUMO gap in evGW iteration', &
+                     count_ev_sc_GW, ' (eV)', E_GAP_GW*evolt
+               ELSE IF (count_sc_GW0 > 1) THEN
+                  WRITE (unit_nr, '(T3,A,I3,A,F38.4)') 'HOMO-LUMO gap in evGW0 iteration', &
+                     count_sc_GW0, ' (eV)', E_GAP_GW*evolt
+               ELSE
+                  WRITE (unit_nr, '(T3,A,F55.4)') 'G0W0 HOMO-LUMO gap (eV)', E_GAP_GW*evolt
+               END IF
             ELSE IF (do_alpha) THEN
                WRITE (unit_nr, '(T3,A)') ' '
                WRITE (unit_nr, '(T3,A,F51.4)') 'Alpha GW HOMO-LUMO gap (eV)', E_GAP_GW*evolt

tests/QS/regtest-gw-cubic/TEST_FILES

@@ -1,5 +1,5 @@
-scGW0_H2O_PBE_default_values.inp                      78      1e-04                         17.188
-evGW_H2O_PBE_default_values.inp                       78      1e-04                         17.474
+scGW0_H2O_PBE_default_values.inp                     111      1e-04                        17.2026
+evGW_H2O_PBE_default_values.inp                      112      1e-04                        17.5131
 G0W0_H2O_PBE_GAPW.inp                                 78      1e-04                         16.699
 G0W0_H2O_PBE0.inp                                     78      1e-04                         16.717
 G0W0_H2O_PBE0_30_pts.inp                              78      1e-04                         15.607

tests/QS/regtest-gw-cubic/evGW_H2O_PBE_default_values.inp

@@ -45,7 +45,7 @@
           &GW
             CORR_MOS_OCC 1
             CORR_MOS_VIRT 1
-            EV_GW_ITER 2
+            EV_GW_ITER 3
             RI_SIGMA_X FALSE
           &END GW
           &HF

tests/QS/regtest-gw-cubic/scGW0_H2O_PBE_default_values.inp

@@ -46,7 +46,7 @@
             CORR_MOS_OCC 1
             CORR_MOS_VIRT 1
             RI_SIGMA_X FALSE
-            SC_GW0_ITER 2
+            SC_GW0_ITER 3
           &END GW
           &HF
             FRACTION 1.0000000

tests/QS/regtest-gw/TEST_FILES

@@ -11,4 +11,6 @@ evGW_OH_PBE_svd.inp                                   11      1e-08            -
 evGW_H2O_PBE_RI_HFX_svd.inp                           11      1e-08            -13.112474185078236
 scGW0_H2O_PBE.inp                                     11      1e-08            -17.118441200805574
 G0W0_H2O_PBE0_Hedin_shift.inp                         78      1e-04                         16.806
+evGW0_H2O_PBE_set_sc_flag.inp                        111      1e-04                        16.5546
+evGW_H2O_PBE_set_sc_flag.inp                         112      1e-04                        17.0863
 #EOF

tests/QS/regtest-gw/evGW0_H2O_PBE_set_sc_flag.inp

@@ -0,0 +1,93 @@
+&GLOBAL
+  PRINT_LEVEL MEDIUM
+  PROJECT G0W0_H2O_PBE_ev_sc
+  RUN_TYPE ENERGY
+  &TIMINGS
+    THRESHOLD 0.01
+  &END TIMINGS
+&END GLOBAL
+
+&FORCE_EVAL
+  METHOD Quickstep
+  &DFT
+    BASIS_SET_FILE_NAME HFX_BASIS
+    POTENTIAL_FILE_NAME GTH_POTENTIALS
+    &MGRID
+      CUTOFF 100
+      REL_CUTOFF 20
+    &END MGRID
+    &POISSON
+      PERIODIC NONE
+      POISSON_SOLVER WAVELET
+    &END POISSON
+    &QS
+      EPS_DEFAULT 1.0E-15
+      EPS_PGF_ORB 1.0E-30
+      METHOD GPW
+    &END QS
+    &SCF
+      EPS_SCF 1.0E-7
+      MAX_SCF 100
+      SCF_GUESS ATOMIC
+      &PRINT
+        &RESTART OFF
+        &END RESTART
+      &END PRINT
+    &END SCF
+    &XC
+      &WF_CORRELATION
+        MEMORY 200.
+        NUMBER_PROC 1
+        &INTEGRALS
+          &WFC_GPW
+            CUTOFF 100
+            REL_CUTOFF 20
+          &END WFC_GPW
+        &END INTEGRALS
+        &RI_RPA
+          RPA_NUM_QUAD_POINTS 10
+          &GW
+            CORR_MOS_OCC 10
+            CORR_MOS_VIRT 10
+            SELF_CONSISTENCY EVGW0
+          &END GW
+          &HF
+            FRACTION 1.0000000
+            &SCREENING
+              EPS_SCHWARZ 1.0E-6
+              SCREEN_ON_INITIAL_P FALSE
+            &END SCREENING
+          &END HF
+        &END RI_RPA
+      &END WF_CORRELATION
+      &XC_FUNCTIONAL PBE
+        &PBE
+          SCALE_C 1.0000000
+          SCALE_X 1.0000000
+        &END PBE
+      &END XC_FUNCTIONAL
+    &END XC
+  &END DFT
+  &SUBSYS
+    &CELL
+      ABC [angstrom] 6.000 6.000 6.000
+      PERIODIC NONE
+    &END CELL
+    &KIND H
+      BASIS_SET DZVP-GTH
+      BASIS_SET RI_AUX RI_DZVP-GTH
+      POTENTIAL GTH-PBE-q1
+    &END KIND
+    &KIND O
+      BASIS_SET DZVP-GTH
+      BASIS_SET RI_AUX RI_DZVP-GTH
+      POTENTIAL GTH-PBE-q6
+    &END KIND
+    &TOPOLOGY
+      COORD_FILE_FORMAT xyz
+      COORD_FILE_NAME H2O_gas.xyz
+      &CENTER_COORDINATES
+      &END CENTER_COORDINATES
+    &END TOPOLOGY
+  &END SUBSYS
+&END FORCE_EVAL