GW full k-points small cell: now also with perturbative SOC

(important for 2d TMDs where SOC can close bandgap by 0.3 eV)
hfp · Jun 27, 2024 · 1370f7d · 1370f7d
1 parent 0eac0bc
commit 1370f7d
Show file tree

Hide file tree

Showing 11 changed files with 751 additions and 622 deletions.
diff --git a/src/gw_large_cell_gamma.F b/src/gw_large_cell_gamma.F
@@ -56,7 +56,6 @@ MODULE gw_large_cell_gamma
                                               local_dbt_to_global_mat
    USE gw_utils,                        ONLY: analyt_conti_and_print,&
                                               de_init_bs_env,&
-                                              get_VBM_CBM_bandgaps,&
                                               time_to_freq
    USE input_section_types,             ONLY: section_vals_type
    USE kinds,                           ONLY: default_string_length,&
@@ -74,7 +73,8 @@ MODULE gw_large_cell_gamma
    USE particle_types,                  ONLY: particle_type
    USE post_scf_bandstructure_types,    ONLY: post_scf_bandstructure_type
    USE post_scf_bandstructure_utils,    ONLY: MIC_contribution_from_ikp,&
-                                              cfm_ikp_from_fm_Gamma
+                                              cfm_ikp_from_fm_Gamma,&
+                                              get_all_VBM_CBM_bandgaps
    USE qs_environment_types,            ONLY: get_qs_env,&
                                               qs_environment_type
    USE qs_kind_types,                   ONLY: qs_kind_type
