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 : smearing_separate_spin,&
22 : smearing_spin_targets
23 :
24 : IMPLICIT NONE
25 :
26 : LOGICAL :: first_fractional, has_weight, last_occupied
27 : LOGICAL, DIMENSION(2) :: last_occupied_spin
28 : REAL(KIND=dp), DIMENSION(2) :: weight
29 : REAL(KIND=dp), DIMENSION(3, 2, 2) :: occupation
30 : TYPE(smear_type) :: smear
31 :
32 2 : CALL add_all_references()
33 2 : smear%do_smear = .TRUE.
34 2 : smear%eps_fermi_dirac = 1.0E-5_dp
35 2 : occupation = 0.0_dp
36 14 : occupation(1, :, :) = 1.0_dp
37 2 : weight = [1.0_dp, 0.0_dp]
38 :
39 2 : occupation(1, 2, 1) = 0.5_dp
40 2 : occupation(3, 2, 2) = 0.5_dp
41 : CALL kpoint_smearing_edge_status(occupation, weight, smear, 2, has_weight, &
42 2 : first_fractional, last_occupied, last_occupied_spin)
43 2 : IF (.NOT. has_weight) ERROR STOP "Positive K-point weight was not detected"
44 2 : IF (first_fractional) ERROR STOP "Zero-weight first-band occupation was included"
45 2 : IF (last_occupied) ERROR STOP "Zero-weight last-band occupation was included"
46 6 : IF (ANY(last_occupied_spin)) ERROR STOP "Zero-weight spin edge was included"
47 :
48 2 : weight(2) = 0.25_dp
49 : CALL kpoint_smearing_edge_status(occupation, weight, smear, 2, has_weight, &
50 2 : first_fractional, last_occupied, last_occupied_spin)
51 2 : IF (.NOT. first_fractional) ERROR STOP "Fractional first band was not detected"
52 2 : IF (.NOT. last_occupied) ERROR STOP "Occupied last band was not detected"
53 2 : IF (last_occupied_spin(1)) ERROR STOP "Wrong occupied spin edge was reported"
54 2 : IF (.NOT. last_occupied_spin(2)) ERROR STOP "Occupied spin edge was not reported"
55 :
56 2 : occupation = 0.0_dp
57 6 : occupation(1, :, 1) = 2.0_dp
58 : CALL kpoint_smearing_edge_status(occupation(:, :, 1:1), weight, smear, 1, has_weight, &
59 2 : first_fractional, last_occupied)
60 2 : IF (first_fractional) ERROR STOP "Restricted full occupation was classified as fractional"
61 :
62 2 : weight = 0.0_dp
63 : CALL kpoint_smearing_edge_status(occupation, weight, smear, 2, has_weight, &
64 2 : first_fractional, last_occupied, last_occupied_spin)
65 2 : IF (has_weight .OR. first_fractional .OR. last_occupied) THEN
66 0 : ERROR STOP "Zero-weight K-point set produced an edge status"
67 : END IF
68 :
69 2 : smear%do_smear = .FALSE.
70 6 : weight = 1.0_dp
71 : CALL kpoint_smearing_edge_status(occupation, weight, smear, 2, has_weight, &
72 2 : first_fractional, last_occupied, last_occupied_spin)
73 2 : IF (has_weight .OR. first_fractional .OR. last_occupied) THEN
74 0 : ERROR STOP "Disabled smearing produced an edge status"
75 : END IF
76 :
77 2 : CALL test_fixed_n_smearing_response()
78 2 : CALL test_fixed_spin_populations()
79 2 : CALL remove_all_references()
80 :
81 : CONTAINS
82 :
83 : ! **************************************************************************************************
84 : !> \brief Check fractional, zero, integer-preserving, and unconstrained spin targets.
85 : ! **************************************************************************************************
86 2 : SUBROUTINE test_fixed_spin_populations()
87 : REAL(KIND=dp), DIMENSION(2) :: targets
88 :
89 2 : IF (smearing_separate_spin(-100.0_dp)) ERROR STOP "Unconstrained spin was fixed"
90 2 : IF (.NOT. smearing_separate_spin(-1.0_dp)) ERROR STOP "Integer spin populations were released"
91 2 : IF (.NOT. smearing_separate_spin(0.0_dp)) ERROR STOP "Zero magnetization was not constrained"
92 2 : IF (.NOT. smearing_separate_spin(0.5_dp)) ERROR STOP "Fractional magnetization was not constrained"
93 2 : targets(:) = smearing_spin_targets([3, 1], -1.0_dp)
94 6 : IF (MAXVAL(ABS(targets - [3.0_dp, 1.0_dp])) > 1.0E-14_dp) THEN
95 0 : ERROR STOP "Integer spin targets changed"
96 : END IF
97 2 : targets(:) = smearing_spin_targets([3, 1], 0.0_dp)
98 6 : IF (MAXVAL(ABS(targets - 2.0_dp)) > 1.0E-14_dp) ERROR STOP "Zero-moment targets are wrong"
99 2 : targets(:) = smearing_spin_targets([3, 1], 0.5_dp)
100 6 : IF (MAXVAL(ABS(targets - [2.25_dp, 1.75_dp])) > 1.0E-14_dp) THEN
101 0 : ERROR STOP "Fractional spin targets were rounded"
102 : END IF
103 2 : END SUBROUTINE test_fixed_spin_populations
104 :
105 : ! **************************************************************************************************
106 : !> \brief Check the fixed-N occupation response for every supported K-point smearing method.
107 : ! **************************************************************************************************
108 2 : SUBROUTINE test_fixed_n_smearing_response()
109 : INTEGER, PARAMETER :: nkpoint = 2, nstate = 5
110 : REAL(KIND=dp), PARAMETER :: fd_step = 1.0E-6_dp, maxocc = 2.0_dp, &
111 : sigma = 0.15_dp, TARGET = 4.3_dp
112 :
113 : INTEGER :: ikpoint, imethod
114 : INTEGER, DIMENSION(4) :: methods
115 : REAL(KIND=dp) :: error, kTS, mu, response_sum
116 : REAL(KIND=dp), DIMENSION(nkpoint) :: kpoint_weight
117 : REAL(KIND=dp), DIMENSION(nstate*nkpoint) :: fd_response, flat_direction, &
118 : predicted_response
119 : REAL(KIND=dp), &
120 : DIMENSION(nstate*nkpoint, nstate*nkpoint) :: hessian
121 : REAL(KIND=dp), DIMENSION(nstate, nkpoint) :: direction, eigenvalue, eigenvalue_minus, &
122 : eigenvalue_plus, occupation, occupation_minus, occupation_plus, response
123 :
124 2 : methods(:) = [smear_fermi_dirac, smear_gaussian, smear_mp, smear_mv]
125 2 : kpoint_weight(:) = [0.4_dp, 0.6_dp]
126 12 : eigenvalue(:, 1) = [-0.38_dp, -0.14_dp, 0.02_dp, 0.23_dp, 0.52_dp]
127 12 : eigenvalue(:, 2) = [-0.31_dp, -0.08_dp, 0.07_dp, 0.28_dp, 0.61_dp]
128 12 : direction(:, 1) = [0.13_dp, -0.07_dp, 0.11_dp, -0.04_dp, 0.09_dp]
129 12 : direction(:, 2) = [-0.08_dp, 0.05_dp, -0.12_dp, 0.06_dp, -0.03_dp]
130 2 : flat_direction(:) = RESHAPE(direction, [nstate*nkpoint])
131 :
132 10 : DO imethod = 1, SIZE(methods)
133 : CALL Smearkp(occupation, mu, kTS, eigenvalue, TARGET, kpoint_weight, sigma, maxocc, &
134 8 : methods(imethod))
135 24 : DO ikpoint = 1, nkpoint
136 : CALL smearing_response_weight( &
137 : response(:, ikpoint), occupation(:, ikpoint), eigenvalue(:, ikpoint), mu, sigma, &
138 16 : maxocc, nstate, methods(imethod))
139 104 : response(:, ikpoint) = kpoint_weight(ikpoint)*response(:, ikpoint)
140 : END DO
141 104 : response_sum = SUM(response)
142 8 : IF (ABS(response_sum) <= EPSILON(response_sum)) THEN
143 0 : ERROR STOP "Singular fixed-N smearing response in unit test"
144 : END IF
145 56 : IF ((methods(imethod) == smear_mp .OR. methods(imethod) == smear_mv) .AND. &
146 : .NOT. ANY(response < 0.0_dp)) THEN
147 0 : ERROR STOP "Signed smearing-response branch was not exercised"
148 : END IF
149 : CALL qs_ot_fixed_n_energy_hessian( &
150 8 : RESHAPE(response, [nstate*nkpoint]), response_sum, hessian)
151 :
152 104 : eigenvalue_plus(:, :) = eigenvalue + fd_step*direction
153 104 : eigenvalue_minus(:, :) = eigenvalue - fd_step*direction
154 : CALL Smearkp(occupation_plus, mu, kTS, eigenvalue_plus, TARGET, kpoint_weight, sigma, &
155 8 : maxocc, methods(imethod))
156 : CALL Smearkp(occupation_minus, mu, kTS, eigenvalue_minus, TARGET, kpoint_weight, sigma, &
157 8 : maxocc, methods(imethod))
158 24 : DO ikpoint = 1, nkpoint
159 : fd_response((ikpoint - 1)*nstate + 1:ikpoint*nstate) = &
160 : kpoint_weight(ikpoint)*(occupation_plus(:, ikpoint) - &
161 104 : occupation_minus(:, ikpoint))/(2.0_dp*fd_step)
162 : END DO
163 976 : predicted_response(:) = -MATMUL(hessian, flat_direction)
164 88 : error = MAXVAL(ABS(fd_response - predicted_response))
165 88 : error = MAX(error, ABS(SUM(fd_response)))
166 976 : error = MAX(error, MAXVAL(ABS(MATMUL(hessian, SPREAD(1.0_dp, 1, nstate*nkpoint)))))
167 10 : IF (error > 2.0E-6_dp) ERROR STOP "Fixed-N smearing response finite difference failed"
168 : END DO
169 :
170 2 : END SUBROUTINE test_fixed_n_smearing_response
171 :
172 : END PROGRAM kpoint_smearing_unittest
|
