kpoint_smearing_unittest.F

Go to the documentation of this file.

!--------------------------------------------------------------------------------------------------!
!   CP2K: A general program to perform molecular dynamics simulations                              !
!   Copyright 2000-2026 CP2K developers group <https://cp2k.org>                                   !
!                                                                                                  !
!   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
!--------------------------------------------------------------------------------------------------!
 
PROGRAM kpoint_smearing_unittest
   USE bibliography,                    ONLY: add_all_references
   USE input_constants,                 ONLY: smear_fermi_dirac,&
                                              smear_gaussian,&
                                              smear_mp,&
                                              smear_mv
   USE kinds,                           ONLY: dp
   USE kpoint_methods,                  ONLY: kpoint_smearing_edge_status
   USE reference_manager,               ONLY: remove_all_references
   USE scf_control_types,               ONLY: smear_type
   USE smearing_utils,                  ONLY: smearkp,&
                                              smearing_response_weight,&
                                              smearing_separate_spin,&
                                              smearing_spin_targets
 
   IMPLICIT NONE
 
   LOGICAL                                            :: first_fractional, has_weight, last_occupied
   LOGICAL, DIMENSION(2)                              :: last_occupied_spin
   REAL(kind=dp), DIMENSION(2)                        :: weight
   REAL(kind=dp), DIMENSION(3, 2, 2)                  :: occupation
   TYPE(smear_type)                                   :: smear
 
   CALL add_all_references()
   smear%do_smear = .true.
   smear%eps_fermi_dirac = 1.0e-5_dp
   occupation = 0.0_dp
   occupation(1, :, :) = 1.0_dp
   weight = [1.0_dp, 0.0_dp]
 
   occupation(1, 2, 1) = 0.5_dp
   occupation(3, 2, 2) = 0.5_dp
   CALL kpoint_smearing_edge_status(occupation, weight, smear, 2, has_weight, &
                                    first_fractional, last_occupied, last_occupied_spin)
   IF (.NOT. has_weight) error stop "Positive K-point weight was not detected"
   IF (first_fractional) error stop "Zero-weight first-band occupation was included"
   IF (last_occupied) error stop "Zero-weight last-band occupation was included"
   IF (any(last_occupied_spin)) error stop "Zero-weight spin edge was included"
 
   weight(2) = 0.25_dp
   CALL kpoint_smearing_edge_status(occupation, weight, smear, 2, has_weight, &
                                    first_fractional, last_occupied, last_occupied_spin)
   IF (.NOT. first_fractional) error stop "Fractional first band was not detected"
   IF (.NOT. last_occupied) error stop "Occupied last band was not detected"
   IF (last_occupied_spin(1)) error stop "Wrong occupied spin edge was reported"
   IF (.NOT. last_occupied_spin(2)) error stop "Occupied spin edge was not reported"
 
   occupation = 0.0_dp
   occupation(1, :, 1) = 2.0_dp
   CALL kpoint_smearing_edge_status(occupation(:, :, 1:1), weight, smear, 1, has_weight, &
                                    first_fractional, last_occupied)
   IF (first_fractional) error stop "Restricted full occupation was classified as fractional"
 
   weight = 0.0_dp
   CALL kpoint_smearing_edge_status(occupation, weight, smear, 2, has_weight, &
                                    first_fractional, last_occupied, last_occupied_spin)
   IF (has_weight .OR. first_fractional .OR. last_occupied) THEN
      error stop "Zero-weight K-point set produced an edge status"
   END IF
 
   smear%do_smear = .false.
   weight = 1.0_dp
   CALL kpoint_smearing_edge_status(occupation, weight, smear, 2, has_weight, &
                                    first_fractional, last_occupied, last_occupied_spin)
   IF (has_weight .OR. first_fractional .OR. last_occupied) THEN
      error stop "Disabled smearing produced an edge status"
   END IF
 
   CALL test_fixed_n_smearing_response()
   CALL test_fixed_spin_populations()
   CALL remove_all_references()
 
CONTAINS
 
! **************************************************************************************************
!> \brief Check fractional, zero, integer-preserving, and unconstrained spin targets.
! **************************************************************************************************
   SUBROUTINE test_fixed_spin_populations()
      REAL(kind=dp), DIMENSION(2)                        :: targets
 
      IF (smearing_separate_spin(-100.0_dp)) error stop "Unconstrained spin was fixed"
      IF (.NOT. smearing_separate_spin(-1.0_dp)) error stop "Integer spin populations were released"
      IF (.NOT. smearing_separate_spin(0.0_dp)) error stop "Zero magnetization was not constrained"
      IF (.NOT. smearing_separate_spin(0.5_dp)) error stop "Fractional magnetization was not constrained"
      targets(:) = smearing_spin_targets([3, 1], -1.0_dp)
      IF (maxval(abs(targets - [3.0_dp, 1.0_dp])) > 1.0e-14_dp) THEN
         error stop "Integer spin targets changed"
      END IF
      targets(:) = smearing_spin_targets([3, 1], 0.0_dp)
      IF (maxval(abs(targets - 2.0_dp)) > 1.0e-14_dp) error stop "Zero-moment targets are wrong"
      targets(:) = smearing_spin_targets([3, 1], 0.5_dp)
      IF (maxval(abs(targets - [2.25_dp, 1.75_dp])) > 1.0e-14_dp) THEN
         error stop "Fractional spin targets were rounded"
      END IF
   END SUBROUTINE test_fixed_spin_populations
 
