Implement k-point PDOS with Lowdin projections (#5299)

This commit is contained in:
SY Wang 2026-05-28 05:09:27 +08:00 committed by GitHub
parent cdd8209268
commit ddfed05c5a
No known key found for this signature in database
GPG key ID: B5690EEEBB952194
6 changed files with 551 additions and 57 deletions

View file

@ -19,6 +19,7 @@ MODULE kpoint_methods
USE cp_blacs_env, ONLY: cp_blacs_env_create,&
cp_blacs_env_type
USE cp_cfm_types, ONLY: cp_cfm_create,&
cp_cfm_get_info,&
cp_cfm_release,&
cp_cfm_to_fm,&
cp_cfm_type,&
@ -104,7 +105,7 @@ MODULE kpoint_methods
PUBLIC :: kpoint_initialize, kpoint_env_initialize, kpoint_initialize_mos, kpoint_initialize_mo_set
PUBLIC :: kpoint_init_cell_index, kpoint_set_mo_occupation
PUBLIC :: kpoint_density_matrices, kpoint_density_transform
PUBLIC :: rskp_transform, lowdin_kp_trans
PUBLIC :: rskp_transform, lowdin_kp_trans, lowdin_kp_mo_coeff
! **************************************************************************************************
@ -1647,6 +1648,71 @@ CONTAINS
END SUBROUTINE lowdin_kp_trans
! **************************************************************************************************
!> \brief Calculate S(k)^1/2 C(k) for real or complex k-point wavefunctions
!> \param kp K-point environment for one local k point
!> \param ispin Spin index
!> \param use_real_wfn Use real k-point wavefunctions
!> \param shalfc Output matrix containing S(k)^1/2 C(k) for real wavefunctions
!> \param cshalfc Output matrix containing S(k)^1/2 C(k) for complex wavefunctions
! **************************************************************************************************
SUBROUTINE lowdin_kp_mo_coeff(kp, ispin, use_real_wfn, shalfc, cshalfc)
TYPE(kpoint_env_type), POINTER :: kp
INTEGER, INTENT(IN) :: ispin
LOGICAL, INTENT(IN) :: use_real_wfn
TYPE(cp_fm_type), INTENT(INOUT), OPTIONAL :: shalfc
TYPE(cp_cfm_type), INTENT(INOUT), OPTIONAL :: cshalfc
INTEGER :: nao, nmo
TYPE(cp_fm_struct_type), POINTER :: matrix_struct_mo, matrix_struct_shalf
TYPE(cp_fm_type) :: cshalf_im, cshalf_re, shalf_im, shalf_re
TYPE(mo_set_type), POINTER :: mo_set, mo_set_im, mo_set_re
IF (use_real_wfn) THEN
CPASSERT(PRESENT(shalfc))
mo_set => kp%mos(1, ispin)
CALL get_mo_set(mo_set, nao=nao, nmo=nmo)
CALL parallel_gemm("N", "N", nao, nmo, nao, 1.0_dp, kp%shalf, &
mo_set%mo_coeff, 0.0_dp, shalfc)
ELSE
CPASSERT(PRESENT(cshalfc))
mo_set_re => kp%mos(1, ispin)
mo_set_im => kp%mos(2, ispin)
CALL get_mo_set(mo_set_re, nao=nao, nmo=nmo)
CALL cp_fm_get_info(mo_set_re%mo_coeff, matrix_struct=matrix_struct_mo)
CALL cp_cfm_get_info(kp%cshalf, matrix_struct=matrix_struct_shalf)
CALL cp_fm_create(shalf_re, matrix_struct_shalf, nrow=nao, ncol=nao)
CALL cp_fm_create(shalf_im, matrix_struct_shalf, nrow=nao, ncol=nao)
CALL cp_fm_create(cshalf_re, matrix_struct_mo, nrow=nao, ncol=nmo)
CALL cp_fm_create(cshalf_im, matrix_struct_mo, nrow=nao, ncol=nmo)
CALL cp_cfm_to_fm(kp%cshalf, mtargetr=shalf_re, mtargeti=shalf_im)
! Re[S(k)^1/2 C(k)] = Re[S(k)^1/2] C_re(k) - Im[S(k)^1/2] C_im(k)
CALL parallel_gemm("N", "N", nao, nmo, nao, 1.0_dp, shalf_re, &
mo_set_re%mo_coeff, 0.0_dp, cshalf_re)
CALL parallel_gemm("N", "N", nao, nmo, nao, -1.0_dp, shalf_im, &
mo_set_im%mo_coeff, 1.0_dp, cshalf_re)
! Im[S(k)^1/2 C(k)] = Re[S(k)^1/2] C_im(k) + Im[S(k)^1/2] C_re(k)
CALL parallel_gemm("N", "N", nao, nmo, nao, 1.0_dp, shalf_re, &
mo_set_im%mo_coeff, 0.0_dp, cshalf_im)
CALL parallel_gemm("N", "N", nao, nmo, nao, 1.0_dp, shalf_im, &
mo_set_re%mo_coeff, 1.0_dp, cshalf_im)
CALL cp_fm_to_cfm(cshalf_re, cshalf_im, cshalfc)
CALL cp_fm_release(shalf_re)
CALL cp_fm_release(shalf_im)
CALL cp_fm_release(cshalf_re)
CALL cp_fm_release(cshalf_im)
END IF
END SUBROUTINE lowdin_kp_mo_coeff
! **************************************************************************************************
!> \brief generate real space density matrices in DBCSR format
!> \param kpoint Kpoint environment

View file

@ -24,7 +24,8 @@ MODULE qs_dos_utils
broadening_lorentzian = 2, &
broadening_pseudo_voigt = 3
PUBLIC :: add_broadened_peak, broadening_cutoff, broadening_function, write_broadening_info
PUBLIC :: add_broadened_peak, add_broadened_value, broadening_cutoff, broadening_function, &
write_broadening_info
CONTAINS
@ -64,6 +65,39 @@ CONTAINS
END SUBROUTINE add_broadened_peak
! **************************************************************************************************
!> \brief Add a broadened spectral line with a scalar weight to a curve.
!> \param curve ...
!> \param emin ...
!> \param de ...
!> \param eig ...
!> \param weight ...
!> \param broaden_type ...
!> \param broaden_width ...
!> \param voigt_mixing ...
! **************************************************************************************************
SUBROUTINE add_broadened_value(curve, emin, de, eig, weight, broaden_type, &
broaden_width, voigt_mixing)
REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: curve
REAL(KIND=dp), INTENT(IN) :: emin, de, eig, weight
INTEGER, INTENT(IN) :: broaden_type
REAL(KIND=dp), INTENT(IN) :: broaden_width, voigt_mixing
INTEGER :: i, iend, istart, nhist
REAL(KIND=dp) :: eval, line_shape
nhist = SIZE(curve)
istart = MAX(1, FLOOR((eig - broadening_cutoff(broaden_type, broaden_width) - emin)/de) + 1)
iend = MIN(nhist, CEILING((eig + broadening_cutoff(broaden_type, broaden_width) - emin)/de) + 1)
DO i = istart, iend
eval = emin + (i - 1)*de
line_shape = broadening_function(eval - eig, broaden_type, broaden_width, voigt_mixing)
curve(i) = curve(i) + weight*line_shape
END DO
END SUBROUTINE add_broadened_value
! **************************************************************************************************
!> \brief Broadening cutoff used for numerical accumulation.
!> \param broaden_type ...

View file

@ -22,6 +22,10 @@ MODULE qs_pdos
pbc
USE cp_array_utils, ONLY: cp_1d_r_p_type
USE cp_blacs_env, ONLY: cp_blacs_env_type
USE cp_cfm_types, ONLY: cp_cfm_create,&
cp_cfm_get_submatrix,&
cp_cfm_release,&
cp_cfm_type
USE cp_control_types, ONLY: dft_control_type
USE cp_dbcsr_api, ONLY: dbcsr_p_type
USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm
@ -48,6 +52,9 @@ MODULE qs_pdos
section_vals_val_get
USE kinds, ONLY: default_string_length,&
dp
USE kpoint_methods, ONLY: lowdin_kp_mo_coeff
USE kpoint_types, ONLY: kpoint_env_type,&
kpoint_type
USE memory_utilities, ONLY: reallocate
USE message_passing, ONLY: mp_para_env_type
USE orbital_pointers, ONLY: nso,&
@ -63,7 +70,8 @@ MODULE qs_pdos
USE pw_types, ONLY: pw_c1d_gs_type,&
pw_r3d_rs_type
USE qs_collocate_density, ONLY: calculate_wavefunction
USE qs_dos_utils, ONLY: broadening_cutoff,&
USE qs_dos_utils, ONLY: add_broadened_value,&
broadening_cutoff,&
broadening_function,&
write_broadening_info
USE qs_environment_types, ONLY: get_qs_env,&
@ -73,6 +81,8 @@ MODULE qs_pdos
qs_kind_type
USE qs_mo_types, ONLY: get_mo_set,&
mo_set_type
USE qs_scf_diagonalization, ONLY: diag_kp_smat
USE qs_scf_types, ONLY: qs_scf_env_type
#include "./base/base_uses.f90"
IMPLICIT NONE
@ -84,7 +94,7 @@ MODULE qs_pdos
! **************************************************************************************************
! *** Public subroutines ***
PUBLIC :: calculate_projected_dos
PUBLIC :: calculate_projected_dos, calculate_projected_dos_kp
TYPE ldos_type
INTEGER :: maxl = -1, nlist = -1
@ -948,6 +958,422 @@ CONTAINS
END SUBROUTINE calculate_projected_dos
! **************************************************************************************************
!> \brief Compute and write broadened projected density of states for k-point calculations.
!> \param qs_env ...
!> \param dft_section ...
! **************************************************************************************************
SUBROUTINE calculate_projected_dos_kp(qs_env, dft_section)
TYPE(qs_environment_type), POINTER :: qs_env
TYPE(section_vals_type), POINTER :: dft_section
CHARACTER(len=*), PARAMETER :: routineN = 'calculate_projected_dos_kp'
CHARACTER(LEN=default_string_length) :: kind_name, my_act, my_mittle, my_pos, &
spin(2)
COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: zvecbuffer
INTEGER :: broaden_type, handle, icomp, ik, ikind, imo, ispin, iterstep, maxl, maxlgto, &
n_r_ldos, nao, ncomp, ndigits, nhist, nkind, nldos, nmo_kp, nspins, output_unit
INTEGER, ALLOCATABLE, DIMENSION(:) :: ao_comp, ao_kind, ao_l, kind_maxl
LOGICAL :: append, separate_components, &
should_output
REAL(KIND=dp) :: broaden_cutoff, broaden_width, de, e1, &
e2, e_fermi, emax, emin, voigt_mixing, &
wkp
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: ao_weight
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: proj_weight, vecbuffer
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: pdos_curve
REAL(KIND=dp), DIMENSION(:), POINTER :: eigenvalues
TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
TYPE(cp_cfm_type) :: cshalfc
TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp
TYPE(cp_fm_type) :: shalfc
TYPE(cp_logger_type), POINTER :: logger
TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s_kp
TYPE(dft_control_type), POINTER :: dft_control
TYPE(gto_basis_set_type), POINTER :: orb_basis_set
TYPE(kpoint_env_type), POINTER :: kp
TYPE(kpoint_type), POINTER :: kpoints
TYPE(mo_set_type), POINTER :: mo_set
TYPE(mp_para_env_type), POINTER :: para_env
TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
TYPE(qs_scf_env_type), POINTER :: scf_env
TYPE(section_vals_type), POINTER :: ldos_section
NULLIFY (logger, kpoints, dft_control, para_env, atomic_kind_set, qs_kind_set, particle_set)
NULLIFY (matrix_s_kp, scf_env)
NULLIFY (kp, mo_set, eigenvalues, fm_struct_tmp, orb_basis_set, ldos_section)
logger => cp_get_default_logger()
should_output = BTEST(cp_print_key_should_output(logger%iter_info, dft_section, &
"PRINT%PDOS"), cp_p_file)
output_unit = cp_logger_get_default_io_unit(logger)
IF ((.NOT. should_output)) RETURN
CALL timeset(routineN, handle)
iterstep = logger%iter_info%iteration(logger%iter_info%n_rlevel)
IF (output_unit > 0) WRITE (UNIT=output_unit, FMT='(/,(T3,A,T61,I10))') &
" Calculate k-point PDOS at iteration step ", iterstep
CALL get_qs_env(qs_env=qs_env, &
kpoints=kpoints, &
dft_control=dft_control, &
matrix_s_kp=matrix_s_kp, &
scf_env=scf_env, &
atomic_kind_set=atomic_kind_set, &
qs_kind_set=qs_kind_set, &
particle_set=particle_set)
para_env => kpoints%para_env_inter_kp
nspins = dft_control%nspins
nkind = SIZE(atomic_kind_set)
CALL get_qs_kind_set(qs_kind_set, maxlgto=maxlgto)
IF (.NOT. ASSOCIATED(kpoints%kp_env)) THEN
CPWARN("No local k points available for k-point PDOS")
CALL timestop(handle)
RETURN
END IF
IF (SIZE(kpoints%kp_env) == 0) THEN
CPWARN("No local k points available for k-point PDOS")
CALL timestop(handle)
RETURN
END IF
CALL section_vals_val_get(dft_section, "PRINT%PDOS%DELTA_E", r_val=de)
CALL section_vals_val_get(dft_section, "PRINT%PDOS%APPEND", l_val=append)
CALL section_vals_val_get(dft_section, "PRINT%PDOS%COMPONENTS", l_val=separate_components)
CALL section_vals_val_get(dft_section, "PRINT%PDOS%NDIGITS", i_val=ndigits)
CALL section_vals_val_get(dft_section, "PRINT%PDOS%BROADEN_TYPE", i_val=broaden_type)
CALL section_vals_val_get(dft_section, "PRINT%PDOS%BROADEN_WIDTH", r_val=broaden_width)
CALL section_vals_val_get(dft_section, "PRINT%PDOS%VOIGT_MIXING", r_val=voigt_mixing)
ndigits = MIN(MAX(ndigits, 1), 10)
de = MAX(de, 0.00001_dp)
ldos_section => section_vals_get_subs_vals(dft_section, "PRINT%PDOS%LDOS")
CALL section_vals_get(ldos_section, n_repetition=nldos)
ldos_section => section_vals_get_subs_vals(dft_section, "PRINT%PDOS%R_LDOS")
CALL section_vals_get(ldos_section, n_repetition=n_r_ldos)
IF (nldos > 0 .OR. n_r_ldos > 0) THEN
CPWARN("LDOS/R_LDOS are not implemented for k-point PDOS and will be ignored")
END IF
IF (broaden_width <= 0.0_dp) THEN
CPWARN("Raw k-point PDOS output is not implemented; use a finite BROADEN_WIDTH")
CALL timestop(handle)
RETURN
END IF
IF (separate_components) THEN
ncomp = nsoset(maxlgto)
ELSE
ncomp = maxlgto + 1
END IF
ALLOCATE (kind_maxl(nkind))
kind_maxl = -1
DO ikind = 1, nkind
NULLIFY (orb_basis_set)
CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set)
CALL get_gto_basis_set(gto_basis_set=orb_basis_set, maxl=maxl)
kind_maxl(ikind) = maxl
END DO
emin = HUGE(0.0_dp)
emax = -HUGE(0.0_dp)
e_fermi = 0.0_dp
IF (kpoints%nkp /= 0) THEN
DO ik = 1, SIZE(kpoints%kp_env)
kp => kpoints%kp_env(ik)%kpoint_env
DO ispin = 1, nspins
mo_set => kp%mos(1, ispin)
CALL get_mo_set(mo_set=mo_set, nmo=nmo_kp, mu=e_fermi)
eigenvalues => mo_set%eigenvalues
e1 = MINVAL(eigenvalues(1:nmo_kp))
e2 = MAXVAL(eigenvalues(1:nmo_kp))
emin = MIN(emin, e1)
emax = MAX(emax, e2)
END DO
END DO
END IF
CALL para_env%min(emin)
CALL para_env%max(emax)
broaden_cutoff = broadening_cutoff(broaden_type, broaden_width)
emin = emin - broaden_cutoff
emax = emax + broaden_cutoff
nhist = NINT((emax - emin)/de) + 1
ALLOCATE (pdos_curve(nhist, ncomp, nkind, nspins))
pdos_curve = 0.0_dp
! Ensure that S(k)^1/2 is available for the Lowdin projection.
! This is normally only constructed for Lowdin population/DFT+U paths.
CALL diag_kp_smat(matrix_s_kp, kpoints, scf_env%scf_work1)
! Use the first local k point to construct the AO -> kind/l/component map.
kp => kpoints%kp_env(1)%kpoint_env
mo_set => kp%mos(1, 1)
CALL get_mo_set(mo_set=mo_set, nao=nao)
CALL build_pdos_ao_map(qs_kind_set, particle_set, nao, ao_kind, ao_l, ao_comp)
ALLOCATE (ao_weight(nao), proj_weight(ncomp, nkind))
ALLOCATE (vecbuffer(1, nao), zvecbuffer(1, nao))
IF (kpoints%nkp /= 0) THEN
DO ik = 1, SIZE(kpoints%kp_env)
kp => kpoints%kp_env(ik)%kpoint_env
wkp = kp%wkp
DO ispin = 1, nspins
mo_set => kp%mos(1, ispin)
CALL get_mo_set(mo_set=mo_set, nao=nao, nmo=nmo_kp)
eigenvalues => mo_set%eigenvalues
CALL cp_fm_get_info(mo_set%mo_coeff, matrix_struct=fm_struct_tmp)
IF (kpoints%use_real_wfn) THEN
CALL cp_fm_create(shalfc, fm_struct_tmp, name="shalfc")
CALL lowdin_kp_mo_coeff(kp, ispin, kpoints%use_real_wfn, shalfc=shalfc)
ELSE
CALL cp_cfm_create(cshalfc, fm_struct_tmp, name="cshalfc")
CALL lowdin_kp_mo_coeff(kp, ispin, kpoints%use_real_wfn, cshalfc=cshalfc)
END IF
DO imo = 1, nmo_kp
IF (kpoints%use_real_wfn) THEN
CALL cp_fm_get_submatrix(shalfc, vecbuffer, 1, imo, nao, 1, transpose=.TRUE.)
ao_weight(:) = vecbuffer(1, 1:nao)**2
ELSE
CALL cp_cfm_get_submatrix(cshalfc, zvecbuffer, 1, imo, nao, 1, transpose=.TRUE.)
ao_weight(:) = REAL(CONJG(zvecbuffer(1, 1:nao))*zvecbuffer(1, 1:nao), KIND=dp)
END IF
proj_weight = 0.0_dp
CALL accumulate_pdos_weights(ao_weight, ao_kind, ao_l, ao_comp, &
separate_components, proj_weight)
DO ikind = 1, nkind
IF (kind_maxl(ikind) < 0) CYCLE
IF (separate_components) THEN
DO icomp = 1, nsoset(kind_maxl(ikind))
CALL add_broadened_value(pdos_curve(:, icomp, ikind, ispin), &
emin, de, eigenvalues(imo), &
wkp*proj_weight(icomp, ikind), &
broaden_type, broaden_width, voigt_mixing)
END DO
ELSE
DO icomp = 1, kind_maxl(ikind) + 1
CALL add_broadened_value(pdos_curve(:, icomp, ikind, ispin), &
emin, de, eigenvalues(imo), &
wkp*proj_weight(icomp, ikind), &
broaden_type, broaden_width, voigt_mixing)
END DO
END IF
END DO
END DO
IF (kpoints%use_real_wfn) THEN
CALL cp_fm_release(shalfc)
ELSE
CALL cp_cfm_release(cshalfc)
END IF
END DO
END DO
END IF
CALL para_env%sum(pdos_curve)
IF (append .AND. iterstep > 1) THEN
my_pos = "APPEND"
ELSE
my_pos = "REWIND"
END IF
my_act = "WRITE"
spin(1) = "ALPHA"
spin(2) = "BETA"
DO ikind = 1, nkind
IF (kind_maxl(ikind) < 0) CYCLE
CALL get_atomic_kind(atomic_kind_set(ikind), name=kind_name)
DO ispin = 1, nspins
IF (nspins == 2) THEN
my_mittle = TRIM(spin(ispin))//"_k"//TRIM(ADJUSTL(cp_to_string(ikind)))
ELSE
my_mittle = "k"//TRIM(ADJUSTL(cp_to_string(ikind)))
END IF
CALL write_broadened_pdos_curve(logger, dft_section, TRIM(my_mittle), my_pos, my_act, &
iterstep, e_fermi, &
"K-point projected DOS for atomic kind "//TRIM(kind_name), &
kind_maxl(ikind), separate_components, &
pdos_curve(:, :, ikind, ispin), emin, de, &
broaden_type, broaden_width, voigt_mixing, ndigits)
END DO
END DO
DEALLOCATE (ao_comp, ao_kind, ao_l, ao_weight, kind_maxl, pdos_curve, proj_weight, &
vecbuffer, zvecbuffer)
CALL timestop(handle)
END SUBROUTINE calculate_projected_dos_kp
! **************************************************************************************************
!> \brief Build AO mapping arrays for PDOS accumulation.
!> \param qs_kind_set ...
!> \param particle_set ...
!> \param nao ...
!> \param ao_kind ...
!> \param ao_l ...
!> \param ao_comp ...
! **************************************************************************************************
SUBROUTINE build_pdos_ao_map(qs_kind_set, particle_set, nao, ao_kind, ao_l, ao_comp)
TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
INTEGER, INTENT(IN) :: nao
INTEGER, ALLOCATABLE, DIMENSION(:), INTENT(OUT) :: ao_kind, ao_l, ao_comp
INTEGER :: iatom, ikind, irow, iset, ishell, iso, &
lshell, maxl, nset
INTEGER, DIMENSION(:), POINTER :: nshell
INTEGER, DIMENSION(:, :), POINTER :: l
TYPE(gto_basis_set_type), POINTER :: orb_basis_set
ALLOCATE (ao_kind(nao), ao_l(nao), ao_comp(nao))
irow = 0
DO iatom = 1, SIZE(particle_set)
NULLIFY (orb_basis_set)
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)
CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
nset=nset, &
nshell=nshell, &
l=l, maxl=maxl)
DO iset = 1, nset
DO ishell = 1, nshell(iset)
lshell = l(ishell, iset)
DO iso = 1, nso(lshell)
irow = irow + 1
CPASSERT(irow <= nao)
ao_kind(irow) = ikind
ao_l(irow) = lshell
ao_comp(irow) = nsoset(lshell - 1) + iso
END DO
END DO
END DO
END DO
CPASSERT(irow == nao)
END SUBROUTINE build_pdos_ao_map
! **************************************************************************************************
!> \brief Accumulate AO weights into kind/l or kind/component projected weights.
!> \param ao_weight ...
!> \param ao_kind ...
!> \param ao_l ...
!> \param ao_comp ...
!> \param separate_components ...
!> \param proj_weight ...
! **************************************************************************************************
SUBROUTINE accumulate_pdos_weights(ao_weight, ao_kind, ao_l, ao_comp, &
separate_components, proj_weight)
REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: ao_weight
INTEGER, DIMENSION(:), INTENT(IN) :: ao_kind, ao_l, ao_comp
LOGICAL, INTENT(IN) :: separate_components
REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: proj_weight
INTEGER :: iao, icomp, ikind
DO iao = 1, SIZE(ao_weight)
ikind = ao_kind(iao)
IF (separate_components) THEN
icomp = ao_comp(iao)
ELSE
icomp = ao_l(iao) + 1
END IF
proj_weight(icomp, ikind) = proj_weight(icomp, ikind) + ao_weight(iao)
END DO
END SUBROUTINE accumulate_pdos_weights
! **************************************************************************************************
!> \brief Write a broadened k-point PDOS curve.
!> \param logger ...
!> \param dft_section ...
!> \param middle_name ...
!> \param file_position ...
!> \param file_action ...
!> \param iterstep ...
!> \param e_fermi ...
!> \param title ...
!> \param maxl ...
!> \param separate_components ...
!> \param pdos_curve ...
!> \param emin ...
!> \param de ...
!> \param broaden_type ...
!> \param broaden_width ...
!> \param voigt_mixing ...
!> \param ndigits ...
! **************************************************************************************************
SUBROUTINE write_broadened_pdos_curve(logger, dft_section, middle_name, file_position, &
file_action, iterstep, e_fermi, title, maxl, &
separate_components, pdos_curve, emin, de, &
broaden_type, broaden_width, voigt_mixing, ndigits)
TYPE(cp_logger_type), POINTER :: logger
TYPE(section_vals_type), POINTER :: dft_section
CHARACTER(LEN=*), INTENT(IN) :: middle_name, file_position, file_action
INTEGER, INTENT(IN) :: iterstep
REAL(KIND=dp), INTENT(IN) :: e_fermi
CHARACTER(LEN=*), INTENT(IN) :: title
INTEGER, INTENT(IN) :: maxl
LOGICAL, INTENT(IN) :: separate_components
REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: pdos_curve
REAL(KIND=dp), INTENT(IN) :: emin, de
INTEGER, INTENT(IN) :: broaden_type
REAL(KIND=dp), INTENT(IN) :: broaden_width, voigt_mixing
INTEGER, INTENT(IN) :: ndigits
CHARACTER(LEN=20) :: fmtstr_data
CHARACTER(LEN=6), ALLOCATABLE, DIMENSION(:, :, :) :: tmp_str
INTEGER :: i, icomp, il, im, iw, ncomponents, nhist
REAL(KIND=dp) :: eval
nhist = SIZE(pdos_curve, 1)
IF (separate_components) THEN
ncomponents = nsoset(maxl)
ELSE
ncomponents = maxl + 1
END IF
iw = cp_print_key_unit_nr(logger, dft_section, "PRINT%PDOS", &
extension=".pdos", file_position=file_position, file_action=file_action, &
file_form="FORMATTED", middle_name=TRIM(middle_name))
IF (iw > 0) THEN
WRITE (UNIT=iw, FMT="(A,I0,A,F12.6,A)") &
"# "//TRIM(title)//" at iteration step i = ", iterstep, ", E(Fermi) = ", e_fermi, " a.u."
CALL write_broadening_info(iw, broaden_type, broaden_width, voigt_mixing)
WRITE (UNIT=iw, FMT="(A)", ADVANCE="NO") "# Energy[a.u.]"
IF (separate_components) THEN
ALLOCATE (tmp_str(0:0, 0:maxl, -maxl:maxl))
tmp_str = ""
DO il = 0, maxl
DO im = -il, il
tmp_str(0, il, im) = sgf_symbol(0, il, im)
WRITE (UNIT=iw, FMT="(2X,A)", ADVANCE="NO") TRIM(tmp_str(0, il, im))
END DO
END DO
DEALLOCATE (tmp_str)
ELSE
DO il = 0, maxl
WRITE (UNIT=iw, FMT="(2X,A)", ADVANCE="NO") TRIM(l_sym(il))
END DO
END IF
WRITE (UNIT=iw, FMT="()")
WRITE (UNIT=fmtstr_data, FMT="(A,I0,A)") "(2X,F20.", ndigits, ")"
DO i = 1, nhist
eval = emin + (i - 1)*de
WRITE (UNIT=iw, FMT="(F15.8)", ADVANCE="NO") eval
DO icomp = 1, ncomponents
WRITE (UNIT=iw, FMT=fmtstr_data, ADVANCE="NO") pdos_curve(i, icomp)
END DO
WRITE (UNIT=iw, FMT="()")
END DO
END IF
CALL cp_print_key_finished_output(iw, logger, dft_section, "PRINT%PDOS")
END SUBROUTINE write_broadened_pdos_curve
! **************************************************************************************************
!> \brief Write a broadened PDOS curve for a projected weight matrix.
!> \param logger ...
@ -990,10 +1416,7 @@ CONTAINS
REAL(KIND=dp), INTENT(IN) :: broaden_width, voigt_mixing
INTEGER, INTENT(IN) :: ndigits
CHARACTER(LEN=20) :: fmtstr_data
CHARACTER(LEN=6), ALLOCATABLE, DIMENSION(:, :, :) :: tmp_str
INTEGER :: i, icomp, il, im, imo, iw, ncomponents, &
nhist
INTEGER :: i, icomp, imo, ncomponents, nhist
REAL(KIND=dp) :: cutoff, emax, emin, eval, line_shape
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: pdos_curve
@ -1019,41 +1442,10 @@ CONTAINS
END DO
END DO
iw = cp_print_key_unit_nr(logger, dft_section, "PRINT%PDOS", &
extension=".pdos", file_position=file_position, file_action=file_action, &
file_form="FORMATTED", middle_name=TRIM(middle_name))
IF (iw > 0) THEN
WRITE (UNIT=iw, FMT="(A,I0,A,F12.6,A)") &
"# "//TRIM(title)//" at iteration step i = ", iterstep, ", E(Fermi) = ", e_fermi, " a.u."
CALL write_broadening_info(iw, broaden_type, broaden_width, voigt_mixing)
WRITE (UNIT=iw, FMT="(A)", ADVANCE="NO") "# Energy[a.u.]"
IF (separate_components) THEN
ALLOCATE (tmp_str(0:0, 0:maxl, -maxl:maxl))
tmp_str = ""
DO il = 0, maxl
DO im = -il, il
tmp_str(0, il, im) = sgf_symbol(0, il, im)
WRITE (UNIT=iw, FMT="(2X,A)", ADVANCE="NO") TRIM(tmp_str(0, il, im))
END DO
END DO
DEALLOCATE (tmp_str)
ELSE
DO il = 0, maxl
WRITE (UNIT=iw, FMT="(2X,A)", ADVANCE="NO") TRIM(l_sym(il))
END DO
END IF
WRITE (UNIT=iw, FMT="()")
WRITE (UNIT=fmtstr_data, FMT="(A,I0,A)") "(2X,F20.", ndigits, ")"
DO i = 1, nhist
eval = emin + (i - 1)*de
WRITE (UNIT=iw, FMT="(F15.8)", ADVANCE="NO") eval
DO icomp = 1, ncomponents
WRITE (UNIT=iw, FMT=fmtstr_data, ADVANCE="NO") pdos_curve(i, icomp)
END DO
WRITE (UNIT=iw, FMT="()")
END DO
END IF
CALL cp_print_key_finished_output(iw, logger, dft_section, "PRINT%PDOS")
CALL write_broadened_pdos_curve(logger, dft_section, middle_name, file_position, file_action, &
iterstep, e_fermi, title, maxl, separate_components, &
pdos_curve, emin, de, broaden_type, broaden_width, &
voigt_mixing, ndigits)
DEALLOCATE (pdos_curve)

