GW2X| Implementation of the GW2X method for core ionization potential…

…s and XAS_TDP correction

src/common/bibliography.F

@@ -91,7 +91,7 @@ MODULE bibliography
                             Toulouse2004, Limpanuparb2011, Martin2003, Yin2017, Goerigk2017, &
                             Wilhelm2016a, Wilhelm2016b, Wilhelm2017, Wilhelm2018, Lass2018, cp2kqs2020, &
                             Behler2007, Behler2011, Schran2020a, Schran2020b, &
-                            Rycroft2009, Thomas2015, Brehm2018, Brehm2020
+                            Rycroft2009, Thomas2015, Brehm2018, Brehm2020, Shigeta2001
 
 CONTAINS
 
@@ -4675,6 +4675,19 @@ SUBROUTINE add_all_references()
                          "ER"), &
                          DOI="10.1063/5.0005078")
 
+      CALL add_reference(key=Shigeta2001, ISI_record=s2a( &
+                         "AU Y. Shigeta", &
+                         "   A. M. Ferreira", &
+                         "   V. G. Zakrzewski", &
+                         "   J. V. Ortiz", &
+                         "TI Electron propagator calculations with Kohn–Sham reference states", &
+                         "SO International Journal of Quantum Chemistry", &
+                         "PY 2001", &
+                         "VL 85", &
+                         "IS 4-5", &
+                         "ER"), &
+                         DOI="10.1002/qua.1543")
+
    END SUBROUTINE add_all_references
 
 END MODULE bibliography

src/common/mathlib.F

@@ -62,7 +62,7 @@ MODULE mathlib
              dotprod_3d, &
              transpose_3d, &
              expint, abnormal_value, &
-             get_diag
+             get_diag, set_diag
 
    INTERFACE det_3x3
       MODULE PROCEDURE det_3x3_1, det_3x3_2

src/fm/cp_fm_basic_linalg.F

@@ -15,7 +15,7 @@ MODULE cp_fm_basic_linalg
    USE cp_blacs_env,                    ONLY: cp_blacs_env_type
    USE cp_fm_struct,                    ONLY: cp_fm_struct_equivalent
    USE cp_fm_types,                     ONLY: &
-        cp_fm_create, cp_fm_get_element, cp_fm_get_info, cp_fm_get_submatrix, cp_fm_p_type, &
+        cp_fm_create, cp_fm_get_diag, cp_fm_get_info, cp_fm_get_submatrix, cp_fm_p_type, &
         cp_fm_release, cp_fm_set_all, cp_fm_set_element, cp_fm_set_submatrix, cp_fm_to_fm, &
         cp_fm_type
    USE cp_log_handling,                 ONLY: cp_to_string
@@ -272,6 +272,7 @@ END SUBROUTINE cp_fm_geadd
 !>      - SERIOUS WARNING (KNOWN BUG) : the sign of the determinant depends on ipivot
 !>      - one should be able to find out if ipivot is an even or an odd permutation...
 !>        if you need the correct sign, just add correct_sign==.TRUE. (fschiff)
+!>      - Use cp_fm_get_diag instead of n times cp_fm_get_element (A. Bussy)
 ! **************************************************************************************************
    SUBROUTINE cp_fm_lu_decompose(matrix_a, almost_determinant, correct_sign)
       TYPE(cp_fm_type), POINTER                :: matrix_a
@@ -304,9 +305,7 @@ SUBROUTINE cp_fm_lu_decompose(matrix_a, almost_determinant, correct_sign)
 
       ALLOCATE (diag(n))
       diag(:) = 0.0_dp
-      DO i = 1, n
-         CALL cp_fm_get_element(matrix_a, i, i, diag(i)) !  not completely optimal in speed i would say
-      ENDDO
+      CALL cp_fm_get_diag(matrix_a, diag)
       determinant = 1.0_dp
       DO i = 1, n
          determinant = determinant*diag(i)
@@ -1452,6 +1451,7 @@ END SUBROUTINE cp_fm_row_scale
 !>      - computation of determinant corrected
 !>      - determinant only computed if det_a is present
 !>      12.2016 added option to use SVD instead of LU [Nico Holmberg]
+!>      - Use cp_fm_get diag instead of n times cp_fm_get_element (A. Bussy)
 !> \author Florian Schiffmann(02.2007)
 ! **************************************************************************************************
    SUBROUTINE cp_fm_invert(matrix_a, matrix_inverse, det_a, eps_svd, eigval)
@@ -1505,9 +1505,7 @@ SUBROUTINE cp_fm_invert(matrix_a, matrix_inverse, det_a, eps_svd, eigval)
 
             ALLOCATE (diag(n))
             diag(:) = 0.0_dp
-            DO i = 1, n
-               CALL cp_fm_get_element(matrix_lu, i, i, diag(i)) !  not completely optimal in speed i would say
-            ENDDO
+            CALL cp_fm_get_diag(matrix_lu, diag)
 
             exponent_of_minus_one = 0
             determinant = 1.0_dp

