Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 2 : PROGRAM kpoint_smearing_unittest
9 2 : USE bibliography, ONLY: add_all_references
10 : USE input_constants, ONLY: smear_fermi_dirac,&
11 : smear_gaussian,&
12 : smear_mp,&
13 : smear_mv
14 : USE kinds, ONLY: dp
15 : USE kpoint_methods, ONLY: kpoint_smearing_edge_status
16 : USE qs_ot, ONLY: qs_ot_fixed_n_energy_hessian
17 : USE reference_manager, ONLY: remove_all_references
18 : USE scf_control_types, ONLY: smear_type
19 : USE smearing_utils, ONLY: Smearkp,&
20 : smearing_response_weight
21 :
22 : IMPLICIT NONE
23 :
24 : LOGICAL :: first_fractional, has_weight, last_occupied
25 : LOGICAL, DIMENSION(2) :: last_occupied_spin
26 : REAL(KIND=dp), DIMENSION(2) :: weight
27 : REAL(KIND=dp), DIMENSION(3, 2, 2) :: occupation
28 : TYPE(smear_type) :: smear
29 :
30 2 : CALL add_all_references()
31 2 : smear%do_smear = .TRUE.
32 2 : smear%eps_fermi_dirac = 1.0E-5_dp
33 2 : occupation = 0.0_dp
34 14 : occupation(1, :, :) = 1.0_dp
35 2 : weight = [1.0_dp, 0.0_dp]
36 :
37 2 : occupation(1, 2, 1) = 0.5_dp
38 2 : occupation(3, 2, 2) = 0.5_dp
39 : CALL kpoint_smearing_edge_status(occupation, weight, smear, 2, has_weight, &
40 2 : first_fractional, last_occupied, last_occupied_spin)
41 2 : IF (.NOT. has_weight) ERROR STOP "Positive K-point weight was not detected"
42 2 : IF (first_fractional) ERROR STOP "Zero-weight first-band occupation was included"
43 2 : IF (last_occupied) ERROR STOP "Zero-weight last-band occupation was included"
44 6 : IF (ANY(last_occupied_spin)) ERROR STOP "Zero-weight spin edge was included"
45 :
46 2 : weight(2) = 0.25_dp
47 : CALL kpoint_smearing_edge_status(occupation, weight, smear, 2, has_weight, &
48 2 : first_fractional, last_occupied, last_occupied_spin)
49 2 : IF (.NOT. first_fractional) ERROR STOP "Fractional first band was not detected"
50 2 : IF (.NOT. last_occupied) ERROR STOP "Occupied last band was not detected"
51 2 : IF (last_occupied_spin(1)) ERROR STOP "Wrong occupied spin edge was reported"
52 2 : IF (.NOT. last_occupied_spin(2)) ERROR STOP "Occupied spin edge was not reported"
53 :
54 2 : occupation = 0.0_dp
55 6 : occupation(1, :, 1) = 2.0_dp
56 : CALL kpoint_smearing_edge_status(occupation(:, :, 1:1), weight, smear, 1, has_weight, &
57 2 : first_fractional, last_occupied)
58 2 : IF (first_fractional) ERROR STOP "Restricted full occupation was classified as fractional"
59 :
60 2 : weight = 0.0_dp
61 : CALL kpoint_smearing_edge_status(occupation, weight, smear, 2, has_weight, &
62 2 : first_fractional, last_occupied, last_occupied_spin)
63 2 : IF (has_weight .OR. first_fractional .OR. last_occupied) THEN
64 0 : ERROR STOP "Zero-weight K-point set produced an edge status"
65 : END IF
66 :
67 2 : smear%do_smear = .FALSE.
68 6 : weight = 1.0_dp
69 : CALL kpoint_smearing_edge_status(occupation, weight, smear, 2, has_weight, &
70 2 : first_fractional, last_occupied, last_occupied_spin)
71 2 : IF (has_weight .OR. first_fractional .OR. last_occupied) THEN
72 0 : ERROR STOP "Disabled smearing produced an edge status"
73 : END IF
74 :
75 2 : CALL test_fixed_n_smearing_response()
76 2 : CALL remove_all_references()
77 :
78 : CONTAINS
79 :
80 : ! **************************************************************************************************
81 : !> \brief Check the fixed-N occupation response for every supported K-point smearing method.
82 : ! **************************************************************************************************
83 2 : SUBROUTINE test_fixed_n_smearing_response()
84 : INTEGER, PARAMETER :: nkpoint = 2, nstate = 5
85 : REAL(KIND=dp), PARAMETER :: fd_step = 1.0E-6_dp, maxocc = 2.0_dp, &
86 : sigma = 0.15_dp, TARGET = 4.3_dp
87 :
88 : INTEGER :: ikpoint, imethod
89 : INTEGER, DIMENSION(4) :: methods
90 : REAL(KIND=dp) :: error, kTS, mu, response_sum
91 : REAL(KIND=dp), DIMENSION(nkpoint) :: kpoint_weight
92 : REAL(KIND=dp), DIMENSION(nstate*nkpoint) :: fd_response, flat_direction, &
93 : predicted_response
94 : REAL(KIND=dp), &
95 : DIMENSION(nstate*nkpoint, nstate*nkpoint) :: hessian
96 : REAL(KIND=dp), DIMENSION(nstate, nkpoint) :: direction, eigenvalue, eigenvalue_minus, &
97 : eigenvalue_plus, occupation, occupation_minus, occupation_plus, response
98 :
99 2 : methods(:) = [smear_fermi_dirac, smear_gaussian, smear_mp, smear_mv]
100 2 : kpoint_weight(:) = [0.4_dp, 0.6_dp]
101 12 : eigenvalue(:, 1) = [-0.38_dp, -0.14_dp, 0.02_dp, 0.23_dp, 0.52_dp]
102 12 : eigenvalue(:, 2) = [-0.31_dp, -0.08_dp, 0.07_dp, 0.28_dp, 0.61_dp]
103 12 : direction(:, 1) = [0.13_dp, -0.07_dp, 0.11_dp, -0.04_dp, 0.09_dp]
104 12 : direction(:, 2) = [-0.08_dp, 0.05_dp, -0.12_dp, 0.06_dp, -0.03_dp]
105 2 : flat_direction(:) = RESHAPE(direction, [nstate*nkpoint])
106 :
107 10 : DO imethod = 1, SIZE(methods)
108 : CALL Smearkp(occupation, mu, kTS, eigenvalue, TARGET, kpoint_weight, sigma, maxocc, &
109 8 : methods(imethod))
110 24 : DO ikpoint = 1, nkpoint
111 : CALL smearing_response_weight( &
112 : response(:, ikpoint), occupation(:, ikpoint), eigenvalue(:, ikpoint), mu, sigma, &
113 16 : maxocc, nstate, methods(imethod))
114 104 : response(:, ikpoint) = kpoint_weight(ikpoint)*response(:, ikpoint)
115 : END DO
116 104 : response_sum = SUM(response)
117 8 : IF (ABS(response_sum) <= EPSILON(response_sum)) THEN
118 0 : ERROR STOP "Singular fixed-N smearing response in unit test"
119 : END IF
120 56 : IF ((methods(imethod) == smear_mp .OR. methods(imethod) == smear_mv) .AND. &
121 : .NOT. ANY(response < 0.0_dp)) THEN
122 0 : ERROR STOP "Signed smearing-response branch was not exercised"
123 : END IF
124 : CALL qs_ot_fixed_n_energy_hessian( &
125 8 : RESHAPE(response, [nstate*nkpoint]), response_sum, hessian)
126 :
127 104 : eigenvalue_plus(:, :) = eigenvalue + fd_step*direction
128 104 : eigenvalue_minus(:, :) = eigenvalue - fd_step*direction
129 : CALL Smearkp(occupation_plus, mu, kTS, eigenvalue_plus, TARGET, kpoint_weight, sigma, &
130 8 : maxocc, methods(imethod))
131 : CALL Smearkp(occupation_minus, mu, kTS, eigenvalue_minus, TARGET, kpoint_weight, sigma, &
132 8 : maxocc, methods(imethod))
133 24 : DO ikpoint = 1, nkpoint
134 : fd_response((ikpoint - 1)*nstate + 1:ikpoint*nstate) = &
135 : kpoint_weight(ikpoint)*(occupation_plus(:, ikpoint) - &
136 104 : occupation_minus(:, ikpoint))/(2.0_dp*fd_step)
137 : END DO
138 976 : predicted_response(:) = -MATMUL(hessian, flat_direction)
139 88 : error = MAXVAL(ABS(fd_response - predicted_response))
140 88 : error = MAX(error, ABS(SUM(fd_response)))
141 976 : error = MAX(error, MAXVAL(ABS(MATMUL(hessian, SPREAD(1.0_dp, 1, nstate*nkpoint)))))
142 10 : IF (error > 2.0E-6_dp) ERROR STOP "Fixed-N smearing response finite difference failed"
143 : END DO
144 :
145 2 : END SUBROUTINE test_fixed_n_smearing_response
146 :
147 : END PROGRAM kpoint_smearing_unittest
|
