From bf6c9a47c984b7c8bb1cd90d6dbd08098a8232e6 Mon Sep 17 00:00:00 2001 From: Jan Wilhelm Date: Wed, 27 Dec 2023 12:21:00 +0200 Subject: [PATCH] =?UTF-8?q?Print=20LDOS=20in=20units=20of=201/(eV*=C3=85^2?= =?UTF-8?q?)=20(LDOS=20is=20integrated=20in=20z=20direction)?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- src/gw_methods.F | 139 ++----------- src/post_scf_bandstructure_types.F | 9 +- src/post_scf_bandstructure_utils.F | 304 ++++++++++++++++++++++------- 3 files changed, 256 insertions(+), 196 deletions(-) diff --git a/src/gw_methods.F b/src/gw_methods.F index 27eee13fad..d80711d0b6 100644 --- a/src/gw_methods.F +++ b/src/gw_methods.F @@ -64,8 +64,7 @@ MODULE gw_methods dp,& int_8 USE kpoint_coulomb_2c, ONLY: build_2c_coulomb_matrix_kp - USE kpoint_types, ONLY: get_kpoint_info,& - kpoint_type + USE kpoint_types, ONLY: kpoint_type USE machine, ONLY: m_walltime USE mathconstants, ONLY: gaussi,& twopi,& @@ -78,15 +77,14 @@ MODULE gw_methods USE particle_types, ONLY: particle_type USE physcon, ONLY: evolt USE post_scf_bandstructure_types, ONLY: post_scf_bandstructure_type - USE post_scf_bandstructure_utils, ONLY: cfm_ikp_from_fm_Gamma,& + USE post_scf_bandstructure_utils, ONLY: MIC_contribution_from_ikp,& + cfm_ikp_from_fm_Gamma,& get_fname USE qs_environment_types, ONLY: get_qs_env,& qs_environment_type USE qs_kind_types, ONLY: qs_kind_type USE qs_tensors, ONLY: build_3c_integrals USE rpa_gw, ONLY: continuation_pade - USE rpa_gw_im_time_util, ONLY: compute_weight_re_im,& - get_atom_index_from_basis_function_index USE rpa_gw_kpoints_util, ONLY: cp_cfm_power,& cp_cfm_upper_to_full #include "./base/base_uses.f90" @@ -1187,10 +1185,10 @@ CONTAINS CALL compute_cfm_W_ikp_freq_j(bs_env, cfm_chi_ikp_freq_j, cfm_V_sqrt_ikp, & cfm_M_inv_V_sqrt_ikp, cfm_W_ikp_freq_j) - ! 5. k-point integration W_PQ(iω_j,k_i) to W_PQ^MIC(iω_j) + ! 5. k-point integration W_PQ(iω_j, k_i) to W_PQ^MIC(iω_j) IF (.NOT. bs_env%approx_kp_extrapol) THEN - CALL compute_fm_W_MIC_freq_j_ikp_contribution(bs_env, qs_env, fm_W_MIC_freq_j, & - cfm_W_ikp_freq_j, ikp) + CALL MIC_contribution_from_ikp(bs_env, qs_env, fm_W_MIC_freq_j, cfm_W_ikp_freq_j, ikp, & + bs_env%kpoints_chi_eps_W, "RI_AUX") END IF ! 6. for approximate kpoint extrapolation: get W_PQ^MIC(iω_1) with and without k-point @@ -1202,24 +1200,23 @@ CONTAINS IF (j_w == 1) THEN ! k-point extrapolated - CALL compute_fm_W_MIC_freq_j_ikp_contribution(bs_env, qs_env, & - bs_env%fm_W_MIC_freq_1_extra, & - cfm_W_ikp_freq_j, ikp) + CALL MIC_contribution_from_ikp(bs_env, qs_env, bs_env%fm_W_MIC_freq_1_extra, & + cfm_W_ikp_freq_j, ikp, bs_env%kpoints_chi_eps_W, & + "RI_AUX") ! non-kpoint extrapolated IF (ikp .LE. bs_env%nkp_chi_eps_W_orig) THEN - CALL compute_fm_W_MIC_freq_j_ikp_contribution(bs_env, qs_env, & - bs_env%fm_W_MIC_freq_1_no_extra, & - cfm_W_ikp_freq_j, ikp, & - wkp_ext=bs_env%wkp_orig) + CALL MIC_contribution_from_ikp(bs_env, qs_env, bs_env%fm_W_MIC_freq_1_no_extra, & + cfm_W_ikp_freq_j, ikp, bs_env%kpoints_chi_eps_W, & + "RI_AUX", wkp_ext=bs_env%wkp_orig) END IF END IF ! for all ω_j, we need to compute W^MIC without k-point extrpolation IF (ikp .LE. bs_env%nkp_chi_eps_W_orig) THEN - CALL compute_fm_W_MIC_freq_j_ikp_contribution(bs_env, qs_env, fm_W_MIC_freq_j, & - cfm_W_ikp_freq_j, ikp, & - wkp_ext=bs_env%wkp_orig) + CALL MIC_contribution_from_ikp(bs_env, qs_env, fm_W_MIC_freq_j, cfm_W_ikp_freq_j, & + ikp, bs_env%kpoints_chi_eps_W, "RI_AUX", & + wkp_ext=bs_env%wkp_orig) END IF END IF @@ -1549,112 +1546,6 @@ CONTAINS END SUBROUTINE cfm_add_on_diag -! ************************************************************************************************** -!> \brief ... -!> \param bs_env ... -!> \param qs_env ... -!> \param fm_W_MIC_freq_j ... -!> \param cfm_W_ikp_freq_j ... -!> \param ikp ... -!> \param wkp_ext ... -! ************************************************************************************************** - SUBROUTINE compute_fm_W_MIC_freq_j_ikp_contribution(bs_env, qs_env, fm_W_MIC_freq_j, & - cfm_W_ikp_freq_j, ikp, wkp_ext) - TYPE(post_scf_bandstructure_type), POINTER :: bs_env - TYPE(qs_environment_type), POINTER :: qs_env - TYPE(cp_fm_type) :: fm_W_MIC_freq_j - TYPE(cp_cfm_type) :: cfm_W_ikp_freq_j - INTEGER :: ikp - REAL(KIND=dp), OPTIONAL :: wkp_ext - - CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_fm_W_MIC_freq_j_ikp_contribution' - - INTEGER :: handle, i_RI, iatom, iatom_old, irow, & - j_RI, jatom, jatom_old, jcol, n_RI, & - ncol_local, nrow_local, num_cells - INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_from_RI_index - INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices - INTEGER, DIMENSION(:, :), POINTER :: index_to_cell - REAL(KIND=dp) :: contribution, weight_im, weight_re, & - wkp_of_ikp - REAL(KIND=dp), DIMENSION(3, 3) :: hmat - REAL(KIND=dp), DIMENSION(:), POINTER :: wkp - REAL(KIND=dp), DIMENSION(:, :), POINTER :: xkp - TYPE(cell_type), POINTER :: cell - TYPE(particle_type), DIMENSION(:), POINTER :: particle_set - - CALL timeset(routineN, handle) - - n_RI = bs_env%n_RI - ALLOCATE (atom_from_RI_index(n_RI)) - CALL get_atom_index_from_basis_function_index(qs_env, atom_from_RI_index, n_RI, "RI_AUX") - - NULLIFY (cell, particle_set) - CALL get_qs_env(qs_env, cell=cell, particle_set=particle_set) - CALL get_cell(cell=cell, h=hmat) - - CALL cp_cfm_get_info(matrix=cfm_W_ikp_freq_j, & - nrow_local=nrow_local, & - ncol_local=ncol_local, & - row_indices=row_indices, & - col_indices=col_indices) - - CALL get_kpoint_info(bs_env%kpoints_chi_eps_W, xkp=xkp, wkp=wkp) - index_to_cell => bs_env%kpoints_chi_eps_W%index_to_cell - num_cells = SIZE(index_to_cell, 2) - - iatom_old = 0 - jatom_old = 0 - - DO irow = 1, nrow_local - DO jcol = 1, ncol_local - - i_RI = row_indices(irow) - j_RI = col_indices(jcol) - - iatom = atom_from_RI_index(i_RI) - jatom = atom_from_RI_index(j_RI) - - IF (PRESENT(wkp_ext)) THEN - wkp_of_ikp = wkp_ext - ELSE - SELECT CASE (bs_env%l_RI(i_RI) + bs_env%l_RI(j_RI)) - CASE (0) - ! both RI functions are s-functions, k-extrapolation for 2D and 3D - wkp_of_ikp = wkp(ikp) - CASE (1) - ! one function is an s-function, the other a p-function, k-extrapolation for 3D - wkp_of_ikp = bs_env%wkp_s_p(ikp) - CASE DEFAULT - ! for any other matrix element of W, there is no need for extrapolation - wkp_of_ikp = bs_env%wkp_no_extra(ikp) - END SELECT - END IF - - IF (iatom .NE. iatom_old .OR. jatom .NE. jatom_old) THEN - - CALL compute_weight_re_im(weight_re, weight_im, & - num_cells, iatom, jatom, xkp(1:3, ikp), wkp_of_ikp, & - cell, index_to_cell, hmat, particle_set) - - iatom_old = iatom - jatom_old = jatom - - END IF - - contribution = weight_re*REAL(cfm_W_ikp_freq_j%local_data(irow, jcol)) + & - weight_im*AIMAG(cfm_W_ikp_freq_j%local_data(irow, jcol)) - - fm_W_MIC_freq_j%local_data(irow, jcol) = fm_W_MIC_freq_j%local_data(irow, jcol) & - + contribution - - END DO - END DO - - CALL timestop(handle) - - END SUBROUTINE compute_fm_W_MIC_freq_j_ikp_contribution - ! ************************************************************************************************** !> \brief ... !> \param bs_env ... diff --git a/src/post_scf_bandstructure_types.F b/src/post_scf_bandstructure_types.F index 43a4d43c6d..b539d0ec29 100644 --- a/src/post_scf_bandstructure_types.F +++ b/src/post_scf_bandstructure_types.F @@ -129,7 +129,7 @@ MODULE post_scf_bandstructure_types hinv_primitive_cell = -1.0_dp, & hmat = -1.0_dp - ! imaginary time and frequency grids + ! imaginary time and imaginary frequency grids INTEGER :: num_time_freq_points = -1, & num_freq_points_fit = -1 REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: imag_time_points, & @@ -181,6 +181,7 @@ MODULE post_scf_bandstructure_types mat_Sigma_from_Gvir_tensor, & mat_W_MIC_time_tensor + ! tensors for sparse matrix-tensor operations TYPE(dbt_type) :: t_G, & t_chi, & t_W, & @@ -246,7 +247,7 @@ MODULE post_scf_bandstructure_types TYPE(libint_potential_type) :: ri_metric, & trunc_coulomb - ! parameters for SOC calculation (SOC is always computed by default because it is cheap) + ! parameters for SOC calculation REAL(KIND=dp) :: energy_window_soc = -1.0_dp TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: mat_V_SOC_xyz => NULL() TYPE(cp_fm_type), DIMENSION(3) :: fm_V_SOC_xyz_mo @@ -255,15 +256,17 @@ MODULE post_scf_bandstructure_types cfm_s_spinor_Gamma TYPE(band_edges_type) :: band_edges_scf_SOC, & band_edges_G0W0_SOC + ! parameters for DOS and PDOS calculation REAL(KIND=dp) :: energy_window_DOS = -1.0_dp, & energy_step_DOS = -1.0_dp, & broadening_DOS = -1.0_dp - ! parameters for local VBM, CBM, gap calculation (from local density of states, ldos) + ! parameters for LDOS calculation (LDOS: local density of states) INTEGER :: int_ldos_xyz = -1 INTEGER, DIMENSION(:), POINTER :: bin_mesh => NULL() INTEGER :: n_bins_max_for_printing = -1 + REAL(KIND=dp) :: unit_ldos_int_z_inv_Ang2_eV END TYPE post_scf_bandstructure_type diff --git a/src/post_scf_bandstructure_utils.F b/src/post_scf_bandstructure_utils.F index 03c7cd0ca6..dcbfadd8fd 100644 --- a/src/post_scf_bandstructure_utils.F +++ b/src/post_scf_bandstructure_utils.F @@ -35,7 +35,6 @@ MODULE post_scf_bandstructure_utils dbcsr_deallocate_matrix_set USE cp_files, ONLY: close_file,& open_file - USE cp_fm_basic_linalg, ONLY: cp_fm_scale_and_add USE cp_fm_diag, ONLY: cp_fm_geeig_canon USE cp_fm_struct, ONLY: cp_fm_struct_create,& cp_fm_struct_release,& @@ -44,6 +43,7 @@ MODULE post_scf_bandstructure_utils cp_fm_get_diag,& cp_fm_get_info,& cp_fm_release,& + cp_fm_set_all,& cp_fm_to_fm,& cp_fm_type USE cp_log_handling, ONLY: cp_logger_get_default_io_unit @@ -62,7 +62,8 @@ MODULE post_scf_bandstructure_utils USE kinds, ONLY: default_string_length,& dp,& max_line_length - USE kpoint_types, ONLY: kpoint_create,& + USE kpoint_types, ONLY: get_kpoint_info,& + kpoint_create,& kpoint_type USE machine, ONLY: m_walltime USE mathconstants, ONLY: gaussi,& @@ -92,7 +93,8 @@ MODULE post_scf_bandstructure_utils USE qs_ks_types, ONLY: qs_ks_env_type USE qs_mo_types, ONLY: get_mo_set,& mo_set_type - USE rpa_gw_im_time_util, ONLY: get_atom_index_from_basis_function_index + 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_utils, ONLY: add_cfm_submat,& @@ -106,7 +108,7 @@ MODULE post_scf_bandstructure_utils PUBLIC :: create_and_init_bs_env, & bandstructure_primitive_cell, bandstructure_primitive_cell_spinor, & - dos_pdos_ldos, cfm_ikp_from_fm_Gamma, get_fname + dos_pdos_ldos, cfm_ikp_from_fm_Gamma, get_fname, MIC_contribution_from_ikp CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'post_scf_bandstructure_utils' @@ -1360,7 +1362,7 @@ CONTAINS band_edges_scf_SOC TYPE(cp_cfm_type) :: cfm_ks_ikp, cfm_ks_ikp_spinor, cfm_mos_ikp_spinor, cfm_s_ikp, & cfm_s_ikp_copy, cfm_s_ikp_spinor, cfm_s_ikp_spinor_copy, cfm_SOC_ikp_spinor, & - cfm_work_ikp, cfm_work_ikp_spinor + cfm_spinor_wf_ikp, cfm_work_ikp, cfm_work_ikp_spinor TYPE(cp_cfm_type), DIMENSION(2) :: cfm_mos_ikp CALL timeset(routineN, handle) @@ -1458,6 +1460,7 @@ CONTAINS 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) homo_spinor = bs_env%n_occ(1) + bs_env%n_occ(bs_env%n_spin) @@ -1515,10 +1518,12 @@ CONTAINS IF (bs_env%do_ldos) THEN CALL add_to_LDOS_2d(LDOS_scf_2d, qs_env, ikp, bs_env, cfm_mos_ikp(ispin), & eigenval(:), band_edges_scf_guess) + IF (bs_env%do_gw) THEN CALL add_to_LDOS_2d(LDOS_G0W0_2d, qs_env, ikp, bs_env, cfm_mos_ikp(ispin), & bs_env%eigenval_G0W0(:, ikp, 1), band_edges_G0W0) END IF + END IF homo = bs_env%n_occ(ispin) @@ -1532,23 +1537,24 @@ CONTAINS ! now the same with spin-orbit coupling IF (bs_env%do_soc) THEN - ! compute G0W0+SOC eigenvalues; based on these, compute band edges and DOS + ! 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) + DOS_scf_SOC, PDOS_scf_SOC, band_edges_scf_SOC, eigenval_spinor, & + cfm_spinor_wf_ikp) IF (.NOT. bs_env%do_gw) CALL write_SOC_eigenvalues(eigenval_spinor, "SCF", ikp, bs_env) IF (bs_env%do_ldos) THEN - CALL add_to_LDOS_2d(LDOS_scf_2d_SOC, qs_env, ikp, bs_env, cfm_mos_ikp_spinor, & + CALL add_to_LDOS_2d(LDOS_scf_2d_SOC, qs_env, ikp, bs_env, cfm_spinor_wf_ikp, & eigenval_spinor, band_edges_scf_guess, .TRUE., cfm_work_ikp) END IF IF (bs_env%do_gw) THEN - ! compute G0W0+SOC eigenvalues; based on these, compute band edges and DOS + ! 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) + band_edges_G0W0_SOC, eigenval_spinor_G0W0, cfm_spinor_wf_ikp) ! write SCF+SOC and G0W0+SOC eigenvalues to file ! SCF_and_G0W0_band_structure_for_kpoint__+_SOC @@ -1621,6 +1627,7 @@ CONTAINS CALL cp_cfm_release(cfm_mos_ikp_spinor) CALL cp_cfm_release(cfm_work_ikp_spinor) CALL cp_cfm_release(cfm_s_ikp_spinor_copy) + CALL cp_cfm_release(cfm_spinor_wf_ikp) CALL timestop(handle) @@ -1644,20 +1651,19 @@ CONTAINS INTEGER :: handle, i_x, i_x_bin, i_x_end, i_x_end_bin, i_x_end_glob, i_x_start, & i_x_start_bin, i_x_start_glob, i_y, i_y_bin, i_y_end, i_y_end_bin, i_y_end_glob, & i_y_start, i_y_start_bin, i_y_start_glob, n_E + INTEGER, ALLOCATABLE, DIMENSION(:, :) :: n_sum_for_bins INTEGER, DIMENSION(2) :: bin_mesh LOGICAL :: do_xy_bins - REAL(KIND=dp) :: E_min, energy_step_ldos, & - energy_window_ldos - REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: avg_LDOS_occ, avg_LDOS_vir + REAL(KIND=dp) :: E_min, energy_step, energy_window REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: LDOS_2d_bins CALL timeset(routineN, handle) n_E = SIZE(LDOS_2d, 3) - energy_window_ldos = bs_env%energy_window_DOS - energy_step_ldos = bs_env%energy_step_DOS - E_min = band_edges%VBM - 0.5_dp*energy_window_ldos + energy_window = bs_env%energy_window_DOS + energy_step = bs_env%energy_step_DOS + E_min = band_edges%VBM - 0.5_dp*energy_window bin_mesh(1:2) = bs_env%bin_mesh(1:2) do_xy_bins = (bin_mesh(1) > 0 .AND. bin_mesh(2) > 0) @@ -1679,8 +1685,6 @@ CONTAINS i_y_end_bin = i_y_end END IF - ALLOCATE (avg_LDOS_occ(i_x_start_bin:i_x_end_bin, i_y_start_bin:i_y_end_bin)) - ALLOCATE (avg_LDOS_vir(i_x_start_bin:i_x_end_bin, i_y_start_bin:i_y_end_bin)) ALLOCATE (LDOS_2d_bins(i_x_start_bin:i_x_end_bin, i_y_start_bin:i_y_end_bin, n_E)) LDOS_2d_bins(:, :, :) = 0.0_dp @@ -1696,6 +1700,9 @@ CONTAINS CALL bs_env%para_env%min(i_y_start_glob) CALL bs_env%para_env%max(i_y_end_glob) + ALLOCATE (n_sum_for_bins(bin_mesh(1), bin_mesh(2))) + n_sum_for_bins(:, :) = 0 + ! transform interval [i_x_start, i_x_end] to [1, bin_mesh(1)] (and same for y) DO i_x = i_x_start, i_x_end DO i_y = i_y_start, i_y_end @@ -1703,10 +1710,20 @@ CONTAINS i_y_bin = bin_mesh(2)*(i_y - i_y_start_glob)/(i_y_end_glob - i_y_start_glob + 1) + 1 LDOS_2d_bins(i_x_bin, i_y_bin, :) = LDOS_2d_bins(i_x_bin, i_y_bin, :) + & LDOS_2d(i_x, i_y, :) + n_sum_for_bins(i_x_bin, i_y_bin) = n_sum_for_bins(i_x_bin, i_y_bin) + 1 END DO END DO CALL bs_env%para_env%sum(LDOS_2d_bins) + CALL bs_env%para_env%sum(n_sum_for_bins) + + ! divide by number of terms in the sum so we have the average LDOS(x,y,E) + DO i_x_bin = 1, bin_mesh(1) + DO i_y_bin = 1, bin_mesh(2) + LDOS_2d_bins(i_x_bin, i_y_bin, :) = LDOS_2d_bins(i_x_bin, i_y_bin, :)/ & + REAL(n_sum_for_bins(i_x_bin, i_y_bin), KIND=dp) + END DO + END DO ELSE @@ -1716,6 +1733,8 @@ CONTAINS IF (bin_mesh(1)*bin_mesh(2) < bs_env%n_bins_max_for_printing) THEN CALL print_LDOS_2d_bins(LDOS_2d_bins, bs_env, E_min, scf_gw_soc) + ELSE + CPWARN("The number of bins for the LDOS is too large. Decrease BIN_MESH.") END IF CALL timestop(handle) @@ -1778,9 +1797,14 @@ CONTAINS CALL open_file(TRIM(fname), unit_number=iunit, file_status="REPLACE", & file_action="WRITE") + WRITE (iunit, "(2A)") " Energy E (eV) average LDOS(x,y,E) (1/(eV*Å^2), ", & + "integrated over z, averaged inside bin)" + DO i_E = 1, n_E energy = E_min + i_E*bs_env%energy_step_DOS - WRITE (iunit, "(2F17.3)") energy*evolt, LDOS_2d_bins(i_x, i_y, i_E)*evolt + WRITE (iunit, "(2F17.3)") energy*evolt, & + LDOS_2d_bins(i_x, i_y, i_E)* & + bs_env%unit_ldos_int_z_inv_Ang2_eV END DO CALL close_file(iunit) @@ -1848,9 +1872,10 @@ CONTAINS !> \param PDOS ... !> \param band_edges ... !> \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) + 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 @@ -1863,13 +1888,15 @@ CONTAINS REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: PDOS TYPE(band_edges_type) :: band_edges REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eigenval_spinor + TYPE(cp_cfm_type) :: cfm_spinor_wf_ikp CHARACTER(LEN=*), PARAMETER :: routineN = 'SOC' INTEGER :: handle, homo_spinor, n_ao, n_E, nkind REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eigenval_spinor_no_SOC REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: proj_mo_on_kind_spinor - TYPE(cp_cfm_type) :: cfm_ks_ikp_spinor, cfm_mos_ikp_spinor, & + TYPE(cp_cfm_type) :: cfm_eigenvec_ikp_spinor, & + cfm_ks_ikp_spinor, cfm_mos_ikp_spinor, & cfm_SOC_ikp_spinor, cfm_work_ikp_spinor CALL timeset(routineN, handle) @@ -1883,6 +1910,7 @@ CONTAINS 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) ALLOCATE (eigenval_spinor_no_SOC(2*n_ao)) ALLOCATE (proj_mo_on_kind_spinor(2*n_ao, nkind)) @@ -1930,7 +1958,7 @@ CONTAINS CALL cfm_add_on_diag(cfm_ks_ikp_spinor, eigenval_spinor_no_SOC) ! 5. diagonalize h^G0W0+SOC_nn',σσ'(k_i) to get eigenvalues - CALL cp_cfm_heevd(cfm_ks_ikp_spinor, cfm_work_ikp_spinor, eigenval_spinor) + CALL cp_cfm_heevd(cfm_ks_ikp_spinor, cfm_eigenvec_ikp_spinor, eigenval_spinor) ! 6. DOS from spinors, no PDOS CALL add_to_DOS_PDOS(DOS, PDOS, eigenval_spinor, & @@ -1942,10 +1970,16 @@ CONTAINS band_edges%DBG = MIN(band_edges%DBG, eigenval_spinor(homo_spinor + 1) & - eigenval_spinor(homo_spinor)) + ! 8. spinor wavefunctions: + CALL parallel_gemm('N', 'N', 2*n_ao, 2*n_ao, 2*n_ao, z_one, & + cfm_mos_ikp_spinor, cfm_eigenvec_ikp_spinor, & + z_zero, cfm_spinor_wf_ikp) + CALL cp_cfm_release(cfm_ks_ikp_spinor) CALL cp_cfm_release(cfm_SOC_ikp_spinor) - CALL cp_cfm_release(cfm_mos_ikp_spinor) CALL cp_cfm_release(cfm_work_ikp_spinor) + CALL cp_cfm_release(cfm_eigenvec_ikp_spinor) + CALL cp_cfm_release(cfm_mos_ikp_spinor) CALL timestop(handle) @@ -2039,17 +2073,16 @@ CONTAINS CHARACTER(LEN=*), PARAMETER :: routineN = 'add_to_LDOS_2d' - INTEGER :: handle, i_E, i_x, i_x_end, i_x_start, i_y, i_y_end, i_y_start, i_z, i_z_end, & - i_z_start, j_col, j_mo, n_E, n_mo, ncol_local, nimages + INTEGER :: handle, i_E, i_x_end, i_x_start, i_y_end, i_y_start, i_z, i_z_end, i_z_start, & + j_col, j_mo, n_E, n_mo, n_z, ncol_local, nimages, z_end_global, z_start_global INTEGER, DIMENSION(:), POINTER :: col_indices LOGICAL :: is_any_weight_non_zero, my_do_spinor - REAL(KIND=dp) :: arg, E_max, E_min, E_total_window, & - energy, energy_step_ldos, & - energy_window_ldos, weight, wkp + REAL(KIND=dp) :: broadening, E_max, E_min, & + E_total_window, energy, energy_step, & + energy_window, spin_degeneracy, weight TYPE(cp_cfm_type) :: cfm_weighted_dm_ikp, cfm_work - TYPE(cp_fm_type) :: fm_non_spinor, fm_non_spinor_2, & - fm_weighted_dm_ikp_re - TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: weighted_dm_ikp_re + TYPE(cp_fm_type) :: fm_non_spinor, fm_weighted_dm_MIC + TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: weighted_dm_MIC TYPE(dft_control_type), POINTER :: dft_control TYPE(pw_env_type), POINTER :: pw_env TYPE(pw_pool_type), POINTER :: auxbas_pw_pool @@ -2067,14 +2100,15 @@ CONTAINS nimages = dft_control%nimages dft_control%nimages = 1 - energy_window_ldos = bs_env%energy_window_DOS - energy_step_ldos = bs_env%energy_step_DOS + energy_window = bs_env%energy_window_DOS + energy_step = bs_env%energy_step_DOS + broadening = bs_env%broadening_DOS - E_min = band_edges%VBM - 0.5_dp*energy_window_ldos - E_max = band_edges%CBM + 0.5_dp*energy_window_ldos + E_min = band_edges%VBM - 0.5_dp*energy_window + E_max = band_edges%CBM + 0.5_dp*energy_window E_total_window = E_max - E_min - n_E = INT(E_total_window/energy_step_ldos) + n_E = INT(E_total_window/energy_step) CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool) @@ -2088,6 +2122,18 @@ CONTAINS i_z_start = LBOUND(LDOS_3d%cr3d, 3) i_z_end = UBOUND(LDOS_3d%cr3d, 3) + z_start_global = i_z_start + z_end_global = i_z_end + + CALL bs_env%para_env%min(z_start_global) + CALL bs_env%para_env%max(z_end_global) + n_z = z_end_global - z_start_global + 1 + + IF (ANY(ABS(bs_env%hmat(1:2, 3)) > 1.0E-6_dp) .OR. ANY(ABS(bs_env%hmat(3, 1:2)) > 1.0E-6_dp)) & + CPABORT("Please choose a cell that has 90° angles to the z-direction.") + ! for integration, we need the dz and the conversion from H -> eV and a_Bohr -> Å + bs_env%unit_ldos_int_z_inv_Ang2_eV = bs_env%hmat(3, 3)/REAL(n_z, KIND=dp)/evolt/angstrom**2 + IF (ikp == 1) THEN ALLOCATE (LDOS_2d(i_x_start:i_x_end, i_y_start:i_y_end, n_E)) LDOS_2d(:, :, :) = 0.0_dp @@ -2095,10 +2141,9 @@ CONTAINS CALL cp_cfm_create(cfm_work, cfm_mos_ikp%matrix_struct) CALL cp_cfm_create(cfm_weighted_dm_ikp, cfm_mos_ikp%matrix_struct) - CALL cp_fm_create(fm_weighted_dm_ikp_re, cfm_mos_ikp%matrix_struct) + CALL cp_fm_create(fm_weighted_dm_MIC, cfm_mos_ikp%matrix_struct) IF (my_do_spinor) THEN CALL cp_fm_create(fm_non_spinor, cfm_non_spinor%matrix_struct) - CALL cp_fm_create(fm_non_spinor_2, cfm_non_spinor%matrix_struct) END IF CALL cp_cfm_get_info(matrix=cfm_mos_ikp, & @@ -2106,28 +2151,15 @@ CONTAINS ncol_local=ncol_local, & col_indices=col_indices) - NULLIFY (weighted_dm_ikp_re) - CALL dbcsr_allocate_matrix_set(weighted_dm_ikp_re, 1) - ALLOCATE (weighted_dm_ikp_re(1)%matrix) - CALL dbcsr_create(weighted_dm_ikp_re(1)%matrix, template=bs_env%mat_ao_ao%matrix, & + NULLIFY (weighted_dm_MIC) + CALL dbcsr_allocate_matrix_set(weighted_dm_MIC, 1) + ALLOCATE (weighted_dm_MIC(1)%matrix) + CALL dbcsr_create(weighted_dm_MIC(1)%matrix, template=bs_env%mat_ao_ao%matrix, & matrix_type=dbcsr_type_symmetric) - wkp = 0 - DO i_x = -bs_env%periodic(1), bs_env%periodic(1) - DO i_y = -bs_env%periodic(2), bs_env%periodic(2) - DO i_z = -bs_env%periodic(3), bs_env%periodic(3) - arg = i_x*bs_env%kpoints_DOS%xkp(1, ikp) + & - i_y*bs_env%kpoints_DOS%xkp(2, ikp) + & - i_z*bs_env%kpoints_DOS%xkp(3, ikp) - ! only cosine weight necessary; sine weight will cancel because of R <-> -R - wkp = wkp + bs_env%kpoints_DOS%wkp(ikp)*COS(twopi*arg) - END DO - END DO - END DO - DO i_E = 1, n_E - energy = E_min + i_E*energy_step_ldos + energy = E_min + i_E*energy_step is_any_weight_non_zero = .FALSE. @@ -2135,7 +2167,13 @@ CONTAINS j_mo = col_indices(j_col) - weight = wkp*Gaussian(energy - eigenval(j_mo), energy_step_ldos) + IF (my_do_spinor) THEN + spin_degeneracy = 1.0_dp + ELSE + spin_degeneracy = bs_env%spin_degeneracy + END IF + + weight = Gaussian(energy - eigenval(j_mo), broadening)*spin_degeneracy cfm_work%local_data(:, j_col) = cfm_mos_ikp%local_data(:, j_col)*weight @@ -2154,22 +2192,33 @@ CONTAINS cfm_mos_ikp, cfm_work, z_zero, cfm_weighted_dm_ikp) IF (my_do_spinor) THEN + + ! contribution from up,up to fm_non_spinor CALL get_cfm_submat(cfm_non_spinor, cfm_weighted_dm_ikp, 1, 1) - CALL cp_cfm_to_fm(cfm_non_spinor, fm_non_spinor) + CALL cp_fm_set_all(fm_non_spinor, 0.0_dp) + CALL MIC_contribution_from_ikp(bs_env, qs_env, fm_non_spinor, & + cfm_non_spinor, ikp, bs_env%kpoints_DOS, & + "ORB", bs_env%kpoints_DOS%wkp(ikp)) + + ! add contribution from down,down to fm_non_spinor CALL get_cfm_submat(cfm_non_spinor, cfm_weighted_dm_ikp, n_mo/2, n_mo/2) - CALL cp_cfm_to_fm(cfm_non_spinor, fm_non_spinor_2) - CALL cp_fm_scale_and_add(1.0_dp, fm_non_spinor, 1.0_dp, fm_non_spinor_2) - CALL copy_fm_to_dbcsr(fm_non_spinor, weighted_dm_ikp_re(1)%matrix, & + CALL MIC_contribution_from_ikp(bs_env, qs_env, fm_non_spinor, & + cfm_non_spinor, ikp, bs_env%kpoints_DOS, & + "ORB", bs_env%kpoints_DOS%wkp(ikp)) + CALL copy_fm_to_dbcsr(fm_non_spinor, weighted_dm_MIC(1)%matrix, & keep_sparsity=.FALSE.) ELSE - CALL cp_cfm_to_fm(cfm_weighted_dm_ikp, fm_weighted_dm_ikp_re) - CALL copy_fm_to_dbcsr(fm_weighted_dm_ikp_re, weighted_dm_ikp_re(1)%matrix, & + CALL cp_fm_set_all(fm_weighted_dm_MIC, 0.0_dp) + CALL MIC_contribution_from_ikp(bs_env, qs_env, fm_weighted_dm_MIC, & + cfm_weighted_dm_ikp, ikp, bs_env%kpoints_DOS, & + "ORB", bs_env%kpoints_DOS%wkp(ikp)) + CALL copy_fm_to_dbcsr(fm_weighted_dm_MIC, weighted_dm_MIC(1)%matrix, & keep_sparsity=.FALSE.) END IF LDOS_3d%cr3d(:, :, :) = 0.0_dp - CALL calculate_rho_elec(matrix_p_kp=weighted_dm_ikp_re, & + CALL calculate_rho_elec(matrix_p_kp=weighted_dm_MIC, & rho=LDOS_3d, & rho_gspace=rho_g, & ks_env=ks_env) @@ -2191,13 +2240,12 @@ CONTAINS CALL cp_cfm_release(cfm_work) CALL cp_cfm_release(cfm_weighted_dm_ikp) - CALL cp_fm_release(fm_weighted_dm_ikp_re) + CALL cp_fm_release(fm_weighted_dm_MIC) - CALL dbcsr_deallocate_matrix_set(weighted_dm_ikp_re) + CALL dbcsr_deallocate_matrix_set(weighted_dm_MIC) IF (my_do_spinor) THEN CALL cp_fm_release(fm_non_spinor) - CALL cp_fm_release(fm_non_spinor_2) END IF CALL timestop(handle) @@ -2473,8 +2521,8 @@ CONTAINS REAL(KIND=dp), INTENT(IN) :: energy, broadening REAL(KIND=dp) :: Gaussian - IF (energy < 5*broadening) THEN - Gaussian = 1/broadening/SQRT(twopi)*EXP(-energy**2/broadening**2) + IF (ABS(energy) < 5*broadening) THEN + Gaussian = 1.0_dp/broadening/SQRT(twopi)*EXP(-0.5_dp*energy**2/broadening**2) ELSE Gaussian = 0.0_dp END IF @@ -2772,4 +2820,122 @@ CONTAINS END SUBROUTINE cfm_ikp_from_fm_Gamma +! ************************************************************************************************** +!> \brief ... +!> \param bs_env ... +!> \param qs_env ... +!> \param fm_W_MIC_freq_j ... +!> \param cfm_W_ikp_freq_j ... +!> \param ikp ... +!> \param kpoints ... +!> \param basis_type ... +!> \param wkp_ext ... +! ************************************************************************************************** + SUBROUTINE MIC_contribution_from_ikp(bs_env, qs_env, fm_W_MIC_freq_j, & + cfm_W_ikp_freq_j, ikp, kpoints, basis_type, wkp_ext) + TYPE(post_scf_bandstructure_type), POINTER :: bs_env + TYPE(qs_environment_type), POINTER :: qs_env + TYPE(cp_fm_type) :: fm_W_MIC_freq_j + TYPE(cp_cfm_type) :: cfm_W_ikp_freq_j + INTEGER :: ikp + TYPE(kpoint_type), POINTER :: kpoints + CHARACTER(LEN=*) :: basis_type + REAL(KIND=dp), OPTIONAL :: wkp_ext + + CHARACTER(LEN=*), PARAMETER :: routineN = 'MIC_contribution_from_ikp' + + INTEGER :: handle, i_bf, iatom, iatom_old, irow, & + j_bf, jatom, jatom_old, jcol, n_bf, & + ncol_local, nrow_local, num_cells + INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_from_bf_index + INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices + INTEGER, DIMENSION(:, :), POINTER :: index_to_cell + REAL(KIND=dp) :: contribution, weight_im, weight_re, & + wkp_of_ikp + REAL(KIND=dp), DIMENSION(3, 3) :: hmat + REAL(KIND=dp), DIMENSION(:), POINTER :: wkp + REAL(KIND=dp), DIMENSION(:, :), POINTER :: xkp + TYPE(cell_type), POINTER :: cell + TYPE(particle_type), DIMENSION(:), POINTER :: particle_set + + CALL timeset(routineN, handle) + + ! get number of basis functions (bf) for different basis sets + IF (basis_type == "ORB") THEN + n_bf = qs_env%bs_env%n_ao + ELSE IF (basis_type == "RI_AUX") THEN + n_bf = qs_env%bs_env%n_RI + ELSE + CPABORT("Only ORB and RI_AUX basis implemented.") + END IF + + ALLOCATE (atom_from_bf_index(n_bf)) + CALL get_atom_index_from_basis_function_index(qs_env, atom_from_bf_index, n_bf, basis_type) + + NULLIFY (cell, particle_set) + CALL get_qs_env(qs_env, cell=cell, particle_set=particle_set) + CALL get_cell(cell=cell, h=hmat) + + CALL cp_cfm_get_info(matrix=cfm_W_ikp_freq_j, & + nrow_local=nrow_local, & + ncol_local=ncol_local, & + row_indices=row_indices, & + col_indices=col_indices) + + CALL get_kpoint_info(kpoints, xkp=xkp, wkp=wkp) + index_to_cell => kpoints%index_to_cell + num_cells = SIZE(index_to_cell, 2) + + iatom_old = 0 + jatom_old = 0 + + DO irow = 1, nrow_local + DO jcol = 1, ncol_local + + i_bf = row_indices(irow) + j_bf = col_indices(jcol) + + iatom = atom_from_bf_index(i_bf) + jatom = atom_from_bf_index(j_bf) + + IF (PRESENT(wkp_ext)) THEN + wkp_of_ikp = wkp_ext + ELSE + SELECT CASE (bs_env%l_RI(i_bf) + bs_env%l_RI(j_bf)) + CASE (0) + ! both RI functions are s-functions, k-extrapolation for 2D and 3D + wkp_of_ikp = wkp(ikp) + CASE (1) + ! one function is an s-function, the other a p-function, k-extrapolation for 3D + wkp_of_ikp = bs_env%wkp_s_p(ikp) + CASE DEFAULT + ! for any other matrix element of W, there is no need for extrapolation + wkp_of_ikp = bs_env%wkp_no_extra(ikp) + END SELECT + END IF + + IF (iatom .NE. iatom_old .OR. jatom .NE. jatom_old) THEN + + CALL compute_weight_re_im(weight_re, weight_im, & + num_cells, iatom, jatom, xkp(1:3, ikp), wkp_of_ikp, & + cell, index_to_cell, hmat, particle_set) + + iatom_old = iatom + jatom_old = jatom + + END IF + + contribution = weight_re*REAL(cfm_W_ikp_freq_j%local_data(irow, jcol)) + & + weight_im*AIMAG(cfm_W_ikp_freq_j%local_data(irow, jcol)) + + fm_W_MIC_freq_j%local_data(irow, jcol) = fm_W_MIC_freq_j%local_data(irow, jcol) & + + contribution + + END DO + END DO + + CALL timestop(handle) + + END SUBROUTINE MIC_contribution_from_ikp + END MODULE post_scf_bandstructure_utils