From c1b10bd6526dca92a6deb944b6159f1d8b1af7ba Mon Sep 17 00:00:00 2001 From: Hossam Elgabarty Date: Fri, 20 Sep 2024 11:42:34 +0200 Subject: [PATCH] Print wannier states coefficients in AO basis * Print Wannier states in AO basis when input section "FORCE_EVAL/DFT/LOCALIZE/PRINT/WANNIER_STATES" is used * Modify behaviour of the above input section so that it is not ignored when input section " ... /PRINT/WANNIER_CENTERS" is missing --- src/input_cp2k_loc.F | 12 +- src/qs_scf_post_gpw.F | 4 +- src/wannier_states.F | 275 ++++++++++++++++++++++++++++++++++++++++-- 3 files changed, 279 insertions(+), 12 deletions(-) diff --git a/src/input_cp2k_loc.F b/src/input_cp2k_loc.F index 92492352fd..6c74dc498f 100644 --- a/src/input_cp2k_loc.F +++ b/src/input_cp2k_loc.F @@ -454,7 +454,7 @@ SUBROUTINE create_molecular_states_section(print_key) CALL keyword_create( & keyword, __LOCATION__, name="CUBE_EVAL_RANGE", & - description="only write cubes if the eigenvalues of the corresponding molecular states lie in the given interval. "// & + description="only write cubes if the energies of the corresponding molecular states lie in the given interval. "// & "Default is all states.", & usage="CUBE_EVAL_RANGE -1.0 1.0", unit_str="hartree", n_var=2, type_of_var=real_t) CALL section_add_keyword(print_key, keyword) @@ -499,12 +499,12 @@ SUBROUTINE create_wannier_states_section(print_key) CPASSERT(.NOT. ASSOCIATED(print_key)) NULLIFY (print_key2, keyword) CALL cp_print_key_section_create(print_key, __LOCATION__, "WANNIER_STATES", & - description="Controls printing of molecular states ", & + description="Controls printing of Wannier states ", & print_level=high_print_level, filename=" ") CALL keyword_create( & keyword, __LOCATION__, name="CUBE_EVAL_RANGE", & - description="only write cubes if the eigenvalues of the corresponding molecular states lie in the given interval. "// & + description="only write cubes if the energies of the corresponding molecular states lie in the given interval. "// & "Default is all states.", & usage="CUBE_EVAL_RANGE -1.0 1.0", unit_str="hartree", n_var=2, type_of_var=real_t) CALL section_add_keyword(print_key, keyword) @@ -522,6 +522,12 @@ SUBROUTINE create_wannier_states_section(print_key) CALL section_add_keyword(print_key, keyword) CALL keyword_release(keyword) + CALL keyword_create(keyword, __LOCATION__, name="CARTESIAN", & + description="Print the Wannier states in the Cartesian basis instead of the default spherical basis.", & + default_l_val=.FALSE., lone_keyword_l_val=.TRUE.) + CALL section_add_keyword(print_key, keyword) + CALL keyword_release(keyword) + CALL cp_print_key_section_create(print_key2, __LOCATION__, "cubes", & description="Controls the printing of cube files", & print_level=high_print_level, filename="") diff --git a/src/qs_scf_post_gpw.F b/src/qs_scf_post_gpw.F index 9860a16fc4..ea80310bd6 100644 --- a/src/qs_scf_post_gpw.F +++ b/src/qs_scf_post_gpw.F @@ -400,6 +400,8 @@ SUBROUTINE scf_post_calculation_gpw(qs_env, wf_type, do_mp2) p_loc = p_loc .OR. BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file) print_key => section_vals_get_subs_vals(loc_print_section, "MOLECULAR_MOMENTS") p_loc = p_loc .OR. BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file) + print_key => section_vals_get_subs_vals(loc_print_section, "WANNIER_STATES") + p_loc = p_loc .OR. BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file) ELSE p_loc = .FALSE. END IF @@ -422,8 +424,6 @@ SUBROUTINE scf_post_calculation_gpw(qs_env, wf_type, do_mp2) "therefore localization of unoccupied states will be skipped!") p_loc_lumo = .FALSE. END IF - print_key => section_vals_get_subs_vals(loc_print_section, "WANNIER_STATES") - p_loc = p_loc .OR. BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file) ! Control for STM stm_section => section_vals_get_subs_vals(input, "DFT%PRINT%STM") diff --git a/src/wannier_states.F b/src/wannier_states.F index 497f4e6d9e..629b1fccd6 100644 --- a/src/wannier_states.F +++ b/src/wannier_states.F @@ -10,6 +10,11 @@ !> \author Alin M Elena ! ************************************************************************************************** MODULE wannier_states + USE atomic_kind_types, ONLY: atomic_kind_type,& + get_atomic_kind,& + get_atomic_kind_set + USE basis_set_types, ONLY: get_gto_basis_set,& + gto_basis_set_type USE cp_dbcsr_api, ONLY: dbcsr_type USE cp_dbcsr_operations, ONLY: cp_dbcsr_sm_fm_multiply USE cp_fm_struct, ONLY: cp_fm_struct_create,& @@ -18,6 +23,7 @@ MODULE wannier_states USE cp_fm_types, ONLY: cp_fm_create,& cp_fm_get_element,& cp_fm_get_info,& + cp_fm_get_submatrix,& cp_fm_release,& cp_fm_to_fm,& cp_fm_type @@ -33,9 +39,21 @@ MODULE wannier_states USE kinds, ONLY: default_string_length,& dp USE message_passing, ONLY: mp_para_env_type + USE orbital_pointers, ONLY: indco,& + nco,& + nso + USE orbital_symbols, ONLY: cgf_symbol,& + sgf_symbol + USE orbital_transformation_matrices, ONLY: orbtramat USE parallel_gemm_api, ONLY: parallel_gemm + USE particle_types, ONLY: particle_type + USE qs_dftb_types, ONLY: qs_dftb_atom_type + USE qs_dftb_utils, ONLY: get_dftb_atom_param USE qs_environment_types, ONLY: get_qs_env,& qs_environment_type + USE qs_kind_types, ONLY: get_qs_kind,& + get_qs_kind_set,& + qs_kind_type USE wannier_states_types, ONLY: wannier_centres_type !!!! this ones are needed to mapping #include "./base/base_uses.f90" @@ -65,6 +83,8 @@ MODULE wannier_states !> \param WannierCentres ... !> \param ns ... !> \param states ... +!> \par History +!> - Printout Wannier states in AO basis (11.09.2025, H. Elgabarty) ! ************************************************************************************************** SUBROUTINE construct_wannier_states(mo_localized, & Hks, qs_env, loc_print_section, WannierCentres, ns, states) @@ -80,14 +100,28 @@ SUBROUTINE construct_wannier_states(mo_localized, & CHARACTER(len=*), PARAMETER :: routineN = 'construct_wannier_states' CHARACTER(default_string_length) :: unit_str - INTEGER :: handle, i, iproc, ncol_global, nproc, & - nrow_global, nstates(2), output_unit, & - unit_mat + CHARACTER(LEN=12) :: symbol + CHARACTER(LEN=12), DIMENSION(:), POINTER :: bcgf_symbol + CHARACTER(LEN=2) :: element_symbol + CHARACTER(LEN=40) :: fmtstr1, fmtstr2, fmtstr3 + CHARACTER(LEN=6), DIMENSION(:), POINTER :: bsgf_symbol + INTEGER :: after = 6, before = 4, from, handle, i, iatom, icgf, ico, icol, ikind, iproc, & + irow, iset, isgf, ishell, iso, jcol, left, lmax, lshell, natom, ncgf, ncol, ncol_global, & + nproc, nrow_global, nset, nsgf, nstates(2), output_unit, right, to, unit_mat + INTEGER, DIMENSION(:), POINTER :: nshell + INTEGER, DIMENSION(:, :), POINTER :: l + LOGICAL :: ionode, print_cartesian REAL(KIND=dp) :: unit_conv + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: cmatrix, smatrix + TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp TYPE(cp_fm_type) :: b, c, d TYPE(cp_logger_type), POINTER :: logger + TYPE(gto_basis_set_type), POINTER :: orb_basis_set TYPE(mp_para_env_type), POINTER :: para_env + TYPE(particle_type), DIMENSION(:), POINTER :: particle_set + TYPE(qs_dftb_atom_type), POINTER :: dftb_parameter + TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set TYPE(section_vals_type), POINTER :: print_key !----------------------------------------------------------------------- @@ -96,10 +130,16 @@ SUBROUTINE construct_wannier_states(mo_localized, & CALL timeset(routineN, handle) NULLIFY (logger, para_env) - CALL get_qs_env(qs_env, para_env=para_env) + CALL get_qs_env(qs_env, para_env=para_env, & + atomic_kind_set=atomic_kind_set, & + qs_kind_set=qs_kind_set, & + particle_set=particle_set) + nproc = para_env%num_pe logger => cp_get_default_logger() + ionode = logger%para_env%is_source() + output_unit = cp_logger_get_default_io_unit(logger) CALL cp_fm_get_info(mo_localized, & ncol_global=ncol_global, & @@ -114,6 +154,8 @@ SUBROUTINE construct_wannier_states(mo_localized, & CALL section_vals_val_get(print_key, "UNIT", c_val=unit_str) unit_conv = cp_unit_from_cp2k(1.0_dp, TRIM(unit_str)) + print_key => section_vals_get_subs_vals(loc_print_section, "WANNIER_STATES") + CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=nrow_global, & ncol_global=1, & para_env=mo_localized%matrix_struct%para_env, & @@ -132,7 +174,6 @@ SUBROUTINE construct_wannier_states(mo_localized, & CALL cp_fm_struct_release(fm_struct_tmp) WannierCentres%WannierHamDiag = 0.0_dp - ! try to print the matrix unit_mat = cp_print_key_unit_nr(logger, loc_print_section, & "WANNIER_STATES", extension=".whks", & @@ -153,8 +194,7 @@ SUBROUTINE construct_wannier_states(mo_localized, & END DO IF (unit_mat > 0) WRITE (unit_mat, *) - CALL cp_print_key_finished_output(unit_mat, logger, loc_print_section, & - "WANNIER_STATES") + IF (output_unit > 0) THEN WRITE (output_unit, *) "" WRITE (output_unit, *) "NUMBER OF Wannier STATES ", ns @@ -163,9 +203,230 @@ SUBROUTINE construct_wannier_states(mo_localized, & WRITE (output_unit, '(f16.8,2x,i0)') WannierCentres%WannierHamDiag(i), states(i) END DO END IF + CALL cp_fm_release(b) CALL cp_fm_release(c) CALL cp_fm_release(d) + + ! Print the states in AO basis + CALL section_vals_val_get(print_key, "CARTESIAN", l_val=print_cartesian) + + IF (unit_mat > 0) THEN + + NULLIFY (nshell) + NULLIFY (bsgf_symbol) + NULLIFY (l) + NULLIFY (bsgf_symbol) + + ALLOCATE (smatrix(nrow_global, ncol_global)) + CALL cp_fm_get_submatrix(mo_localized, smatrix(1:nrow_global, 1:ncol_global)) + + IF (.NOT. ionode) THEN + DEALLOCATE (smatrix) + END IF + + CALL get_atomic_kind_set(atomic_kind_set, natom=natom) + CALL get_qs_kind_set(qs_kind_set, ncgf=ncgf, nsgf=nsgf) + + ncol = INT(56/(before + after + 3)) + + fmtstr1 = "(T2,A,21X, ( X,I5, X))" + fmtstr2 = "(T2,A,9X, (1X,F . ))" + fmtstr3 = "(T2,A,I5,1X,I5,1X,A,1X,A6, (1X,F . ))" + + right = MAX((after - 2), 1) + left = (before + after + 3) - right - 5 + + IF (print_cartesian) THEN + WRITE (UNIT=unit_mat, FMT="(T2,A,16X,A)") "WS|", "Wannier states in the cartesian AO basis" + ELSE + WRITE (UNIT=unit_mat, FMT="(T2,A,16X,A)") "WS|", "Wannier states in the spherical AO basis" + END IF + WRITE (UNIT=fmtstr1(11:12), FMT="(I2)") ncol + WRITE (UNIT=fmtstr1(14:15), FMT="(I2)") left + WRITE (UNIT=fmtstr1(21:22), FMT="(I2)") right + + WRITE (UNIT=fmtstr2(10:11), FMT="(I2)") ncol + WRITE (UNIT=fmtstr2(17:18), FMT="(I2)") before + after + 2 + WRITE (UNIT=fmtstr2(20:21), FMT="(I2)") after + + WRITE (UNIT=fmtstr3(27:28), FMT="(I2)") ncol + WRITE (UNIT=fmtstr3(34:35), FMT="(I2)") before + after + 2 + WRITE (UNIT=fmtstr3(37:38), FMT="(I2)") after + + ! get MO coefficients in terms of contracted cartesian functions + IF (print_cartesian) THEN + + ALLOCATE (cmatrix(ncgf, ncgf)) + cmatrix = 0.0_dp + + ! Transform spherical MOs to Cartesian MOs + icgf = 1 + isgf = 1 + DO iatom = 1, natom + NULLIFY (orb_basis_set, dftb_parameter) + CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind) + CALL get_qs_kind(qs_kind_set(ikind), & + basis_set=orb_basis_set, & + dftb_parameter=dftb_parameter) + IF (ASSOCIATED(orb_basis_set)) THEN + CALL get_gto_basis_set(gto_basis_set=orb_basis_set, & + nset=nset, & + nshell=nshell, & + l=l) + DO iset = 1, nset + DO ishell = 1, nshell(iset) + lshell = l(ishell, iset) + CALL dgemm("T", "N", nco(lshell), ncol_global, nso(lshell), 1.0_dp, & + orbtramat(lshell)%c2s, nso(lshell), & + smatrix(isgf, 1), nsgf, 0.0_dp, & + cmatrix(icgf, 1), ncgf) + icgf = icgf + nco(lshell) + isgf = isgf + nso(lshell) + END DO + END DO + ELSE IF (ASSOCIATED(dftb_parameter)) THEN + CALL get_dftb_atom_param(dftb_parameter, lmax=lmax) + DO ishell = 1, lmax + 1 + lshell = ishell - 1 + CALL dgemm("T", "N", nco(lshell), nsgf, nso(lshell), 1.0_dp, & + orbtramat(lshell)%c2s, nso(lshell), & + smatrix(isgf, 1), nsgf, 0.0_dp, & + cmatrix(icgf, 1), ncgf) + icgf = icgf + nco(lshell) + isgf = isgf + nso(lshell) + END DO + ELSE + ! assume atom without basis set + END IF + END DO ! iatom + + END IF + + ! Print to file + DO icol = 1, ncol_global, ncol + + from = icol + to = MIN((from + ncol - 1), ncol_global) + + WRITE (UNIT=unit_mat, FMT="(T2,A)") "WS|" + WRITE (UNIT=unit_mat, FMT=fmtstr1) "WS|", (jcol, jcol=from, to) + WRITE (UNIT=unit_mat, FMT=fmtstr2) "WS| Energies", & + (WannierCentres%WannierHamDiag(states(jcol)), jcol=from, to) + WRITE (UNIT=unit_mat, FMT="(T2,A)") "WS|" + + irow = 1 + + DO iatom = 1, natom + + IF (iatom /= 1) WRITE (UNIT=unit_mat, FMT="(T2,A)") "WS|" + + NULLIFY (orb_basis_set, dftb_parameter) + CALL get_atomic_kind(particle_set(iatom)%atomic_kind, & + element_symbol=element_symbol, kind_number=ikind) + CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set, & + dftb_parameter=dftb_parameter) + + IF (print_cartesian) THEN + + NULLIFY (bcgf_symbol) + IF (ASSOCIATED(orb_basis_set)) THEN + CALL get_gto_basis_set(gto_basis_set=orb_basis_set, & + nset=nset, & + nshell=nshell, & + l=l, & + cgf_symbol=bcgf_symbol) + + icgf = 1 + DO iset = 1, nset + DO ishell = 1, nshell(iset) + lshell = l(ishell, iset) + DO ico = 1, nco(lshell) + WRITE (UNIT=unit_mat, FMT=fmtstr3) & + "WS|", irow, iatom, ADJUSTR(element_symbol), bcgf_symbol(icgf), & + (cmatrix(irow, jcol), jcol=from, to) + icgf = icgf + 1 + irow = irow + 1 + END DO + END DO + END DO + ELSE IF (ASSOCIATED(dftb_parameter)) THEN + CALL get_dftb_atom_param(dftb_parameter, lmax=lmax) + icgf = 1 + DO ishell = 1, lmax + 1 + lshell = ishell - 1 + DO ico = 1, nco(lshell) + symbol = cgf_symbol(1, indco(1:3, icgf)) + symbol(1:2) = " " + WRITE (UNIT=unit_mat, FMT=fmtstr3) & + "WS|", irow, iatom, ADJUSTR(element_symbol), symbol, & + (cmatrix(irow, jcol), jcol=from, to) + icgf = icgf + 1 + irow = irow + 1 + END DO + END DO + ELSE + ! assume atom without basis set + END IF + + ELSE !print in spherical AO basis + + IF (ASSOCIATED(orb_basis_set)) THEN + CALL get_gto_basis_set(gto_basis_set=orb_basis_set, & + nset=nset, & + nshell=nshell, & + l=l, & + sgf_symbol=bsgf_symbol) + isgf = 1 + DO iset = 1, nset + DO ishell = 1, nshell(iset) + lshell = l(ishell, iset) + DO iso = 1, nso(lshell) + WRITE (UNIT=unit_mat, FMT=fmtstr3) & + "WS|", irow, iatom, ADJUSTR(element_symbol), bsgf_symbol(isgf), & + (smatrix(irow, jcol), jcol=from, to) + isgf = isgf + 1 + irow = irow + 1 + END DO + END DO + END DO + ELSE IF (ASSOCIATED(dftb_parameter)) THEN + CALL get_dftb_atom_param(dftb_parameter, lmax=lmax) + isgf = 1 + DO ishell = 1, lmax + 1 + lshell = ishell - 1 + DO iso = 1, nso(lshell) + symbol = sgf_symbol(1, lshell, -lshell + iso - 1) + symbol(1:2) = " " + WRITE (UNIT=unit_mat, FMT=fmtstr3) & + "WS|", irow, iatom, ADJUSTR(element_symbol), symbol, & + (smatrix(irow, jcol), jcol=from, to) + isgf = isgf + 1 + irow = irow + 1 + END DO + END DO + ELSE + ! assume atom without basis set + END IF + + END IF ! print cartesian + + END DO ! iatom + + END DO ! icol + + WRITE (UNIT=unit_mat, FMT="(T2,A)") "MO|" + + DEALLOCATE (smatrix) + IF (print_cartesian) THEN + DEALLOCATE (cmatrix) + END IF + + END IF ! output Wannier states in AO + + CALL cp_print_key_finished_output(unit_mat, logger, loc_print_section, & + "WANNIER_STATES") + CALL timestop(handle) END SUBROUTINE construct_wannier_states