! **************************************************************************************************
!> \brief Check the fixed-N occupation response for every supported K-point smearing method.
! **************************************************************************************************
   SUBROUTINE test_fixed_n_smearing_response()
      INTEGER, PARAMETER                                 :: nkpoint = 2, nstate = 5
      REAL(kind=dp), PARAMETER                           :: fd_step = 1.0e-6_dp, maxocc = 2.0_dp, &
                                                            sigma = 0.15_dp, TARGET = 4.3_dp
 
      INTEGER                                            :: ikpoint, imethod
      INTEGER, DIMENSION(4)                              :: methods
      REAL(kind=dp)                                      :: error, kts, mu, response_sum
      REAL(kind=dp), DIMENSION(nkpoint)                  :: kpoint_weight
      REAL(kind=dp), DIMENSION(nstate*nkpoint)           :: fd_response, flat_direction, &
                                                            flat_response, predicted_response
      REAL(kind=dp), DIMENSION(nstate, nkpoint) :: direction, eigenvalue, eigenvalue_minus, &
         eigenvalue_plus, occupation, occupation_minus, occupation_plus, response
 
      methods(:) = [smear_fermi_dirac, smear_gaussian, smear_mp, smear_mv]
      kpoint_weight(:) = [0.4_dp, 0.6_dp]
      eigenvalue(:, 1) = [-0.38_dp, -0.14_dp, 0.02_dp, 0.23_dp, 0.52_dp]
      eigenvalue(:, 2) = [-0.31_dp, -0.08_dp, 0.07_dp, 0.28_dp, 0.61_dp]
      direction(:, 1) = [0.13_dp, -0.07_dp, 0.11_dp, -0.04_dp, 0.09_dp]
      direction(:, 2) = [-0.08_dp, 0.05_dp, -0.12_dp, 0.06_dp, -0.03_dp]
      flat_direction(:) = reshape(direction, [nstate*nkpoint])
 
      DO imethod = 1, SIZE(methods)
         CALL smearkp(occupation, mu, kts, eigenvalue, TARGET, kpoint_weight, sigma, maxocc, &
                      methods(imethod))
         DO ikpoint = 1, nkpoint
            CALL smearing_response_weight( &
               response(:, ikpoint), occupation(:, ikpoint), eigenvalue(:, ikpoint), mu, sigma, &
               maxocc, nstate, methods(imethod))
            response(:, ikpoint) = kpoint_weight(ikpoint)*response(:, ikpoint)
         END DO
         response_sum = sum(response)
         IF (abs(response_sum) <= epsilon(response_sum)) THEN
            error stop "Singular fixed-N smearing response in unit test"
         END IF
         IF ((methods(imethod) == smear_mp .OR. methods(imethod) == smear_mv) .AND. &
             .NOT. any(response < 0.0_dp)) THEN
            error stop "Signed smearing-response branch was not exercised"
         END IF
         flat_response = reshape(response, [nstate*nkpoint])
 
         eigenvalue_plus(:, :) = eigenvalue + fd_step*direction
         eigenvalue_minus(:, :) = eigenvalue - fd_step*direction
         CALL smearkp(occupation_plus, mu, kts, eigenvalue_plus, TARGET, kpoint_weight, sigma, &
                      maxocc, methods(imethod))
         CALL smearkp(occupation_minus, mu, kts, eigenvalue_minus, TARGET, kpoint_weight, sigma, &
                      maxocc, methods(imethod))
         DO ikpoint = 1, nkpoint
            fd_response((ikpoint - 1)*nstate + 1:ikpoint*nstate) = &
               kpoint_weight(ikpoint)*(occupation_plus(:, ikpoint) - &
                                       occupation_minus(:, ikpoint))/(2.0_dp*fd_step)
         END DO
         predicted_response = -flat_response* &
                              (flat_direction - dot_product(flat_response, flat_direction)/response_sum)
         error = maxval(abs(fd_response - predicted_response))
         error = max(error, abs(sum(fd_response)))
         ! A common energy shift leaves fixed-N occupations unchanged.
         CALL smearkp(occupation_plus, mu, kts, eigenvalue + 0.37_dp, TARGET, kpoint_weight, &
                      sigma, maxocc, methods(imethod))
         error = max(error, maxval(abs(occupation_plus - occupation)))
         IF (error > 2.0e-6_dp) error stop "Fixed-N smearing response finite difference failed"
      END DO
 
   END SUBROUTINE test_fixed_n_smearing_response
 
END PROGRAM kpoint_smearing_unittest