(git:6344619)
Loading...
Searching...
No Matches
qs_ot_complex_ref_unittest.F
Go to the documentation of this file.
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! **************************************************************************************************
9!> \brief Focused regressions for complex k-point OT derivatives.
10! **************************************************************************************************
11PROGRAM qs_ot_complex_ref_unittest
12 USE bibliography, ONLY: add_all_references
13 USE cp_blacs_env, ONLY: cp_blacs_env_create,&
14 cp_blacs_env_release,&
15 cp_blacs_env_type
16 USE cp_dbcsr_api, ONLY: &
17 dbcsr_create, dbcsr_distribution_new, dbcsr_distribution_release, dbcsr_distribution_type, &
18 dbcsr_finalize, dbcsr_finalize_lib, dbcsr_get_block_p, dbcsr_init_lib, dbcsr_put_block, &
19 dbcsr_release, dbcsr_reserve_blocks, dbcsr_type, dbcsr_type_no_symmetry
20 USE cp_log_handling, ONLY: cp_add_default_logger,&
21 cp_logger_create,&
22 cp_logger_release,&
23 cp_logger_type,&
24 cp_rm_default_logger
25 USE ieee_arithmetic, ONLY: ieee_is_finite
26 USE kinds, ONLY: dp
27 USE machine, ONLY: default_output_unit
28 USE mathlib, ONLY: diag_complex
29 USE message_passing, ONLY: mp_comm_self,&
30 mp_comm_type,&
31 mp_para_env_release,&
32 mp_para_env_type,&
33 mp_world_finalize,&
34 mp_world_init
35 USE qs_ot, ONLY: qs_ot_apply_complex_frechet_dbcsr,&
36 qs_ot_complex_exp_frechet_kernel,&
37 qs_ot_density_tangent,&
38 qs_ot_get_derivative_complex,&
39 qs_ot_get_derivative_ref_complex,&
40 qs_ot_get_orbitals_ref_complex,&
41 qs_ot_get_p_complex
42 USE qs_ot_types, ONLY: qs_ot_type
43 USE reference_manager, ONLY: remove_all_references
44
45 IMPLICIT NONE
46
47 INTEGER, PARAMETER :: k = 2, n = 3
48
49 INTEGER :: io_unit, mynode, nfail
50 INTEGER, DIMENSION(:), POINTER :: col_dist, col_k, row_dist, row_k, row_n
51 INTEGER, DIMENSION(:, :), POINTER :: pgrid
52 TYPE(cp_logger_type), POINTER :: logger
53 TYPE(dbcsr_distribution_type) :: dist
54 TYPE(mp_comm_type) :: mp_comm
55 TYPE(mp_para_env_type), POINTER :: para_env
56
57 nfail = 0
58 NULLIFY (logger, para_env)
59 CALL mp_world_init(mp_comm)
60 mynode = mp_comm%mepos
61 io_unit = -1
62 IF (mynode == 0) io_unit = default_output_unit
63 CALL test_density_tangent(mynode, nfail)
64 ALLOCATE (para_env)
65 CALL para_env%from_dup(mp_comm)
66 CALL cp_logger_create(logger, para_env=para_env, default_global_unit_nr=io_unit, &
67 close_global_unit_on_dealloc=.false.)
68 CALL cp_add_default_logger(logger)
69 CALL add_all_references()
70 CALL dbcsr_init_lib(mp_comm%get_handle(), io_unit)
71 ! Sparse matrix rotations and derivative pullbacks.
72 CALL test_sparse_frechet_patterns(para_env, nfail)
73
74 ALLOCATE (pgrid(0:0, 0:0), row_dist(1), col_dist(1), row_n(1), row_k(1), col_k(1))
75 pgrid(:, :) = 0
76 row_dist(:) = 0
77 col_dist(:) = 0
78 row_n(:) = n
79 row_k(:) = k
80 col_k(:) = k
81 CALL dbcsr_distribution_new(dist, group=mp_comm_self%get_handle(), pgrid=pgrid, &
82 row_dist=row_dist, col_dist=col_dist)
83 ! Both derivative fixtures use one block on MPI_COMM_SELF.
84 CALL test_complex_ref_rotation(dist, row_n, row_k, col_k, para_env, nfail)
85 CALL test_complex_strict_derivative(dist, row_n, row_k, col_k, para_env, nfail)
86 CALL dbcsr_distribution_release(dist)
87 DEALLOCATE (pgrid, row_dist, col_dist, row_n, row_k, col_k)
88
89 CALL dbcsr_finalize_lib()
90 CALL remove_all_references()
91 CALL cp_rm_default_logger()
92 CALL cp_logger_release(logger)
93 CALL mp_para_env_release(para_env)
94 CALL mp_world_finalize()
95
96 IF (nfail > 0) error stop "qs_ot_complex_ref_unittest failed"
97
98CONTAINS
99
100! **************************************************************************************************
101!> \brief Check the coupled finite-chart density tangent against a central finite difference.
102!> \param mynode MPI rank
103!> \param nfail accumulated number of failures
104! **************************************************************************************************
105 SUBROUTINE test_density_tangent(mynode, nfail)
106 INTEGER, INTENT(IN) :: mynode
107 INTEGER, INTENT(INOUT) :: nfail
108
109 INTEGER, PARAMETER :: nbands = 3, nrotation = 6
110 REAL(KIND=dp), PARAMETER :: fd_step = 2.0e-5_dp, &
111 kpoint_weight = 0.625_dp
112
113 COMPLEX(KIND=dp), DIMENSION(nbands, nbands) :: density_minus, density_plus, density_tangent, &
114 direction_generator, generator, rotation, rotation_minus, rotation_plus, weighted_rotation
115 INTEGER :: i, j, r
116 REAL(KIND=dp) :: error, trace_error
117 REAL(KIND=dp), DIMENSION(nbands) :: occupation, weighted_occupation_step
118 REAL(KIND=dp), DIMENSION(nrotation) :: rotation_step
119
120 generator(:, :) = cmplx(0.0_dp, 0.0_dp, kind=dp)
121 generator(1, 2) = cmplx(0.17_dp, -0.09_dp, kind=dp)
122 generator(2, 1) = -conjg(generator(1, 2))
123 generator(1, 3) = cmplx(-0.06_dp, 0.11_dp, kind=dp)
124 generator(3, 1) = -conjg(generator(1, 3))
125 generator(2, 3) = cmplx(0.08_dp, 0.04_dp, kind=dp)
126 generator(3, 2) = -conjg(generator(2, 3))
127 occupation(:) = [1.72_dp, 0.83_dp, 0.31_dp]
128 rotation_step(:) = [0.21_dp, -0.14_dp, -0.08_dp, 0.17_dp, 0.11_dp, 0.06_dp]
129 weighted_occupation_step(:) = [-0.037_dp, 0.052_dp, -0.015_dp]
130
131 CALL qs_ot_density_tangent( &
132 generator, occupation, kpoint_weight, rotation_step, weighted_occupation_step, &
133 density_tangent, difference_step=fd_step)
134
135 direction_generator(:, :) = cmplx(0.0_dp, 0.0_dp, kind=dp)
136 r = 0
137 DO i = 1, nbands - 1
138 DO j = i + 1, nbands
139 r = r + 1
140 direction_generator(i, j) = cmplx(rotation_step(r), 0.0_dp, kind=dp)
141 direction_generator(j, i) = -direction_generator(i, j)
142 r = r + 1
143 direction_generator(i, j) = direction_generator(i, j) + &
144 cmplx(0.0_dp, rotation_step(r), kind=dp)
145 direction_generator(j, i) = direction_generator(j, i) + &
146 cmplx(0.0_dp, rotation_step(r), kind=dp)
147 END DO
148 END DO
149 rotation = dense_antihermitian_exp(generator)
150 rotation_plus = dense_antihermitian_exp(generator + fd_step*direction_generator)
151 rotation_minus = dense_antihermitian_exp(generator - fd_step*direction_generator)
152 weighted_rotation(:, :) = rotation_plus
153 DO j = 1, nbands
154 weighted_rotation(:, j) = &
155 (kpoint_weight*occupation(j) + fd_step*weighted_occupation_step(j))* &
156 weighted_rotation(:, j)
157 END DO
158 density_plus(:, :) = matmul(weighted_rotation, conjg(transpose(rotation_plus)))
159 weighted_rotation(:, :) = rotation_minus
160 DO j = 1, nbands
161 weighted_rotation(:, j) = &
162 (kpoint_weight*occupation(j) - fd_step*weighted_occupation_step(j))* &
163 weighted_rotation(:, j)
164 END DO
165 density_minus(:, :) = matmul(weighted_rotation, conjg(transpose(rotation_minus)))
166 density_plus(:, :) = matmul(conjg(transpose(rotation)), &
167 matmul((density_plus - density_minus)/(2.0_dp*fd_step), rotation))
168 error = maxval(abs(density_tangent - density_plus))
169 trace_error = abs(real(sum([(density_tangent(i, i), i=1, nbands)]), kind=dp) - &
170 sum(weighted_occupation_step))
171
172 IF (error > 3.0e-9_dp .OR. trace_error > 3.0e-10_dp) nfail = nfail + 1
173 IF (mynode == 0) WRITE (io_unit, '(A,2(1X,ES13.6))') &
174 "finite complex density tangent errors:", error, trace_error
175
176 END SUBROUTINE test_density_tangent
177
178! **************************************************************************************************
179!> \brief Check the combined finite complex REF and rotation derivatives.
180!> \param dist single-block matrix distribution
181!> \param row_n AO block sizes
182!> \param row_k occupied row block sizes
183!> \param col_k occupied column block sizes
184!> \param para_env parallel environment
185!> \param nfail accumulated number of failures
186! **************************************************************************************************
187 SUBROUTINE test_complex_ref_rotation(dist, row_n, row_k, col_k, para_env, nfail)
188 TYPE(dbcsr_distribution_type), INTENT(IN) :: dist
189 INTEGER, DIMENSION(:), POINTER :: row_n, row_k, col_k
190 TYPE(mp_para_env_type), POINTER :: para_env
191 INTEGER, INTENT(INOUT) :: nfail
192
193 COMPLEX(KIND=dp), DIMENSION(k, k) :: direction_a, generator, gradient_a, &
194 inverse_sqrt, rotation
195 COMPLEX(KIND=dp), DIMENSION(n, k) :: c_dense, c_reference, direction_x, &
196 gradient_x, hc_dense, q_dense, x_dense
197 COMPLEX(KIND=dp), DIMENSION(n, n) :: hamiltonian, overlap
198 INTEGER :: i, j
199 REAL(KIND=dp) :: error_c, error_ref, error_rot, fd_ref, &
200 fd_rot, hstep, predicted_ref, &
201 predicted_rot
202 REAL(KIND=dp), DIMENSION(k) :: occupation
203 REAL(KIND=dp), DIMENSION(k), TARGET :: rotation_evals
204 TYPE(cp_blacs_env_type), POINTER :: blacs_env
205 TYPE(dbcsr_type), DIMENSION(10), TARGET :: rot_matrix
206 TYPE(dbcsr_type), DIMENSION(2), TARGET :: nk_buffer
207 TYPE(dbcsr_type), DIMENSION(3), TARGET :: nonsym_buffer
208 TYPE(dbcsr_type), DIMENSION(4), TARGET :: sym_buffer
209 TYPE(dbcsr_type), POINTER :: c_im_p, c_re_p, hc_im_p, hc_re_p, &
210 s_im_p, s_re_p
211 TYPE(dbcsr_type), TARGET :: c_im, c_re, f_im, f_re, g_im, g_re, &
212 hc_im, hc_re, s_im, s_re, sx_im, &
213 sx_re, x_im, x_re
214 TYPE(mp_para_env_type), POINTER :: local_para_env
215 TYPE(qs_ot_type) :: rotation_env
216
217 NULLIFY (blacs_env, c_im_p, c_re_p, hc_im_p, hc_re_p, local_para_env, &
218 s_im_p, s_re_p)
219
220 x_dense(:, 1) = [cmplx(1.05_dp, 0.12_dp, kind=dp), cmplx(0.18_dp, -0.27_dp, kind=dp), &
221 cmplx(-0.14_dp, 0.16_dp, kind=dp)]
222 x_dense(:, 2) = [cmplx(0.26_dp, 0.31_dp, kind=dp), cmplx(0.94_dp, -0.17_dp, kind=dp), &
223 cmplx(0.21_dp, 0.13_dp, kind=dp)]
224 direction_x(:, 1) = [cmplx(0.11_dp, -0.07_dp, kind=dp), cmplx(-0.16_dp, 0.09_dp, kind=dp), &
225 cmplx(0.05_dp, 0.14_dp, kind=dp)]
226 direction_x(:, 2) = [cmplx(-0.08_dp, 0.15_dp, kind=dp), cmplx(0.13_dp, 0.04_dp, kind=dp), &
227 cmplx(-0.10_dp, -0.06_dp, kind=dp)]
228
229 generator(:, :) = cmplx(0.0_dp, 0.0_dp, kind=dp)
230 generator(1, 1) = cmplx(0.0_dp, 0.19_dp, kind=dp)
231 generator(2, 2) = cmplx(0.0_dp, -0.08_dp, kind=dp)
232 generator(1, 2) = cmplx(0.23_dp, -0.11_dp, kind=dp)
233 generator(2, 1) = -conjg(generator(1, 2))
234 direction_a(:, :) = cmplx(0.0_dp, 0.0_dp, kind=dp)
235 direction_a(1, 1) = cmplx(0.0_dp, -0.07_dp, kind=dp)
236 direction_a(2, 2) = cmplx(0.0_dp, 0.05_dp, kind=dp)
237 direction_a(1, 2) = cmplx(-0.09_dp, 0.06_dp, kind=dp)
238 direction_a(2, 1) = -conjg(direction_a(1, 2))
239
240 hamiltonian(:, :) = cmplx(0.0_dp, 0.0_dp, kind=dp)
241 hamiltonian(1, 1) = cmplx(0.65_dp, 0.0_dp, kind=dp)
242 hamiltonian(2, 2) = cmplx(-0.35_dp, 0.0_dp, kind=dp)
243 hamiltonian(3, 3) = cmplx(1.15_dp, 0.0_dp, kind=dp)
244 hamiltonian(1, 2) = cmplx(0.18_dp, 0.12_dp, kind=dp)
245 hamiltonian(1, 3) = cmplx(-0.09_dp, 0.07_dp, kind=dp)
246 hamiltonian(2, 3) = cmplx(0.14_dp, -0.17_dp, kind=dp)
247 DO i = 1, n
248 DO j = i + 1, n
249 hamiltonian(j, i) = conjg(hamiltonian(i, j))
250 END DO
251 END DO
252 overlap(:, :) = cmplx(0.0_dp, 0.0_dp, kind=dp)
253 DO i = 1, n
254 overlap(i, i) = cmplx(1.0_dp, 0.0_dp, kind=dp)
255 END DO
256 occupation = [1.0_dp, 0.37_dp]
257
258 CALL dense_inverse_sqrt(matmul(conjg(transpose(x_dense)), x_dense), inverse_sqrt)
259 q_dense = matmul(x_dense, inverse_sqrt)
260 rotation = dense_antihermitian_exp(generator)
261 c_reference = matmul(q_dense, rotation)
262
263 ALLOCATE (local_para_env)
264 CALL local_para_env%from_dup(mp_comm_self)
265 CALL cp_blacs_env_create(blacs_env=blacs_env, para_env=local_para_env)
266
267 CALL create_one_block(x_re, 'rotation_ref_x_re', dist, row_n, col_k)
268 CALL create_one_block(x_im, 'rotation_ref_x_im', dist, row_n, col_k)
269 CALL create_one_block(sx_re, 'rotation_ref_sx_re', dist, row_n, col_k)
270 CALL create_one_block(sx_im, 'rotation_ref_sx_im', dist, row_n, col_k)
271 CALL create_one_block(c_re, 'rotation_ref_c_re', dist, row_n, col_k)
272 CALL create_one_block(c_im, 'rotation_ref_c_im', dist, row_n, col_k)
273 CALL create_one_block(hc_re, 'rotation_ref_hc_re', dist, row_n, col_k)
274 CALL create_one_block(hc_im, 'rotation_ref_hc_im', dist, row_n, col_k)
275 CALL create_one_block(g_re, 'rotation_ref_g_re', dist, row_n, col_k)
276 CALL create_one_block(g_im, 'rotation_ref_g_im', dist, row_n, col_k)
277 CALL create_one_block(s_re, 'rotation_ref_s_re', dist, row_n, row_n)
278 CALL create_one_block(s_im, 'rotation_ref_s_im', dist, row_n, row_n)
279 CALL create_one_block(f_re, 'rotation_ref_f_re', dist, row_k, col_k)
280 CALL create_one_block(f_im, 'rotation_ref_f_im', dist, row_k, col_k)
281 DO i = 1, SIZE(sym_buffer)
282 CALL create_one_block(sym_buffer(i), 'rotation_ref_sym', dist, row_k, col_k)
283 END DO
284 DO i = 1, SIZE(nonsym_buffer)
285 CALL create_one_block(nonsym_buffer(i), 'rotation_ref_nonsym', dist, row_k, col_k)
286 END DO
287 DO i = 1, SIZE(nk_buffer)
288 CALL create_one_block(nk_buffer(i), 'rotation_ref_nk', dist, row_n, col_k)
289 END DO
290 DO i = 1, SIZE(rot_matrix)
291 CALL create_one_block(rot_matrix(i), 'rotation_ref_rot', dist, row_k, col_k)
292 END DO
293
294 CALL put_complex_pair(x_re, x_im, x_dense)
295 CALL put_complex_pair(s_re, s_im, overlap)
296 CALL put_complex_pair(rot_matrix(3), rot_matrix(4), generator)
297
298 rotation_env%has_complex_kpoint_state = .true.
299 rotation_env%settings%do_rotation = .true.
300 rotation_env%settings%eps_irac_filter_matrix = 0.0_dp
301 rotation_env%settings%ot_method = 'LBFG'
302 rotation_env%para_env => local_para_env
303 rotation_env%blacs_env => blacs_env
304 rotation_env%matrix_x => x_re
305 rotation_env%matrix_x_im => x_im
306 rotation_env%matrix_sx => sx_re
307 rotation_env%matrix_sx_im => sx_im
308 rotation_env%matrix_gx => g_re
309 rotation_env%matrix_gx_im => g_im
310 rotation_env%matrix_ref_inv_sqrt => f_re
311 rotation_env%matrix_ref_inv_sqrt_im => f_im
312 rotation_env%buf1_k_k_sym => sym_buffer(1)
313 rotation_env%buf2_k_k_sym => sym_buffer(2)
314 rotation_env%buf3_k_k_sym => sym_buffer(3)
315 rotation_env%buf4_k_k_sym => sym_buffer(4)
316 rotation_env%buf1_k_k_nosym => nonsym_buffer(1)
317 rotation_env%buf2_k_k_nosym => nonsym_buffer(2)
318 rotation_env%buf3_k_k_nosym => nonsym_buffer(3)
319 rotation_env%buf1_n_k => nk_buffer(1)
320 rotation_env%buf1_n_k_dp => nk_buffer(2)
321 rotation_env%rot_mat_u => rot_matrix(1)
322 rotation_env%rot_mat_u_im => rot_matrix(2)
323 rotation_env%rot_mat_x => rot_matrix(3)
324 rotation_env%rot_mat_x_im => rot_matrix(4)
325 rotation_env%rot_mat_dedu => rot_matrix(5)
326 rotation_env%rot_mat_dedu_im => rot_matrix(6)
327 rotation_env%rot_mat_gx => rot_matrix(7)
328 rotation_env%rot_mat_gx_im => rot_matrix(8)
329 rotation_env%rot_mat_evec_re => rot_matrix(9)
330 rotation_env%rot_mat_evec_im => rot_matrix(10)
331 rotation_env%rot_mat_evals => rotation_evals
332
333 c_re_p => c_re
334 c_im_p => c_im
335 s_re_p => s_re
336 s_im_p => s_im
337 CALL qs_ot_get_orbitals_ref_complex(c_re_p, c_im_p, s_re_p, s_im_p, rotation_env)
338 CALL get_complex_pair(c_re, c_im, c_dense)
339 error_c = maxval(abs(c_dense - c_reference))
340
341 hc_dense = 2.0_dp*matmul(hamiltonian, c_dense)
342 DO j = 1, k
343 hc_dense(:, j) = occupation(j)*hc_dense(:, j)
344 END DO
345 CALL put_complex_pair(hc_re, hc_im, hc_dense)
346 hc_re_p => hc_re
347 hc_im_p => hc_im
348 CALL qs_ot_get_derivative_ref_complex(hc_re_p, hc_im_p, rotation_env)
349 CALL get_complex_pair(g_re, g_im, gradient_x)
350 CALL get_complex_pair(rot_matrix(7), rot_matrix(8), gradient_a)
351
352 hstep = 1.0e-6_dp
353 fd_ref = (rotated_polar_energy(x_dense + hstep*direction_x, generator, &
354 hamiltonian, occupation) - &
355 rotated_polar_energy(x_dense - hstep*direction_x, generator, &
356 hamiltonian, occupation))/(2.0_dp*hstep)
357 predicted_ref = real(sum(conjg(gradient_x)*direction_x), kind=dp)
358 error_ref = abs(predicted_ref - fd_ref)
359
360 fd_rot = (rotated_polar_energy(x_dense, generator + hstep*direction_a, &
361 hamiltonian, occupation) - &
362 rotated_polar_energy(x_dense, generator - hstep*direction_a, &
363 hamiltonian, occupation))/(2.0_dp*hstep)
364 predicted_rot = 0.5_dp*(sum(real(gradient_a, kind=dp)*real(direction_a, kind=dp)) + &
365 sum(aimag(gradient_a)*aimag(direction_a)))
366 error_rot = abs(predicted_rot - fd_rot)
367
368 IF (error_c > 5.0e-11_dp .OR. error_ref > 5.0e-8_dp .OR. &
369 error_rot > 5.0e-8_dp) nfail = nfail + 1
370 IF (para_env%is_source()) THEN
371 WRITE (io_unit, '(A,1X,ES13.6)') 'complex finite REF rotation orbital error', error_c
372 WRITE (io_unit, '(A,3(1X,ES13.6))') 'complex finite REF slope fd/predicted/error', &
373 fd_ref, predicted_ref, error_ref
374 WRITE (io_unit, '(A,3(1X,ES13.6))') 'complex finite rotation slope fd/predicted/error', &
375 fd_rot, predicted_rot, error_rot
376 END IF
377
378 DO i = 1, SIZE(rot_matrix)
379 CALL dbcsr_release(rot_matrix(i))
380 END DO
381 DO i = 1, SIZE(nk_buffer)
382 CALL dbcsr_release(nk_buffer(i))
383 END DO
384 DO i = 1, SIZE(nonsym_buffer)
385 CALL dbcsr_release(nonsym_buffer(i))
386 END DO
387 DO i = 1, SIZE(sym_buffer)
388 CALL dbcsr_release(sym_buffer(i))
389 END DO
390 CALL dbcsr_release(f_im)
391 CALL dbcsr_release(f_re)
392 CALL dbcsr_release(s_im)
393 CALL dbcsr_release(s_re)
394 CALL dbcsr_release(g_im)
395 CALL dbcsr_release(g_re)
396 CALL dbcsr_release(hc_im)
397 CALL dbcsr_release(hc_re)
398 CALL dbcsr_release(c_im)
399 CALL dbcsr_release(c_re)
400 CALL dbcsr_release(sx_im)
401 CALL dbcsr_release(sx_re)
402 CALL dbcsr_release(x_im)
403 CALL dbcsr_release(x_re)
404 CALL cp_blacs_env_release(blacs_env)
405 CALL mp_para_env_release(local_para_env)
406
407 END SUBROUTINE test_complex_ref_rotation
408
409! **************************************************************************************************
410!> \brief Check union handling for distinct real/imaginary DBCSR block patterns.
411!> \param para_env parallel environment
412!> \param nfail accumulated number of failures
413! **************************************************************************************************
414 SUBROUTINE test_sparse_frechet_patterns(para_env, nfail)
415 TYPE(mp_para_env_type), POINTER :: para_env
416 INTEGER, INTENT(INOUT) :: nfail
417
418 COMPLEX(KIND=dp) :: expected, value
419 INTEGER, DIMENSION(:), POINTER :: col_dist, col_size, row_dist, row_size
420 INTEGER, DIMENSION(:, :), POINTER :: pgrid
421 LOGICAL :: found
422 REAL(KIND=dp) :: error
423 REAL(KIND=dp), DIMENSION(2) :: evals
424 REAL(KIND=dp), DIMENSION(:, :), POINTER :: block
425 TYPE(dbcsr_distribution_type) :: dist
426 TYPE(dbcsr_type) :: inner_im, inner_re, outer_im, outer_re
427
428 ALLOCATE (pgrid(0:0, 0:0), row_dist(2), col_dist(2), row_size(2), col_size(2))
429 pgrid(:, :) = 0
430 row_dist(:) = 0
431 col_dist(:) = 0
432 row_size(:) = 1
433 col_size(:) = 1
434 CALL dbcsr_distribution_new(dist, group=mp_comm_self%get_handle(), pgrid=pgrid, &
435 row_dist=row_dist, col_dist=col_dist)
436 CALL dbcsr_create(inner_re, 'frechet_sparse_re', dist, dbcsr_type_no_symmetry, row_size, col_size)
437 CALL dbcsr_reserve_blocks(inner_re, rows=[1], cols=[2])
438 CALL dbcsr_finalize(inner_re)
439 CALL dbcsr_create(inner_im, 'frechet_sparse_im', dist, dbcsr_type_no_symmetry, row_size, col_size)
440 CALL dbcsr_reserve_blocks(inner_im, rows=[2], cols=[1])
441 CALL dbcsr_finalize(inner_im)
442 CALL dbcsr_put_block(inner_re, 1, 2, reshape([0.37_dp], [1, 1]))
443 CALL dbcsr_put_block(inner_im, 2, 1, reshape([-0.21_dp], [1, 1]))
444
445 evals = [-0.4_dp, 0.9_dp]
446 CALL qs_ot_apply_complex_frechet_dbcsr(evals, inner_re, inner_im, outer_re, outer_im)
447 error = 0.0_dp
448 expected = qs_ot_complex_exp_frechet_kernel(evals(1), evals(2))*cmplx(0.37_dp, 0.0_dp, kind=dp)
449 value = cmplx(0.0_dp, 0.0_dp, kind=dp)
450 CALL dbcsr_get_block_p(outer_re, 1, 2, block, found)
451 IF (found) value = value + cmplx(block(1, 1), 0.0_dp, kind=dp)
452 CALL dbcsr_get_block_p(outer_im, 1, 2, block, found)
453 IF (found) value = value + cmplx(0.0_dp, block(1, 1), kind=dp)
454 error = max(error, abs(value - expected))
455
456 expected = qs_ot_complex_exp_frechet_kernel(evals(2), evals(1))*cmplx(0.0_dp, -0.21_dp, kind=dp)
457 value = cmplx(0.0_dp, 0.0_dp, kind=dp)
458 CALL dbcsr_get_block_p(outer_re, 2, 1, block, found)
459 IF (found) value = value + cmplx(block(1, 1), 0.0_dp, kind=dp)
460 CALL dbcsr_get_block_p(outer_im, 2, 1, block, found)
461 IF (found) value = value + cmplx(0.0_dp, block(1, 1), kind=dp)
462 error = max(error, abs(value - expected))
463 CALL check_error("sparse_frechet_patterns", error, 5.0e-14_dp, nfail)
464 IF (para_env%is_source()) THEN
465 WRITE (io_unit, '(A,1X,ES13.6)') 'sparse complex Frechet block-union error', error
466 END IF
467
468 CALL dbcsr_release(outer_im)
469 CALL dbcsr_release(outer_re)
470 CALL dbcsr_release(inner_im)
471 CALL dbcsr_release(inner_re)
472 CALL dbcsr_distribution_release(dist)
473 DEALLOCATE (pgrid, row_dist, col_dist, row_size, col_size)
474
475 END SUBROUTINE test_sparse_frechet_patterns
476
477! **************************************************************************************************
478!> \brief Check the finite complex STRICT chart and its analytic derivative.
479!> \param dist single-block matrix distribution
480!> \param row_n AO block sizes
481!> \param row_k occupied-space row block sizes
482!> \param col_k occupied-space column block sizes
483!> \param para_env parallel environment
484!> \param nfail accumulated number of failures
485! **************************************************************************************************
486 SUBROUTINE test_complex_strict_derivative(dist, row_n, row_k, col_k, para_env, nfail)
487 TYPE(dbcsr_distribution_type), INTENT(IN) :: dist
488 INTEGER, DIMENSION(:), INTENT(INOUT), POINTER :: row_n, row_k, col_k
489 TYPE(mp_para_env_type), POINTER :: para_env
490 INTEGER, INTENT(INOUT) :: nfail
491
492 CHARACTER(LEN=32) :: name
493 COMPLEX(KIND=dp), DIMENSION(k, k) :: direction_a, generator, gradient_a, &
494 inverse_sqrt, rotation
495 COMPLEX(KIND=dp), DIMENSION(n, k) :: c0, direction, gradient, hc, orbitals, &
496 orbitals_base, sc0, sx, x
497 COMPLEX(KIND=dp), DIMENSION(n, n) :: hamiltonian, overlap
498 INTEGER :: i
499 REAL(KIND=dp) :: error, error_rot, fd_slope, &
500 fd_slope_rot, hstep, slope, slope_rot
501 REAL(KIND=dp), DIMENSION(k) :: occupation
502 REAL(KIND=dp), DIMENSION(k), TARGET :: rotation_evals
503 TYPE(cp_blacs_env_type), POINTER :: blacs_env
504 TYPE(dbcsr_type), DIMENSION(10), TARGET :: rot_matrix
505 TYPE(dbcsr_type), DIMENSION(15), TARGET :: nk
506 TYPE(dbcsr_type), DIMENSION(25), TARGET :: kk
507 TYPE(dbcsr_type), POINTER :: hc_im_p, hc_re_p, sx_im_p, sx_re_p, &
508 x_im_p, x_re_p
509 TYPE(mp_para_env_type), POINTER :: strict_para_env
510 TYPE(qs_ot_type) :: strict_env
511
512 NULLIFY (blacs_env, hc_im_p, hc_re_p, strict_para_env, x_im_p, x_re_p, sx_im_p, sx_re_p)
513 ALLOCATE (strict_para_env)
514 CALL strict_para_env%from_dup(mp_comm_self)
515 CALL cp_blacs_env_create(blacs_env=blacs_env, para_env=strict_para_env)
516 DO i = 1, SIZE(nk)
517 WRITE (name, '(A,I0)') 'complex_strict_nk_', i
518 CALL create_one_block(nk(i), trim(name), dist, row_n, col_k)
519 END DO
520 DO i = 1, SIZE(kk)
521 WRITE (name, '(A,I0)') 'complex_strict_kk_', i
522 CALL create_one_block(kk(i), trim(name), dist, row_k, col_k)
523 END DO
524 DO i = 1, SIZE(rot_matrix)
525 WRITE (name, '(A,I0)') 'complex_strict_rot_', i
526 CALL create_one_block(rot_matrix(i), trim(name), dist, row_k, col_k)
527 END DO
528
529 overlap(:, :) = cmplx(0.0_dp, 0.0_dp, kind=dp)
530 overlap(1, 1) = cmplx(1.0_dp, 0.0_dp, kind=dp)
531 overlap(2, 2) = cmplx(1.2_dp, 0.0_dp, kind=dp)
532 overlap(3, 3) = cmplx(1.4_dp, 0.0_dp, kind=dp)
533 overlap(1, 2) = cmplx(0.06_dp, 0.03_dp, kind=dp)
534 overlap(2, 1) = conjg(overlap(1, 2))
535 overlap(1, 3) = cmplx(-0.04_dp, 0.02_dp, kind=dp)
536 overlap(3, 1) = conjg(overlap(1, 3))
537 overlap(2, 3) = cmplx(0.05_dp, -0.01_dp, kind=dp)
538 overlap(3, 2) = conjg(overlap(2, 3))
539
540 c0(:, 1) = [cmplx(1.0_dp, 0.10_dp, kind=dp), cmplx(0.12_dp, -0.08_dp, kind=dp), &
541 cmplx(-0.07_dp, 0.04_dp, kind=dp)]
542 c0(:, 2) = [cmplx(0.09_dp, 0.05_dp, kind=dp), cmplx(0.85_dp, -0.12_dp, kind=dp), &
543 cmplx(0.11_dp, 0.06_dp, kind=dp)]
544 CALL dense_inverse_sqrt(matmul(conjg(transpose(c0)), matmul(overlap, c0)), inverse_sqrt)
545 c0 = matmul(c0, inverse_sqrt)
546
547 x(:, :) = cmplx(0.0_dp, 0.0_dp, kind=dp)
548 x(3, 1) = cmplx(0.23_dp, 0.17_dp, kind=dp)
549 x(3, 2) = cmplx(-0.14_dp, 0.21_dp, kind=dp)
550 x = x - matmul(c0, matmul(conjg(transpose(c0)), matmul(overlap, x)))
551 direction(:, :) = cmplx(0.0_dp, 0.0_dp, kind=dp)
552 direction(3, 1) = cmplx(0.11_dp, -0.09_dp, kind=dp)
553 direction(3, 2) = cmplx(0.07_dp, 0.13_dp, kind=dp)
554 direction = direction - matmul(c0, matmul(conjg(transpose(c0)), matmul(overlap, direction)))
555
556 hamiltonian(:, :) = cmplx(0.0_dp, 0.0_dp, kind=dp)
557 hamiltonian(1, 1) = cmplx(-0.70_dp, 0.0_dp, kind=dp)
558 hamiltonian(2, 2) = cmplx(-0.20_dp, 0.0_dp, kind=dp)
559 hamiltonian(3, 3) = cmplx(0.90_dp, 0.0_dp, kind=dp)
560 hamiltonian(1, 2) = cmplx(0.08_dp, 0.04_dp, kind=dp)
561 hamiltonian(2, 1) = conjg(hamiltonian(1, 2))
562 hamiltonian(1, 3) = cmplx(-0.12_dp, 0.07_dp, kind=dp)
563 hamiltonian(3, 1) = conjg(hamiltonian(1, 3))
564 hamiltonian(2, 3) = cmplx(0.05_dp, -0.09_dp, kind=dp)
565 hamiltonian(3, 2) = conjg(hamiltonian(2, 3))
566
567 generator(:, :) = cmplx(0.0_dp, 0.0_dp, kind=dp)
568 generator(1, 1) = cmplx(0.0_dp, 0.16_dp, kind=dp)
569 generator(2, 2) = cmplx(0.0_dp, -0.09_dp, kind=dp)
570 generator(1, 2) = cmplx(0.18_dp, -0.07_dp, kind=dp)
571 generator(2, 1) = -conjg(generator(1, 2))
572 direction_a(:, :) = cmplx(0.0_dp, 0.0_dp, kind=dp)
573 direction_a(1, 1) = cmplx(0.0_dp, -0.05_dp, kind=dp)
574 direction_a(2, 2) = cmplx(0.0_dp, 0.04_dp, kind=dp)
575 direction_a(1, 2) = cmplx(-0.08_dp, 0.06_dp, kind=dp)
576 direction_a(2, 1) = -conjg(direction_a(1, 2))
577 occupation = [1.0_dp, 0.34_dp]
578
579 CALL strict_orbitals(c0, x, overlap, orbitals_base)
580 rotation = dense_antihermitian_exp(generator)
581 orbitals = matmul(orbitals_base, rotation)
582 hc = 2.0_dp*matmul(hamiltonian, orbitals)
583 DO i = 1, k
584 hc(:, i) = occupation(i)*hc(:, i)
585 END DO
586 sx = matmul(overlap, x)
587 sc0 = matmul(overlap, c0)
588 CALL put_complex_pair(nk(1), nk(2), x)
589 CALL put_complex_pair(nk(3), nk(4), sx)
590 CALL put_complex_pair(nk(5), nk(6), hc)
591 CALL put_complex_pair(nk(9), nk(10), c0)
592 CALL put_complex_pair(nk(11), nk(12), sc0)
593 CALL put_complex_pair(rot_matrix(3), rot_matrix(4), generator)
594
595 strict_env%has_complex_kpoint_state = .true.
596 strict_env%settings%ot_algorithm = 'TOD'
597 strict_env%settings%do_rotation = .true.
598 strict_env%para_env => strict_para_env
599 strict_env%blacs_env => blacs_env
600 strict_env%matrix_x => nk(1)
601 strict_env%matrix_x_im => nk(2)
602 strict_env%matrix_sx => nk(3)
603 strict_env%matrix_sx_im => nk(4)
604 strict_env%matrix_gx => nk(7)
605 strict_env%matrix_gx_im => nk(8)
606 strict_env%matrix_c0 => nk(9)
607 strict_env%matrix_c0_im => nk(10)
608 strict_env%matrix_sc0 => nk(11)
609 strict_env%matrix_sc0_im => nk(12)
610 strict_env%matrix_buf_nk => nk(13)
611 strict_env%matrix_buf_nk_im => nk(14)
612 strict_env%matrix_tmp_nk => nk(15)
613 strict_env%matrix_p => kk(1)
614 strict_env%matrix_p_im => kk(2)
615 strict_env%matrix_r => kk(3)
616 strict_env%matrix_r_im => kk(4)
617 strict_env%matrix_sinp => kk(5)
618 strict_env%matrix_sinp_im => kk(6)
619 strict_env%matrix_cosp => kk(7)
620 strict_env%matrix_cosp_im => kk(8)
621 strict_env%matrix_sinp_b => kk(9)
622 strict_env%matrix_cosp_b => kk(10)
623 strict_env%matrix_buf1 => kk(11)
624 strict_env%matrix_buf1_im => kk(12)
625 strict_env%matrix_buf2 => kk(13)
626 strict_env%matrix_buf2_im => kk(14)
627 strict_env%matrix_buf3 => kk(15)
628 strict_env%matrix_buf3_im => kk(16)
629 strict_env%matrix_buf4 => kk(17)
630 strict_env%matrix_buf4_im => kk(18)
631 strict_env%matrix_os => kk(19)
632 strict_env%matrix_os_im => kk(20)
633 strict_env%matrix_buf1_ortho => kk(21)
634 strict_env%matrix_buf1_ortho_im => kk(22)
635 strict_env%matrix_buf2_ortho => kk(23)
636 strict_env%matrix_buf2_ortho_im => kk(24)
637 strict_env%matrix_tmp_ortho => kk(25)
638 strict_env%rot_mat_u => rot_matrix(1)
639 strict_env%rot_mat_u_im => rot_matrix(2)
640 strict_env%rot_mat_x => rot_matrix(3)
641 strict_env%rot_mat_x_im => rot_matrix(4)
642 strict_env%rot_mat_dedu => rot_matrix(5)
643 strict_env%rot_mat_dedu_im => rot_matrix(6)
644 strict_env%rot_mat_gx => rot_matrix(7)
645 strict_env%rot_mat_gx_im => rot_matrix(8)
646 strict_env%rot_mat_evec_re => rot_matrix(9)
647 strict_env%rot_mat_evec_im => rot_matrix(10)
648 strict_env%rot_mat_evals => rotation_evals
649 ALLOCATE (strict_env%evals(k), strict_env%dum(k))
650
651 x_re_p => nk(1)
652 x_im_p => nk(2)
653 sx_re_p => nk(3)
654 sx_im_p => nk(4)
655 hc_re_p => nk(5)
656 hc_im_p => nk(6)
657 CALL qs_ot_get_p_complex(x_re_p, x_im_p, sx_re_p, sx_im_p, strict_env)
658 CALL qs_ot_get_derivative_complex(hc_re_p, hc_im_p, strict_env)
659
660 CALL get_complex_pair(nk(7), nk(8), gradient)
661 CALL get_complex_pair(rot_matrix(7), rot_matrix(8), gradient_a)
662
663 hstep = 1.0e-6_dp
664 fd_slope = (strict_rotated_energy(c0, x + hstep*direction, generator, overlap, &
665 hamiltonian, occupation) - &
666 strict_rotated_energy(c0, x - hstep*direction, generator, overlap, &
667 hamiltonian, occupation))/(2.0_dp*hstep)
668 slope = real(sum(conjg(gradient)*direction), kind=dp)
669 error = abs(slope - fd_slope)
670 fd_slope_rot = (strict_rotated_energy(c0, x, generator + hstep*direction_a, overlap, &
671 hamiltonian, occupation) - &
672 strict_rotated_energy(c0, x, generator - hstep*direction_a, overlap, &
673 hamiltonian, occupation))/(2.0_dp*hstep)
674 slope_rot = 0.5_dp*(sum(real(gradient_a, kind=dp)*real(direction_a, kind=dp)) + &
675 sum(aimag(gradient_a)*aimag(direction_a)))
676 error_rot = abs(slope_rot - fd_slope_rot)
677 IF (error > 5.0e-8_dp .OR. error_rot > 5.0e-8_dp) nfail = nfail + 1
678 IF (para_env%is_source()) THEN
679 WRITE (io_unit, '(A,3(1X,ES13.6))') 'complex finite-STRICT derivative: fd/predicted/error', &
680 fd_slope, slope, error
681 WRITE (io_unit, '(A,3(1X,ES13.6))') 'complex finite-STRICT rotation: fd/predicted/error', &
682 fd_slope_rot, slope_rot, error_rot
683 END IF
684
685 DEALLOCATE (strict_env%evals, strict_env%dum)
686 DO i = SIZE(rot_matrix), 1, -1
687 CALL dbcsr_release(rot_matrix(i))
688 END DO
689 DO i = SIZE(kk), 1, -1
690 CALL dbcsr_release(kk(i))
691 END DO
692 DO i = SIZE(nk), 1, -1
693 CALL dbcsr_release(nk(i))
694 END DO
695 CALL cp_blacs_env_release(blacs_env)
696 CALL mp_para_env_release(strict_para_env)
697
698 END SUBROUTINE test_complex_strict_derivative
699
700! **************************************************************************************************
701!> \brief Create and reserve a one-block DBCSR matrix.
702!> \param matrix ...
703!> \param name ...
704!> \param dist ...
705!> \param row_size ...
706!> \param col_size ...
707! **************************************************************************************************
708 SUBROUTINE create_one_block(matrix, name, dist, row_size, col_size)
709 TYPE(dbcsr_type), INTENT(INOUT) :: matrix
710 CHARACTER(LEN=*), INTENT(IN) :: name
711 TYPE(dbcsr_distribution_type), INTENT(IN) :: dist
712 INTEGER, DIMENSION(:), INTENT(INOUT), POINTER :: row_size, col_size
713
714 CALL dbcsr_create(matrix, name, dist, dbcsr_type_no_symmetry, row_size, col_size)
715 CALL dbcsr_reserve_blocks(matrix, rows=[1], cols=[1])
716 CALL dbcsr_finalize(matrix)
717
718 END SUBROUTINE create_one_block
719
720! **************************************************************************************************
721!> \brief Store a dense complex matrix in a real/imaginary DBCSR pair.
722!> \param matrix_re ...
723!> \param matrix_im ...
724!> \param matrix ...
725! **************************************************************************************************
726 SUBROUTINE put_complex_pair(matrix_re, matrix_im, matrix)
727 TYPE(dbcsr_type), INTENT(INOUT) :: matrix_re, matrix_im
728 COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(IN) :: matrix
729
730 CALL dbcsr_put_block(matrix_re, 1, 1, real(matrix, kind=dp))
731 CALL dbcsr_put_block(matrix_im, 1, 1, aimag(matrix))
732
733 END SUBROUTINE put_complex_pair
734
735! **************************************************************************************************
736!> \brief Retrieve a dense complex matrix from a one-block DBCSR pair.
737!> \param matrix_re ...
738!> \param matrix_im ...
739!> \param matrix ...
740! **************************************************************************************************
741 SUBROUTINE get_complex_pair(matrix_re, matrix_im, matrix)
742 TYPE(dbcsr_type), INTENT(INOUT) :: matrix_re, matrix_im
743 COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: matrix
744
745 LOGICAL :: found
746 REAL(KIND=dp), DIMENSION(:, :), POINTER :: block
747
748 matrix(:, :) = cmplx(0.0_dp, 0.0_dp, kind=dp)
749 CALL dbcsr_get_block_p(matrix_re, 1, 1, block, found)
750 IF (found) matrix = matrix + cmplx(block, 0.0_dp, kind=dp)
751 CALL dbcsr_get_block_p(matrix_im, 1, 1, block, found)
752 IF (found) matrix = matrix + cmplx(0.0_dp, block, kind=dp)
753
754 END SUBROUTINE get_complex_pair
755
756! **************************************************************************************************
757!> \brief Dense fixed-occupation energy after finite REF orthogonalization and rotation.
758!> \param x REF coordinate
759!> \param generator anti-Hermitian rotation generator
760!> \param hamiltonian Hermitian one-particle Hamiltonian
761!> \param occupation orbital occupations
762!> \return energy
763! **************************************************************************************************
764 FUNCTION rotated_polar_energy(x, generator, hamiltonian, occupation) RESULT(energy)
765 COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(IN) :: x, generator, hamiltonian
766 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: occupation
767 REAL(KIND=dp) :: energy
768
769 COMPLEX(KIND=dp), &
770 DIMENSION(SIZE(x, 1), SIZE(x, 2)) :: orbitals
771 COMPLEX(KIND=dp), &
772 DIMENSION(SIZE(x, 2), SIZE(x, 2)) :: inverse_sqrt, rotation
773 INTEGER :: i
774
775 CALL dense_inverse_sqrt(matmul(conjg(transpose(x)), x), inverse_sqrt)
776 rotation = dense_antihermitian_exp(generator)
777 orbitals = matmul(matmul(x, inverse_sqrt), rotation)
778 energy = 0.0_dp
779 DO i = 1, SIZE(occupation)
780 energy = energy + occupation(i)* &
781 REAL(DOT_PRODUCT(orbitals(:, i), MATMUL(hamiltonian, orbitals(:, i))), KIND=dp)
782 END DO
783
784 END FUNCTION rotated_polar_energy
785
786! **************************************************************************************************
787!> \brief Dense reference exponential for an anti-Hermitian matrix.
788!> \param generator anti-Hermitian generator
789!> \return its unitary exponential
790! **************************************************************************************************
791 FUNCTION dense_antihermitian_exp(generator) RESULT(rotation)
792 COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(IN) :: generator
793 COMPLEX(KIND=dp), DIMENSION(SIZE(generator, 1), & SIZE(generator, 2)) :: rotation
794
795 COMPLEX(KIND=dp), DIMENSION(SIZE(generator, 1), & SIZE(generator, 2)) :: hermitian, vectors
796 INTEGER :: i
797 REAL(KIND=dp), DIMENSION(SIZE(generator, 1)) :: eigenvalues
798
799 hermitian = cmplx(0.0_dp, 1.0_dp, kind=dp)*generator
800 CALL diag_complex(hermitian, vectors, eigenvalues)
801 rotation(:, :) = cmplx(0.0_dp, 0.0_dp, kind=dp)
802 DO i = 1, SIZE(eigenvalues)
803 rotation = rotation + exp(cmplx(0.0_dp, -eigenvalues(i), kind=dp))* &
804 matmul(reshape(vectors(:, i), [SIZE(generator, 1), 1]), &
805 reshape(conjg(vectors(:, i)), [1, SIZE(generator, 1)]))
806 END DO
807
808 END FUNCTION dense_antihermitian_exp
809
810! **************************************************************************************************
811!> \brief Hermitian inverse square root for the independent dense reference.
812!> \param matrix ...
813!> \param inverse ...
814! **************************************************************************************************
815 SUBROUTINE dense_inverse_sqrt(matrix, inverse)
816 COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(IN) :: matrix
817 COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: inverse
818
819 COMPLEX(KIND=dp), &
820 DIMENSION(SIZE(matrix, 1), SIZE(matrix, 2)) :: vectors
821 INTEGER :: i
822 REAL(KIND=dp), DIMENSION(SIZE(matrix, 1)) :: eigenvalues
823
824 CALL diag_complex(matrix, vectors, eigenvalues)
825 inverse(:, :) = cmplx(0.0_dp, 0.0_dp, kind=dp)
826 DO i = 1, SIZE(eigenvalues)
827 inverse = inverse + matmul(reshape(vectors(:, i), [SIZE(matrix, 1), 1]), &
828 reshape(conjg(vectors(:, i)), [1, SIZE(matrix, 1)]))/ &
829 sqrt(eigenvalues(i))
830 END DO
831
832 END SUBROUTINE dense_inverse_sqrt
833
834! **************************************************************************************************
835!> \brief Independent dense finite STRICT transformation for a Hermitian overlap.
836!> \param c0 reference orbitals
837!> \param coordinate tangent coordinate
838!> \param overlap Hermitian overlap matrix
839!> \param orbitals transformed orbitals
840! **************************************************************************************************
841 SUBROUTINE strict_orbitals(c0, coordinate, overlap, orbitals)
842 COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(IN) :: c0, coordinate, overlap
843 COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: orbitals
844
845 COMPLEX(KIND=dp), DIMENSION(SIZE(coordinate, 2), & SIZE(coordinate, 2)) :: cosp, sinp, vectors
846 INTEGER :: i
847 REAL(KIND=dp) :: root
848 REAL(KIND=dp), DIMENSION(SIZE(coordinate, 2)) :: eigenvalues
849
850 CALL diag_complex(matmul(conjg(transpose(coordinate)), matmul(overlap, coordinate)), &
851 vectors, eigenvalues)
852 cosp(:, :) = cmplx(0.0_dp, 0.0_dp, kind=dp)
853 sinp(:, :) = cmplx(0.0_dp, 0.0_dp, kind=dp)
854 DO i = 1, SIZE(eigenvalues)
855 root = sqrt(max(0.0_dp, eigenvalues(i)))
856 cosp = cosp + cos(root)*matmul(reshape(vectors(:, i), [SIZE(vectors, 1), 1]), &
857 reshape(conjg(vectors(:, i)), [1, SIZE(vectors, 1)]))
858 IF (root < 1.0e-12_dp) THEN
859 sinp = sinp + matmul(reshape(vectors(:, i), [SIZE(vectors, 1), 1]), &
860 reshape(conjg(vectors(:, i)), [1, SIZE(vectors, 1)]))
861 ELSE
862 sinp = sinp + sin(root)/root* &
863 matmul(reshape(vectors(:, i), [SIZE(vectors, 1), 1]), &
864 reshape(conjg(vectors(:, i)), [1, SIZE(vectors, 1)]))
865 END IF
866 END DO
867 orbitals = matmul(c0, cosp) + matmul(coordinate, sinp)
868
869 END SUBROUTINE strict_orbitals
870
871! **************************************************************************************************
872!> \brief Energy of finite STRICT orbitals followed by a complex occupied-space rotation.
873!> \param c0 reference orbitals
874!> \param coordinate tangent coordinate
875!> \param generator anti-Hermitian rotation generator
876!> \param overlap Hermitian overlap matrix
877!> \param hamiltonian Hermitian Hamiltonian
878!> \param occupation orbital occupations
879!> \return energy
880! **************************************************************************************************
881 FUNCTION strict_rotated_energy(c0, coordinate, generator, overlap, hamiltonian, &
882 occupation) RESULT(energy)
883 COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(IN) :: c0, coordinate, generator, overlap, &
884 hamiltonian
885 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: occupation
886 REAL(KIND=dp) :: energy
887
888 COMPLEX(KIND=dp), DIMENSION(SIZE(coordinate, 1), & SIZE(coordinate, 2)) :: orbitals
889 INTEGER :: i
890
891 CALL strict_orbitals(c0, coordinate, overlap, orbitals)
892 orbitals = matmul(orbitals, dense_antihermitian_exp(generator))
893 energy = 0.0_dp
894 DO i = 1, SIZE(occupation)
895 energy = energy + occupation(i)* &
896 REAL(DOT_PRODUCT(orbitals(:, i), MATMUL(hamiltonian, orbitals(:, i))), KIND=dp)
897 END DO
898
899 END FUNCTION strict_rotated_energy
900
901! **************************************************************************************************
902!> \brief Report a failed numerical check, including non-finite results.
903!> \param label check name
904!> \param error measured error
905!> \param tolerance maximum error
906!> \param nfail accumulated failures
907! **************************************************************************************************
908 SUBROUTINE check_error(label, error, tolerance, nfail)
909 CHARACTER(LEN=*), INTENT(IN) :: label
910 REAL(KIND=dp), INTENT(IN) :: error, tolerance
911 INTEGER, INTENT(INOUT) :: nfail
912
913 IF (ieee_is_finite(error)) THEN
914 IF (error <= tolerance) RETURN
915 END IF
916 nfail = nfail + 1
917 IF (io_unit >= 0) WRITE (io_unit, '(A,2(1X,ES13.6))') trim(label)//': error/tolerance', error, tolerance
918 END SUBROUTINE check_error
919
920END PROGRAM qs_ot_complex_ref_unittest
921
cp_dbcsr_api::dbcsr_create
Definition cp_dbcsr_api.F:194
bibliography
collects all references to literature in CP2K as new algorithms / method are included from literature...
Definition bibliography.F:21
bibliography::add_all_references
subroutine, public add_all_references()
adds references that can later be cited / printed using the key
Definition bibliography.F:126
cp_blacs_env
methods related to the blacs parallel environment
Definition cp_blacs_env.F:15
cp_blacs_env::cp_blacs_env_release
subroutine, public cp_blacs_env_release(blacs_env)
releases the given blacs_env
Definition cp_blacs_env.F:274
cp_blacs_env::cp_blacs_env_create
subroutine, public cp_blacs_env_create(blacs_env, para_env, blacs_grid_layout, blacs_repeatable, row_major, grid_2d)
allocates and initializes a type that represent a blacs context
Definition cp_blacs_env.F:123
cp_dbcsr_api
Definition cp_dbcsr_api.F:8
cp_dbcsr_api::dbcsr_distribution_release
subroutine, public dbcsr_distribution_release(dist)
...
Definition cp_dbcsr_api.F:643
cp_dbcsr_api::dbcsr_distribution_new
subroutine, public dbcsr_distribution_new(dist, template, group, pgrid, row_dist, col_dist, reuse_arrays)
...
Definition cp_dbcsr_api.F:616
cp_dbcsr_api::dbcsr_get_block_p
subroutine, public dbcsr_get_block_p(matrix, row, col, block, found, row_size, col_size)
...
Definition cp_dbcsr_api.F:693
cp_dbcsr_api::dbcsr_reserve_blocks
subroutine, public dbcsr_reserve_blocks(matrix, rows, cols)
...
Definition cp_dbcsr_api.F:1172
cp_dbcsr_api::dbcsr_finalize
subroutine, public dbcsr_finalize(matrix)
...
Definition cp_dbcsr_api.F:673
cp_dbcsr_api::dbcsr_release
subroutine, public dbcsr_release(matrix)
...
Definition cp_dbcsr_api.F:1142
cp_dbcsr_api::dbcsr_put_block
subroutine, public dbcsr_put_block(matrix, row, col, block, summation)
...
Definition cp_dbcsr_api.F:1125
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition cp_log_handling.F:41
cp_log_handling::cp_rm_default_logger
subroutine, public cp_rm_default_logger()
the cousin of cp_add_default_logger, decrements the stack, so that the default logger is what it has ...
Definition cp_log_handling.F:213
cp_log_handling::cp_logger_release
subroutine, public cp_logger_release(logger)
releases this logger
Definition cp_log_handling.F:429
cp_log_handling::cp_logger_create
subroutine, public cp_logger_create(logger, para_env, print_level, default_global_unit_nr, default_local_unit_nr, global_filename, local_filename, close_global_unit_on_dealloc, iter_info, close_local_unit_on_dealloc, suffix, template_logger)
initializes a logger
Definition cp_log_handling.F:285
cp_log_handling::cp_add_default_logger
subroutine, public cp_add_default_logger(logger)
adds a default logger. MUST be called before logging occours
Definition cp_log_handling.F:188
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
machine
Machine interface based on Fortran 2003 and POSIX.
Definition machine.F:17
machine::default_output_unit
integer, parameter, public default_output_unit
Definition machine.F:48
mathlib
Collection of simple mathematical functions and subroutines.
Definition mathlib.F:15
mathlib::diag_complex
subroutine, public diag_complex(matrix, eigenvectors, eigenvalues)
Diagonalizes a local complex Hermitian matrix using LAPACK. Based on cp_cfm_heevd.
Definition mathlib.F:1907
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:27
message_passing::mp_world_init
subroutine, public mp_world_init(mp_comm)
initializes the system default communicator
Definition message_passing.F:1389
message_passing::mp_para_env_release
subroutine, public mp_para_env_release(para_env)
releases the para object (to be called when you don't want anymore the shared copy of this object)
Definition message_passing.F:2163
message_passing::mp_comm_self
type(mp_comm_type), parameter, public mp_comm_self
Definition message_passing.F:817
message_passing::mp_world_finalize
subroutine, public mp_world_finalize()
Finalize the system default communicator and MPI when CP2K owns MPI.
Definition message_passing.F:1481
qs_ot_types
orbital transformations
Definition qs_ot_types.F:15
qs_ot
orbital transformations
Definition qs_ot.F:15
qs_ot::qs_ot_get_derivative_ref_complex
subroutine, public qs_ot_get_derivative_ref_complex(matrix_hc, matrix_hc_im, qs_ot_env, matrix_hc_rotation, matrix_hc_rotation_im)
complex k-point REF derivative dE/dX from H(k)C(k), S(k)C(k), and C(k)
Definition qs_ot.F:1997
qs_ot::qs_ot_get_orbitals_ref_complex
subroutine, public qs_ot_get_orbitals_ref_complex(matrix_c, matrix_c_im, matrix_s, matrix_s_im, qs_ot_env, qs_ot_env1)
update complex REF k-point orbitals and their S(k)C(k) images
Definition qs_ot.F:1689
qs_ot::qs_ot_get_derivative_complex
subroutine, public qs_ot_get_derivative_complex(matrix_hc, matrix_hc_im, qs_ot_env, matrix_hc_rotation, matrix_hc_rotation_im)
finite complex STRICT derivative, projected onto C0^H*S*X=0
Definition qs_ot.F:3254
qs_ot::qs_ot_complex_exp_frechet_kernel
pure complex(kind=dp) function, public qs_ot_complex_exp_frechet_kernel(e1, e2)
Frechet divided-difference kernel for exp(-i*evals).
Definition qs_ot.F:910
qs_ot::qs_ot_apply_complex_frechet_dbcsr
subroutine, public qs_ot_apply_complex_frechet_dbcsr(evals, inner_deriv_re, inner_deriv_im, outer_deriv_re, outer_deriv_im, adjoint)
apply the complex exponential Frechet kernel to sparse DBCSR Re/Im matrices
Definition qs_ot.F:943
qs_ot::qs_ot_density_tangent
subroutine, public qs_ot_density_tangent(rotation_generator, occupation, kpoint_weight, rotation_step, weighted_occupation_step, density_tangent, difference_step)
finite-chart density tangent for coupled complex rotations and fixed-N occupations
Definition qs_ot.F:461
qs_ot::qs_ot_get_p_complex
subroutine, public qs_ot_get_p_complex(matrix_x, matrix_x_im, matrix_sx, matrix_sx_im, qs_ot_env)
compute P=X^H*S*X and the STRICT matrix functions for a complex K-point channel
Definition qs_ot.F:2595
reference_manager
provides a uniform framework to add references to CP2K cite and output these
Definition reference_manager.F:25
reference_manager::remove_all_references
subroutine, public remove_all_references()
deallocate the bibliography
Definition reference_manager.F:227
qs_ot_complex_ref_unittest
program qs_ot_complex_ref_unittest
Focused regressions for complex k-point OT derivatives.
Definition qs_ot_complex_ref_unittest.F:11
test_density_tangent
subroutine test_density_tangent(mynode, nfail)
Check the coupled finite-chart density tangent against a central finite difference.
Definition qs_ot_complex_ref_unittest.F:106
cp_blacs_env::cp_blacs_env_type
represent a blacs multidimensional parallel environment (for the mpi corrispective see cp_paratypes/m...
Definition cp_blacs_env.F:53
cp_dbcsr_api::dbcsr_distribution_type
Definition cp_dbcsr_api.F:182
cp_dbcsr_api::dbcsr_type
Definition cp_dbcsr_api.F:176
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
message_passing::mp_comm_type
Definition message_passing.F:179
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:768
qs_ot_types::qs_ot_type
Definition qs_ot_types.F:140