comp-chem/cp2k
1
0
Fork 0
mirror of https://github.com/cp2k/cp2k.git synced 2026-07-21 06:25:15 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

HP-DFT modules and regtests (#4138)

Browse source
This commit is contained in:
MaggieB092 2025-05-15 13:35:20 +01:00 • committed by GitHub
parent 5ea6a88af0
commit 61d4209430
No known key found for this signature in database
GPG key ID: B5690EEEBB952194
19 changed files with 1342 additions and 90 deletions
Download patch file Download diff file Expand all files Collapse all files

1
src/CMakeLists.txt
View file

@ -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

30
src/cp_control_types.F
View file

@ -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)

86
src/cp_control_utils.F
View file

@ -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

621
src/hairy_probes.F Normal file
View file

@ -0,0 +1,621 @@
!--------------------------------------------------------------------------------------------------!
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
! !
! 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

71
src/input_cp2k_dft.F
View file

@ -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 <INTEGER> .. <INTEGER>", &
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 <REAL>", &
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 <REAL>", &
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 <REAL>", &
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 <REAL>", &
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 ...

215
src/kpoint_methods.F
View file

@ -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)

21
src/qs_mo_methods.F
View file

@ -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

79
src/qs_mo_occupation.F
View file

@ -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

33
src/qs_scf.F
View file

@ -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)

33
src/qs_scf_diagonalization.F
View file

@ -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

29
src/qs_scf_initialization.F
View file

@ -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

36
src/qs_scf_loop_utils.F
View file

@ -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."

6
src/qs_scf_output.F
View file

@ -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

9
src/qs_scf_post_gpw.F
View file

@ -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)

10
tests/HP-DFT/TEST_FILES.toml Normal file
View file

@ -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

74
tests/HP-DFT/hp-dft_0eV.inp Normal file
View file

@ -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

74
tests/HP-DFT/hp-dft_1eV.inp Normal file
View file

@ -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

3
tests/TEST_DIRS
View file

@ -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

1
tests/matchers.py
View file

@ -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