(git:5e7fe52)
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 Regressions for complex k-point OT derivatives and preconditioner guards.
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_cfm_types, ONLY: cp_cfm_create,&
17 cp_cfm_get_element,&
18 cp_cfm_release,&
19 cp_cfm_set_all,&
20 cp_cfm_set_element,&
21 cp_cfm_type
22 USE cp_dbcsr_api, ONLY: &
23 dbcsr_create, dbcsr_distribution_new, dbcsr_distribution_release, dbcsr_distribution_type, &
24 dbcsr_finalize, dbcsr_finalize_lib, dbcsr_get_block_p, dbcsr_init_lib, dbcsr_put_block, &
25 dbcsr_release, dbcsr_reserve_blocks, dbcsr_type, dbcsr_type_no_symmetry
26 USE cp_fm_basic_linalg, ONLY: cp_complex_fm_gemm
27 USE cp_fm_struct, ONLY: cp_fm_struct_create,&
28 cp_fm_struct_release,&
29 cp_fm_struct_type
30 USE cp_fm_types, ONLY: cp_fm_create,&
31 cp_fm_get_submatrix,&
32 cp_fm_release,&
33 cp_fm_set_submatrix,&
34 cp_fm_type
35 USE cp_log_handling, ONLY: cp_add_default_logger,&
36 cp_logger_create,&
37 cp_logger_release,&
38 cp_logger_type,&
39 cp_rm_default_logger
40 USE input_constants, ONLY: &
41 ot_precond_fermi_low_rank, ot_precond_full_all, ot_precond_full_all_covariant, &
42 ot_precond_full_kinetic, ot_precond_full_single, ot_precond_full_single_inverse, &
43 ot_precond_none, ot_precond_s_inverse, ot_precond_solver_chebyshev, &
44 ot_precond_solver_default, ot_precond_solver_direct, ot_precond_solver_inv_chol
45 USE kinds, ONLY: dp
46 USE kpoint_methods, ONLY: kpoint_ot_energy_weighted_density
47 USE machine, ONLY: default_output_unit
48 USE mathlib, ONLY: diag_complex
49 USE message_passing, ONLY: mp_comm_self,&
50 mp_comm_type,&
51 mp_para_env_release,&
52 mp_para_env_type,&
53 mp_world_finalize,&
54 mp_world_init
55 USE preconditioner, ONLY: apply_preconditioner
56 USE preconditioner_makes, ONLY: make_complex_fermi_low_rank,&
57 make_complex_full_all,&
58 make_complex_full_all_covariant,&
59 make_complex_full_kinetic,&
60 make_complex_full_s_inverse,&
61 make_complex_full_single,&
62 make_complex_full_single_inverse
63 USE preconditioner_types, ONLY: destroy_preconditioner,&
64 init_preconditioner,&
65 preconditioner_type
66 USE qs_ot, ONLY: &
67 qs_ot_antihermitian_spectral_norm, qs_ot_apply_complex_frechet_dbcsr, &
68 qs_ot_complex_exp_frechet_kernel, qs_ot_density_secant_hessian, &
69 qs_ot_density_secant_orbital_overlaps, qs_ot_density_secant_projected_hessian, &
70 qs_ot_density_tangent, qs_ot_finite_rotation_response, qs_ot_fixed_n_energy_gradient, &
71 qs_ot_fixed_n_energy_hessian, qs_ot_fixed_n_multigroup_schur_block, &
72 qs_ot_fixed_n_projector_frechet, qs_ot_fixed_n_schur_block, qs_ot_generate_rotation, &
73 qs_ot_get_derivative_complex, qs_ot_get_derivative_ref_complex, &
74 qs_ot_get_orbitals_ref_complex, qs_ot_get_p_complex, qs_ot_projected_response_update, &
75 qs_ot_rot_mat_derivative, qs_ot_symmetric_abs_solve, qs_ot_symmetric_sr1_update
76 USE qs_ot_types, ONLY: qs_ot_kpoint_preconditioner_scale,&
77 qs_ot_kpoint_preconditioner_solver_supported,&
78 qs_ot_kpoint_preconditioner_supported,&
79 qs_ot_type
80 USE reference_manager, ONLY: remove_all_references
81
82 IMPLICIT NONE
83
84 INTEGER, PARAMETER :: k = 2, n = 3
85
86 COMPLEX(KIND=dp), DIMENSION(n, k) :: direction, gradient, gradient_old, hc, q, x
87 COMPLEX(KIND=dp), DIMENSION(n, n) :: hamiltonian
88 COMPLEX(KIND=dp), DIMENSION(k, k) :: b, inverse_sqrt
89 INTEGER :: io_unit, mynode, nfail
90 INTEGER, DIMENSION(:), POINTER :: col_dist, col_k, row_dist, row_k, row_n
91 INTEGER, DIMENSION(:, :), POINTER :: pgrid
92 LOGICAL :: found
93 REAL(kind=dp) :: error, fd_slope, hstep, old_error, slope
94 REAL(kind=dp), DIMENSION(:, :), POINTER :: block_p
95 TYPE(cp_logger_type), POINTER :: logger
96 TYPE(dbcsr_distribution_type) :: dist_kk, dist_nk
97 TYPE(dbcsr_type), TARGET :: b_im, b_re, f_im, f_re, g_im, g_re, &
98 hc_im, hc_re, q_im, q_re, sx_im, sx_re, &
99 tmp_kk, x_im, x_re
100 TYPE(dbcsr_type), POINTER :: hc_im_p, hc_re_p
101 TYPE(mp_comm_type) :: mp_comm
102 TYPE(mp_para_env_type), POINTER :: para_env
103 TYPE(qs_ot_type) :: ot_env
104
105 nfail = 0
106 CALL test_preconditioner_support(nfail)
107 NULLIFY (hc_im_p, hc_re_p, logger, para_env)
108 CALL mp_world_init(mp_comm)
109 mynode = mp_comm%mepos
110 io_unit = -1
111 IF (mynode == 0) io_unit = default_output_unit
112 CALL test_antihermitian_spectral_norm(io_unit, nfail)
113 CALL test_fixed_n_mermin_energy(mynode, nfail)
114 CALL test_fixed_n_projector_frechet(mynode, nfail)
115 CALL test_finite_rotation_response(mynode, nfail)
116 CALL test_fixed_n_rotation_schur(mynode, nfail)
117 CALL test_fixed_n_multigroup_schur(mynode, nfail)
118 CALL test_roks_fixed_n_rotation_schur(mynode, nfail)
119 CALL test_symmetric_abs_solve(mynode, nfail)
120 CALL test_projected_response_update(mynode, nfail)
121 CALL test_symmetric_sr1_update(mynode, nfail)
122 CALL test_density_secant_hessian(mynode, nfail)
123 CALL test_joint_spin_density_secant(mynode, nfail)
124 CALL test_density_tangent(mynode, nfail)
125 CALL test_density_secant_moving_subspace(mynode, nfail)
126 ALLOCATE (para_env)
127 CALL para_env%from_dup(mp_comm)
128 CALL cp_logger_create(logger, para_env=para_env, default_global_unit_nr=io_unit, &
129 close_global_unit_on_dealloc=.false.)
130 CALL cp_add_default_logger(logger)
131 CALL add_all_references()
132 CALL dbcsr_init_lib(mp_comm%get_handle(), io_unit)
133 CALL test_real_rotation_frechet(para_env, nfail)
134 CALL test_complex_rotation_frechet(para_env, nfail)
135 CALL test_sparse_frechet_patterns(para_env, nfail)
136 CALL test_complex_ref_rotation(para_env, nfail)
137 CALL test_split_complex_fm_gemm(para_env, io_unit, nfail)
138 CALL test_complex_preconditioner_gauge(para_env, nfail)
139 CALL test_complex_occupation_block_preconditioner(para_env, nfail)
140 CALL test_kpoint_ot_energy_weighted_density(para_env, io_unit, nfail)
141
142 x(:, 1) = [cmplx(1.10_dp, 0.10_dp, kind=dp), cmplx(0.20_dp, -0.30_dp, kind=dp), &
143 cmplx(-0.10_dp, 0.20_dp, kind=dp)]
144 x(:, 2) = [cmplx(0.30_dp, 0.40_dp, kind=dp), cmplx(0.90_dp, -0.20_dp, kind=dp), &
145 cmplx(0.25_dp, 0.10_dp, kind=dp)]
146 direction(:, 1) = [cmplx(0.17_dp, -0.09_dp, kind=dp), cmplx(-0.21_dp, 0.13_dp, kind=dp), &
147 cmplx(0.08_dp, 0.19_dp, kind=dp)]
148 direction(:, 2) = [cmplx(-0.12_dp, 0.23_dp, kind=dp), cmplx(0.16_dp, 0.07_dp, kind=dp), &
149 cmplx(-0.14_dp, -0.11_dp, kind=dp)]
150
151 hamiltonian(:, :) = cmplx(0.0_dp, 0.0_dp, kind=dp)
152 hamiltonian(1, 1) = cmplx(0.70_dp, 0.0_dp, kind=dp)
153 hamiltonian(2, 2) = cmplx(-0.40_dp, 0.0_dp, kind=dp)
154 hamiltonian(3, 3) = cmplx(1.10_dp, 0.0_dp, kind=dp)
155 hamiltonian(1, 2) = cmplx(0.20_dp, 0.10_dp, kind=dp)
156 hamiltonian(2, 1) = conjg(hamiltonian(1, 2))
157 hamiltonian(1, 3) = cmplx(-0.10_dp, 0.05_dp, kind=dp)
158 hamiltonian(3, 1) = conjg(hamiltonian(1, 3))
159 hamiltonian(2, 3) = cmplx(0.15_dp, 0.20_dp, kind=dp)
160 hamiltonian(3, 2) = conjg(hamiltonian(2, 3))
161
162 CALL dense_inverse_sqrt(matmul(conjg(transpose(x)), x), inverse_sqrt)
163 q = matmul(x, inverse_sqrt)
164 hc = 2.0_dp*matmul(hamiltonian, q)
165
166 ALLOCATE (pgrid(0:0, 0:0), row_dist(1), col_dist(1), row_n(1), row_k(1), col_k(1))
167 pgrid(:, :) = 0
168 row_dist(:) = 0
169 col_dist(:) = 0
170 row_n(:) = n
171 row_k(:) = k
172 col_k(:) = k
173 CALL dbcsr_distribution_new(dist_nk, group=mp_comm_self%get_handle(), pgrid=pgrid, &
174 row_dist=row_dist, col_dist=col_dist)
175 CALL dbcsr_distribution_new(dist_kk, group=mp_comm_self%get_handle(), pgrid=pgrid, &
176 row_dist=row_dist, col_dist=col_dist)
177 CALL test_complex_strict_derivative(dist_nk, dist_kk, row_n, row_k, col_k, para_env, nfail)
178
179 CALL create_one_block(x_re, "complex_ref_x_re", dist_nk, row_n, col_k)
180 CALL create_one_block(x_im, "complex_ref_x_im", dist_nk, row_n, col_k)
181 CALL create_one_block(sx_re, "complex_ref_sx_re", dist_nk, row_n, col_k)
182 CALL create_one_block(sx_im, "complex_ref_sx_im", dist_nk, row_n, col_k)
183 CALL create_one_block(hc_re, "complex_ref_hc_re", dist_nk, row_n, col_k)
184 CALL create_one_block(hc_im, "complex_ref_hc_im", dist_nk, row_n, col_k)
185 CALL create_one_block(g_re, "complex_ref_g_re", dist_nk, row_n, col_k)
186 CALL create_one_block(g_im, "complex_ref_g_im", dist_nk, row_n, col_k)
187 CALL create_one_block(q_re, "complex_ref_q_re", dist_nk, row_n, col_k)
188 CALL create_one_block(q_im, "complex_ref_q_im", dist_nk, row_n, col_k)
189 CALL create_one_block(f_re, "complex_ref_f_re", dist_kk, row_k, col_k)
190 CALL create_one_block(f_im, "complex_ref_f_im", dist_kk, row_k, col_k)
191 CALL create_one_block(b_re, "complex_ref_b_re", dist_kk, row_k, col_k)
192 CALL create_one_block(b_im, "complex_ref_b_im", dist_kk, row_k, col_k)
193 CALL create_one_block(tmp_kk, "complex_ref_tmp_kk", dist_kk, row_k, col_k)
194
195 CALL put_complex_pair(x_re, x_im, x)
196 CALL put_complex_pair(sx_re, sx_im, x)
197 CALL put_complex_pair(hc_re, hc_im, hc)
198 CALL put_complex_pair(f_re, f_im, inverse_sqrt)
199
200 ot_env%has_complex_kpoint_state = .true.
201 ot_env%settings%eps_irac_filter_matrix = 0.0_dp
202 ot_env%matrix_x => x_re
203 ot_env%matrix_x_im => x_im
204 ot_env%matrix_sx => sx_re
205 ot_env%matrix_sx_im => sx_im
206 ot_env%matrix_gx => g_re
207 ot_env%matrix_gx_im => g_im
208 ot_env%matrix_ref_inv_sqrt => f_re
209 ot_env%matrix_ref_inv_sqrt_im => f_im
210 ot_env%buf1_k_k_sym => b_re
211 ot_env%buf2_k_k_sym => b_im
212 ot_env%buf3_k_k_sym => tmp_kk
213 ot_env%buf1_n_k => q_re
214 ot_env%buf1_n_k_dp => q_im
215
216 hc_re_p => hc_re
217 hc_im_p => hc_im
218 CALL qs_ot_get_derivative_ref_complex(hc_re_p, hc_im_p, ot_env)
219
220 gradient(:, :) = cmplx(0.0_dp, 0.0_dp, kind=dp)
221 CALL dbcsr_get_block_p(g_re, 1, 1, block_p, found)
222 IF (found) gradient = gradient + cmplx(block_p, 0.0_dp, kind=dp)
223 CALL dbcsr_get_block_p(g_im, 1, 1, block_p, found)
224 IF (found) gradient = gradient + cmplx(0.0_dp, block_p, kind=dp)
225
226 hstep = 1.0e-6_dp
227 fd_slope = (polar_energy(x + hstep*direction, hamiltonian) - &
228 polar_energy(x - hstep*direction, hamiltonian))/(2.0_dp*hstep)
229 slope = real(sum(conjg(gradient)*direction), kind=dp)
230 error = abs(slope - fd_slope)
231
232 b = matmul(conjg(transpose(x)), hc)
233 gradient_old = hc - matmul(x, b)
234 old_error = abs(real(sum(conjg(gradient_old)*direction), kind=dp) - fd_slope)
235 IF (error > 5.0e-8_dp .OR. old_error < 1.0e-3_dp) nfail = nfail + 1
236
237 IF (mynode == 0) THEN
238 WRITE (io_unit, '(A,3(1X,ES13.6))') "complex finite-REF pullback: fd/predicted/error", &
239 fd_slope, slope, error
240 WRITE (io_unit, '(A,1X,ES13.6)') "discarded projected-gradient error", old_error
241 END IF
242
243 CALL dbcsr_release(b_im)
244 CALL dbcsr_release(b_re)
245 CALL dbcsr_release(tmp_kk)
246 CALL dbcsr_release(f_im)
247 CALL dbcsr_release(f_re)
248 CALL dbcsr_release(q_im)
249 CALL dbcsr_release(q_re)
250 CALL dbcsr_release(g_im)
251 CALL dbcsr_release(g_re)
252 CALL dbcsr_release(hc_im)
253 CALL dbcsr_release(hc_re)
254 CALL dbcsr_release(sx_im)
255 CALL dbcsr_release(sx_re)
256 CALL dbcsr_release(x_im)
257 CALL dbcsr_release(x_re)
258 CALL dbcsr_distribution_release(dist_kk)
259 CALL dbcsr_distribution_release(dist_nk)
260 DEALLOCATE (pgrid, row_dist, col_dist, row_n, row_k, col_k)
261
262 CALL dbcsr_finalize_lib()
263 CALL remove_all_references()
264 CALL cp_rm_default_logger()
265 CALL cp_logger_release(logger)
266 CALL mp_para_env_release(para_env)
267 CALL mp_world_finalize()
268
269 IF (nfail > 0) error stop "qs_ot_complex_ref_unittest failed"
270
271CONTAINS
272
273! **************************************************************************************************
274!> \brief Check split-complex full-matrix multiplication including adjoints.
275!> \param para_env parallel environment
276!> \param io_unit output unit
277!> \param nfail accumulated number of failures
278! **************************************************************************************************
279 SUBROUTINE test_split_complex_fm_gemm(para_env, io_unit, nfail)
280 TYPE(mp_para_env_type), POINTER :: para_env
281 INTEGER, INTENT(IN) :: io_unit
282 INTEGER, INTENT(INOUT) :: nfail
283
284 INTEGER, PARAMETER :: ndim = 3
285
286 COMPLEX(KIND=dp), DIMENSION(ndim, ndim) :: a, actual, b, expected
287 INTEGER :: i, j
288 REAL(KIND=dp) :: error_cn, error_nc, error_nn
289 REAL(KIND=dp), DIMENSION(ndim, ndim) :: actual_im, actual_re
290 TYPE(cp_blacs_env_type), POINTER :: blacs_env
291 TYPE(cp_fm_struct_type), POINTER :: matrix_struct
292 TYPE(cp_fm_type) :: a_im, a_re, b_im, b_re, c_im, c_re
293
294 NULLIFY (blacs_env, matrix_struct)
295 DO j = 1, ndim
296 DO i = 1, ndim
297 a(i, j) = cmplx(0.13_dp*i - 0.07_dp*j, 0.05_dp*i*j - 0.11_dp*j, kind=dp)
298 b(i, j) = cmplx(-0.09_dp*i + 0.17_dp*j, 0.08_dp*i*j + 0.04_dp*i, kind=dp)
299 END DO
300 END DO
301
302 CALL cp_blacs_env_create(blacs_env=blacs_env, para_env=para_env)
303 CALL cp_fm_struct_create(matrix_struct, nrow_global=ndim, ncol_global=ndim, &
304 context=blacs_env, para_env=para_env)
305 CALL cp_fm_create(a_re, matrix_struct)
306 CALL cp_fm_create(a_im, matrix_struct)
307 CALL cp_fm_create(b_re, matrix_struct)
308 CALL cp_fm_create(b_im, matrix_struct)
309 CALL cp_fm_create(c_re, matrix_struct)
310 CALL cp_fm_create(c_im, matrix_struct)
311 CALL cp_fm_set_submatrix(a_re, real(a, kind=dp))
312 CALL cp_fm_set_submatrix(a_im, aimag(a))
313 CALL cp_fm_set_submatrix(b_re, real(b, kind=dp))
314 CALL cp_fm_set_submatrix(b_im, aimag(b))
315
316 CALL cp_complex_fm_gemm('N', 'N', ndim, ndim, ndim, 1.0_dp, &
317 a_re, a_im, b_re, b_im, 0.0_dp, c_re, c_im)
318 CALL cp_fm_get_submatrix(c_re, actual_re)
319 CALL cp_fm_get_submatrix(c_im, actual_im)
320 actual = cmplx(actual_re, actual_im, kind=dp)
321 expected = matmul(a, b)
322 error_nn = maxval(abs(actual - expected))
323
324 CALL cp_complex_fm_gemm('C', 'N', ndim, ndim, ndim, 1.0_dp, &
325 a_re, a_im, b_re, b_im, 0.0_dp, c_re, c_im)
326 CALL cp_fm_get_submatrix(c_re, actual_re)
327 CALL cp_fm_get_submatrix(c_im, actual_im)
328 actual = cmplx(actual_re, actual_im, kind=dp)
329 expected = matmul(conjg(transpose(a)), b)
330 error_cn = maxval(abs(actual - expected))
331
332 CALL cp_complex_fm_gemm('N', 'C', ndim, ndim, ndim, 1.0_dp, &
333 a_re, a_im, b_re, b_im, 0.0_dp, c_re, c_im)
334 CALL cp_fm_get_submatrix(c_re, actual_re)
335 CALL cp_fm_get_submatrix(c_im, actual_im)
336 actual = cmplx(actual_re, actual_im, kind=dp)
337 expected = matmul(a, conjg(transpose(b)))
338 error_nc = maxval(abs(actual - expected))
339
340 IF (max(error_nn, error_cn, error_nc) > 1.0e-12_dp) nfail = nfail + 1
341 IF (io_unit >= 0) WRITE (io_unit, '(A,3(1X,ES13.6))') &
342 'split-complex GEMM N/N, C/N, N/C errors:', error_nn, error_cn, error_nc
343
344 CALL cp_fm_release(c_im)
345 CALL cp_fm_release(c_re)
346 CALL cp_fm_release(b_im)
347 CALL cp_fm_release(b_re)
348 CALL cp_fm_release(a_im)
349 CALL cp_fm_release(a_re)
350 CALL cp_fm_struct_release(matrix_struct)
351 CALL cp_blacs_env_release(blacs_env)
352
353 END SUBROUTINE test_split_complex_fm_gemm
354
355! **************************************************************************************************
356!> \brief Check the noncanonical complex OT energy-weighted density and its gauge covariance.
357!> \param para_env parallel environment
358!> \param io_unit output unit
359!> \param nfail accumulated number of failures
360! **************************************************************************************************
361 SUBROUTINE test_kpoint_ot_energy_weighted_density(para_env, io_unit, nfail)
362 TYPE(mp_para_env_type), POINTER :: para_env
363 INTEGER, INTENT(IN) :: io_unit
364 INTEGER, INTENT(INOUT) :: nfail
365
366 INTEGER, PARAMETER :: nao = 4, nmo = 3
367
368 COMPLEX(KIND=dp), DIMENSION(nao, nao) :: actual, expected, hmat, rotated
369 COMPLEX(KIND=dp), DIMENSION(nao, nmo) :: coeff, coeff_rot, hc_dense, hc_rot
370 COMPLEX(KIND=dp), DIMENSION(nmo, nmo) :: hblock, lagrange, rotation
371 INTEGER :: i, j
372 REAL(KIND=dp) :: error, gauge_error, hermitian_error, &
373 phase, pi
374 REAL(KIND=dp), DIMENSION(nao, nao) :: actual_im, actual_re, rotated_im, &
375 rotated_re
376 REAL(KIND=dp), DIMENSION(nmo) :: occupation, uniform_occupation
377 TYPE(cp_blacs_env_type), POINTER :: blacs_env
378 TYPE(cp_fm_struct_type), POINTER :: coeff_struct, hmat_struct
379 TYPE(cp_fm_type) :: coeff_im_fm, coeff_re_fm, hc_im_fm, &
380 hc_re_fm, wmat_im_fm, wmat_re_fm
381
382 NULLIFY (blacs_env, coeff_struct, hmat_struct)
383 pi = acos(-1.0_dp)
384 hmat(:, :) = cmplx(0.0_dp, 0.0_dp, kind=dp)
385 hmat(1, 1) = cmplx(-0.8_dp, 0.0_dp, kind=dp)
386 hmat(2, 2) = cmplx(0.1_dp, 0.0_dp, kind=dp)
387 hmat(3, 3) = cmplx(0.7_dp, 0.0_dp, kind=dp)
388 hmat(4, 4) = cmplx(1.4_dp, 0.0_dp, kind=dp)
389 hmat(1, 2) = cmplx(0.21_dp, 0.13_dp, kind=dp)
390 hmat(1, 3) = cmplx(-0.17_dp, 0.09_dp, kind=dp)
391 hmat(1, 4) = cmplx(0.08_dp, -0.14_dp, kind=dp)
392 hmat(2, 3) = cmplx(0.31_dp, -0.11_dp, kind=dp)
393 hmat(2, 4) = cmplx(-0.19_dp, -0.07_dp, kind=dp)
394 hmat(3, 4) = cmplx(0.27_dp, 0.16_dp, kind=dp)
395 DO i = 1, nao
396 DO j = i + 1, nao
397 hmat(j, i) = conjg(hmat(i, j))
398 END DO
399 END DO
400 DO j = 1, nmo
401 DO i = 1, nao
402 phase = 2.0_dp*pi*real((i - 1)*(j - 1), kind=dp)/real(nao, kind=dp)
403 coeff(i, j) = cmplx(cos(phase), sin(phase), kind=dp)/sqrt(real(nao, kind=dp))
404 END DO
405 END DO
406 hc_dense = matmul(hmat, coeff)
407 occupation = [1.0_dp, 0.63_dp, 0.11_dp]
408
409 hblock = matmul(conjg(transpose(coeff)), hc_dense)
410 DO j = 1, nmo
411 lagrange(:, j) = hblock(:, j)*occupation(j)
412 END DO
413 lagrange = 0.5_dp*(lagrange + conjg(transpose(lagrange)))
414 expected = matmul(matmul(coeff, lagrange), conjg(transpose(coeff)))
415
416 CALL cp_blacs_env_create(blacs_env=blacs_env, para_env=para_env)
417 CALL cp_fm_struct_create(hmat_struct, nrow_global=nao, ncol_global=nao, &
418 context=blacs_env, para_env=para_env)
419 CALL cp_fm_struct_create(coeff_struct, nrow_global=nao, ncol_global=nmo, &
420 context=blacs_env, para_env=para_env)
421 CALL cp_fm_create(coeff_re_fm, coeff_struct)
422 CALL cp_fm_create(coeff_im_fm, coeff_struct)
423 CALL cp_fm_create(hc_re_fm, coeff_struct)
424 CALL cp_fm_create(hc_im_fm, coeff_struct)
425 CALL cp_fm_create(wmat_re_fm, hmat_struct)
426 CALL cp_fm_create(wmat_im_fm, hmat_struct)
427 CALL cp_fm_set_submatrix(coeff_re_fm, real(coeff, kind=dp))
428 CALL cp_fm_set_submatrix(coeff_im_fm, aimag(coeff))
429 CALL cp_fm_set_submatrix(hc_re_fm, real(hc_dense, kind=dp))
430 CALL cp_fm_set_submatrix(hc_im_fm, aimag(hc_dense))
431
432 CALL kpoint_ot_energy_weighted_density( &
433 coeff_re_fm, coeff_im_fm, hc_re_fm, hc_im_fm, occupation, &
434 wmat_re_fm, wmat_im_fm)
435 CALL cp_fm_get_submatrix(wmat_re_fm, actual_re)
436 CALL cp_fm_get_submatrix(wmat_im_fm, actual_im)
437 actual = cmplx(actual_re, actual_im, kind=dp)
438 error = maxval(abs(actual - expected))
439 hermitian_error = maxval(abs(actual - conjg(transpose(actual))))
440
441 rotation(:, :) = cmplx(0.0_dp, 0.0_dp, kind=dp)
442 rotation(1, 1) = cmplx(sqrt(0.5_dp), 0.0_dp, kind=dp)
443 rotation(1, 2) = cmplx(0.0_dp, sqrt(0.5_dp), kind=dp)
444 rotation(2, 1) = cmplx(0.0_dp, sqrt(0.5_dp), kind=dp)
445 rotation(2, 2) = cmplx(sqrt(0.5_dp), 0.0_dp, kind=dp)
446 rotation(3, 3) = cmplx(0.0_dp, 1.0_dp, kind=dp)
447 coeff_rot = matmul(coeff, rotation)
448 hc_rot = matmul(hc_dense, rotation)
449 uniform_occupation(:) = 1.0_dp
450 CALL cp_fm_set_submatrix(coeff_re_fm, real(coeff_rot, kind=dp))
451 CALL cp_fm_set_submatrix(coeff_im_fm, aimag(coeff_rot))
452 CALL cp_fm_set_submatrix(hc_re_fm, real(hc_rot, kind=dp))
453 CALL cp_fm_set_submatrix(hc_im_fm, aimag(hc_rot))
454 CALL kpoint_ot_energy_weighted_density( &
455 coeff_re_fm, coeff_im_fm, hc_re_fm, hc_im_fm, uniform_occupation, &
456 wmat_re_fm, wmat_im_fm)
457 CALL cp_fm_get_submatrix(wmat_re_fm, rotated_re)
458 CALL cp_fm_get_submatrix(wmat_im_fm, rotated_im)
459 rotated = cmplx(rotated_re, rotated_im, kind=dp)
460 expected = matmul(matmul(coeff, hblock), conjg(transpose(coeff)))
461 gauge_error = maxval(abs(rotated - expected))
462
463 IF (error > 5.0e-12_dp .OR. hermitian_error > 5.0e-12_dp .OR. &
464 gauge_error > 5.0e-12_dp) nfail = nfail + 1
465 IF (io_unit >= 0) THEN
466 WRITE (io_unit, '(A,3(1X,ES13.6))') 'complex OT W error/hermitian/gauge', &
467 error, hermitian_error, gauge_error
468 END IF
469
470 CALL cp_fm_release(wmat_im_fm)
471 CALL cp_fm_release(wmat_re_fm)
472 CALL cp_fm_release(hc_im_fm)
473 CALL cp_fm_release(hc_re_fm)
474 CALL cp_fm_release(coeff_im_fm)
475 CALL cp_fm_release(coeff_re_fm)
476 CALL cp_fm_struct_release(coeff_struct)
477 CALL cp_fm_struct_release(hmat_struct)
478 CALL cp_blacs_env_release(blacs_env)
479
480 END SUBROUTINE test_kpoint_ot_energy_weighted_density
481! **************************************************************************************************
482!> \brief Check the complex rotation norm against an analytic pair and a unitary gauge change.
483!> \param io_unit output unit
484!> \param nfail accumulated failures
485! **************************************************************************************************
486 SUBROUTINE test_antihermitian_spectral_norm(io_unit, nfail)
487 INTEGER, INTENT(IN) :: io_unit
488 INTEGER, INTENT(INOUT) :: nfail
489
490 COMPLEX(KIND=dp), DIMENSION(3) :: phase
491 COMPLEX(KIND=dp), DIMENSION(3, 3) :: generator, transformed
492 INTEGER :: i, j
493 REAL(KIND=dp) :: error, norm, transformed_norm
494
495 generator(:, :) = cmplx(0.0_dp, 0.0_dp, kind=dp)
496 generator(1, 2) = cmplx(0.3_dp, 0.4_dp, kind=dp)
497 generator(2, 1) = -conjg(generator(1, 2))
498 phase = exp(cmplx(0.0_dp, 1.0_dp, kind=dp)*[0.37_dp, -0.51_dp, 0.83_dp])
499 DO j = 1, 3
500 DO i = 1, 3
501 transformed(i, j) = conjg(phase(i))*generator(i, j)*phase(j)
502 END DO
503 END DO
504
505 norm = qs_ot_antihermitian_spectral_norm(generator)
506 transformed_norm = qs_ot_antihermitian_spectral_norm(transformed)
507 error = max(abs(norm - 0.5_dp), abs(transformed_norm - norm))
508 IF (io_unit >= 0) WRITE (io_unit, '(A,ES14.6)') &
509 "complex rotation spectral-norm error: ", error
510 IF (error > 1.0e-13_dp) nfail = nfail + 1
511
512 END SUBROUTINE test_antihermitian_spectral_norm
513
514! **************************************************************************************************
515!> \brief Check the fixed-N Mermin gradient, dense gauge projector, and occupation Hessian.
516!> \param mynode MPI rank
517!> \param nfail accumulated number of failures
518! **************************************************************************************************
519 SUBROUTINE test_fixed_n_mermin_energy(mynode, nfail)
520 INTEGER, INTENT(IN) :: mynode
521 INTEGER, INTENT(INOUT) :: nfail
522
523 INTEGER, PARAMETER :: nbands = 4
524 REAL(KIND=dp), PARAMETER :: fd_step = 1.0e-6_dp, maxocc = 2.0_dp, &
525 target_electrons = 2.25_dp, &
526 temperature = 0.07_dp
527
528 REAL(KIND=dp) :: fd_slope, gradient_error, hessian_error, &
529 mu, predicted_slope, response_sum, &
530 weighted_residual
531 REAL(KIND=dp), DIMENSION(nbands) :: direction, energies, gradient, gradient_minus, &
532 gradient_plus, hessian_action, occupation, rayleigh, response, response_minus, &
533 response_plus, weights
534 REAL(KIND=dp), DIMENSION(nbands, nbands) :: hessian
535
536 energies(:) = [-0.24_dp, -0.05_dp, 0.08_dp, 0.31_dp]
537 rayleigh(:) = [-0.19_dp, -0.02_dp, 0.04_dp, 0.27_dp]
538 weights(:) = [0.35_dp, 0.35_dp, 0.65_dp, 0.65_dp]
539 direction(:) = [0.17_dp, -0.11_dp, 0.08_dp, -0.05_dp]
540
541 CALL fixed_n_fermi_occupations(energies, weights, target_electrons, temperature, &
542 maxocc, occupation, mu)
543 response(:) = weights(:)*occupation(:)*(maxocc - occupation(:))/(maxocc*temperature)
544 response_sum = sum(response)
545 weighted_residual = dot_product(response, rayleigh - energies)
546 CALL qs_ot_fixed_n_energy_gradient(rayleigh, energies, response, response_sum, &
547 weighted_residual, gradient)
548
549 fd_slope = (fixed_n_mermin_value(rayleigh, energies + fd_step*direction, weights, &
550 target_electrons, temperature, maxocc) - &
551 fixed_n_mermin_value(rayleigh, energies - fd_step*direction, weights, &
552 target_electrons, temperature, maxocc))/(2.0_dp*fd_step)
553 predicted_slope = dot_product(gradient, direction)
554 gradient_error = max(abs(fd_slope - predicted_slope), abs(sum(gradient)))
555
556 ! At the canonical point h_i=e_i, the exact local Hessian is the dense fixed-N projector.
557 rayleigh(:) = energies(:)
558 CALL qs_ot_fixed_n_energy_hessian(response, response_sum, hessian)
559 hessian_action(:) = matmul(hessian, direction)
560
561 CALL fixed_n_fermi_occupations(energies + fd_step*direction, weights, target_electrons, &
562 temperature, maxocc, occupation, mu)
563 response_plus(:) = weights(:)*occupation(:)*(maxocc - occupation(:))/(maxocc*temperature)
564 weighted_residual = dot_product(response_plus, rayleigh - energies - fd_step*direction)
565 CALL qs_ot_fixed_n_energy_gradient(rayleigh, energies + fd_step*direction, response_plus, &
566 sum(response_plus), weighted_residual, gradient_plus)
567
568 CALL fixed_n_fermi_occupations(energies - fd_step*direction, weights, target_electrons, &
569 temperature, maxocc, occupation, mu)
570 response_minus(:) = weights(:)*occupation(:)*(maxocc - occupation(:))/(maxocc*temperature)
571 weighted_residual = dot_product(response_minus, rayleigh - energies + fd_step*direction)
572 CALL qs_ot_fixed_n_energy_gradient(rayleigh, energies - fd_step*direction, response_minus, &
573 sum(response_minus), weighted_residual, gradient_minus)
574
575 hessian_error = maxval(abs((gradient_plus - gradient_minus)/(2.0_dp*fd_step) - &
576 hessian_action))
577 hessian_error = max(hessian_error, maxval(abs(hessian - transpose(hessian))))
578 hessian_error = max(hessian_error, maxval(abs(matmul(hessian, spread(1.0_dp, 1, nbands)))))
579
580 IF (gradient_error > 2.0e-8_dp .OR. hessian_error > 2.0e-7_dp) nfail = nfail + 1
581 IF (mynode == 0) THEN
582 WRITE (io_unit, '(A,2(1X,ES13.6))') "fixed-N Mermin gradient/Hessian errors:", &
583 gradient_error, hessian_error
584 END IF
585
586 END SUBROUTINE test_fixed_n_mermin_energy
587
588! **************************************************************************************************
589!> \brief Check the complex fixed-N spectral occupation-projector derivative by finite differences.
590!> \param mynode MPI rank
591!> \param nfail accumulated number of failures
592! **************************************************************************************************
593 SUBROUTINE test_fixed_n_projector_frechet(mynode, nfail)
594 INTEGER, INTENT(IN) :: mynode
595 INTEGER, INTENT(INOUT) :: nfail
596
597 INTEGER, PARAMETER :: nbands = 3
598 REAL(KIND=dp), PARAMETER :: fd_step = 1.0e-6_dp, kpoint_weight = 0.75_dp, maxocc = 2.0_dp, &
599 target_electrons = 2.10_dp, temperature = 0.09_dp
600
601 COMPLEX(KIND=dp), DIMENSION(nbands, nbands) :: dchc, derivative, derivative_fd, &
602 hamiltonian, projector_minus, &
603 projector_plus, vectors
604 INTEGER :: i
605 REAL(KIND=dp) :: error, mu, trace_error
606 REAL(KIND=dp), DIMENSION(nbands) :: eigenvalues, occupation, response, &
607 weights
608
609 hamiltonian(:, :) = cmplx(0.0_dp, 0.0_dp, kind=dp)
610 hamiltonian(1, 1) = cmplx(-0.22_dp, 0.0_dp, kind=dp)
611 hamiltonian(2, 2) = cmplx(0.03_dp, 0.0_dp, kind=dp)
612 hamiltonian(3, 3) = cmplx(0.28_dp, 0.0_dp, kind=dp)
613 hamiltonian(1, 2) = cmplx(0.04_dp, -0.03_dp, kind=dp)
614 hamiltonian(2, 1) = conjg(hamiltonian(1, 2))
615 hamiltonian(2, 3) = cmplx(-0.02_dp, 0.05_dp, kind=dp)
616 hamiltonian(3, 2) = conjg(hamiltonian(2, 3))
617
618 dchc(:, :) = cmplx(0.0_dp, 0.0_dp, kind=dp)
619 dchc(1, 1) = cmplx(0.13_dp, 0.0_dp, kind=dp)
620 dchc(2, 2) = cmplx(-0.08_dp, 0.0_dp, kind=dp)
621 dchc(3, 3) = cmplx(0.05_dp, 0.0_dp, kind=dp)
622 dchc(1, 2) = cmplx(0.17_dp, -0.11_dp, kind=dp)
623 dchc(2, 1) = conjg(dchc(1, 2))
624 dchc(2, 3) = cmplx(0.12_dp, 0.09_dp, kind=dp)
625 dchc(3, 2) = conjg(dchc(2, 3))
626
627 CALL diag_complex(hamiltonian, vectors, eigenvalues)
628 weights(:) = kpoint_weight
629 CALL fixed_n_fermi_occupations(eigenvalues, weights, target_electrons, temperature, &
630 maxocc, occupation, mu)
631 response(:) = weights(:)*occupation(:)*(maxocc - occupation(:))/(maxocc*temperature)
632 CALL qs_ot_fixed_n_projector_frechet(hamiltonian, dchc, occupation, kpoint_weight, &
633 response, sum(response), derivative)
634
635 CALL fixed_n_spectral_projector(hamiltonian + fd_step*dchc, weights, target_electrons, &
636 temperature, maxocc, projector_plus)
637 CALL fixed_n_spectral_projector(hamiltonian - fd_step*dchc, weights, target_electrons, &
638 temperature, maxocc, projector_minus)
639 derivative_fd(:, :) = (projector_plus - projector_minus)/(2.0_dp*fd_step)
640 error = maxval(abs(derivative - derivative_fd))
641 trace_error = 0.0_dp
642 DO i = 1, nbands
643 trace_error = trace_error + real(derivative(i, i), kind=dp)
644 END DO
645 error = max(error, abs(trace_error))
646
647 IF (error > 3.0e-7_dp) nfail = nfail + 1
648 IF (mynode == 0) WRITE (io_unit, '(A,1X,ES13.6)') &
649 "fixed-N complex occupation-projector error:", error
650
651 END SUBROUTINE test_fixed_n_projector_frechet
652
653! **************************************************************************************************
654!> \brief Check the finite complex REF rotation Hessian against the scalar band energy.
655!> \param mynode MPI rank
656!> \param nfail accumulated number of failures
657! **************************************************************************************************
658 SUBROUTINE test_finite_rotation_response(mynode, nfail)
659 INTEGER, INTENT(IN) :: mynode
660 INTEGER, INTENT(INOUT) :: nfail
661
662 INTEGER, PARAMETER :: nbands = 3, nrotation = 6
663 REAL(KIND=dp), PARAMETER :: fd_step = 2.0e-4_dp, &
664 kpoint_weight = 0.625_dp
665
666 COMPLEX(KIND=dp), DIMENSION(nbands, nbands) :: base_hamiltonian, chc, &
667 direction_generator, generator, &
668 rotation, rotation_minus, rotation_plus
669 INTEGER :: i
670 REAL(KIND=dp) :: energy, energy_minus, energy_plus, &
671 gradient_error, hessian_error, &
672 rayleigh_error, symmetry_error
673 REAL(KIND=dp), DIMENSION(nbands) :: occupation, rayleigh_minus, rayleigh_plus
674 REAL(KIND=dp), DIMENSION(nbands, nrotation) :: rayleigh_response
675 REAL(KIND=dp), DIMENSION(nrotation) :: direction, gradient
676 REAL(KIND=dp), DIMENSION(nrotation, nrotation) :: hessian
677
678 base_hamiltonian(:, :) = cmplx(0.0_dp, 0.0_dp, kind=dp)
679 base_hamiltonian(1, 1) = cmplx(-0.31_dp, 0.0_dp, kind=dp)
680 base_hamiltonian(2, 2) = cmplx(0.07_dp, 0.0_dp, kind=dp)
681 base_hamiltonian(3, 3) = cmplx(0.42_dp, 0.0_dp, kind=dp)
682 base_hamiltonian(1, 2) = cmplx(0.09_dp, -0.06_dp, kind=dp)
683 base_hamiltonian(1, 3) = cmplx(-0.04_dp, 0.08_dp, kind=dp)
684 base_hamiltonian(2, 3) = cmplx(0.11_dp, 0.05_dp, kind=dp)
685 DO i = 1, nbands
686 base_hamiltonian(i, 1:i - 1) = conjg(base_hamiltonian(1:i - 1, i))
687 END DO
688
689 generator(:, :) = cmplx(0.0_dp, 0.0_dp, kind=dp)
690 generator(1, 2) = cmplx(0.17_dp, -0.08_dp, kind=dp)
691 generator(1, 3) = cmplx(-0.09_dp, 0.04_dp, kind=dp)
692 generator(2, 3) = cmplx(0.12_dp, 0.07_dp, kind=dp)
693 DO i = 1, nbands
694 generator(i, 1:i - 1) = -conjg(generator(1:i - 1, i))
695 END DO
696
697 direction(:) = [0.14_dp, -0.07_dp, -0.11_dp, 0.09_dp, 0.05_dp, 0.13_dp]
698 direction_generator(:, :) = cmplx(0.0_dp, 0.0_dp, kind=dp)
699 direction_generator(1, 2) = cmplx(direction(1), direction(2), kind=dp)
700 direction_generator(2, 1) = -conjg(direction_generator(1, 2))
701 direction_generator(1, 3) = cmplx(direction(3), direction(4), kind=dp)
702 direction_generator(3, 1) = -conjg(direction_generator(1, 3))
703 direction_generator(2, 3) = cmplx(direction(5), direction(6), kind=dp)
704 direction_generator(3, 2) = -conjg(direction_generator(2, 3))
705 occupation(:) = [1.73_dp, 0.88_dp, 0.19_dp]
706
707 rotation(:, :) = dense_antihermitian_exp(generator)
708 chc(:, :) = matmul(conjg(transpose(rotation)), &
709 matmul(base_hamiltonian, rotation))
710 CALL qs_ot_finite_rotation_response(chc, generator, occupation, kpoint_weight, &
711 gradient, hessian, rayleigh_response)
712
713 rotation_plus(:, :) = dense_antihermitian_exp(generator + fd_step*direction_generator)
714 rotation_minus(:, :) = dense_antihermitian_exp(generator - fd_step*direction_generator)
715 energy = finite_rotation_energy(rotation, base_hamiltonian, occupation, kpoint_weight)
716 energy_plus = finite_rotation_energy(rotation_plus, base_hamiltonian, occupation, kpoint_weight)
717 energy_minus = finite_rotation_energy(rotation_minus, base_hamiltonian, occupation, kpoint_weight)
718 gradient_error = abs((energy_plus - energy_minus)/(2.0_dp*fd_step) - &
719 dot_product(gradient, direction))
720 hessian_error = abs((energy_plus - 2.0_dp*energy + energy_minus)/fd_step**2 - &
721 dot_product(direction, matmul(hessian, direction)))
722 symmetry_error = maxval(abs(hessian - transpose(hessian)))
723
724 DO i = 1, nbands
725 rayleigh_plus(i) = real(dot_product(rotation_plus(:, i), &
726 matmul(base_hamiltonian, rotation_plus(:, i))), kind=dp)
727 rayleigh_minus(i) = real(dot_product(rotation_minus(:, i), &
728 matmul(base_hamiltonian, rotation_minus(:, i))), kind=dp)
729 END DO
730 rayleigh_error = maxval(abs((rayleigh_plus - rayleigh_minus)/(2.0_dp*fd_step) - &
731 matmul(rayleigh_response, direction)))
732
733 IF (gradient_error > 5.0e-8_dp .OR. hessian_error > 2.0e-6_dp .OR. &
734 rayleigh_error > 2.0e-7_dp .OR. symmetry_error > 1.0e-12_dp) nfail = nfail + 1
735 IF (mynode == 0) THEN
736 WRITE (io_unit, '(A,4(1X,ES13.6))') "finite complex rotation response errors:", &
737 gradient_error, hessian_error, rayleigh_error, symmetry_error
738 END IF
739
740 END SUBROUTINE test_finite_rotation_response
741
742! **************************************************************************************************
743!> \brief Check the coupled fixed-N rotation/energy Hessian and its Schur decomposition.
744!> \param mynode MPI rank
745!> \param nfail accumulated number of failures
746! **************************************************************************************************
747 SUBROUTINE test_fixed_n_rotation_schur(mynode, nfail)
748 INTEGER, INTENT(IN) :: mynode
749 INTEGER, INTENT(INOUT) :: nfail
750
751 INTEGER, PARAMETER :: nbands = 3, nrotation = 6
752 REAL(KIND=dp), PARAMETER :: fd_step = 1.0e-4_dp, &
753 kpoint_weight = 0.625_dp, &
754 maxocc = 2.0_dp, temperature = 0.08_dp
755
756 COMPLEX(KIND=dp), DIMENSION(nbands, nbands) :: base_hamiltonian, chc, &
757 direction_generator, generator, &
758 rotation, rotation_minus, rotation_plus
759 INTEGER :: i
760 REAL(KIND=dp) :: curvature, curvature_error, &
761 energy_minus, energy_plus, &
762 energy_zero, mu, response_sum, &
763 schur_error
764 REAL(KIND=dp), DIMENSION(nbands) :: energy_coordinate, energy_direction, energy_gradient, &
765 energy_minus_coordinate, energy_plus_coordinate, occupation, rayleigh, rayleigh_minus, &
766 rayleigh_plus, response, weights
767 REAL(KIND=dp), DIMENSION(nbands, nbands) :: energy_hessian
768 REAL(KIND=dp), DIMENSION(nbands, nrotation) :: rayleigh_response
769 REAL(KIND=dp), DIMENSION(nrotation) :: coupling, rotation_direction, &
770 rotation_gradient, schur_rhs
771 REAL(KIND=dp), DIMENSION(nrotation, nrotation) :: rotation_hessian, schur, schur_block
772
773 base_hamiltonian(:, :) = cmplx(0.0_dp, 0.0_dp, kind=dp)
774 base_hamiltonian(1, 1) = cmplx(-0.31_dp, 0.0_dp, kind=dp)
775 base_hamiltonian(2, 2) = cmplx(0.07_dp, 0.0_dp, kind=dp)
776 base_hamiltonian(3, 3) = cmplx(0.42_dp, 0.0_dp, kind=dp)
777 base_hamiltonian(1, 2) = cmplx(0.09_dp, -0.06_dp, kind=dp)
778 base_hamiltonian(1, 3) = cmplx(-0.04_dp, 0.08_dp, kind=dp)
779 base_hamiltonian(2, 3) = cmplx(0.11_dp, 0.05_dp, kind=dp)
780 DO i = 1, nbands
781 base_hamiltonian(i, 1:i - 1) = conjg(base_hamiltonian(1:i - 1, i))
782 END DO
783
784 generator(:, :) = cmplx(0.0_dp, 0.0_dp, kind=dp)
785 generator(1, 2) = cmplx(0.17_dp, -0.08_dp, kind=dp)
786 generator(1, 3) = cmplx(-0.09_dp, 0.04_dp, kind=dp)
787 generator(2, 3) = cmplx(0.12_dp, 0.07_dp, kind=dp)
788 DO i = 1, nbands
789 generator(i, 1:i - 1) = -conjg(generator(1:i - 1, i))
790 END DO
791 rotation_direction(:) = [0.14_dp, -0.07_dp, -0.11_dp, 0.09_dp, 0.05_dp, 0.13_dp]
792 energy_direction(:) = [0.08_dp, -0.11_dp, 0.03_dp]
793 direction_generator(:, :) = cmplx(0.0_dp, 0.0_dp, kind=dp)
794 direction_generator(1, 2) = cmplx(rotation_direction(1), rotation_direction(2), kind=dp)
795 direction_generator(2, 1) = -conjg(direction_generator(1, 2))
796 direction_generator(1, 3) = cmplx(rotation_direction(3), rotation_direction(4), kind=dp)
797 direction_generator(3, 1) = -conjg(direction_generator(1, 3))
798 direction_generator(2, 3) = cmplx(rotation_direction(5), rotation_direction(6), kind=dp)
799 direction_generator(3, 2) = -conjg(direction_generator(2, 3))
800
801 rotation(:, :) = dense_antihermitian_exp(generator)
802 chc(:, :) = matmul(conjg(transpose(rotation)), matmul(base_hamiltonian, rotation))
803 DO i = 1, nbands
804 rayleigh(i) = real(chc(i, i), kind=dp)
805 END DO
806 energy_coordinate(:) = rayleigh(:)
807 weights(:) = kpoint_weight
808 CALL fixed_n_fermi_occupations(energy_coordinate, weights, 1.7_dp, temperature, &
809 maxocc, occupation, mu)
810 response(:) = weights*occupation*(maxocc - occupation)/(maxocc*temperature)
811 response_sum = sum(response)
812 CALL qs_ot_finite_rotation_response(chc, generator, occupation, kpoint_weight, &
813 rotation_gradient, rotation_hessian, rayleigh_response)
814 CALL qs_ot_fixed_n_energy_gradient(rayleigh, energy_coordinate, response, response_sum, &
815 0.0_dp, energy_gradient)
816 CALL qs_ot_fixed_n_energy_hessian(response, response_sum, energy_hessian)
817 CALL qs_ot_fixed_n_schur_block( &
818 rotation_hessian, rayleigh_response, response, rotation_gradient, energy_gradient, &
819 schur_block, coupling, schur_rhs)
820 schur(:, :) = schur_block + spread(coupling, dim=2, ncopies=nrotation)* &
821 spread(coupling, dim=1, ncopies=nrotation)/response_sum
822 schur_error = maxval(abs(schur - rotation_hessian + &
823 matmul(transpose(rayleigh_response), &
824 matmul(energy_hessian, rayleigh_response))))
825
826 rotation_plus(:, :) = dense_antihermitian_exp(generator + fd_step*direction_generator)
827 rotation_minus(:, :) = dense_antihermitian_exp(generator - fd_step*direction_generator)
828 DO i = 1, nbands
829 rayleigh_plus(i) = real(dot_product(rotation_plus(:, i), &
830 matmul(base_hamiltonian, rotation_plus(:, i))), kind=dp)
831 rayleigh_minus(i) = real(dot_product(rotation_minus(:, i), &
832 matmul(base_hamiltonian, rotation_minus(:, i))), kind=dp)
833 END DO
834 energy_plus_coordinate(:) = energy_coordinate + fd_step*energy_direction
835 energy_minus_coordinate(:) = energy_coordinate - fd_step*energy_direction
836 energy_zero = fixed_n_mermin_value(rayleigh, energy_coordinate, weights, 1.7_dp, &
837 temperature, maxocc)
838 energy_plus = fixed_n_mermin_value(rayleigh_plus, energy_plus_coordinate, weights, 1.7_dp, &
839 temperature, maxocc)
840 energy_minus = fixed_n_mermin_value(rayleigh_minus, energy_minus_coordinate, weights, 1.7_dp, &
841 temperature, maxocc)
842 curvature = dot_product(rotation_direction, matmul(rotation_hessian, rotation_direction)) - &
843 2.0_dp*dot_product(energy_direction, &
844 matmul(energy_hessian, &
845 matmul(rayleigh_response, rotation_direction))) + &
846 dot_product(energy_direction, matmul(energy_hessian, energy_direction))
847 curvature_error = abs((energy_plus - 2.0_dp*energy_zero + energy_minus)/fd_step**2 - curvature)
848
849 IF (schur_error > 2.0e-11_dp .OR. curvature_error > 3.0e-6_dp) nfail = nfail + 1
850 IF (mynode == 0) THEN
851 WRITE (io_unit, '(A,2(1X,ES13.6))') "fixed-N rotation Schur errors:", &
852 schur_error, curvature_error
853 END IF
854
855 END SUBROUTINE test_fixed_n_rotation_schur
856
857! **************************************************************************************************
858!> \brief verify the two-chemical-potential Schur reduction used by shared ROKS orbitals
859!> \param mynode MPI rank
860!> \param nfail accumulated failures
861! **************************************************************************************************
862 SUBROUTINE test_fixed_n_multigroup_schur(mynode, nfail)
863 INTEGER, INTENT(IN) :: mynode
864 INTEGER, INTENT(INOUT) :: nfail
865
866 INTEGER, PARAMETER :: nenergy = 4, ngroups = 2, nrotation = 3
867
868 INTEGER, DIMENSION(nenergy) :: response_group
869 REAL(KIND=dp) :: coupling_error, rhs_error, schur_error
870 REAL(KIND=dp), DIMENSION(nenergy) :: energy_gradient, response_weight
871 REAL(KIND=dp), DIMENSION(nenergy, nrotation) :: rayleigh_response
872 REAL(KIND=dp), DIMENSION(nrotation) :: rotation_gradient, schur_rhs
873 REAL(KIND=dp), DIMENSION(nrotation, ngroups) :: coupling
874 REAL(KIND=dp), DIMENSION(nrotation, nrotation) :: expected, rotation_hessian, schur_block
875
876 rotation_hessian(:, :) = reshape( &
877 [2.3_dp, 0.2_dp, -0.1_dp, 0.2_dp, 1.8_dp, 0.3_dp, -0.1_dp, 0.3_dp, 1.5_dp], &
878 [nrotation, nrotation])
879 rayleigh_response(:, :) = reshape( &
880 [0.4_dp, -0.2_dp, 0.1_dp, 0.3_dp, &
881 -0.1_dp, 0.5_dp, 0.2_dp, -0.4_dp, &
882 0.3_dp, 0.1_dp, -0.5_dp, 0.2_dp], &
883 [nenergy, nrotation])
884 response_weight(:) = [0.7_dp, 0.5_dp, 0.4_dp, 0.2_dp]
885 response_group(:) = [1, 1, 2, 2]
886 rotation_gradient(:) = [0.2_dp, -0.1_dp, 0.3_dp]
887 energy_gradient(:) = [0.04_dp, -0.03_dp, 0.02_dp, -0.01_dp]
888
889 CALL qs_ot_fixed_n_multigroup_schur_block( &
890 rotation_hessian, rayleigh_response, response_weight, response_group, &
891 rotation_gradient, energy_gradient, schur_block, coupling, schur_rhs)
892
893 expected(:, :) = rotation_hessian - &
894 matmul(transpose(rayleigh_response), &
895 matmul(reshape( &
896 [0.2916666666666667_dp, -0.2916666666666667_dp, 0.0_dp, 0.0_dp, &
897 -0.2916666666666667_dp, 0.2916666666666667_dp, 0.0_dp, 0.0_dp, &
898 0.0_dp, 0.0_dp, 0.1333333333333333_dp, -0.1333333333333333_dp, &
899 0.0_dp, 0.0_dp, -0.1333333333333333_dp, 0.1333333333333333_dp], &
900 [nenergy, nenergy]), rayleigh_response))
901 schur_error = maxval(abs( &
902 schur_block + &
903 matmul(coupling(:, 1:1), transpose(coupling(:, 1:1)))/1.2_dp + &
904 matmul(coupling(:, 2:2), transpose(coupling(:, 2:2)))/0.6_dp - expected))
905 coupling_error = maxval(abs(coupling(:, 1) - &
906 matmul(transpose(rayleigh_response(1:2, :)), &
907 response_weight(1:2))))
908 coupling_error = max(coupling_error, maxval(abs(coupling(:, 2) - &
909 matmul(transpose(rayleigh_response(3:4, :)), &
910 response_weight(3:4)))))
911 rhs_error = maxval(abs(schur_rhs - rotation_gradient - &
912 matmul(transpose(rayleigh_response), energy_gradient)))
913
914 IF (max(schur_error, coupling_error, rhs_error) > 2.0e-14_dp) nfail = nfail + 1
915 IF (mynode == 0) THEN
916 WRITE (io_unit, '(A,3(1X,ES13.6))') "fixed-N multigroup Schur errors:", &
917 schur_error, coupling_error, rhs_error
918 END IF
919
920 END SUBROUTINE test_fixed_n_multigroup_schur
921
922! **************************************************************************************************
923!> \brief finite-difference the shared-orbital, two-fixed-N ROKS Mermin Hessian
924!> \param mynode MPI rank
925!> \param nfail accumulated failures
926! **************************************************************************************************
927 SUBROUTINE test_roks_fixed_n_rotation_schur(mynode, nfail)
928 INTEGER, INTENT(IN) :: mynode
929 INTEGER, INTENT(INOUT) :: nfail
930
931 INTEGER, PARAMETER :: nbands = 3, nenergy = 6, ngroups = 2, &
932 nrotation = 6
933 REAL(KIND=dp), PARAMETER :: fd_step = 1.0e-4_dp, &
934 kpoint_weight = 0.7_dp, &
935 maxocc = 1.0_dp, temperature = 0.09_dp
936
937 COMPLEX(KIND=dp), DIMENSION(nbands, nbands) :: chc_a, chc_b, direction_generator, generator, &
938 hamiltonian_a, hamiltonian_b, rotation, rotation_minus, rotation_plus
939 INTEGER :: i
940 INTEGER, DIMENSION(nenergy) :: response_group
941 REAL(KIND=dp) :: curvature, curvature_a, curvature_b, curvature_error, energy_minus, &
942 energy_minus_a, energy_minus_b, energy_plus, energy_plus_a, energy_plus_b, energy_zero, &
943 energy_zero_a, energy_zero_b, mu_a, mu_b, response_sum_a, response_sum_b, schur_error
944 REAL(KIND=dp), DIMENSION(nbands) :: energy_a, energy_a_minus, energy_a_plus, energy_b, &
945 energy_b_minus, energy_b_plus, occupation_a, occupation_b, rayleigh_a, rayleigh_a_minus, &
946 rayleigh_a_plus, rayleigh_b, rayleigh_b_minus, rayleigh_b_plus, response_a, response_b, &
947 weights
948 REAL(KIND=dp), DIMENSION(nbands, nrotation) :: spin_rayleigh_response
949 REAL(KIND=dp), DIMENSION(nenergy) :: energy_direction, energy_gradient, &
950 response
951 REAL(KIND=dp), DIMENSION(nenergy, nenergy) :: energy_hessian
952 REAL(KIND=dp), DIMENSION(nenergy, nrotation) :: rayleigh_response
953 REAL(KIND=dp), DIMENSION(nrotation) :: rotation_direction, rotation_gradient, &
954 schur_rhs, spin_gradient
955 REAL(KIND=dp), DIMENSION(nrotation, ngroups) :: coupling
956 REAL(KIND=dp), DIMENSION(nrotation, nrotation) :: rotation_hessian, rotation_hessian_a, &
957 rotation_hessian_b, schur, &
958 schur_block, spin_hessian
959
960 hamiltonian_a(:, :) = cmplx(0.0_dp, 0.0_dp, kind=dp)
961 hamiltonian_b(:, :) = cmplx(0.0_dp, 0.0_dp, kind=dp)
962 hamiltonian_a(1, 1) = -0.42_dp
963 hamiltonian_a(2, 2) = 0.03_dp
964 hamiltonian_a(3, 3) = 0.37_dp
965 hamiltonian_a(1, 2) = cmplx(0.08_dp, -0.05_dp, kind=dp)
966 hamiltonian_a(1, 3) = cmplx(-0.03_dp, 0.07_dp, kind=dp)
967 hamiltonian_a(2, 3) = cmplx(0.09_dp, 0.04_dp, kind=dp)
968 hamiltonian_b(1, 1) = -0.31_dp
969 hamiltonian_b(2, 2) = 0.11_dp
970 hamiltonian_b(3, 3) = 0.46_dp
971 hamiltonian_b(1, 2) = cmplx(-0.06_dp, 0.02_dp, kind=dp)
972 hamiltonian_b(1, 3) = cmplx(0.05_dp, 0.04_dp, kind=dp)
973 hamiltonian_b(2, 3) = cmplx(-0.07_dp, 0.06_dp, kind=dp)
974 DO i = 1, nbands
975 hamiltonian_a(i, 1:i - 1) = conjg(hamiltonian_a(1:i - 1, i))
976 hamiltonian_b(i, 1:i - 1) = conjg(hamiltonian_b(1:i - 1, i))
977 END DO
978
979 generator(:, :) = cmplx(0.0_dp, 0.0_dp, kind=dp)
980 generator(1, 2) = cmplx(0.14_dp, -0.06_dp, kind=dp)
981 generator(1, 3) = cmplx(-0.08_dp, 0.05_dp, kind=dp)
982 generator(2, 3) = cmplx(0.10_dp, 0.03_dp, kind=dp)
983 DO i = 1, nbands
984 generator(i, 1:i - 1) = -conjg(generator(1:i - 1, i))
985 END DO
986 rotation_direction(:) = [0.09_dp, -0.04_dp, -0.07_dp, 0.06_dp, 0.03_dp, 0.08_dp]
987 energy_direction(:) = [0.05_dp, -0.04_dp, 0.02_dp, -0.03_dp, 0.06_dp, -0.01_dp]
988 direction_generator(:, :) = cmplx(0.0_dp, 0.0_dp, kind=dp)
989 direction_generator(1, 2) = cmplx(rotation_direction(1), rotation_direction(2), kind=dp)
990 direction_generator(2, 1) = -conjg(direction_generator(1, 2))
991 direction_generator(1, 3) = cmplx(rotation_direction(3), rotation_direction(4), kind=dp)
992 direction_generator(3, 1) = -conjg(direction_generator(1, 3))
993 direction_generator(2, 3) = cmplx(rotation_direction(5), rotation_direction(6), kind=dp)
994 direction_generator(3, 2) = -conjg(direction_generator(2, 3))
995
996 rotation(:, :) = dense_antihermitian_exp(generator)
997 chc_a(:, :) = matmul(conjg(transpose(rotation)), matmul(hamiltonian_a, rotation))
998 chc_b(:, :) = matmul(conjg(transpose(rotation)), matmul(hamiltonian_b, rotation))
999 DO i = 1, nbands
1000 rayleigh_a(i) = real(chc_a(i, i), kind=dp)
1001 rayleigh_b(i) = real(chc_b(i, i), kind=dp)
1002 END DO
1003 energy_a(:) = rayleigh_a
1004 energy_b(:) = rayleigh_b
1005 weights(:) = kpoint_weight
1006 CALL fixed_n_fermi_occupations(energy_a, weights, 1.65_dp, temperature, &
1007 maxocc, occupation_a, mu_a)
1008 CALL fixed_n_fermi_occupations(energy_b, weights, 0.95_dp, temperature, &
1009 maxocc, occupation_b, mu_b)
1010 response_a(:) = weights*occupation_a*(maxocc - occupation_a)/(maxocc*temperature)
1011 response_b(:) = weights*occupation_b*(maxocc - occupation_b)/(maxocc*temperature)
1012 response_sum_a = sum(response_a)
1013 response_sum_b = sum(response_b)
1014 response(1:nbands) = response_a
1015 response(nbands + 1:nenergy) = response_b
1016 response_group(:) = [1, 1, 1, 2, 2, 2]
1017
1018 CALL qs_ot_finite_rotation_response( &
1019 chc_a, generator, occupation_a, kpoint_weight, spin_gradient, spin_hessian, &
1020 spin_rayleigh_response)
1021 rotation_gradient(:) = spin_gradient
1022 rotation_hessian_a(:, :) = spin_hessian
1023 rayleigh_response(1:nbands, :) = spin_rayleigh_response
1024 CALL qs_ot_finite_rotation_response( &
1025 chc_b, generator, occupation_b, kpoint_weight, spin_gradient, spin_hessian, &
1026 spin_rayleigh_response)
1027 rotation_gradient(:) = rotation_gradient + spin_gradient
1028 rotation_hessian_b(:, :) = spin_hessian
1029 rotation_hessian(:, :) = rotation_hessian_a + rotation_hessian_b
1030 rayleigh_response(nbands + 1:nenergy, :) = spin_rayleigh_response
1031
1032 CALL qs_ot_fixed_n_energy_gradient( &
1033 rayleigh_a, energy_a, response_a, response_sum_a, 0.0_dp, &
1034 energy_gradient(1:nbands))
1035 CALL qs_ot_fixed_n_energy_gradient( &
1036 rayleigh_b, energy_b, response_b, response_sum_b, 0.0_dp, &
1037 energy_gradient(nbands + 1:nenergy))
1038 energy_hessian(:, :) = 0.0_dp
1039 CALL qs_ot_fixed_n_energy_hessian( &
1040 response_a, response_sum_a, energy_hessian(1:nbands, 1:nbands))
1041 CALL qs_ot_fixed_n_energy_hessian( &
1042 response_b, response_sum_b, energy_hessian(nbands + 1:nenergy, nbands + 1:nenergy))
1043 CALL qs_ot_fixed_n_multigroup_schur_block( &
1044 rotation_hessian, rayleigh_response, response, response_group, rotation_gradient, &
1045 energy_gradient, schur_block, coupling, schur_rhs)
1046 schur(:, :) = schur_block + &
1047 matmul(coupling(:, 1:1), transpose(coupling(:, 1:1)))/response_sum_a + &
1048 matmul(coupling(:, 2:2), transpose(coupling(:, 2:2)))/response_sum_b
1049 schur_error = maxval(abs(schur - rotation_hessian + &
1050 matmul(transpose(rayleigh_response), &
1051 matmul(energy_hessian, rayleigh_response))))
1052
1053 rotation_plus(:, :) = dense_antihermitian_exp(generator + fd_step*direction_generator)
1054 rotation_minus(:, :) = dense_antihermitian_exp(generator - fd_step*direction_generator)
1055 DO i = 1, nbands
1056 rayleigh_a_plus(i) = real(dot_product(rotation_plus(:, i), &
1057 matmul(hamiltonian_a, rotation_plus(:, i))), kind=dp)
1058 rayleigh_a_minus(i) = real(dot_product(rotation_minus(:, i), &
1059 matmul(hamiltonian_a, rotation_minus(:, i))), kind=dp)
1060 rayleigh_b_plus(i) = real(dot_product(rotation_plus(:, i), &
1061 matmul(hamiltonian_b, rotation_plus(:, i))), kind=dp)
1062 rayleigh_b_minus(i) = real(dot_product(rotation_minus(:, i), &
1063 matmul(hamiltonian_b, rotation_minus(:, i))), kind=dp)
1064 END DO
1065 energy_a_plus(:) = energy_a + fd_step*energy_direction(1:nbands)
1066 energy_a_minus(:) = energy_a - fd_step*energy_direction(1:nbands)
1067 energy_b_plus(:) = energy_b + fd_step*energy_direction(nbands + 1:nenergy)
1068 energy_b_minus(:) = energy_b - fd_step*energy_direction(nbands + 1:nenergy)
1069 energy_zero_a = fixed_n_mermin_value( &
1070 rayleigh_a, energy_a, weights, 1.65_dp, temperature, maxocc)
1071 energy_zero_b = fixed_n_mermin_value( &
1072 rayleigh_b, energy_b, weights, 0.95_dp, temperature, maxocc)
1073 energy_plus_a = fixed_n_mermin_value( &
1074 rayleigh_a_plus, energy_a_plus, weights, 1.65_dp, temperature, maxocc)
1075 energy_plus_b = fixed_n_mermin_value( &
1076 rayleigh_b_plus, energy_b_plus, weights, 0.95_dp, temperature, maxocc)
1077 energy_minus_a = fixed_n_mermin_value( &
1078 rayleigh_a_minus, energy_a_minus, weights, 1.65_dp, temperature, maxocc)
1079 energy_minus_b = fixed_n_mermin_value( &
1080 rayleigh_b_minus, energy_b_minus, weights, 0.95_dp, temperature, maxocc)
1081 energy_zero = energy_zero_a + energy_zero_b
1082 energy_plus = energy_plus_a + energy_plus_b
1083 energy_minus = energy_minus_a + energy_minus_b
1084 curvature_a = dot_product(rotation_direction, &
1085 matmul(rotation_hessian_a, rotation_direction)) - &
1086 2.0_dp*dot_product(energy_direction(1:nbands), &
1087 matmul(energy_hessian(1:nbands, 1:nbands), &
1088 matmul(rayleigh_response(1:nbands, :), &
1089 rotation_direction))) + &
1090 dot_product(energy_direction(1:nbands), &
1091 matmul(energy_hessian(1:nbands, 1:nbands), &
1092 energy_direction(1:nbands)))
1093 curvature_b = dot_product(rotation_direction, &
1094 matmul(rotation_hessian_b, rotation_direction)) - &
1095 2.0_dp*dot_product(energy_direction(nbands + 1:nenergy), &
1096 matmul(energy_hessian(nbands + 1:nenergy, &
1097 nbands + 1:nenergy), &
1098 matmul(rayleigh_response(nbands + 1:nenergy, :), &
1099 rotation_direction))) + &
1100 dot_product(energy_direction(nbands + 1:nenergy), &
1101 matmul(energy_hessian(nbands + 1:nenergy, nbands + 1:nenergy), &
1102 energy_direction(nbands + 1:nenergy)))
1103 curvature = curvature_a + curvature_b
1104 curvature_error = abs((energy_plus_a - 2.0_dp*energy_zero_a + energy_minus_a)/ &
1105 fd_step**2 - curvature_a)
1106 curvature_error = max( &
1107 curvature_error, &
1108 abs((energy_plus_b - 2.0_dp*energy_zero_b + energy_minus_b)/fd_step**2 - curvature_b), &
1109 abs((energy_plus - 2.0_dp*energy_zero + energy_minus)/fd_step**2 - curvature))
1110
1111 IF (schur_error > 3.0e-11_dp .OR. curvature_error > 5.0e-6_dp) nfail = nfail + 1
1112 IF (mynode == 0) THEN
1113 WRITE (io_unit, '(A,2(1X,ES13.6))') "ROKS fixed-N rotation Schur errors:", &
1114 schur_error, curvature_error
1115 END IF
1116
1117 END SUBROUTINE test_roks_fixed_n_rotation_schur
1118
1119! **************************************************************************************************
1120!> \brief Check that positive and negative response modes retain their spectral magnitude.
1121!> \param mynode MPI rank
1122!> \param nfail accumulated number of failures
1123! **************************************************************************************************
1124 SUBROUTINE test_symmetric_abs_solve(mynode, nfail)
1125 INTEGER, INTENT(IN) :: mynode
1126 INTEGER, INTENT(INOUT) :: nfail
1127
1128 LOGICAL :: valid
1129 REAL(KIND=dp) :: error
1130 REAL(KIND=dp), DIMENSION(3) :: eigenvalues
1131 REAL(KIND=dp), DIMENSION(3, 2) :: expected, rhs, solution
1132 REAL(KIND=dp), DIMENSION(3, 3) :: absolute_matrix, matrix, vectors
1133
1134 vectors(:, 1) = [0.8_dp, 0.6_dp, 0.0_dp]
1135 vectors(:, 2) = [-0.6_dp, 0.8_dp, 0.0_dp]
1136 vectors(:, 3) = [0.0_dp, 0.0_dp, 1.0_dp]
1137 eigenvalues(:) = [-0.4_dp, 1.7_dp, 0.25_dp]
1138 matrix(:, :) = matmul(vectors, matmul( &
1139 reshape([eigenvalues(1), 0.0_dp, 0.0_dp, 0.0_dp, eigenvalues(2), 0.0_dp, &
1140 0.0_dp, 0.0_dp, eigenvalues(3)], [3, 3]), transpose(vectors)))
1141 absolute_matrix(:, :) = matmul(vectors, matmul( &
1142 reshape([abs(eigenvalues(1)), 0.0_dp, 0.0_dp, 0.0_dp, abs(eigenvalues(2)), 0.0_dp, &
1143 0.0_dp, 0.0_dp, abs(eigenvalues(3))], [3, 3]), transpose(vectors)))
1144 rhs(:, 1) = [0.3_dp, -0.5_dp, 0.7_dp]
1145 rhs(:, 2) = [-0.2_dp, 0.4_dp, 0.1_dp]
1146 CALL qs_ot_symmetric_abs_solve(matrix, rhs, solution, valid)
1147 expected(:, :) = matmul(absolute_matrix, solution)
1148 error = maxval(abs(expected - rhs))
1149 IF (.NOT. valid .OR. error > 2.0e-13_dp) nfail = nfail + 1
1150
1151 eigenvalues(:) = [-0.4_dp, 1.7_dp, 1.0e-12_dp]
1152 matrix(:, :) = matmul(vectors, matmul( &
1153 reshape([eigenvalues(1), 0.0_dp, 0.0_dp, 0.0_dp, eigenvalues(2), 0.0_dp, &
1154 0.0_dp, 0.0_dp, eigenvalues(3)], [3, 3]), transpose(vectors)))
1155 rhs(:, 1) = [0.3_dp, -0.5_dp, 0.7_dp]
1156 rhs(:, 2) = [-0.2_dp, 0.4_dp, 0.1_dp]
1157 CALL qs_ot_symmetric_abs_solve(matrix, rhs, solution, valid)
1158 expected(:, :) = rhs(:, :)
1159 expected(3, :) = 0.0_dp
1160 absolute_matrix(:, :) = matmul(vectors(:, 1:2), matmul( &
1161 reshape([abs(eigenvalues(1)), 0.0_dp, 0.0_dp, abs(eigenvalues(2))], [2, 2]), &
1162 transpose(vectors(:, 1:2))))
1163 error = max(error, maxval(abs(matmul(absolute_matrix, solution) - expected)))
1164 IF (.NOT. valid .OR. error > 2.0e-13_dp .OR. maxval(abs(solution(3, :))) > 2.0e-13_dp) THEN
1165 nfail = nfail + 1
1166 END IF
1167 IF (mynode == 0) WRITE (io_unit, '(A,1X,ES13.6)') "symmetric absolute solve error:", error
1168
1169 END SUBROUTINE test_symmetric_abs_solve
1170
1171! **************************************************************************************************
1172!> \brief Check that the projected response update preserves the baseline and rejects indefiniteness.
1173!> \param mynode MPI rank
1174!> \param nfail accumulated number of failures
1175! **************************************************************************************************
1176 SUBROUTINE test_projected_response_update(mynode, nfail)
1177 INTEGER, INTENT(IN) :: mynode
1178 INTEGER, INTENT(INOUT) :: nfail
1179
1180 LOGICAL :: valid
1181 REAL(KIND=dp) :: error, preservation_error
1182 REAL(KIND=dp), DIMENSION(2) :: coefficients, gradient, residual
1183 REAL(KIND=dp), DIMENSION(2, 2) :: correction, reference, total
1184
1185 reference(:, :) = reshape([4.0_dp, 0.7_dp, 0.7_dp, 2.5_dp], [2, 2])
1186 correction(:, :) = 0.0_dp
1187 CALL qs_ot_projected_response_update(reference, correction, coefficients, valid)
1188 preservation_error = maxval(abs(coefficients - [1.0_dp, 0.0_dp]))
1189 IF (.NOT. valid .OR. preservation_error > 1.0e-12_dp) nfail = nfail + 1
1190
1191 correction(:, :) = reshape([0.8_dp, -0.3_dp, -0.3_dp, 0.4_dp], [2, 2])
1192 gradient(:) = [1.2_dp, -0.4_dp]
1193 CALL qs_ot_projected_response_update( &
1194 reference, correction, coefficients, valid, projected_gradient=gradient)
1195 total(:, :) = reference + correction
1196 residual(:) = matmul(total, coefficients) - gradient
1197 error = maxval(abs(residual))
1198 IF (.NOT. valid .OR. error > 1.0e-12_dp) nfail = nfail + 1
1199
1200 correction(:, :) = reshape([-8.0_dp, 0.0_dp, 0.0_dp, 0.0_dp], [2, 2])
1201 CALL qs_ot_projected_response_update(reference, correction, coefficients, valid)
1202 IF (valid) nfail = nfail + 1
1203 IF (mynode == 0) WRITE (io_unit, '(A,2(1X,ES14.6))') &
1204 "projected response preservation/residual errors:", preservation_error, error
1205
1206 END SUBROUTINE test_projected_response_update
1207
1208! **************************************************************************************************
1209!> \brief Check exact signed SR1 secants and rejection of an unresolved orthogonal residual.
1210!> \param mynode MPI rank
1211!> \param nfail accumulated number of failures
1212! **************************************************************************************************
1213 SUBROUTINE test_symmetric_sr1_update(mynode, nfail)
1214 INTEGER, INTENT(IN) :: mynode
1215 INTEGER, INTENT(INOUT) :: nfail
1216
1217 LOGICAL :: valid
1218 REAL(KIND=dp) :: rejected_error, secant_error, &
1219 symmetry_error
1220 REAL(KIND=dp), DIMENSION(3) :: orthogonal_residual, response, step
1221 REAL(KIND=dp), DIMENSION(3, 3) :: matrix, updated
1222
1223 matrix(:, :) = 0.0_dp
1224 matrix(1, 1) = 1.0_dp
1225 matrix(2, 2) = 2.0_dp
1226 matrix(3, 3) = 3.0_dp
1227 step(:) = [0.3_dp, -0.4_dp, 0.2_dp]
1228 response(:) = [0.0_dp, -0.3_dp, 0.8_dp]
1229 CALL qs_ot_symmetric_sr1_update(matrix, step, response, updated, valid)
1230 secant_error = maxval(abs(matmul(updated, step) - response))
1231 symmetry_error = maxval(abs(updated - transpose(updated)))
1232 IF (.NOT. valid .OR. secant_error > 2.0e-14_dp .OR. symmetry_error > 2.0e-14_dp) THEN
1233 nfail = nfail + 1
1234 END IF
1235
1236 orthogonal_residual(:) = [0.4_dp, 0.3_dp, 0.0_dp]
1237 response(:) = matmul(matrix, step) + orthogonal_residual
1238 CALL qs_ot_symmetric_sr1_update(matrix, step, response, updated, valid)
1239 rejected_error = maxval(abs(updated - matrix))
1240 IF (valid .OR. rejected_error > 2.0e-14_dp) nfail = nfail + 1
1241 IF (mynode == 0) WRITE (io_unit, '(A,3(1X,ES13.6))') &
1242 "symmetric SR1 secant/symmetry/rejection errors:", &
1243 secant_error, symmetry_error, rejected_error
1244
1245 END SUBROUTINE test_symmetric_sr1_update
1246
1247! **************************************************************************************************
1248!> \brief Check the complex gauge-invariant, sign-aware density-space response secant.
1249!> \param mynode MPI rank
1250!> \param nfail accumulated number of failures
1251! **************************************************************************************************
1252 SUBROUTINE test_density_secant_hessian(mynode, nfail)
1253 INTEGER, INTENT(IN) :: mynode
1254 INTEGER, INTENT(INOUT) :: nfail
1255
1256 INTEGER, PARAMETER :: n = 3, nmode = 2
1257 REAL(KIND=dp), PARAMETER :: secant_position = 0.4_dp
1258
1259 COMPLEX(KIND=dp), DIMENSION(n, n) :: density_step, hamiltonian_step, unitary
1260 COMPLEX(KIND=dp), DIMENSION(n, n, nmode) :: density_modes, transformed_modes
1261 LOGICAL :: calibrated_valid, rejected_valid, valid
1262 REAL(KIND=dp) :: calibrated_error, density_norm, direct_error, gauge_error, overlap_error, &
1263 rejected_error, response_work, sign_error
1264 REAL(KIND=dp), DIMENSION(nmode) :: density_overlap, response_overlap
1265 REAL(KIND=dp), DIMENSION(nmode, nmode) :: calibrated, correction, &
1266 correction_gauge, correction_overlap, &
1267 expected, rejected
1268
1269 density_step(:, :) = cmplx(0.0_dp, 0.0_dp, kind=dp)
1270 density_step(1, 1) = cmplx(0.30_dp, 0.0_dp, kind=dp)
1271 density_step(2, 2) = cmplx(-0.18_dp, 0.0_dp, kind=dp)
1272 density_step(3, 3) = cmplx(-0.12_dp, 0.0_dp, kind=dp)
1273 density_step(1, 2) = cmplx(0.07_dp, -0.11_dp, kind=dp)
1274 density_step(2, 1) = conjg(density_step(1, 2))
1275 density_step(1, 3) = cmplx(-0.04_dp, 0.05_dp, kind=dp)
1276 density_step(3, 1) = conjg(density_step(1, 3))
1277
1278 hamiltonian_step(:, :) = -0.45_dp*density_step
1279 hamiltonian_step(2, 3) = cmplx(0.06_dp, 0.03_dp, kind=dp)
1280 hamiltonian_step(3, 2) = conjg(hamiltonian_step(2, 3))
1281 density_modes(:, :, 1) = density_step
1282 density_modes(:, :, 2) = cmplx(0.0_dp, 0.0_dp, kind=dp)
1283 density_modes(1, 1, 2) = cmplx(-0.2_dp, 0.0_dp, kind=dp)
1284 density_modes(2, 2, 2) = cmplx(0.1_dp, 0.0_dp, kind=dp)
1285 density_modes(3, 3, 2) = cmplx(0.1_dp, 0.0_dp, kind=dp)
1286 density_modes(1, 3, 2) = cmplx(0.03_dp, 0.08_dp, kind=dp)
1287 density_modes(3, 1, 2) = conjg(density_modes(1, 3, 2))
1288
1289 CALL qs_ot_density_secant_hessian( &
1290 density_step, hamiltonian_step, density_modes, correction, valid, &
1291 density_norm, response_work)
1292 density_overlap(1) = density_norm
1293 density_overlap(2) = sum(real(conjg(density_step)*density_modes(:, :, 2), kind=dp))
1294 response_overlap(1) = response_work
1295 response_overlap(2) = sum(real(conjg(hamiltonian_step)*density_modes(:, :, 2), kind=dp))
1296 CALL qs_ot_density_secant_projected_hessian( &
1297 density_norm, response_work, density_overlap, response_overlap, correction_overlap, valid)
1298
1299 expected(1, 1) = response_work
1300 expected(1, 2) = response_overlap(2)
1301 expected(2, 1) = response_overlap(2)
1302 expected(2, 2) = (2.0_dp*density_overlap(2)*response_overlap(2) - &
1303 response_work*density_overlap(2)**2/density_norm)/density_norm
1304 direct_error = maxval(abs(correction - expected))
1305 overlap_error = maxval(abs(correction_overlap - correction))
1306 sign_error = max(0.0_dp, correction(1, 1))
1307
1308 CALL qs_ot_density_secant_projected_hessian( &
1309 density_norm, response_work, density_overlap, response_overlap, calibrated, &
1310 calibrated_valid, secant_mode=2, secant_position=secant_position)
1311 expected(1, 1) = response_work
1312 expected(1, 2) = response_work/secant_position
1313 expected(2, 1) = expected(1, 2)
1314 expected(2, 2) = response_work/secant_position**2
1315 calibrated_error = maxval(abs(calibrated - expected))
1316 CALL qs_ot_density_secant_projected_hessian( &
1317 density_norm, response_work, density_overlap, response_overlap, rejected, &
1318 rejected_valid, secant_mode=2, secant_position=0.0_dp)
1319 rejected_error = maxval(abs(rejected))
1320
1321 unitary(:, :) = cmplx(0.0_dp, 0.0_dp, kind=dp)
1322 unitary(1, 1) = cmplx(sqrt(0.5_dp), 0.0_dp, kind=dp)
1323 unitary(2, 1) = cmplx(0.0_dp, sqrt(0.5_dp), kind=dp)
1324 unitary(1, 2) = cmplx(0.0_dp, sqrt(0.5_dp), kind=dp)
1325 unitary(2, 2) = cmplx(sqrt(0.5_dp), 0.0_dp, kind=dp)
1326 unitary(3, 3) = cmplx(1.0_dp, 0.0_dp, kind=dp)
1327 density_step = matmul(conjg(transpose(unitary)), matmul(density_step, unitary))
1328 hamiltonian_step = matmul(conjg(transpose(unitary)), matmul(hamiltonian_step, unitary))
1329 transformed_modes(:, :, 1) = &
1330 matmul(conjg(transpose(unitary)), matmul(density_modes(:, :, 1), unitary))
1331 transformed_modes(:, :, 2) = &
1332 matmul(conjg(transpose(unitary)), matmul(density_modes(:, :, 2), unitary))
1333 CALL qs_ot_density_secant_hessian( &
1334 density_step, hamiltonian_step, transformed_modes, correction_gauge, valid)
1335 gauge_error = maxval(abs(correction_gauge - correction))
1336
1337 IF (.NOT. valid .OR. .NOT. calibrated_valid .OR. rejected_valid .OR. &
1338 direct_error > 2.0e-14_dp .OR. overlap_error > 2.0e-14_dp .OR. &
1339 gauge_error > 2.0e-14_dp .OR. sign_error > 2.0e-14_dp .OR. &
1340 calibrated_error > 2.0e-14_dp .OR. rejected_error > 2.0e-14_dp) nfail = nfail + 1
1341 IF (mynode == 0) WRITE (io_unit, '(A,6(1X,ES13.6))') &
1342 "density secant direct/overlap/gauge/sign/calibrated/rejected errors:", &
1343 direct_error, overlap_error, gauge_error, sign_error, calibrated_error, rejected_error
1344
1345 END SUBROUTINE test_density_secant_hessian
1346
1347! **************************************************************************************************
1348!> \brief Check that the Hxc secant is formed in the joint direct sum of both spin channels.
1349!> \param mynode MPI rank
1350!> \param nfail accumulated number of failures
1351! **************************************************************************************************
1352 SUBROUTINE test_joint_spin_density_secant(mynode, nfail)
1353 INTEGER, INTENT(IN) :: mynode
1354 INTEGER, INTENT(INOUT) :: nfail
1355
1356 INTEGER, PARAMETER :: n = 2, nmode = 2, nspin = 2, &
1357 ntotal = n*nspin
1358
1359 COMPLEX(KIND=dp), DIMENSION(n, n, nmode, nspin) :: density_modes
1360 COMPLEX(KIND=dp), DIMENSION(n, n, nspin) :: density_step, hamiltonian_step
1361 COMPLEX(KIND=dp), DIMENSION(ntotal, ntotal) :: direct_density_step, &
1362 direct_hamiltonian_step
1363 COMPLEX(KIND=dp), DIMENSION(ntotal, ntotal, nmode) :: direct_density_modes
1364 INTEGER :: first, ispin, last, mode
1365 LOGICAL :: direct_valid, joint_valid, spin_valid
1366 REAL(KIND=dp) :: coupling_size, density_norm, &
1367 direct_error, response_work
1368 REAL(KIND=dp), DIMENSION(nmode) :: density_overlap, response_overlap
1369 REAL(KIND=dp), DIMENSION(nmode, nmode) :: direct_correction, joint_correction, &
1370 separate_correction, spin_correction
1371
1372 density_step(:, :, :) = cmplx(0.0_dp, 0.0_dp, kind=dp)
1373 density_step(1, 1, 1) = 0.31_dp
1374 density_step(2, 2, 1) = -0.31_dp
1375 density_step(1, 2, 1) = cmplx(0.08_dp, -0.04_dp, kind=dp)
1376 density_step(2, 1, 1) = conjg(density_step(1, 2, 1))
1377 density_step(1, 1, 2) = -0.19_dp
1378 density_step(2, 2, 2) = 0.19_dp
1379 density_step(1, 2, 2) = cmplx(-0.05_dp, 0.06_dp, kind=dp)
1380 density_step(2, 1, 2) = conjg(density_step(1, 2, 2))
1381
1382 hamiltonian_step(:, :, 1) = 0.67_dp*density_step(:, :, 1)
1383 hamiltonian_step(1, 2, 1) = hamiltonian_step(1, 2, 1) + cmplx(-0.03_dp, 0.09_dp, kind=dp)
1384 hamiltonian_step(2, 1, 1) = conjg(hamiltonian_step(1, 2, 1))
1385 hamiltonian_step(:, :, 2) = -0.43_dp*density_step(:, :, 2)
1386 hamiltonian_step(1, 2, 2) = hamiltonian_step(1, 2, 2) + cmplx(0.07_dp, 0.02_dp, kind=dp)
1387 hamiltonian_step(2, 1, 2) = conjg(hamiltonian_step(1, 2, 2))
1388
1389 density_modes(:, :, :, :) = cmplx(0.0_dp, 0.0_dp, kind=dp)
1390 density_modes(1, 1, 1, 1) = 0.21_dp
1391 density_modes(2, 2, 1, 1) = -0.21_dp
1392 density_modes(1, 2, 1, 1) = cmplx(0.04_dp, 0.03_dp, kind=dp)
1393 density_modes(2, 1, 1, 1) = conjg(density_modes(1, 2, 1, 1))
1394 density_modes(1, 1, 2, 1) = -0.13_dp
1395 density_modes(2, 2, 2, 1) = 0.13_dp
1396 density_modes(1, 2, 2, 1) = cmplx(-0.02_dp, 0.05_dp, kind=dp)
1397 density_modes(2, 1, 2, 1) = conjg(density_modes(1, 2, 2, 1))
1398 density_modes(1, 1, 1, 2) = -0.17_dp
1399 density_modes(2, 2, 1, 2) = 0.17_dp
1400 density_modes(1, 2, 1, 2) = cmplx(0.06_dp, -0.01_dp, kind=dp)
1401 density_modes(2, 1, 1, 2) = conjg(density_modes(1, 2, 1, 2))
1402 density_modes(1, 1, 2, 2) = 0.09_dp
1403 density_modes(2, 2, 2, 2) = -0.09_dp
1404 density_modes(1, 2, 2, 2) = cmplx(0.03_dp, 0.04_dp, kind=dp)
1405 density_modes(2, 1, 2, 2) = conjg(density_modes(1, 2, 2, 2))
1406
1407 direct_density_step(:, :) = cmplx(0.0_dp, 0.0_dp, kind=dp)
1408 direct_hamiltonian_step(:, :) = cmplx(0.0_dp, 0.0_dp, kind=dp)
1409 direct_density_modes(:, :, :) = cmplx(0.0_dp, 0.0_dp, kind=dp)
1410 density_norm = 0.0_dp
1411 response_work = 0.0_dp
1412 density_overlap(:) = 0.0_dp
1413 response_overlap(:) = 0.0_dp
1414 separate_correction(:, :) = 0.0_dp
1415 spin_valid = .true.
1416 DO ispin = 1, nspin
1417 first = (ispin - 1)*n + 1
1418 last = ispin*n
1419 direct_density_step(first:last, first:last) = density_step(:, :, ispin)
1420 direct_hamiltonian_step(first:last, first:last) = hamiltonian_step(:, :, ispin)
1421 direct_density_modes(first:last, first:last, :) = density_modes(:, :, :, ispin)
1422 density_norm = density_norm + &
1423 sum(real(conjg(density_step(:, :, ispin))* &
1424 density_step(:, :, ispin), kind=dp))
1425 response_work = response_work + &
1426 sum(real(conjg(density_step(:, :, ispin))* &
1427 hamiltonian_step(:, :, ispin), kind=dp))
1428 DO mode = 1, nmode
1429 density_overlap(mode) = density_overlap(mode) + &
1430 sum(real(conjg(density_step(:, :, ispin))* &
1431 density_modes(:, :, mode, ispin), kind=dp))
1432 response_overlap(mode) = response_overlap(mode) + &
1433 sum(real(conjg(hamiltonian_step(:, :, ispin))* &
1434 density_modes(:, :, mode, ispin), kind=dp))
1435 END DO
1436 CALL qs_ot_density_secant_hessian( &
1437 density_step(:, :, ispin), hamiltonian_step(:, :, ispin), &
1438 density_modes(:, :, :, ispin), spin_correction, direct_valid)
1439 spin_valid = spin_valid .AND. direct_valid
1440 separate_correction = separate_correction + spin_correction
1441 END DO
1442
1443 CALL qs_ot_density_secant_hessian( &
1444 direct_density_step, direct_hamiltonian_step, direct_density_modes, &
1445 direct_correction, direct_valid)
1446 CALL qs_ot_density_secant_projected_hessian( &
1447 density_norm, response_work, density_overlap, response_overlap, &
1448 joint_correction, joint_valid)
1449 direct_error = maxval(abs(joint_correction - direct_correction))
1450 coupling_size = maxval(abs(joint_correction - separate_correction))
1451
1452 IF (.NOT. direct_valid .OR. .NOT. joint_valid .OR. .NOT. spin_valid .OR. &
1453 direct_error > 3.0e-14_dp .OR. coupling_size < 1.0e-5_dp) nfail = nfail + 1
1454 IF (mynode == 0) WRITE (io_unit, '(A,2(1X,ES13.6))') &
1455 "joint-spin density secant direct error/coupling:", direct_error, coupling_size
1456
1457 END SUBROUTINE test_joint_spin_density_secant
1458
1459! **************************************************************************************************
1460!> \brief Check the coupled finite-chart density tangent against a central finite difference.
1461!> \param mynode MPI rank
1462!> \param nfail accumulated number of failures
1463! **************************************************************************************************
1464 SUBROUTINE test_density_tangent(mynode, nfail)
1465 INTEGER, INTENT(IN) :: mynode
1466 INTEGER, INTENT(INOUT) :: nfail
1467
1468 INTEGER, PARAMETER :: nbands = 3, nrotation = 6
1469 REAL(KIND=dp), PARAMETER :: fd_step = 2.0e-6_dp, &
1470 kpoint_weight = 0.625_dp
1471
1472 COMPLEX(KIND=dp), DIMENSION(nbands, nbands) :: density_minus, density_plus, density_tangent, &
1473 direction_generator, generator, rotation, rotation_minus, rotation_plus, weighted_rotation
1474 INTEGER :: i, j, r
1475 REAL(KIND=dp) :: error, trace_error
1476 REAL(KIND=dp), DIMENSION(nbands) :: occupation, weighted_occupation_step
1477 REAL(KIND=dp), DIMENSION(nrotation) :: rotation_step
1478
1479 generator(:, :) = cmplx(0.0_dp, 0.0_dp, kind=dp)
1480 generator(1, 2) = cmplx(0.17_dp, -0.09_dp, kind=dp)
1481 generator(2, 1) = -conjg(generator(1, 2))
1482 generator(1, 3) = cmplx(-0.06_dp, 0.11_dp, kind=dp)
1483 generator(3, 1) = -conjg(generator(1, 3))
1484 generator(2, 3) = cmplx(0.08_dp, 0.04_dp, kind=dp)
1485 generator(3, 2) = -conjg(generator(2, 3))
1486 occupation(:) = [1.72_dp, 0.83_dp, 0.31_dp]
1487 rotation_step(:) = [0.21_dp, -0.14_dp, -0.08_dp, 0.17_dp, 0.11_dp, 0.06_dp]
1488 weighted_occupation_step(:) = [-0.037_dp, 0.052_dp, -0.015_dp]
1489
1490 CALL qs_ot_density_tangent( &
1491 generator, occupation, kpoint_weight, rotation_step, weighted_occupation_step, &
1492 density_tangent, difference_step=fd_step)
1493
1494 direction_generator(:, :) = cmplx(0.0_dp, 0.0_dp, kind=dp)
1495 r = 0
1496 DO i = 1, nbands - 1
1497 DO j = i + 1, nbands
1498 r = r + 1
1499 direction_generator(i, j) = cmplx(rotation_step(r), 0.0_dp, kind=dp)
1500 direction_generator(j, i) = -direction_generator(i, j)
1501 r = r + 1
1502 direction_generator(i, j) = direction_generator(i, j) + &
1503 cmplx(0.0_dp, rotation_step(r), kind=dp)
1504 direction_generator(j, i) = direction_generator(j, i) + &
1505 cmplx(0.0_dp, rotation_step(r), kind=dp)
1506 END DO
1507 END DO
1508 rotation = dense_antihermitian_exp(generator)
1509 rotation_plus = dense_antihermitian_exp(generator + fd_step*direction_generator)
1510 rotation_minus = dense_antihermitian_exp(generator - fd_step*direction_generator)
1511 weighted_rotation(:, :) = rotation_plus
1512 DO j = 1, nbands
1513 weighted_rotation(:, j) = &
1514 (kpoint_weight*occupation(j) + fd_step*weighted_occupation_step(j))* &
1515 weighted_rotation(:, j)
1516 END DO
1517 density_plus(:, :) = matmul(weighted_rotation, conjg(transpose(rotation_plus)))
1518 weighted_rotation(:, :) = rotation_minus
1519 DO j = 1, nbands
1520 weighted_rotation(:, j) = &
1521 (kpoint_weight*occupation(j) - fd_step*weighted_occupation_step(j))* &
1522 weighted_rotation(:, j)
1523 END DO
1524 density_minus(:, :) = matmul(weighted_rotation, conjg(transpose(rotation_minus)))
1525 density_plus(:, :) = matmul(conjg(transpose(rotation)), &
1526 matmul((density_plus - density_minus)/(2.0_dp*fd_step), rotation))
1527 error = maxval(abs(density_tangent - density_plus))
1528 trace_error = abs(real(sum([(density_tangent(i, i), i=1, nbands)]), kind=dp) - &
1529 sum(weighted_occupation_step))
1530
1531 IF (error > 3.0e-9_dp .OR. trace_error > 3.0e-10_dp) nfail = nfail + 1
1532 IF (mynode == 0) WRITE (io_unit, '(A,2(1X,ES13.6))') &
1533 "finite complex density tangent errors:", error, trace_error
1534
1535 END SUBROUTINE test_density_tangent
1536
1537! **************************************************************************************************
1538!> \brief Check density-secants when accepted endpoint orbital subspaces differ.
1539!> \param mynode MPI rank
1540!> \param nfail accumulated number of failures
1541! **************************************************************************************************
1542 SUBROUTINE test_density_secant_moving_subspace(mynode, nfail)
1543 INTEGER, INTENT(IN) :: mynode
1544 INTEGER, INTENT(INOUT) :: nfail
1545
1546 INTEGER, PARAMETER :: nao = 4, nbands = 3, nmode = 2
1547 REAL(KIND=dp), PARAMETER :: weight = 0.625_dp
1548
1549 COMPLEX(KIND=dp), DIMENSION(nao, nao) :: density_step, hamiltonian_step
1550 COMPLEX(KIND=dp), DIMENSION(nao, nao, nmode) :: density_modes_ao
1551 COMPLEX(KIND=dp), DIMENSION(nao, nbands) :: c0, c1
1552 COMPLEX(KIND=dp), DIMENSION(nbands, nbands) :: h0, h1, overlap
1553 COMPLEX(KIND=dp), DIMENSION(nbands, nbands, nmode) :: density_modes
1554 INTEGER :: i
1555 LOGICAL :: valid
1556 REAL(KIND=dp) :: density_norm, density_norm_ref, error, &
1557 response_work, response_work_ref
1558 REAL(KIND=dp), DIMENSION(nbands) :: occupation0, occupation1
1559 REAL(KIND=dp), DIMENSION(nmode) :: density_overlap, density_overlap_ref, &
1560 response_overlap, response_overlap_ref
1561
1562 c0(:, :) = cmplx(0.0_dp, 0.0_dp, kind=dp)
1563 c0(1, 1) = 1.0_dp
1564 c0(2, 2) = 1.0_dp
1565 c0(3, 3) = 1.0_dp
1566 c1(:, :) = c0
1567 c1(1, 1) = cos(0.31_dp)
1568 c1(4, 1) = sin(0.31_dp)
1569 c1(2, 2) = cos(0.23_dp)
1570 c1(3, 2) = cmplx(0.0_dp, sin(0.23_dp), kind=dp)
1571 c1(2, 3) = cmplx(0.0_dp, sin(0.23_dp), kind=dp)
1572 c1(3, 3) = cos(0.23_dp)
1573 occupation0(:) = [1.73_dp, 0.82_dp, 0.82_dp]
1574 occupation1(:) = [1.61_dp, 0.91_dp, 0.76_dp]
1575
1576 hamiltonian_step(:, :) = cmplx(0.0_dp, 0.0_dp, kind=dp)
1577 hamiltonian_step(1, 1) = 0.21_dp
1578 hamiltonian_step(2, 2) = -0.17_dp
1579 hamiltonian_step(3, 3) = 0.09_dp
1580 hamiltonian_step(4, 4) = -0.04_dp
1581 hamiltonian_step(1, 4) = cmplx(0.07_dp, -0.03_dp, kind=dp)
1582 hamiltonian_step(4, 1) = conjg(hamiltonian_step(1, 4))
1583 hamiltonian_step(2, 3) = cmplx(-0.05_dp, 0.08_dp, kind=dp)
1584 hamiltonian_step(3, 2) = conjg(hamiltonian_step(2, 3))
1585 density_modes(:, :, :) = cmplx(0.0_dp, 0.0_dp, kind=dp)
1586 density_modes(1, 1, 1) = 0.13_dp
1587 density_modes(2, 2, 1) = -0.08_dp
1588 density_modes(3, 3, 1) = -0.05_dp
1589 density_modes(1, 2, 2) = cmplx(0.04_dp, -0.06_dp, kind=dp)
1590 density_modes(2, 1, 2) = conjg(density_modes(1, 2, 2))
1591
1592 overlap(:, :) = matmul(conjg(transpose(c0)), c1)
1593 ! Keep the reference path independent of optimized complex MATMUL expressions.
1594 CALL reference_project_operator(c0, hamiltonian_step, h0)
1595 CALL reference_project_operator(c1, hamiltonian_step, h1)
1596 CALL qs_ot_density_secant_orbital_overlaps( &
1597 overlap, occupation0, occupation1, h0, h1, density_modes, weight, density_norm, &
1598 response_work, density_overlap, response_overlap, valid)
1599
1600 density_step(:, :) = cmplx(0.0_dp, 0.0_dp, kind=dp)
1601 DO i = 1, nbands
1602 density_step = density_step + weight*occupation1(i)* &
1603 spread(c1(:, i), dim=2, ncopies=nao)* &
1604 spread(conjg(c1(:, i)), dim=1, ncopies=nao) - &
1605 weight*occupation0(i)* &
1606 spread(c0(:, i), dim=2, ncopies=nao)* &
1607 spread(conjg(c0(:, i)), dim=1, ncopies=nao)
1608 END DO
1609 DO i = 1, nmode
1610 CALL reference_transform_density_mode(c1, density_modes(:, :, i), density_modes_ao(:, :, i))
1611 density_overlap_ref(i) = &
1612 reference_frobenius_product_real(density_step, density_modes_ao(:, :, i))/weight
1613 response_overlap_ref(i) = &
1614 reference_frobenius_product_real(hamiltonian_step, density_modes_ao(:, :, i))
1615 END DO
1616 density_norm_ref = reference_frobenius_product_real(density_step, density_step)/weight
1617 response_work_ref = reference_frobenius_product_real(density_step, hamiltonian_step)
1618 error = max(abs(density_norm - density_norm_ref), abs(response_work - response_work_ref))
1619 error = max(error, maxval(abs(density_overlap - density_overlap_ref)))
1620 error = max(error, maxval(abs(response_overlap - response_overlap_ref)))
1621
1622 IF (.NOT. valid .OR. error > 3.0e-13_dp) nfail = nfail + 1
1623 IF (mynode == 0) WRITE (io_unit, '(A,1X,ES13.6)') &
1624 "moving-subspace density secant error:", error
1625
1626 END SUBROUTINE test_density_secant_moving_subspace
1627
1628! **************************************************************************************************
1629!> \brief Reference projection C^H H C using scalar operations.
1630!> \param orbitals orbital coefficient matrix C
1631!> \param hamiltonian operator matrix H
1632!> \param projected projected operator
1633! **************************************************************************************************
1634 SUBROUTINE reference_project_operator(orbitals, hamiltonian, projected)
1635 COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(IN) :: orbitals, hamiltonian
1636 COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: projected
1637
1638 COMPLEX(KIND=dp) :: value
1639 INTEGER :: i, j, mu, nu
1640
1641 DO j = 1, SIZE(orbitals, 2)
1642 DO i = 1, SIZE(orbitals, 2)
1643 value = cmplx(0.0_dp, 0.0_dp, kind=dp)
1644 DO nu = 1, SIZE(orbitals, 1)
1645 DO mu = 1, SIZE(orbitals, 1)
1646 value = value + conjg(orbitals(mu, i))*hamiltonian(mu, nu)*orbitals(nu, j)
1647 END DO
1648 END DO
1649 projected(i, j) = value
1650 END DO
1651 END DO
1652
1653 END SUBROUTINE reference_project_operator
1654
1655! **************************************************************************************************
1656!> \brief Reference transformation C M C^H using scalar operations.
1657!> \param orbitals orbital coefficient matrix C
1658!> \param mode density mode M
1659!> \param transformed transformed density mode
1660! **************************************************************************************************
1661 SUBROUTINE reference_transform_density_mode(orbitals, mode, transformed)
1662 COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(IN) :: orbitals, mode
1663 COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: transformed
1664
1665 COMPLEX(KIND=dp) :: value
1666 INTEGER :: i, j, m, n
1667
1668 DO j = 1, SIZE(orbitals, 1)
1669 DO i = 1, SIZE(orbitals, 1)
1670 value = cmplx(0.0_dp, 0.0_dp, kind=dp)
1671 DO n = 1, SIZE(orbitals, 2)
1672 DO m = 1, SIZE(orbitals, 2)
1673 value = value + orbitals(i, m)*mode(m, n)*conjg(orbitals(j, n))
1674 END DO
1675 END DO
1676 transformed(i, j) = value
1677 END DO
1678 END DO
1679
1680 END SUBROUTINE reference_transform_density_mode
1681
1682! **************************************************************************************************
1683!> \brief Real part of the complex Frobenius inner product.
1684!> \param matrix_a first matrix
1685!> \param matrix_b second matrix
1686!> \return real Frobenius product
1687! **************************************************************************************************
1688 FUNCTION reference_frobenius_product_real(matrix_a, matrix_b) RESULT(value)
1689 COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(IN) :: matrix_a, matrix_b
1690 REAL(KIND=dp) :: value
1691
1692 value = sum(real(matrix_a, kind=dp)*real(matrix_b, kind=dp) + &
1693 aimag(matrix_a)*aimag(matrix_b))
1694
1695 END FUNCTION reference_frobenius_product_real
1696
1697! **************************************************************************************************
1698!> \brief fixed-occupation band energy for a dense finite rotation
1699!> \param rotation unitary REF rotation
1700!> \param hamiltonian fixed Hermitian Hamiltonian in the REF basis
1701!> \param occupation band occupations
1702!> \param kpoint_weight irreducible-k-point weight
1703!> \return weighted band energy
1704! **************************************************************************************************
1705 FUNCTION finite_rotation_energy(rotation, hamiltonian, occupation, kpoint_weight) RESULT(energy)
1706 COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(IN) :: rotation, hamiltonian
1707 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: occupation
1708 REAL(KIND=dp), INTENT(IN) :: kpoint_weight
1709 REAL(KIND=dp) :: energy
1710
1711 INTEGER :: i
1712
1713 energy = 0.0_dp
1714 DO i = 1, SIZE(occupation)
1715 energy = energy + kpoint_weight*occupation(i)* &
1716 REAL(DOT_PRODUCT(rotation(:, i), &
1717 MATMUL(hamiltonian, rotation(:, i))), KIND=dp)
1718 END DO
1719
1720 END FUNCTION finite_rotation_energy
1721
1722! **************************************************************************************************
1723!> \brief Solve fixed-N Fermi occupations by bisection.
1724!> \param energies band energies
1725!> \param weights integration weights
1726!> \param target_electrons requested electron number
1727!> \param temperature electronic temperature in energy units
1728!> \param maxocc maximum band occupation
1729!> \param occupation resulting occupations
1730!> \param mu resulting chemical potential
1731! **************************************************************************************************
1732 SUBROUTINE fixed_n_fermi_occupations(energies, weights, target_electrons, temperature, &
1733 maxocc, occupation, mu)
1734 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: energies, weights
1735 REAL(KIND=dp), INTENT(IN) :: target_electrons, temperature, maxocc
1736 REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: occupation
1737 REAL(KIND=dp), INTENT(OUT) :: mu
1738
1739 INTEGER :: i, iteration
1740 REAL(KIND=dp) :: high, low, number
1741
1742 low = minval(energies) - 80.0_dp*temperature
1743 high = maxval(energies) + 80.0_dp*temperature
1744 DO iteration = 1, 160
1745 mu = 0.5_dp*(low + high)
1746 number = 0.0_dp
1747 DO i = 1, SIZE(energies)
1748 occupation(i) = fermi_occupation(energies(i), mu, temperature, maxocc)
1749 number = number + weights(i)*occupation(i)
1750 END DO
1751 IF (number < target_electrons) THEN
1752 low = mu
1753 ELSE
1754 high = mu
1755 END IF
1756 END DO
1757 mu = 0.5_dp*(low + high)
1758 DO i = 1, SIZE(energies)
1759 occupation(i) = fermi_occupation(energies(i), mu, temperature, maxocc)
1760 END DO
1761
1762 END SUBROUTINE fixed_n_fermi_occupations
1763
1764! **************************************************************************************************
1765!> \brief Numerically stable Fermi occupation.
1766!> \param energy band energy
1767!> \param mu chemical potential
1768!> \param temperature electronic temperature in energy units
1769!> \param maxocc maximum band occupation
1770!> \return Fermi occupation
1771! **************************************************************************************************
1772 PURE FUNCTION fermi_occupation(energy, mu, temperature, maxocc) RESULT(occupation)
1773 REAL(KIND=dp), INTENT(IN) :: energy, mu, temperature, maxocc
1774 REAL(KIND=dp) :: occupation
1775
1776 REAL(KIND=dp) :: x
1777
1778 x = (energy - mu)/temperature
1779 IF (x > 40.0_dp) THEN
1780 occupation = maxocc*exp(-x)
1781 ELSE IF (x < -40.0_dp) THEN
1782 occupation = maxocc*(1.0_dp - exp(x))
1783 ELSE
1784 occupation = maxocc/(1.0_dp + exp(x))
1785 END IF
1786
1787 END FUNCTION fermi_occupation
1788
1789! **************************************************************************************************
1790!> \brief Fixed-H fixed-N Mermin free energy used by the finite-difference test.
1791!> \param rayleigh Hamiltonian expectation values
1792!> \param energies auxiliary band energies
1793!> \param weights integration weights
1794!> \param target_electrons requested electron number
1795!> \param temperature electronic temperature in energy units
1796!> \param maxocc maximum band occupation
1797!> \return fixed-N Mermin free energy
1798! **************************************************************************************************
1799 FUNCTION fixed_n_mermin_value(rayleigh, energies, weights, target_electrons, &
1800 temperature, maxocc) RESULT(value)
1801 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: rayleigh, energies, weights
1802 REAL(KIND=dp), INTENT(IN) :: target_electrons, temperature, maxocc
1803 REAL(KIND=dp) :: value
1804
1805 INTEGER :: i
1806 REAL(KIND=dp) :: mu, p
1807 REAL(KIND=dp), DIMENSION(SIZE(energies)) :: occupation
1808
1809 CALL fixed_n_fermi_occupations(energies, weights, target_electrons, temperature, &
1810 maxocc, occupation, mu)
1811 value = 0.0_dp
1812 DO i = 1, SIZE(energies)
1813 p = max(epsilon(1.0_dp), min(1.0_dp - epsilon(1.0_dp), occupation(i)/maxocc))
1814 value = value + weights(i)*(occupation(i)*rayleigh(i) + &
1815 temperature*maxocc*(p*log(p) + (1.0_dp - p)*log(1.0_dp - p)))
1816 END DO
1817
1818 END FUNCTION fixed_n_mermin_value
1819
1820! **************************************************************************************************
1821!> \brief Canonical weighted occupation projector used by the finite-difference test.
1822!> \param hamiltonian Hermitian Hamiltonian
1823!> \param weights integration weights
1824!> \param target_electrons requested electron number
1825!> \param temperature electronic temperature in energy units
1826!> \param maxocc maximum band occupation
1827!> \param projector weighted spectral occupation projector
1828! **************************************************************************************************
1829 SUBROUTINE fixed_n_spectral_projector(hamiltonian, weights, target_electrons, temperature, &
1830 maxocc, projector)
1831 COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(IN) :: hamiltonian
1832 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: weights
1833 REAL(KIND=dp), INTENT(IN) :: target_electrons, temperature, maxocc
1834 COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: projector
1835
1836 COMPLEX(KIND=dp), DIMENSION(SIZE(hamiltonian, 1), & SIZE(hamiltonian, 1)) :: vectors, weighted_vectors
1837 INTEGER :: i
1838 REAL(KIND=dp) :: mu
1839 REAL(KIND=dp), DIMENSION(SIZE(hamiltonian, 1)) :: eigenvalues, occupation
1840
1841 CALL diag_complex(hamiltonian, vectors, eigenvalues)
1842 CALL fixed_n_fermi_occupations(eigenvalues, weights, target_electrons, temperature, &
1843 maxocc, occupation, mu)
1844 weighted_vectors(:, :) = vectors(:, :)
1845 DO i = 1, SIZE(eigenvalues)
1846 weighted_vectors(:, i) = weights(i)*occupation(i)*weighted_vectors(:, i)
1847 END DO
1848 projector(:, :) = matmul(weighted_vectors, conjg(transpose(vectors)))
1849
1850 END SUBROUTINE fixed_n_spectral_projector
1851
1852! **************************************************************************************************
1853!> \brief Check the real finite rotation pullback used by Gamma-point REF OT.
1854!> \param para_env parallel environment
1855!> \param nfail accumulated number of failures
1856! **************************************************************************************************
1857 SUBROUTINE test_real_rotation_frechet(para_env, nfail)
1858 TYPE(mp_para_env_type), POINTER :: para_env
1859 INTEGER, INTENT(INOUT) :: nfail
1860
1861 INTEGER, PARAMETER :: m = 3
1862
1863 INTEGER :: i
1864 INTEGER, DIMENSION(:), POINTER :: col_dist, col_size, row_dist, row_size
1865 INTEGER, DIMENSION(:, :), POINTER :: pgrid
1866 LOGICAL :: found
1867 REAL(kind=dp) :: error, fd_slope, hstep, predicted
1868 REAL(kind=dp), DIMENSION(:, :), POINTER :: block
1869 REAL(kind=dp), DIMENSION(m), TARGET :: rotation_evals
1870 REAL(kind=dp), DIMENSION(m, m) :: dedu, direction, generator, gradient, &
1871 rotation_minus, rotation_plus
1872 TYPE(cp_blacs_env_type), POINTER :: blacs_env
1873 TYPE(dbcsr_distribution_type) :: dist
1874 TYPE(dbcsr_type), DIMENSION(7), TARGET :: matrices
1875 TYPE(mp_para_env_type), POINTER :: local_para_env
1876 TYPE(qs_ot_type) :: rotation_env
1877
1878 NULLIFY (blacs_env, local_para_env)
1879 generator(:, :) = 0.0_dp
1880 generator(1, 2) = 0.23_dp
1881 generator(2, 1) = -generator(1, 2)
1882 generator(1, 3) = -0.11_dp
1883 generator(3, 1) = -generator(1, 3)
1884 generator(2, 3) = 0.17_dp
1885 generator(3, 2) = -generator(2, 3)
1886
1887 direction(:, :) = 0.0_dp
1888 direction(1, 2) = -0.09_dp
1889 direction(2, 1) = -direction(1, 2)
1890 direction(1, 3) = 0.14_dp
1891 direction(3, 1) = -direction(1, 3)
1892 direction(2, 3) = 0.06_dp
1893 direction(3, 2) = -direction(2, 3)
1894
1895 dedu = reshape([0.31_dp, -0.27_dp, 0.08_dp, &
1896 0.19_dp, 0.42_dp, -0.16_dp, &
1897 -0.07_dp, 0.21_dp, 0.35_dp], [m, m])
1898
1899 ALLOCATE (local_para_env)
1900 CALL local_para_env%from_dup(mp_comm_self)
1901 CALL cp_blacs_env_create(blacs_env=blacs_env, para_env=local_para_env)
1902 ALLOCATE (pgrid(0:0, 0:0), row_dist(1), col_dist(1), row_size(1), col_size(1))
1903 pgrid(:, :) = 0
1904 row_dist(:) = 0
1905 col_dist(:) = 0
1906 row_size(:) = m
1907 col_size(:) = m
1908 CALL dbcsr_distribution_new(dist, group=mp_comm_self%get_handle(), pgrid=pgrid, &
1909 row_dist=row_dist, col_dist=col_dist)
1910 DO i = 1, SIZE(matrices)
1911 CALL create_one_block(matrices(i), 'real_rotation', dist, row_size, col_size)
1912 END DO
1913 CALL dbcsr_put_block(matrices(2), 1, 1, generator)
1914 CALL dbcsr_put_block(matrices(3), 1, 1, dedu)
1915
1916 rotation_env%para_env => local_para_env
1917 rotation_env%blacs_env => blacs_env
1918 rotation_env%rot_mat_u => matrices(1)
1919 rotation_env%rot_mat_x => matrices(2)
1920 rotation_env%rot_mat_dedu => matrices(3)
1921 rotation_env%rot_mat_gx => matrices(4)
1922 rotation_env%rot_mat_evec_re => matrices(5)
1923 rotation_env%rot_mat_evec_im => matrices(6)
1924 rotation_env%matrix_buf1 => matrices(7)
1925 rotation_env%rot_mat_evals => rotation_evals
1926
1927 CALL qs_ot_generate_rotation(rotation_env)
1928 CALL qs_ot_rot_mat_derivative(rotation_env)
1929 gradient(:, :) = 0.0_dp
1930 CALL dbcsr_get_block_p(matrices(4), 1, 1, block, found)
1931 IF (found) gradient = block
1932
1933 hstep = 1.0e-6_dp
1934 rotation_plus = real(dense_antihermitian_exp( &
1935 cmplx(generator + hstep*direction, 0.0_dp, kind=dp)), kind=dp)
1936 rotation_minus = real(dense_antihermitian_exp( &
1937 cmplx(generator - hstep*direction, 0.0_dp, kind=dp)), kind=dp)
1938 fd_slope = sum(dedu*(rotation_plus - rotation_minus))/(2.0_dp*hstep)
1939 predicted = 0.5_dp*sum(gradient*direction)
1940 error = abs(predicted - fd_slope)
1941 IF (error > 5.0e-8_dp) nfail = nfail + 1
1942 IF (para_env%is_source()) THEN
1943 WRITE (io_unit, '(A,3(1X,ES13.6))') 'real finite rotation slope fd/predicted/error', &
1944 fd_slope, predicted, error
1945 END IF
1946
1947 DO i = 1, SIZE(matrices)
1948 CALL dbcsr_release(matrices(i))
1949 END DO
1950 CALL dbcsr_distribution_release(dist)
1951 DEALLOCATE (pgrid, row_dist, col_dist, row_size, col_size)
1952 CALL cp_blacs_env_release(blacs_env)
1953 CALL mp_para_env_release(local_para_env)
1954
1955 END SUBROUTINE test_real_rotation_frechet
1956
1957! **************************************************************************************************
1958!> \brief Check the combined finite complex REF and rotation derivatives.
1959!> \param para_env parallel environment
1960!> \param nfail accumulated number of failures
1961! **************************************************************************************************
1962 SUBROUTINE test_complex_ref_rotation(para_env, nfail)
1963 TYPE(mp_para_env_type), POINTER :: para_env
1964 INTEGER, INTENT(INOUT) :: nfail
1965
1966 COMPLEX(KIND=dp), DIMENSION(k, k) :: direction_a, generator, gradient_a, &
1967 inverse_sqrt, rotation
1968 COMPLEX(KIND=dp), DIMENSION(n, k) :: c_dense, c_reference, direction_x, &
1969 gradient_x, hc_dense, q_dense, x_dense
1970 COMPLEX(KIND=dp), DIMENSION(n, n) :: hamiltonian, overlap
1971 INTEGER :: i, j
1972 INTEGER, DIMENSION(:), POINTER :: col_dist, col_k, col_n, row_dist, row_k, &
1973 row_n
1974 INTEGER, DIMENSION(:, :), POINTER :: pgrid
1975 REAL(KIND=dp) :: error_c, error_ref, error_rot, fd_ref, &
1976 fd_rot, hstep, predicted_ref, &
1977 predicted_rot
1978 REAL(KIND=dp), DIMENSION(k) :: occupation
1979 REAL(KIND=dp), DIMENSION(k), TARGET :: rotation_evals
1980 TYPE(cp_blacs_env_type), POINTER :: blacs_env
1981 TYPE(dbcsr_distribution_type) :: dist_kk, dist_nk, dist_nn
1982 TYPE(dbcsr_type), DIMENSION(10), TARGET :: rot_matrix
1983 TYPE(dbcsr_type), DIMENSION(2), TARGET :: nk_buffer
1984 TYPE(dbcsr_type), DIMENSION(3), TARGET :: nonsym_buffer
1985 TYPE(dbcsr_type), DIMENSION(4), TARGET :: sym_buffer
1986 TYPE(dbcsr_type), POINTER :: c_im_p, c_re_p, hc_im_p, hc_re_p, &
1987 s_im_p, s_re_p
1988 TYPE(dbcsr_type), TARGET :: c_im, c_re, f_im, f_re, g_im, g_re, &
1989 hc_im, hc_re, s_im, s_re, sx_im, &
1990 sx_re, x_im, x_re
1991 TYPE(mp_para_env_type), POINTER :: local_para_env
1992 TYPE(qs_ot_type) :: rotation_env
1993
1994 NULLIFY (blacs_env, c_im_p, c_re_p, hc_im_p, hc_re_p, local_para_env, &
1995 s_im_p, s_re_p)
1996
1997 x_dense(:, 1) = [cmplx(1.05_dp, 0.12_dp, kind=dp), cmplx(0.18_dp, -0.27_dp, kind=dp), &
1998 cmplx(-0.14_dp, 0.16_dp, kind=dp)]
1999 x_dense(:, 2) = [cmplx(0.26_dp, 0.31_dp, kind=dp), cmplx(0.94_dp, -0.17_dp, kind=dp), &
2000 cmplx(0.21_dp, 0.13_dp, kind=dp)]
2001 direction_x(:, 1) = [cmplx(0.11_dp, -0.07_dp, kind=dp), cmplx(-0.16_dp, 0.09_dp, kind=dp), &
2002 cmplx(0.05_dp, 0.14_dp, kind=dp)]
2003 direction_x(:, 2) = [cmplx(-0.08_dp, 0.15_dp, kind=dp), cmplx(0.13_dp, 0.04_dp, kind=dp), &
2004 cmplx(-0.10_dp, -0.06_dp, kind=dp)]
2005
2006 generator(:, :) = cmplx(0.0_dp, 0.0_dp, kind=dp)
2007 generator(1, 1) = cmplx(0.0_dp, 0.19_dp, kind=dp)
2008 generator(2, 2) = cmplx(0.0_dp, -0.08_dp, kind=dp)
2009 generator(1, 2) = cmplx(0.23_dp, -0.11_dp, kind=dp)
2010 generator(2, 1) = -conjg(generator(1, 2))
2011 direction_a(:, :) = cmplx(0.0_dp, 0.0_dp, kind=dp)
2012 direction_a(1, 1) = cmplx(0.0_dp, -0.07_dp, kind=dp)
2013 direction_a(2, 2) = cmplx(0.0_dp, 0.05_dp, kind=dp)
2014 direction_a(1, 2) = cmplx(-0.09_dp, 0.06_dp, kind=dp)
2015 direction_a(2, 1) = -conjg(direction_a(1, 2))
2016
2017 hamiltonian(:, :) = cmplx(0.0_dp, 0.0_dp, kind=dp)
2018 hamiltonian(1, 1) = cmplx(0.65_dp, 0.0_dp, kind=dp)
2019 hamiltonian(2, 2) = cmplx(-0.35_dp, 0.0_dp, kind=dp)
2020 hamiltonian(3, 3) = cmplx(1.15_dp, 0.0_dp, kind=dp)
2021 hamiltonian(1, 2) = cmplx(0.18_dp, 0.12_dp, kind=dp)
2022 hamiltonian(1, 3) = cmplx(-0.09_dp, 0.07_dp, kind=dp)
2023 hamiltonian(2, 3) = cmplx(0.14_dp, -0.17_dp, kind=dp)
2024 DO i = 1, n
2025 DO j = i + 1, n
2026 hamiltonian(j, i) = conjg(hamiltonian(i, j))
2027 END DO
2028 END DO
2029 overlap(:, :) = cmplx(0.0_dp, 0.0_dp, kind=dp)
2030 DO i = 1, n
2031 overlap(i, i) = cmplx(1.0_dp, 0.0_dp, kind=dp)
2032 END DO
2033 occupation = [1.0_dp, 0.37_dp]
2034
2035 CALL dense_inverse_sqrt(matmul(conjg(transpose(x_dense)), x_dense), inverse_sqrt)
2036 q_dense = matmul(x_dense, inverse_sqrt)
2037 rotation = dense_antihermitian_exp(generator)
2038 c_reference = matmul(q_dense, rotation)
2039
2040 ALLOCATE (local_para_env)
2041 CALL local_para_env%from_dup(mp_comm_self)
2042 CALL cp_blacs_env_create(blacs_env=blacs_env, para_env=local_para_env)
2043
2044 ALLOCATE (pgrid(0:0, 0:0), row_dist(1), col_dist(1), row_n(1), row_k(1), &
2045 col_n(1), col_k(1))
2046 pgrid(:, :) = 0
2047 row_dist(:) = 0
2048 col_dist(:) = 0
2049 row_n(:) = n
2050 col_n(:) = n
2051 row_k(:) = k
2052 col_k(:) = k
2053 CALL dbcsr_distribution_new(dist_nk, group=mp_comm_self%get_handle(), pgrid=pgrid, &
2054 row_dist=row_dist, col_dist=col_dist)
2055 CALL dbcsr_distribution_new(dist_nn, group=mp_comm_self%get_handle(), pgrid=pgrid, &
2056 row_dist=row_dist, col_dist=col_dist)
2057 CALL dbcsr_distribution_new(dist_kk, group=mp_comm_self%get_handle(), pgrid=pgrid, &
2058 row_dist=row_dist, col_dist=col_dist)
2059
2060 CALL create_one_block(x_re, 'rotation_ref_x_re', dist_nk, row_n, col_k)
2061 CALL create_one_block(x_im, 'rotation_ref_x_im', dist_nk, row_n, col_k)
2062 CALL create_one_block(sx_re, 'rotation_ref_sx_re', dist_nk, row_n, col_k)
2063 CALL create_one_block(sx_im, 'rotation_ref_sx_im', dist_nk, row_n, col_k)
2064 CALL create_one_block(c_re, 'rotation_ref_c_re', dist_nk, row_n, col_k)
2065 CALL create_one_block(c_im, 'rotation_ref_c_im', dist_nk, row_n, col_k)
2066 CALL create_one_block(hc_re, 'rotation_ref_hc_re', dist_nk, row_n, col_k)
2067 CALL create_one_block(hc_im, 'rotation_ref_hc_im', dist_nk, row_n, col_k)
2068 CALL create_one_block(g_re, 'rotation_ref_g_re', dist_nk, row_n, col_k)
2069 CALL create_one_block(g_im, 'rotation_ref_g_im', dist_nk, row_n, col_k)
2070 CALL create_one_block(s_re, 'rotation_ref_s_re', dist_nn, row_n, col_n)
2071 CALL create_one_block(s_im, 'rotation_ref_s_im', dist_nn, row_n, col_n)
2072 CALL create_one_block(f_re, 'rotation_ref_f_re', dist_kk, row_k, col_k)
2073 CALL create_one_block(f_im, 'rotation_ref_f_im', dist_kk, row_k, col_k)
2074 DO i = 1, SIZE(sym_buffer)
2075 CALL create_one_block(sym_buffer(i), 'rotation_ref_sym', dist_kk, row_k, col_k)
2076 END DO
2077 DO i = 1, SIZE(nonsym_buffer)
2078 CALL create_one_block(nonsym_buffer(i), 'rotation_ref_nonsym', dist_kk, row_k, col_k)
2079 END DO
2080 DO i = 1, SIZE(nk_buffer)
2081 CALL create_one_block(nk_buffer(i), 'rotation_ref_nk', dist_nk, row_n, col_k)
2082 END DO
2083 DO i = 1, SIZE(rot_matrix)
2084 CALL create_one_block(rot_matrix(i), 'rotation_ref_rot', dist_kk, row_k, col_k)
2085 END DO
2086
2087 CALL put_complex_pair(x_re, x_im, x_dense)
2088 CALL put_complex_pair(s_re, s_im, overlap)
2089 CALL put_complex_pair(rot_matrix(3), rot_matrix(4), generator)
2090
2091 rotation_env%has_complex_kpoint_state = .true.
2092 rotation_env%settings%do_rotation = .true.
2093 rotation_env%settings%eps_irac_filter_matrix = 0.0_dp
2094 rotation_env%settings%ot_method = 'LBFG'
2095 rotation_env%para_env => local_para_env
2096 rotation_env%blacs_env => blacs_env
2097 rotation_env%matrix_x => x_re
2098 rotation_env%matrix_x_im => x_im
2099 rotation_env%matrix_sx => sx_re
2100 rotation_env%matrix_sx_im => sx_im
2101 rotation_env%matrix_gx => g_re
2102 rotation_env%matrix_gx_im => g_im
2103 rotation_env%matrix_ref_inv_sqrt => f_re
2104 rotation_env%matrix_ref_inv_sqrt_im => f_im
2105 rotation_env%buf1_k_k_sym => sym_buffer(1)
2106 rotation_env%buf2_k_k_sym => sym_buffer(2)
2107 rotation_env%buf3_k_k_sym => sym_buffer(3)
2108 rotation_env%buf4_k_k_sym => sym_buffer(4)
2109 rotation_env%buf1_k_k_nosym => nonsym_buffer(1)
2110 rotation_env%buf2_k_k_nosym => nonsym_buffer(2)
2111 rotation_env%buf3_k_k_nosym => nonsym_buffer(3)
2112 rotation_env%buf1_n_k => nk_buffer(1)
2113 rotation_env%buf1_n_k_dp => nk_buffer(2)
2114 rotation_env%rot_mat_u => rot_matrix(1)
2115 rotation_env%rot_mat_u_im => rot_matrix(2)
2116 rotation_env%rot_mat_x => rot_matrix(3)
2117 rotation_env%rot_mat_x_im => rot_matrix(4)
2118 rotation_env%rot_mat_dedu => rot_matrix(5)
2119 rotation_env%rot_mat_dedu_im => rot_matrix(6)
2120 rotation_env%rot_mat_gx => rot_matrix(7)
2121 rotation_env%rot_mat_gx_im => rot_matrix(8)
2122 rotation_env%rot_mat_evec_re => rot_matrix(9)
2123 rotation_env%rot_mat_evec_im => rot_matrix(10)
2124 rotation_env%rot_mat_evals => rotation_evals
2125
2126 c_re_p => c_re
2127 c_im_p => c_im
2128 s_re_p => s_re
2129 s_im_p => s_im
2130 CALL qs_ot_get_orbitals_ref_complex(c_re_p, c_im_p, s_re_p, s_im_p, rotation_env)
2131 CALL get_complex_pair(c_re, c_im, c_dense)
2132 error_c = maxval(abs(c_dense - c_reference))
2133
2134 hc_dense = 2.0_dp*matmul(hamiltonian, c_dense)
2135 DO j = 1, k
2136 hc_dense(:, j) = occupation(j)*hc_dense(:, j)
2137 END DO
2138 CALL put_complex_pair(hc_re, hc_im, hc_dense)
2139 hc_re_p => hc_re
2140 hc_im_p => hc_im
2141 CALL qs_ot_get_derivative_ref_complex(hc_re_p, hc_im_p, rotation_env)
2142 CALL get_complex_pair(g_re, g_im, gradient_x)
2143 CALL get_complex_pair(rot_matrix(7), rot_matrix(8), gradient_a)
2144
2145 hstep = 1.0e-6_dp
2146 fd_ref = (rotated_polar_energy(x_dense + hstep*direction_x, generator, &
2147 hamiltonian, occupation) - &
2148 rotated_polar_energy(x_dense - hstep*direction_x, generator, &
2149 hamiltonian, occupation))/(2.0_dp*hstep)
2150 predicted_ref = real(sum(conjg(gradient_x)*direction_x), kind=dp)
2151 error_ref = abs(predicted_ref - fd_ref)
2152
2153 fd_rot = (rotated_polar_energy(x_dense, generator + hstep*direction_a, &
2154 hamiltonian, occupation) - &
2155 rotated_polar_energy(x_dense, generator - hstep*direction_a, &
2156 hamiltonian, occupation))/(2.0_dp*hstep)
2157 predicted_rot = 0.5_dp*(sum(real(gradient_a, kind=dp)*real(direction_a, kind=dp)) + &
2158 sum(aimag(gradient_a)*aimag(direction_a)))
2159 error_rot = abs(predicted_rot - fd_rot)
2160
2161 IF (error_c > 5.0e-11_dp .OR. error_ref > 5.0e-8_dp .OR. &
2162 error_rot > 5.0e-8_dp) nfail = nfail + 1
2163 IF (para_env%is_source()) THEN
2164 WRITE (io_unit, '(A,1X,ES13.6)') 'complex finite REF rotation orbital error', error_c
2165 WRITE (io_unit, '(A,3(1X,ES13.6))') 'complex finite REF slope fd/predicted/error', &
2166 fd_ref, predicted_ref, error_ref
2167 WRITE (io_unit, '(A,3(1X,ES13.6))') 'complex finite rotation slope fd/predicted/error', &
2168 fd_rot, predicted_rot, error_rot
2169 END IF
2170
2171 DO i = 1, SIZE(rot_matrix)
2172 CALL dbcsr_release(rot_matrix(i))
2173 END DO
2174 DO i = 1, SIZE(nk_buffer)
2175 CALL dbcsr_release(nk_buffer(i))
2176 END DO
2177 DO i = 1, SIZE(nonsym_buffer)
2178 CALL dbcsr_release(nonsym_buffer(i))
2179 END DO
2180 DO i = 1, SIZE(sym_buffer)
2181 CALL dbcsr_release(sym_buffer(i))
2182 END DO
2183 CALL dbcsr_release(f_im)
2184 CALL dbcsr_release(f_re)
2185 CALL dbcsr_release(s_im)
2186 CALL dbcsr_release(s_re)
2187 CALL dbcsr_release(g_im)
2188 CALL dbcsr_release(g_re)
2189 CALL dbcsr_release(hc_im)
2190 CALL dbcsr_release(hc_re)
2191 CALL dbcsr_release(c_im)
2192 CALL dbcsr_release(c_re)
2193 CALL dbcsr_release(sx_im)
2194 CALL dbcsr_release(sx_re)
2195 CALL dbcsr_release(x_im)
2196 CALL dbcsr_release(x_re)
2197 CALL dbcsr_distribution_release(dist_kk)
2198 CALL dbcsr_distribution_release(dist_nn)
2199 CALL dbcsr_distribution_release(dist_nk)
2200 DEALLOCATE (pgrid, row_dist, col_dist, row_n, row_k, col_n, col_k)
2201 CALL cp_blacs_env_release(blacs_env)
2202 CALL mp_para_env_release(local_para_env)
2203
2204 END SUBROUTINE test_complex_ref_rotation
2205
2206! **************************************************************************************************
2207!> \brief Check the finite complex rotation Frechet derivative against central differences.
2208!> \param para_env parallel environment
2209!> \param nfail accumulated number of failures
2210! **************************************************************************************************
2211 SUBROUTINE test_complex_rotation_frechet(para_env, nfail)
2212 TYPE(mp_para_env_type), POINTER :: para_env
2213 INTEGER, INTENT(INOUT) :: nfail
2214
2215 INTEGER, PARAMETER :: m = 3
2216
2217 COMPLEX(KIND=dp), DIMENSION(m, m) :: direction, fd, generator, predicted
2218 INTEGER :: i, j
2219 REAL(KIND=dp) :: error, hstep
2220 REAL(KIND=dp), DIMENSION(m) :: evals
2221
2222 evals = [-0.7_dp, 0.2000003_dp, 1.1_dp]
2223 generator(:, :) = cmplx(0.0_dp, 0.0_dp, kind=dp)
2224 DO i = 1, m
2225 generator(i, i) = cmplx(0.0_dp, -evals(i), kind=dp)
2226 END DO
2227 direction(:, :) = cmplx(0.0_dp, 0.0_dp, kind=dp)
2228 direction(1, 2) = cmplx(0.17_dp, -0.09_dp, kind=dp)
2229 direction(2, 1) = -conjg(direction(1, 2))
2230 direction(1, 3) = cmplx(-0.06_dp, 0.13_dp, kind=dp)
2231 direction(3, 1) = -conjg(direction(1, 3))
2232 direction(2, 3) = cmplx(0.11_dp, 0.04_dp, kind=dp)
2233 direction(3, 2) = -conjg(direction(2, 3))
2234 direction(1, 1) = cmplx(0.0_dp, 0.08_dp, kind=dp)
2235 direction(2, 2) = cmplx(0.0_dp, -0.03_dp, kind=dp)
2236 direction(3, 3) = cmplx(0.0_dp, 0.05_dp, kind=dp)
2237
2238 DO i = 1, m
2239 DO j = 1, m
2240 predicted(i, j) = qs_ot_complex_exp_frechet_kernel(evals(i), evals(j))* &
2241 direction(i, j)
2242 END DO
2243 END DO
2244 hstep = 1.0e-6_dp
2245 fd = (dense_antihermitian_exp(generator + hstep*direction) - &
2246 dense_antihermitian_exp(generator - hstep*direction))/(2.0_dp*hstep)
2247 error = maxval(abs(fd - predicted))
2248 IF (error > 5.0e-9_dp) nfail = nfail + 1
2249 IF (para_env%is_source()) THEN
2250 WRITE (io_unit, '(A,1X,ES13.6)') 'complex rotation Frechet finite-difference error', error
2251 END IF
2252
2253 END SUBROUTINE test_complex_rotation_frechet
2254
2255! **************************************************************************************************
2256!> \brief Check union handling for distinct real/imaginary DBCSR block patterns.
2257!> \param para_env parallel environment
2258!> \param nfail accumulated number of failures
2259! **************************************************************************************************
2260 SUBROUTINE test_sparse_frechet_patterns(para_env, nfail)
2261 TYPE(mp_para_env_type), POINTER :: para_env
2262 INTEGER, INTENT(INOUT) :: nfail
2263
2264 COMPLEX(KIND=dp) :: expected, value
2265 INTEGER, DIMENSION(:), POINTER :: col_dist, col_size, row_dist, row_size
2266 INTEGER, DIMENSION(:, :), POINTER :: pgrid
2267 LOGICAL :: found
2268 REAL(kind=dp) :: error
2269 REAL(kind=dp), DIMENSION(2) :: evals
2270 REAL(kind=dp), DIMENSION(:, :), POINTER :: block
2271 TYPE(dbcsr_distribution_type) :: dist
2272 TYPE(dbcsr_type) :: inner_im, inner_re, outer_im, outer_re
2273
2274 ALLOCATE (pgrid(0:0, 0:0), row_dist(2), col_dist(2), row_size(2), col_size(2))
2275 pgrid(:, :) = 0
2276 row_dist(:) = 0
2277 col_dist(:) = 0
2278 row_size(:) = 1
2279 col_size(:) = 1
2280 CALL dbcsr_distribution_new(dist, group=mp_comm_self%get_handle(), pgrid=pgrid, &
2281 row_dist=row_dist, col_dist=col_dist)
2282 CALL dbcsr_create(inner_re, 'frechet_sparse_re', dist, dbcsr_type_no_symmetry, row_size, col_size)
2283 CALL dbcsr_reserve_blocks(inner_re, rows=[1], cols=[2])
2284 CALL dbcsr_finalize(inner_re)
2285 CALL dbcsr_create(inner_im, 'frechet_sparse_im', dist, dbcsr_type_no_symmetry, row_size, col_size)
2286 CALL dbcsr_reserve_blocks(inner_im, rows=[2], cols=[1])
2287 CALL dbcsr_finalize(inner_im)
2288 CALL dbcsr_put_block(inner_re, 1, 2, reshape([0.37_dp], [1, 1]))
2289 CALL dbcsr_put_block(inner_im, 2, 1, reshape([-0.21_dp], [1, 1]))
2290
2291 evals = [-0.4_dp, 0.9_dp]
2292 CALL qs_ot_apply_complex_frechet_dbcsr(evals, inner_re, inner_im, outer_re, outer_im)
2293 error = 0.0_dp
2294 expected = qs_ot_complex_exp_frechet_kernel(evals(1), evals(2))*cmplx(0.37_dp, 0.0_dp, kind=dp)
2295 value = cmplx(0.0_dp, 0.0_dp, kind=dp)
2296 CALL dbcsr_get_block_p(outer_re, 1, 2, block, found)
2297 IF (found) value = value + cmplx(block(1, 1), 0.0_dp, kind=dp)
2298 CALL dbcsr_get_block_p(outer_im, 1, 2, block, found)
2299 IF (found) value = value + cmplx(0.0_dp, block(1, 1), kind=dp)
2300 error = max(error, abs(value - expected))
2301
2302 expected = qs_ot_complex_exp_frechet_kernel(evals(2), evals(1))*cmplx(0.0_dp, -0.21_dp, kind=dp)
2303 value = cmplx(0.0_dp, 0.0_dp, kind=dp)
2304 CALL dbcsr_get_block_p(outer_re, 2, 1, block, found)
2305 IF (found) value = value + cmplx(block(1, 1), 0.0_dp, kind=dp)
2306 CALL dbcsr_get_block_p(outer_im, 2, 1, block, found)
2307 IF (found) value = value + cmplx(0.0_dp, block(1, 1), kind=dp)
2308 error = max(error, abs(value - expected))
2309 IF (error > 5.0e-14_dp) nfail = nfail + 1
2310 IF (para_env%is_source()) THEN
2311 WRITE (io_unit, '(A,1X,ES13.6)') 'sparse complex Frechet block-union error', error
2312 END IF
2313
2314 CALL dbcsr_release(outer_im)
2315 CALL dbcsr_release(outer_re)
2316 CALL dbcsr_release(inner_im)
2317 CALL dbcsr_release(inner_re)
2318 CALL dbcsr_distribution_release(dist)
2319 DEALLOCATE (pgrid, row_dist, col_dist, row_size, col_size)
2320
2321 END SUBROUTINE test_sparse_frechet_patterns
2322! **************************************************************************************************
2323!> \brief Check the finite complex STRICT chart and its analytic derivative.
2324!> \param dist_nk distribution for orbital-shaped matrices
2325!> \param dist_kk distribution for occupied-space matrices
2326!> \param row_n AO block sizes
2327!> \param row_k occupied-space row block sizes
2328!> \param col_k occupied-space column block sizes
2329!> \param para_env parallel environment
2330!> \param nfail accumulated number of failures
2331! **************************************************************************************************
2332 SUBROUTINE test_complex_strict_derivative(dist_nk, dist_kk, row_n, row_k, col_k, para_env, nfail)
2333 TYPE(dbcsr_distribution_type), INTENT(IN) :: dist_nk, dist_kk
2334 INTEGER, DIMENSION(:), INTENT(INOUT), POINTER :: row_n, row_k, col_k
2335 TYPE(mp_para_env_type), POINTER :: para_env
2336 INTEGER, INTENT(INOUT) :: nfail
2337
2338 CHARACTER(LEN=32) :: name
2339 COMPLEX(KIND=dp), DIMENSION(k, k) :: direction_a, generator, gradient_a, &
2340 inverse_sqrt, rotation
2341 COMPLEX(KIND=dp), DIMENSION(n, k) :: c0, direction, gradient, hc, orbitals, &
2342 orbitals_base, sc0, sx, x
2343 COMPLEX(KIND=dp), DIMENSION(n, n) :: hamiltonian, overlap
2344 INTEGER :: i
2345 LOGICAL :: found
2346 REAL(KIND=dp) :: error, error_rot, fd_slope, &
2347 fd_slope_rot, hstep, ortho_error, &
2348 slope, slope_rot
2349 REAL(KIND=dp), DIMENSION(:, :), POINTER :: block
2350 REAL(KIND=dp), DIMENSION(k) :: occupation
2351 REAL(KIND=dp), DIMENSION(k), TARGET :: rotation_evals
2352 TYPE(cp_blacs_env_type), POINTER :: blacs_env
2353 TYPE(dbcsr_type), DIMENSION(10), TARGET :: rot_matrix
2354 TYPE(dbcsr_type), DIMENSION(15), TARGET :: nk
2355 TYPE(dbcsr_type), DIMENSION(25), TARGET :: kk
2356 TYPE(dbcsr_type), POINTER :: hc_im_p, hc_re_p, sx_im_p, sx_re_p, &
2357 x_im_p, x_re_p
2358 TYPE(mp_para_env_type), POINTER :: strict_para_env
2359 TYPE(qs_ot_type) :: strict_env
2360
2361 NULLIFY (blacs_env, hc_im_p, hc_re_p, strict_para_env, x_im_p, x_re_p, sx_im_p, sx_re_p)
2362 ALLOCATE (strict_para_env)
2363 CALL strict_para_env%from_dup(mp_comm_self)
2364 CALL cp_blacs_env_create(blacs_env=blacs_env, para_env=strict_para_env)
2365 DO i = 1, SIZE(nk)
2366 WRITE (name, '(A,I0)') 'complex_strict_nk_', i
2367 CALL create_one_block(nk(i), trim(name), dist_nk, row_n, col_k)
2368 END DO
2369 DO i = 1, SIZE(kk)
2370 WRITE (name, '(A,I0)') 'complex_strict_kk_', i
2371 CALL create_one_block(kk(i), trim(name), dist_kk, row_k, col_k)
2372 END DO
2373 DO i = 1, SIZE(rot_matrix)
2374 WRITE (name, '(A,I0)') 'complex_strict_rot_', i
2375 CALL create_one_block(rot_matrix(i), trim(name), dist_kk, row_k, col_k)
2376 END DO
2377
2378 overlap(:, :) = cmplx(0.0_dp, 0.0_dp, kind=dp)
2379 overlap(1, 1) = cmplx(1.0_dp, 0.0_dp, kind=dp)
2380 overlap(2, 2) = cmplx(1.2_dp, 0.0_dp, kind=dp)
2381 overlap(3, 3) = cmplx(1.4_dp, 0.0_dp, kind=dp)
2382 overlap(1, 2) = cmplx(0.06_dp, 0.03_dp, kind=dp)
2383 overlap(2, 1) = conjg(overlap(1, 2))
2384 overlap(1, 3) = cmplx(-0.04_dp, 0.02_dp, kind=dp)
2385 overlap(3, 1) = conjg(overlap(1, 3))
2386 overlap(2, 3) = cmplx(0.05_dp, -0.01_dp, kind=dp)
2387 overlap(3, 2) = conjg(overlap(2, 3))
2388
2389 c0(:, 1) = [cmplx(1.0_dp, 0.10_dp, kind=dp), cmplx(0.12_dp, -0.08_dp, kind=dp), &
2390 cmplx(-0.07_dp, 0.04_dp, kind=dp)]
2391 c0(:, 2) = [cmplx(0.09_dp, 0.05_dp, kind=dp), cmplx(0.85_dp, -0.12_dp, kind=dp), &
2392 cmplx(0.11_dp, 0.06_dp, kind=dp)]
2393 CALL dense_inverse_sqrt(matmul(conjg(transpose(c0)), matmul(overlap, c0)), inverse_sqrt)
2394 c0 = matmul(c0, inverse_sqrt)
2395
2396 x(:, :) = cmplx(0.0_dp, 0.0_dp, kind=dp)
2397 x(3, 1) = cmplx(0.23_dp, 0.17_dp, kind=dp)
2398 x(3, 2) = cmplx(-0.14_dp, 0.21_dp, kind=dp)
2399 x = x - matmul(c0, matmul(conjg(transpose(c0)), matmul(overlap, x)))
2400 direction(:, :) = cmplx(0.0_dp, 0.0_dp, kind=dp)
2401 direction(3, 1) = cmplx(0.11_dp, -0.09_dp, kind=dp)
2402 direction(3, 2) = cmplx(0.07_dp, 0.13_dp, kind=dp)
2403 direction = direction - matmul(c0, matmul(conjg(transpose(c0)), matmul(overlap, direction)))
2404
2405 hamiltonian(:, :) = cmplx(0.0_dp, 0.0_dp, kind=dp)
2406 hamiltonian(1, 1) = cmplx(-0.70_dp, 0.0_dp, kind=dp)
2407 hamiltonian(2, 2) = cmplx(-0.20_dp, 0.0_dp, kind=dp)
2408 hamiltonian(3, 3) = cmplx(0.90_dp, 0.0_dp, kind=dp)
2409 hamiltonian(1, 2) = cmplx(0.08_dp, 0.04_dp, kind=dp)
2410 hamiltonian(2, 1) = conjg(hamiltonian(1, 2))
2411 hamiltonian(1, 3) = cmplx(-0.12_dp, 0.07_dp, kind=dp)
2412 hamiltonian(3, 1) = conjg(hamiltonian(1, 3))
2413 hamiltonian(2, 3) = cmplx(0.05_dp, -0.09_dp, kind=dp)
2414 hamiltonian(3, 2) = conjg(hamiltonian(2, 3))
2415
2416 generator(:, :) = cmplx(0.0_dp, 0.0_dp, kind=dp)
2417 generator(1, 1) = cmplx(0.0_dp, 0.16_dp, kind=dp)
2418 generator(2, 2) = cmplx(0.0_dp, -0.09_dp, kind=dp)
2419 generator(1, 2) = cmplx(0.18_dp, -0.07_dp, kind=dp)
2420 generator(2, 1) = -conjg(generator(1, 2))
2421 direction_a(:, :) = cmplx(0.0_dp, 0.0_dp, kind=dp)
2422 direction_a(1, 1) = cmplx(0.0_dp, -0.05_dp, kind=dp)
2423 direction_a(2, 2) = cmplx(0.0_dp, 0.04_dp, kind=dp)
2424 direction_a(1, 2) = cmplx(-0.08_dp, 0.06_dp, kind=dp)
2425 direction_a(2, 1) = -conjg(direction_a(1, 2))
2426 occupation = [1.0_dp, 0.34_dp]
2427
2428 CALL strict_orbitals(c0, x, overlap, orbitals_base)
2429 rotation = dense_antihermitian_exp(generator)
2430 orbitals = matmul(orbitals_base, rotation)
2431 hc = 2.0_dp*matmul(hamiltonian, orbitals)
2432 DO i = 1, k
2433 hc(:, i) = occupation(i)*hc(:, i)
2434 END DO
2435 sx = matmul(overlap, x)
2436 sc0 = matmul(overlap, c0)
2437 CALL put_complex_pair(nk(1), nk(2), x)
2438 CALL put_complex_pair(nk(3), nk(4), sx)
2439 CALL put_complex_pair(nk(5), nk(6), hc)
2440 CALL put_complex_pair(nk(9), nk(10), c0)
2441 CALL put_complex_pair(nk(11), nk(12), sc0)
2442 CALL put_complex_pair(rot_matrix(3), rot_matrix(4), generator)
2443
2444 strict_env%has_complex_kpoint_state = .true.
2445 strict_env%settings%ot_algorithm = 'TOD'
2446 strict_env%settings%do_rotation = .true.
2447 strict_env%para_env => strict_para_env
2448 strict_env%blacs_env => blacs_env
2449 strict_env%matrix_x => nk(1)
2450 strict_env%matrix_x_im => nk(2)
2451 strict_env%matrix_sx => nk(3)
2452 strict_env%matrix_sx_im => nk(4)
2453 strict_env%matrix_gx => nk(7)
2454 strict_env%matrix_gx_im => nk(8)
2455 strict_env%matrix_c0 => nk(9)
2456 strict_env%matrix_c0_im => nk(10)
2457 strict_env%matrix_sc0 => nk(11)
2458 strict_env%matrix_sc0_im => nk(12)
2459 strict_env%matrix_buf_nk => nk(13)
2460 strict_env%matrix_buf_nk_im => nk(14)
2461 strict_env%matrix_tmp_nk => nk(15)
2462 strict_env%matrix_p => kk(1)
2463 strict_env%matrix_p_im => kk(2)
2464 strict_env%matrix_r => kk(3)
2465 strict_env%matrix_r_im => kk(4)
2466 strict_env%matrix_sinp => kk(5)
2467 strict_env%matrix_sinp_im => kk(6)
2468 strict_env%matrix_cosp => kk(7)
2469 strict_env%matrix_cosp_im => kk(8)
2470 strict_env%matrix_sinp_b => kk(9)
2471 strict_env%matrix_cosp_b => kk(10)
2472 strict_env%matrix_buf1 => kk(11)
2473 strict_env%matrix_buf1_im => kk(12)
2474 strict_env%matrix_buf2 => kk(13)
2475 strict_env%matrix_buf2_im => kk(14)
2476 strict_env%matrix_buf3 => kk(15)
2477 strict_env%matrix_buf3_im => kk(16)
2478 strict_env%matrix_buf4 => kk(17)
2479 strict_env%matrix_buf4_im => kk(18)
2480 strict_env%matrix_os => kk(19)
2481 strict_env%matrix_os_im => kk(20)
2482 strict_env%matrix_buf1_ortho => kk(21)
2483 strict_env%matrix_buf1_ortho_im => kk(22)
2484 strict_env%matrix_buf2_ortho => kk(23)
2485 strict_env%matrix_buf2_ortho_im => kk(24)
2486 strict_env%matrix_tmp_ortho => kk(25)
2487 strict_env%rot_mat_u => rot_matrix(1)
2488 strict_env%rot_mat_u_im => rot_matrix(2)
2489 strict_env%rot_mat_x => rot_matrix(3)
2490 strict_env%rot_mat_x_im => rot_matrix(4)
2491 strict_env%rot_mat_dedu => rot_matrix(5)
2492 strict_env%rot_mat_dedu_im => rot_matrix(6)
2493 strict_env%rot_mat_gx => rot_matrix(7)
2494 strict_env%rot_mat_gx_im => rot_matrix(8)
2495 strict_env%rot_mat_evec_re => rot_matrix(9)
2496 strict_env%rot_mat_evec_im => rot_matrix(10)
2497 strict_env%rot_mat_evals => rotation_evals
2498 ALLOCATE (strict_env%evals(k), strict_env%dum(k))
2499
2500 x_re_p => nk(1)
2501 x_im_p => nk(2)
2502 sx_re_p => nk(3)
2503 sx_im_p => nk(4)
2504 hc_re_p => nk(5)
2505 hc_im_p => nk(6)
2506 CALL qs_ot_get_p_complex(x_re_p, x_im_p, sx_re_p, sx_im_p, strict_env)
2507 CALL qs_ot_get_derivative_complex(hc_re_p, hc_im_p, strict_env)
2508
2509 gradient(:, :) = cmplx(0.0_dp, 0.0_dp, kind=dp)
2510 CALL dbcsr_get_block_p(nk(7), 1, 1, block, found)
2511 IF (found) gradient = gradient + cmplx(block, 0.0_dp, kind=dp)
2512 CALL dbcsr_get_block_p(nk(8), 1, 1, block, found)
2513 IF (found) gradient = gradient + cmplx(0.0_dp, block, kind=dp)
2514 CALL get_complex_pair(rot_matrix(7), rot_matrix(8), gradient_a)
2515
2516 hstep = 1.0e-6_dp
2517 fd_slope = (strict_rotated_energy(c0, x + hstep*direction, generator, overlap, &
2518 hamiltonian, occupation) - &
2519 strict_rotated_energy(c0, x - hstep*direction, generator, overlap, &
2520 hamiltonian, occupation))/(2.0_dp*hstep)
2521 slope = real(sum(conjg(gradient)*direction), kind=dp)
2522 error = abs(slope - fd_slope)
2523 fd_slope_rot = (strict_rotated_energy(c0, x, generator + hstep*direction_a, overlap, &
2524 hamiltonian, occupation) - &
2525 strict_rotated_energy(c0, x, generator - hstep*direction_a, overlap, &
2526 hamiltonian, occupation))/(2.0_dp*hstep)
2527 slope_rot = 0.5_dp*(sum(real(gradient_a, kind=dp)*real(direction_a, kind=dp)) + &
2528 sum(aimag(gradient_a)*aimag(direction_a)))
2529 error_rot = abs(slope_rot - fd_slope_rot)
2530 ortho_error = maxval(abs(matmul(conjg(transpose(orbitals)), matmul(overlap, orbitals)) - &
2531 reshape([cmplx(1.0_dp, 0.0_dp, kind=dp), &
2532 cmplx(0.0_dp, 0.0_dp, kind=dp), &
2533 cmplx(0.0_dp, 0.0_dp, kind=dp), &
2534 cmplx(1.0_dp, 0.0_dp, kind=dp)], [k, k])))
2535 IF (error > 5.0e-8_dp .OR. error_rot > 5.0e-8_dp .OR. &
2536 ortho_error > 5.0e-13_dp) nfail = nfail + 1
2537 IF (para_env%is_source()) THEN
2538 WRITE (io_unit, '(A,3(1X,ES13.6))') 'complex finite-STRICT derivative: fd/predicted/error', &
2539 fd_slope, slope, error
2540 WRITE (io_unit, '(A,3(1X,ES13.6))') 'complex finite-STRICT rotation: fd/predicted/error', &
2541 fd_slope_rot, slope_rot, error_rot
2542 WRITE (io_unit, '(A,1X,ES13.6)') 'complex finite-STRICT orthonormality error', ortho_error
2543 END IF
2544
2545 DEALLOCATE (strict_env%evals, strict_env%dum)
2546 DO i = SIZE(rot_matrix), 1, -1
2547 CALL dbcsr_release(rot_matrix(i))
2548 END DO
2549 DO i = SIZE(kk), 1, -1
2550 CALL dbcsr_release(kk(i))
2551 END DO
2552 DO i = SIZE(nk), 1, -1
2553 CALL dbcsr_release(nk(i))
2554 END DO
2555 CALL cp_blacs_env_release(blacs_env)
2556 CALL mp_para_env_release(strict_para_env)
2557
2558 END SUBROUTINE test_complex_strict_derivative
2559! **************************************************************************************************
2560!> \brief Check gauge invariance, Hermiticity, and positivity of the complex inverse operator.
2561!> \param para_env parallel environment
2562!> \param nfail accumulated number of failures
2563! **************************************************************************************************
2564 SUBROUTINE test_complex_preconditioner_gauge(para_env, nfail)
2565 TYPE(mp_para_env_type), POINTER :: para_env
2566 INTEGER, INTENT(INOUT) :: nfail
2567
2568 COMPLEX(KIND=dp) :: value
2569 COMPLEX(KIND=dp), DIMENSION(k) :: phases
2570 COMPLEX(KIND=dp), DIMENSION(n, k) :: c0, c0_phase, c0_rot, gradient, &
2571 gradient_phase, gradient_rot, output, &
2572 output_phase, output_rot
2573 COMPLEX(KIND=dp), DIMENSION(n, n) :: eigenvectors, full_all_vectors, h, inverse, &
2574 inverse_rot, inverse_s, inverse_single, inverse_t, reference_s, reference_t, s, &
2575 spectral_metric, t
2576 INTEGER :: i, j
2577 REAL(KIND=dp) :: covariant_descent, covariant_gauge_error, full_all_descent, &
2578 full_all_gauge_error, full_all_orthogonality_error, full_single_hermitian_error, &
2579 gauge_error, hermitian_error, low_rank_descent, low_rank_gauge_error, s_error, t_error
2580 REAL(KIND=dp), DIMENSION(k) :: c0_evals
2581 REAL(KIND=dp), DIMENSION(n) :: eigenvalues, eigenvalues_single
2582 TYPE(cp_blacs_env_type), POINTER :: blacs_env
2583 TYPE(cp_cfm_type) :: matrix_c0, matrix_c0_phase, matrix_c0_rot, matrix_gradient, &
2584 matrix_gradient_phase, matrix_gradient_rot, matrix_h, matrix_output, matrix_output_phase, &
2585 matrix_output_rot, matrix_s, matrix_t
2586 TYPE(cp_fm_struct_type), POINTER :: fmstruct_nk, fmstruct_nn
2587 TYPE(preconditioner_type) :: covariant, covariant_rot, full_all, full_all_phase, low_rank, &
2588 low_rank_rot, preconditioner, preconditioner_rot, preconditioner_s, &
2589 preconditioner_single, preconditioner_t
2590
2591 NULLIFY (blacs_env, fmstruct_nk, fmstruct_nn)
2592 CALL cp_blacs_env_create(blacs_env=blacs_env, para_env=para_env)
2593 CALL cp_fm_struct_create(fmstruct=fmstruct_nn, para_env=para_env, context=blacs_env, &
2594 nrow_global=n, ncol_global=n)
2595 CALL cp_fm_struct_create(fmstruct=fmstruct_nk, para_env=para_env, context=blacs_env, &
2596 nrow_global=n, ncol_global=k)
2597 CALL cp_cfm_create(matrix_h, fmstruct_nn, name='complex preconditioner H')
2598 CALL cp_cfm_create(matrix_s, fmstruct_nn, name='complex preconditioner S')
2599 CALL cp_cfm_create(matrix_t, fmstruct_nn, name='complex preconditioner T')
2600 CALL cp_cfm_create(matrix_c0, fmstruct_nk, name='complex preconditioner C')
2601 CALL cp_cfm_create(matrix_c0_phase, fmstruct_nk, name='complex preconditioner phased C')
2602 CALL cp_cfm_create(matrix_c0_rot, fmstruct_nk, name='complex preconditioner C U')
2603 CALL cp_cfm_create(matrix_gradient, fmstruct_nk, name='complex preconditioner gradient')
2604 CALL cp_cfm_create(matrix_gradient_phase, fmstruct_nk, &
2605 name='complex preconditioner phased gradient')
2606 CALL cp_cfm_create(matrix_gradient_rot, fmstruct_nk, &
2607 name='complex preconditioner rotated gradient')
2608 CALL cp_cfm_create(matrix_output, fmstruct_nk, name='complex preconditioner output')
2609 CALL cp_cfm_create(matrix_output_phase, fmstruct_nk, &
2610 name='complex preconditioner phased output')
2611 CALL cp_cfm_create(matrix_output_rot, fmstruct_nk, &
2612 name='complex preconditioner rotated output')
2613 CALL cp_cfm_set_all(matrix_h, cmplx(0.0_dp, 0.0_dp, kind=dp))
2614 CALL cp_cfm_set_all(matrix_s, cmplx(0.0_dp, 0.0_dp, kind=dp))
2615 CALL cp_cfm_set_all(matrix_t, cmplx(0.0_dp, 0.0_dp, kind=dp))
2616 CALL cp_cfm_set_all(matrix_c0, cmplx(0.0_dp, 0.0_dp, kind=dp))
2617 CALL cp_cfm_set_all(matrix_c0_phase, cmplx(0.0_dp, 0.0_dp, kind=dp))
2618 CALL cp_cfm_set_all(matrix_c0_rot, cmplx(0.0_dp, 0.0_dp, kind=dp))
2619
2620 h(:, :) = cmplx(0.0_dp, 0.0_dp, kind=dp)
2621 h(1, 1) = cmplx(-0.35_dp, 0.0_dp, kind=dp)
2622 h(2, 2) = cmplx(0.15_dp, 0.0_dp, kind=dp)
2623 h(3, 3) = cmplx(0.90_dp, 0.0_dp, kind=dp)
2624 h(1, 2) = cmplx(0.08_dp, 0.04_dp, kind=dp)
2625 h(2, 1) = conjg(h(1, 2))
2626 h(1, 3) = cmplx(-0.06_dp, 0.03_dp, kind=dp)
2627 h(3, 1) = conjg(h(1, 3))
2628 h(2, 3) = cmplx(0.05_dp, -0.07_dp, kind=dp)
2629 h(3, 2) = conjg(h(2, 3))
2630 s(:, :) = cmplx(0.0_dp, 0.0_dp, kind=dp)
2631 s(1, 1) = cmplx(1.0_dp, 0.0_dp, kind=dp)
2632 s(2, 2) = cmplx(1.0_dp, 0.0_dp, kind=dp)
2633 s(3, 3) = cmplx(1.4_dp, 0.0_dp, kind=dp)
2634 s(1, 3) = cmplx(0.08_dp, 0.03_dp, kind=dp)
2635 s(3, 1) = conjg(s(1, 3))
2636 s(2, 3) = cmplx(-0.04_dp, 0.02_dp, kind=dp)
2637 s(3, 2) = conjg(s(2, 3))
2638 t = h
2639 DO i = 1, n
2640 t(i, i) = t(i, i) + cmplx(1.50_dp, 0.0_dp, kind=dp)
2641 END DO
2642 c0(:, :) = cmplx(0.0_dp, 0.0_dp, kind=dp)
2643 c0(1, 1) = cmplx(1.0_dp, 0.0_dp, kind=dp)
2644 c0(2, 2) = cmplx(1.0_dp, 0.0_dp, kind=dp)
2645 c0_rot(:, :) = cmplx(0.0_dp, 0.0_dp, kind=dp)
2646 c0_rot(1, 1) = cmplx(sqrt(0.5_dp), 0.0_dp, kind=dp)
2647 c0_rot(1, 2) = cmplx(0.0_dp, sqrt(0.5_dp), kind=dp)
2648 c0_rot(2, 1) = cmplx(0.0_dp, sqrt(0.5_dp), kind=dp)
2649 c0_rot(2, 2) = cmplx(sqrt(0.5_dp), 0.0_dp, kind=dp)
2650 phases = [exp(cmplx(0.0_dp, 0.43_dp, kind=dp)), &
2651 exp(cmplx(0.0_dp, -0.67_dp, kind=dp))]
2652 DO j = 1, k
2653 c0_phase(:, j) = phases(j)*c0(:, j)
2654 END DO
2655 gradient(:, 1) = [cmplx(0.17_dp, 0.02_dp, kind=dp), &
2656 cmplx(-0.08_dp, 0.05_dp, kind=dp), &
2657 cmplx(0.11_dp, -0.04_dp, kind=dp)]
2658 gradient(:, 2) = [cmplx(-0.03_dp, 0.06_dp, kind=dp), &
2659 cmplx(0.09_dp, -0.02_dp, kind=dp), &
2660 cmplx(-0.14_dp, 0.05_dp, kind=dp)]
2661 DO j = 1, k
2662 gradient_phase(:, j) = phases(j)*gradient(:, j)
2663 END DO
2664 gradient_rot = matmul(gradient, c0_rot(1:k, :))
2665 c0_evals = [real(h(1, 1), kind=dp), real(h(2, 2), kind=dp)]
2666
2667 DO j = 1, n
2668 DO i = 1, n
2669 CALL cp_cfm_set_element(matrix_h, i, j, h(i, j))
2670 CALL cp_cfm_set_element(matrix_s, i, j, s(i, j))
2671 CALL cp_cfm_set_element(matrix_t, i, j, t(i, j))
2672 END DO
2673 END DO
2674 DO j = 1, k
2675 DO i = 1, n
2676 CALL cp_cfm_set_element(matrix_c0, i, j, c0(i, j))
2677 CALL cp_cfm_set_element(matrix_c0_phase, i, j, c0_phase(i, j))
2678 CALL cp_cfm_set_element(matrix_c0_rot, i, j, c0_rot(i, j))
2679 CALL cp_cfm_set_element(matrix_gradient, i, j, gradient(i, j))
2680 CALL cp_cfm_set_element(matrix_gradient_phase, i, j, gradient_phase(i, j))
2681 CALL cp_cfm_set_element(matrix_gradient_rot, i, j, gradient_rot(i, j))
2682 END DO
2683 END DO
2684
2685 CALL init_preconditioner(full_all, para_env, blacs_env)
2686 CALL init_preconditioner(full_all_phase, para_env, blacs_env)
2687 CALL init_preconditioner(low_rank, para_env, blacs_env)
2688 CALL init_preconditioner(low_rank_rot, para_env, blacs_env)
2689 CALL init_preconditioner(covariant, para_env, blacs_env)
2690 CALL init_preconditioner(covariant_rot, para_env, blacs_env)
2691 CALL init_preconditioner(preconditioner, para_env, blacs_env)
2692 CALL init_preconditioner(preconditioner_rot, para_env, blacs_env)
2693 CALL init_preconditioner(preconditioner_s, para_env, blacs_env)
2694 CALL init_preconditioner(preconditioner_single, para_env, blacs_env)
2695 CALL init_preconditioner(preconditioner_t, para_env, blacs_env)
2696 CALL make_complex_full_single_inverse(preconditioner, matrix_c0, matrix_h, matrix_s, 0.20_dp)
2697 CALL make_complex_full_single_inverse(preconditioner_rot, matrix_c0_rot, matrix_h, matrix_s, 0.20_dp)
2698 CALL make_complex_full_s_inverse(preconditioner_s, matrix_s)
2699 CALL make_complex_full_single(preconditioner_single, matrix_h, matrix_s, 0.15_dp, 0.20_dp)
2700 CALL make_complex_full_kinetic(preconditioner_t, matrix_t, matrix_s, 0.20_dp)
2701 CALL make_complex_full_all(full_all, matrix_c0, matrix_h, matrix_s, c0_evals, 0.08_dp)
2702 CALL make_complex_full_all(full_all_phase, matrix_c0_phase, matrix_h, matrix_s, &
2703 c0_evals, 0.08_dp)
2704 CALL apply_preconditioner(full_all, matrix_gradient, matrix_output)
2705 CALL apply_preconditioner(full_all_phase, matrix_gradient_phase, matrix_output_phase)
2706 DO j = 1, n
2707 DO i = 1, n
2708 CALL cp_cfm_get_element(preconditioner%complex_fm, i, j, value)
2709 inverse(i, j) = value
2710 CALL cp_cfm_get_element(preconditioner_rot%complex_fm, i, j, value)
2711 inverse_rot(i, j) = value
2712 CALL cp_cfm_get_element(preconditioner_s%complex_fm, i, j, value)
2713 inverse_s(i, j) = value
2714 CALL cp_cfm_get_element(preconditioner_single%complex_fm, i, j, value)
2715 inverse_single(i, j) = value
2716 CALL cp_cfm_get_element(preconditioner_t%complex_fm, i, j, value)
2717 inverse_t(i, j) = value
2718 END DO
2719 END DO
2720 gauge_error = maxval(abs(inverse - inverse_rot))
2721 hermitian_error = maxval(abs(inverse - conjg(transpose(inverse))))
2722 CALL diag_complex(inverse, eigenvectors, eigenvalues)
2723 full_single_hermitian_error = maxval(abs(inverse_single - conjg(transpose(inverse_single))))
2724 CALL diag_complex(inverse_single, eigenvectors, eigenvalues_single)
2725 CALL dense_hermitian_inverse(s, reference_s)
2726 CALL dense_hermitian_inverse(t + 0.20_dp*s, reference_t)
2727 s_error = maxval(abs(inverse_s - reference_s))
2728 t_error = maxval(abs(inverse_t - reference_t))
2729 DO j = 1, k
2730 DO i = 1, n
2731 CALL cp_cfm_get_element(matrix_output, i, j, value)
2732 output(i, j) = value
2733 CALL cp_cfm_get_element(matrix_output_phase, i, j, value)
2734 output_phase(i, j) = value
2735 END DO
2736 END DO
2737 DO j = 1, n
2738 DO i = 1, n
2739 CALL cp_cfm_get_element(full_all%complex_fm, i, j, value)
2740 full_all_vectors(i, j) = value
2741 END DO
2742 END DO
2743 full_all_gauge_error = 0.0_dp
2744 DO j = 1, k
2745 full_all_gauge_error = max(full_all_gauge_error, &
2746 maxval(abs(output_phase(:, j) - phases(j)*output(:, j))))
2747 END DO
2748 spectral_metric = matmul(conjg(transpose(full_all_vectors)), &
2749 matmul(s, full_all_vectors))
2750 full_all_orthogonality_error = maxval(abs(spectral_metric - &
2751 reshape([cmplx(1.0_dp, 0.0_dp, kind=dp), &
2752 cmplx(0.0_dp, 0.0_dp, kind=dp), &
2753 cmplx(0.0_dp, 0.0_dp, kind=dp), &
2754 cmplx(0.0_dp, 0.0_dp, kind=dp), &
2755 cmplx(1.0_dp, 0.0_dp, kind=dp), &
2756 cmplx(0.0_dp, 0.0_dp, kind=dp), &
2757 cmplx(0.0_dp, 0.0_dp, kind=dp), &
2758 cmplx(0.0_dp, 0.0_dp, kind=dp), &
2759 cmplx(1.0_dp, 0.0_dp, kind=dp)], [n, n])))
2760 full_all_descent = real(sum(conjg(gradient)*output), kind=dp)
2761
2762 CALL make_complex_fermi_low_rank(low_rank, matrix_c0, matrix_h, matrix_s, &
2763 0.08_dp, 1, 1.0_dp, 1.0e-8_dp)
2764 CALL make_complex_fermi_low_rank(low_rank_rot, matrix_c0_rot, matrix_h, matrix_s, &
2765 0.08_dp, 1, 1.0_dp, 1.0e-8_dp)
2766 CALL apply_preconditioner(low_rank, matrix_gradient, matrix_output)
2767 CALL apply_preconditioner(low_rank_rot, matrix_gradient_rot, matrix_output_rot)
2768 DO j = 1, k
2769 DO i = 1, n
2770 CALL cp_cfm_get_element(matrix_output, i, j, value)
2771 output(i, j) = value
2772 CALL cp_cfm_get_element(matrix_output_rot, i, j, value)
2773 output_rot(i, j) = value
2774 END DO
2775 END DO
2776 low_rank_gauge_error = maxval(abs(output_rot - matmul(output, c0_rot(1:k, :))))
2777 low_rank_descent = real(sum(conjg(gradient)*output), kind=dp)
2778
2779 CALL make_complex_full_all_covariant(covariant, matrix_c0, matrix_h, matrix_s, 0.08_dp)
2780 CALL make_complex_full_all_covariant(covariant_rot, matrix_c0_rot, matrix_h, matrix_s, 0.08_dp)
2781 CALL apply_preconditioner(covariant, matrix_gradient, matrix_output)
2782 CALL apply_preconditioner(covariant_rot, matrix_gradient_rot, matrix_output_rot)
2783 DO j = 1, k
2784 DO i = 1, n
2785 CALL cp_cfm_get_element(matrix_output, i, j, value)
2786 output(i, j) = value
2787 CALL cp_cfm_get_element(matrix_output_rot, i, j, value)
2788 output_rot(i, j) = value
2789 END DO
2790 END DO
2791 covariant_gauge_error = maxval(abs(output_rot - matmul(output, c0_rot(1:k, :))))
2792 covariant_descent = real(sum(conjg(gradient)*output), kind=dp)
2793 IF (gauge_error > 1.0e-11_dp .OR. hermitian_error > 1.0e-12_dp .OR. &
2794 minval(eigenvalues) <= 0.0_dp) nfail = nfail + 1
2795 IF (s_error > 1.0e-11_dp .OR. t_error > 1.0e-11_dp) nfail = nfail + 1
2796 IF (full_single_hermitian_error > 1.0e-12_dp .OR. &
2797 minval(eigenvalues_single) <= 0.0_dp) nfail = nfail + 1
2798 IF (full_all_gauge_error > 1.0e-11_dp .OR. &
2799 full_all_orthogonality_error > 1.0e-11_dp .OR. full_all_descent <= 0.0_dp) nfail = nfail + 1
2800 IF (low_rank_gauge_error > 1.0e-11_dp .OR. low_rank_descent <= 0.0_dp) nfail = nfail + 1
2801 IF (covariant_gauge_error > 1.0e-11_dp .OR. covariant_descent <= 0.0_dp) nfail = nfail + 1
2802 IF (para_env%is_source()) THEN
2803 WRITE (io_unit, '(A,3(1X,ES13.6))') 'complex FULL_SINGLE gauge/hermitian/min-eigenvalue', &
2804 gauge_error, hermitian_error, minval(eigenvalues)
2805 WRITE (io_unit, '(A,2(1X,ES13.6))') 'complex S/KINETIC inverse errors', s_error, t_error
2806 WRITE (io_unit, '(A,2(1X,ES13.6))') 'complex FULL_SINGLE hermitian/min-eigenvalue', &
2807 full_single_hermitian_error, minval(eigenvalues_single)
2808 WRITE (io_unit, '(A,3(1X,ES13.6))') 'complex FULL_ALL phase/orthogonality/descent', &
2809 full_all_gauge_error, full_all_orthogonality_error, full_all_descent
2810 WRITE (io_unit, '(A,2(1X,ES13.6))') 'complex FERMI_LOW_RANK covariance/descent', &
2811 low_rank_gauge_error, low_rank_descent
2812 WRITE (io_unit, '(A,2(1X,ES13.6))') 'complex FULL_ALL_COVARIANT covariance/descent', &
2813 covariant_gauge_error, covariant_descent
2814 END IF
2815
2816 CALL destroy_preconditioner(covariant_rot)
2817 CALL destroy_preconditioner(covariant)
2818 CALL destroy_preconditioner(low_rank_rot)
2819 CALL destroy_preconditioner(low_rank)
2820 CALL destroy_preconditioner(full_all_phase)
2821 CALL destroy_preconditioner(full_all)
2822 CALL destroy_preconditioner(preconditioner_t)
2823 CALL destroy_preconditioner(preconditioner_single)
2824 CALL destroy_preconditioner(preconditioner_s)
2825 CALL destroy_preconditioner(preconditioner_rot)
2826 CALL destroy_preconditioner(preconditioner)
2827 CALL cp_cfm_release(matrix_c0_rot)
2828 CALL cp_cfm_release(matrix_c0_phase)
2829 CALL cp_cfm_release(matrix_c0)
2830 CALL cp_cfm_release(matrix_output_phase)
2831 CALL cp_cfm_release(matrix_output)
2832 CALL cp_cfm_release(matrix_output_rot)
2833 CALL cp_cfm_release(matrix_gradient_rot)
2834 CALL cp_cfm_release(matrix_gradient_phase)
2835 CALL cp_cfm_release(matrix_gradient)
2836 CALL cp_cfm_release(matrix_t)
2837 CALL cp_cfm_release(matrix_s)
2838 CALL cp_cfm_release(matrix_h)
2839 CALL cp_fm_struct_release(fmstruct_nk)
2840 CALL cp_fm_struct_release(fmstruct_nn)
2841 CALL cp_blacs_env_release(blacs_env)
2842
2843 END SUBROUTINE test_complex_preconditioner_gauge
2844
2845! **************************************************************************************************
2846!> \brief Check complex covariant preconditioning within equal-occupation ROKS blocks.
2847!> \param para_env parallel environment
2848!> \param nfail accumulated number of failures
2849! **************************************************************************************************
2850 SUBROUTINE test_complex_occupation_block_preconditioner(para_env, nfail)
2851 TYPE(mp_para_env_type), POINTER :: para_env
2852 INTEGER, INTENT(INOUT) :: nfail
2853
2854 INTEGER, PARAMETER :: nbas = 4, norb = 3
2855
2856 COMPLEX(KIND=dp) :: value
2857 COMPLEX(KIND=dp), DIMENSION(nbas, nbas) :: hamiltonian, overlap
2858 COMPLEX(KIND=dp), DIMENSION(nbas, norb) :: coeff, coeff_rot, gradient, &
2859 gradient_rot, output, output_rot
2860 COMPLEX(KIND=dp), DIMENSION(norb, norb) :: block_rotation, private_rotation
2861 INTEGER :: i, j
2862 REAL(KIND=dp) :: covariance_error, cross_block_error, &
2863 descent, shifted_eigenvalue_error
2864 REAL(KIND=dp), DIMENSION(norb, 2) :: occupation_signature
2865 TYPE(cp_blacs_env_type), POINTER :: blacs_env
2866 TYPE(cp_cfm_type) :: matrix_c, matrix_c_rot, matrix_gradient, matrix_gradient_rot, matrix_h, &
2867 matrix_output, matrix_output_rot, matrix_s
2868 TYPE(cp_fm_struct_type), POINTER :: fmstruct_nk, fmstruct_nn
2869 TYPE(preconditioner_type) :: preconditioner, preconditioner_rot
2870
2871 NULLIFY (blacs_env, fmstruct_nk, fmstruct_nn)
2872 CALL cp_blacs_env_create(blacs_env=blacs_env, para_env=para_env)
2873 CALL cp_fm_struct_create(fmstruct=fmstruct_nn, para_env=para_env, context=blacs_env, &
2874 nrow_global=nbas, ncol_global=nbas)
2875 CALL cp_fm_struct_create(fmstruct=fmstruct_nk, para_env=para_env, context=blacs_env, &
2876 nrow_global=nbas, ncol_global=norb)
2877 CALL cp_cfm_create(matrix_h, fmstruct_nn, name='occupation-block H')
2878 CALL cp_cfm_create(matrix_s, fmstruct_nn, name='occupation-block S')
2879 CALL cp_cfm_create(matrix_c, fmstruct_nk, name='occupation-block C')
2880 CALL cp_cfm_create(matrix_c_rot, fmstruct_nk, name='occupation-block C U')
2881 CALL cp_cfm_create(matrix_gradient, fmstruct_nk, name='occupation-block gradient')
2882 CALL cp_cfm_create(matrix_gradient_rot, fmstruct_nk, &
2883 name='occupation-block rotated gradient')
2884 CALL cp_cfm_create(matrix_output, fmstruct_nk, name='occupation-block output')
2885 CALL cp_cfm_create(matrix_output_rot, fmstruct_nk, name='occupation-block rotated output')
2886
2887 hamiltonian(:, :) = cmplx(0.0_dp, 0.0_dp, kind=dp)
2888 hamiltonian(1, 1) = cmplx(-0.42_dp, 0.0_dp, kind=dp)
2889 hamiltonian(2, 2) = cmplx(-0.18_dp, 0.0_dp, kind=dp)
2890 hamiltonian(3, 3) = cmplx(0.07_dp, 0.0_dp, kind=dp)
2891 hamiltonian(4, 4) = cmplx(0.83_dp, 0.0_dp, kind=dp)
2892 hamiltonian(1, 2) = cmplx(0.09_dp, 0.04_dp, kind=dp)
2893 hamiltonian(2, 1) = conjg(hamiltonian(1, 2))
2894 hamiltonian(1, 3) = cmplx(-0.06_dp, 0.03_dp, kind=dp)
2895 hamiltonian(3, 1) = conjg(hamiltonian(1, 3))
2896 hamiltonian(2, 3) = cmplx(0.05_dp, -0.02_dp, kind=dp)
2897 hamiltonian(3, 2) = conjg(hamiltonian(2, 3))
2898 hamiltonian(1, 4) = cmplx(0.04_dp, -0.03_dp, kind=dp)
2899 hamiltonian(4, 1) = conjg(hamiltonian(1, 4))
2900 hamiltonian(2, 4) = cmplx(-0.02_dp, 0.05_dp, kind=dp)
2901 hamiltonian(4, 2) = conjg(hamiltonian(2, 4))
2902 hamiltonian(3, 4) = cmplx(0.07_dp, 0.01_dp, kind=dp)
2903 hamiltonian(4, 3) = conjg(hamiltonian(3, 4))
2904 overlap(:, :) = cmplx(0.0_dp, 0.0_dp, kind=dp)
2905 DO i = 1, nbas
2906 overlap(i, i) = cmplx(1.0_dp, 0.0_dp, kind=dp)
2907 END DO
2908 coeff(:, :) = cmplx(0.0_dp, 0.0_dp, kind=dp)
2909 DO i = 1, norb
2910 coeff(i, i) = cmplx(1.0_dp, 0.0_dp, kind=dp)
2911 END DO
2912 block_rotation(:, :) = cmplx(0.0_dp, 0.0_dp, kind=dp)
2913 block_rotation(1, 1) = cmplx(sqrt(0.7_dp), 0.0_dp, kind=dp)
2914 block_rotation(1, 3) = cmplx(0.0_dp, sqrt(0.3_dp), kind=dp)
2915 block_rotation(2, 2) = exp(cmplx(0.0_dp, 0.37_dp, kind=dp))
2916 block_rotation(3, 1) = cmplx(0.0_dp, sqrt(0.3_dp), kind=dp)
2917 block_rotation(3, 3) = cmplx(sqrt(0.7_dp), 0.0_dp, kind=dp)
2918 coeff_rot = matmul(coeff, block_rotation)
2919 gradient(:, 1) = [cmplx(0.12_dp, 0.03_dp, kind=dp), &
2920 cmplx(-0.04_dp, 0.07_dp, kind=dp), &
2921 cmplx(0.08_dp, -0.02_dp, kind=dp), &
2922 cmplx(-0.05_dp, 0.06_dp, kind=dp)]
2923 gradient(:, 2) = [cmplx(-0.06_dp, 0.02_dp, kind=dp), &
2924 cmplx(0.11_dp, -0.05_dp, kind=dp), &
2925 cmplx(0.03_dp, 0.09_dp, kind=dp), &
2926 cmplx(0.07_dp, -0.01_dp, kind=dp)]
2927 gradient(:, 3) = [cmplx(0.04_dp, -0.08_dp, kind=dp), &
2928 cmplx(-0.02_dp, 0.05_dp, kind=dp), &
2929 cmplx(0.10_dp, 0.01_dp, kind=dp), &
2930 cmplx(-0.03_dp, -0.06_dp, kind=dp)]
2931 gradient_rot = matmul(gradient, block_rotation)
2932 occupation_signature(:, 1) = [1.0_dp, 1.0_dp, 1.0_dp]
2933 occupation_signature(:, 2) = [1.0_dp, 0.0_dp, 1.0_dp]
2934
2935 DO j = 1, nbas
2936 DO i = 1, nbas
2937 CALL cp_cfm_set_element(matrix_h, i, j, hamiltonian(i, j))
2938 CALL cp_cfm_set_element(matrix_s, i, j, overlap(i, j))
2939 END DO
2940 END DO
2941 DO j = 1, norb
2942 DO i = 1, nbas
2943 CALL cp_cfm_set_element(matrix_c, i, j, coeff(i, j))
2944 CALL cp_cfm_set_element(matrix_c_rot, i, j, coeff_rot(i, j))
2945 CALL cp_cfm_set_element(matrix_gradient, i, j, gradient(i, j))
2946 CALL cp_cfm_set_element(matrix_gradient_rot, i, j, gradient_rot(i, j))
2947 END DO
2948 END DO
2949
2950 CALL init_preconditioner(preconditioner, para_env, blacs_env)
2951 CALL init_preconditioner(preconditioner_rot, para_env, blacs_env)
2952 CALL make_complex_full_all_covariant(preconditioner, matrix_c, matrix_h, matrix_s, &
2953 0.08_dp, occupation_signature)
2954 CALL make_complex_full_all_covariant(preconditioner_rot, matrix_c_rot, matrix_h, matrix_s, &
2955 0.08_dp, occupation_signature)
2956 CALL apply_preconditioner(preconditioner, matrix_gradient, matrix_output)
2957 CALL apply_preconditioner(preconditioner_rot, matrix_gradient_rot, matrix_output_rot)
2958 DO j = 1, norb
2959 DO i = 1, nbas
2960 CALL cp_cfm_get_element(matrix_output, i, j, value)
2961 output(i, j) = value
2962 CALL cp_cfm_get_element(matrix_output_rot, i, j, value)
2963 output_rot(i, j) = value
2964 END DO
2965 END DO
2966 DO j = 1, norb
2967 DO i = 1, norb
2968 CALL cp_cfm_get_element(preconditioner%occ_rotation_complex, i, j, value)
2969 private_rotation(i, j) = value
2970 END DO
2971 END DO
2972
2973 covariance_error = maxval(abs(output_rot - matmul(output, block_rotation)))
2974 cross_block_error = max(maxval(abs(private_rotation([1, 3], 3))), &
2975 maxval(abs(private_rotation(2, 1:2))))
2976 shifted_eigenvalue_error = maxval(abs(preconditioner%occ_evals))
2977 descent = real(sum(conjg(gradient)*output), kind=dp)
2978 IF (covariance_error > 1.0e-11_dp .OR. cross_block_error > 1.0e-12_dp .OR. &
2979 shifted_eigenvalue_error > 1.0_dp .OR. descent <= 0.0_dp) nfail = nfail + 1
2980 IF (para_env%is_source()) THEN
2981 WRITE (io_unit, '(A,4(1X,ES13.6))') &
2982 'complex occupation-block covariance/cross-block/max-eigenvalue/descent', &
2983 covariance_error, cross_block_error, shifted_eigenvalue_error, descent
2984 END IF
2985
2986 CALL destroy_preconditioner(preconditioner_rot)
2987 CALL destroy_preconditioner(preconditioner)
2988 CALL cp_cfm_release(matrix_output_rot)
2989 CALL cp_cfm_release(matrix_output)
2990 CALL cp_cfm_release(matrix_gradient_rot)
2991 CALL cp_cfm_release(matrix_gradient)
2992 CALL cp_cfm_release(matrix_c_rot)
2993 CALL cp_cfm_release(matrix_c)
2994 CALL cp_cfm_release(matrix_s)
2995 CALL cp_cfm_release(matrix_h)
2996 CALL cp_fm_struct_release(fmstruct_nk)
2997 CALL cp_fm_struct_release(fmstruct_nn)
2998 CALL cp_blacs_env_release(blacs_env)
2999
3000 END SUBROUTINE test_complex_occupation_block_preconditioner
3001
3002! **************************************************************************************************
3003!> \brief Check the intentionally narrow complex k-point preconditioner support matrix.
3004!> \param nfail accumulated number of failures
3005! **************************************************************************************************
3006 SUBROUTINE test_preconditioner_support(nfail)
3007 INTEGER, INTENT(INOUT) :: nfail
3008
3009 IF (.NOT. qs_ot_kpoint_preconditioner_supported(ot_precond_none, .false.)) nfail = nfail + 1
3010 IF (.NOT. qs_ot_kpoint_preconditioner_supported(ot_precond_fermi_low_rank, .false.)) nfail = nfail + 1
3011 IF (.NOT. qs_ot_kpoint_preconditioner_supported(ot_precond_full_single_inverse, .false.)) THEN
3012 nfail = nfail + 1
3013 END IF
3014 IF (.NOT. qs_ot_kpoint_preconditioner_supported(ot_precond_full_all, .false.)) nfail = nfail + 1
3015 IF (.NOT. qs_ot_kpoint_preconditioner_supported(ot_precond_full_single, .false.)) nfail = nfail + 1
3016 IF (.NOT. qs_ot_kpoint_preconditioner_supported( &
3017 ot_precond_full_all_covariant, .false.)) nfail = nfail + 1
3018 IF (.NOT. qs_ot_kpoint_preconditioner_supported(ot_precond_full_kinetic, .false.)) nfail = nfail + 1
3019 IF (.NOT. qs_ot_kpoint_preconditioner_supported(ot_precond_s_inverse, .false.)) nfail = nfail + 1
3020 IF (qs_ot_kpoint_preconditioner_supported(ot_precond_full_single_inverse, .true.)) nfail = nfail + 1
3021
3022 IF (.NOT. qs_ot_kpoint_preconditioner_solver_supported( &
3023 ot_precond_full_single, ot_precond_solver_default, .false.)) nfail = nfail + 1
3024 IF (qs_ot_kpoint_preconditioner_solver_supported( &
3025 ot_precond_full_single, ot_precond_solver_inv_chol, .false.)) nfail = nfail + 1
3026
3027 IF (.NOT. qs_ot_kpoint_preconditioner_solver_supported( &
3028 ot_precond_full_single_inverse, ot_precond_solver_default, .false.)) nfail = nfail + 1
3029 IF (.NOT. qs_ot_kpoint_preconditioner_solver_supported( &
3030 ot_precond_full_single_inverse, ot_precond_solver_inv_chol, .false.)) nfail = nfail + 1
3031 IF (qs_ot_kpoint_preconditioner_solver_supported( &
3032 ot_precond_full_single_inverse, ot_precond_solver_direct, .false.)) nfail = nfail + 1
3033 IF (qs_ot_kpoint_preconditioner_solver_supported( &
3034 ot_precond_full_single_inverse, ot_precond_solver_chebyshev, .false.)) nfail = nfail + 1
3035 IF (.NOT. qs_ot_kpoint_preconditioner_solver_supported( &
3036 ot_precond_full_all, ot_precond_solver_default, .false.)) nfail = nfail + 1
3037 IF (qs_ot_kpoint_preconditioner_solver_supported( &
3038 ot_precond_full_all, ot_precond_solver_inv_chol, .false.)) nfail = nfail + 1
3039 IF (.NOT. qs_ot_kpoint_preconditioner_solver_supported( &
3040 ot_precond_fermi_low_rank, ot_precond_solver_default, .false.)) nfail = nfail + 1
3041 IF (qs_ot_kpoint_preconditioner_solver_supported( &
3042 ot_precond_fermi_low_rank, ot_precond_solver_inv_chol, .false.)) nfail = nfail + 1
3043 IF (.NOT. qs_ot_kpoint_preconditioner_solver_supported( &
3044 ot_precond_full_all_covariant, ot_precond_solver_default, .false.)) nfail = nfail + 1
3045 IF (qs_ot_kpoint_preconditioner_solver_supported( &
3046 ot_precond_full_all_covariant, ot_precond_solver_inv_chol, .false.)) nfail = nfail + 1
3047 IF (.NOT. qs_ot_kpoint_preconditioner_solver_supported( &
3048 ot_precond_full_kinetic, ot_precond_solver_default, .false.)) nfail = nfail + 1
3049 IF (.NOT. qs_ot_kpoint_preconditioner_solver_supported( &
3050 ot_precond_full_kinetic, ot_precond_solver_inv_chol, .false.)) nfail = nfail + 1
3051 IF (qs_ot_kpoint_preconditioner_solver_supported( &
3052 ot_precond_full_kinetic, ot_precond_solver_direct, .false.)) nfail = nfail + 1
3053 IF (.NOT. qs_ot_kpoint_preconditioner_solver_supported( &
3054 ot_precond_s_inverse, ot_precond_solver_default, .false.)) nfail = nfail + 1
3055 IF (.NOT. qs_ot_kpoint_preconditioner_solver_supported( &
3056 ot_precond_s_inverse, ot_precond_solver_inv_chol, .false.)) nfail = nfail + 1
3057 IF (qs_ot_kpoint_preconditioner_solver_supported( &
3058 ot_precond_s_inverse, ot_precond_solver_chebyshev, .false.)) nfail = nfail + 1
3059 IF (abs(qs_ot_kpoint_preconditioner_scale(0.25_dp) - 4.0_dp) > epsilon(1.0_dp)) THEN
3060 nfail = nfail + 1
3061 END IF
3062
3063 END SUBROUTINE test_preconditioner_support
3064
3065! **************************************************************************************************
3066!> \brief Create and reserve a one-block DBCSR matrix.
3067!> \param matrix ...
3068!> \param name ...
3069!> \param dist ...
3070!> \param row_size ...
3071!> \param col_size ...
3072! **************************************************************************************************
3073 SUBROUTINE create_one_block(matrix, name, dist, row_size, col_size)
3074 TYPE(dbcsr_type), INTENT(INOUT) :: matrix
3075 CHARACTER(LEN=*), INTENT(IN) :: name
3076 TYPE(dbcsr_distribution_type), INTENT(IN) :: dist
3077 INTEGER, DIMENSION(:), INTENT(INOUT), POINTER :: row_size, col_size
3078
3079 CALL dbcsr_create(matrix, name, dist, dbcsr_type_no_symmetry, row_size, col_size)
3080 CALL dbcsr_reserve_blocks(matrix, rows=[1], cols=[1])
3081 CALL dbcsr_finalize(matrix)
3082
3083 END SUBROUTINE create_one_block
3084
3085! **************************************************************************************************
3086!> \brief Store a dense complex matrix in a real/imaginary DBCSR pair.
3087!> \param matrix_re ...
3088!> \param matrix_im ...
3089!> \param matrix ...
3090! **************************************************************************************************
3091 SUBROUTINE put_complex_pair(matrix_re, matrix_im, matrix)
3092 TYPE(dbcsr_type), INTENT(INOUT) :: matrix_re, matrix_im
3093 COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(IN) :: matrix
3094
3095 CALL dbcsr_put_block(matrix_re, 1, 1, real(matrix, kind=dp))
3096 CALL dbcsr_put_block(matrix_im, 1, 1, aimag(matrix))
3097
3098 END SUBROUTINE put_complex_pair
3099
3100! **************************************************************************************************
3101!> \brief Retrieve a dense complex matrix from a one-block DBCSR pair.
3102!> \param matrix_re ...
3103!> \param matrix_im ...
3104!> \param matrix ...
3105! **************************************************************************************************
3106 SUBROUTINE get_complex_pair(matrix_re, matrix_im, matrix)
3107 TYPE(dbcsr_type), INTENT(INOUT) :: matrix_re, matrix_im
3108 COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: matrix
3109
3110 LOGICAL :: found
3111 REAL(KIND=dp), DIMENSION(:, :), POINTER :: block
3112
3113 matrix(:, :) = cmplx(0.0_dp, 0.0_dp, kind=dp)
3114 CALL dbcsr_get_block_p(matrix_re, 1, 1, block, found)
3115 IF (found) matrix = matrix + cmplx(block, 0.0_dp, kind=dp)
3116 CALL dbcsr_get_block_p(matrix_im, 1, 1, block, found)
3117 IF (found) matrix = matrix + cmplx(0.0_dp, block, kind=dp)
3118
3119 END SUBROUTINE get_complex_pair
3120
3121! **************************************************************************************************
3122!> \brief Dense fixed-occupation energy after finite REF orthogonalization and rotation.
3123!> \param x REF coordinate
3124!> \param generator anti-Hermitian rotation generator
3125!> \param hamiltonian Hermitian one-particle Hamiltonian
3126!> \param occupation orbital occupations
3127!> \return energy
3128! **************************************************************************************************
3129 FUNCTION rotated_polar_energy(x, generator, hamiltonian, occupation) RESULT(energy)
3130 COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(IN) :: x, generator, hamiltonian
3131 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: occupation
3132 REAL(KIND=dp) :: energy
3133
3134 COMPLEX(KIND=dp), &
3135 DIMENSION(SIZE(x, 1), SIZE(x, 2)) :: orbitals
3136 COMPLEX(KIND=dp), &
3137 DIMENSION(SIZE(x, 2), SIZE(x, 2)) :: inverse_sqrt, rotation
3138 INTEGER :: i
3139
3140 CALL dense_inverse_sqrt(matmul(conjg(transpose(x)), x), inverse_sqrt)
3141 rotation = dense_antihermitian_exp(generator)
3142 orbitals = matmul(matmul(x, inverse_sqrt), rotation)
3143 energy = 0.0_dp
3144 DO i = 1, SIZE(occupation)
3145 energy = energy + occupation(i)* &
3146 REAL(DOT_PRODUCT(orbitals(:, i), MATMUL(hamiltonian, orbitals(:, i))), KIND=dp)
3147 END DO
3148
3149 END FUNCTION rotated_polar_energy
3150
3151! **************************************************************************************************
3152!> \brief Dense reference exponential for an anti-Hermitian matrix.
3153!> \param generator anti-Hermitian generator
3154!> \return its unitary exponential
3155! **************************************************************************************************
3156 FUNCTION dense_antihermitian_exp(generator) RESULT(rotation)
3157 COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(IN) :: generator
3158 COMPLEX(KIND=dp), DIMENSION(SIZE(generator, 1), & SIZE(generator, 2)) :: rotation
3159
3160 COMPLEX(KIND=dp), DIMENSION(SIZE(generator, 1), & SIZE(generator, 2)) :: hermitian, vectors
3161 INTEGER :: i
3162 REAL(KIND=dp), DIMENSION(SIZE(generator, 1)) :: eigenvalues
3163
3164 hermitian = cmplx(0.0_dp, 1.0_dp, kind=dp)*generator
3165 CALL diag_complex(hermitian, vectors, eigenvalues)
3166 rotation(:, :) = cmplx(0.0_dp, 0.0_dp, kind=dp)
3167 DO i = 1, SIZE(eigenvalues)
3168 rotation = rotation + exp(cmplx(0.0_dp, -eigenvalues(i), kind=dp))* &
3169 matmul(reshape(vectors(:, i), [SIZE(generator, 1), 1]), &
3170 reshape(conjg(vectors(:, i)), [1, SIZE(generator, 1)]))
3171 END DO
3172
3173 END FUNCTION dense_antihermitian_exp
3174
3175! **************************************************************************************************
3176!> \brief Hermitian inverse for an independent dense reference.
3177!> \param matrix ...
3178!> \param inverse ...
3179! **************************************************************************************************
3180 SUBROUTINE dense_hermitian_inverse(matrix, inverse)
3181 COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(IN) :: matrix
3182 COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: inverse
3183
3184 COMPLEX(KIND=dp), &
3185 DIMENSION(SIZE(matrix, 1), SIZE(matrix, 2)) :: vectors
3186 INTEGER :: i
3187 REAL(KIND=dp), DIMENSION(SIZE(matrix, 1)) :: eigenvalues
3188
3189 CALL diag_complex(matrix, vectors, eigenvalues)
3190 inverse(:, :) = cmplx(0.0_dp, 0.0_dp, kind=dp)
3191 DO i = 1, SIZE(eigenvalues)
3192 inverse = inverse + matmul(reshape(vectors(:, i), [SIZE(matrix, 1), 1]), &
3193 reshape(conjg(vectors(:, i)), [1, SIZE(matrix, 1)]))/ &
3194 eigenvalues(i)
3195 END DO
3196
3197 END SUBROUTINE dense_hermitian_inverse
3198
3199! **************************************************************************************************
3200!> \brief Hermitian inverse square root for the independent dense reference.
3201!> \param matrix ...
3202!> \param inverse ...
3203! **************************************************************************************************
3204 SUBROUTINE dense_inverse_sqrt(matrix, inverse)
3205 COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(IN) :: matrix
3206 COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: inverse
3207
3208 COMPLEX(KIND=dp), &
3209 DIMENSION(SIZE(matrix, 1), SIZE(matrix, 2)) :: vectors
3210 INTEGER :: i
3211 REAL(KIND=dp), DIMENSION(SIZE(matrix, 1)) :: eigenvalues
3212
3213 CALL diag_complex(matrix, vectors, eigenvalues)
3214 inverse(:, :) = cmplx(0.0_dp, 0.0_dp, kind=dp)
3215 DO i = 1, SIZE(eigenvalues)
3216 inverse = inverse + matmul(reshape(vectors(:, i), [SIZE(matrix, 1), 1]), &
3217 reshape(conjg(vectors(:, i)), [1, SIZE(matrix, 1)]))/ &
3218 sqrt(eigenvalues(i))
3219 END DO
3220
3221 END SUBROUTINE dense_inverse_sqrt
3222
3223! **************************************************************************************************
3224!> \brief Energy of the polar-orthonormalized occupied subspace.
3225!> \param coordinate ...
3226!> \param hamiltonian ...
3227!> \return ...
3228! **************************************************************************************************
3229 FUNCTION polar_energy(coordinate, hamiltonian) RESULT(energy)
3230 COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(IN) :: coordinate, hamiltonian
3231 REAL(KIND=dp) :: energy
3232
3233 COMPLEX(KIND=dp), DIMENSION(SIZE(coordinate, 1), & SIZE(coordinate, 2)) :: orbitals
3234 COMPLEX(KIND=dp), DIMENSION(SIZE(coordinate, 2), & SIZE(coordinate, 2)) :: inverse
3235
3236 CALL dense_inverse_sqrt(matmul(conjg(transpose(coordinate)), coordinate), inverse)
3237 orbitals = matmul(coordinate, inverse)
3238 energy = real(sum(conjg(orbitals)*matmul(hamiltonian, orbitals)), kind=dp)
3239
3240 END FUNCTION polar_energy
3241
3242! **************************************************************************************************
3243!> \brief Independent dense finite STRICT transformation for a Hermitian overlap.
3244!> \param c0 reference orbitals
3245!> \param coordinate tangent coordinate
3246!> \param overlap Hermitian overlap matrix
3247!> \param orbitals transformed orbitals
3248! **************************************************************************************************
3249 SUBROUTINE strict_orbitals(c0, coordinate, overlap, orbitals)
3250 COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(IN) :: c0, coordinate, overlap
3251 COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: orbitals
3252
3253 COMPLEX(KIND=dp), DIMENSION(SIZE(coordinate, 2), & SIZE(coordinate, 2)) :: cosp, sinp, vectors
3254 INTEGER :: i
3255 REAL(KIND=dp) :: root
3256 REAL(KIND=dp), DIMENSION(SIZE(coordinate, 2)) :: eigenvalues
3257
3258 CALL diag_complex(matmul(conjg(transpose(coordinate)), matmul(overlap, coordinate)), &
3259 vectors, eigenvalues)
3260 cosp(:, :) = cmplx(0.0_dp, 0.0_dp, kind=dp)
3261 sinp(:, :) = cmplx(0.0_dp, 0.0_dp, kind=dp)
3262 DO i = 1, SIZE(eigenvalues)
3263 root = sqrt(max(0.0_dp, eigenvalues(i)))
3264 cosp = cosp + cos(root)*matmul(reshape(vectors(:, i), [SIZE(vectors, 1), 1]), &
3265 reshape(conjg(vectors(:, i)), [1, SIZE(vectors, 1)]))
3266 IF (root < 1.0e-12_dp) THEN
3267 sinp = sinp + matmul(reshape(vectors(:, i), [SIZE(vectors, 1), 1]), &
3268 reshape(conjg(vectors(:, i)), [1, SIZE(vectors, 1)]))
3269 ELSE
3270 sinp = sinp + sin(root)/root* &
3271 matmul(reshape(vectors(:, i), [SIZE(vectors, 1), 1]), &
3272 reshape(conjg(vectors(:, i)), [1, SIZE(vectors, 1)]))
3273 END IF
3274 END DO
3275 orbitals = matmul(c0, cosp) + matmul(coordinate, sinp)
3276
3277 END SUBROUTINE strict_orbitals
3278
3279! **************************************************************************************************
3280!> \brief Energy of the independently transformed finite STRICT orbitals.
3281!> \param c0 reference orbitals
3282!> \param coordinate tangent coordinate
3283!> \param overlap Hermitian overlap matrix
3284!> \param hamiltonian Hermitian Hamiltonian
3285!> \return energy
3286! **************************************************************************************************
3287 FUNCTION strict_energy(c0, coordinate, overlap, hamiltonian) RESULT(energy)
3288 COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(IN) :: c0, coordinate, overlap, hamiltonian
3289 REAL(KIND=dp) :: energy
3290
3291 COMPLEX(KIND=dp), DIMENSION(SIZE(coordinate, 1), & SIZE(coordinate, 2)) :: orbitals
3292
3293 CALL strict_orbitals(c0, coordinate, overlap, orbitals)
3294 energy = real(sum(conjg(orbitals)*matmul(hamiltonian, orbitals)), kind=dp)
3295
3296 END FUNCTION strict_energy
3297
3298! **************************************************************************************************
3299!> \brief Energy of finite STRICT orbitals followed by a complex occupied-space rotation.
3300!> \param c0 reference orbitals
3301!> \param coordinate tangent coordinate
3302!> \param generator anti-Hermitian rotation generator
3303!> \param overlap Hermitian overlap matrix
3304!> \param hamiltonian Hermitian Hamiltonian
3305!> \param occupation orbital occupations
3306!> \return energy
3307! **************************************************************************************************
3308 FUNCTION strict_rotated_energy(c0, coordinate, generator, overlap, hamiltonian, &
3309 occupation) RESULT(energy)
3310 COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(IN) :: c0, coordinate, generator, overlap, &
3311 hamiltonian
3312 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: occupation
3313 REAL(KIND=dp) :: energy
3314
3315 COMPLEX(KIND=dp), DIMENSION(SIZE(coordinate, 1), & SIZE(coordinate, 2)) :: orbitals
3316 INTEGER :: i
3317
3318 CALL strict_orbitals(c0, coordinate, overlap, orbitals)
3319 orbitals = matmul(orbitals, dense_antihermitian_exp(generator))
3320 energy = 0.0_dp
3321 DO i = 1, SIZE(occupation)
3322 energy = energy + occupation(i)* &
3323 REAL(DOT_PRODUCT(orbitals(:, i), MATMUL(hamiltonian, orbitals(:, i))), KIND=dp)
3324 END DO
3325
3326 END FUNCTION strict_rotated_energy
3327
3328END PROGRAM qs_ot_complex_ref_unittest
3329
cp_dbcsr_api::dbcsr_create
Definition cp_dbcsr_api.F:194
cp_fm_types::cp_fm_release
Definition cp_fm_types.F:88
preconditioner::apply_preconditioner
Definition preconditioner.F:115
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:122
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_cfm_types
Represents a complex full matrix distributed on many processors.
Definition cp_cfm_types.F:12
cp_cfm_types::cp_cfm_get_element
subroutine, public cp_cfm_get_element(matrix, irow_global, icol_global, alpha)
Get the matrix element by its global index.
Definition cp_cfm_types.F:255
cp_cfm_types::cp_cfm_release
subroutine, public cp_cfm_release(matrix)
Releases a full matrix.
Definition cp_cfm_types.F:182
cp_cfm_types::cp_cfm_set_element
subroutine, public cp_cfm_set_element(matrix, irow_global, icol_global, alpha)
Set the matrix element (irow_global,icol_global) of the full matrix to alpha.
Definition cp_cfm_types.F:302
cp_cfm_types::cp_cfm_create
subroutine, public cp_cfm_create(matrix, matrix_struct, name, nrow, ncol, set_zero)
Creates a new full matrix with the given structure.
Definition cp_cfm_types.F:125
cp_cfm_types::cp_cfm_set_all
subroutine, public cp_cfm_set_all(matrix, alpha, beta)
Set all elements of the full matrix to alpha. Besides, set all diagonal matrix elements to beta (if g...
Definition cp_cfm_types.F:202
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_fm_basic_linalg
Basic linear algebra operations for full matrices.
Definition cp_fm_basic_linalg.F:14
cp_fm_basic_linalg::cp_complex_fm_gemm
subroutine, public cp_complex_fm_gemm(transa, transb, m, n, k, alpha, a_re, a_im, b_re, b_im, beta, c_re, c_im, a_first_col, a_first_row, b_first_col, b_first_row, c_first_col, c_first_row)
Convenience function. Computes the matrix multiplications needed for the multiplication of complex ma...
Definition cp_fm_basic_linalg.F:2181
cp_fm_struct
represent the structure of a full matrix
Definition cp_fm_struct.F:14
cp_fm_struct::cp_fm_struct_create
subroutine, public cp_fm_struct_create(fmstruct, para_env, context, nrow_global, ncol_global, nrow_block, ncol_block, descriptor, first_p_pos, local_leading_dimension, template_fmstruct, square_blocks, force_block)
allocates and initializes a full matrix structure
Definition cp_fm_struct.F:132
cp_fm_struct::cp_fm_struct_release
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
Definition cp_fm_struct.F:355
cp_fm_types
represent a full matrix distributed on many processors
Definition cp_fm_types.F:15
cp_fm_types::cp_fm_set_submatrix
subroutine, public cp_fm_set_submatrix(fm, new_values, start_row, start_col, n_rows, n_cols, alpha, beta, transpose)
sets a submatrix of a full matrix fm(start_row:start_row+n_rows,start_col:start_col+n_cols) = alpha*o...
Definition cp_fm_types.F:750
cp_fm_types::cp_fm_create
subroutine, public cp_fm_create(matrix, matrix_struct, name, nrow, ncol, set_zero)
creates a new full matrix with the given structure
Definition cp_fm_types.F:165
cp_fm_types::cp_fm_get_submatrix
subroutine, public cp_fm_get_submatrix(fm, target_m, start_row, start_col, n_rows, n_cols, transpose)
gets a submatrix of a full matrix op(target_m)(1:n_rows,1:n_cols) =fm(start_row:start_row+n_rows,...
Definition cp_fm_types.F:925
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
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:18
input_constants::ot_precond_solver_chebyshev
integer, parameter, public ot_precond_solver_chebyshev
Definition input_constants.F:547
input_constants::ot_precond_full_all_covariant
integer, parameter, public ot_precond_full_all_covariant
Definition input_constants.F:533
input_constants::ot_precond_full_kinetic
integer, parameter, public ot_precond_full_kinetic
Definition input_constants.F:533
input_constants::ot_precond_solver_default
integer, parameter, public ot_precond_solver_default
Definition input_constants.F:547
input_constants::ot_precond_full_single
integer, parameter, public ot_precond_full_single
Definition input_constants.F:533
input_constants::ot_precond_solver_inv_chol
integer, parameter, public ot_precond_solver_inv_chol
Definition input_constants.F:547
input_constants::ot_precond_none
integer, parameter, public ot_precond_none
Definition input_constants.F:533
input_constants::ot_precond_full_single_inverse
integer, parameter, public ot_precond_full_single_inverse
Definition input_constants.F:533
input_constants::ot_precond_fermi_low_rank
integer, parameter, public ot_precond_fermi_low_rank
Definition input_constants.F:533
input_constants::ot_precond_s_inverse
integer, parameter, public ot_precond_s_inverse
Definition input_constants.F:533
input_constants::ot_precond_full_all
integer, parameter, public ot_precond_full_all
Definition input_constants.F:533
input_constants::ot_precond_solver_direct
integer, parameter, public ot_precond_solver_direct
Definition input_constants.F:547
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
kpoint_methods
Routines needed for kpoint calculation.
Definition kpoint_methods.F:15
kpoint_methods::kpoint_ot_energy_weighted_density
subroutine, public kpoint_ot_energy_weighted_density(coeff_re, coeff_im, hc_re, hc_im, occupation, wmat_re, wmat_im)
Build W(k) for noncanonical complex OT orbitals from H(k) C(k). The occupied-space Lagrange multiplie...
Definition kpoint_methods.F:2059
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:46
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:1882
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:1388
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:2162
message_passing::mp_comm_self
type(mp_comm_type), parameter, public mp_comm_self
Definition message_passing.F:816
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:1480
preconditioner_makes
computes preconditioners, and implements methods to apply them currently used in qs_ot
Definition preconditioner_makes.F:15
preconditioner_makes::make_complex_full_single_inverse
subroutine, public make_complex_full_single_inverse(preconditioner_env, matrix_c0, matrix_h, matrix_s, energy_gap)
Build a gauge-covariant FULL_SINGLE_INVERSE operator for a complex k-point channel.
Definition preconditioner_makes.F:1318
preconditioner_makes::make_complex_full_kinetic
subroutine, public make_complex_full_kinetic(preconditioner_env, matrix_t, matrix_s, energy_gap)
Build the inverse complex kinetic-plus-overlap preconditioner.
Definition preconditioner_makes.F:1441
preconditioner_makes::make_complex_full_s_inverse
subroutine, public make_complex_full_s_inverse(preconditioner_env, matrix_s)
Build the inverse complex overlap preconditioner.
Definition preconditioner_makes.F:1425
preconditioner_makes::make_complex_full_all_covariant
subroutine, public make_complex_full_all_covariant(preconditioner_env, matrix_c0, matrix_h, matrix_s, energy_gap, occupation_signature)
Build the gauge-covariant FULL_ALL inverse for a complex k-point channel.
Definition preconditioner_makes.F:385
preconditioner_makes::make_complex_full_single
subroutine, public make_complex_full_single(preconditioner_env, matrix_h, matrix_s, energy_homo, energy_gap)
Build the complex spectral FULL_SINGLE preconditioner.
Definition preconditioner_makes.F:236
preconditioner_makes::make_complex_fermi_low_rank
subroutine, public make_complex_fermi_low_rank(preconditioner_env, matrix_c0, matrix_h, matrix_s, energy_gap, max_rank, spectral_window, degeneracy_tolerance)
Build the bounded low-rank spectral inverse for a complex k-point channel.
Definition preconditioner_makes.F:320
preconditioner_makes::make_complex_full_all
subroutine, public make_complex_full_all(preconditioner_env, matrix_c0, matrix_h, matrix_s, c0_evals, energy_gap)
Build the state-selective FULL_ALL operator for a complex k-point channel. The occupied/reference sub...
Definition preconditioner_makes.F:104
preconditioner_types
types of preconditioners
Definition preconditioner_types.F:14
preconditioner_types::init_preconditioner
subroutine, public init_preconditioner(preconditioner_env, para_env, blacs_env)
...
Definition preconditioner_types.F:137
preconditioner_types::destroy_preconditioner
subroutine, public destroy_preconditioner(preconditioner_env)
...
Definition preconditioner_types.F:193
preconditioner
computes preconditioners, and implements methods to apply them currently used in qs_ot
Definition preconditioner.F:15
qs_ot_types
orbital transformations
Definition qs_ot_types.F:15
qs_ot_types::qs_ot_kpoint_preconditioner_supported
pure elemental logical function, public qs_ot_kpoint_preconditioner_supported(preconditioner_type, use_real_wfn)
Return whether a preconditioner is implemented for complex k-point OT.
Definition qs_ot_types.F:410
qs_ot_types::qs_ot_kpoint_preconditioner_scale
pure elemental real(kind=dp) function, public qs_ot_kpoint_preconditioner_scale(kpoint_weight)
Scale an inverse k-point Hessian block consistently with its irreducible weight.
Definition qs_ot_types.F:463
qs_ot_types::qs_ot_kpoint_preconditioner_solver_supported
pure elemental logical function, public qs_ot_kpoint_preconditioner_solver_supported(preconditioner_type, solver_type, use_real_wfn)
Return whether a solver is implemented for a complex k-point OT preconditioner.
Definition qs_ot_types.F:437
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:2294
qs_ot::qs_ot_density_secant_hessian
subroutine, public qs_ot_density_secant_hessian(density_step, hamiltonian_step, density_modes, correction, valid, density_norm_sq, response_work)
project a self-adjoint density/Hamiltonian secant onto density-response modes
Definition qs_ot.F:527
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:1986
qs_ot::qs_ot_symmetric_abs_solve
subroutine, public qs_ot_symmetric_abs_solve(matrix, rhs, solution, valid, relative_floor)
apply a positive spectral inverse of a real symmetric response matrix
Definition qs_ot.F:324
qs_ot::qs_ot_symmetric_sr1_update
pure subroutine, public qs_ot_symmetric_sr1_update(matrix, step, response, updated_matrix, valid, relative_tolerance)
add one accepted symmetric response secant to a reference Hessian
Definition qs_ot.F:467
qs_ot::qs_ot_fixed_n_projector_frechet
subroutine, public qs_ot_fixed_n_projector_frechet(chc, dchc, occupation, kpoint_weight, response_weight, fixed_n_weight_sum, projector_derivative, density_factor)
fixed-N Frechet derivative of a smooth occupation projector
Definition qs_ot.F:872
qs_ot::qs_ot_fixed_n_schur_block
subroutine, public qs_ot_fixed_n_schur_block(rotation_hessian, rayleigh_response, response_weight, rotation_gradient, energy_gradient, schur_block, coupling_vector, schur_rhs)
local block of the fixed-N rotation/energy Schur complement
Definition qs_ot.F:229
qs_ot::qs_ot_finite_rotation_response
subroutine, public qs_ot_finite_rotation_response(chc, rotation_generator, occupation, kpoint_weight, rotation_gradient, rotation_hessian, rayleigh_response, difference_step)
finite complex REF rotation Hessian and Rayleigh-energy response
Definition qs_ot.F:964
qs_ot::qs_ot_fixed_n_energy_gradient
pure subroutine, public qs_ot_fixed_n_energy_gradient(rayleigh_energy, energy_coordinate, response_weight, fixed_n_weight_sum, fixed_n_weighted_residual, gradient)
fixed-N Mermin gradient in auxiliary-energy coordinates
Definition qs_ot.F:156
qs_ot::qs_ot_antihermitian_spectral_norm
real(kind=dp) function, public qs_ot_antihermitian_spectral_norm(rotation_generator)
spectral norm of a dense anti-Hermitian rotation generator
Definition qs_ot.F:102
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:3491
qs_ot::qs_ot_projected_response_update
subroutine, public qs_ot_projected_response_update(reference_hessian, response_correction, coefficients, valid, projected_gradient, relative_floor)
update a baseline response direction in a small positive physical-response subspace
Definition qs_ot.F:389
qs_ot::qs_ot_fixed_n_energy_hessian
pure subroutine, public qs_ot_fixed_n_energy_hessian(response_weight, fixed_n_weight_sum, hessian)
dense fixed-N occupation Hessian in auxiliary-energy coordinates
Definition qs_ot.F:183
qs_ot::qs_ot_density_secant_orbital_overlaps
subroutine, public qs_ot_density_secant_orbital_overlaps(overlap_start_current, occupation_start, occupation_current, hamiltonian_step_start, hamiltonian_step_current, density_modes, kpoint_weight, density_norm_sq, response_work, density_overlap, response_overlap, valid)
project a physical density/Hamiltonian secant between moving orbital subspaces
Definition qs_ot.F:768
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:1208
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:1241
qs_ot::qs_ot_density_secant_projected_hessian
subroutine, public qs_ot_density_secant_projected_hessian(density_norm_sq, response_work, density_overlap, response_overlap, correction, valid, secant_mode, secant_position)
form a projected self-adjoint Hxc response from distributed density-space overlaps
Definition qs_ot.F:592
qs_ot::qs_ot_rot_mat_derivative
subroutine, public qs_ot_rot_mat_derivative(qs_ot_env)
computes the derivative fields with respect to rot_mat_x
Definition qs_ot.F:2988
qs_ot::qs_ot_fixed_n_multigroup_schur_block
subroutine, public qs_ot_fixed_n_multigroup_schur_block(rotation_hessian, rayleigh_response, response_weight, response_group, rotation_gradient, energy_gradient, schur_block, coupling_matrix, schur_rhs)
eliminate spin-resolved auxiliary energies while retaining every fixed-N constraint
Definition qs_ot.F:270
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:668
qs_ot::qs_ot_generate_rotation
subroutine, public qs_ot_generate_rotation(qs_ot_env)
computes the rotation matrix rot_mat_u that is associated to a given rot_mat_x using rot_mat_u=exp(ro...
Definition qs_ot.F:2923
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:2892
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
test_split_complex_fm_gemm
subroutine test_split_complex_fm_gemm(para_env, io_unit, nfail)
Check split-complex full-matrix multiplication including adjoints.
Definition qs_ot_complex_ref_unittest.F:280
qs_ot_complex_ref_unittest
program qs_ot_complex_ref_unittest
Regressions for complex k-point OT derivatives and preconditioner guards.
Definition qs_ot_complex_ref_unittest.F:11
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_cfm_types::cp_cfm_type
Represent a complex full matrix.
Definition cp_cfm_types.F:69
cp_dbcsr_api::dbcsr_distribution_type
Definition cp_dbcsr_api.F:182
cp_dbcsr_api::dbcsr_type
Definition cp_dbcsr_api.F:176
cp_fm_struct::cp_fm_struct_type
keeps the information about the structure of a full matrix
Definition cp_fm_struct.F:82
cp_fm_types::cp_fm_type
represent a full matrix
Definition cp_fm_types.F:114
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:178
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:767
preconditioner_types::preconditioner_type
Definition preconditioner_types.F:60
qs_ot_types::qs_ot_type
Definition qs_ot_types.F:139