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 2 : PROGRAM ot_covariant_preconditioner_unittest
8 2 : USE cp_blacs_env, ONLY: cp_blacs_env_create,&
9 : cp_blacs_env_release,&
10 : cp_blacs_env_type
11 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
12 : cp_fm_struct_release,&
13 : cp_fm_struct_type
14 : USE cp_fm_types, ONLY: cp_fm_create,&
15 : cp_fm_get_submatrix,&
16 : cp_fm_release,&
17 : cp_fm_set_submatrix,&
18 : cp_fm_type
19 : USE input_constants, ONLY: ot_precond_full_all_covariant
20 : USE iso_fortran_env, ONLY: error_unit
21 : USE kinds, ONLY: dp
22 : USE message_passing, ONLY: mp_para_env_type,&
23 : mp_world_finalize,&
24 : mp_world_init
25 : USE preconditioner_apply, ONLY: apply_preconditioner_fm
26 : USE preconditioner_types, ONLY: destroy_preconditioner,&
27 : init_preconditioner,&
28 : preconditioner_type
29 :
30 : IMPLICIT NONE
31 :
32 : INTEGER, PARAMETER :: k = 3, n = 5
33 : REAL(KIND=dp), PARAMETER :: eps = 2.0E-11_dp
34 :
35 : INTEGER :: i, j
36 : REAL(KIND=dp), DIMENSION(n, k) :: gradient, gradient_rotated, hessian_x, &
37 : output, output_reference, &
38 : output_rotated, x_reference
39 : REAL(KIND=dp), DIMENSION(n) :: full_evals
40 : REAL(KIND=dp), DIMENSION(k) :: occ_evals
41 : REAL(KIND=dp), DIMENSION(k, k) :: occupied_h, q, rotation, rotation_2, &
42 : rotation_gauge
43 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
44 : TYPE(cp_fm_struct_type), POINTER :: fm_struct
45 : TYPE(cp_fm_type) :: matrix_in, matrix_out
46 : TYPE(mp_para_env_type), POINTER :: para_env
47 8 : TYPE(preconditioner_type) :: preconditioner_env
48 :
49 2 : NULLIFY (blacs_env, fm_struct, para_env)
50 2 : ALLOCATE (para_env)
51 2 : CALL mp_world_init(para_env)
52 2 : CALL cp_blacs_env_create(blacs_env=blacs_env, para_env=para_env)
53 2 : CALL init_preconditioner(preconditioner_env, para_env, blacs_env)
54 :
55 : CALL cp_fm_struct_create(fm_struct, nrow_global=n, ncol_global=k, &
56 2 : context=blacs_env, para_env=para_env)
57 2 : CALL cp_fm_create(matrix_in, fm_struct, name="covariant test input")
58 2 : CALL cp_fm_create(matrix_out, fm_struct, name="covariant test output")
59 2 : CALL cp_fm_struct_release(fm_struct)
60 :
61 : CALL cp_fm_struct_create(fm_struct, nrow_global=n, ncol_global=n, &
62 2 : context=blacs_env, para_env=para_env)
63 2 : ALLOCATE (preconditioner_env%fm)
64 2 : CALL cp_fm_create(preconditioner_env%fm, fm_struct, name="spectral eigenvectors")
65 2 : CALL cp_fm_struct_release(fm_struct)
66 :
67 : CALL cp_fm_struct_create(fm_struct, nrow_global=k, ncol_global=k, &
68 2 : context=blacs_env, para_env=para_env)
69 2 : ALLOCATE (preconditioner_env%occ_rotation)
70 2 : CALL cp_fm_create(preconditioner_env%occ_rotation, fm_struct, name="occupied rotation")
71 2 : CALL cp_fm_struct_release(fm_struct)
72 :
73 2 : preconditioner_env%in_use = ot_precond_full_all_covariant
74 2 : preconditioner_env%energy_gap = 0.10_dp
75 2 : ALLOCATE (preconditioner_env%full_evals(n), preconditioner_env%occ_evals(k))
76 2 : full_evals = [0.8_dp, 1.3_dp, 2.1_dp, 3.0_dp, 4.4_dp]
77 2 : occ_evals = [-0.9_dp, -0.2_dp, 0.35_dp]
78 12 : preconditioner_env%full_evals = full_evals
79 8 : preconditioner_env%occ_evals = occ_evals
80 :
81 2 : CALL set_identity(preconditioner_env%fm, n)
82 2 : CALL make_rotation(0.43_dp, -0.31_dp, rotation)
83 2 : CALL cp_fm_set_submatrix(preconditioner_env%occ_rotation, rotation)
84 :
85 8 : DO j = 1, k
86 38 : DO i = 1, n
87 : gradient(i, j) = SIN(0.37_dp*REAL(2*i + j, dp)) + &
88 30 : 0.2_dp*COS(0.19_dp*REAL(i - 3*j, dp))
89 : x_reference(i, j) = COS(0.23_dp*REAL(i + 2*j, dp)) - &
90 36 : 0.1_dp*SIN(0.41_dp*REAL(3*i - j, dp))
91 : END DO
92 : END DO
93 :
94 : ! Check the production FM wrapper against its explicit spectral formula.
95 : CALL apply_model(gradient, rotation, full_evals, occ_evals, &
96 2 : preconditioner_env%energy_gap, output_reference)
97 2 : CALL apply_fm(preconditioner_env, matrix_in, matrix_out, gradient, output)
98 2 : CALL assert_close(output, output_reference, eps, "Explicit covariant spectral inverse")
99 :
100 : ! Invert the frozen-H Sylvester operator A X - X B exactly.
101 2 : occupied_h = 0.0_dp
102 8 : DO j = 1, k
103 80 : occupied_h = occupied_h + occ_evals(j)*outer_product(rotation(:, j), rotation(:, j))
104 : END DO
105 152 : hessian_x = -MATMUL(x_reference, occupied_h)
106 12 : DO i = 1, n
107 42 : hessian_x(i, :) = hessian_x(i, :) + full_evals(i)*x_reference(i, :)
108 : END DO
109 2 : CALL apply_fm(preconditioner_env, matrix_in, matrix_out, hessian_x, output)
110 2 : CALL assert_close(output, x_reference, eps, "Frozen-H Sylvester inverse")
111 :
112 : ! Rotate the occupied gauge and verify P_(C R)(G R) = P_C(G) R.
113 2 : CALL make_rotation(-0.37_dp, 0.28_dp, rotation_gauge)
114 80 : rotation_2 = MATMUL(TRANSPOSE(rotation_gauge), rotation)
115 116 : gradient_rotated = MATMUL(gradient, rotation_gauge)
116 2 : CALL cp_fm_set_submatrix(preconditioner_env%occ_rotation, rotation_2)
117 2 : CALL apply_fm(preconditioner_env, matrix_in, matrix_out, gradient_rotated, output_rotated)
118 : CALL assert_close(output_rotated, MATMUL(output_reference, rotation_gauge), eps, &
119 116 : "Occupied-gauge covariance")
120 :
121 : ! Eigenvectors may rotate freely inside a degenerate occupied eigenspace.
122 2 : occ_evals = [-0.4_dp, -0.4_dp, 0.2_dp]
123 8 : preconditioner_env%occ_evals = occ_evals
124 2 : CALL cp_fm_set_submatrix(preconditioner_env%occ_rotation, rotation)
125 2 : CALL apply_fm(preconditioner_env, matrix_in, matrix_out, gradient, output)
126 2 : CALL plane_rotation(0.61_dp, 1, 2, q)
127 80 : rotation_2 = MATMUL(rotation, q)
128 2 : CALL cp_fm_set_submatrix(preconditioner_env%occ_rotation, rotation_2)
129 2 : CALL apply_fm(preconditioner_env, matrix_in, matrix_out, gradient, output_rotated)
130 2 : CALL assert_close(output_rotated, output, eps, "Degenerate occupied subspace")
131 :
132 : ! The gap floor must leave the inverse positive definite.
133 2 : preconditioner_env%energy_gap = 0.75_dp
134 8 : preconditioner_env%occ_evals = [0.7_dp, 0.9_dp, 1.1_dp]
135 2 : CALL cp_fm_set_submatrix(preconditioner_env%occ_rotation, rotation)
136 2 : CALL apply_fm(preconditioner_env, matrix_in, matrix_out, gradient, output)
137 38 : IF (SUM(gradient*output) <= 0.0_dp) THEN
138 0 : ERROR STOP "Gap-floor inverse is not positive definite"
139 : END IF
140 :
141 2 : CALL cp_fm_release(matrix_in)
142 2 : CALL cp_fm_release(matrix_out)
143 2 : CALL destroy_preconditioner(preconditioner_env)
144 2 : CALL cp_blacs_env_release(blacs_env)
145 2 : CALL mp_world_finalize()
146 4 : DEALLOCATE (para_env)
147 :
148 : CONTAINS
149 :
150 : ! **************************************************************************************************
151 : !> \brief Apply the production FM path to a replicated test matrix.
152 : !> \param preconditioner_env ...
153 : !> \param matrix_in ...
154 : !> \param matrix_out ...
155 : !> \param input ...
156 : !> \param RESULT ...
157 : ! **************************************************************************************************
158 12 : SUBROUTINE apply_fm(preconditioner_env, matrix_in, matrix_out, input, RESULT)
159 :
160 : TYPE(preconditioner_type) :: preconditioner_env
161 : TYPE(cp_fm_type), INTENT(IN) :: matrix_in, matrix_out
162 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: input
163 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: result
164 :
165 12 : CALL cp_fm_set_submatrix(matrix_in, input)
166 12 : CALL apply_preconditioner_fm(preconditioner_env, matrix_in, matrix_out)
167 12 : CALL cp_fm_get_submatrix(matrix_out, RESULT)
168 :
169 2 : END SUBROUTINE apply_fm
170 :
171 : ! **************************************************************************************************
172 : !> \brief Explicit reference for the covariant state-selective spectral inverse.
173 : !> \param input ...
174 : !> \param rotation ...
175 : !> \param full_evals ...
176 : !> \param occ_evals ...
177 : !> \param gap ...
178 : !> \param RESULT ...
179 : ! **************************************************************************************************
180 2 : SUBROUTINE apply_model(input, rotation, full_evals, occ_evals, gap, RESULT)
181 :
182 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: input, rotation
183 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: full_evals, occ_evals
184 : REAL(KIND=dp), INTENT(IN) :: gap
185 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: result
186 :
187 : INTEGER :: i, j
188 : REAL(KIND=dp), &
189 4 : DIMENSION(SIZE(input, 1), SIZE(input, 2)) :: canonical
190 :
191 156 : canonical = MATMUL(input, rotation)
192 8 : DO j = 1, SIZE(input, 2)
193 38 : DO i = 1, SIZE(input, 1)
194 36 : canonical(i, j) = canonical(i, j)/MAX(gap, full_evals(i) - occ_evals(j))
195 : END DO
196 : END DO
197 154 : RESULT = MATMUL(canonical, TRANSPOSE(rotation))
198 :
199 2 : END SUBROUTINE apply_model
200 :
201 : ! **************************************************************************************************
202 : !> \brief Create a product of two plane rotations.
203 : !> \param angle_12 ...
204 : !> \param angle_23 ...
205 : !> \param rotation ...
206 : ! **************************************************************************************************
207 4 : SUBROUTINE make_rotation(angle_12, angle_23, rotation)
208 :
209 : REAL(KIND=dp), INTENT(IN) :: angle_12, angle_23
210 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: rotation
211 :
212 : REAL(KIND=dp), &
213 8 : DIMENSION(SIZE(rotation, 1), SIZE(rotation, 2)) :: r12, r23
214 :
215 4 : CALL plane_rotation(angle_12, 1, 2, r12)
216 4 : CALL plane_rotation(angle_23, 2, 3, r23)
217 216 : rotation = MATMUL(r12, r23)
218 :
219 4 : END SUBROUTINE make_rotation
220 :
221 : ! **************************************************************************************************
222 : !> \brief Create an orthogonal plane rotation.
223 : !> \param angle ...
224 : !> \param axis_1 ...
225 : !> \param axis_2 ...
226 : !> \param rotation ...
227 : ! **************************************************************************************************
228 10 : SUBROUTINE plane_rotation(angle, axis_1, axis_2, rotation)
229 :
230 : REAL(KIND=dp), INTENT(IN) :: angle
231 : INTEGER, INTENT(IN) :: axis_1, axis_2
232 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: rotation
233 :
234 : INTEGER :: i
235 :
236 130 : rotation = 0.0_dp
237 40 : DO i = 1, SIZE(rotation, 1)
238 40 : rotation(i, i) = 1.0_dp
239 : END DO
240 10 : rotation(axis_1, axis_1) = COS(angle)
241 10 : rotation(axis_2, axis_2) = COS(angle)
242 10 : rotation(axis_1, axis_2) = -SIN(angle)
243 10 : rotation(axis_2, axis_1) = SIN(angle)
244 :
245 10 : END SUBROUTINE plane_rotation
246 :
247 : ! **************************************************************************************************
248 : !> \brief Set a distributed full matrix to the identity.
249 : !> \param matrix ...
250 : !> \param n ...
251 : ! **************************************************************************************************
252 2 : SUBROUTINE set_identity(matrix, n)
253 :
254 : TYPE(cp_fm_type), INTENT(IN) :: matrix
255 : INTEGER, INTENT(IN) :: n
256 :
257 : INTEGER :: i
258 2 : REAL(KIND=dp), DIMENSION(n, n) :: identity
259 :
260 62 : identity = 0.0_dp
261 12 : DO i = 1, n
262 12 : identity(i, i) = 1.0_dp
263 : END DO
264 2 : CALL cp_fm_set_submatrix(matrix, identity)
265 :
266 2 : END SUBROUTINE set_identity
267 :
268 : ! **************************************************************************************************
269 : !> \brief Form a real outer product.
270 : !> \param left ...
271 : !> \param right ...
272 : !> \return ...
273 : ! **************************************************************************************************
274 6 : PURE FUNCTION outer_product(left, right) RESULT(product)
275 :
276 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: left, right
277 : REAL(KIND=dp), DIMENSION(SIZE(left), SIZE(right)) :: product
278 :
279 : INTEGER :: i, j
280 :
281 24 : DO j = 1, SIZE(right)
282 78 : DO i = 1, SIZE(left)
283 72 : PRODUCT(i, j) = left(i)*right(j)
284 : END DO
285 : END DO
286 :
287 6 : END FUNCTION outer_product
288 :
289 : ! **************************************************************************************************
290 : !> \brief Abort when two matrices differ beyond the requested tolerance.
291 : !> \param actual ...
292 : !> \param reference ...
293 : !> \param tolerance ...
294 : !> \param label ...
295 : ! **************************************************************************************************
296 8 : SUBROUTINE assert_close(actual, reference, tolerance, label)
297 :
298 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: actual, reference
299 : REAL(KIND=dp), INTENT(IN) :: tolerance
300 : CHARACTER(LEN=*), INTENT(IN) :: label
301 :
302 : INTEGER :: io_unit
303 :
304 152 : IF (MAXVAL(ABS(actual - reference)) > tolerance) THEN
305 0 : io_unit = error_unit
306 0 : WRITE (io_unit, '(A)') TRIM(label)
307 0 : ERROR STOP "Matrix comparison failed"
308 : END IF
309 :
310 8 : END SUBROUTINE assert_close
311 :
312 : END PROGRAM ot_covariant_preconditioner_unittest
|
