From 61d4209430e287abcf9f6edf52047f033dc846b6 Mon Sep 17 00:00:00 2001 From: MaggieB092 <76487868+MaggieB092@users.noreply.github.com> Date: Thu, 15 May 2025 13:35:20 +0100 Subject: [PATCH] HP-DFT modules and regtests (#4138) --- src/CMakeLists.txt | 1 + src/cp_control_types.F | 30 +- src/cp_control_utils.F | 86 ++++- src/hairy_probes.F | 621 +++++++++++++++++++++++++++++++++++ src/input_cp2k_dft.F | 71 ++++ src/kpoint_methods.F | 215 +++++++++--- src/qs_mo_methods.F | 21 +- src/qs_mo_occupation.F | 79 ++++- src/qs_scf.F | 33 +- src/qs_scf_diagonalization.F | 33 +- src/qs_scf_initialization.F | 29 +- src/qs_scf_loop_utils.F | 36 +- src/qs_scf_output.F | 6 + src/qs_scf_post_gpw.F | 9 +- tests/HP-DFT/TEST_FILES.toml | 10 + tests/HP-DFT/hp-dft_0eV.inp | 74 +++++ tests/HP-DFT/hp-dft_1eV.inp | 74 +++++ tests/TEST_DIRS | 3 +- tests/matchers.py | 1 + 19 files changed, 1342 insertions(+), 90 deletions(-) create mode 100644 src/hairy_probes.F create mode 100644 tests/HP-DFT/TEST_FILES.toml create mode 100644 tests/HP-DFT/hp-dft_0eV.inp create mode 100644 tests/HP-DFT/hp-dft_1eV.inp diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index f161438162..669ee9afce 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -225,6 +225,7 @@ list( gw_utils.F gw_integrals.F gw_kp_to_real_space_and_back.F + hairy_probes.F hartree_local_methods.F hartree_local_types.F header.F diff --git a/src/cp_control_types.F b/src/cp_control_types.F index ff1a6b23c9..146be5187e 100644 --- a/src/cp_control_types.F +++ b/src/cp_control_types.F @@ -35,6 +35,18 @@ MODULE cp_control_types PRIVATE +!*************************************************************************************************** +!\brief Control parameters for hairy-probes DFT +!*************************************************************************************************** + TYPE hairy_probes_type + REAL(kind=dp) :: alpha !solution probes parameter + REAL(kind=dp) :: mu !chemical potenatial of electrons in reservoir + REAL(kind=dp) :: T !temperature of electrons in reservoir + REAL(kind=dp) :: eps_hp !tolerance for accuracy checks on occupation numbers + INTEGER :: natoms, last_ao, first_ao + INTEGER, DIMENSION(:), POINTER :: atom_ids !atom ids to which the probes are attached + END TYPE hairy_probes_type + ! ************************************************************************************************** ! \brief Control parameters for pw grids ! ************************************************************************************************** @@ -582,6 +594,8 @@ MODULE cp_control_types TYPE(smeagol_control_type), POINTER :: smeagol_control => NULL() TYPE(efield_p_type), POINTER, & DIMENSION(:) :: efield_fields => NULL() + TYPE(hairy_probes_type), POINTER, & + DIMENSION(:) :: probe => NULL() INTEGER :: nspins = 0, & charge = 0, & multiplicity = 0, & @@ -631,7 +645,8 @@ MODULE cp_control_types correct_el_density_dip = .FALSE., & do_sccs = .FALSE., & apply_embed_pot = .FALSE., & - apply_dmfet_pot = .FALSE. + apply_dmfet_pot = .FALSE., & + hairy_probes = .FALSE. END TYPE dft_control_type CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'cp_control_types' @@ -656,7 +671,8 @@ MODULE cp_control_types rtp_control_type, & sccs_control_type, & stda_control_type, & - smeared_type + smeared_type, & + hairy_probes_type ! Public subroutines @@ -788,9 +804,11 @@ CONTAINS NULLIFY (dft_control%smeagol_control) NULLIFY (dft_control%rtp_control) NULLIFY (dft_control%sccs_control) + NULLIFY (dft_control%probe) dft_control%do_sccs = .FALSE. dft_control%apply_embed_pot = .FALSE. dft_control%apply_dmfet_pot = .FALSE. + dft_control%hairy_probes = .FALSE. CALL qs_control_create(dft_control%qs_control) CALL tddfpt2_control_create(dft_control%tddfpt2_control) CALL smeagol_control_create(dft_control%smeagol_control) @@ -806,6 +824,8 @@ CONTAINS SUBROUTINE dft_control_release(dft_control) TYPE(dft_control_type), INTENT(INOUT) :: dft_control + INTEGER :: i + CALL qs_control_release(dft_control%qs_control) CALL tddfpt2_control_release(dft_control%tddfpt2_control) IF (ASSOCIATED(dft_control%xas_control)) THEN @@ -817,6 +837,12 @@ CONTAINS CALL maxwell_control_release(dft_control%maxwell_control) CALL smeagol_control_release(dft_control%smeagol_control) CALL efield_fields_release(dft_control%efield_fields) + IF (ASSOCIATED(dft_control%probe)) THEN + DO i = 1, SIZE(dft_control%probe) + DEALLOCATE (dft_control%probe(i)%atom_ids) + END DO + DEALLOCATE (dft_control%probe) + END IF IF (ASSOCIATED(dft_control%sccs_control)) DEALLOCATE (dft_control%sccs_control) IF (ASSOCIATED(dft_control%period_efield)) THEN DEALLOCATE (dft_control%period_efield) diff --git a/src/cp_control_utils.F b/src/cp_control_utils.F index 8b8b518a07..8d1a6ffb9c 100644 --- a/src/cp_control_utils.F +++ b/src/cp_control_utils.F @@ -122,9 +122,8 @@ CONTAINS REAL(KIND=dp), DIMENSION(:), POINTER :: pol TYPE(cp_logger_type), POINTER :: logger TYPE(cp_parser_type) :: parser - TYPE(section_vals_type), POINTER :: hfx_section, maxwell_section, & - sccs_section, scf_section, & - tmp_section, xc_fun_section, xc_section + TYPE(section_vals_type), POINTER :: hairy_probes_section, hfx_section, maxwell_section, & + sccs_section, scf_section, tmp_section, xc_fun_section, xc_section was_present = .FALSE. @@ -396,6 +395,16 @@ CONTAINS "unrestricted Kohn-Sham (UKS) calculations") END IF + !Read the HAIR PROBES input section if present + hairy_probes_section => section_vals_get_subs_vals(dft_section, "HAIRY_PROBES") + CALL section_vals_get(hairy_probes_section, n_repetition=nrep, explicit=is_present) + + IF (is_present) THEN + dft_control%hairy_probes = .TRUE. + ALLOCATE (dft_control%probe(nrep)) + CALL read_hairy_probes_sections(dft_control, hairy_probes_section) + END IF + ! check for the presence of the low spin roks section tmp_section => section_vals_get_subs_vals(dft_section, "LOW_SPIN_ROKS") CALL section_vals_get(tmp_section, explicit=dft_control%low_spin_roks) @@ -1843,7 +1852,7 @@ CONTAINS CHARACTER(len=*), PARAMETER :: routineN = 'write_dft_control' CHARACTER(LEN=20) :: tmpStr - INTEGER :: handle, i, output_unit + INTEGER :: handle, i, i_rep, n_rep, output_unit REAL(kind=dp) :: density_cut, density_smooth_cut_range, & gradient_cut, tau_cut TYPE(cp_logger_type), POINTER :: logger @@ -2084,6 +2093,25 @@ CONTAINS END IF + IF (dft_control%hairy_probes .EQV. .TRUE.) THEN + n_rep = SIZE(dft_control%probe) + IF (output_unit > 0) THEN + DO i_rep = 1, n_rep + WRITE (UNIT=output_unit, FMT="(T2,A,I5)") & + "HP | hair probe set", i_rep + WRITE (UNIT=output_unit, FMT="(T2,A,T61,*(I5))") & + "HP| atom indexes", & + (dft_control%probe(i_rep)%atom_ids(i), i=1, dft_control%probe(i_rep)%natoms) + WRITE (UNIT=output_unit, FMT="(T2,A,T61,ES20.6)") & + "HP| potential", dft_control%probe(i_rep)%mu + WRITE (UNIT=output_unit, FMT="(T2,A,T61,F20.2)") & + "HP| temperature", dft_control%probe(i_rep)%T + WRITE (UNIT=output_unit, FMT="(T2,A,T61,ES20.6)") & + "HP| eps_hp", dft_control%probe(i_rep)%eps_hp + END DO + END IF + END IF + CALL cp_print_key_finished_output(output_unit, logger, dft_section, & "PRINT%DFT_CONTROL_PARAMETERS") @@ -2816,4 +2844,54 @@ CONTAINS END SUBROUTINE read_admm_block_list +! ************************************************************************************************** +!> \brief ... +!> \param dft_control ... +!> \param hairy_probes_section ... +!> \param +!> \param +! ************************************************************************************************** + SUBROUTINE read_hairy_probes_sections(dft_control, hairy_probes_section) + TYPE(dft_control_type), POINTER :: dft_control + TYPE(section_vals_type), POINTER :: hairy_probes_section + + INTEGER :: i, j, jj, kk, n_rep + INTEGER, DIMENSION(:), POINTER :: tmplist + + DO i = 1, SIZE(dft_control%probe) + NULLIFY (dft_control%probe(i)%atom_ids) + + CALL section_vals_val_get(hairy_probes_section, "ATOM_IDS", i_rep_section=i, n_rep_val=n_rep) + jj = 0 + DO kk = 1, n_rep + CALL section_vals_val_get(hairy_probes_section, "ATOM_IDS", i_rep_section=i, i_rep_val=kk, i_vals=tmplist) + jj = jj + SIZE(tmplist) + END DO + + dft_control%probe(i)%natoms = jj + IF (dft_control%probe(i)%natoms < 1) & + CPABORT("Need at least 1 atom to use hair probes formalism") + ALLOCATE (dft_control%probe(i)%atom_ids(dft_control%probe(i)%natoms)) + + jj = 0 + DO kk = 1, n_rep + CALL section_vals_val_get(hairy_probes_section, "ATOM_IDS", i_rep_section=i, i_rep_val=kk, i_vals=tmplist) + DO j = 1, SIZE(tmplist) + jj = jj + 1 + dft_control%probe(i)%atom_ids(jj) = tmplist(j) + END DO + END DO + + CALL section_vals_val_get(hairy_probes_section, "MU", i_rep_section=i, r_val=dft_control%probe(i)%mu) + + CALL section_vals_val_get(hairy_probes_section, "T", i_rep_section=i, r_val=dft_control%probe(i)%T) + + CALL section_vals_val_get(hairy_probes_section, "ALPHA", i_rep_section=i, r_val=dft_control%probe(i)%alpha) + + CALL section_vals_val_get(hairy_probes_section, "eps_hp", i_rep_section=i, r_val=dft_control%probe(i)%eps_hp) + END DO + + END SUBROUTINE read_hairy_probes_sections +! ************************************************************************************************** + END MODULE cp_control_utils diff --git a/src/hairy_probes.F b/src/hairy_probes.F new file mode 100644 index 0000000000..9d9d6af9bb --- /dev/null +++ b/src/hairy_probes.F @@ -0,0 +1,621 @@ +!--------------------------------------------------------------------------------------------------! +! CP2K: A general program to perform molecular dynamics simulations ! +! Copyright 2000-2025 CP2K developers group ! +! ! +! SPDX-License-Identifier: GPL-2.0-or-later ! +!--------------------------------------------------------------------------------------------------! + +MODULE hairy_probes + + 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_control_types, ONLY: hairy_probes_type + USE cp_fm_types, ONLY: cp_fm_get_info,& + cp_fm_get_submatrix,& + cp_fm_type + USE kahan_sum, ONLY: accurate_sum + USE kinds, ONLY: dp + USE orbital_pointers, ONLY: nso + USE particle_types, ONLY: particle_type + USE qs_kind_types, ONLY: get_qs_kind,& + get_qs_kind_set,& + qs_kind_type +#include "./base/base_uses.f90" + + IMPLICIT NONE + PRIVATE + + CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'hairy_probes' + + INTEGER, PARAMETER, PRIVATE :: BISECT_MAX_ITER = 400 + PUBLIC :: probe_occupancy, probe_occupancy_kp, AO_boundaries + +CONTAINS + +!************************************************************************************** +!> \brief subroutine to calculate occupation number and 'Fermi' level using the +!> \brief HAIR PROBE approach; gamma point calculation. +!> \param occ occupation numbers +!> \param fermi fermi level +!> \param kTS entropic energy contribution +!> \param energies MOs eigenvalues +!> \param coeff MOs coefficient +!> \param maxocc maximum allowed occupation number of an MO (1 or 2) +!> \param probe hairy probe +!> \param N number of electrons +! ************************************************************************************************** + + SUBROUTINE probe_occupancy(occ, fermi, kTS, energies, coeff, maxocc, probe, N) + + !i/o variables and arrays + REAL(KIND=dp), INTENT(out) :: occ(:), fermi, kTS + REAL(KIND=dp), INTENT(IN) :: energies(:) + TYPE(cp_fm_type), INTENT(IN), POINTER :: coeff + REAL(KIND=dp), INTENT(IN) :: maxocc + TYPE(hairy_probes_type), INTENT(INOUT) :: probe(:) + REAL(KIND=dp), INTENT(IN) :: N + + REAL(KIND=dp), PARAMETER :: epsocc = 1.0e-12_dp + + INTEGER :: iter, ncol_global, nrow_global + REAL(KIND=dp) :: de, delta_fermi, fermi_fit, fermi_half, & + fermi_max, fermi_min, h, N_fit, & + N_half, N_max, N_min, N_now, y0, y1, y2 + REAL(KIND=dp), ALLOCATABLE :: smatrix_squared(:, :) + REAL(KIND=dp), POINTER :: smatrix(:, :) + +!subroutine variables and arrays +!smatrix: Spherical MOs matrix +!squared Spherical MOs matrix + + CALL cp_fm_get_info(coeff, & + nrow_global=nrow_global, & + ncol_global=ncol_global) + ALLOCATE (smatrix(nrow_global, ncol_global)) + CALL cp_fm_get_submatrix(coeff, smatrix) + + ALLOCATE (smatrix_squared(nrow_global, ncol_global)) + smatrix_squared(:, :) = smatrix(:, :)**2.0d0 + +!************************************************************************************** +!1)calculate Fermi energy using the hairy probe formula +!************************************************************************************** + de = probe(1)%T*LOG((1.0_dp - epsocc)/epsocc) + de = MAX(de, 0.5_dp) + fermi_max = MAXVAL(probe%mu) + de + fermi_min = MINVAL(probe%mu) - de + + CALL HP_occupancy(probe=probe, matrix=smatrix_squared, energies=energies, & + maxocc=maxocc, fermi=fermi_max, occupancy=occ, kTS=kTS) + N_max = accurate_sum(occ) + + CALL HP_occupancy(probe=probe, matrix=smatrix_squared, energies=energies, & + maxocc=maxocc, fermi=fermi_min, occupancy=occ, kTS=kTS) + N_min = accurate_sum(occ) + + iter = 0 + DO WHILE (ABS(N_max - N_min) > N*epsocc) + iter = iter + 1 + + fermi_half = (fermi_max + fermi_min)/2.0_dp + CALL HP_occupancy(probe=probe, matrix=smatrix_squared, energies=energies, & + maxocc=maxocc, fermi=fermi_half, occupancy=occ, kTS=kTS) + N_half = accurate_sum(occ) + + h = fermi_half - fermi_min + IF (h .GT. N*epsocc*100) THEN + y0 = N_min - N + y1 = N_half - N + y2 = N_max - N + + CALL three_point_zero(y0, y1, y2, h, delta_fermi) + fermi_fit = fermi_min + delta_fermi + + CALL HP_occupancy(probe=probe, matrix=smatrix_squared, energies=energies, & + maxocc=maxocc, fermi=fermi_fit, occupancy=occ, kTS=kTS) + N_fit = accurate_sum(occ) + END IF + + !define 1st bracked using fermi_half + IF (N_half < N) THEN + fermi_min = fermi_half + N_min = N_half + ELSE IF (N_half > N) THEN + fermi_max = fermi_half + N_max = N_half + ELSE + fermi_min = fermi_half + N_min = N_half + fermi_max = fermi_half + N_max = N_half + h = 0.0d0 + END IF + + !define 2nd bracker using fermi_fit + IF (h .GT. N*epsocc*100) THEN + IF (fermi_fit .GE. fermi_min .AND. fermi_fit .LE. fermi_max) THEN + IF (N_fit < N) THEN + fermi_min = fermi_fit + N_min = N_fit + ELSE IF (N_fit > N) THEN + fermi_max = fermi_fit + N_max = N_fit + END IF + END IF + END IF + + IF (ABS(N_max - N) < N*epsocc) THEN + fermi = fermi_max + EXIT + ELSE IF (ABS(N_min - N) < N*epsocc) THEN + fermi = fermi_min + EXIT + END IF + + IF (iter > BISECT_MAX_ITER) THEN + CPWARN("Maximum number of iterations reached while finding the Fermi energy") + EXIT + END IF + END DO + +!************************************************************************************** +!2)calculate occupation numbers according to hairy probe formula +!************************************************************************************** + occ(:) = 0.0_dp + N_now = 0.0_dp + CALL HP_occupancy(probe=probe, matrix=smatrix_squared, energies=energies, & + maxocc=maxocc, fermi=fermi, occupancy=occ, kTS=kTS) + N_now = accurate_sum(occ) + + IF (ABS(N_now - N) > N*epsocc) CPWARN("Total number of electrons is not accurate - HP") + + DEALLOCATE (smatrix, smatrix_squared) + + END SUBROUTINE probe_occupancy + +!************************************************************************************** +!> \brief subroutine to calculate occupation number and 'Fermi' level using the +!> \brief HAIR PROBE approach; kpoints calculation. +!> \param occ occupation numbers +!> \param fermi fermi level +!> \param kTS entropic energy contribution +!> \param energies eigenvalues +!> \param rcoeff ... +!> \param icoeff ... +!> \param maxocc maximum allowed occupation number of an MO (1 or 2) +!> \param probe hairy probe +!> \param N number of electrons +!> \param wk weight of kpoints +! ************************************************************************************************** + SUBROUTINE probe_occupancy_kp(occ, fermi, kTS, energies, rcoeff, icoeff, maxocc, probe, N, wk) + + REAL(KIND=dp), INTENT(OUT) :: occ(:, :, :), fermi, kTS + REAL(KIND=dp), INTENT(IN) :: energies(:, :, :), rcoeff(:, :, :, :), & + icoeff(:, :, :, :), maxocc + TYPE(hairy_probes_type), INTENT(IN) :: probe(:) + REAL(KIND=dp), INTENT(IN) :: N, wk(:) + + CHARACTER(LEN=*), PARAMETER :: routineN = 'probe_occupancy_kp' + REAL(KIND=dp), PARAMETER :: epsocc = 1.0e-12_dp + + INTEGER :: handle, ikp, ispin, iter, nAO, nkp, nMO, & + nspin + REAL(KIND=dp) :: de, delta_fermi, fermi_fit, fermi_half, & + fermi_max, fermi_min, h, kTS_kp, & + N_fit, N_half, N_max, N_min, N_now, & + y0, y1, y2 + REAL(KIND=dp), ALLOCATABLE :: coeff_squared(:, :, :, :) + + CALL timeset(routineN, handle) + +!************************************************************************************** +!1)calculate Fermi energy using the hairy probe formula +!************************************************************************************** + nAO = SIZE(rcoeff, 1) + nMO = SIZE(rcoeff, 2) + nkp = SIZE(rcoeff, 3) + nspin = SIZE(rcoeff, 4) + + ALLOCATE (coeff_squared(nAO, nMO, nkp, nspin)) + coeff_squared(:, :, :, :) = rcoeff(:, :, :, :)**2.0d0 + icoeff(:, :, :, :)**2.0d0 + + occ(:, :, :) = 0.0_dp + + !define initial brackets + de = probe(1)%T*LOG((1.0_dp - epsocc)/epsocc) + de = MAX(de, 0.5_dp) + fermi_max = MAXVAL(probe%mu) + de + fermi_min = MINVAL(probe%mu) - de + + N_max = 0.0_dp + kTS = 0.0_dp + !***HP loop + DO ispin = 1, nspin + DO ikp = 1, nkp + CALL HP_occupancy(probe, coeff_squared(:, :, ikp, ispin), energies(:, ikp, ispin), & + maxocc, fermi_max, occ(:, ikp, ispin), kTS_kp) + + kTS = kTS + kTS_kp*wk(ikp) !entropic contribution + N_max = N_max + accurate_sum(occ(1:nMo, ikp, ispin))*wk(ikp) + END DO + END DO + !***HP loop + + N_min = 0.0_dp + kTS = 0.0_dp + !***HP loop + DO ispin = 1, nspin + DO ikp = 1, nkp + CALL HP_occupancy(probe, coeff_squared(:, :, ikp, ispin), energies(:, ikp, ispin), & + maxocc, fermi_min, occ(:, ikp, ispin), kTS_kp) + + kTS = kTS + kTS_kp*wk(ikp) !entropic contribution + N_min = N_min + accurate_sum(occ(1:nMo, ikp, ispin))*wk(ikp) + END DO + END DO + !***HP loop + + iter = 0 + DO WHILE (ABS(N_max - N_min) > N*epsocc) + iter = iter + 1 + fermi_half = (fermi_max + fermi_min)/2.0_dp + N_half = 0.0_dp + kTS = 0.0_dp + !***HP loop + DO ispin = 1, nspin + DO ikp = 1, nkp + CALL HP_occupancy(probe, coeff_squared(:, :, ikp, ispin), energies(:, ikp, ispin), & + maxocc, fermi_half, occ(:, ikp, ispin), kTS_kp) + + kTS = kTS + kTS_kp*wk(ikp) !entropic contribution + N_half = N_half + accurate_sum(occ(1:nMo, ikp, ispin))*wk(ikp) + END DO + END DO + !***HP loop + + h = fermi_half - fermi_min + IF (h .GT. N*epsocc*100) THEN + y0 = N_min - N + y1 = N_half - N + y2 = N_max - N + + CALL three_point_zero(y0, y1, y2, h, delta_fermi) + fermi_fit = fermi_min + delta_fermi + N_fit = 0.0_dp + kTS = 0.0_dp + + !***HP loop + DO ispin = 1, nspin + DO ikp = 1, nkp + CALL HP_occupancy(probe, coeff_squared(:, :, ikp, ispin), energies(:, ikp, ispin), & + maxocc, fermi_fit, occ(:, ikp, ispin), kTS_kp) + + kTS = kTS + kTS_kp*wk(ikp) !entropic contribution + N_fit = N_fit + accurate_sum(occ(1:nMo, ikp, ispin))*wk(ikp) + END DO + END DO + !***HP loop + + END IF + + !define 1st bracked using fermi_half + IF (N_half < N) THEN + fermi_min = fermi_half + N_min = N_half + ELSE IF (N_half > N) THEN + fermi_max = fermi_half + N_max = N_half + ELSE + fermi_min = fermi_half + N_min = N_half + fermi_max = fermi_half + N_max = N_half + h = 0.0d0 + END IF + + !define 2nd bracker using fermi_fit + IF (h .GT. N*epsocc*100) THEN + IF (fermi_fit .GE. fermi_min .AND. fermi_fit .LE. fermi_max) THEN + IF (N_fit < N) THEN + fermi_min = fermi_fit + N_min = N_fit + ELSE IF (N_fit > N) THEN + fermi_max = fermi_fit + N_max = N_fit + END IF + END IF + END IF + + IF (ABS(N_max - N) < N*epsocc) THEN + fermi = fermi_max + EXIT + ELSE IF (ABS(N_min - N) < N*epsocc) THEN + fermi = fermi_min + EXIT + END IF + + IF (iter > BISECT_MAX_ITER) THEN + CPWARN("Maximum number of iterations reached while finding the Fermi energy") + EXIT + END IF + END DO + +!************************************************************************************** +!3)calculate occupation numbers using the hairy probe formula +!************************************************************************************** + N_now = 0.0_dp + kTS = 0.0_dp + + !***HP loop + DO ispin = 1, nspin + DO ikp = 1, nkp + CALL HP_occupancy(probe, coeff_squared(:, :, ikp, ispin), energies(:, ikp, ispin), & + maxocc, fermi, occ(:, ikp, ispin), kTS_kp) + + kTS = kTS + kTS_kp*wk(ikp) !entropic contribution + N_now = N_now + accurate_sum(occ(1:nMo, ikp, ispin))*wk(ikp) + END DO + END DO + !***HP loop + + DEALLOCATE (coeff_squared) + + CALL timestop(handle) + END SUBROUTINE probe_occupancy_kp + +!************************************************************************************** +! +! +! +!************************************************************************************** +! ************************************************************************************************** +!> \brief ... +!> \param probe ... +!> \param matrix ... +!> \param energies ... +!> \param maxocc ... +!> \param fermi ... +!> \param occupancy ... +!> \param kTS ... +! ************************************************************************************************** + SUBROUTINE HP_occupancy(probe, matrix, energies, maxocc, fermi, occupancy, kTS) + + TYPE(hairy_probes_type), INTENT(IN) :: probe(:) + REAL(KIND=dp), INTENT(IN) :: matrix(:, :), energies(:), maxocc, fermi + REAL(KIND=dp), INTENT(OUT) :: occupancy(:), kTS + + INTEGER :: iMO, ip, nMO, np + REAL(KIND=dp) :: alpha, C, f, fermi_fun, fermi_fun_sol, & + mu, s, sum_coeff, sum_coeff_sol, & + sum_fermi_fun, sum_fermi_fun_sol + +!squared coefficient matrix + + nMO = SIZE(matrix, 2) + np = SIZE(probe) + kTS = 0.0_dp + alpha = 1.0_dp + ! kTS is the entropic contribution to the electronic energy + + mos2: DO iMO = 1, nMO + + sum_fermi_fun = 0.0_dp + sum_coeff = 0.0_dp + + sum_fermi_fun_sol = 0.0_dp + sum_coeff_sol = 0.0_dp + fermi_fun_sol = 0.0_dp + + probes2: DO ip = 1, np + IF (probe(ip)%alpha .LT. 1.0_dp) THEN + alpha = probe(ip)%alpha + !fermi distribution, solution probes + CALL fermi_distribution(energies(iMO), fermi, probe(ip)%T, fermi_fun_sol) + !sum of coefficients, solution probes + sum_coeff_sol = sum_coeff_sol + SUM(matrix(probe(ip)%first_ao:probe(ip)%last_ao, iMO)) + ELSE + C = SUM(matrix(probe(ip)%first_ao:probe(ip)%last_ao, iMO)) + !bias probes + mu = fermi - probe(ip)%mu + !fermi distribution, main probes + CALL fermi_distribution(energies(iMO), mu, probe(ip)%T, fermi_fun) + !sum fermi distribution * coefficients + sum_fermi_fun = sum_fermi_fun + (fermi_fun*C) + !sum cofficients + sum_coeff = sum_coeff + C + END IF + END DO probes2 + + sum_fermi_fun_sol = alpha*fermi_fun_sol*sum_coeff_sol + sum_coeff_sol = alpha*sum_coeff_sol + f = (sum_fermi_fun_sol + sum_fermi_fun)/(sum_coeff_sol + sum_coeff) + occupancy(iMO) = f*maxocc + + !entropy kTS= kT*[f ln f + (1-f) ln (1-f)] + IF (f .EQ. 0.0d0 .OR. f .EQ. 1.0d0) THEN + s = 0.0d0 + ELSE + s = f*LOG(f) + (1.0d0 - f)*LOG(1.0d0 - f) + END IF + kTS = kTS + probe(np)%T*maxocc*s + END DO MOs2 + + END SUBROUTINE HP_occupancy + +!************************************************************************************** +! +! +! +!************************************************************************************** +! ************************************************************************************************** +!> \brief ... +!> \param probe ... +!> \param atomic_kind_set ... +!> \param qs_kind_set ... +!> \param particle_set ... +!> \param nAO ... +! ************************************************************************************************** + SUBROUTINE AO_boundaries(probe, atomic_kind_set, qs_kind_set, particle_set, nAO) + + TYPE(hairy_probes_type), INTENT(INOUT) :: probe + TYPE(atomic_kind_type), INTENT(IN), POINTER :: atomic_kind_set(:) + TYPE(qs_kind_type), INTENT(IN), POINTER :: qs_kind_set(:) + TYPE(particle_type), INTENT(IN), POINTER :: particle_set(:) + INTEGER, INTENT(IN) :: nAO + + INTEGER :: iatom, ii, ikind, iset, isgf, ishell, & + lshell, natom, nset, nsgf, p_atom + INTEGER, DIMENSION(:), POINTER :: nshell + INTEGER, DIMENSION(:, :), POINTER :: l + TYPE(gto_basis_set_type), POINTER :: orb_basis_set + +!from get_atomic_kind_set +!from get_atomic_kind +!from get_qs_kind_set +!subroutine variables and arrays +!from get_qs_kind +!from get_gto_basis_set +!from get_gto_basis_set +!from get_gto_basis_set + + !get sets from different modules: + CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, natom=natom) + CALL get_qs_kind_set(qs_kind_set=qs_kind_set, nsgf=nsgf) !indexes for orbital symbols + + !get iAO boundaries: + p_atom = SIZE(probe%atom_ids) + probe%first_ao = nsgf !nAO + probe%last_ao = 1 + isgf = 0 + + atoms: DO iatom = 1, natom + NULLIFY (orb_basis_set) + !1) iatom is used to find the correct atom kind and its index (ikind) is the particle set + CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind) + + !2) ikind is used to find the basis set associate to that atomic kind + CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set) + + !3) orb_basis_set is used to get the gto basis set variables + 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 ) + END IF + + !4) get iAO boundaries + sets: DO iset = 1, nset + + shells: DO ishell = 1, nshell(iset) + lshell = l(ishell, iset) + + isgf = isgf + nso(lshell) + + boundaries: DO ii = 1, p_atom + IF (iatom .NE. probe%atom_ids(ii)) THEN + CYCLE boundaries + ELSE + !defines iAO boundaries*********** + probe%first_ao = MIN(probe%first_ao, isgf) + probe%last_ao = MAX(probe%last_ao, isgf) + END IF + END DO boundaries + + END DO shells + END DO sets + END DO atoms + + IF (isgf .NE. nAO) CPWARN("row count does not correspond to nAO, number of rows in mo_coeff") + + END SUBROUTINE AO_boundaries + +!******************************************************************************************************* +!******************************************************************************************************* + +! ************************************************************************************************** +!> \brief ... +!> \param E ... +!> \param pot ... +!> \param temp ... +!> \param f_p ... +! ************************************************************************************************** + SUBROUTINE fermi_distribution(E, pot, temp, f_p) + + REAL(kind=dp), INTENT(IN) :: E, pot, temp + REAL(kind=dp), INTENT(OUT) :: f_p + + REAL(KIND=dp) :: arg, exponential, exponential_plus_1, f, & + one_minus_f + + ! have the result of exp go to zero instead of overflowing + IF (E > pot) THEN + arg = -(E - pot)/temp + ! exponential is smaller than 1 + exponential = EXP(arg) + exponential_plus_1 = exponential + 1.0_dp + + one_minus_f = exponential/exponential_plus_1 + f = 1.0_dp/exponential_plus_1 + f_p = one_minus_f + ELSE + arg = (E - pot)/temp + ! exponential is smaller than 1 + exponential = EXP(arg) + exponential_plus_1 = exponential + 1.0_dp + + f = 1.0_dp/exponential_plus_1 + one_minus_f = exponential/exponential_plus_1 + f_p = f + END IF + + END SUBROUTINE fermi_distribution + +!******************************************************************************************************* +!******************************************************************************************************* +! ************************************************************************************************** +!> \brief ... +!> \param y0 ... +!> \param y1 ... +!> \param y2 ... +!> \param h ... +!> \param delta ... +! ************************************************************************************************** + SUBROUTINE three_point_zero(y0, y1, y2, h, delta) + + REAL(kind=dp), INTENT(IN) :: y0, y1, y2, h + REAL(kind=dp), INTENT(OUT) :: delta + + REAL(kind=dp) :: a, b, c, d, u0 + + a = (y2 - 2.0_dp*y1 + y0)/(2.0_dp*h*h) + b = (4.0_dp*y1 - 3.0_dp*y0 - y2)/(2.0_dp*h) + c = y0 + + IF (ABS(a) .LT. 1.0e-15_dp) THEN + delta = 0.0_dp + RETURN + END IF + + d = SQRT(b*b - 4.0_dp*a*c) + + u0 = (-b + d)/(2.0_dp*a) + + IF (u0 .GE. 0.0_dp .AND. u0 .LE. 2.0_dp*h) THEN + delta = u0 + !RETURN + ELSE + u0 = (-b - d)/(2.0_dp*a) + IF (u0 .GE. 0.0_dp .AND. u0 .LE. 2.0_dp*h) THEN + delta = u0 + !RETURN + ELSE + IF (y1 .LT. 0.0_dp) delta = 2.0_dp*h + IF (y1 .GE. 0.0_dp) delta = 0.0_dp + END IF + END IF + + END SUBROUTINE three_point_zero + +END MODULE hairy_probes diff --git a/src/input_cp2k_dft.F b/src/input_cp2k_dft.F index 4e71f32c93..bbe6ae6df1 100644 --- a/src/input_cp2k_dft.F +++ b/src/input_cp2k_dft.F @@ -445,8 +445,79 @@ CONTAINS CALL section_add_subsection(section, subsection) CALL section_release(subsection) + CALL create_hairy_probes_section(subsection) + CALL section_add_subsection(section, subsection) + CALL section_release(subsection) + END SUBROUTINE create_dft_section +! ************************************************************************************************** +!> \brief Hairy Probe DFT Model +!> \param section ... +!> \author Margherita Buraschi +! ************************************************************************************************** + + SUBROUTINE create_hairy_probes_section(section) + TYPE(section_type), POINTER :: section + + TYPE(keyword_type), POINTER :: keyword + + NULLIFY (keyword) + CPASSERT(.NOT. ASSOCIATED(section)) + CALL section_create(section, __LOCATION__, & + name="HAIRY_PROBES", & + description="Sets up a Hairy Probe calculation. ", & + n_keywords=0, n_subsections=0, repeats=.TRUE.) + + CALL keyword_create(keyword, __LOCATION__, & + name="_SECTION_PARAMETERS_", & + description="Controls the activation of hairy probe", & + usage="&HAIRY_PROBES ON", & + default_l_val=.FALSE., & + lone_keyword_l_val=.TRUE.) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, __LOCATION__, name="ATOM_IDS", & + description="Indexes of the atoms to which the probes are attached.", & + usage="ATOM_IDS .. ", & + type_of_var=integer_t, n_var=-1) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, __LOCATION__, name="T", & + description="Electronic temperature [K]", & + usage="T ", & + default_r_val=cp_unit_to_cp2k(value=300.0_dp, unit_str="K"), & + unit_str="K") + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, __LOCATION__, name="MU", & + description="Chemical potential of the electrons in the probes [eV] ", & + usage="MU ", & + default_r_val=cp_unit_to_cp2k(value=0.0_dp, unit_str="eV"), & + unit_str="eV") + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, __LOCATION__, name="ALPHA", & + description="Parameter for solution probes ", & + usage="ALPHA ", & + default_r_val=1.0_dp) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, __LOCATION__, name="EPS_HP", & + description=" Tolerance for accuracy checks on occupation numbers "// & + "calculated using hair-probes. ", & + usage="EPS_HP ", & + default_r_val=1.0E-5_dp) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + END SUBROUTINE create_hairy_probes_section +!#################################################################################### + ! ************************************************************************************************** !> \brief Implicit Solvation Model !> \param section ... diff --git a/src/kpoint_methods.F b/src/kpoint_methods.F index 13ffc5711f..50b0acc2c9 100644 --- a/src/kpoint_methods.F +++ b/src/kpoint_methods.F @@ -18,7 +18,8 @@ MODULE kpoint_methods real_to_scaled USE cp_blacs_env, ONLY: cp_blacs_env_create,& cp_blacs_env_type - USE cp_control_types, ONLY: dft_control_type + USE cp_control_types, ONLY: dft_control_type,& + hairy_probes_type USE cp_dbcsr_api, ONLY: & dbcsr_copy, dbcsr_create, dbcsr_deallocate_matrix, dbcsr_distribution_get, & dbcsr_distribution_type, dbcsr_get_block_p, dbcsr_get_info, dbcsr_iterator_blocks_left, & @@ -34,7 +35,8 @@ MODULE kpoint_methods USE cp_fm_struct, ONLY: cp_fm_struct_type USE cp_fm_types, ONLY: & copy_info_type, cp_fm_cleanup_copy_general, cp_fm_create, cp_fm_finish_copy_general, & - cp_fm_get_info, cp_fm_release, cp_fm_start_copy_general, cp_fm_to_fm, cp_fm_type + cp_fm_get_info, cp_fm_get_submatrix, cp_fm_release, cp_fm_start_copy_general, cp_fm_to_fm, & + cp_fm_type USE cp_log_handling, ONLY: cp_logger_get_default_io_unit USE cryssym, ONLY: crys_sym_gen,& csym_type,& @@ -44,6 +46,7 @@ MODULE kpoint_methods release_csym_type USE fermi_utils, ONLY: fermikp,& fermikp2 + USE hairy_probes, ONLY: probe_occupancy_kp USE input_constants, ONLY: smear_fermi_dirac USE kinds, ONLY: dp USE kpoint_types, ONLY: get_kpoint_info,& @@ -903,20 +906,27 @@ CONTAINS !> \brief Given the eigenvalues of all kpoints, calculates the occupation numbers !> \param kpoint Kpoint environment !> \param smear Smearing information +!> \param probe ... ! ************************************************************************************************** - SUBROUTINE kpoint_set_mo_occupation(kpoint, smear) + SUBROUTINE kpoint_set_mo_occupation(kpoint, smear, probe) TYPE(kpoint_type), POINTER :: kpoint TYPE(smear_type), POINTER :: smear + TYPE(hairy_probes_type), DIMENSION(:), OPTIONAL, & + POINTER :: probe CHARACTER(LEN=*), PARAMETER :: routineN = 'kpoint_set_mo_occupation' - INTEGER :: handle, ik, ikpgr, ispin, kplocal, nb, & - ne_a, ne_b, nelectron, nkp, nmo, nspin + INTEGER :: handle, ik, ikpgr, ispin, kplocal, nao, & + nb, ncol_global, ne_a, ne_b, & + nelectron, nkp, nmo, nrow_global, nspin INTEGER, DIMENSION(2) :: kp_range REAL(KIND=dp) :: kTS, mu, mus(2), nel + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: smatrix REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: weig, wocc + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: icoeff, rcoeff REAL(KIND=dp), DIMENSION(:), POINTER :: eigenvalues, occupation, wkp + TYPE(cp_fm_type), POINTER :: mo_coeff TYPE(kpoint_env_type), POINTER :: kp TYPE(mo_set_type), POINTER :: mo_set TYPE(mp_para_env_type), POINTER :: para_env_inter_kp @@ -928,7 +938,7 @@ CONTAINS kp => kpoint%kp_env(1)%kpoint_env nspin = SIZE(kp%mos, 2) mo_set => kp%mos(1, 1) - CALL get_mo_set(mo_set, nmo=nmo, nelectron=nelectron) + CALL get_mo_set(mo_set, nmo=nmo, nao=nao, nelectron=nelectron) ne_a = nelectron IF (nspin == 2) THEN CALL get_mo_set(kp%mos(1, 2), nmo=nb, nelectron=ne_b) @@ -937,6 +947,11 @@ CONTAINS ALLOCATE (weig(nmo, nkp, nspin), wocc(nmo, nkp, nspin)) weig = 0.0_dp wocc = 0.0_dp + IF (PRESENT(probe)) THEN + ALLOCATE (rcoeff(nao, nmo, nkp, nspin), icoeff(nao, nmo, nkp, nspin)) + rcoeff = 0.0_dp !coeff, real part + icoeff = 0.0_dp !coeff, imaginary part + END IF CALL get_kpoint_info(kpoint, kp_range=kp_range) kplocal = kp_range(2) - kp_range(1) + 1 DO ikpgr = 1, kplocal @@ -946,63 +961,159 @@ CONTAINS mo_set => kp%mos(1, ispin) CALL get_mo_set(mo_set, eigenvalues=eigenvalues) weig(1:nmo, ik, ispin) = eigenvalues(1:nmo) + IF (PRESENT(probe)) THEN + CALL get_mo_set(mo_set, mo_coeff=mo_coeff) + CALL cp_fm_get_info(mo_coeff, & + nrow_global=nrow_global, & + ncol_global=ncol_global) + ALLOCATE (smatrix(nrow_global, ncol_global)) + CALL cp_fm_get_submatrix(mo_coeff, smatrix) + + rcoeff(1:nao, 1:nmo, ik, ispin) = smatrix(1:nrow_global, 1:ncol_global) + + DEALLOCATE (smatrix) + + mo_set => kp%mos(2, ispin) + + CALL get_mo_set(mo_set, mo_coeff=mo_coeff) + CALL cp_fm_get_info(mo_coeff, & + nrow_global=nrow_global, & + ncol_global=ncol_global) + ALLOCATE (smatrix(nrow_global, ncol_global)) + CALL cp_fm_get_submatrix(mo_coeff, smatrix) + + icoeff(1:nao, 1:nmo, ik, ispin) = smatrix(1:nrow_global, 1:ncol_global) + + mo_set => kp%mos(1, ispin) + + DEALLOCATE (smatrix) + END IF END DO END DO CALL get_kpoint_info(kpoint, para_env_inter_kp=para_env_inter_kp) CALL para_env_inter_kp%sum(weig) + IF (PRESENT(probe)) THEN + CALL para_env_inter_kp%sum(rcoeff) + CALL para_env_inter_kp%sum(icoeff) + END IF + CALL get_kpoint_info(kpoint, wkp=wkp) - IF (smear%do_smear) THEN - ! finite electronic temperature - SELECT CASE (smear%method) - CASE (smear_fermi_dirac) - IF (nspin == 1) THEN - nel = REAL(nelectron, KIND=dp) - CALL Fermikp(wocc(:, :, 1), mus(1), kTS, weig(:, :, 1), nel, wkp, & - smear%electronic_temperature, 2.0_dp) - ELSE IF (smear%fixed_mag_mom > 0.0_dp) THEN - nel = REAL(ne_a, KIND=dp) - CALL Fermikp(wocc(:, :, 1), mus(1), kTS, weig(:, :, 1), nel, wkp, & - smear%electronic_temperature, 1.0_dp) - nel = REAL(ne_b, KIND=dp) - CALL Fermikp(wocc(:, :, 2), mus(2), kTS, weig(:, :, 2), nel, wkp, & - smear%electronic_temperature, 1.0_dp) - ELSE - nel = REAL(ne_a, KIND=dp) + REAL(ne_b, KIND=dp) - CALL Fermikp2(wocc(:, :, :), mu, kTS, weig(:, :, :), nel, wkp, & - smear%electronic_temperature) - kTS = kTS/2._dp - mus(1:2) = mu - END IF - CASE DEFAULT - CPABORT("kpoints: Selected smearing not (yet) supported") - END SELECT - ELSE - ! fixed occupations (2/1) + +!calling of HP module HERE, before smear + IF (PRESENT(probe)) THEN + smear%do_smear = .FALSE. !ensures smearing is switched off + IF (nspin == 1) THEN nel = REAL(nelectron, KIND=dp) - CALL Fermikp(wocc(:, :, 1), mus(1), kTS, weig(:, :, 1), nel, wkp, 0.0_dp, 2.0_dp) + CALL probe_occupancy_kp(wocc(:, :, :), mus(1), kTS, weig(:, :, :), rcoeff(:, :, :, :), icoeff(:, :, :, :), 2.0d0, & + probe, nel, wkp) ELSE - nel = REAL(ne_a, KIND=dp) - CALL Fermikp(wocc(:, :, 1), mus(1), kTS, weig(:, :, 1), nel, wkp, 0.0_dp, 1.0_dp) - nel = REAL(ne_b, KIND=dp) - CALL Fermikp(wocc(:, :, 2), mus(2), kTS, weig(:, :, 2), nel, wkp, 0.0_dp, 1.0_dp) + nel = REAL(ne_a, KIND=dp) + REAL(ne_b, KIND=dp) + CALL probe_occupancy_kp(wocc(:, :, :), mu, kTS, weig(:, :, :), rcoeff(:, :, :, :), icoeff(:, :, :, :), 1.0d0, & + probe, nel, wkp) + kTS = kTS/2._dp + mus(1:2) = mu END IF - END IF - DO ikpgr = 1, kplocal - ik = kp_range(1) + ikpgr - 1 - kp => kpoint%kp_env(ikpgr)%kpoint_env - DO ispin = 1, nspin - mo_set => kp%mos(1, ispin) - CALL get_mo_set(mo_set, eigenvalues=eigenvalues, occupation_numbers=occupation) - eigenvalues(1:nmo) = weig(1:nmo, ik, ispin) - occupation(1:nmo) = wocc(1:nmo, ik, ispin) - mo_set%kTS = kTS - mo_set%mu = mus(ispin) - END DO - END DO - DEALLOCATE (weig, wocc) + DO ikpgr = 1, kplocal + ik = kp_range(1) + ikpgr - 1 + kp => kpoint%kp_env(ikpgr)%kpoint_env + DO ispin = 1, nspin + mo_set => kp%mos(1, ispin) + CALL get_mo_set(mo_set, eigenvalues=eigenvalues, occupation_numbers=occupation) + eigenvalues(1:nmo) = weig(1:nmo, ik, ispin) + occupation(1:nmo) = wocc(1:nmo, ik, ispin) + mo_set%kTS = kTS + mo_set%mu = mus(ispin) + + CALL get_mo_set(mo_set, mo_coeff=mo_coeff) + !get smatrix for kpoint_env ikp + CALL cp_fm_get_info(mo_coeff, & + nrow_global=nrow_global, & + ncol_global=ncol_global) + ALLOCATE (smatrix(nrow_global, ncol_global)) + CALL cp_fm_get_submatrix(mo_coeff, smatrix) + + smatrix(1:nrow_global, 1:ncol_global) = rcoeff(1:nao, 1:nmo, ik, ispin) + DEALLOCATE (smatrix) + + mo_set => kp%mos(2, ispin) + + CALL get_mo_set(mo_set, mo_coeff=mo_coeff) + !get smatrix for kpoint_env ikp + CALL cp_fm_get_info(mo_coeff, & + nrow_global=nrow_global, & + ncol_global=ncol_global) + ALLOCATE (smatrix(nrow_global, ncol_global)) + CALL cp_fm_get_submatrix(mo_coeff, smatrix) + + smatrix(1:nrow_global, 1:ncol_global) = icoeff(1:nao, 1:nmo, ik, ispin) + DEALLOCATE (smatrix) + + mo_set => kp%mos(1, ispin) + + END DO + END DO + + DEALLOCATE (weig, wocc, rcoeff, icoeff) + + END IF + + IF (PRESENT(probe) .EQV. .FALSE.) THEN + IF (smear%do_smear) THEN + ! finite electronic temperature + SELECT CASE (smear%method) + CASE (smear_fermi_dirac) + IF (nspin == 1) THEN + nel = REAL(nelectron, KIND=dp) + CALL Fermikp(wocc(:, :, 1), mus(1), kTS, weig(:, :, 1), nel, wkp, & + smear%electronic_temperature, 2.0_dp) + ELSE IF (smear%fixed_mag_mom > 0.0_dp) THEN + nel = REAL(ne_a, KIND=dp) + CALL Fermikp(wocc(:, :, 1), mus(1), kTS, weig(:, :, 1), nel, wkp, & + smear%electronic_temperature, 1.0_dp) + nel = REAL(ne_b, KIND=dp) + CALL Fermikp(wocc(:, :, 2), mus(2), kTS, weig(:, :, 2), nel, wkp, & + smear%electronic_temperature, 1.0_dp) + ELSE + nel = REAL(ne_a, KIND=dp) + REAL(ne_b, KIND=dp) + CALL Fermikp2(wocc(:, :, :), mu, kTS, weig(:, :, :), nel, wkp, & + smear%electronic_temperature) + kTS = kTS/2._dp + mus(1:2) = mu + END IF + CASE DEFAULT + CPABORT("kpoints: Selected smearing not (yet) supported") + END SELECT + ELSE + ! fixed occupations (2/1) + IF (nspin == 1) THEN + nel = REAL(nelectron, KIND=dp) + CALL Fermikp(wocc(:, :, 1), mus(1), kTS, weig(:, :, 1), nel, wkp, 0.0_dp, 2.0_dp) + ELSE + nel = REAL(ne_a, KIND=dp) + CALL Fermikp(wocc(:, :, 1), mus(1), kTS, weig(:, :, 1), nel, wkp, 0.0_dp, 1.0_dp) + nel = REAL(ne_b, KIND=dp) + CALL Fermikp(wocc(:, :, 2), mus(2), kTS, weig(:, :, 2), nel, wkp, 0.0_dp, 1.0_dp) + END IF + END IF + DO ikpgr = 1, kplocal + ik = kp_range(1) + ikpgr - 1 + kp => kpoint%kp_env(ikpgr)%kpoint_env + DO ispin = 1, nspin + mo_set => kp%mos(1, ispin) + CALL get_mo_set(mo_set, eigenvalues=eigenvalues, occupation_numbers=occupation) + eigenvalues(1:nmo) = weig(1:nmo, ik, ispin) + occupation(1:nmo) = wocc(1:nmo, ik, ispin) + mo_set%kTS = kTS + mo_set%mu = mus(ispin) + END DO + END DO + + DEALLOCATE (weig, wocc) + + END IF CALL timestop(handle) diff --git a/src/qs_mo_methods.F b/src/qs_mo_methods.F index abb5d98a0b..284ac6138d 100644 --- a/src/qs_mo_methods.F +++ b/src/qs_mo_methods.F @@ -16,6 +16,7 @@ MODULE qs_mo_methods USE admm_utils, ONLY: admm_correct_for_eigenvalues,& admm_uncorrect_for_eigenvalues USE cp_blacs_env, ONLY: cp_blacs_env_type + USE cp_control_types, ONLY: hairy_probes_type USE cp_dbcsr_api, ONLY: dbcsr_copy,& dbcsr_get_info,& dbcsr_init_p,& @@ -809,11 +810,14 @@ CONTAINS !> \param scf_control ... !> \param mo_derivs ... !> \param admm_env ... +!> \param hairy_probes ... +!> \param probe ... !> \par History !> 02.2013 moved from qs_scf_post_gpw !> ! ************************************************************************************************** - SUBROUTINE make_mo_eig(mos, nspins, ks_rmpv, scf_control, mo_derivs, admm_env) + SUBROUTINE make_mo_eig(mos, nspins, ks_rmpv, scf_control, mo_derivs, admm_env, & + hairy_probes, probe) TYPE(mo_set_type), DIMENSION(:), INTENT(INOUT) :: mos INTEGER, INTENT(IN) :: nspins @@ -821,6 +825,9 @@ CONTAINS TYPE(scf_control_type), POINTER :: scf_control TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mo_derivs TYPE(admm_type), OPTIONAL, POINTER :: admm_env + LOGICAL, OPTIONAL :: hairy_probes + TYPE(hairy_probes_type), DIMENSION(:), OPTIONAL, & + POINTER :: probe CHARACTER(len=*), PARAMETER :: routineN = 'make_mo_eig' @@ -868,8 +875,16 @@ CONTAINS ELSE IF (output_unit > 0) WRITE (output_unit, '(4(1X,1F16.8))') mo_eigenvalues(1:homo) END IF - IF (.NOT. scf_control%diagonalization%mom) & - CALL set_mo_occupation(mo_set=mos(ispin), smear=scf_control%smear) + IF (.NOT. scf_control%diagonalization%mom) THEN + IF (PRESENT(hairy_probes) .EQV. .TRUE.) THEN + scf_control%smear%do_smear = .FALSE. + CALL set_mo_occupation(mo_set=mos(ispin), & + smear=scf_control%smear, & + probe=probe) + ELSE + CALL set_mo_occupation(mo_set=mos(ispin), smear=scf_control%smear) + END IF + END IF IF (output_unit > 0) WRITE (output_unit, '(T2,A,F12.6)') & "Fermi Energy [eV] :", mos(ispin)%mu*evolt END DO diff --git a/src/qs_mo_occupation.F b/src/qs_mo_occupation.F index 123821107a..3f14d5228d 100644 --- a/src/qs_mo_occupation.F +++ b/src/qs_mo_occupation.F @@ -14,9 +14,11 @@ MODULE qs_mo_occupation + USE cp_control_types, ONLY: hairy_probes_type USE cp_log_handling, ONLY: cp_to_string USE fermi_utils, ONLY: FermiFixed,& FermiFixedDeriv + USE hairy_probes, ONLY: probe_occupancy USE input_constants, ONLY: smear_energy_window,& smear_fermi_dirac,& smear_list @@ -181,6 +183,7 @@ CONTAINS !> \param smear ... !> \param eval_deriv ... !> \param tot_zeff_corr ... +!> \param probe ... !> \date 25.01.2010 (MK) !> \par History !> 10.2019 Added functionality to adjust mo occupation if the core @@ -191,12 +194,14 @@ CONTAINS !> \author Matthias Krack (MK) !> \version 1.0 ! ************************************************************************************************** - SUBROUTINE set_mo_occupation_2(mo_array, smear, eval_deriv, tot_zeff_corr) + SUBROUTINE set_mo_occupation_2(mo_array, smear, eval_deriv, tot_zeff_corr, probe) TYPE(mo_set_type), DIMENSION(:), INTENT(INOUT) :: mo_array TYPE(smear_type), POINTER :: smear REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER :: eval_deriv REAL(KIND=dp), OPTIONAL :: tot_zeff_corr + TYPE(hairy_probes_type), DIMENSION(:), OPTIONAL, & + POINTER :: probe CHARACTER(LEN=*), PARAMETER :: routineN = 'set_mo_occupation_2' @@ -210,7 +215,9 @@ CONTAINS ! Fall back for the case that we have only one MO set IF (SIZE(mo_array) == 1) THEN - IF (PRESENT(eval_deriv)) THEN + IF (PRESENT(probe)) THEN + CALL set_mo_occupation_1(mo_array(1), smear=smear, probe=probe) + ELSE IF (PRESENT(eval_deriv)) THEN ! Change of MO occupancy to account for CORE_CORRECTION is not yet implemented CALL set_mo_occupation_1(mo_array(1), smear=smear, eval_deriv=eval_deriv) ELSE @@ -224,6 +231,11 @@ CONTAINS RETURN END IF + IF (PRESENT(probe)) THEN + CALL set_mo_occupation_1(mo_array(1), smear=smear, probe=probe) + CALL set_mo_occupation_1(mo_array(2), smear=smear, probe=probe) + END IF + IF (smear%do_smear) THEN IF (smear%fixed_mag_mom < 0.0_dp) THEN IF (PRESENT(tot_zeff_corr)) THEN @@ -259,7 +271,10 @@ CONTAINS IF (.NOT. ((mo_array(1)%flexible_electron_count > 0.0_dp) .AND. & (mo_array(2)%flexible_electron_count > 0.0_dp))) THEN - IF (PRESENT(eval_deriv)) THEN + IF (PRESENT(probe)) THEN + CALL set_mo_occupation_1(mo_array(1), smear=smear, probe=probe) + CALL set_mo_occupation_1(mo_array(2), smear=smear, probe=probe) + ELSE IF (PRESENT(eval_deriv)) THEN CALL set_mo_occupation_1(mo_array(1), smear=smear, eval_deriv=eval_deriv) CALL set_mo_occupation_1(mo_array(2), smear=smear, eval_deriv=eval_deriv) ELSE @@ -349,7 +364,10 @@ CONTAINS TRIM(ADJUSTL(cp_to_string(multiplicity_old)))//" to "// & TRIM(ADJUSTL(cp_to_string(multiplicity_new)))) - IF (PRESENT(eval_deriv)) THEN + IF (PRESENT(probe)) THEN + CALL set_mo_occupation_1(mo_array(1), smear=smear, probe=probe) + CALL set_mo_occupation_1(mo_array(2), smear=smear, probe=probe) + ELSE IF (PRESENT(eval_deriv)) THEN CALL set_mo_occupation_1(mo_array(1), smear=smear, eval_deriv=eval_deriv) CALL set_mo_occupation_1(mo_array(2), smear=smear, eval_deriv=eval_deriv) ELSE @@ -374,6 +392,7 @@ CONTAINS !> on exit the derivative of the full free energy (i.e. KS and entropy) wrt to the eigenvalue !> \param xas_env ... !> \param tot_zeff_corr ... +!> \param probe ... !> \date 17.04.2002 (v1.0), 26.08.2008 (v1.1) !> \par History !> 10.2019 Added functionality to adjust mo occupation if the core @@ -384,13 +403,15 @@ CONTAINS !> \author Matthias Krack !> \version 1.1 ! ************************************************************************************************** - SUBROUTINE set_mo_occupation_1(mo_set, smear, eval_deriv, xas_env, tot_zeff_corr) + SUBROUTINE set_mo_occupation_1(mo_set, smear, eval_deriv, xas_env, tot_zeff_corr, probe) TYPE(mo_set_type), INTENT(INOUT) :: mo_set TYPE(smear_type), OPTIONAL, POINTER :: smear REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER :: eval_deriv TYPE(xas_environment_type), OPTIONAL, POINTER :: xas_env REAL(KIND=dp), OPTIONAL :: tot_zeff_corr + TYPE(hairy_probes_type), DIMENSION(:), OPTIONAL, & + POINTER :: probe CHARACTER(LEN=*), PARAMETER :: routineN = 'set_mo_occupation_1' @@ -510,6 +531,54 @@ CONTAINS CPASSERT(equal_size) END IF +!calling of HP module HERE, before smear + IF (PRESENT(probe)) THEN + i_first = 1 + IF (smear%fixed_mag_mom == -1.0_dp) THEN + nelec = REAL(mo_set%nelectron, dp) + ELSE + nelec = mo_set%n_el_f + END IF + + mo_set%occupation_numbers(:) = 0.0_dp + + CALL probe_occupancy(mo_set%occupation_numbers, mo_set%mu, mo_set%kTS, & + mo_set%eigenvalues, mo_set%mo_coeff, mo_set%maxocc, & + probe, N=nelec) + !NB: mu and T are taken from the hairy_probe type (defined in cp_control_types.F); these values are set in the input + + ! Find the lowest fractional occupied MO (LFOMO) + DO imo = i_first, nmo + IF (mo_set%occupation_numbers(imo) < mo_set%maxocc) THEN + mo_set%lfomo = imo + EXIT + END IF + END DO + is_large = ABS(MAXVAL(mo_set%occupation_numbers) - mo_set%maxocc) > probe(1)%eps_hp + ! this is not a real problem, but the temperature might be a bit large + IF (is_large) & + CPWARN("Hair-probes occupancy distribution includes the first MO") + + ! Find the highest (fractional) occupied MO which will be now the HOMO + DO imo = nmo, mo_set%lfomo, -1 + IF (mo_set%occupation_numbers(imo) > probe(1)%eps_hp) THEN + mo_set%homo = imo + EXIT + END IF + END DO + is_large = ABS(MINVAL(mo_set%occupation_numbers)) > probe(1)%eps_hp + IF (is_large) & + CALL cp_warn(__LOCATION__, & + "Hair-probes occupancy distribution includes the last MO => "// & + "Add more MOs for proper smearing.") + + ! check that the total electron count is accurate + is_large = (ABS(nelec - accurate_sum(mo_set%occupation_numbers(:))) > probe(1)%eps_hp*nelec) + IF (is_large) & + CPWARN("Total number of electrons is not accurate") + + END IF + ! Quick return, if no smearing information is supplied (TO BE FIXED, smear should become non-optional...) IF (.NOT. PRESENT(smear)) THEN ! there is no dependence of the energy on the eigenvalues diff --git a/src/qs_scf.F b/src/qs_scf.F index 4f4d361f38..7d2fa565a0 100644 --- a/src/qs_scf.F +++ b/src/qs_scf.F @@ -501,10 +501,21 @@ CONTAINS IF (do_kpoints) THEN ! kpoints - CALL qs_scf_new_mos_kp(qs_env, scf_env, scf_control, diis_step) + IF (dft_control%hairy_probes .EQV. .TRUE.) THEN + scf_control%smear%do_smear = .FALSE. + CALL qs_scf_new_mos_kp(qs_env, scf_env, scf_control, diis_step, dft_control%probe) + ELSE + CALL qs_scf_new_mos_kp(qs_env, scf_env, scf_control, diis_step) + END IF ELSE ! Gamma points only - CALL qs_scf_new_mos(qs_env, scf_env, scf_control, scf_section, diis_step, energy_only) + IF (dft_control%hairy_probes .EQV. .TRUE.) THEN + scf_control%smear%do_smear = .FALSE. + CALL qs_scf_new_mos(qs_env, scf_env, scf_control, scf_section, diis_step, energy_only, & + dft_control%probe) + ELSE + CALL qs_scf_new_mos(qs_env, scf_env, scf_control, scf_section, diis_step, energy_only) + END IF END IF ! Print requested MO information (can be computationally expensive with OT) @@ -699,9 +710,21 @@ CONTAINS ! this just guarantees that all mo_occupations match the eigenvalues, if smear DO ispin = 1, dft_control%nspins ! do not reset mo_occupations if the maximum overlap method is in use - IF (.NOT. scf_control%diagonalization%mom) & - CALL set_mo_occupation(mo_set=mos(ispin), & - smear=scf_control%smear) + IF (.NOT. scf_control%diagonalization%mom) THEN + !if the hair probes section is present, this sends hairy_probes to set_mo_occupation subroutine + !and switches off the standard smearing + IF (dft_control%hairy_probes .EQV. .TRUE.) THEN + IF (scf_env%outer_scf%iter_count > 0) THEN + scf_control%smear%do_smear = .FALSE. + CALL set_mo_occupation(mo_set=mos(ispin), & + smear=scf_control%smear, & + probe=dft_control%probe) + END IF + ELSE + CALL set_mo_occupation(mo_set=mos(ispin), & + smear=scf_control%smear) + END IF + END IF END DO SELECT CASE (scf_env%method) diff --git a/src/qs_scf_diagonalization.F b/src/qs_scf_diagonalization.F index 12d376fe6b..a4c605e7c8 100644 --- a/src/qs_scf_diagonalization.F +++ b/src/qs_scf_diagonalization.F @@ -24,7 +24,8 @@ MODULE qs_scf_diagonalization cp_cfm_to_cfm,& cp_cfm_to_fm,& cp_cfm_type - USE cp_control_types, ONLY: dft_control_type + USE cp_control_types, ONLY: dft_control_type,& + hairy_probes_type USE cp_dbcsr_api, ONLY: & dbcsr_copy, dbcsr_create, dbcsr_deallocate_matrix, dbcsr_desymmetrize, dbcsr_p_type, & dbcsr_set, dbcsr_type, dbcsr_type_antisymmetric, dbcsr_type_no_symmetry, & @@ -346,10 +347,11 @@ CONTAINS !> \param scf_control ... !> \param scf_section ... !> \param diis_step ... +!> \param probe ... ! ************************************************************************************************** SUBROUTINE do_general_diag(scf_env, mos, matrix_ks, & matrix_s, scf_control, scf_section, & - diis_step) + diis_step, probe) TYPE(qs_scf_env_type), POINTER :: scf_env TYPE(mo_set_type), DIMENSION(:), INTENT(INOUT) :: mos @@ -357,6 +359,8 @@ CONTAINS TYPE(scf_control_type), POINTER :: scf_control TYPE(section_vals_type), POINTER :: scf_section LOGICAL, INTENT(INOUT) :: diis_step + TYPE(hairy_probes_type), DIMENSION(:), OPTIONAL, & + POINTER :: probe INTEGER :: ispin, nspin REAL(KIND=dp) :: total_zeff_corr @@ -373,8 +377,15 @@ CONTAINS CALL set_mo_occupation(mo_array=mos, & smear=scf_control%smear, tot_zeff_corr=total_zeff_corr) ELSE - CALL set_mo_occupation(mo_array=mos, & - smear=scf_control%smear) + IF (PRESENT(probe) .EQV. .TRUE.) THEN + scf_control%smear%do_smear = .FALSE. + CALL set_mo_occupation(mo_array=mos, & + smear=scf_control%smear, & + probe=probe) + ELSE + CALL set_mo_occupation(mo_array=mos, & + smear=scf_control%smear) + END IF END IF DO ispin = 1, nspin @@ -399,11 +410,12 @@ CONTAINS !> \param diis_step ... !> \param diis_error ... !> \param qs_env ... +!> \param probe ... !> \par History !> 08.2014 created [JGH] ! ************************************************************************************************** SUBROUTINE do_general_diag_kp(matrix_ks, matrix_s, kpoints, scf_env, scf_control, update_p, & - diis_step, diis_error, qs_env) + diis_step, diis_error, qs_env, probe) TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks, matrix_s TYPE(kpoint_type), POINTER :: kpoints @@ -413,6 +425,8 @@ CONTAINS LOGICAL, INTENT(INOUT) :: diis_step REAL(dp), INTENT(INOUT), OPTIONAL :: diis_error TYPE(qs_environment_type), OPTIONAL, POINTER :: qs_env + TYPE(hairy_probes_type), DIMENSION(:), OPTIONAL, & + POINTER :: probe CHARACTER(len=*), PARAMETER :: routineN = 'do_general_diag_kp' COMPLEX(KIND=dp), PARAMETER :: cone = CMPLX(1.0_dp, 0.0_dp, KIND=dp), & @@ -711,8 +725,13 @@ CONTAINS IF (update_p) THEN ! MO occupations - CALL kpoint_set_mo_occupation(kpoints, scf_control%smear) - + IF (PRESENT(probe) .EQV. .TRUE.) THEN + scf_control%smear%do_smear = .FALSE. + CALL kpoint_set_mo_occupation(kpoints, scf_control%smear, & + probe=probe) + ELSE + CALL kpoint_set_mo_occupation(kpoints, scf_control%smear) + END IF ! density matrices CALL kpoint_density_matrices(kpoints) ! density matrices in real space diff --git a/src/qs_scf_initialization.F b/src/qs_scf_initialization.F index b54bb5d3e9..fe51cf35f9 100644 --- a/src/qs_scf_initialization.F +++ b/src/qs_scf_initialization.F @@ -9,6 +9,7 @@ !> \brief Utility routines for qs_scf ! ************************************************************************************************** MODULE qs_scf_initialization + USE atomic_kind_types, ONLY: atomic_kind_type USE cp_control_types, ONLY: dft_control_type USE cp_dbcsr_api, ONLY: dbcsr_create,& dbcsr_init_p,& @@ -47,6 +48,7 @@ MODULE qs_scf_initialization cp_print_key_finished_output,& cp_print_key_should_output,& cp_print_key_unit_nr + USE hairy_probes, ONLY: AO_boundaries USE input_constants, ONLY: & broy_mix, cholesky_dbcsr, cholesky_inverse, cholesky_off, diag_block_davidson, & diag_block_krylov, diag_filter_matrix, diag_ot, diag_standard, direct_p_mix, kerker_mix, & @@ -60,6 +62,7 @@ MODULE qs_scf_initialization USE kpoint_types, ONLY: kpoint_type USE message_passing, ONLY: mp_para_env_type USE parallel_gemm_api, ONLY: parallel_gemm + USE particle_types, ONLY: particle_type USE pw_types, ONLY: pw_c1d_gs_type USE qmmm_image_charge, ONLY: conditional_calc_image_matrix USE qs_block_davidson_types, ONLY: block_davidson_allocate,& @@ -145,12 +148,28 @@ CONTAINS TYPE(scf_control_type), OPTIONAL, POINTER :: scf_control TYPE(section_vals_type), OPTIONAL, POINTER :: scf_section + INTEGER :: ip, np + TYPE(atomic_kind_type), POINTER :: atomic_kind_set(:) TYPE(dft_control_type), POINTER :: dft_control + TYPE(mo_set_type), DIMENSION(:), POINTER :: mos + TYPE(particle_type), POINTER :: particle_set(:) + TYPE(qs_kind_type), POINTER :: qs_kind_set(:) TYPE(scf_control_type), POINTER :: my_scf_control TYPE(section_vals_type), POINTER :: dft_section, input, my_scf_section CALL get_qs_env(qs_env, input=input, dft_control=dft_control) + !Initialize Hairy Probe calculation + IF (dft_control%hairy_probes .EQV. .TRUE.) THEN + CALL get_qs_env(qs_env, mos=mos, & + atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set, particle_set=particle_set) + np = SIZE(dft_control%probe) + DO ip = 1, np + CALL AO_boundaries(probe=dft_control%probe(ip), atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set, & + particle_set=particle_set, nAO=mos(1)%nao) !FIX THIS! + END DO + END IF + IF (PRESENT(scf_control)) THEN my_scf_control => scf_control ELSE @@ -1202,8 +1221,14 @@ CONTAINS DO ispin = 1, SIZE(mos) CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff, nmo=nmo) CALL reorthogonalize_vectors(qs_env, v_matrix=mo_coeff, n_col=nmo) - CALL set_mo_occupation(mo_set=mos(ispin), & - smear=scf_control%smear) + IF (dft_control%hairy_probes .EQV. .TRUE.) THEN + scf_control%smear%do_smear = .FALSE. + CALL set_mo_occupation(mo_set=mos(ispin), & + smear=scf_control%smear, probe=dft_control%probe) + ELSE + CALL set_mo_occupation(mo_set=mos(ispin), & + smear=scf_control%smear) + END IF END DO END IF END IF diff --git a/src/qs_scf_loop_utils.F b/src/qs_scf_loop_utils.F index 7aa1310ae1..38c52b89bc 100644 --- a/src/qs_scf_loop_utils.F +++ b/src/qs_scf_loop_utils.F @@ -8,7 +8,8 @@ !> \brief Utility routines for qs_scf ! ************************************************************************************************** MODULE qs_scf_loop_utils - USE cp_control_types, ONLY: dft_control_type + USE cp_control_types, ONLY: dft_control_type,& + hairy_probes_type USE cp_dbcsr_api, ONLY: dbcsr_copy,& dbcsr_get_info,& dbcsr_p_type,& @@ -114,14 +115,17 @@ CONTAINS !> \param scf_section ... !> \param diis_step ... !> \param energy_only ... +!> \param probe ... ! ************************************************************************************************** SUBROUTINE qs_scf_new_mos(qs_env, scf_env, scf_control, scf_section, diis_step, & - energy_only) + energy_only, probe) TYPE(qs_environment_type), POINTER :: qs_env TYPE(qs_scf_env_type), POINTER :: scf_env TYPE(scf_control_type), POINTER :: scf_control TYPE(section_vals_type), POINTER :: scf_section LOGICAL :: diis_step, energy_only + TYPE(hairy_probes_type), DIMENSION(:), OPTIONAL, & + POINTER :: probe CHARACTER(LEN=*), PARAMETER :: routineN = 'qs_scf_new_mos' @@ -204,9 +208,16 @@ CONTAINS matrix_s, scf_control, scf_section, & diis_step) ELSE - CALL do_general_diag(scf_env, mos, matrix_ks, & - matrix_s, scf_control, scf_section, & - diis_step) + IF (dft_control%hairy_probes .EQV. .TRUE.) THEN + CALL do_general_diag(scf_env, mos, matrix_ks, & + matrix_s, scf_control, scf_section, & + diis_step, & + probe) + ELSE + CALL do_general_diag(scf_env, mos, matrix_ks, & + matrix_s, scf_control, scf_section, & + diis_step) + END IF END IF IF (scf_control%do_diag_sub) THEN skip_diag_sub = (scf_env%subspace_env%eps_diag_sub > 0.0_dp) .AND. & @@ -285,12 +296,15 @@ CONTAINS !> \param scf_env ... !> \param scf_control ... !> \param diis_step ... +!> \param probe ... ! ************************************************************************************************** - SUBROUTINE qs_scf_new_mos_kp(qs_env, scf_env, scf_control, diis_step) + SUBROUTINE qs_scf_new_mos_kp(qs_env, scf_env, scf_control, diis_step, probe) TYPE(qs_environment_type), POINTER :: qs_env TYPE(qs_scf_env_type), POINTER :: scf_env TYPE(scf_control_type), POINTER :: scf_control LOGICAL :: diis_step + TYPE(hairy_probes_type), DIMENSION(:), OPTIONAL, & + POINTER :: probe CHARACTER(LEN=*), PARAMETER :: routineN = 'qs_scf_new_mos_kp' @@ -321,8 +335,14 @@ CONTAINS CASE (general_diag_method_nr) ! Diagonlization in non orthonormal case CALL get_qs_env(qs_env, matrix_ks_kp=matrix_ks, matrix_s_kp=matrix_s) - CALL do_general_diag_kp(matrix_ks, matrix_s, kpoints, scf_env, scf_control, .TRUE., & - diis_step, diis_error, qs_env) + IF (dft_control%hairy_probes .EQV. .TRUE.) THEN + scf_control%smear%do_smear = .FALSE. + CALL do_general_diag_kp(matrix_ks, matrix_s, kpoints, scf_env, scf_control, .TRUE., & + diis_step, diis_error, qs_env, probe) + ELSE + CALL do_general_diag_kp(matrix_ks, matrix_s, kpoints, scf_env, scf_control, .TRUE., & + diis_step, diis_error, qs_env) + END IF IF (diis_step) THEN scf_env%iter_param = diis_error scf_env%iter_method = "DIIS/Diag." diff --git a/src/qs_scf_output.F b/src/qs_scf_output.F index 1607194d16..1154fcd734 100644 --- a/src/qs_scf_output.F +++ b/src/qs_scf_output.F @@ -760,6 +760,12 @@ CONTAINS "GAPW_XC| Exc from hard and soft atomic rho1: ", exc1_energy END IF END IF + IF (dft_control%hairy_probes .EQV. .TRUE.) THEN + WRITE (UNIT=output_unit, FMT="((T3,A,T56,F25.14))") & + "Electronic entropic energy:", energy%kTS + WRITE (UNIT=output_unit, FMT="((T3,A,T56,F25.14))") & + "Fermi energy:", energy%efermi + END IF IF (dft_control%smear) THEN WRITE (UNIT=output_unit, FMT="((T3,A,T56,F25.14))") & "Electronic entropic energy:", energy%kTS diff --git a/src/qs_scf_post_gpw.F b/src/qs_scf_post_gpw.F index 93be3d73fc..ba55a6a19b 100644 --- a/src/qs_scf_post_gpw.F +++ b/src/qs_scf_post_gpw.F @@ -478,7 +478,14 @@ CONTAINS CALL get_qs_env(qs_env, admm_env=admm_env) CALL make_mo_eig(mos, nspins, ks_rmpv, scf_control, mo_derivs, admm_env=admm_env) ELSE - CALL make_mo_eig(mos, nspins, ks_rmpv, scf_control, mo_derivs) + IF (dft_control%hairy_probes) THEN + scf_control%smear%do_smear = .FALSE. + CALL make_mo_eig(mos, dft_control%nspins, ks_rmpv, scf_control, mo_derivs, & + hairy_probes=dft_control%hairy_probes, & + probe=dft_control%probe) + ELSE + CALL make_mo_eig(mos, dft_control%nspins, ks_rmpv, scf_control, mo_derivs) + END IF END IF DO ispin = 1, dft_control%nspins CALL get_mo_set(mo_set=mos(ispin), eigenvalues=mo_eigenvalues, homo=homo) diff --git a/tests/HP-DFT/TEST_FILES.toml b/tests/HP-DFT/TEST_FILES.toml new file mode 100644 index 0000000000..6166ae3da2 --- /dev/null +++ b/tests/HP-DFT/TEST_FILES.toml @@ -0,0 +1,10 @@ +# runs are executed in the same order as in this file +# the second field tells which test should be run in order to compare with the last available output +# e.g. 0 means do not compare anything, running is enough +# 1 compares the last total energy in the file +# for details see cp2k/tests/do_regtest +"hp-dft_0eV.inp" = [{matcher="M011", tol=5.0E-2, ref=-236.439459678130504}, + {matcher="M126", tol=1.0E-6, ref=0.0}] +"hp-dft_1eV.inp" = [{matcher="M011", tol=5.0E-2, ref=-238.228758936954392}, + {matcher="M126", tol=1.0E-6, ref=0.0}] +#EOF \ No newline at end of file diff --git a/tests/HP-DFT/hp-dft_0eV.inp b/tests/HP-DFT/hp-dft_0eV.inp new file mode 100644 index 0000000000..3708f5ea49 --- /dev/null +++ b/tests/HP-DFT/hp-dft_0eV.inp @@ -0,0 +1,74 @@ +&GLOBAL + PRINT_LEVEL medium + PROJECT_NAME Pt_capacitor + RUN_TYPE ENERGY +&END GLOBAL + +&FORCE_EVAL + METHOD QS + &DFT + BASIS_SET_FILE_NAME BASIS_MOLOPT + POTENTIAL_FILE_NAME GTH_POTENTIALS + &HAIRY_PROBES + ATOM_IDS 1 + EPS_HP 1.0E-5 + MU 0.0 + T 300 + &END HAIRY_PROBES + &HAIRY_PROBES + ATOM_IDS 2 + EPS_HP 1.0E-5 + MU 0.0 + T 300 + &END HAIRY_PROBES + &MGRID + CUTOFF 250 + NGRIDS 4 + REL_CUTOFF 50 + &END MGRID + &QS + EPS_DEFAULT 1.0E-14 + &END QS + &SCF + ADDED_MOS -1 + CHOLESKY INVERSE + EPS_DIIS 1.0E-01 + EPS_SCF 1.0E-05 + MAX_SCF 500 + SCF_GUESS ATOMIC + &DIAGONALIZATION + ALGORITHM STANDARD + &END DIAGONALIZATION + &MIXING + ALPHA 0.1 + BETA 0.5 + METHOD BROYDEN_MIXING + NBROYDEN 8 + &END MIXING + &OUTER_SCF + EPS_SCF 1.0E-05 + MAX_SCF 60 + &END OUTER_SCF + &END SCF + &XC + &XC_FUNCTIONAL PBE + &END XC_FUNCTIONAL + &END XC + &END DFT + &SUBSYS + &CELL + ABC 2.807213921310594 2.807213921310594 20.0 + ALPHA_BETA_GAMMA 90 90 60 + PERIODIC XY + &END CELL + &COORD + Pt 0.00000000 0.00000000 5.00 + Pt 1.40360696 0.81037286 15.0 + &END COORD + &KIND Pt + BASIS_SET SZV-MOLOPT-SR-GTH + ELEMENT Pt + POTENTIAL GTH-PBE-q18 + &END KIND + &END SUBSYS +&END FORCE_EVAL diff --git a/tests/HP-DFT/hp-dft_1eV.inp b/tests/HP-DFT/hp-dft_1eV.inp new file mode 100644 index 0000000000..7902d721bb --- /dev/null +++ b/tests/HP-DFT/hp-dft_1eV.inp @@ -0,0 +1,74 @@ +&GLOBAL + PRINT_LEVEL medium + PROJECT_NAME Pt_capacitor + RUN_TYPE ENERGY +&END GLOBAL + +&FORCE_EVAL + METHOD QS + &DFT + BASIS_SET_FILE_NAME BASIS_MOLOPT + POTENTIAL_FILE_NAME GTH_POTENTIALS + &HAIRY_PROBES + ATOM_IDS 1 + EPS_HP 1.0E-5 + MU 0.5 + T 300 + &END HAIRY_PROBES + &HAIRY_PROBES + ATOM_IDS 2 + EPS_HP 1.0E-5 + MU -0.5 + T 300 + &END HAIRY_PROBES + &MGRID + CUTOFF 250 + NGRIDS 4 + REL_CUTOFF 50 + &END MGRID + &QS + EPS_DEFAULT 1.0E-14 + &END QS + &SCF + ADDED_MOS -1 + CHOLESKY INVERSE + EPS_DIIS 1.0E-01 + EPS_SCF 1.0E-05 + MAX_SCF 500 + SCF_GUESS ATOMIC + &DIAGONALIZATION + ALGORITHM STANDARD + &END DIAGONALIZATION + &MIXING + ALPHA 0.1 + BETA 0.5 + METHOD BROYDEN_MIXING + NBROYDEN 8 + &END MIXING + &OUTER_SCF + EPS_SCF 1.0E-05 + MAX_SCF 60 + &END OUTER_SCF + &END SCF + &XC + &XC_FUNCTIONAL PBE + &END XC_FUNCTIONAL + &END XC + &END DFT + &SUBSYS + &CELL + ABC 2.807213921310594 2.807213921310594 20.0 + ALPHA_BETA_GAMMA 90 90 60 + PERIODIC XY + &END CELL + &COORD + Pt 0.00000000 0.00000000 5.00 + Pt 1.40360696 0.81037286 15.0 + &END COORD + &KIND Pt + BASIS_SET SZV-MOLOPT-SR-GTH + ELEMENT Pt + POTENTIAL GTH-PBE-q18 + &END KIND + &END SUBSYS +&END FORCE_EVAL diff --git a/tests/TEST_DIRS b/tests/TEST_DIRS index 3c006eba74..ec193e99ee 100644 --- a/tests/TEST_DIRS +++ b/tests/TEST_DIRS @@ -3,6 +3,7 @@ # Directories have been reordered according the execution time needed for a gfortran pdbg run using 2 MPI tasks # in case a new directory is added just add it at the top of the list.. # the order will be regularly checked and modified... +HP-DFT xTB/regtest-tblite-ipea1 tblite xTB/regtest-tblite-gfn2 tblite xTB/regtest-tblite-gfn1 tblite @@ -388,4 +389,4 @@ Fist/regtest-quip quip QS/regtest-eht-guess libdftd4 QS/regtest-trexio trexio QS/regtest-trexio-2 trexio libgrpp -QS/regtest-rtbse-gxac libint greenx +QS/regtest-rtbse-gxac libint greenx \ No newline at end of file diff --git a/tests/matchers.py b/tests/matchers.py index 23ac459821..01eb8191fa 100644 --- a/tests/matchers.py +++ b/tests/matchers.py @@ -250,4 +250,5 @@ registry["RTBSE_GXAC_H2_pol"] = GenericMatcher( r"POLARIZABILITY_PADE| 0.30450000E+002", col=4 ) +registry["M126"] = GenericMatcher(r" # Total charge ", col=5) # EOF