src/input_cp2k_dft.F

@@ -21,7 +21,7 @@ MODULE input_cp2k_dft
         Heinzmann1976, Holmberg2017, Holmberg2018, Hunt2003, Iannuzzi2005, Iannuzzi2006, &
         Iannuzzi2007, Kolafa2004, Krack2000, Krack2002, Kuhne2007, Kunert2003, Lippert1997, &
         Lippert1999, Lu2004, Perdew1981, Repasky2002, Rocha2006, Rycroft2009, Schenter2008, &
-        Schiffmann2015, Stewart1982, Stewart1989, Stewart2007, Thiel1992, Thomas2015, &
+        Schiffmann2015, Shigeta2001, Stewart1982, Stewart1989, Stewart2007, Thiel1992, Thomas2015, &
         VanVoorhis2015, VandeVondele2003, VandeVondele2005a, VandeVondele2005b, VandeVondele2006, &
         Weber2008, Yin2017
    USE cp_output_handling,              ONLY: add_last_numeric,&
@@ -7818,7 +7818,7 @@ SUBROUTINE create_xas_tdp_section(section)
                           description="XAS simulations using linear-response TDDFT. Excitation from "// &
                           "specified core orbitals are considered one at a time. In case of high "// &
                           "symmetry structures, donor core orbitals should be localized.", &
-                          n_keywords=15, n_subsections=4, repeats=.FALSE.)
+                          n_keywords=19, n_subsections=4, repeats=.FALSE.)
 
       NULLIFY (keyword, subsection, print_key)
 
@@ -7921,6 +7921,77 @@ SUBROUTINE create_xas_tdp_section(section)
       CALL section_add_keyword(section, keyword)
       CALL keyword_release(keyword)
 