View file

@ -182,7 +182,8 @@ MODULE qs_scf_post_gpw
neighbor_list_iterator_release,&
neighbor_list_set_p_type
USE qs_ot_eigensolver, ONLY: ot_eigensolver
USE qs_pdos, ONLY: calculate_projected_dos
USE qs_pdos, ONLY: calculate_projected_dos,&
calculate_projected_dos_kp
USE qs_resp, ONLY: resp_fit
USE qs_rho0_types, ONLY: get_rho0_mpole,&
mpole_rho_atom,&
@ -2354,7 +2355,7 @@ CONTAINS
! Print the projected density of states (pDOS) for each atomic kind
IF (do_pdos) THEN
IF (do_kpoints) THEN
CPWARN("Projected density of states (pDOS) is not implemented for k points")
CALL calculate_projected_dos_kp(qs_env, dft_section)
ELSE
CALL get_qs_env(qs_env, &
mos=mos, &

View file

@ -14,8 +14,8 @@
"cc5.inp" = [{matcher="M011", tol=5.0E-10, ref=-45.389916356646268}]
"cc6.inp" = [{matcher="M011", tol=5.0E-10, ref=-45.389916356646268}]
"cc7.inp" = [{matcher="M011", tol=5.0E-10, ref=-45.389916355960565},
{matcher="gext", tol=5e-08, ref=0.0700127407}]
{matcher="gext", tol=5.0E-08, ref=0.0700127407}]
"cc8.inp" = [{matcher="M011", tol=5.0E-10, ref=-45.389916356641265},
{matcher="gext", tol=5e-08, ref=0.7099959941}]
"cc9.inp" = [{matcher="M011", tol=5.0E-10, ref=-45.270019242293365}]
{matcher="gext", tol=5.0E-08, ref=0.7099959941}]
"cc9.inp" = [{matcher="M011", tol=5.0E-10, ref=-45.270016601865670}]
#EOF

View file

@ -9,29 +9,30 @@
BASIS_SET_FILE_NAME GTH_BASIS_SETS
POTENTIAL_FILE_NAME POTENTIAL
&KPOINTS
EPS_GEO 1.e-8
FULL_GRID OFF
PARALLEL_GROUP_SIZE 0
SCHEME MONKHORST-PACK 2 2 2
SYMMETRY ON
VERBOSE T
&END KPOINTS
&MGRID
CUTOFF 280
REL_CUTOFF 50
CUTOFF 250
REL_CUTOFF 40
&END MGRID
&PRINT
&PDOS
COMPONENTS T
NLUMO 5
&END PDOS
&END PRINT
&QS
EPS_DEFAULT 1.0E-12
METHOD GPW
&END QS
&SCF
CHOLESKY OFF
EPS_EIGVAL 1.e-8
EPS_SCF 1.0E-8
EPS_SCF 1.0E-6
MAX_SCF 20
SCF_GUESS ATOMIC
&MIXING
ALPHA 0.50
ALPHA 0.40
METHOD BROYDEN_MIXING
&END MIXING
&PRINT