@@ -2293,7 +2293,7 @@ SUBROUTINE compute_QP_energies(bs_env, qs_env, fm_Sigma_x_Gamma, fm_Sigma_c_Gamm
 
       END DO ! ispin
 
-      CALL get_VBM_CBM_bandgaps(bs_env)
+      CALL get_all_VBM_CBM_bandgaps(bs_env)
 
       CALL cp_fm_release(fm_Sigma_x_Gamma)
       CALL cp_fm_release(fm_Sigma_c_Gamma_time)

diff --git a/src/gw_small_cell_full_kp.F b/src/gw_small_cell_full_kp.F
@@ -15,9 +15,9 @@ MODULE gw_small_cell_full_kp
    USE cp_cfm_types,                    ONLY: cp_cfm_create,&
                                               cp_cfm_get_info,&
                                               cp_cfm_release,&
+                                              cp_cfm_to_cfm,&
                                               cp_cfm_to_fm,&
-                                              cp_cfm_type,&
-                                              cp_fm_to_cfm
+                                              cp_cfm_type
    USE cp_dbcsr_api,                    ONLY: dbcsr_create,&
                                               dbcsr_distribution_release,&
                                               dbcsr_distribution_type,&
@@ -50,7 +50,6 @@ MODULE gw_small_cell_full_kp
    USE gw_utils,                        ONLY: add_R,&
                                               analyt_conti_and_print,&
                                               de_init_bs_env,&
-                                              get_VBM_CBM_bandgaps,&
                                               is_cell_in_index_to_cell,&
                                               time_to_freq
    USE kinds,                           ONLY: dp
@@ -63,6 +62,7 @@ MODULE gw_small_cell_full_kp
    USE particle_methods,                ONLY: get_particle_set
    USE particle_types,                  ONLY: particle_type
    USE post_scf_bandstructure_types,    ONLY: post_scf_bandstructure_type
+   USE post_scf_bandstructure_utils,    ONLY: get_all_VBM_CBM_bandgaps
    USE qs_environment_types,            ONLY: get_qs_env,&
                                               qs_environment_type
    USE qs_kind_types,                   ONLY: qs_kind_type
@@ -417,8 +417,7 @@ SUBROUTINE G_occ_vir(bs_env, tau, G_S, ispin, occ, vir)
       DO ikp = 1, nkp
 
          ! get C_µn(k)
-         CALL cp_fm_to_cfm(bs_env%fm_mo_coeff_kp(ikp, ispin, 1), &
-                           bs_env%fm_mo_coeff_kp(ikp, ispin, 2), bs_env%cfm_work_mo)
+         CALL cp_cfm_to_cfm(bs_env%cfm_mo_coeff_kp(ikp, ispin), bs_env%cfm_work_mo)
 
          ! G^occ/vir_µλ(i|τ|,k) = sum_n^occ/vir C_µn(k)^* e^(-|(ϵ_nk-ϵ_F)τ|) C_λn(k)
          DO j_col_local = 1, ncol_local
@@ -1763,8 +1762,7 @@ SUBROUTINE compute_QP_energies(bs_env)
          DO ikp = 1, bs_env%nkp_bs_and_DOS
 
             ! 1. get C_µn(k)
-            CALL cp_fm_to_cfm(bs_env%fm_mo_coeff_kp(ikp, ispin, 1), &
-                              bs_env%fm_mo_coeff_kp(ikp, ispin, 2), cfm_mo_coeff)
+            CALL cp_cfm_to_cfm(bs_env%cfm_mo_coeff_kp(ikp, ispin), cfm_mo_coeff)
 
             ! 2. Σ^x_µν(k) = sum_R Σ^x_µν^R e^ikR
             !    Σ^x_nn(k) = sum_µν C^*_µn(k) Σ^x_µν(k) C_νn(k)
@@ -1792,7 +1790,7 @@ SUBROUTINE compute_QP_energies(bs_env)
 
       END DO ! ispin
 
-      CALL get_VBM_CBM_bandgaps(bs_env)
+      CALL get_all_VBM_CBM_bandgaps(bs_env)
 
       CALL cp_cfm_release(cfm_mo_coeff)
 

diff --git a/src/gw_utils.F b/src/gw_utils.F
diff --git a/src/input_cp2k_properties_dft.F b/src/input_cp2k_properties_dft.F
@@ -2064,6 +2064,15 @@ SUBROUTINE create_gw_section(section)
       CALL section_add_keyword(section, keyword)
       CALL keyword_release(keyword)
 
+      CALL keyword_create(keyword, __LOCATION__, name="REGULARIZATION_RI", &
+                          description="Parameter for RI regularization, setting a negative "// &
+                          "value triggers the default value. Affects RI basis set convergence "// &
+                          "but in any case large RI basis will give RI basis set convergence.", &
+                          usage="REGULARIZATION_RI 1.0E-4", &
+                          default_r_val=-1.0_dp)
+      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 "// &

diff --git a/src/post_scf_bandstructure_methods.F b/src/post_scf_bandstructure_methods.F
@@ -8,13 +8,11 @@
 MODULE post_scf_bandstructure_methods
    USE gw_main,                         ONLY: gw
    USE input_section_types,             ONLY: section_vals_type
-   USE post_scf_bandstructure_types,    ONLY: post_scf_bandstructure_type
    USE post_scf_bandstructure_utils,    ONLY: create_and_init_bs_env,&
-                                              dos_pdos_ldos
+                                              dos_pdos_ldos,&
+                                              soc
    USE qs_environment_types,            ONLY: qs_environment_type
    USE qs_scf,                          ONLY: scf
-   USE soc_pseudopotential_methods,     ONLY: H_KS_spinor,&
-                                              V_SOC_xyz_from_pseudopotential
 #include "./base/base_uses.f90"
 
    IMPLICIT NONE
@@ -64,33 +62,4 @@ SUBROUTINE post_scf_bandstructure(qs_env, post_scf_bandstructure_section)
 
    END SUBROUTINE post_scf_bandstructure
 
-! **************************************************************************************************
-!> \brief ...
-!> \param qs_env ...
-!> \param bs_env ...
-! **************************************************************************************************
-   SUBROUTINE soc(qs_env, bs_env)
-      TYPE(qs_environment_type), POINTER                 :: qs_env
-      TYPE(post_scf_bandstructure_type), POINTER         :: bs_env
-
-      CHARACTER(LEN=*), PARAMETER                        :: routineN = 'soc'
-
-      INTEGER                                            :: handle
-
-      CALL timeset(routineN, handle)
-
-      ! Compute V^SOC_µν^(α) = ħ/2 < ϕ_µ | sum_ℓ ΔV_ℓ^SO(r,r') L^(α) | ϕ_ν >, α = x, y, z, see
-      ! Hartwigsen, Goedecker, Hutter, Eq.(18), (19) (doi.org/10.1103/PhysRevB.58.3641)
-      CALL V_SOC_xyz_from_pseudopotential(qs_env, bs_env%mat_V_SOC_xyz)
-
-      ! Spinor KS-matrix H_µν,σσ' = h^SCF_µν*δ_σσ' + sum_α V^SOC_µν^(α)*Pauli-matrix^(α)_σσ', see
-      ! Hartwigsen, Goedecker, Hutter, Eq.(18) (doi.org/10.1103/PhysRevB.58.3641)
-      CALL H_KS_spinor(bs_env%cfm_ks_spinor_ao_Gamma, bs_env%fm_ks_Gamma(1:2), bs_env%n_spin, &
-                       bs_env%mat_V_SOC_xyz(:, 1), bs_env%cfm_s_spinor_Gamma, bs_env%fm_s_Gamma, &
-                       bs_env%cfm_SOC_spinor_ao_Gamma)
-
-      CALL timestop(handle)
-
-   END SUBROUTINE soc
-
 END MODULE post_scf_bandstructure_methods
diff --git a/src/post_scf_bandstructure_types.F b/src/post_scf_bandstructure_types.F
@@ -142,6 +142,7 @@ MODULE post_scf_bandstructure_types
       ! threshold for inverting ao overlap matrix, RI cfm_1d
       REAL(KIND=dp)                                :: eps_eigval_mat_s = -1.0_dp, &
                                                       eps_eigval_mat_RI = -1.0_dp, &
+                                                      input_regularization_RI = -1.0_dp, &
                                                       regularization_RI = -1.0_dp
 
       ! global full cfm_1d used in GW
@@ -231,11 +232,11 @@ MODULE post_scf_bandstructure_types
 
       ! parameters for SOC calculation
       REAL(KIND=dp)                                :: energy_window_soc = -1.0_dp
+      ! sizes: mat_V_SOC_xyz: xyz, img
       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: mat_V_SOC_xyz => NULL()
       TYPE(cp_fm_type), DIMENSION(3)               :: fm_V_SOC_xyz_mo
-      TYPE(cp_cfm_type)                            :: cfm_ks_spinor_ao_Gamma, &
-                                                      cfm_SOC_spinor_ao_Gamma, &
-                                                      cfm_s_spinor_Gamma
+      ! small-cell GW: dimension = number of kpoints; large-cell GW: Gamma-point, dimension = 1
+      TYPE(cp_cfm_type), DIMENSION(:), ALLOCATABLE :: cfm_SOC_spinor_ao
       TYPE(band_edges_type)                        :: band_edges_scf_SOC, &
                                                       band_edges_G0W0_SOC
 
@@ -269,10 +270,10 @@ MODULE post_scf_bandstructure_types
       INTEGER, DIMENSION(:), ALLOCATABLE           :: task_Delta_R
       INTEGER, DIMENSION(:, :), ALLOCATABLE        :: nblocks_3c
 
-      ! full complex fm for k-dependent mo coefficients C_μn(k,spin,re/im) from SCF
-      TYPE(cp_fm_type), DIMENSION(:, :, :), ALLOCATABLE :: fm_mo_coeff_kp, &
-                                                           fm_ks_kp, &
-                                                           fm_s_kp
+      ! cfm for k-dep overl mat S_µν(k), KS mat h_µν(k,spin) and mo coeff C_μn(k,spin) from SCF
+      TYPE(cp_cfm_type), DIMENSION(:), ALLOCATABLE      :: cfm_s_kp
+      TYPE(cp_cfm_type), DIMENSION(:, :), ALLOCATABLE   :: cfm_mo_coeff_kp, &
+                                                           cfm_ks_kp
       TYPE(cp_fm_type), DIMENSION(:), ALLOCATABLE       :: fm_G_S, &
                                                            fm_Sigma_x_R
       TYPE(cp_fm_type), DIMENSION(:, :), ALLOCATABLE    :: fm_V_xc_R, &
@@ -371,9 +372,7 @@ SUBROUTINE bs_env_release(bs_env)
       CALL safe_fm_destroy_3d(bs_env%fm_Sigma_c_R_neg_tau)
       CALL safe_fm_destroy_3d(bs_env%fm_Sigma_c_R_pos_tau)
 
-      IF (bs_env%small_cell_full_kp_or_large_cell_Gamma == small_cell_full_kp) THEN
-         CALL t_destroy_2d(bs_env%t_3c_int)
-      END IF
+      CALL t_destroy_2d(bs_env%t_3c_int)
 
       CALL release_dbcsr_p_type(bs_env%mat_ao_ao)
       CALL release_dbcsr_p_type(bs_env%mat_RI_RI)
@@ -382,9 +381,9 @@ SUBROUTINE bs_env_release(bs_env)
       CALL release_dbcsr_p_type(bs_env%mat_ao_ao_tensor)
       CALL release_dbcsr_p_type(bs_env%mat_RI_RI_tensor)
 
-      CALL safe_fm_destroy_3d(bs_env%fm_ks_kp)
-      CALL safe_fm_destroy_3d(bs_env%fm_s_kp)
-      CALL safe_fm_destroy_3d(bs_env%fm_mo_coeff_kp)
+      CALL safe_cfm_destroy_1d(bs_env%cfm_s_kp)
+      CALL safe_cfm_destroy_2d(bs_env%cfm_ks_kp)
+      CALL safe_cfm_destroy_2d(bs_env%cfm_mo_coeff_kp)
 
       CALL mp_para_env_release(bs_env%para_env)
       IF (ASSOCIATED(bs_env%para_env_tensor)) CALL mp_para_env_release(bs_env%para_env_tensor)
@@ -403,9 +402,7 @@ SUBROUTINE bs_env_release(bs_env)
       CALL cp_fm_release(bs_env%fm_V_SOC_xyz_mo(1))
       CALL cp_fm_release(bs_env%fm_V_SOC_xyz_mo(2))
       CALL cp_fm_release(bs_env%fm_V_SOC_xyz_mo(3))
-      CALL cp_cfm_release(bs_env%cfm_ks_spinor_ao_Gamma)
-      CALL cp_cfm_release(bs_env%cfm_SOC_spinor_ao_Gamma)
-      CALL cp_cfm_release(bs_env%cfm_s_spinor_Gamma)
+      CALL safe_cfm_destroy_1d(bs_env%cfm_SOC_spinor_ao)
 
       DEALLOCATE (bs_env)
 
@@ -531,6 +528,44 @@ SUBROUTINE safe_fm_destroy_3d(fm_3d)
 
    END SUBROUTINE safe_fm_destroy_3d
 
+! **************************************************************************************************
+!> \brief ...
+!> \param cfm_1d ...
+! **************************************************************************************************
+   SUBROUTINE safe_cfm_destroy_1d(cfm_1d)
+      TYPE(cp_cfm_type), ALLOCATABLE, DIMENSION(:)       :: cfm_1d
+
+      INTEGER                                            :: i
+
+      IF (ALLOCATED(cfm_1d)) THEN
+         DO i = 1, SIZE(cfm_1d, 1)
+            CALL cp_cfm_release(cfm_1d(i))
+         END DO
+         DEALLOCATE (cfm_1d)
+      END IF
+
+   END SUBROUTINE safe_cfm_destroy_1d
+
+! **************************************************************************************************
+!> \brief ...
+!> \param cfm_2d ...
+! **************************************************************************************************
+   SUBROUTINE safe_cfm_destroy_2d(cfm_2d)
+      TYPE(cp_cfm_type), ALLOCATABLE, DIMENSION(:, :)    :: cfm_2d
+
+      INTEGER                                            :: i, j
+
+      IF (ALLOCATED(cfm_2d)) THEN
+         DO i = 1, SIZE(cfm_2d, 1)
+         DO j = 1, SIZE(cfm_2d, 2)
+            CALL cp_cfm_release(cfm_2d(i, j))
+         END DO
+         END DO
+         DEALLOCATE (cfm_2d)
+      END IF
+
+   END SUBROUTINE safe_cfm_destroy_2d
+
 ! **************************************************************************************************
 !> \brief ...
 !> \param t_2d ...
@@ -546,8 +581,8 @@ SUBROUTINE t_destroy_2d(t_2d)
                CALL dbt_destroy(t_2d(i, j))
             END DO
          END DO
+         DEALLOCATE (t_2d)
       END IF
-      DEALLOCATE (t_2d)
 
    END SUBROUTINE t_destroy_2d