+!  the GW2X correction subsection
+      CALL section_create(subsection, __LOCATION__, name="GW2X", &
+                          description="Specifications for the GW2X calculation of core "// &
+                          "ionization potentials. Note that donor states need to be actively "// &
+                          "localized using the LOCALIZE keyword in DONOR_STATES. N_SEARCH"// &
+                          "should be kept to the minimum, such that only core states are localized.", &
+                          citations=(/Shigeta2001/), &
+                          n_keywords=8, &
+                          n_subsections=0, &
+                          repeats=.FALSE.)
+
+      CALL keyword_create(keyword, __LOCATION__, name="_SECTION_PARAMETERS_", &
+                          description="Enables the GW2X correction of the core ionization potentials", &
+                          usage="&GW2X {logical}", &
+                          default_l_val=.FALSE., &
+                          lone_keyword_l_val=.TRUE.)
+      CALL section_add_keyword(subsection, keyword)
+      CALL keyword_release(keyword)
+
+      CALL keyword_create(keyword, __LOCATION__, name="XPS_ONLY", &
+                          description="If set to .TRUE., only run GW2X calculations for XPS "// &
+                          "spectroscopy and ignore all XAS calculations. It is still "// &
+                          "required to define the DONOR_STATES and KERNEL%EXACT_EXCHANGE "// &
+                          "subsections.", &
+                          default_l_val=.FALSE., &
+                          lone_keyword_l_val=.TRUE.)
+      CALL section_add_keyword(subsection, keyword)
+      CALL keyword_release(keyword)
+
+      CALL keyword_create(keyword, __LOCATION__, name="BATCH_SIZE", &
+                          description="MOs batch size for batched tensor contraction. Larger "// &
+                          "size is faster, but uses more memory. Default should be safe "// &
+                          "in most cases.", &
+                          default_i_val=64)
+      CALL section_add_keyword(subsection, keyword)
+      CALL keyword_release(keyword)
+
+      CALL keyword_create(keyword, __LOCATION__, name="EPS_GW2X", &
+                          description="Convergence threshold for GW2X iterations (in eV)", &
+                          default_r_val=1.E-2_dp)
+      CALL section_add_keyword(subsection, keyword)
+      CALL keyword_release(keyword)
+
+      CALL keyword_create(keyword, __LOCATION__, name="C_OS", &
+                          description="Opposite-spin scling factor. SCS => 6/5, SOS => 1.3", &
+                          default_r_val=1.0_dp)
+      CALL section_add_keyword(subsection, keyword)
+      CALL keyword_release(keyword)
+
+      CALL keyword_create(keyword, __LOCATION__, name="C_SS", &
+                          description="Same-spin scling factor. SCS => 1/3, SOS => 0.0", &
+                          default_r_val=1.0_dp)
+      CALL section_add_keyword(subsection, keyword)
+      CALL keyword_release(keyword)
+
+      CALL keyword_create(keyword, __LOCATION__, name="MAX_GW2X_ITER", &
+                          description="Maximum number of iterations for GW2X", &
+                          default_i_val=10)
+      CALL section_add_keyword(subsection, keyword)
+      CALL keyword_release(keyword)
+
+      CALL keyword_create(keyword, __LOCATION__, name="PSEUDO_CANONICAL", &
+                          description="Whether the pseudo-canonical version of GW2X should be used "// &
+                          "(versus only using the diagonal of the generalized Fock matrix)", &
+                          default_l_val=.TRUE.)
+      CALL section_add_keyword(subsection, keyword)
+      CALL keyword_release(keyword)
+
+      CALL section_add_subsection(section, subsection)
+      CALL section_release(subsection)
+
 !  The donor state subsection
 
       CALL section_create(subsection, __LOCATION__, name="DONOR_STATES", &
@@ -7989,9 +8060,9 @@ SUBROUTINE create_xas_tdp_section(section)
 
       CALL keyword_create(keyword, __LOCATION__, name="LOCALIZE", &
                           variants=s2a("LOC", "DO_LOC"), &
-                          description="Whether the N_SEARCH potential donor states are should be "// &
+                          description="Whether the N_SEARCH potential donor states should be "// &
                           "actively localized. Necessary in case of excited atoms "// &
-                          "equivalent under symmetry", &
+                          "equivalent under symmetry or GW2X correction.", &
                           usage="LOCALIZE {logical}", &
                           default_l_val=.FALSE., &
                           lone_keyword_l_val=.TRUE.)

src/qs_o3c_methods.F

@@ -396,7 +396,7 @@ END SUBROUTINE contract12_o3c
    SUBROUTINE contract3_o3c(o3c, vec, matrix)
       TYPE(o3c_container_type), POINTER                  :: o3c
       TYPE(o3c_vec_type), DIMENSION(:), POINTER          :: vec
-      TYPE(dbcsr_type), POINTER                          :: matrix
+      TYPE(dbcsr_type)                                   :: matrix
 
       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'contract3_o3c'
 

src/xas_tdp_atom.F

@@ -159,7 +159,7 @@ SUBROUTINE init_xas_atom_env(xas_atom_env, xas_tdp_env, xas_tdp_control, qs_env)
 
       INTEGER                                            :: handle, ikind, natom, nex_atoms, &
                                                             nex_kinds, nkind, nspins
-      LOGICAL                                            :: do_soc, do_xc
+      LOGICAL                                            :: do_xc
       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
 
@@ -174,7 +174,6 @@ SUBROUTINE init_xas_atom_env(xas_atom_env, xas_tdp_env, xas_tdp_control, qs_env)
       nex_kinds = xas_tdp_env%nex_kinds
       nex_atoms = xas_tdp_env%nex_atoms
       do_xc = xas_tdp_control%do_xc
-      do_soc = xas_tdp_control%do_soc
       nspins = 1; IF (xas_tdp_control%do_uks .OR. xas_tdp_control%do_roks) nspins = 2
 
       xas_atom_env%nspins = nspins
@@ -201,16 +200,14 @@ SUBROUTINE init_xas_atom_env(xas_atom_env, xas_tdp_env, xas_tdp_control, qs_env)
 !  Compute the contraction coefficients for spherical orbitals
       DO ikind = 1, nkind
          NULLIFY (xas_atom_env%orb_sphi_so(ikind)%array, xas_atom_env%ri_sphi_so(ikind)%array)
-         IF (do_soc) THEN
-            CALL compute_sphi_so(ikind, "ORB", xas_atom_env%orb_sphi_so(ikind)%array, qs_env)
-         END IF
-         IF (ANY(xas_atom_env%excited_kinds == ikind) .AND. do_xc) THEN
+         CALL compute_sphi_so(ikind, "ORB", xas_atom_env%orb_sphi_so(ikind)%array, qs_env)
+         IF (ANY(xas_atom_env%excited_kinds == ikind)) THEN
             CALL compute_sphi_so(ikind, "RI_XAS", xas_atom_env%ri_sphi_so(ikind)%array, qs_env)
          END IF
       END DO !ikind
 
 !  Compute the coefficients to expand the density in the RI_XAS basis, if requested
-      IF (do_xc) CALL calculate_density_coeffs(xas_atom_env, qs_env)
+      IF (do_xc .AND. (.NOT. xas_tdp_control%xps_only)) CALL calculate_density_coeffs(xas_atom_env, qs_env)
 
       CALL timestop(handle)
 
@@ -3383,7 +3380,6 @@ SUBROUTINE integrate_soc_atoms(matrix_soc, xas_atom_env, qs_env)
          nr = grid%nr
          na = grid%ng_sphere
          ALLOCATE (Vr(nr))
-         !TODO: do we need potential contribution from neighboring atoms ?!?
          CALL calculate_model_potential(Vr, grid, zeff)
 
          !Compute V/(4c^2-2V) and weight it