diff --git a/src/gw_large_cell_gamma.F b/src/gw_large_cell_gamma.F index 10c166af61..06e74b449b 100644 --- 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 index aba3fc0431..282f01a2b0 100644 --- 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 index a976dd1745..e0a5e1ecba 100644 --- a/src/gw_utils.F +++ b/src/gw_utils.F @@ -22,17 +22,15 @@ MODULE gw_utils USE cp_blacs_env, ONLY: cp_blacs_env_create,& cp_blacs_env_release,& cp_blacs_env_type - USE cp_cfm_diag, ONLY: cp_cfm_geeig_canon USE cp_cfm_types, ONLY: cp_cfm_create,& cp_cfm_release,& + cp_cfm_to_cfm,& cp_cfm_to_fm,& - cp_cfm_type,& - cp_fm_to_cfm + cp_cfm_type USE cp_control_types, ONLY: dft_control_type - USE cp_dbcsr_api, ONLY: & - dbcsr_create, dbcsr_deallocate_matrix, dbcsr_desymmetrize, dbcsr_p_type, dbcsr_set, & - dbcsr_type, dbcsr_type_antisymmetric, dbcsr_type_no_symmetry, dbcsr_type_symmetric - USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl + USE cp_dbcsr_api, ONLY: dbcsr_create,& + dbcsr_p_type,& + dbcsr_set USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,& copy_fm_to_dbcsr,& dbcsr_allocate_matrix_set,& @@ -68,8 +66,7 @@ MODULE gw_utils USE kinds, ONLY: default_string_length,& dp,& int_8 - USE kpoint_methods, ONLY: kpoint_init_cell_index,& - rskp_transform + USE kpoint_methods, ONLY: kpoint_init_cell_index USE kpoint_types, ONLY: get_kpoint_info,& kpoint_create,& kpoint_type @@ -97,9 +94,8 @@ MODULE gw_utils USE particle_types, ONLY: particle_type USE physcon, ONLY: angstrom,& evolt - USE post_scf_bandstructure_types, ONLY: band_edges_type,& - post_scf_bandstructure_type - USE post_scf_bandstructure_utils, ONLY: get_fname + USE post_scf_bandstructure_types, ONLY: post_scf_bandstructure_type + USE post_scf_bandstructure_utils, ONLY: rsmat_to_kp USE qs_energy_types, ONLY: qs_energy_type USE qs_environment_types, ONLY: get_qs_env,& qs_env_part_release,& @@ -126,7 +122,7 @@ MODULE gw_utils PRIVATE - PUBLIC :: create_and_init_bs_env_for_gw, de_init_bs_env, get_i_j_atoms, get_VBM_CBM_bandgaps, & + PUBLIC :: create_and_init_bs_env_for_gw, de_init_bs_env, get_i_j_atoms, & kpoint_init_cell_index_simple, compute_xkp, time_to_freq, analyt_conti_and_print, & add_R, is_cell_in_index_to_cell @@ -198,8 +194,6 @@ SUBROUTINE create_and_init_bs_env_for_gw(qs_env, bs_env, bs_sec) CALL setup_parallelization_Delta_R(bs_env) - CALL compute_cfm_mo_coeff_kp_and_eigenval_scf_kp(qs_env, bs_env) - CALL allocate_matrices_small_cell_full_kp(qs_env, bs_env) CALL trafo_V_xc_R_to_kp(qs_env, bs_env) @@ -267,6 +261,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, "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) @@ -1494,11 +1489,23 @@ SUBROUTINE set_heuristic_parameters(bs_env, qs_env) bs_env%eps_eigval_mat_RI = 0.0_dp - ! for periodic systems, we use the regularized resolution of the identity - IF (SUM(bs_env%periodic) == 0) THEN - bs_env%regularization_RI = 0.0_dp + IF (bs_env%input_regularization_RI > 0.0_dp) THEN + bs_env%regularization_RI = bs_env%input_regularization_RI ELSE - bs_env%regularization_RI = 1.0E-2_dp + + ! default case: + + ! 1. for periodic systems, we use the regularized resolution of the identity per default + SELECT CASE (bs_env%small_cell_full_kp_or_large_cell_Gamma) + CASE (large_cell_Gamma) + bs_env%regularization_RI = 1.0E-2_dp + CASE (small_cell_full_kp) + bs_env%regularization_RI = 1.0E-3_dp + END SELECT + + ! 2. for molecules, no regularization is necessary + IF (SUM(bs_env%periodic) == 0) bs_env%regularization_RI = 0.0_dp + END IF ! truncated Coulomb operator for exchange self-energy @@ -1517,6 +1524,7 @@ 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,A,ES48.1)") "Regularization parameter of RI ", bs_env%regularization_RI END IF CALL timestop(handle) @@ -1949,37 +1957,37 @@ SUBROUTINE setup_cells_3c(qs_env, bs_env) IF (MODULO(cell_pair_count, bs_env%para_env%num_pe) .NE. bs_env%para_env%mepos) CYCLE DO atom_j = 1, bs_env%n_atom - DO atom_k = 1, bs_env%n_atom - DO atom_i = 1, bs_env%n_atom + DO atom_k = 1, bs_env%n_atom + DO atom_i = 1, bs_env%n_atom - j_size = bs_env%i_ao_end_from_atom(atom_j) - bs_env%i_ao_start_from_atom(atom_j) + 1 - k_size = bs_env%i_ao_end_from_atom(atom_k) - bs_env%i_ao_start_from_atom(atom_k) + 1 - i_size = bs_env%i_RI_end_from_atom(atom_i) - bs_env%i_RI_start_from_atom(atom_i) + 1 + j_size = bs_env%i_ao_end_from_atom(atom_j) - bs_env%i_ao_start_from_atom(atom_j) + 1 + k_size = bs_env%i_ao_end_from_atom(atom_k) - bs_env%i_ao_start_from_atom(atom_k) + 1 + i_size = bs_env%i_RI_end_from_atom(atom_i) - bs_env%i_RI_start_from_atom(atom_i) + 1 - ALLOCATE (int_3c(j_size, k_size, i_size)) + ALLOCATE (int_3c(j_size, k_size, i_size)) - ! compute 3-c int. ( μ(atom j) R , ν (atom k) S | P (atom i) 0 ) - ! ("|": truncated Coulomb operator), inside build_3c_integrals: (j k | i) - CALL build_3c_integral_block(int_3c, qs_env, bs_env%ri_metric, & - basis_j=bs_env%basis_set_AO, & - basis_k=bs_env%basis_set_AO, & - basis_i=bs_env%basis_set_RI, & - cell_j=index_to_cell_3c_max(j_cell, 1:3), & - cell_k=index_to_cell_3c_max(k_cell, 1:3), & - atom_k=atom_k, atom_j=atom_j, atom_i=atom_i) + ! compute 3-c int. ( μ(atom j) R , ν (atom k) S | P (atom i) 0 ) + ! ("|": truncated Coulomb operator), inside build_3c_integrals: (j k | i) + CALL build_3c_integral_block(int_3c, qs_env, bs_env%ri_metric, & + basis_j=bs_env%basis_set_AO, & + basis_k=bs_env%basis_set_AO, & + basis_i=bs_env%basis_set_RI, & + cell_j=index_to_cell_3c_max(j_cell, 1:3), & + cell_k=index_to_cell_3c_max(k_cell, 1:3), & + atom_k=atom_k, atom_j=atom_j, atom_i=atom_i) - frobenius_norm = SQRT(SUM(int_3c(:, :, :)**2)) + frobenius_norm = SQRT(SUM(int_3c(:, :, :)**2)) - DEALLOCATE (int_3c) + DEALLOCATE (int_3c) - ! we use a higher threshold here to safe memory when storing the 3c integrals - ! in every tensor group - IF (frobenius_norm > eps) THEN - nblocks_3c_max(j_cell, k_cell) = nblocks_3c_max(j_cell, k_cell) + 1 - END IF + ! we use a higher threshold here to safe memory when storing the 3c integrals + ! in every tensor group + IF (frobenius_norm > eps) THEN + nblocks_3c_max(j_cell, k_cell) = nblocks_3c_max(j_cell, k_cell) + 1 + END IF - END DO - END DO + END DO + END DO END DO END DO @@ -2610,164 +2618,6 @@ SUBROUTINE is_cell_in_index_to_cell(cell, index_to_cell, cell_found) END SUBROUTINE is_cell_in_index_to_cell -! ************************************************************************************************** -!> \brief ... -!> \param qs_env ... -!> \param bs_env ... -! ************************************************************************************************** - SUBROUTINE compute_cfm_mo_coeff_kp_and_eigenval_scf_kp(qs_env, bs_env) - TYPE(qs_environment_type), POINTER :: qs_env - TYPE(post_scf_bandstructure_type), POINTER :: bs_env - - CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_cfm_mo_coeff_kp_and_eigenval_scf_kp' - - INTEGER :: handle, ikp, ispin, nkp_bs_and_DOS, re_im - INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index_scf - REAL(KIND=dp) :: CBM, VBM - REAL(KIND=dp), DIMENSION(3) :: xkp - TYPE(cp_cfm_type) :: cfm_ks, cfm_mos, cfm_s - TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks, matrix_s - TYPE(kpoint_type), POINTER :: kpoints_scf - TYPE(neighbor_list_set_p_type), DIMENSION(:), & - POINTER :: sab_nl - - CALL timeset(routineN, handle) - - CALL get_qs_env(qs_env, & - matrix_ks_kp=matrix_ks, & - matrix_s_kp=matrix_s, & - kpoints=kpoints_scf) - - NULLIFY (sab_nl) - CALL get_kpoint_info(kpoints_scf, sab_nl=sab_nl, cell_to_index=cell_to_index_scf) - - CALL cp_cfm_create(cfm_ks, bs_env%cfm_work_mo%matrix_struct) - CALL cp_cfm_create(cfm_s, bs_env%cfm_work_mo%matrix_struct) - CALL cp_cfm_create(cfm_mos, bs_env%cfm_work_mo%matrix_struct) - - ! nkp_bs_and_DOS contains desymmetrized k-point mesh from SCF and k-points from GW bandstructure - nkp_bs_and_DOS = bs_env%nkp_bs_and_DOS - - ALLOCATE (bs_env%eigenval_G0W0(bs_env%n_ao, nkp_bs_and_DOS, bs_env%n_spin)) - ALLOCATE (bs_env%eigenval_HF(bs_env%n_ao, nkp_bs_and_DOS, bs_env%n_spin)) - ALLOCATE (bs_env%fm_mo_coeff_kp(nkp_bs_and_DOS, bs_env%n_spin, 2)) - ALLOCATE (bs_env%fm_ks_kp(nkp_bs_and_DOS, bs_env%n_spin, 2)) - ALLOCATE (bs_env%fm_s_kp(nkp_bs_and_DOS, bs_env%n_spin, 2)) - DO ikp = 1, nkp_bs_and_DOS - DO ispin = 1, bs_env%n_spin - DO re_im = 1, 2 - CALL cp_fm_create(bs_env%fm_mo_coeff_kp(ikp, ispin, re_im), & - bs_env%cfm_work_mo%matrix_struct) - CALL cp_fm_create(bs_env%fm_ks_kp(ikp, ispin, re_im), & - bs_env%cfm_work_mo%matrix_struct) - CALL cp_fm_create(bs_env%fm_s_kp(ikp, ispin, re_im), & - bs_env%cfm_work_mo%matrix_struct) - END DO - END DO - END DO - - DO ispin = 1, bs_env%n_spin - DO ikp = 1, nkp_bs_and_DOS - - xkp(1:3) = bs_env%kpoints_DOS%xkp(1:3, ikp) - - ! h^KS^R -> h^KS(k) - CALL rsmat_to_kp(matrix_ks, ispin, xkp, cell_to_index_scf, sab_nl, bs_env, cfm_ks) - - ! S^R -> S(k) - CALL rsmat_to_kp(matrix_s, 1, xkp, cell_to_index_scf, sab_nl, bs_env, cfm_s) - - ! JW comment: one might remove h^KS(k) again later - ! we store the complex KS matrix as fm matrix because the infrastructure for fm is - ! much nicer compared to cfm - CALL cp_cfm_to_fm(cfm_ks, & - bs_env%fm_ks_kp(ikp, ispin, 1), & - bs_env%fm_ks_kp(ikp, ispin, 2)) - CALL cp_cfm_to_fm(cfm_s, & - bs_env%fm_s_kp(ikp, ispin, 1), & - bs_env%fm_s_kp(ikp, ispin, 2)) - - ! Diagonalize KS-matrix via Rothaan-Hall equation: - ! H^KS(k) C(k) = S(k) C(k) ε(k) - CALL cp_cfm_geeig_canon(cfm_ks, cfm_s, cfm_mos, & - bs_env%eigenval_scf(:, ikp, ispin), & - bs_env%cfm_work_mo, bs_env%eps_eigval_mat_s) - - ! we store the complex MO coeff as fm matrix because the infrastructure for fm is - ! much nicer compared to cfm - CALL cp_cfm_to_fm(cfm_mos, & - bs_env%fm_mo_coeff_kp(ikp, ispin, 1), & - bs_env%fm_mo_coeff_kp(ikp, ispin, 2)) - - END DO - - VBM = MAXVAL(bs_env%eigenval_scf(bs_env%n_occ(ispin), :, ispin)) - CBM = MINVAL(bs_env%eigenval_scf(bs_env%n_occ(ispin) + 1, :, ispin)) - - bs_env%e_fermi(ispin) = 0.5_dp*(VBM + CBM) - - END DO - - CALL cp_cfm_release(cfm_ks) - CALL cp_cfm_release(cfm_s) - CALL cp_cfm_release(cfm_mos) - - CALL timestop(handle) - - END SUBROUTINE compute_cfm_mo_coeff_kp_and_eigenval_scf_kp - -! ************************************************************************************************** -!> \brief ... -!> \param matrix_ks ... -!> \param ispin ... -!> \param xkp ... -!> \param cell_to_index_scf ... -!> \param sab_nl ... -!> \param bs_env ... -!> \param cfm_ks ... -! ************************************************************************************************** - SUBROUTINE rsmat_to_kp(matrix_ks, ispin, xkp, cell_to_index_scf, sab_nl, bs_env, cfm_ks) - TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks - INTEGER :: ispin - REAL(KIND=dp), DIMENSION(3) :: xkp - INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index_scf - TYPE(neighbor_list_set_p_type), DIMENSION(:), & - POINTER :: sab_nl - TYPE(post_scf_bandstructure_type), POINTER :: bs_env - TYPE(cp_cfm_type) :: cfm_ks - - CHARACTER(LEN=*), PARAMETER :: routineN = 'rsmat_to_kp' - - INTEGER :: handle - TYPE(dbcsr_type), POINTER :: cmat, nsmat, rmat - - CALL timeset(routineN, handle) - - ALLOCATE (rmat, cmat, nsmat) - CALL dbcsr_create(rmat, template=matrix_ks(1, 1)%matrix, matrix_type=dbcsr_type_symmetric) - CALL dbcsr_create(cmat, template=matrix_ks(1, 1)%matrix, matrix_type=dbcsr_type_antisymmetric) - CALL dbcsr_create(nsmat, template=matrix_ks(1, 1)%matrix, matrix_type=dbcsr_type_no_symmetry) - CALL cp_dbcsr_alloc_block_from_nbl(rmat, sab_nl) - CALL cp_dbcsr_alloc_block_from_nbl(cmat, sab_nl) - - CALL dbcsr_set(rmat, 0.0_dp) - CALL dbcsr_set(cmat, 0.0_dp) - CALL rskp_transform(rmatrix=rmat, cmatrix=cmat, rsmat=matrix_ks, ispin=ispin, & - xkp=xkp, cell_to_index=cell_to_index_scf, sab_nl=sab_nl) - CALL dbcsr_desymmetrize(rmat, nsmat) - CALL copy_dbcsr_to_fm(nsmat, bs_env%fm_work_mo(1)) - CALL dbcsr_desymmetrize(cmat, nsmat) - CALL copy_dbcsr_to_fm(nsmat, bs_env%fm_work_mo(2)) - CALL cp_fm_to_cfm(bs_env%fm_work_mo(1), bs_env%fm_work_mo(2), cfm_ks) - - CALL dbcsr_deallocate_matrix(rmat) - CALL dbcsr_deallocate_matrix(cmat) - CALL dbcsr_deallocate_matrix(nsmat) - - CALL timestop(handle) - - END SUBROUTINE rsmat_to_kp - ! ************************************************************************************************** !> \brief ... !> \param qs_env ... @@ -2872,8 +2722,7 @@ SUBROUTINE trafo_V_xc_R_to_kp(qs_env, bs_env) cell_to_index_scf, sab_nl, bs_env, cfm_V_xc) ! 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) ! v^xc_nm(k_i) = sum_µν C^*_µn(k_i) v^xc_µν(k_i) C_νn(k_i) CALL parallel_gemm('N', 'N', n_ao, n_ao, n_ao, z_one, cfm_V_xc, cfm_mo_coeff, & @@ -2983,7 +2832,7 @@ SUBROUTINE analyt_conti_and_print(bs_env, Sigma_c_ikp_n_freq, Sigma_x_ikp_n, V_x CHARACTER(len=3) :: occ_vir CHARACTER(len=default_string_length) :: fname - INTEGER :: handle, i_mo, ikp_for_filename, iunit, & + INTEGER :: handle, i_mo, ikp_for_print, iunit, & n_mo, nkp LOGICAL :: is_bandstruc_kpoint, print_DOS_kpoints, & print_ikp @@ -3007,8 +2856,8 @@ SUBROUTINE analyt_conti_and_print(bs_env, Sigma_c_ikp_n_freq, Sigma_x_ikp_n, V_x Sigma_x_ikp_n(:) - V_xc_ikp_n(:), & eigenval_scf(:), & eigenval_scf(:), & - i_mo, bs_env%n_occ(ispin), bs_env%nparam_pade, & - bs_env%num_freq_points_fit, & + i_mo, bs_env%n_occ(ispin), bs_env%n_vir(ispin), & + bs_env%nparam_pade, bs_env%num_freq_points_fit, & ri_rpa_g0w0_crossing_newton, bs_env%n_occ(ispin), & 0.0_dp, .TRUE., .FALSE., 1, e_fermi_ext=bs_env%e_fermi(ispin)) END DO @@ -3034,28 +2883,34 @@ SUBROUTINE analyt_conti_and_print(bs_env, Sigma_c_ikp_n_freq, Sigma_x_ikp_n, V_x IF (print_DOS_kpoints) THEN nkp = bs_env%nkp_only_DOS - ikp_for_filename = ikp + ikp_for_print = ikp ELSE nkp = bs_env%nkp_only_bs - ikp_for_filename = ikp - bs_env%nkp_only_DOS + ikp_for_print = ikp - bs_env%nkp_only_DOS END IF - CALL get_fname(fname, bs_env, ikp_for_filename, nkp, "SCF_and_G0W0", ispin=ispin) + fname = "bandstructure_SCF_and_G0W0" - CALL open_file(TRIM(fname), unit_number=iunit, file_status="REPLACE", file_action="WRITE") + IF (ikp_for_print == 1) THEN + CALL open_file(TRIM(fname), unit_number=iunit, file_status="REPLACE", & + file_action="WRITE") + ELSE + CALL open_file(TRIM(fname), unit_number=iunit, file_status="OLD", & + file_action="WRITE", file_position="APPEND") + END IF WRITE (iunit, "(A)") " " - - WRITE (iunit, "(A10,3F10.4)") "kpoint: ", bs_env%kpoints_DOS%xkp(:, ikp) + WRITE (iunit, "(A10,I7,A25,3F10.4)") "kpoint: ", ikp_for_print, "coordinate: ", & + bs_env%kpoints_DOS%xkp(:, ikp) WRITE (iunit, "(A)") " " - WRITE (iunit, "(A5,A24,2A17,A16,A18)") "n", "ϵ_nk^DFT (eV)", "Σ^c_nk (eV)", & + WRITE (iunit, "(A5,A12,3A17,A16,A18)") "n", "k", "ϵ_nk^DFT (eV)", "Σ^c_nk (eV)", & "Σ^x_nk (eV)", "v_nk^xc (eV)", "ϵ_nk^G0W0 (eV)" WRITE (iunit, "(A)") " " DO i_mo = 1, n_mo IF (i_mo .LE. bs_env%n_occ(ispin)) occ_vir = 'occ' IF (i_mo > bs_env%n_occ(ispin)) occ_vir = 'vir' - WRITE (iunit, "(I5,3A,4F16.3,F17.3)") i_mo, ' (', occ_vir, ') ', & + WRITE (iunit, "(I5,3A,I5,4F16.3,F17.3)") i_mo, ' (', occ_vir, ') ', ikp_for_print, & eigenval_scf(i_mo)*evolt, & Sigma_c_ikp_n_qp(i_mo)*evolt, & Sigma_x_ikp_n(i_mo)*evolt, & @@ -3063,6 +2918,8 @@ SUBROUTINE analyt_conti_and_print(bs_env, Sigma_c_ikp_n_freq, Sigma_x_ikp_n, V_x bs_env%eigenval_G0W0(i_mo, ikp, ispin)*evolt END DO + WRITE (iunit, "(A)") " " + CALL close_file(iunit) END IF @@ -3130,86 +2987,4 @@ SUBROUTINE correct_obvious_fitting_fails(Sigma_c_ikp_n_qp, ispin, bs_env) END SUBROUTINE correct_obvious_fitting_fails -! ************************************************************************************************** -!> \brief ... -!> \param bs_env ... -! ************************************************************************************************** - SUBROUTINE get_VBM_CBM_bandgaps(bs_env) - - TYPE(post_scf_bandstructure_type), POINTER :: bs_env - - CHARACTER(LEN=*), PARAMETER :: routineN = 'get_VBM_CBM_bandgaps' - - INTEGER :: handle - - CALL timeset(routineN, handle) - - CALL get_VBM_CBM_bandgaps_single(bs_env%band_edges_scf, bs_env%eigenval_scf, bs_env) - CALL get_VBM_CBM_bandgaps_single(bs_env%band_edges_G0W0, bs_env%eigenval_G0W0, bs_env) - CALL get_VBM_CBM_bandgaps_single(bs_env%band_edges_HF, bs_env%eigenval_HF, bs_env) - - CALL timestop(handle) - - END SUBROUTINE get_VBM_CBM_bandgaps - -! ************************************************************************************************** -!> \brief ... -!> \param band_edges ... -!> \param ev ... -!> \param bs_env ... -! ************************************************************************************************** - SUBROUTINE get_VBM_CBM_bandgaps_single(band_edges, ev, bs_env) - TYPE(band_edges_type) :: band_edges - REAL(KIND=dp), DIMENSION(:, :, :) :: ev - TYPE(post_scf_bandstructure_type), POINTER :: bs_env - - CHARACTER(LEN=*), PARAMETER :: routineN = 'get_VBM_CBM_bandgaps_single' - - INTEGER :: handle, homo, homo_1, homo_2, ikp, & - ispin, lumo, lumo_1, lumo_2, n_mo - REAL(KIND=dp) :: E_DBG_at_ikp - - CALL timeset(routineN, handle) - - n_mo = bs_env%n_ao - - band_edges%DBG = 1000.0_dp - - SELECT CASE (bs_env%n_spin) - CASE (1) - homo = bs_env%n_occ(1) - lumo = homo + 1 - band_edges%VBM = MAXVAL(ev(1:homo, :, 1)) - band_edges%CBM = MINVAL(ev(homo + 1:n_mo, :, 1)) - CASE (2) - homo_1 = bs_env%n_occ(1) - lumo_1 = homo_1 + 1 - homo_2 = bs_env%n_occ(2) - lumo_2 = homo_2 + 1 - band_edges%VBM = MAX(MAXVAL(ev(1:homo_1, :, 1)), MAXVAL(ev(1:homo_2, :, 2))) - band_edges%CBM = MIN(MINVAL(ev(homo_1 + 1:n_mo, :, 1)), MINVAL(ev(homo_2 + 1:n_mo, :, 2))) - CASE DEFAULT - CPABORT("Error with number of spins.") - END SELECT - - band_edges%IDBG = band_edges%CBM - band_edges%VBM - - DO ispin = 1, bs_env%n_spin - - homo = bs_env%n_occ(ispin) - - DO ikp = 1, bs_env%nkp_bs_and_DOS - - E_DBG_at_ikp = -MAXVAL(ev(1:homo, ikp, ispin)) + MINVAL(ev(homo + 1:n_mo, ikp, ispin)) - - IF (E_DBG_at_ikp < band_edges%DBG) band_edges%DBG = E_DBG_at_ikp - - END DO - - END DO - - CALL timestop(handle) - - END SUBROUTINE get_VBM_CBM_bandgaps_single - END MODULE gw_utils diff --git a/src/input_cp2k_properties_dft.F b/src/input_cp2k_properties_dft.F index ea3250b1f6..272f3c28e7 100644 --- 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 index b30720151d..8ab305b484 100644 --- 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 index 17410d0153..0be67ae5f3 100644 --- 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 diff --git a/src/post_scf_bandstructure_utils.F b/src/post_scf_bandstructure_utils.F index a8cfd95955..9443eca2c5 100644 --- a/src/post_scf_bandstructure_utils.F +++ b/src/post_scf_bandstructure_utils.F @@ -18,8 +18,10 @@ MODULE post_scf_bandstructure_utils get_cell,& pbc USE cp_blacs_env, ONLY: cp_blacs_env_type - USE cp_cfm_basic_linalg, ONLY: cp_cfm_cholesky_decompose + USE cp_cfm_basic_linalg, ONLY: cp_cfm_cholesky_decompose,& + cp_cfm_scale USE cp_cfm_diag, ONLY: cp_cfm_geeig,& + cp_cfm_geeig_canon,& cp_cfm_heevd USE cp_cfm_types, ONLY: cp_cfm_create,& cp_cfm_get_info,& @@ -30,10 +32,10 @@ MODULE post_scf_bandstructure_utils cp_cfm_type,& cp_fm_to_cfm USE cp_control_types, ONLY: dft_control_type - USE cp_dbcsr_api, ONLY: dbcsr_create,& - dbcsr_p_type,& - dbcsr_type_no_symmetry,& - dbcsr_type_symmetric + USE cp_dbcsr_api, ONLY: & + dbcsr_create, dbcsr_deallocate_matrix, dbcsr_desymmetrize, dbcsr_p_type, dbcsr_set, & + dbcsr_type, dbcsr_type_antisymmetric, dbcsr_type_no_symmetry, dbcsr_type_symmetric + USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,& copy_fm_to_dbcsr,& dbcsr_allocate_matrix_set,& @@ -63,7 +65,8 @@ MODULE post_scf_bandstructure_utils USE kinds, ONLY: default_string_length,& dp,& max_line_length - USE kpoint_methods, ONLY: kpoint_init_cell_index + USE kpoint_methods, ONLY: kpoint_init_cell_index,& + rskp_transform USE kpoint_types, ONLY: get_kpoint_info,& kpoint_create,& kpoint_type @@ -94,9 +97,12 @@ MODULE post_scf_bandstructure_utils USE rpa_gw_im_time_util, ONLY: compute_weight_re_im,& get_atom_index_from_basis_function_index USE scf_control_types, ONLY: scf_control_type - USE soc_pseudopotential_methods, ONLY: remove_soc_outside_energy_window_mo + USE soc_pseudopotential_methods, ONLY: V_SOC_xyz_from_pseudopotential,& + remove_soc_outside_energy_window_mo USE soc_pseudopotential_utils, ONLY: add_cfm_submat,& + add_dbcsr_submat,& cfm_add_on_diag,& + create_cfm_double,& get_cfm_submat #include "base/base_uses.f90" @@ -105,8 +111,9 @@ MODULE post_scf_bandstructure_utils PRIVATE PUBLIC :: create_and_init_bs_env, & - dos_pdos_ldos, cfm_ikp_from_fm_Gamma, get_fname, MIC_contribution_from_ikp, & - compute_xkp, kpoint_init_cell_index_simple + dos_pdos_ldos, cfm_ikp_from_fm_Gamma, MIC_contribution_from_ikp, & + compute_xkp, kpoint_init_cell_index_simple, rsmat_to_kp, soc, & + get_VBM_CBM_bandgaps, get_all_VBM_CBM_bandgaps CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'post_scf_bandstructure_utils' @@ -144,6 +151,12 @@ SUBROUTINE create_and_init_bs_env(qs_env, bs_env, post_scf_bandstructure_section CALL setup_kpoints_DOS_large_cell_Gamma(qs_env, bs_env, bs_env%kpoints_DOS) + CALL allocate_and_fill_fm_ks_fm_s(qs_env, bs_env) + + CALL diagonalize_ks_matrix(bs_env) + + CALL check_positive_definite_overlap_mat(bs_env, qs_env) + CASE (small_cell_full_kp) CALL setup_kpoints_scf_desymm(qs_env, bs_env, bs_env%kpoints_scf_desymm, .TRUE.) @@ -151,13 +164,11 @@ SUBROUTINE create_and_init_bs_env(qs_env, bs_env, post_scf_bandstructure_section CALL setup_kpoints_DOS_small_cell_full_kp(bs_env, bs_env%kpoints_DOS) - END SELECT - - CALL allocate_and_fill_fm_ks_fm_s(qs_env, bs_env) + CALL allocate_and_fill_fm_ks_fm_s(qs_env, bs_env) - CALL diagonalize_ks_matrix(bs_env) + CALL compute_cfm_mo_coeff_kp_and_eigenval_scf_kp(qs_env, bs_env) - CALL check_positive_definite_overlap_mat(bs_env, qs_env) + END SELECT CALL timestop(handle) @@ -542,6 +553,166 @@ SUBROUTINE setup_kpoints_DOS_small_cell_full_kp(bs_env, kpoints) END SUBROUTINE setup_kpoints_DOS_small_cell_full_kp +! ************************************************************************************************** +!> \brief ... +!> \param qs_env ... +!> \param bs_env ... +! ************************************************************************************************** + SUBROUTINE compute_cfm_mo_coeff_kp_and_eigenval_scf_kp(qs_env, bs_env) + TYPE(qs_environment_type), POINTER :: qs_env + TYPE(post_scf_bandstructure_type), POINTER :: bs_env + + CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_cfm_mo_coeff_kp_and_eigenval_scf_kp' + + INTEGER :: handle, ikp, ispin, nkp_bs_and_DOS + INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index_scf + REAL(KIND=dp) :: CBM, VBM + REAL(KIND=dp), DIMENSION(3) :: xkp + TYPE(cp_cfm_type) :: cfm_ks, cfm_mos, cfm_s + TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks, matrix_s + TYPE(kpoint_type), POINTER :: kpoints_scf + TYPE(neighbor_list_set_p_type), DIMENSION(:), & + POINTER :: sab_nl + + CALL timeset(routineN, handle) + + CALL get_qs_env(qs_env, & + matrix_ks_kp=matrix_ks, & + matrix_s_kp=matrix_s, & + kpoints=kpoints_scf) + + NULLIFY (sab_nl) + CALL get_kpoint_info(kpoints_scf, sab_nl=sab_nl, cell_to_index=cell_to_index_scf) + + CALL cp_cfm_create(cfm_ks, bs_env%cfm_work_mo%matrix_struct) + CALL cp_cfm_create(cfm_s, bs_env%cfm_work_mo%matrix_struct) + CALL cp_cfm_create(cfm_mos, bs_env%cfm_work_mo%matrix_struct) + + ! nkp_bs_and_DOS contains desymmetrized k-point mesh from SCF and k-points from GW bandstructure + nkp_bs_and_DOS = bs_env%nkp_bs_and_DOS + + ALLOCATE (bs_env%eigenval_G0W0(bs_env%n_ao, nkp_bs_and_DOS, bs_env%n_spin)) + ALLOCATE (bs_env%eigenval_HF(bs_env%n_ao, nkp_bs_and_DOS, bs_env%n_spin)) + ALLOCATE (bs_env%cfm_mo_coeff_kp(nkp_bs_and_DOS, bs_env%n_spin)) + ALLOCATE (bs_env%cfm_ks_kp(nkp_bs_and_DOS, bs_env%n_spin)) + ALLOCATE (bs_env%cfm_s_kp(nkp_bs_and_DOS)) + DO ikp = 1, nkp_bs_and_DOS + DO ispin = 1, bs_env%n_spin + CALL cp_cfm_create(bs_env%cfm_mo_coeff_kp(ikp, ispin), bs_env%cfm_work_mo%matrix_struct) + CALL cp_cfm_create(bs_env%cfm_ks_kp(ikp, ispin), bs_env%cfm_work_mo%matrix_struct) + END DO + CALL cp_cfm_create(bs_env%cfm_s_kp(ikp), bs_env%cfm_work_mo%matrix_struct) + END DO + + DO ispin = 1, bs_env%n_spin + DO ikp = 1, nkp_bs_and_DOS + + xkp(1:3) = bs_env%kpoints_DOS%xkp(1:3, ikp) + + ! h^KS^R -> h^KS(k) + CALL rsmat_to_kp(matrix_ks, ispin, xkp, cell_to_index_scf, sab_nl, bs_env, cfm_ks) + + ! S^R -> S(k) + CALL rsmat_to_kp(matrix_s, 1, xkp, cell_to_index_scf, sab_nl, bs_env, cfm_s) + + ! we store the complex KS matrix as fm matrix because the infrastructure for fm is + ! much nicer compared to cfm + CALL cp_cfm_to_cfm(cfm_ks, bs_env%cfm_ks_kp(ikp, ispin)) + CALL cp_cfm_to_cfm(cfm_s, bs_env%cfm_s_kp(ikp)) + + ! Diagonalize KS-matrix via Rothaan-Hall equation: + ! H^KS(k) C(k) = S(k) C(k) ε(k) + CALL cp_cfm_geeig_canon(cfm_ks, cfm_s, cfm_mos, & + bs_env%eigenval_scf(:, ikp, ispin), & + bs_env%cfm_work_mo, bs_env%eps_eigval_mat_s) + + ! we store the complex MO coeff as fm matrix because the infrastructure for fm is + ! much nicer compared to cfm + CALL cp_cfm_to_cfm(cfm_mos, bs_env%cfm_mo_coeff_kp(ikp, ispin)) + + END DO + + VBM = MAXVAL(bs_env%eigenval_scf(bs_env%n_occ(ispin), :, ispin)) + CBM = MINVAL(bs_env%eigenval_scf(bs_env%n_occ(ispin) + 1, :, ispin)) + + bs_env%e_fermi(ispin) = 0.5_dp*(VBM + CBM) + + END DO + + CALL get_VBM_CBM_bandgaps(bs_env%band_edges_scf, bs_env%eigenval_scf, bs_env) + + CALL cp_cfm_release(cfm_ks) + CALL cp_cfm_release(cfm_s) + CALL cp_cfm_release(cfm_mos) + + CALL timestop(handle) + + END SUBROUTINE compute_cfm_mo_coeff_kp_and_eigenval_scf_kp + +! ************************************************************************************************** +!> \brief ... +!> \param matrix_ks ... +!> \param ispin ... +!> \param xkp ... +!> \param cell_to_index_scf ... +!> \param sab_nl ... +!> \param bs_env ... +!> \param cfm_ks ... +!> \param imag_rs_mat ... +! ************************************************************************************************** + SUBROUTINE rsmat_to_kp(matrix_ks, ispin, xkp, cell_to_index_scf, sab_nl, bs_env, cfm_ks, imag_rs_mat) + TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks + INTEGER :: ispin + REAL(KIND=dp), DIMENSION(3) :: xkp + INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index_scf + TYPE(neighbor_list_set_p_type), DIMENSION(:), & + POINTER :: sab_nl + TYPE(post_scf_bandstructure_type), POINTER :: bs_env + TYPE(cp_cfm_type) :: cfm_ks + LOGICAL, OPTIONAL :: imag_rs_mat + + CHARACTER(LEN=*), PARAMETER :: routineN = 'rsmat_to_kp' + + INTEGER :: handle + LOGICAL :: imag_rs_mat_private + TYPE(dbcsr_type), POINTER :: cmat, nsmat, rmat + + CALL timeset(routineN, handle) + + ALLOCATE (rmat, cmat, nsmat) + + imag_rs_mat_private = .FALSE. + IF (PRESENT(imag_rs_mat)) imag_rs_mat_private = imag_rs_mat + + IF (imag_rs_mat_private) THEN + CALL dbcsr_create(rmat, template=matrix_ks(1, 1)%matrix, matrix_type=dbcsr_type_antisymmetric) + CALL dbcsr_create(cmat, template=matrix_ks(1, 1)%matrix, matrix_type=dbcsr_type_symmetric) + ELSE + CALL dbcsr_create(rmat, template=matrix_ks(1, 1)%matrix, matrix_type=dbcsr_type_symmetric) + CALL dbcsr_create(cmat, template=matrix_ks(1, 1)%matrix, matrix_type=dbcsr_type_antisymmetric) + END IF + CALL dbcsr_create(nsmat, template=matrix_ks(1, 1)%matrix, matrix_type=dbcsr_type_no_symmetry) + CALL cp_dbcsr_alloc_block_from_nbl(rmat, sab_nl) + CALL cp_dbcsr_alloc_block_from_nbl(cmat, sab_nl) + + CALL dbcsr_set(rmat, 0.0_dp) + CALL dbcsr_set(cmat, 0.0_dp) + CALL rskp_transform(rmatrix=rmat, cmatrix=cmat, rsmat=matrix_ks, ispin=ispin, & + xkp=xkp, cell_to_index=cell_to_index_scf, sab_nl=sab_nl) + CALL dbcsr_desymmetrize(rmat, nsmat) + CALL copy_dbcsr_to_fm(nsmat, bs_env%fm_work_mo(1)) + CALL dbcsr_desymmetrize(cmat, nsmat) + CALL copy_dbcsr_to_fm(nsmat, bs_env%fm_work_mo(2)) + CALL cp_fm_to_cfm(bs_env%fm_work_mo(1), bs_env%fm_work_mo(2), cfm_ks) + + CALL dbcsr_deallocate_matrix(rmat) + CALL dbcsr_deallocate_matrix(cmat) + CALL dbcsr_deallocate_matrix(nsmat) + + CALL timestop(handle) + + END SUBROUTINE rsmat_to_kp + ! ************************************************************************************************** !> \brief ... !> \param bs_env ... @@ -816,8 +987,8 @@ SUBROUTINE dos_pdos_ldos(qs_env, bs_env) CHARACTER(LEN=*), PARAMETER :: routineN = 'dos_pdos_ldos' INTEGER :: handle, homo, homo_1, homo_2, & - homo_spinor, ikp, ikp_for_filename, & - ispin, n_ao, n_E, nkind, nkp + homo_spinor, ikp, ikp_for_file, ispin, & + n_ao, n_E, nkind, nkp LOGICAL :: is_bandstruc_kpoint, print_DOS_kpoints, & print_ikp REAL(KIND=dp) :: broadening, E_max, E_max_G0W0, E_min, & @@ -845,8 +1016,9 @@ SUBROUTINE dos_pdos_ldos(qs_env, bs_env) energy_step_DOS = bs_env%energy_step_DOS broadening = bs_env%broadening_DOS - ! if we have done GW, we already have the band edges - IF (bs_env%do_gw) THEN + ! if we have done GW or a full kpoint SCF, we already have the band edges + IF (bs_env%do_gw .OR. & + bs_env%small_cell_full_kp_or_large_cell_Gamma == small_cell_full_kp) THEN band_edges_scf = bs_env%band_edges_scf band_edges_scf_guess = band_edges_scf ELSE @@ -926,13 +1098,13 @@ SUBROUTINE dos_pdos_ldos(qs_env, bs_env) IF (bs_env%do_soc) THEN - CALL cp_cfm_create(cfm_mos_ikp_spinor, bs_env%cfm_ks_spinor_ao_Gamma%matrix_struct) - CALL cp_cfm_create(cfm_work_ikp_spinor, bs_env%cfm_ks_spinor_ao_Gamma%matrix_struct) - CALL cp_cfm_create(cfm_s_ikp_spinor_copy, bs_env%cfm_ks_spinor_ao_Gamma%matrix_struct) - CALL cp_cfm_create(cfm_ks_ikp_spinor, bs_env%cfm_ks_spinor_ao_Gamma%matrix_struct) - CALL cp_cfm_create(cfm_SOC_ikp_spinor, bs_env%cfm_ks_spinor_ao_Gamma%matrix_struct) - CALL cp_cfm_create(cfm_s_ikp_spinor, bs_env%cfm_ks_spinor_ao_Gamma%matrix_struct) - CALL cp_cfm_create(cfm_spinor_wf_ikp, bs_env%cfm_ks_spinor_ao_Gamma%matrix_struct) + CALL cp_cfm_create(cfm_mos_ikp_spinor, bs_env%cfm_SOC_spinor_ao(1)%matrix_struct) + CALL cp_cfm_create(cfm_work_ikp_spinor, bs_env%cfm_SOC_spinor_ao(1)%matrix_struct) + CALL cp_cfm_create(cfm_s_ikp_spinor_copy, bs_env%cfm_SOC_spinor_ao(1)%matrix_struct) + CALL cp_cfm_create(cfm_ks_ikp_spinor, bs_env%cfm_SOC_spinor_ao(1)%matrix_struct) + CALL cp_cfm_create(cfm_SOC_ikp_spinor, bs_env%cfm_SOC_spinor_ao(1)%matrix_struct) + CALL cp_cfm_create(cfm_s_ikp_spinor, bs_env%cfm_SOC_spinor_ao(1)%matrix_struct) + CALL cp_cfm_create(cfm_spinor_wf_ikp, bs_env%cfm_SOC_spinor_ao(1)%matrix_struct) homo_spinor = bs_env%n_occ(1) + bs_env%n_occ(bs_env%n_spin) @@ -957,8 +1129,8 @@ SUBROUTINE dos_pdos_ldos(qs_env, bs_env) END IF IF (bs_env%small_cell_full_kp_or_large_cell_Gamma == small_cell_full_kp) THEN - CALL cp_cfm_create(cfm_ks_ikp, bs_env%fm_ks_kp(1, 1, 1)%matrix_struct) - CALL cp_cfm_create(cfm_s_ikp, bs_env%fm_ks_kp(1, 1, 1)%matrix_struct) + CALL cp_cfm_create(cfm_ks_ikp, bs_env%cfm_ks_kp(1, 1)%matrix_struct) + CALL cp_cfm_create(cfm_s_ikp, bs_env%cfm_ks_kp(1, 1)%matrix_struct) END IF DO ikp = 1, bs_env%nkp_bs_and_DOS @@ -986,14 +1158,13 @@ SUBROUTINE dos_pdos_ldos(qs_env, bs_env) CASE (small_cell_full_kp) ! 1. get H^KS_µν(k_i) - CALL cp_fm_to_cfm(bs_env%fm_ks_kp(ikp, ispin, 1), & - bs_env%fm_ks_kp(ikp, ispin, 2), cfm_ks_ikp) + CALL cp_cfm_to_cfm(bs_env%cfm_ks_kp(ikp, ispin), cfm_ks_ikp) + ! 2. get S_µν(k_i) - CALL cp_fm_to_cfm(bs_env%fm_s_kp(ikp, ispin, 1), & - bs_env%fm_s_kp(ikp, ispin, 2), cfm_s_ikp) + CALL cp_cfm_to_cfm(bs_env%cfm_s_kp(ikp), cfm_s_ikp) + ! 3. get C_µn(k_i) and ϵ_n(k_i) - CALL cp_fm_to_cfm(bs_env%fm_mo_coeff_kp(ikp, ispin, 1), & - bs_env%fm_mo_coeff_kp(ikp, ispin, 2), cfm_mos_ikp(ispin)) + CALL cp_cfm_to_cfm(bs_env%cfm_mo_coeff_kp(ikp, ispin), cfm_mos_ikp(ispin)) eigenval(:) = bs_env%eigenval_scf(:, ikp, ispin) END SELECT @@ -1040,19 +1211,19 @@ SUBROUTINE dos_pdos_ldos(qs_env, bs_env) IF (print_DOS_kpoints) THEN nkp = bs_env%nkp_only_DOS - ikp_for_filename = ikp + ikp_for_file = ikp ELSE nkp = bs_env%nkp_only_bs - ikp_for_filename = ikp - bs_env%nkp_only_DOS + ikp_for_file = ikp - bs_env%nkp_only_DOS END IF ! compute DFT+SOC eigenvalues; based on these, compute band edges, DOS and LDOS - CALL SOC(bs_env, qs_env, ikp, bs_env%eigenval_scf, band_edges_scf, E_min, cfm_mos_ikp, & - DOS_scf_SOC, PDOS_scf_SOC, band_edges_scf_SOC, eigenval_spinor, & - cfm_spinor_wf_ikp) + CALL SOC_ev(bs_env, qs_env, ikp, bs_env%eigenval_scf, band_edges_scf, & + E_min, cfm_mos_ikp, DOS_scf_SOC, PDOS_scf_SOC, & + band_edges_scf_SOC, eigenval_spinor, cfm_spinor_wf_ikp) IF (.NOT. bs_env%do_gw .AND. print_ikp) THEN - CALL write_SOC_eigenvalues(eigenval_spinor, "SCF", ikp_for_filename, nkp, bs_env) + CALL write_SOC_eigenvalues(eigenval_spinor, ikp_for_file, ikp, bs_env) END IF IF (bs_env%do_ldos) THEN @@ -1063,15 +1234,15 @@ SUBROUTINE dos_pdos_ldos(qs_env, bs_env) IF (bs_env%do_gw) THEN ! compute G0W0+SOC eigenvalues; based on these, compute band edges, DOS and LDOS - CALL SOC(bs_env, qs_env, ikp, bs_env%eigenval_G0W0, band_edges_G0W0, & - E_min, cfm_mos_ikp, DOS_G0W0_SOC, PDOS_G0W0_SOC, & - band_edges_G0W0_SOC, eigenval_spinor_G0W0, cfm_spinor_wf_ikp) + CALL SOC_ev(bs_env, qs_env, ikp, bs_env%eigenval_G0W0, band_edges_G0W0, & + E_min, cfm_mos_ikp, DOS_G0W0_SOC, PDOS_G0W0_SOC, & + band_edges_G0W0_SOC, eigenval_spinor_G0W0, cfm_spinor_wf_ikp) IF (print_ikp) THEN ! write SCF+SOC and G0W0+SOC eigenvalues to file ! SCF_and_G0W0_band_structure_for_kpoint__+_SOC - CALL write_SOC_eigenvalues(eigenval_spinor, "SCF_and_G0W0", ikp_for_filename, & - nkp, bs_env, eigenval_spinor_G0W0) + CALL write_SOC_eigenvalues(eigenval_spinor, ikp_for_file, ikp, bs_env, & + eigenval_spinor_G0W0) END IF END IF ! do_gw @@ -1388,8 +1559,8 @@ END SUBROUTINE get_print_format !> \param eigenval_spinor ... !> \param cfm_spinor_wf_ikp ... ! ************************************************************************************************** - SUBROUTINE SOC(bs_env, qs_env, ikp, eigenval_no_SOC, band_edges_no_SOC, E_min, cfm_mos_ikp, & - DOS, PDOS, band_edges, eigenval_spinor, cfm_spinor_wf_ikp) + SUBROUTINE SOC_ev(bs_env, qs_env, ikp, eigenval_no_SOC, band_edges_no_SOC, E_min, cfm_mos_ikp, & + DOS, PDOS, band_edges, eigenval_spinor, cfm_spinor_wf_ikp) TYPE(post_scf_bandstructure_type), POINTER :: bs_env TYPE(qs_environment_type), POINTER :: qs_env @@ -1404,7 +1575,7 @@ SUBROUTINE SOC(bs_env, qs_env, ikp, eigenval_no_SOC, band_edges_no_SOC, E_min, c REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eigenval_spinor TYPE(cp_cfm_type) :: cfm_spinor_wf_ikp - CHARACTER(LEN=*), PARAMETER :: routineN = 'SOC' + CHARACTER(LEN=*), PARAMETER :: routineN = 'SOC_ev' INTEGER :: handle, homo_spinor, n_ao, n_E, nkind REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eigenval_spinor_no_SOC @@ -1420,22 +1591,33 @@ SUBROUTINE SOC(bs_env, qs_env, ikp, eigenval_no_SOC, band_edges_no_SOC, E_min, c n_E = SIZE(DOS) nkind = SIZE(PDOS, 2) - CALL cp_cfm_create(cfm_ks_ikp_spinor, bs_env%cfm_ks_spinor_ao_Gamma%matrix_struct) - CALL cp_cfm_create(cfm_SOC_ikp_spinor, bs_env%cfm_ks_spinor_ao_Gamma%matrix_struct) - CALL cp_cfm_create(cfm_mos_ikp_spinor, bs_env%cfm_ks_spinor_ao_Gamma%matrix_struct) - CALL cp_cfm_create(cfm_work_ikp_spinor, bs_env%cfm_ks_spinor_ao_Gamma%matrix_struct) - CALL cp_cfm_create(cfm_eigenvec_ikp_spinor, bs_env%cfm_ks_spinor_ao_Gamma%matrix_struct) + CALL cp_cfm_create(cfm_ks_ikp_spinor, bs_env%cfm_SOC_spinor_ao(1)%matrix_struct) + CALL cp_cfm_create(cfm_SOC_ikp_spinor, bs_env%cfm_SOC_spinor_ao(1)%matrix_struct) + CALL cp_cfm_create(cfm_mos_ikp_spinor, bs_env%cfm_SOC_spinor_ao(1)%matrix_struct) + CALL cp_cfm_create(cfm_work_ikp_spinor, bs_env%cfm_SOC_spinor_ao(1)%matrix_struct) + CALL cp_cfm_create(cfm_eigenvec_ikp_spinor, bs_env%cfm_SOC_spinor_ao(1)%matrix_struct) ALLOCATE (eigenval_spinor_no_SOC(2*n_ao)) ALLOCATE (proj_mo_on_kind_spinor(2*n_ao, nkind)) ! PDOS not yet implemented -> projection is just zero -> PDOS is zero proj_mo_on_kind_spinor(:, :) = 0.0_dp - ! 1. get V^SOC_µν,σσ'(k_i) from V^SOC_µν,σσ'(k=0) - CALL cfm_ikp_from_cfm_spinor_Gamma(cfm_SOC_ikp_spinor, & - bs_env%cfm_SOC_spinor_ao_Gamma, & - bs_env%fm_s_Gamma%matrix_struct, & - ikp, qs_env, bs_env%kpoints_DOS, "ORB") + ! 1. get V^SOC_µν,σσ'(k_i) + SELECT CASE (bs_env%small_cell_full_kp_or_large_cell_Gamma) + CASE (large_cell_Gamma) + + ! 1. get V^SOC_µν,σσ'(k_i) from V^SOC_µν,σσ'(k=0) + CALL cfm_ikp_from_cfm_spinor_Gamma(cfm_SOC_ikp_spinor, & + bs_env%cfm_SOC_spinor_ao(1), & + bs_env%fm_s_Gamma%matrix_struct, & + ikp, qs_env, bs_env%kpoints_DOS, "ORB") + + CASE (small_cell_full_kp) + + ! 1. V^SOC_µν,σσ'(k_i) already there + CALL cp_cfm_to_cfm(bs_env%cfm_SOC_spinor_ao(ikp), cfm_SOC_ikp_spinor) + + END SELECT ! 2. V^SOC_nn',σσ'(k_i) = sum_µν C^*_µn,σ(k_i) V^SOC_µν,σσ'(k_i) C_νn'(k_i), ! C_µn,σ(k_i): MO coefficiencts from diagonalizing KS-matrix h^KS_nn',σσ'(k_i) @@ -1497,7 +1679,7 @@ SUBROUTINE SOC(bs_env, qs_env, ikp, eigenval_no_SOC, band_edges_no_SOC, E_min, c CALL timestop(handle) - END SUBROUTINE SOC + END SUBROUTINE SOC_ev ! ************************************************************************************************** !> \brief ... @@ -1766,17 +1948,15 @@ END SUBROUTINE add_to_LDOS_2d ! ************************************************************************************************** !> \brief ... !> \param eigenval_spinor ... -!> \param scf_gw ... +!> \param ikp_for_file ... !> \param ikp ... -!> \param nkp ... !> \param bs_env ... !> \param eigenval_spinor_G0W0 ... ! ************************************************************************************************** - SUBROUTINE write_SOC_eigenvalues(eigenval_spinor, scf_gw, ikp, nkp, bs_env, eigenval_spinor_G0W0) + SUBROUTINE write_SOC_eigenvalues(eigenval_spinor, ikp_for_file, ikp, bs_env, eigenval_spinor_G0W0) REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eigenval_spinor - CHARACTER(LEN=*) :: scf_gw - INTEGER :: ikp, nkp + INTEGER :: ikp_for_file, ikp TYPE(post_scf_bandstructure_type), POINTER :: bs_env REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), OPTIONAL :: eigenval_spinor_G0W0 @@ -1788,16 +1968,31 @@ SUBROUTINE write_SOC_eigenvalues(eigenval_spinor, scf_gw, ikp, nkp, bs_env, eige CALL timeset(routineN, handle) - CALL get_fname(fname, bs_env, ikp, nkp, scf_gw, SOC=.TRUE.) + fname = "bandstructure_SCF_and_G0W0_plus_SOC" IF (bs_env%para_env%is_source()) THEN - CALL open_file(TRIM(fname), unit_number=iunit, file_status="REPLACE", file_action="WRITE") + IF (ikp_for_file == 1) THEN + CALL open_file(TRIM(fname), unit_number=iunit, file_status="REPLACE", & + file_action="WRITE") + ELSE + CALL open_file(TRIM(fname), unit_number=iunit, file_status="OLD", & + file_action="WRITE", file_position="APPEND") + END IF WRITE (iunit, "(A)") " " - WRITE (iunit, "(A10,3F10.4)") "kpoint: ", bs_env%kpoints_DOS%xkp(:, ikp) + WRITE (iunit, "(A10,I7,A25,3F10.4)") "kpoint: ", ikp_for_file, "coordinate: ", & + bs_env%kpoints_DOS%xkp(:, ikp) WRITE (iunit, "(A)") " " + IF (PRESENT(eigenval_spinor_G0W0)) THEN + ! SCF+SOC and G0W0+SOC eigenvalues + WRITE (iunit, "(A5,A12,2A22)") "n", "k", "ϵ_nk^DFT+SOC (eV)", "ϵ_nk^G0W0+SOC (eV)" + ELSE + ! SCF+SOC eigenvalues only + WRITE (iunit, "(A5,A12,A22)") "n", "k", "ϵ_nk^DFT+SOC (eV)" + END IF + n_occ_spinor = bs_env%n_occ(1) + bs_env%n_occ(bs_env%n_spin) DO i_mo = 1, SIZE(eigenval_spinor) @@ -1805,12 +2000,12 @@ SUBROUTINE write_SOC_eigenvalues(eigenval_spinor, scf_gw, ikp, nkp, bs_env, eige IF (i_mo > n_occ_spinor) occ_vir = 'vir' IF (PRESENT(eigenval_spinor_G0W0)) THEN ! SCF+SOC and G0W0+SOC eigenvalues - WRITE (iunit, "(I5,3A,4F16.3,2F17.3)") i_mo, ' (', occ_vir, ') ', & - eigenval_spinor(i_mo)*evolt, eigenval_spinor_G0W0(i_mo)*evolt + WRITE (iunit, "(I5,3A,I5,4F16.3,2F17.3)") i_mo, ' (', occ_vir, ') ', & + ikp_for_file, eigenval_spinor(i_mo)*evolt, eigenval_spinor_G0W0(i_mo)*evolt ELSE ! SCF+SOC eigenvalues only - WRITE (iunit, "(I5,3A,4F16.3,F17.3)") i_mo, ' (', occ_vir, ') ', & - eigenval_spinor(i_mo)*evolt + WRITE (iunit, "(I5,3A,I5,4F16.3,F17.3)") i_mo, ' (', occ_vir, ') ', & + ikp_for_file, eigenval_spinor(i_mo)*evolt END IF END DO @@ -1822,79 +2017,6 @@ SUBROUTINE write_SOC_eigenvalues(eigenval_spinor, scf_gw, ikp, nkp, bs_env, eige END SUBROUTINE write_SOC_eigenvalues -! ************************************************************************************************** -!> \brief ... -!> \param fname ... -!> \param bs_env ... -!> \param ikp ... -!> \param nkp ... -!> \param scf_gw ... -!> \param ispin ... -!> \param SOC ... -! ************************************************************************************************** - SUBROUTINE get_fname(fname, bs_env, ikp, nkp, scf_gw, ispin, SOC) - CHARACTER(len=default_string_length) :: fname - TYPE(post_scf_bandstructure_type), POINTER :: bs_env - INTEGER :: ikp, nkp - CHARACTER(len=*) :: scf_gw - INTEGER, OPTIONAL :: ispin - LOGICAL, OPTIONAL :: SOC - - CHARACTER(LEN=*), PARAMETER :: routineN = 'get_fname' - - CHARACTER(len=1) :: digits_ikp - CHARACTER(len=default_string_length) :: core_name - INTEGER :: handle, n_zeros - LOGICAL :: my_SOC - - CALL timeset(routineN, handle) - - my_SOC = .FALSE. - IF (PRESENT(SOC)) my_SOC = SOC - - n_zeros = count_digits(nkp) - count_digits(ikp) - - WRITE (digits_ikp, '(I1)') count_digits(ikp) - - IF (bs_env%kpoints_DOS%nkp == 1) THEN - ! for molecules - core_name = TRIM(scf_gw)//"_eigenvalues" - ELSE - ! for periodic systems - core_name = TRIM(scf_gw)//"_band_struct_kp" - END IF - - SELECT CASE (n_zeros) - CASE (0) - WRITE (fname, "(2A,I"//digits_ikp//")") TRIM(core_name), "_", ikp - CASE (1) - WRITE (fname, "(2A,I"//digits_ikp//")") TRIM(core_name), "_0", ikp - CASE (2) - WRITE (fname, "(2A,I"//digits_ikp//")") TRIM(core_name), "_00", ikp - CASE (3) - WRITE (fname, "(2A,I"//digits_ikp//")") TRIM(core_name), "_000", ikp - CASE DEFAULT - CPABORT("Too many zeros.") - END SELECT - - ! for molecules, we don't have any k-points; overwrite fname - IF (ikp == 1 .AND. bs_env%kpoints_DOS%nkp == 1) THEN - WRITE (fname, "(A)") TRIM(core_name) - END IF - - IF (my_SOC) THEN - WRITE (fname, "(2A)") TRIM(fname), "_+_SOC" - END IF - - IF (bs_env%n_spin == 2 .AND. .NOT. my_SOC) THEN - CPASSERT(PRESENT(ispin)) - WRITE (fname, "(2A,I1)") TRIM(fname), "_spin_", ispin - END IF - - CALL timestop(handle) - - END SUBROUTINE get_fname - ! ************************************************************************************************** !> \brief ... !> \param int_number ... @@ -2522,4 +2644,269 @@ SUBROUTINE kpoint_init_cell_index_simple(kpoints, qs_env) END SUBROUTINE kpoint_init_cell_index_simple +! ************************************************************************************************** +!> \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) + + ! V^SOC_µν^(α),R = ħ/2 < ϕ_µ cell O | sum_ℓ ΔV_ℓ^SO(r,r') L^(α) | ϕ_ν cell R>, α = 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) + + ! Calculate H^SOC_µν,σσ'(k) = sum_α V^SOC_µν^(α)(k)*Pauli-matrix^(α)_σσ' + ! see Hartwigsen, Goedecker, Hutter, Eq.(18) (doi.org/10.1103/PhysRevB.58.3641) + SELECT CASE (bs_env%small_cell_full_kp_or_large_cell_Gamma) + CASE (large_cell_Gamma) + + ! H^SOC_µν,σσ' = sum_α V^SOC_µν^(α)*Pauli-matrix^(α)_σσ' + CALL H_KS_spinor_Gamma(bs_env) + + CASE (small_cell_full_kp) + + ! V^SOC_µν^(α),R -> V^SOC_µν^(α)(k); then calculate spinor H^SOC_µν,σσ'(k) (see above) + CALL H_KS_spinor_kp(qs_env, bs_env) + + END SELECT + + CALL timestop(handle) + + END SUBROUTINE soc + +! ************************************************************************************************** +!> \brief ... +!> \param bs_env ... +! ************************************************************************************************** + SUBROUTINE H_KS_spinor_Gamma(bs_env) + + TYPE(post_scf_bandstructure_type), POINTER :: bs_env + + CHARACTER(LEN=*), PARAMETER :: routineN = 'H_KS_spinor_Gamma' + + INTEGER :: handle, nao, s + TYPE(cp_fm_struct_type), POINTER :: str + + CALL timeset(routineN, handle) + + CALL cp_fm_get_info(bs_env%fm_ks_Gamma(1), nrow_global=nao) + + ALLOCATE (bs_env%cfm_SOC_spinor_ao(1)) + CALL create_cfm_double(bs_env%cfm_SOC_spinor_ao(1), fm_orig=bs_env%fm_ks_Gamma(1)) + CALL cp_cfm_set_all(bs_env%cfm_SOC_spinor_ao(1), z_zero) + + str => bs_env%fm_ks_Gamma(1)%matrix_struct + + s = nao + 1 + + ! careful: inside add_dbcsr_submat, mat_V_SOC_xyz is multiplied by i because the real matrix + ! mat_V_SOC_xyz is antisymmetric as V_SOC matrix is purely imaginary and Hermitian + CALL add_dbcsr_submat(bs_env%cfm_SOC_spinor_ao(1), bs_env%mat_V_SOC_xyz(1, 1)%matrix, & + str, s, 1, z_one, .TRUE.) + CALL add_dbcsr_submat(bs_env%cfm_SOC_spinor_ao(1), bs_env%mat_V_SOC_xyz(2, 1)%matrix, & + str, s, 1, gaussi, .TRUE.) + CALL add_dbcsr_submat(bs_env%cfm_SOC_spinor_ao(1), bs_env%mat_V_SOC_xyz(3, 1)%matrix, & + str, 1, 1, z_one, .FALSE.) + CALL add_dbcsr_submat(bs_env%cfm_SOC_spinor_ao(1), bs_env%mat_V_SOC_xyz(3, 1)%matrix, & + str, s, s, -z_one, .FALSE.) + + CALL timestop(handle) + + END SUBROUTINE H_KS_spinor_Gamma + +! ************************************************************************************************** +!> \brief ... +!> \param qs_env ... +!> \param bs_env ... +! ************************************************************************************************** + SUBROUTINE H_KS_spinor_kp(qs_env, bs_env) + TYPE(qs_environment_type), POINTER :: qs_env + TYPE(post_scf_bandstructure_type), POINTER :: bs_env + + CHARACTER(LEN=*), PARAMETER :: routineN = 'H_KS_spinor_kp' + + INTEGER :: handle, i_dim, ikp, n_spin, & + nkp_bs_and_DOS, s + INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index_scf + REAL(KIND=dp), DIMENSION(3) :: xkp + TYPE(cp_cfm_type) :: cfm_V_SOC_xyz + TYPE(cp_fm_struct_type), POINTER :: str + TYPE(kpoint_type), POINTER :: kpoints_scf + TYPE(neighbor_list_set_p_type), DIMENSION(:), & + POINTER :: sab_nl + + CALL timeset(routineN, handle) + + nkp_bs_and_DOS = bs_env%nkp_bs_and_DOS + n_spin = bs_env%n_spin + s = bs_env%n_ao + 1 + str => bs_env%cfm_ks_kp(1, 1)%matrix_struct + + CALL cp_cfm_create(cfm_V_SOC_xyz, bs_env%cfm_work_mo%matrix_struct) + + CALL alloc_cfm_double_array_1d(bs_env%cfm_SOC_spinor_ao, bs_env%cfm_ks_kp(1, 1), nkp_bs_and_DOS) + + CALL get_qs_env(qs_env, kpoints=kpoints_scf) + + NULLIFY (sab_nl) + CALL get_kpoint_info(kpoints_scf, sab_nl=sab_nl, cell_to_index=cell_to_index_scf) + + DO ikp = 1, nkp_bs_and_DOS + + xkp(1:3) = bs_env%kpoints_DOS%xkp(1:3, ikp) + + DO i_dim = 1, 3 + + CALL cp_cfm_set_all(cfm_V_SOC_xyz, z_zero) + CALL rsmat_to_kp(bs_env%mat_V_SOC_xyz, i_dim, xkp, cell_to_index_scf, & + sab_nl, bs_env, cfm_V_SOC_xyz, imag_rs_mat=.TRUE.) + + ! multiply V_SOC with i because bs_env%mat_V_SOC_xyz stores imag. part (real part = 0) + CALL cp_cfm_scale(gaussi, cfm_V_SOC_xyz) + + SELECT CASE (i_dim) + CASE (1) + ! add V^SOC_x * σ_x for σ_x = ( (0,1) (1,0) ) + CALL add_cfm_submat(bs_env%cfm_SOC_spinor_ao(ikp), cfm_V_SOC_xyz, 1, s) + CALL add_cfm_submat(bs_env%cfm_SOC_spinor_ao(ikp), cfm_V_SOC_xyz, s, 1) + CASE (2) + ! add V^SOC_y * σ_y for σ_y = ( (0,-i) (i,0) ) + CALL cp_cfm_scale(gaussi, cfm_V_SOC_xyz) + CALL add_cfm_submat(bs_env%cfm_SOC_spinor_ao(ikp), cfm_V_SOC_xyz, 1, s) + CALL cp_cfm_scale(-z_one, cfm_V_SOC_xyz) + CALL add_cfm_submat(bs_env%cfm_SOC_spinor_ao(ikp), cfm_V_SOC_xyz, s, 1) + CASE (3) + ! add V^SOC_z * σ_z for σ_z = ( (1,0) (0,1) ) + CALL add_cfm_submat(bs_env%cfm_SOC_spinor_ao(ikp), cfm_V_SOC_xyz, 1, 1) + CALL cp_cfm_scale(-z_one, cfm_V_SOC_xyz) + CALL add_cfm_submat(bs_env%cfm_SOC_spinor_ao(ikp), cfm_V_SOC_xyz, s, s) + END SELECT + + END DO + + END DO ! ikp + + CALL cp_cfm_release(cfm_V_SOC_xyz) + + CALL timestop(handle) + + END SUBROUTINE H_KS_spinor_kp + +! ************************************************************************************************** +!> \brief ... +!> \param cfm_array ... +!> \param cfm_template ... +!> \param n ... +! ************************************************************************************************** + SUBROUTINE alloc_cfm_double_array_1d(cfm_array, cfm_template, n) + TYPE(cp_cfm_type), ALLOCATABLE, DIMENSION(:) :: cfm_array + TYPE(cp_cfm_type) :: cfm_template + INTEGER :: n + + CHARACTER(LEN=*), PARAMETER :: routineN = 'alloc_cfm_double_array_1d' + + INTEGER :: handle, i + + CALL timeset(routineN, handle) + + ALLOCATE (cfm_array(n)) + DO i = 1, n + CALL create_cfm_double(cfm_array(i), cfm_orig=cfm_template) + CALL cp_cfm_set_all(cfm_array(i), z_zero) + END DO + + CALL timestop(handle) + + END SUBROUTINE alloc_cfm_double_array_1d + +! ************************************************************************************************** +!> \brief ... +!> \param bs_env ... +! ************************************************************************************************** + SUBROUTINE get_all_VBM_CBM_bandgaps(bs_env) + + TYPE(post_scf_bandstructure_type), POINTER :: bs_env + + CHARACTER(LEN=*), PARAMETER :: routineN = 'get_all_VBM_CBM_bandgaps' + + INTEGER :: handle + + CALL timeset(routineN, handle) + + CALL get_VBM_CBM_bandgaps(bs_env%band_edges_scf, bs_env%eigenval_scf, bs_env) + CALL get_VBM_CBM_bandgaps(bs_env%band_edges_G0W0, bs_env%eigenval_G0W0, bs_env) + CALL get_VBM_CBM_bandgaps(bs_env%band_edges_HF, bs_env%eigenval_HF, bs_env) + + CALL timestop(handle) + + END SUBROUTINE get_all_VBM_CBM_bandgaps + +! ************************************************************************************************** +!> \brief ... +!> \param band_edges ... +!> \param ev ... +!> \param bs_env ... +! ************************************************************************************************** + SUBROUTINE get_VBM_CBM_bandgaps(band_edges, ev, bs_env) + TYPE(band_edges_type) :: band_edges + REAL(KIND=dp), DIMENSION(:, :, :) :: ev + TYPE(post_scf_bandstructure_type), POINTER :: bs_env + + CHARACTER(LEN=*), PARAMETER :: routineN = 'get_VBM_CBM_bandgaps' + + INTEGER :: handle, homo, homo_1, homo_2, ikp, & + ispin, lumo, lumo_1, lumo_2, n_mo + REAL(KIND=dp) :: E_DBG_at_ikp + + CALL timeset(routineN, handle) + + n_mo = bs_env%n_ao + + band_edges%DBG = 1000.0_dp + + SELECT CASE (bs_env%n_spin) + CASE (1) + homo = bs_env%n_occ(1) + lumo = homo + 1 + band_edges%VBM = MAXVAL(ev(1:homo, :, 1)) + band_edges%CBM = MINVAL(ev(homo + 1:n_mo, :, 1)) + CASE (2) + homo_1 = bs_env%n_occ(1) + lumo_1 = homo_1 + 1 + homo_2 = bs_env%n_occ(2) + lumo_2 = homo_2 + 1 + band_edges%VBM = MAX(MAXVAL(ev(1:homo_1, :, 1)), MAXVAL(ev(1:homo_2, :, 2))) + band_edges%CBM = MIN(MINVAL(ev(homo_1 + 1:n_mo, :, 1)), MINVAL(ev(homo_2 + 1:n_mo, :, 2))) + CASE DEFAULT + CPABORT("Error with number of spins.") + END SELECT + + band_edges%IDBG = band_edges%CBM - band_edges%VBM + + DO ispin = 1, bs_env%n_spin + + homo = bs_env%n_occ(ispin) + + DO ikp = 1, bs_env%nkp_bs_and_DOS + + E_DBG_at_ikp = -MAXVAL(ev(1:homo, ikp, ispin)) + MINVAL(ev(homo + 1:n_mo, ikp, ispin)) + + IF (E_DBG_at_ikp < band_edges%DBG) band_edges%DBG = E_DBG_at_ikp + + END DO + + END DO + + CALL timestop(handle) + + END SUBROUTINE get_VBM_CBM_bandgaps + END MODULE post_scf_bandstructure_utils diff --git a/src/rpa_gw.F b/src/rpa_gw.F index b6df1323d5..b638b0275e 100644 --- a/src/rpa_gw.F +++ b/src/rpa_gw.F @@ -897,8 +897,10 @@ SUBROUTINE compute_GW_self_energy(vec_Sigma_c_gw, dimen_nm_gw, dimen_RI, gw_corr DO ispin = 1, nspins CALL compute_GW_self_energy_deep(vec_Sigma_c_gw(:, :, :, ispin), dimen_nm_gw, dimen_RI, & - gw_corr_lev_occ(ispin), homo(ispin), jquad, nmo, & - num_fit_points, do_periodic, fermi_level_offset, omega, Eigenval(:, ispin), delta_corr, & + gw_corr_lev_occ(ispin), gw_corr_lev_virt(ispin), & + homo(ispin), jquad, nmo, & + num_fit_points, do_periodic, fermi_level_offset, omega, & + Eigenval(:, ispin), delta_corr, & vec_omega_fit_gw, vec_W_gw(:, ispin), wj, fm_mat_Q, & fm_mat_S_gw(ispin), fm_mat_S_gw_work(ispin)) END DO @@ -1037,6 +1039,7 @@ END SUBROUTINE compute_W_cubic_GW !> \param dimen_nm_gw ... !> \param dimen_RI ... !> \param gw_corr_lev_occ ... +!> \param gw_corr_lev_virt ... !> \param homo ... !> \param jquad ... !> \param nmo ... @@ -1053,14 +1056,18 @@ END SUBROUTINE compute_W_cubic_GW !> \param fm_mat_S_gw ... !> \param fm_mat_S_gw_work ... ! ************************************************************************************************** - SUBROUTINE compute_GW_self_energy_deep(vec_Sigma_c_gw, dimen_nm_gw, dimen_RI, gw_corr_lev_occ, homo, jquad, nmo, num_fit_points, & - do_periodic, fermi_level_offset, omega, Eigenval, delta_corr, vec_omega_fit_gw, vec_W_gw, & + SUBROUTINE compute_GW_self_energy_deep(vec_Sigma_c_gw, dimen_nm_gw, dimen_RI, & + gw_corr_lev_occ, gw_corr_lev_virt, & + homo, jquad, nmo, num_fit_points, & + do_periodic, fermi_level_offset, omega, Eigenval, & + delta_corr, vec_omega_fit_gw, vec_W_gw, & wj, fm_mat_Q, fm_mat_S_gw, fm_mat_S_gw_work) COMPLEX(KIND=dp), DIMENSION(:, :, :), & INTENT(INOUT) :: vec_Sigma_c_gw INTEGER, INTENT(IN) :: dimen_nm_gw, dimen_RI, gw_corr_lev_occ, & - homo, jquad, nmo, num_fit_points + gw_corr_lev_virt, homo, jquad, nmo, & + num_fit_points LOGICAL, INTENT(IN) :: do_periodic REAL(KIND=dp), INTENT(IN) :: fermi_level_offset REAL(KIND=dp), INTENT(INOUT) :: omega @@ -1119,9 +1126,9 @@ SUBROUTINE compute_GW_self_energy_deep(vec_Sigma_c_gw, dimen_nm_gw, dimen_RI, gw ! set the Fermi energy for occ orbitals slightly above the HOMO and ! for virt orbitals slightly below the LUMO IF (n_global <= homo) THEN - e_fermi = Eigenval(homo) + fermi_level_offset + e_fermi = MAXVAL(Eigenval(homo - gw_corr_lev_occ + 1:homo)) + fermi_level_offset ELSE - e_fermi = Eigenval(homo + 1) - fermi_level_offset + e_fermi = MINVAL(Eigenval(homo + 1:homo + gw_corr_lev_virt)) - fermi_level_offset END IF ! add here the periodic correction @@ -1420,7 +1427,7 @@ SUBROUTINE compute_QP_energies(vec_Sigma_c_gw, count_ev_sc_GW, gw_corr_lev_occ, z_value(:, ikp, 1), m_value(:, ikp, 1), vec_Sigma_c_gw(:, :, ikp, 1), & mp2_env%ri_g0w0%vec_Sigma_x_minus_vxc_gw(:, 1, ikp), & Eigenval(:, ikp, 1), Eigenval_scf(:, ikp, 1), n_level_gw, & - gw_corr_lev_occ(1), num_poles, & + gw_corr_lev_occ(1), gw_corr_lev_virt(1), num_poles, & num_fit_points, crossing_search, homo(1), stop_crit, & fermi_level_offset, do_im_time) @@ -1429,7 +1436,7 @@ SUBROUTINE compute_QP_energies(vec_Sigma_c_gw, count_ev_sc_GW, gw_corr_lev_occ, z_value(:, ikp, 1), m_value(:, ikp, 1), vec_Sigma_c_gw(:, :, ikp, 1), & mp2_env%ri_g0w0%vec_Sigma_x_minus_vxc_gw(:, 1, ikp), & Eigenval(:, ikp, 1), Eigenval_scf(:, ikp, 1), n_level_gw, & - gw_corr_lev_occ(1), mp2_env%ri_g0w0%nparam_pade, & + gw_corr_lev_occ(1), gw_corr_lev_virt(1), mp2_env%ri_g0w0%nparam_pade, & num_fit_points, crossing_search, homo(1), fermi_level_offset, & do_im_time, mp2_env%ri_g0w0%print_self_energy, count_ev_sc_GW, & vec_gw_dos, dos_lower_bound, dos_precision, ndos, & @@ -1449,7 +1456,7 @@ SUBROUTINE compute_QP_energies(vec_Sigma_c_gw, count_ev_sc_GW, gw_corr_lev_occ, z_value(:, ikp, 2), m_value(:, ikp, 2), vec_Sigma_c_gw(:, :, ikp, 2), & mp2_env%ri_g0w0%vec_Sigma_x_minus_vxc_gw(:, 2, ikp), & Eigenval(:, ikp, 2), Eigenval_scf(:, ikp, 2), n_level_gw, & - gw_corr_lev_occ(2), num_poles, & + gw_corr_lev_occ(2), gw_corr_lev_virt(2), num_poles, & num_fit_points, crossing_search, homo(2), stop_crit, & fermi_level_offset, do_im_time) CASE (gw_pade_approx) @@ -1457,7 +1464,7 @@ SUBROUTINE compute_QP_energies(vec_Sigma_c_gw, count_ev_sc_GW, gw_corr_lev_occ, z_value(:, ikp, 2), m_value(:, ikp, 2), vec_Sigma_c_gw(:, :, ikp, 2), & mp2_env%ri_g0w0%vec_Sigma_x_minus_vxc_gw(:, 2, ikp), & Eigenval(:, ikp, 2), Eigenval_scf(:, ikp, 2), n_level_gw, & - gw_corr_lev_occ(2), mp2_env%ri_g0w0%nparam_pade, & + gw_corr_lev_occ(2), gw_corr_lev_virt(2), mp2_env%ri_g0w0%nparam_pade, & num_fit_points, crossing_search, homo(2), & fermi_level_offset, do_im_time, & mp2_env%ri_g0w0%print_self_energy, count_ev_sc_GW, & @@ -3856,6 +3863,7 @@ END SUBROUTINE average_degenerate_levels !> \param Eigenval_scf ... !> \param n_level_gw ... !> \param gw_corr_lev_occ ... +!> \param gw_corr_lev_vir ... !> \param num_poles ... !> \param num_fit_points ... !> \param crossing_search ... @@ -3866,7 +3874,8 @@ END SUBROUTINE average_degenerate_levels ! ************************************************************************************************** SUBROUTINE fit_and_continuation_2pole(vec_gw_energ, vec_omega_fit_gw, & z_value, m_value, vec_Sigma_c_gw, vec_Sigma_x_minus_vxc_gw, & - Eigenval, Eigenval_scf, n_level_gw, gw_corr_lev_occ, num_poles, & + Eigenval, Eigenval_scf, n_level_gw, & + gw_corr_lev_occ, gw_corr_lev_vir, num_poles, & num_fit_points, crossing_search, homo, stop_crit, & fermi_level_offset, do_gw_im_time) @@ -3875,7 +3884,8 @@ SUBROUTINE fit_and_continuation_2pole(vec_gw_energ, vec_omega_fit_gw, & COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(IN) :: vec_Sigma_c_gw REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: vec_Sigma_x_minus_vxc_gw, Eigenval, & Eigenval_scf - INTEGER, INTENT(IN) :: n_level_gw, gw_corr_lev_occ, num_poles, & + INTEGER, INTENT(IN) :: n_level_gw, gw_corr_lev_occ, & + gw_corr_lev_vir, num_poles, & num_fit_points, crossing_search, homo REAL(KIND=dp), INTENT(IN) :: stop_crit, fermi_level_offset LOGICAL, INTENT(IN) :: do_gw_im_time @@ -4120,9 +4130,9 @@ SUBROUTINE fit_and_continuation_2pole(vec_gw_energ, vec_omega_fit_gw, & ! in case of O(N^4) GW, we have the Fermi level differently for occ and virt states, see ! Fig. 1 in JCTC 12, 3623-3635 (2016) IF (n_level_gw <= gw_corr_lev_occ) THEN - e_fermi = Eigenval(homo) + fermi_level_offset + e_fermi = MAXVAL(Eigenval(homo - gw_corr_lev_occ + 1:homo)) + fermi_level_offset ELSE - e_fermi = Eigenval(homo + 1) - fermi_level_offset + e_fermi = MINVAL(Eigenval(homo + 1:homo + gw_corr_lev_vir)) - fermi_level_offset END IF END IF @@ -4271,6 +4281,7 @@ END SUBROUTINE fit_and_continuation_2pole !> \param Eigenval_scf KS/HF eigenvalue !> \param n_level_gw ... !> \param gw_corr_lev_occ ... +!> \param gw_corr_lev_vir ... !> \param nparam_pade number of pade parameters !> \param num_fit_points number of fit points for Sigma_c(iomega) !> \param crossing_search type ofr cross search to find quasiparticle energies @@ -4292,8 +4303,8 @@ END SUBROUTINE fit_and_continuation_2pole ! ************************************************************************************************** SUBROUTINE continuation_pade(vec_gw_energ, vec_omega_fit_gw, & z_value, m_value, vec_Sigma_c_gw, vec_Sigma_x_minus_vxc_gw, & - Eigenval, Eigenval_scf, n_level_gw, gw_corr_lev_occ, nparam_pade, & - num_fit_points, crossing_search, homo, & + Eigenval, Eigenval_scf, n_level_gw, gw_corr_lev_occ, gw_corr_lev_vir, & + nparam_pade, num_fit_points, crossing_search, homo, & fermi_level_offset, do_gw_im_time, print_self_energy, count_ev_sc_GW, & vec_gw_dos, dos_lower_bound, dos_precision, ndos, & min_level_self_energy, max_level_self_energy, & @@ -4307,8 +4318,8 @@ SUBROUTINE continuation_pade(vec_gw_energ, vec_omega_fit_gw, & REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: vec_Sigma_x_minus_vxc_gw, Eigenval, & Eigenval_scf INTEGER, INTENT(IN) :: n_level_gw, gw_corr_lev_occ, & - nparam_pade, num_fit_points, & - crossing_search, homo + gw_corr_lev_vir, nparam_pade, & + num_fit_points, crossing_search, homo REAL(KIND=dp), INTENT(IN) :: fermi_level_offset LOGICAL, INTENT(IN) :: do_gw_im_time, print_self_energy INTEGER, INTENT(IN) :: count_ev_sc_GW @@ -4358,9 +4369,9 @@ SUBROUTINE continuation_pade(vec_gw_energ, vec_omega_fit_gw, & ! in case of O(N^4) GW, we have the Fermi level differently for occ and virt states, see ! Fig. 1 in JCTC 12, 3623-3635 (2016) IF (n_level_gw <= gw_corr_lev_occ) THEN - e_fermi = Eigenval(homo) + fermi_level_offset + e_fermi = MAXVAL(Eigenval(homo - gw_corr_lev_occ + 1:homo)) + fermi_level_offset ELSE - e_fermi = Eigenval(homo + 1) - fermi_level_offset + e_fermi = MINVAL(Eigenval(homo + 1:homo + gw_corr_lev_vir)) - fermi_level_offset END IF END IF diff --git a/src/soc_pseudopotential_methods.F b/src/soc_pseudopotential_methods.F index 009ad5aac9..2986c2c5d7 100644 --- a/src/soc_pseudopotential_methods.F +++ b/src/soc_pseudopotential_methods.F @@ -8,10 +8,7 @@ MODULE soc_pseudopotential_methods USE atomic_kind_types, ONLY: atomic_kind_type USE core_ppnl, ONLY: build_core_ppnl - USE cp_cfm_types, ONLY: cp_cfm_create,& - cp_cfm_get_info,& - cp_cfm_set_all,& - cp_cfm_to_cfm,& + USE cp_cfm_types, ONLY: cp_cfm_get_info,& cp_cfm_type USE cp_control_types, ONLY: dft_control_type USE cp_dbcsr_api, ONLY: dbcsr_add,& @@ -26,7 +23,6 @@ MODULE soc_pseudopotential_methods copy_fm_to_dbcsr,& dbcsr_allocate_matrix_set,& dbcsr_deallocate_matrix_set - USE cp_fm_struct, ONLY: cp_fm_struct_type USE cp_fm_types, ONLY: cp_fm_create,& cp_fm_get_info,& cp_fm_release,& @@ -34,9 +30,6 @@ MODULE soc_pseudopotential_methods USE kinds, ONLY: dp USE kpoint_types, ONLY: get_kpoint_info,& kpoint_type - USE mathconstants, ONLY: gaussi,& - z_one,& - z_zero USE parallel_gemm_api, ONLY: parallel_gemm USE particle_types, ONLY: particle_type USE qs_environment_types, ONLY: get_qs_env,& @@ -45,9 +38,6 @@ MODULE soc_pseudopotential_methods USE qs_kind_types, ONLY: qs_kind_type USE qs_neighbor_list_types, ONLY: get_neighbor_list_set_p,& neighbor_list_set_p_type - USE soc_pseudopotential_utils, ONLY: add_dbcsr_submat,& - add_fm_submat,& - create_cfm_double USE virial_types, ONLY: virial_type #include "./base/base_uses.f90" @@ -57,14 +47,15 @@ MODULE soc_pseudopotential_methods CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'soc_pseudopotential_methods' - PUBLIC :: V_SOC_xyz_from_pseudopotential, H_KS_spinor, remove_soc_outside_energy_window_ao, & + PUBLIC :: V_SOC_xyz_from_pseudopotential, & + remove_soc_outside_energy_window_ao, & remove_soc_outside_energy_window_mo CONTAINS ! ************************************************************************************************** -!> \brief 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) +!> \brief V^SOC_µν^(α),R = ħ/2 < ϕ_µ cell O | sum_ℓ ΔV_ℓ^SO(r,r') L^(α) | ϕ_ν cell R>, α = x,y,z +!> see Hartwigsen, Goedecker, Hutter, Eq.(18), (19) (doi.org/10.1103/PhysRevB.58.3641) !> Caution: V^SOC_µν^(α) is purely imaginary and Hermitian; V^SOC_µν^(α) is stored as real !> dbcsr matrix mat_V_SOC_xyz without symmetry; V^SOC_µν^(α) is stored without !> the imaginary unit, i.e. mat_V_SOC_xyz is real and antisymmetric @@ -72,6 +63,7 @@ MODULE soc_pseudopotential_methods !> \param mat_V_SOC_xyz ... !> \par History !> * 09.2023 created +!> \author Jan Wilhelm ! ************************************************************************************************** SUBROUTINE V_SOC_xyz_from_pseudopotential(qs_env, mat_V_SOC_xyz) TYPE(qs_environment_type), POINTER :: qs_env @@ -143,16 +135,16 @@ SUBROUTINE V_SOC_xyz_from_pseudopotential(qs_env, mat_V_SOC_xyz) eps_ppnl, nimages=nimages, cell_to_index=cell_to_index, & basis_type="ORB", matrix_l=mat_l) -! JW TODO: 1 -> nimages and TEST THIS!!! - NULLIFY (mat_l_nosym) CALL dbcsr_allocate_matrix_set(mat_l_nosym, 3, nimages) DO xyz = 1, 3 DO img = 1, nimages + ALLOCATE (mat_l_nosym(xyz, img)%matrix) CALL dbcsr_create(mat_l_nosym(xyz, img)%matrix, template=matrix_s(1, 1)%matrix, & matrix_type=dbcsr_type_no_symmetry) CALL dbcsr_desymmetrize(mat_l(xyz, img)%matrix, mat_l_nosym(xyz, img)%matrix) + END DO END DO @@ -178,77 +170,6 @@ SUBROUTINE V_SOC_xyz_from_pseudopotential(qs_env, mat_V_SOC_xyz) END SUBROUTINE V_SOC_xyz_from_pseudopotential -! ************************************************************************************************** -!> \brief Spinor KS-matrix H_µν,σσ' = h_µν,σ*δ_σσ' + sum_α V^SOC_µν^(α)*Pauli-matrix^(α)_σσ', see -!> Hartwigsen, Goedecker, Hutter, Eq.(18) (doi.org/10.1103/PhysRevB.58.3641) -!> \param cfm_ks_spinor_ao ... -!> \param fm_ks ... -!> \param n_spin ... -!> \param mat_V_SOC_xyz ... -!> \param cfm_s_double ... -!> \param fm_s ... -!> \param cfm_SOC_spinor_ao ... -! ************************************************************************************************** - SUBROUTINE H_KS_spinor(cfm_ks_spinor_ao, fm_ks, n_spin, mat_V_SOC_xyz, cfm_s_double, fm_s, & - cfm_SOC_spinor_ao) - TYPE(cp_cfm_type) :: cfm_ks_spinor_ao - TYPE(cp_fm_type), DIMENSION(:) :: fm_ks - INTEGER :: n_spin - TYPE(dbcsr_p_type), DIMENSION(:) :: mat_V_SOC_xyz - TYPE(cp_cfm_type), OPTIONAL :: cfm_s_double - TYPE(cp_fm_type), OPTIONAL :: fm_s - TYPE(cp_cfm_type), OPTIONAL :: cfm_SOC_spinor_ao - - CHARACTER(LEN=*), PARAMETER :: routineN = 'H_KS_spinor' - - INTEGER :: handle, nao, s - TYPE(cp_fm_struct_type), POINTER :: str - - CALL timeset(routineN, handle) - - CALL cp_fm_get_info(fm_ks(1), nrow_global=nao) - - CALL create_cfm_double(fm_ks(1), cfm_ks_spinor_ao) - CALL cp_cfm_set_all(cfm_ks_spinor_ao, z_zero) - - str => fm_ks(1)%matrix_struct - - s = nao + 1 - - CALL add_dbcsr_submat(cfm_ks_spinor_ao, mat_V_SOC_xyz(1)%matrix, str, s, 1, z_one, .TRUE.) - CALL add_dbcsr_submat(cfm_ks_spinor_ao, mat_V_SOC_xyz(2)%matrix, str, s, 1, gaussi, .TRUE.) - CALL add_dbcsr_submat(cfm_ks_spinor_ao, mat_V_SOC_xyz(3)%matrix, str, 1, 1, z_one, .FALSE.) - CALL add_dbcsr_submat(cfm_ks_spinor_ao, mat_V_SOC_xyz(3)%matrix, str, s, s, -z_one, .FALSE.) - - IF (PRESENT(cfm_SOC_spinor_ao)) THEN - CALL cp_cfm_create(cfm_SOC_spinor_ao, cfm_ks_spinor_ao%matrix_struct) - CALL cp_cfm_to_cfm(cfm_ks_spinor_ao, cfm_SOC_spinor_ao) - END IF - - CALL add_fm_submat(cfm_ks_spinor_ao, fm_ks(1), 1, 1) - - SELECT CASE (n_spin) - CASE (1) - CALL add_fm_submat(cfm_ks_spinor_ao, fm_ks(1), s, s) - CASE (2) - CPASSERT(SIZE(fm_ks) == 2) - CALL add_fm_submat(cfm_ks_spinor_ao, fm_ks(2), s, s) - CASE DEFAULT - CPABORT("We have either one or two spin channels.") - END SELECT - - IF (PRESENT(cfm_s_double)) THEN - CPASSERT(PRESENT(fm_s)) - CALL create_cfm_double(fm_s, cfm_s_double) - CALL cp_cfm_set_all(cfm_s_double, z_zero) - CALL add_fm_submat(cfm_s_double, fm_s, 1, 1) - CALL add_fm_submat(cfm_s_double, fm_s, s, s) - END IF - - CALL timestop(handle) - - END SUBROUTINE H_KS_spinor - ! ************************************************************************************************** !> \brief Remove SOC outside of energy window (otherwise, numerical problems arise !> because energetically low semicore states and energetically very high diff --git a/src/soc_pseudopotential_utils.F b/src/soc_pseudopotential_utils.F index f900c7cc90..612d98a23d 100644 --- a/src/soc_pseudopotential_utils.F +++ b/src/soc_pseudopotential_utils.F @@ -40,8 +40,8 @@ MODULE soc_pseudopotential_utils CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'soc_pseudopotential_utils' - PUBLIC :: add_dbcsr_submat, cfm_add_on_diag, add_fm_submat, get_cfm_submat, create_cfm_double, & - add_cfm_submat + PUBLIC :: add_dbcsr_submat, cfm_add_on_diag, add_fm_submat, add_cfm_submat, & + get_cfm_submat, create_cfm_double CONTAINS @@ -91,7 +91,8 @@ SUBROUTINE add_dbcsr_submat(cfm_mat_target, mat_source, fm_struct_source, & nrow=nao, ncol=nao, & s_firstrow=1, s_firstcol=1, & t_firstrow=nstart_row, t_firstcol=nstart_col) - + ! careful: inside add_dbcsr_submat, mat_V_SOC_xyz is multiplied by i because the real matrix + ! mat_V_SOC_xyz is antisymmetric as V_SOC matrix is purely imaginary and Hermitian CALL cp_cfm_scale_and_add_fm(z_zero, cfm_mat_work_double, gaussi, fm_mat_work_double_im) CALL cp_cfm_scale(factor, cfm_mat_work_double) @@ -144,9 +145,11 @@ SUBROUTINE cfm_add_on_diag(cfm, alpha) i_global = row_indices(i_row) IF (j_global == i_global) THEN IF (i_global .LE. nao) THEN - cfm%local_data(i_row, j_col) = cfm%local_data(i_row, j_col) + alpha(i_global)*z_one + cfm%local_data(i_row, j_col) = cfm%local_data(i_row, j_col) + & + alpha(i_global)*z_one ELSE - cfm%local_data(i_row, j_col) = cfm%local_data(i_row, j_col) + alpha(i_global - nao)*z_one + cfm%local_data(i_row, j_col) = cfm%local_data(i_row, j_col) + & + alpha(i_global - nao)*z_one END IF END IF END DO @@ -314,31 +317,51 @@ END SUBROUTINE get_cfm_submat ! ************************************************************************************************** !> \brief ... -!> \param fm_orig ... !> \param cfm_double ... +!> \param fm_orig ... +!> \param cfm_orig ... ! ************************************************************************************************** - SUBROUTINE create_cfm_double(fm_orig, cfm_double) - TYPE(cp_fm_type) :: fm_orig + SUBROUTINE create_cfm_double(cfm_double, fm_orig, cfm_orig) TYPE(cp_cfm_type) :: cfm_double + TYPE(cp_fm_type), OPTIONAL :: fm_orig + TYPE(cp_cfm_type), OPTIONAL :: cfm_orig CHARACTER(LEN=*), PARAMETER :: routineN = 'create_cfm_double' INTEGER :: handle, ncol_global_orig, & nrow_global_orig - TYPE(cp_fm_struct_type), POINTER :: fm_struct_double + LOGICAL :: do_cfm_templ, do_fm_templ + TYPE(cp_fm_struct_type), POINTER :: matrix_struct, matrix_struct_double CALL timeset(routineN, handle) - CALL cp_fm_get_info(matrix=fm_orig, nrow_global=nrow_global_orig, ncol_global=ncol_global_orig) + do_fm_templ = PRESENT(fm_orig) + do_cfm_templ = PRESENT(cfm_orig) + + ! either fm template or cfm template + CPASSERT(do_fm_templ .NEQV. do_cfm_templ) - CALL cp_fm_struct_create(fm_struct_double, & + IF (do_fm_templ) THEN + CALL cp_fm_get_info(matrix=fm_orig, nrow_global=nrow_global_orig, & + ncol_global=ncol_global_orig) + matrix_struct => fm_orig%matrix_struct + END IF + IF (do_cfm_templ) THEN + CALL cp_cfm_get_info(matrix=cfm_orig, nrow_global=nrow_global_orig, & + ncol_global=ncol_global_orig) + matrix_struct => cfm_orig%matrix_struct + END IF + + CALL cp_fm_struct_create(matrix_struct_double, & nrow_global=2*nrow_global_orig, & ncol_global=2*ncol_global_orig, & - template_fmstruct=fm_orig%matrix_struct) + template_fmstruct=matrix_struct) + + CALL cp_cfm_create(cfm_double, matrix_struct_double) - CALL cp_cfm_create(cfm_double, fm_struct_double) + CALL cp_cfm_set_all(cfm_double, z_zero) - CALL cp_fm_struct_release(fm_struct_double) + CALL cp_fm_struct_release(matrix_struct_double) CALL timestop(handle) diff --git a/tests/QS/regtest-scalable-gw/TEST_FILES b/tests/QS/regtest-scalable-gw/TEST_FILES index a605c1fbdf..d5f9337dbd 100644 --- a/tests/QS/regtest-scalable-gw/TEST_FILES +++ b/tests/QS/regtest-scalable-gw/TEST_FILES @@ -2,6 +2,7 @@ 02_G0W0_IH_SOC_LDOS.inp 106 1e-03 17.826 02_G0W0_IH_SOC_LDOS.inp 107 1e-03 11.344 03_G0W0_bandstructure_IH_chain.inp 106 1e-03 21.977 +03_G0W0_bandstructure_IH_chain.inp 108 1e-03 21.692 04_G0W0_SOC_TeH_chain_open_shell_kp_extrapol.inp 106 1e-03 5.174 05_G0W0_SOC_TeH_chain_open_shell.inp 106 1e-03 5.178 05_G0W0_SOC_TeH_chain_open_shell.inp 107 1e-03 0.306