Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Regressions for complex k-point OT derivatives and preconditioner guards.
10 : ! **************************************************************************************************
11 320 : PROGRAM qs_ot_complex_ref_unittest
12 2 : 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 4 : INTEGER, DIMENSION(:), POINTER :: col_dist, col_k, row_dist, row_k, row_n
91 2 : INTEGER, DIMENSION(:, :), POINTER :: pgrid
92 : LOGICAL :: found
93 : REAL(KIND=dp) :: error, fd_slope, hstep, old_error, slope
94 2 : 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 324 : TYPE(qs_ot_type) :: ot_env
104 :
105 2 : nfail = 0
106 2 : CALL test_preconditioner_support(nfail)
107 2 : NULLIFY (hc_im_p, hc_re_p, logger, para_env)
108 2 : CALL mp_world_init(mp_comm)
109 2 : mynode = mp_comm%mepos
110 2 : io_unit = -1
111 2 : IF (mynode == 0) io_unit = default_output_unit
112 2 : CALL test_antihermitian_spectral_norm(io_unit, nfail)
113 2 : CALL test_fixed_n_mermin_energy(mynode, nfail)
114 2 : CALL test_fixed_n_projector_frechet(mynode, nfail)
115 2 : CALL test_finite_rotation_response(mynode, nfail)
116 2 : CALL test_fixed_n_rotation_schur(mynode, nfail)
117 2 : CALL test_fixed_n_multigroup_schur(mynode, nfail)
118 2 : CALL test_roks_fixed_n_rotation_schur(mynode, nfail)
119 2 : CALL test_symmetric_abs_solve(mynode, nfail)
120 2 : CALL test_projected_response_update(mynode, nfail)
121 2 : CALL test_symmetric_sr1_update(mynode, nfail)
122 2 : CALL test_density_secant_hessian(mynode, nfail)
123 2 : CALL test_joint_spin_density_secant(mynode, nfail)
124 2 : CALL test_density_tangent(mynode, nfail)
125 2 : CALL test_density_secant_moving_subspace(mynode, nfail)
126 2 : ALLOCATE (para_env)
127 2 : CALL para_env%from_dup(mp_comm)
128 : CALL cp_logger_create(logger, para_env=para_env, default_global_unit_nr=io_unit, &
129 2 : close_global_unit_on_dealloc=.FALSE.)
130 2 : CALL cp_add_default_logger(logger)
131 2 : CALL add_all_references()
132 2 : CALL dbcsr_init_lib(mp_comm%get_handle(), io_unit)
133 2 : CALL test_real_rotation_frechet(para_env, nfail)
134 2 : CALL test_complex_rotation_frechet(para_env, nfail)
135 2 : CALL test_sparse_frechet_patterns(para_env, nfail)
136 2 : CALL test_complex_ref_rotation(para_env, nfail)
137 2 : CALL test_split_complex_fm_gemm(para_env, io_unit, nfail)
138 2 : CALL test_complex_preconditioner_gauge(para_env, nfail)
139 2 : CALL test_complex_occupation_block_preconditioner(para_env, nfail)
140 2 : 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 8 : 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 8 : 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 8 : 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 8 : CMPLX(-0.14_dp, -0.11_dp, KIND=dp)]
150 :
151 : hamiltonian(:, :) = CMPLX(0.0_dp, 0.0_dp, KIND=dp)
152 2 : hamiltonian(1, 1) = CMPLX(0.70_dp, 0.0_dp, KIND=dp)
153 2 : hamiltonian(2, 2) = CMPLX(-0.40_dp, 0.0_dp, KIND=dp)
154 2 : hamiltonian(3, 3) = CMPLX(1.10_dp, 0.0_dp, KIND=dp)
155 2 : hamiltonian(1, 2) = CMPLX(0.20_dp, 0.10_dp, KIND=dp)
156 2 : hamiltonian(2, 1) = CONJG(hamiltonian(1, 2))
157 2 : hamiltonian(1, 3) = CMPLX(-0.10_dp, 0.05_dp, KIND=dp)
158 2 : hamiltonian(3, 1) = CONJG(hamiltonian(1, 3))
159 2 : hamiltonian(2, 3) = CMPLX(0.15_dp, 0.20_dp, KIND=dp)
160 2 : hamiltonian(3, 2) = CONJG(hamiltonian(2, 3))
161 :
162 38 : CALL dense_inverse_sqrt(MATMUL(CONJG(TRANSPOSE(x)), x), inverse_sqrt)
163 38 : q = MATMUL(x, inverse_sqrt)
164 70 : hc = 2.0_dp*MATMUL(hamiltonian, q)
165 :
166 2 : ALLOCATE (pgrid(0:0, 0:0), row_dist(1), col_dist(1), row_n(1), row_k(1), col_k(1))
167 6 : pgrid(:, :) = 0
168 4 : row_dist(:) = 0
169 4 : col_dist(:) = 0
170 4 : row_n(:) = n
171 4 : row_k(:) = k
172 4 : col_k(:) = k
173 : CALL dbcsr_distribution_new(dist_nk, group=mp_comm_self%get_handle(), pgrid=pgrid, &
174 2 : row_dist=row_dist, col_dist=col_dist)
175 : CALL dbcsr_distribution_new(dist_kk, group=mp_comm_self%get_handle(), pgrid=pgrid, &
176 2 : row_dist=row_dist, col_dist=col_dist)
177 2 : CALL test_complex_strict_derivative(dist_nk, dist_kk, row_n, row_k, col_k, para_env, nfail)
178 :
179 2 : CALL create_one_block(x_re, "complex_ref_x_re", dist_nk, row_n, col_k)
180 2 : CALL create_one_block(x_im, "complex_ref_x_im", dist_nk, row_n, col_k)
181 2 : CALL create_one_block(sx_re, "complex_ref_sx_re", dist_nk, row_n, col_k)
182 2 : CALL create_one_block(sx_im, "complex_ref_sx_im", dist_nk, row_n, col_k)
183 2 : CALL create_one_block(hc_re, "complex_ref_hc_re", dist_nk, row_n, col_k)
184 2 : CALL create_one_block(hc_im, "complex_ref_hc_im", dist_nk, row_n, col_k)
185 2 : CALL create_one_block(g_re, "complex_ref_g_re", dist_nk, row_n, col_k)
186 2 : CALL create_one_block(g_im, "complex_ref_g_im", dist_nk, row_n, col_k)
187 2 : CALL create_one_block(q_re, "complex_ref_q_re", dist_nk, row_n, col_k)
188 2 : CALL create_one_block(q_im, "complex_ref_q_im", dist_nk, row_n, col_k)
189 2 : CALL create_one_block(f_re, "complex_ref_f_re", dist_kk, row_k, col_k)
190 2 : CALL create_one_block(f_im, "complex_ref_f_im", dist_kk, row_k, col_k)
191 2 : CALL create_one_block(b_re, "complex_ref_b_re", dist_kk, row_k, col_k)
192 2 : CALL create_one_block(b_im, "complex_ref_b_im", dist_kk, row_k, col_k)
193 2 : CALL create_one_block(tmp_kk, "complex_ref_tmp_kk", dist_kk, row_k, col_k)
194 :
195 2 : CALL put_complex_pair(x_re, x_im, x)
196 2 : CALL put_complex_pair(sx_re, sx_im, x)
197 2 : CALL put_complex_pair(hc_re, hc_im, hc)
198 2 : CALL put_complex_pair(f_re, f_im, inverse_sqrt)
199 :
200 2 : ot_env%has_complex_kpoint_state = .TRUE.
201 2 : ot_env%settings%eps_irac_filter_matrix = 0.0_dp
202 2 : ot_env%matrix_x => x_re
203 2 : ot_env%matrix_x_im => x_im
204 2 : ot_env%matrix_sx => sx_re
205 2 : ot_env%matrix_sx_im => sx_im
206 2 : ot_env%matrix_gx => g_re
207 2 : ot_env%matrix_gx_im => g_im
208 2 : ot_env%matrix_ref_inv_sqrt => f_re
209 2 : ot_env%matrix_ref_inv_sqrt_im => f_im
210 2 : ot_env%buf1_k_k_sym => b_re
211 2 : ot_env%buf2_k_k_sym => b_im
212 2 : ot_env%buf3_k_k_sym => tmp_kk
213 2 : ot_env%buf1_n_k => q_re
214 2 : ot_env%buf1_n_k_dp => q_im
215 :
216 2 : hc_re_p => hc_re
217 2 : hc_im_p => hc_im
218 2 : CALL qs_ot_get_derivative_ref_complex(hc_re_p, hc_im_p, ot_env)
219 :
220 2 : gradient(:, :) = CMPLX(0.0_dp, 0.0_dp, KIND=dp)
221 2 : CALL dbcsr_get_block_p(g_re, 1, 1, block_p, found)
222 18 : IF (found) gradient = gradient + CMPLX(block_p, 0.0_dp, KIND=dp)
223 2 : CALL dbcsr_get_block_p(g_im, 1, 1, block_p, found)
224 18 : IF (found) gradient = gradient + CMPLX(0.0_dp, block_p, KIND=dp)
225 :
226 2 : hstep = 1.0E-6_dp
227 : fd_slope = (polar_energy(x + hstep*direction, hamiltonian) - &
228 34 : polar_energy(x - hstep*direction, hamiltonian))/(2.0_dp*hstep)
229 18 : slope = REAL(SUM(CONJG(gradient)*direction), KIND=dp)
230 2 : error = ABS(slope - fd_slope)
231 :
232 38 : b = MATMUL(CONJG(TRANSPOSE(x)), hc)
233 54 : gradient_old = hc - MATMUL(x, b)
234 18 : old_error = ABS(REAL(SUM(CONJG(gradient_old)*direction), KIND=dp) - fd_slope)
235 2 : IF (error > 5.0E-8_dp .OR. old_error < 1.0E-3_dp) nfail = nfail + 1
236 :
237 2 : IF (mynode == 0) THEN
238 1 : WRITE (io_unit, '(A,3(1X,ES13.6))') "complex finite-REF pullback: fd/predicted/error", &
239 2 : fd_slope, slope, error
240 1 : WRITE (io_unit, '(A,1X,ES13.6)') "discarded projected-gradient error", old_error
241 : END IF
242 :
243 2 : CALL dbcsr_release(b_im)
244 2 : CALL dbcsr_release(b_re)
245 2 : CALL dbcsr_release(tmp_kk)
246 2 : CALL dbcsr_release(f_im)
247 2 : CALL dbcsr_release(f_re)
248 2 : CALL dbcsr_release(q_im)
249 2 : CALL dbcsr_release(q_re)
250 2 : CALL dbcsr_release(g_im)
251 2 : CALL dbcsr_release(g_re)
252 2 : CALL dbcsr_release(hc_im)
253 2 : CALL dbcsr_release(hc_re)
254 2 : CALL dbcsr_release(sx_im)
255 2 : CALL dbcsr_release(sx_re)
256 2 : CALL dbcsr_release(x_im)
257 2 : CALL dbcsr_release(x_re)
258 2 : CALL dbcsr_distribution_release(dist_kk)
259 2 : CALL dbcsr_distribution_release(dist_nk)
260 2 : DEALLOCATE (pgrid, row_dist, col_dist, row_n, row_k, col_k)
261 :
262 2 : CALL dbcsr_finalize_lib()
263 2 : CALL remove_all_references()
264 2 : CALL cp_rm_default_logger()
265 2 : CALL cp_logger_release(logger)
266 2 : CALL mp_para_env_release(para_env)
267 2 : CALL mp_world_finalize()
268 :
269 6 : IF (nfail > 0) ERROR STOP "qs_ot_complex_ref_unittest failed"
270 :
271 : CONTAINS
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 2 : 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 2 : NULLIFY (blacs_env, matrix_struct)
295 8 : DO j = 1, ndim
296 26 : DO i = 1, ndim
297 18 : a(i, j) = CMPLX(0.13_dp*i - 0.07_dp*j, 0.05_dp*i*j - 0.11_dp*j, KIND=dp)
298 24 : 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 2 : 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 2 : context=blacs_env, para_env=para_env)
305 2 : CALL cp_fm_create(a_re, matrix_struct)
306 2 : CALL cp_fm_create(a_im, matrix_struct)
307 2 : CALL cp_fm_create(b_re, matrix_struct)
308 2 : CALL cp_fm_create(b_im, matrix_struct)
309 2 : CALL cp_fm_create(c_re, matrix_struct)
310 2 : CALL cp_fm_create(c_im, matrix_struct)
311 26 : CALL cp_fm_set_submatrix(a_re, REAL(a, KIND=dp))
312 26 : CALL cp_fm_set_submatrix(a_im, AIMAG(a))
313 26 : CALL cp_fm_set_submatrix(b_re, REAL(b, KIND=dp))
314 26 : 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 2 : a_re, a_im, b_re, b_im, 0.0_dp, c_re, c_im)
318 2 : CALL cp_fm_get_submatrix(c_re, actual_re)
319 2 : CALL cp_fm_get_submatrix(c_im, actual_im)
320 26 : actual = CMPLX(actual_re, actual_im, KIND=dp)
321 80 : expected = MATMUL(a, b)
322 26 : error_nn = MAXVAL(ABS(actual - expected))
323 :
324 : CALL cp_complex_fm_gemm('C', 'N', ndim, ndim, ndim, 1.0_dp, &
325 2 : a_re, a_im, b_re, b_im, 0.0_dp, c_re, c_im)
326 2 : CALL cp_fm_get_submatrix(c_re, actual_re)
327 2 : CALL cp_fm_get_submatrix(c_im, actual_im)
328 26 : actual = CMPLX(actual_re, actual_im, KIND=dp)
329 80 : expected = MATMUL(CONJG(TRANSPOSE(a)), b)
330 26 : error_cn = MAXVAL(ABS(actual - expected))
331 :
332 : CALL cp_complex_fm_gemm('N', 'C', ndim, ndim, ndim, 1.0_dp, &
333 2 : a_re, a_im, b_re, b_im, 0.0_dp, c_re, c_im)
334 2 : CALL cp_fm_get_submatrix(c_re, actual_re)
335 2 : CALL cp_fm_get_submatrix(c_im, actual_im)
336 26 : actual = CMPLX(actual_re, actual_im, KIND=dp)
337 80 : expected = MATMUL(a, CONJG(TRANSPOSE(b)))
338 26 : error_nc = MAXVAL(ABS(actual - expected))
339 :
340 2 : IF (MAX(error_nn, error_cn, error_nc) > 1.0E-12_dp) nfail = nfail + 1
341 2 : IF (io_unit >= 0) WRITE (io_unit, '(A,3(1X,ES13.6))') &
342 1 : 'split-complex GEMM N/N, C/N, N/C errors:', error_nn, error_cn, error_nc
343 :
344 2 : CALL cp_fm_release(c_im)
345 2 : CALL cp_fm_release(c_re)
346 2 : CALL cp_fm_release(b_im)
347 2 : CALL cp_fm_release(b_re)
348 2 : CALL cp_fm_release(a_im)
349 2 : CALL cp_fm_release(a_re)
350 2 : CALL cp_fm_struct_release(matrix_struct)
351 2 : CALL cp_blacs_env_release(blacs_env)
352 :
353 6 : 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 2 : 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 2 : NULLIFY (blacs_env, coeff_struct, hmat_struct)
383 2 : pi = ACOS(-1.0_dp)
384 2 : hmat(:, :) = CMPLX(0.0_dp, 0.0_dp, KIND=dp)
385 2 : hmat(1, 1) = CMPLX(-0.8_dp, 0.0_dp, KIND=dp)
386 2 : hmat(2, 2) = CMPLX(0.1_dp, 0.0_dp, KIND=dp)
387 2 : hmat(3, 3) = CMPLX(0.7_dp, 0.0_dp, KIND=dp)
388 2 : hmat(4, 4) = CMPLX(1.4_dp, 0.0_dp, KIND=dp)
389 2 : hmat(1, 2) = CMPLX(0.21_dp, 0.13_dp, KIND=dp)
390 2 : hmat(1, 3) = CMPLX(-0.17_dp, 0.09_dp, KIND=dp)
391 2 : hmat(1, 4) = CMPLX(0.08_dp, -0.14_dp, KIND=dp)
392 2 : hmat(2, 3) = CMPLX(0.31_dp, -0.11_dp, KIND=dp)
393 2 : hmat(2, 4) = CMPLX(-0.19_dp, -0.07_dp, KIND=dp)
394 2 : hmat(3, 4) = CMPLX(0.27_dp, 0.16_dp, KIND=dp)
395 10 : DO i = 1, nao
396 22 : DO j = i + 1, nao
397 20 : hmat(j, i) = CONJG(hmat(i, j))
398 : END DO
399 : END DO
400 8 : DO j = 1, nmo
401 32 : DO i = 1, nao
402 24 : phase = 2.0_dp*pi*REAL((i - 1)*(j - 1), KIND=dp)/REAL(nao, KIND=dp)
403 30 : coeff(i, j) = CMPLX(COS(phase), SIN(phase), KIND=dp)/SQRT(REAL(nao, KIND=dp))
404 : END DO
405 : END DO
406 128 : hc_dense = MATMUL(hmat, coeff)
407 2 : occupation = [1.0_dp, 0.63_dp, 0.11_dp]
408 :
409 98 : hblock = MATMUL(CONJG(TRANSPOSE(coeff)), hc_dense)
410 8 : DO j = 1, nmo
411 26 : lagrange(:, j) = hblock(:, j)*occupation(j)
412 : END DO
413 50 : lagrange = 0.5_dp*(lagrange + CONJG(TRANSPOSE(lagrange)))
414 226 : expected = MATMUL(MATMUL(coeff, lagrange), CONJG(TRANSPOSE(coeff)))
415 :
416 2 : 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 2 : context=blacs_env, para_env=para_env)
419 : CALL cp_fm_struct_create(coeff_struct, nrow_global=nao, ncol_global=nmo, &
420 2 : context=blacs_env, para_env=para_env)
421 2 : CALL cp_fm_create(coeff_re_fm, coeff_struct)
422 2 : CALL cp_fm_create(coeff_im_fm, coeff_struct)
423 2 : CALL cp_fm_create(hc_re_fm, coeff_struct)
424 2 : CALL cp_fm_create(hc_im_fm, coeff_struct)
425 2 : CALL cp_fm_create(wmat_re_fm, hmat_struct)
426 2 : CALL cp_fm_create(wmat_im_fm, hmat_struct)
427 32 : CALL cp_fm_set_submatrix(coeff_re_fm, REAL(coeff, KIND=dp))
428 32 : CALL cp_fm_set_submatrix(coeff_im_fm, AIMAG(coeff))
429 32 : CALL cp_fm_set_submatrix(hc_re_fm, REAL(hc_dense, KIND=dp))
430 32 : 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 2 : wmat_re_fm, wmat_im_fm)
435 2 : CALL cp_fm_get_submatrix(wmat_re_fm, actual_re)
436 2 : CALL cp_fm_get_submatrix(wmat_im_fm, actual_im)
437 42 : actual = CMPLX(actual_re, actual_im, KIND=dp)
438 42 : error = MAXVAL(ABS(actual - expected))
439 42 : hermitian_error = MAXVAL(ABS(actual - CONJG(TRANSPOSE(actual))))
440 :
441 2 : rotation(:, :) = CMPLX(0.0_dp, 0.0_dp, KIND=dp)
442 2 : rotation(1, 1) = CMPLX(SQRT(0.5_dp), 0.0_dp, KIND=dp)
443 2 : rotation(1, 2) = CMPLX(0.0_dp, SQRT(0.5_dp), KIND=dp)
444 2 : rotation(2, 1) = CMPLX(0.0_dp, SQRT(0.5_dp), KIND=dp)
445 2 : rotation(2, 2) = CMPLX(SQRT(0.5_dp), 0.0_dp, KIND=dp)
446 2 : rotation(3, 3) = CMPLX(0.0_dp, 1.0_dp, KIND=dp)
447 98 : coeff_rot = MATMUL(coeff, rotation)
448 98 : hc_rot = MATMUL(hc_dense, rotation)
449 8 : uniform_occupation(:) = 1.0_dp
450 32 : CALL cp_fm_set_submatrix(coeff_re_fm, REAL(coeff_rot, KIND=dp))
451 32 : CALL cp_fm_set_submatrix(coeff_im_fm, AIMAG(coeff_rot))
452 32 : CALL cp_fm_set_submatrix(hc_re_fm, REAL(hc_rot, KIND=dp))
453 32 : 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 2 : wmat_re_fm, wmat_im_fm)
457 2 : CALL cp_fm_get_submatrix(wmat_re_fm, rotated_re)
458 2 : CALL cp_fm_get_submatrix(wmat_im_fm, rotated_im)
459 42 : rotated = CMPLX(rotated_re, rotated_im, KIND=dp)
460 226 : expected = MATMUL(MATMUL(coeff, hblock), CONJG(TRANSPOSE(coeff)))
461 42 : gauge_error = MAXVAL(ABS(rotated - expected))
462 :
463 2 : IF (error > 5.0E-12_dp .OR. hermitian_error > 5.0E-12_dp .OR. &
464 0 : gauge_error > 5.0E-12_dp) nfail = nfail + 1
465 2 : IF (io_unit >= 0) THEN
466 1 : WRITE (io_unit, '(A,3(1X,ES13.6))') 'complex OT W error/hermitian/gauge', &
467 2 : error, hermitian_error, gauge_error
468 : END IF
469 :
470 2 : CALL cp_fm_release(wmat_im_fm)
471 2 : CALL cp_fm_release(wmat_re_fm)
472 2 : CALL cp_fm_release(hc_im_fm)
473 2 : CALL cp_fm_release(hc_re_fm)
474 2 : CALL cp_fm_release(coeff_im_fm)
475 2 : CALL cp_fm_release(coeff_re_fm)
476 2 : CALL cp_fm_struct_release(coeff_struct)
477 2 : CALL cp_fm_struct_release(hmat_struct)
478 2 : CALL cp_blacs_env_release(blacs_env)
479 :
480 4 : 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 2 : 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 2 : generator(:, :) = CMPLX(0.0_dp, 0.0_dp, KIND=dp)
496 2 : generator(1, 2) = CMPLX(0.3_dp, 0.4_dp, KIND=dp)
497 2 : generator(2, 1) = -CONJG(generator(1, 2))
498 2 : phase = EXP(CMPLX(0.0_dp, 1.0_dp, KIND=dp)*[0.37_dp, -0.51_dp, 0.83_dp])
499 8 : DO j = 1, 3
500 26 : DO i = 1, 3
501 24 : transformed(i, j) = CONJG(phase(i))*generator(i, j)*phase(j)
502 : END DO
503 : END DO
504 :
505 2 : norm = qs_ot_antihermitian_spectral_norm(generator)
506 2 : transformed_norm = qs_ot_antihermitian_spectral_norm(transformed)
507 2 : error = MAX(ABS(norm - 0.5_dp), ABS(transformed_norm - norm))
508 2 : IF (io_unit >= 0) WRITE (io_unit, '(A,ES14.6)') &
509 1 : "complex rotation spectral-norm error: ", error
510 2 : IF (error > 1.0E-13_dp) nfail = nfail + 1
511 :
512 2 : 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 2 : 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 2 : energies(:) = [-0.24_dp, -0.05_dp, 0.08_dp, 0.31_dp]
537 2 : rayleigh(:) = [-0.19_dp, -0.02_dp, 0.04_dp, 0.27_dp]
538 2 : weights(:) = [0.35_dp, 0.35_dp, 0.65_dp, 0.65_dp]
539 2 : 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 2 : maxocc, occupation, mu)
543 10 : response(:) = weights(:)*occupation(:)*(maxocc - occupation(:))/(maxocc*temperature)
544 10 : response_sum = SUM(response)
545 10 : weighted_residual = DOT_PRODUCT(response, rayleigh - energies)
546 : CALL qs_ot_fixed_n_energy_gradient(rayleigh, energies, response, response_sum, &
547 2 : 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 18 : target_electrons, temperature, maxocc))/(2.0_dp*fd_step)
553 10 : predicted_slope = DOT_PRODUCT(gradient, direction)
554 10 : 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 2 : rayleigh(:) = energies(:)
558 2 : CALL qs_ot_fixed_n_energy_hessian(response, response_sum, hessian)
559 42 : hessian_action(:) = MATMUL(hessian, direction)
560 :
561 : CALL fixed_n_fermi_occupations(energies + fd_step*direction, weights, target_electrons, &
562 10 : temperature, maxocc, occupation, mu)
563 10 : response_plus(:) = weights(:)*occupation(:)*(maxocc - occupation(:))/(maxocc*temperature)
564 10 : 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 18 : SUM(response_plus), weighted_residual, gradient_plus)
567 :
568 : CALL fixed_n_fermi_occupations(energies - fd_step*direction, weights, target_electrons, &
569 10 : temperature, maxocc, occupation, mu)
570 10 : response_minus(:) = weights(:)*occupation(:)*(maxocc - occupation(:))/(maxocc*temperature)
571 10 : 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 18 : SUM(response_minus), weighted_residual, gradient_minus)
574 :
575 : hessian_error = MAXVAL(ABS((gradient_plus - gradient_minus)/(2.0_dp*fd_step) - &
576 10 : hessian_action))
577 42 : hessian_error = MAX(hessian_error, MAXVAL(ABS(hessian - TRANSPOSE(hessian))))
578 52 : hessian_error = MAX(hessian_error, MAXVAL(ABS(MATMUL(hessian, SPREAD(1.0_dp, 1, nbands)))))
579 :
580 2 : IF (gradient_error > 2.0E-8_dp .OR. hessian_error > 2.0E-7_dp) nfail = nfail + 1
581 2 : IF (mynode == 0) THEN
582 1 : WRITE (io_unit, '(A,2(1X,ES13.6))') "fixed-N Mermin gradient/Hessian errors:", &
583 2 : gradient_error, hessian_error
584 : END IF
585 :
586 2 : 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 2 : 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 2 : hamiltonian(:, :) = CMPLX(0.0_dp, 0.0_dp, KIND=dp)
610 2 : hamiltonian(1, 1) = CMPLX(-0.22_dp, 0.0_dp, KIND=dp)
611 2 : hamiltonian(2, 2) = CMPLX(0.03_dp, 0.0_dp, KIND=dp)
612 2 : hamiltonian(3, 3) = CMPLX(0.28_dp, 0.0_dp, KIND=dp)
613 2 : hamiltonian(1, 2) = CMPLX(0.04_dp, -0.03_dp, KIND=dp)
614 2 : hamiltonian(2, 1) = CONJG(hamiltonian(1, 2))
615 2 : hamiltonian(2, 3) = CMPLX(-0.02_dp, 0.05_dp, KIND=dp)
616 2 : hamiltonian(3, 2) = CONJG(hamiltonian(2, 3))
617 :
618 2 : dchc(:, :) = CMPLX(0.0_dp, 0.0_dp, KIND=dp)
619 2 : dchc(1, 1) = CMPLX(0.13_dp, 0.0_dp, KIND=dp)
620 2 : dchc(2, 2) = CMPLX(-0.08_dp, 0.0_dp, KIND=dp)
621 2 : dchc(3, 3) = CMPLX(0.05_dp, 0.0_dp, KIND=dp)
622 2 : dchc(1, 2) = CMPLX(0.17_dp, -0.11_dp, KIND=dp)
623 2 : dchc(2, 1) = CONJG(dchc(1, 2))
624 2 : dchc(2, 3) = CMPLX(0.12_dp, 0.09_dp, KIND=dp)
625 2 : dchc(3, 2) = CONJG(dchc(2, 3))
626 :
627 2 : CALL diag_complex(hamiltonian, vectors, eigenvalues)
628 8 : weights(:) = kpoint_weight
629 : CALL fixed_n_fermi_occupations(eigenvalues, weights, target_electrons, temperature, &
630 2 : maxocc, occupation, mu)
631 8 : response(:) = weights(:)*occupation(:)*(maxocc - occupation(:))/(maxocc*temperature)
632 : CALL qs_ot_fixed_n_projector_frechet(hamiltonian, dchc, occupation, kpoint_weight, &
633 8 : response, SUM(response), derivative)
634 :
635 : CALL fixed_n_spectral_projector(hamiltonian + fd_step*dchc, weights, target_electrons, &
636 26 : temperature, maxocc, projector_plus)
637 : CALL fixed_n_spectral_projector(hamiltonian - fd_step*dchc, weights, target_electrons, &
638 26 : temperature, maxocc, projector_minus)
639 26 : derivative_fd(:, :) = (projector_plus - projector_minus)/(2.0_dp*fd_step)
640 26 : error = MAXVAL(ABS(derivative - derivative_fd))
641 2 : trace_error = 0.0_dp
642 8 : DO i = 1, nbands
643 8 : trace_error = trace_error + REAL(derivative(i, i), KIND=dp)
644 : END DO
645 2 : error = MAX(error, ABS(trace_error))
646 :
647 2 : IF (error > 3.0E-7_dp) nfail = nfail + 1
648 2 : IF (mynode == 0) WRITE (io_unit, '(A,1X,ES13.6)') &
649 1 : "fixed-N complex occupation-projector error:", error
650 :
651 2 : 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 2 : 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 2 : base_hamiltonian(:, :) = CMPLX(0.0_dp, 0.0_dp, KIND=dp)
679 2 : base_hamiltonian(1, 1) = CMPLX(-0.31_dp, 0.0_dp, KIND=dp)
680 2 : base_hamiltonian(2, 2) = CMPLX(0.07_dp, 0.0_dp, KIND=dp)
681 2 : base_hamiltonian(3, 3) = CMPLX(0.42_dp, 0.0_dp, KIND=dp)
682 2 : base_hamiltonian(1, 2) = CMPLX(0.09_dp, -0.06_dp, KIND=dp)
683 2 : base_hamiltonian(1, 3) = CMPLX(-0.04_dp, 0.08_dp, KIND=dp)
684 2 : base_hamiltonian(2, 3) = CMPLX(0.11_dp, 0.05_dp, KIND=dp)
685 8 : DO i = 1, nbands
686 14 : base_hamiltonian(i, 1:i - 1) = CONJG(base_hamiltonian(1:i - 1, i))
687 : END DO
688 :
689 2 : generator(:, :) = CMPLX(0.0_dp, 0.0_dp, KIND=dp)
690 2 : generator(1, 2) = CMPLX(0.17_dp, -0.08_dp, KIND=dp)
691 2 : generator(1, 3) = CMPLX(-0.09_dp, 0.04_dp, KIND=dp)
692 2 : generator(2, 3) = CMPLX(0.12_dp, 0.07_dp, KIND=dp)
693 8 : DO i = 1, nbands
694 14 : generator(i, 1:i - 1) = -CONJG(generator(1:i - 1, i))
695 : END DO
696 :
697 2 : direction(:) = [0.14_dp, -0.07_dp, -0.11_dp, 0.09_dp, 0.05_dp, 0.13_dp]
698 2 : direction_generator(:, :) = CMPLX(0.0_dp, 0.0_dp, KIND=dp)
699 2 : direction_generator(1, 2) = CMPLX(direction(1), direction(2), KIND=dp)
700 2 : direction_generator(2, 1) = -CONJG(direction_generator(1, 2))
701 2 : direction_generator(1, 3) = CMPLX(direction(3), direction(4), KIND=dp)
702 2 : direction_generator(3, 1) = -CONJG(direction_generator(1, 3))
703 2 : direction_generator(2, 3) = CMPLX(direction(5), direction(6), KIND=dp)
704 2 : direction_generator(3, 2) = -CONJG(direction_generator(2, 3))
705 2 : occupation(:) = [1.73_dp, 0.88_dp, 0.19_dp]
706 :
707 2 : rotation(:, :) = dense_antihermitian_exp(generator)
708 : chc(:, :) = MATMUL(CONJG(TRANSPOSE(rotation)), &
709 158 : MATMUL(base_hamiltonian, rotation))
710 : CALL qs_ot_finite_rotation_response(chc, generator, occupation, kpoint_weight, &
711 2 : gradient, hessian, rayleigh_response)
712 :
713 26 : rotation_plus(:, :) = dense_antihermitian_exp(generator + fd_step*direction_generator)
714 26 : rotation_minus(:, :) = dense_antihermitian_exp(generator - fd_step*direction_generator)
715 2 : energy = finite_rotation_energy(rotation, base_hamiltonian, occupation, kpoint_weight)
716 2 : energy_plus = finite_rotation_energy(rotation_plus, base_hamiltonian, occupation, kpoint_weight)
717 2 : 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 14 : DOT_PRODUCT(gradient, direction))
720 : hessian_error = ABS((energy_plus - 2.0_dp*energy + energy_minus)/fd_step**2 - &
721 98 : DOT_PRODUCT(direction, MATMUL(hessian, direction)))
722 86 : symmetry_error = MAXVAL(ABS(hessian - TRANSPOSE(hessian)))
723 :
724 8 : DO i = 1, nbands
725 : rayleigh_plus(i) = REAL(DOT_PRODUCT(rotation_plus(:, i), &
726 96 : MATMUL(base_hamiltonian, rotation_plus(:, i))), KIND=dp)
727 : rayleigh_minus(i) = REAL(DOT_PRODUCT(rotation_minus(:, i), &
728 98 : 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 56 : MATMUL(rayleigh_response, direction)))
732 :
733 : IF (gradient_error > 5.0E-8_dp .OR. hessian_error > 2.0E-6_dp .OR. &
734 2 : rayleigh_error > 2.0E-7_dp .OR. symmetry_error > 1.0E-12_dp) nfail = nfail + 1
735 2 : IF (mynode == 0) THEN
736 1 : WRITE (io_unit, '(A,4(1X,ES13.6))') "finite complex rotation response errors:", &
737 2 : gradient_error, hessian_error, rayleigh_error, symmetry_error
738 : END IF
739 :
740 2 : 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 2 : 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 2 : base_hamiltonian(:, :) = CMPLX(0.0_dp, 0.0_dp, KIND=dp)
774 2 : base_hamiltonian(1, 1) = CMPLX(-0.31_dp, 0.0_dp, KIND=dp)
775 2 : base_hamiltonian(2, 2) = CMPLX(0.07_dp, 0.0_dp, KIND=dp)
776 2 : base_hamiltonian(3, 3) = CMPLX(0.42_dp, 0.0_dp, KIND=dp)
777 2 : base_hamiltonian(1, 2) = CMPLX(0.09_dp, -0.06_dp, KIND=dp)
778 2 : base_hamiltonian(1, 3) = CMPLX(-0.04_dp, 0.08_dp, KIND=dp)
779 2 : base_hamiltonian(2, 3) = CMPLX(0.11_dp, 0.05_dp, KIND=dp)
780 8 : DO i = 1, nbands
781 14 : base_hamiltonian(i, 1:i - 1) = CONJG(base_hamiltonian(1:i - 1, i))
782 : END DO
783 :
784 2 : generator(:, :) = CMPLX(0.0_dp, 0.0_dp, KIND=dp)
785 2 : generator(1, 2) = CMPLX(0.17_dp, -0.08_dp, KIND=dp)
786 2 : generator(1, 3) = CMPLX(-0.09_dp, 0.04_dp, KIND=dp)
787 2 : generator(2, 3) = CMPLX(0.12_dp, 0.07_dp, KIND=dp)
788 8 : DO i = 1, nbands
789 14 : generator(i, 1:i - 1) = -CONJG(generator(1:i - 1, i))
790 : END DO
791 2 : rotation_direction(:) = [0.14_dp, -0.07_dp, -0.11_dp, 0.09_dp, 0.05_dp, 0.13_dp]
792 2 : energy_direction(:) = [0.08_dp, -0.11_dp, 0.03_dp]
793 2 : direction_generator(:, :) = CMPLX(0.0_dp, 0.0_dp, KIND=dp)
794 2 : direction_generator(1, 2) = CMPLX(rotation_direction(1), rotation_direction(2), KIND=dp)
795 2 : direction_generator(2, 1) = -CONJG(direction_generator(1, 2))
796 2 : direction_generator(1, 3) = CMPLX(rotation_direction(3), rotation_direction(4), KIND=dp)
797 2 : direction_generator(3, 1) = -CONJG(direction_generator(1, 3))
798 2 : direction_generator(2, 3) = CMPLX(rotation_direction(5), rotation_direction(6), KIND=dp)
799 2 : direction_generator(3, 2) = -CONJG(direction_generator(2, 3))
800 :
801 2 : rotation(:, :) = dense_antihermitian_exp(generator)
802 158 : chc(:, :) = MATMUL(CONJG(TRANSPOSE(rotation)), MATMUL(base_hamiltonian, rotation))
803 8 : DO i = 1, nbands
804 8 : rayleigh(i) = REAL(chc(i, i), KIND=dp)
805 : END DO
806 2 : energy_coordinate(:) = rayleigh(:)
807 8 : weights(:) = kpoint_weight
808 : CALL fixed_n_fermi_occupations(energy_coordinate, weights, 1.7_dp, temperature, &
809 2 : maxocc, occupation, mu)
810 8 : response(:) = weights*occupation*(maxocc - occupation)/(maxocc*temperature)
811 8 : response_sum = SUM(response)
812 : CALL qs_ot_finite_rotation_response(chc, generator, occupation, kpoint_weight, &
813 2 : rotation_gradient, rotation_hessian, rayleigh_response)
814 : CALL qs_ot_fixed_n_energy_gradient(rayleigh, energy_coordinate, response, response_sum, &
815 2 : 0.0_dp, energy_gradient)
816 2 : 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 2 : schur_block, coupling, schur_rhs)
820 : schur(:, :) = schur_block + SPREAD(coupling, DIM=2, NCOPIES=nrotation)* &
821 86 : SPREAD(coupling, DIM=1, NCOPIES=nrotation)/response_sum
822 : schur_error = MAXVAL(ABS(schur - rotation_hessian + &
823 : MATMUL(TRANSPOSE(rayleigh_response), &
824 544 : MATMUL(energy_hessian, rayleigh_response))))
825 :
826 26 : rotation_plus(:, :) = dense_antihermitian_exp(generator + fd_step*direction_generator)
827 26 : rotation_minus(:, :) = dense_antihermitian_exp(generator - fd_step*direction_generator)
828 8 : DO i = 1, nbands
829 : rayleigh_plus(i) = REAL(DOT_PRODUCT(rotation_plus(:, i), &
830 96 : MATMUL(base_hamiltonian, rotation_plus(:, i))), KIND=dp)
831 : rayleigh_minus(i) = REAL(DOT_PRODUCT(rotation_minus(:, i), &
832 104 : MATMUL(base_hamiltonian, rotation_minus(:, i))), KIND=dp)
833 : END DO
834 8 : energy_plus_coordinate(:) = energy_coordinate + fd_step*energy_direction
835 8 : energy_minus_coordinate(:) = energy_coordinate - fd_step*energy_direction
836 : energy_zero = fixed_n_mermin_value(rayleigh, energy_coordinate, weights, 1.7_dp, &
837 2 : temperature, maxocc)
838 : energy_plus = fixed_n_mermin_value(rayleigh_plus, energy_plus_coordinate, weights, 1.7_dp, &
839 2 : temperature, maxocc)
840 : energy_minus = fixed_n_mermin_value(rayleigh_minus, energy_minus_coordinate, weights, 1.7_dp, &
841 2 : 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 212 : DOT_PRODUCT(energy_direction, MATMUL(energy_hessian, energy_direction))
847 2 : curvature_error = ABS((energy_plus - 2.0_dp*energy_zero + energy_minus)/fd_step**2 - curvature)
848 :
849 2 : IF (schur_error > 2.0E-11_dp .OR. curvature_error > 3.0E-6_dp) nfail = nfail + 1
850 2 : IF (mynode == 0) THEN
851 1 : WRITE (io_unit, '(A,2(1X,ES13.6))') "fixed-N rotation Schur errors:", &
852 2 : schur_error, curvature_error
853 : END IF
854 :
855 2 : 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 2 : 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 2 : [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 2 : [nenergy, nrotation])
884 2 : response_weight(:) = [0.7_dp, 0.5_dp, 0.4_dp, 0.2_dp]
885 2 : response_group(:) = [1, 1, 2, 2]
886 2 : rotation_gradient(:) = [0.2_dp, -0.1_dp, 0.3_dp]
887 2 : 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 2 : 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 250 : [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 88 : 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 10 : response_weight(1:2))))
908 : coupling_error = MAX(coupling_error, MAXVAL(ABS(coupling(:, 2) - &
909 : MATMUL(TRANSPOSE(rayleigh_response(3:4, :)), &
910 10 : response_weight(3:4)))))
911 : rhs_error = MAXVAL(ABS(schur_rhs - rotation_gradient - &
912 10 : MATMUL(TRANSPOSE(rayleigh_response), energy_gradient)))
913 :
914 2 : IF (MAX(schur_error, coupling_error, rhs_error) > 2.0E-14_dp) nfail = nfail + 1
915 2 : IF (mynode == 0) THEN
916 1 : WRITE (io_unit, '(A,3(1X,ES13.6))') "fixed-N multigroup Schur errors:", &
917 2 : schur_error, coupling_error, rhs_error
918 : END IF
919 :
920 2 : 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 2 : 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 2 : hamiltonian_a(:, :) = CMPLX(0.0_dp, 0.0_dp, KIND=dp)
961 2 : hamiltonian_b(:, :) = CMPLX(0.0_dp, 0.0_dp, KIND=dp)
962 2 : hamiltonian_a(1, 1) = -0.42_dp
963 2 : hamiltonian_a(2, 2) = 0.03_dp
964 2 : hamiltonian_a(3, 3) = 0.37_dp
965 2 : hamiltonian_a(1, 2) = CMPLX(0.08_dp, -0.05_dp, KIND=dp)
966 2 : hamiltonian_a(1, 3) = CMPLX(-0.03_dp, 0.07_dp, KIND=dp)
967 2 : hamiltonian_a(2, 3) = CMPLX(0.09_dp, 0.04_dp, KIND=dp)
968 2 : hamiltonian_b(1, 1) = -0.31_dp
969 2 : hamiltonian_b(2, 2) = 0.11_dp
970 2 : hamiltonian_b(3, 3) = 0.46_dp
971 2 : hamiltonian_b(1, 2) = CMPLX(-0.06_dp, 0.02_dp, KIND=dp)
972 2 : hamiltonian_b(1, 3) = CMPLX(0.05_dp, 0.04_dp, KIND=dp)
973 2 : hamiltonian_b(2, 3) = CMPLX(-0.07_dp, 0.06_dp, KIND=dp)
974 8 : DO i = 1, nbands
975 12 : hamiltonian_a(i, 1:i - 1) = CONJG(hamiltonian_a(1:i - 1, i))
976 14 : hamiltonian_b(i, 1:i - 1) = CONJG(hamiltonian_b(1:i - 1, i))
977 : END DO
978 :
979 2 : generator(:, :) = CMPLX(0.0_dp, 0.0_dp, KIND=dp)
980 2 : generator(1, 2) = CMPLX(0.14_dp, -0.06_dp, KIND=dp)
981 2 : generator(1, 3) = CMPLX(-0.08_dp, 0.05_dp, KIND=dp)
982 2 : generator(2, 3) = CMPLX(0.10_dp, 0.03_dp, KIND=dp)
983 8 : DO i = 1, nbands
984 14 : generator(i, 1:i - 1) = -CONJG(generator(1:i - 1, i))
985 : END DO
986 2 : rotation_direction(:) = [0.09_dp, -0.04_dp, -0.07_dp, 0.06_dp, 0.03_dp, 0.08_dp]
987 2 : energy_direction(:) = [0.05_dp, -0.04_dp, 0.02_dp, -0.03_dp, 0.06_dp, -0.01_dp]
988 2 : direction_generator(:, :) = CMPLX(0.0_dp, 0.0_dp, KIND=dp)
989 2 : direction_generator(1, 2) = CMPLX(rotation_direction(1), rotation_direction(2), KIND=dp)
990 2 : direction_generator(2, 1) = -CONJG(direction_generator(1, 2))
991 2 : direction_generator(1, 3) = CMPLX(rotation_direction(3), rotation_direction(4), KIND=dp)
992 2 : direction_generator(3, 1) = -CONJG(direction_generator(1, 3))
993 2 : direction_generator(2, 3) = CMPLX(rotation_direction(5), rotation_direction(6), KIND=dp)
994 2 : direction_generator(3, 2) = -CONJG(direction_generator(2, 3))
995 :
996 2 : rotation(:, :) = dense_antihermitian_exp(generator)
997 158 : chc_a(:, :) = MATMUL(CONJG(TRANSPOSE(rotation)), MATMUL(hamiltonian_a, rotation))
998 158 : chc_b(:, :) = MATMUL(CONJG(TRANSPOSE(rotation)), MATMUL(hamiltonian_b, rotation))
999 8 : DO i = 1, nbands
1000 6 : rayleigh_a(i) = REAL(chc_a(i, i), KIND=dp)
1001 8 : rayleigh_b(i) = REAL(chc_b(i, i), KIND=dp)
1002 : END DO
1003 2 : energy_a(:) = rayleigh_a
1004 2 : energy_b(:) = rayleigh_b
1005 8 : weights(:) = kpoint_weight
1006 : CALL fixed_n_fermi_occupations(energy_a, weights, 1.65_dp, temperature, &
1007 2 : maxocc, occupation_a, mu_a)
1008 : CALL fixed_n_fermi_occupations(energy_b, weights, 0.95_dp, temperature, &
1009 2 : maxocc, occupation_b, mu_b)
1010 8 : response_a(:) = weights*occupation_a*(maxocc - occupation_a)/(maxocc*temperature)
1011 8 : response_b(:) = weights*occupation_b*(maxocc - occupation_b)/(maxocc*temperature)
1012 8 : response_sum_a = SUM(response_a)
1013 8 : response_sum_b = SUM(response_b)
1014 8 : response(1:nbands) = response_a
1015 8 : response(nbands + 1:nenergy) = response_b
1016 2 : 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 2 : spin_rayleigh_response)
1021 2 : rotation_gradient(:) = spin_gradient
1022 2 : rotation_hessian_a(:, :) = spin_hessian
1023 50 : 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 2 : spin_rayleigh_response)
1027 14 : rotation_gradient(:) = rotation_gradient + spin_gradient
1028 2 : rotation_hessian_b(:, :) = spin_hessian
1029 86 : rotation_hessian(:, :) = rotation_hessian_a + rotation_hessian_b
1030 50 : 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 2 : 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 2 : energy_gradient(nbands + 1:nenergy))
1038 2 : energy_hessian(:, :) = 0.0_dp
1039 : CALL qs_ot_fixed_n_energy_hessian( &
1040 2 : response_a, response_sum_a, energy_hessian(1:nbands, 1:nbands))
1041 : CALL qs_ot_fixed_n_energy_hessian( &
1042 2 : 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 2 : energy_gradient, schur_block, coupling, schur_rhs)
1046 : schur(:, :) = schur_block + &
1047 : MATMUL(coupling(:, 1:1), TRANSPOSE(coupling(:, 1:1)))/response_sum_a + &
1048 280 : 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 1120 : MATMUL(energy_hessian, rayleigh_response))))
1052 :
1053 26 : rotation_plus(:, :) = dense_antihermitian_exp(generator + fd_step*direction_generator)
1054 26 : rotation_minus(:, :) = dense_antihermitian_exp(generator - fd_step*direction_generator)
1055 8 : DO i = 1, nbands
1056 : rayleigh_a_plus(i) = REAL(DOT_PRODUCT(rotation_plus(:, i), &
1057 96 : MATMUL(hamiltonian_a, rotation_plus(:, i))), KIND=dp)
1058 : rayleigh_a_minus(i) = REAL(DOT_PRODUCT(rotation_minus(:, i), &
1059 102 : MATMUL(hamiltonian_a, rotation_minus(:, i))), KIND=dp)
1060 : rayleigh_b_plus(i) = REAL(DOT_PRODUCT(rotation_plus(:, i), &
1061 102 : MATMUL(hamiltonian_b, rotation_plus(:, i))), KIND=dp)
1062 : rayleigh_b_minus(i) = REAL(DOT_PRODUCT(rotation_minus(:, i), &
1063 104 : MATMUL(hamiltonian_b, rotation_minus(:, i))), KIND=dp)
1064 : END DO
1065 8 : energy_a_plus(:) = energy_a + fd_step*energy_direction(1:nbands)
1066 8 : energy_a_minus(:) = energy_a - fd_step*energy_direction(1:nbands)
1067 8 : energy_b_plus(:) = energy_b + fd_step*energy_direction(nbands + 1:nenergy)
1068 8 : energy_b_minus(:) = energy_b - fd_step*energy_direction(nbands + 1:nenergy)
1069 : energy_zero_a = fixed_n_mermin_value( &
1070 2 : rayleigh_a, energy_a, weights, 1.65_dp, temperature, maxocc)
1071 : energy_zero_b = fixed_n_mermin_value( &
1072 2 : rayleigh_b, energy_b, weights, 0.95_dp, temperature, maxocc)
1073 : energy_plus_a = fixed_n_mermin_value( &
1074 2 : rayleigh_a_plus, energy_a_plus, weights, 1.65_dp, temperature, maxocc)
1075 : energy_plus_b = fixed_n_mermin_value( &
1076 2 : rayleigh_b_plus, energy_b_plus, weights, 0.95_dp, temperature, maxocc)
1077 : energy_minus_a = fixed_n_mermin_value( &
1078 2 : rayleigh_a_minus, energy_a_minus, weights, 1.65_dp, temperature, maxocc)
1079 : energy_minus_b = fixed_n_mermin_value( &
1080 2 : rayleigh_b_minus, energy_b_minus, weights, 0.95_dp, temperature, maxocc)
1081 2 : energy_zero = energy_zero_a + energy_zero_b
1082 2 : energy_plus = energy_plus_a + energy_plus_b
1083 2 : 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 212 : 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 212 : energy_direction(nbands + 1:nenergy)))
1103 2 : curvature = curvature_a + curvature_b
1104 : curvature_error = ABS((energy_plus_a - 2.0_dp*energy_zero_a + energy_minus_a)/ &
1105 2 : 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 2 : ABS((energy_plus - 2.0_dp*energy_zero + energy_minus)/fd_step**2 - curvature))
1110 :
1111 2 : IF (schur_error > 3.0E-11_dp .OR. curvature_error > 5.0E-6_dp) nfail = nfail + 1
1112 2 : IF (mynode == 0) THEN
1113 1 : WRITE (io_unit, '(A,2(1X,ES13.6))') "ROKS fixed-N rotation Schur errors:", &
1114 2 : schur_error, curvature_error
1115 : END IF
1116 :
1117 2 : 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 2 : 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 8 : vectors(:, 1) = [0.8_dp, 0.6_dp, 0.0_dp]
1135 8 : vectors(:, 2) = [-0.6_dp, 0.8_dp, 0.0_dp]
1136 8 : vectors(:, 3) = [0.0_dp, 0.0_dp, 1.0_dp]
1137 2 : 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 176 : 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 176 : 0.0_dp, 0.0_dp, ABS(eigenvalues(3))], [3, 3]), TRANSPOSE(vectors)))
1144 8 : rhs(:, 1) = [0.3_dp, -0.5_dp, 0.7_dp]
1145 8 : rhs(:, 2) = [-0.2_dp, 0.4_dp, 0.1_dp]
1146 2 : CALL qs_ot_symmetric_abs_solve(matrix, rhs, solution, valid)
1147 54 : expected(:, :) = MATMUL(absolute_matrix, solution)
1148 18 : error = MAXVAL(ABS(expected - rhs))
1149 2 : IF (.NOT. valid .OR. error > 2.0E-13_dp) nfail = nfail + 1
1150 :
1151 2 : 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 176 : 0.0_dp, 0.0_dp, eigenvalues(3)], [3, 3]), TRANSPOSE(vectors)))
1155 8 : rhs(:, 1) = [0.3_dp, -0.5_dp, 0.7_dp]
1156 8 : rhs(:, 2) = [-0.2_dp, 0.4_dp, 0.1_dp]
1157 2 : CALL qs_ot_symmetric_abs_solve(matrix, rhs, solution, valid)
1158 2 : expected(:, :) = rhs(:, :)
1159 6 : 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 106 : TRANSPOSE(vectors(:, 1:2))))
1163 72 : error = MAX(error, MAXVAL(ABS(MATMUL(absolute_matrix, solution) - expected)))
1164 6 : IF (.NOT. valid .OR. error > 2.0E-13_dp .OR. MAXVAL(ABS(solution(3, :))) > 2.0E-13_dp) THEN
1165 0 : nfail = nfail + 1
1166 : END IF
1167 2 : IF (mynode == 0) WRITE (io_unit, '(A,1X,ES13.6)') "symmetric absolute solve error:", error
1168 :
1169 2 : 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 2 : 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 2 : reference(:, :) = RESHAPE([4.0_dp, 0.7_dp, 0.7_dp, 2.5_dp], [2, 2])
1186 2 : correction(:, :) = 0.0_dp
1187 2 : CALL qs_ot_projected_response_update(reference, correction, coefficients, valid)
1188 6 : preservation_error = MAXVAL(ABS(coefficients - [1.0_dp, 0.0_dp]))
1189 2 : IF (.NOT. valid .OR. preservation_error > 1.0E-12_dp) nfail = nfail + 1
1190 :
1191 2 : correction(:, :) = RESHAPE([0.8_dp, -0.3_dp, -0.3_dp, 0.4_dp], [2, 2])
1192 2 : gradient(:) = [1.2_dp, -0.4_dp]
1193 : CALL qs_ot_projected_response_update( &
1194 2 : reference, correction, coefficients, valid, projected_gradient=gradient)
1195 14 : total(:, :) = reference + correction
1196 20 : residual(:) = MATMUL(total, coefficients) - gradient
1197 6 : error = MAXVAL(ABS(residual))
1198 2 : IF (.NOT. valid .OR. error > 1.0E-12_dp) nfail = nfail + 1
1199 :
1200 2 : correction(:, :) = RESHAPE([-8.0_dp, 0.0_dp, 0.0_dp, 0.0_dp], [2, 2])
1201 2 : CALL qs_ot_projected_response_update(reference, correction, coefficients, valid)
1202 2 : IF (valid) nfail = nfail + 1
1203 2 : IF (mynode == 0) WRITE (io_unit, '(A,2(1X,ES14.6))') &
1204 1 : "projected response preservation/residual errors:", preservation_error, error
1205 :
1206 2 : 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 2 : 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 2 : matrix(:, :) = 0.0_dp
1224 2 : matrix(1, 1) = 1.0_dp
1225 2 : matrix(2, 2) = 2.0_dp
1226 2 : matrix(3, 3) = 3.0_dp
1227 2 : step(:) = [0.3_dp, -0.4_dp, 0.2_dp]
1228 2 : response(:) = [0.0_dp, -0.3_dp, 0.8_dp]
1229 2 : CALL qs_ot_symmetric_sr1_update(matrix, step, response, updated, valid)
1230 34 : secant_error = MAXVAL(ABS(MATMUL(updated, step) - response))
1231 26 : symmetry_error = MAXVAL(ABS(updated - TRANSPOSE(updated)))
1232 2 : IF (.NOT. valid .OR. secant_error > 2.0E-14_dp .OR. symmetry_error > 2.0E-14_dp) THEN
1233 0 : nfail = nfail + 1
1234 : END IF
1235 :
1236 2 : orthogonal_residual(:) = [0.4_dp, 0.3_dp, 0.0_dp]
1237 34 : response(:) = MATMUL(matrix, step) + orthogonal_residual
1238 2 : CALL qs_ot_symmetric_sr1_update(matrix, step, response, updated, valid)
1239 26 : rejected_error = MAXVAL(ABS(updated - matrix))
1240 2 : IF (valid .OR. rejected_error > 2.0E-14_dp) nfail = nfail + 1
1241 2 : IF (mynode == 0) WRITE (io_unit, '(A,3(1X,ES13.6))') &
1242 1 : "symmetric SR1 secant/symmetry/rejection errors:", &
1243 2 : secant_error, symmetry_error, rejected_error
1244 :
1245 2 : 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 2 : 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 2 : density_step(:, :) = CMPLX(0.0_dp, 0.0_dp, KIND=dp)
1270 2 : density_step(1, 1) = CMPLX(0.30_dp, 0.0_dp, KIND=dp)
1271 2 : density_step(2, 2) = CMPLX(-0.18_dp, 0.0_dp, KIND=dp)
1272 2 : density_step(3, 3) = CMPLX(-0.12_dp, 0.0_dp, KIND=dp)
1273 2 : density_step(1, 2) = CMPLX(0.07_dp, -0.11_dp, KIND=dp)
1274 2 : density_step(2, 1) = CONJG(density_step(1, 2))
1275 2 : density_step(1, 3) = CMPLX(-0.04_dp, 0.05_dp, KIND=dp)
1276 2 : density_step(3, 1) = CONJG(density_step(1, 3))
1277 :
1278 26 : hamiltonian_step(:, :) = -0.45_dp*density_step
1279 2 : hamiltonian_step(2, 3) = CMPLX(0.06_dp, 0.03_dp, KIND=dp)
1280 2 : hamiltonian_step(3, 2) = CONJG(hamiltonian_step(2, 3))
1281 26 : density_modes(:, :, 1) = density_step
1282 26 : density_modes(:, :, 2) = CMPLX(0.0_dp, 0.0_dp, KIND=dp)
1283 2 : density_modes(1, 1, 2) = CMPLX(-0.2_dp, 0.0_dp, KIND=dp)
1284 2 : density_modes(2, 2, 2) = CMPLX(0.1_dp, 0.0_dp, KIND=dp)
1285 2 : density_modes(3, 3, 2) = CMPLX(0.1_dp, 0.0_dp, KIND=dp)
1286 2 : density_modes(1, 3, 2) = CMPLX(0.03_dp, 0.08_dp, KIND=dp)
1287 2 : 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 2 : density_norm, response_work)
1292 2 : density_overlap(1) = density_norm
1293 26 : density_overlap(2) = SUM(REAL(CONJG(density_step)*density_modes(:, :, 2), KIND=dp))
1294 2 : response_overlap(1) = response_work
1295 26 : response_overlap(2) = SUM(REAL(CONJG(hamiltonian_step)*density_modes(:, :, 2), KIND=dp))
1296 : CALL qs_ot_density_secant_projected_hessian( &
1297 2 : density_norm, response_work, density_overlap, response_overlap, correction_overlap, valid)
1298 :
1299 2 : expected(1, 1) = response_work
1300 2 : expected(1, 2) = response_overlap(2)
1301 2 : expected(2, 1) = response_overlap(2)
1302 : expected(2, 2) = (2.0_dp*density_overlap(2)*response_overlap(2) - &
1303 2 : response_work*density_overlap(2)**2/density_norm)/density_norm
1304 14 : direct_error = MAXVAL(ABS(correction - expected))
1305 14 : overlap_error = MAXVAL(ABS(correction_overlap - correction))
1306 2 : 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 2 : calibrated_valid, secant_mode=2, secant_position=secant_position)
1311 : expected(1, 1) = response_work
1312 2 : expected(1, 2) = response_work/secant_position
1313 2 : expected(2, 1) = expected(1, 2)
1314 2 : expected(2, 2) = response_work/secant_position**2
1315 14 : 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 2 : rejected_valid, secant_mode=2, secant_position=0.0_dp)
1319 14 : rejected_error = MAXVAL(ABS(rejected))
1320 :
1321 2 : unitary(:, :) = CMPLX(0.0_dp, 0.0_dp, KIND=dp)
1322 2 : unitary(1, 1) = CMPLX(SQRT(0.5_dp), 0.0_dp, KIND=dp)
1323 2 : unitary(2, 1) = CMPLX(0.0_dp, SQRT(0.5_dp), KIND=dp)
1324 2 : unitary(1, 2) = CMPLX(0.0_dp, SQRT(0.5_dp), KIND=dp)
1325 2 : unitary(2, 2) = CMPLX(SQRT(0.5_dp), 0.0_dp, KIND=dp)
1326 2 : unitary(3, 3) = CMPLX(1.0_dp, 0.0_dp, KIND=dp)
1327 158 : density_step = MATMUL(CONJG(TRANSPOSE(unitary)), MATMUL(density_step, unitary))
1328 158 : hamiltonian_step = MATMUL(CONJG(TRANSPOSE(unitary)), MATMUL(hamiltonian_step, unitary))
1329 : transformed_modes(:, :, 1) = &
1330 186 : MATMUL(CONJG(TRANSPOSE(unitary)), MATMUL(density_modes(:, :, 1), unitary))
1331 : transformed_modes(:, :, 2) = &
1332 186 : MATMUL(CONJG(TRANSPOSE(unitary)), MATMUL(density_modes(:, :, 2), unitary))
1333 : CALL qs_ot_density_secant_hessian( &
1334 2 : density_step, hamiltonian_step, transformed_modes, correction_gauge, valid)
1335 14 : 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 2 : calibrated_error > 2.0E-14_dp .OR. rejected_error > 2.0E-14_dp) nfail = nfail + 1
1341 2 : IF (mynode == 0) WRITE (io_unit, '(A,6(1X,ES13.6))') &
1342 1 : "density secant direct/overlap/gauge/sign/calibrated/rejected errors:", &
1343 2 : direct_error, overlap_error, gauge_error, sign_error, calibrated_error, rejected_error
1344 :
1345 2 : 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 2 : 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 2 : density_step(1, 1, 1) = 0.31_dp
1374 2 : density_step(2, 2, 1) = -0.31_dp
1375 2 : density_step(1, 2, 1) = CMPLX(0.08_dp, -0.04_dp, KIND=dp)
1376 2 : density_step(2, 1, 1) = CONJG(density_step(1, 2, 1))
1377 2 : density_step(1, 1, 2) = -0.19_dp
1378 2 : density_step(2, 2, 2) = 0.19_dp
1379 2 : density_step(1, 2, 2) = CMPLX(-0.05_dp, 0.06_dp, KIND=dp)
1380 2 : density_step(2, 1, 2) = CONJG(density_step(1, 2, 2))
1381 :
1382 14 : hamiltonian_step(:, :, 1) = 0.67_dp*density_step(:, :, 1)
1383 2 : hamiltonian_step(1, 2, 1) = hamiltonian_step(1, 2, 1) + CMPLX(-0.03_dp, 0.09_dp, KIND=dp)
1384 2 : hamiltonian_step(2, 1, 1) = CONJG(hamiltonian_step(1, 2, 1))
1385 14 : hamiltonian_step(:, :, 2) = -0.43_dp*density_step(:, :, 2)
1386 2 : hamiltonian_step(1, 2, 2) = hamiltonian_step(1, 2, 2) + CMPLX(0.07_dp, 0.02_dp, KIND=dp)
1387 2 : 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 2 : density_modes(1, 1, 1, 1) = 0.21_dp
1391 2 : density_modes(2, 2, 1, 1) = -0.21_dp
1392 2 : density_modes(1, 2, 1, 1) = CMPLX(0.04_dp, 0.03_dp, KIND=dp)
1393 2 : density_modes(2, 1, 1, 1) = CONJG(density_modes(1, 2, 1, 1))
1394 2 : density_modes(1, 1, 2, 1) = -0.13_dp
1395 2 : density_modes(2, 2, 2, 1) = 0.13_dp
1396 2 : density_modes(1, 2, 2, 1) = CMPLX(-0.02_dp, 0.05_dp, KIND=dp)
1397 2 : density_modes(2, 1, 2, 1) = CONJG(density_modes(1, 2, 2, 1))
1398 2 : density_modes(1, 1, 1, 2) = -0.17_dp
1399 2 : density_modes(2, 2, 1, 2) = 0.17_dp
1400 2 : density_modes(1, 2, 1, 2) = CMPLX(0.06_dp, -0.01_dp, KIND=dp)
1401 2 : density_modes(2, 1, 1, 2) = CONJG(density_modes(1, 2, 1, 2))
1402 2 : density_modes(1, 1, 2, 2) = 0.09_dp
1403 2 : density_modes(2, 2, 2, 2) = -0.09_dp
1404 2 : density_modes(1, 2, 2, 2) = CMPLX(0.03_dp, 0.04_dp, KIND=dp)
1405 2 : density_modes(2, 1, 2, 2) = CONJG(density_modes(1, 2, 2, 2))
1406 :
1407 2 : direct_density_step(:, :) = CMPLX(0.0_dp, 0.0_dp, KIND=dp)
1408 2 : direct_hamiltonian_step(:, :) = CMPLX(0.0_dp, 0.0_dp, KIND=dp)
1409 2 : direct_density_modes(:, :, :) = CMPLX(0.0_dp, 0.0_dp, KIND=dp)
1410 2 : density_norm = 0.0_dp
1411 2 : response_work = 0.0_dp
1412 2 : density_overlap(:) = 0.0_dp
1413 2 : response_overlap(:) = 0.0_dp
1414 2 : separate_correction(:, :) = 0.0_dp
1415 2 : spin_valid = .TRUE.
1416 6 : DO ispin = 1, nspin
1417 4 : first = (ispin - 1)*n + 1
1418 4 : last = ispin*n
1419 28 : direct_density_step(first:last, first:last) = density_step(:, :, ispin)
1420 28 : direct_hamiltonian_step(first:last, first:last) = hamiltonian_step(:, :, ispin)
1421 60 : direct_density_modes(first:last, first:last, :) = density_modes(:, :, :, ispin)
1422 : density_norm = density_norm + &
1423 : SUM(REAL(CONJG(density_step(:, :, ispin))* &
1424 28 : density_step(:, :, ispin), KIND=dp))
1425 : response_work = response_work + &
1426 : SUM(REAL(CONJG(density_step(:, :, ispin))* &
1427 28 : hamiltonian_step(:, :, ispin), KIND=dp))
1428 12 : DO mode = 1, nmode
1429 : density_overlap(mode) = density_overlap(mode) + &
1430 : SUM(REAL(CONJG(density_step(:, :, ispin))* &
1431 56 : density_modes(:, :, mode, ispin), KIND=dp))
1432 : response_overlap(mode) = response_overlap(mode) + &
1433 : SUM(REAL(CONJG(hamiltonian_step(:, :, ispin))* &
1434 60 : density_modes(:, :, mode, ispin), KIND=dp))
1435 : END DO
1436 : CALL qs_ot_density_secant_hessian( &
1437 : density_step(:, :, ispin), hamiltonian_step(:, :, ispin), &
1438 4 : density_modes(:, :, :, ispin), spin_correction, direct_valid)
1439 4 : spin_valid = spin_valid .AND. direct_valid
1440 30 : 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 2 : direct_correction, direct_valid)
1446 : CALL qs_ot_density_secant_projected_hessian( &
1447 : density_norm, response_work, density_overlap, response_overlap, &
1448 2 : joint_correction, joint_valid)
1449 14 : direct_error = MAXVAL(ABS(joint_correction - direct_correction))
1450 14 : coupling_size = MAXVAL(ABS(joint_correction - separate_correction))
1451 :
1452 : IF (.NOT. direct_valid .OR. .NOT. joint_valid .OR. .NOT. spin_valid .OR. &
1453 2 : direct_error > 3.0E-14_dp .OR. coupling_size < 1.0E-5_dp) nfail = nfail + 1
1454 2 : IF (mynode == 0) WRITE (io_unit, '(A,2(1X,ES13.6))') &
1455 1 : "joint-spin density secant direct error/coupling:", direct_error, coupling_size
1456 :
1457 2 : 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 2 : 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 2 : generator(:, :) = CMPLX(0.0_dp, 0.0_dp, KIND=dp)
1480 2 : generator(1, 2) = CMPLX(0.17_dp, -0.09_dp, KIND=dp)
1481 2 : generator(2, 1) = -CONJG(generator(1, 2))
1482 2 : generator(1, 3) = CMPLX(-0.06_dp, 0.11_dp, KIND=dp)
1483 2 : generator(3, 1) = -CONJG(generator(1, 3))
1484 2 : generator(2, 3) = CMPLX(0.08_dp, 0.04_dp, KIND=dp)
1485 2 : generator(3, 2) = -CONJG(generator(2, 3))
1486 2 : occupation(:) = [1.72_dp, 0.83_dp, 0.31_dp]
1487 2 : rotation_step(:) = [0.21_dp, -0.14_dp, -0.08_dp, 0.17_dp, 0.11_dp, 0.06_dp]
1488 2 : 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 2 : density_tangent, difference_step=fd_step)
1493 :
1494 2 : direction_generator(:, :) = CMPLX(0.0_dp, 0.0_dp, KIND=dp)
1495 2 : r = 0
1496 6 : DO i = 1, nbands - 1
1497 12 : DO j = i + 1, nbands
1498 6 : r = r + 1
1499 6 : direction_generator(i, j) = CMPLX(rotation_step(r), 0.0_dp, KIND=dp)
1500 6 : direction_generator(j, i) = -direction_generator(i, j)
1501 6 : r = r + 1
1502 : direction_generator(i, j) = direction_generator(i, j) + &
1503 6 : CMPLX(0.0_dp, rotation_step(r), KIND=dp)
1504 : direction_generator(j, i) = direction_generator(j, i) + &
1505 10 : CMPLX(0.0_dp, rotation_step(r), KIND=dp)
1506 : END DO
1507 : END DO
1508 2 : rotation = dense_antihermitian_exp(generator)
1509 26 : rotation_plus = dense_antihermitian_exp(generator + fd_step*direction_generator)
1510 26 : rotation_minus = dense_antihermitian_exp(generator - fd_step*direction_generator)
1511 2 : weighted_rotation(:, :) = rotation_plus
1512 8 : DO j = 1, nbands
1513 : weighted_rotation(:, j) = &
1514 : (kpoint_weight*occupation(j) + fd_step*weighted_occupation_step(j))* &
1515 26 : weighted_rotation(:, j)
1516 : END DO
1517 80 : density_plus(:, :) = MATMUL(weighted_rotation, CONJG(TRANSPOSE(rotation_plus)))
1518 2 : weighted_rotation(:, :) = rotation_minus
1519 8 : DO j = 1, nbands
1520 : weighted_rotation(:, j) = &
1521 : (kpoint_weight*occupation(j) - fd_step*weighted_occupation_step(j))* &
1522 26 : weighted_rotation(:, j)
1523 : END DO
1524 80 : density_minus(:, :) = MATMUL(weighted_rotation, CONJG(TRANSPOSE(rotation_minus)))
1525 : density_plus(:, :) = MATMUL(CONJG(TRANSPOSE(rotation)), &
1526 182 : MATMUL((density_plus - density_minus)/(2.0_dp*fd_step), rotation))
1527 26 : error = MAXVAL(ABS(density_tangent - density_plus))
1528 : trace_error = ABS(REAL(SUM([(density_tangent(i, i), i=1, nbands)]), KIND=dp) - &
1529 8 : SUM(weighted_occupation_step))
1530 :
1531 2 : IF (error > 3.0E-9_dp .OR. trace_error > 3.0E-10_dp) nfail = nfail + 1
1532 2 : IF (mynode == 0) WRITE (io_unit, '(A,2(1X,ES13.6))') &
1533 1 : "finite complex density tangent errors:", error, trace_error
1534 :
1535 2 : 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 2 : 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 2 : c0(:, :) = CMPLX(0.0_dp, 0.0_dp, KIND=dp)
1563 2 : c0(1, 1) = 1.0_dp
1564 2 : c0(2, 2) = 1.0_dp
1565 2 : c0(3, 3) = 1.0_dp
1566 2 : c1(:, :) = c0
1567 2 : c1(1, 1) = COS(0.31_dp)
1568 2 : c1(4, 1) = SIN(0.31_dp)
1569 2 : c1(2, 2) = COS(0.23_dp)
1570 2 : c1(3, 2) = CMPLX(0.0_dp, SIN(0.23_dp), KIND=dp)
1571 2 : c1(2, 3) = CMPLX(0.0_dp, SIN(0.23_dp), KIND=dp)
1572 2 : c1(3, 3) = COS(0.23_dp)
1573 2 : occupation0(:) = [1.73_dp, 0.82_dp, 0.82_dp]
1574 2 : occupation1(:) = [1.61_dp, 0.91_dp, 0.76_dp]
1575 :
1576 2 : hamiltonian_step(:, :) = CMPLX(0.0_dp, 0.0_dp, KIND=dp)
1577 2 : hamiltonian_step(1, 1) = 0.21_dp
1578 2 : hamiltonian_step(2, 2) = -0.17_dp
1579 2 : hamiltonian_step(3, 3) = 0.09_dp
1580 2 : hamiltonian_step(4, 4) = -0.04_dp
1581 2 : hamiltonian_step(1, 4) = CMPLX(0.07_dp, -0.03_dp, KIND=dp)
1582 2 : hamiltonian_step(4, 1) = CONJG(hamiltonian_step(1, 4))
1583 2 : hamiltonian_step(2, 3) = CMPLX(-0.05_dp, 0.08_dp, KIND=dp)
1584 2 : hamiltonian_step(3, 2) = CONJG(hamiltonian_step(2, 3))
1585 2 : density_modes(:, :, :) = CMPLX(0.0_dp, 0.0_dp, KIND=dp)
1586 2 : density_modes(1, 1, 1) = 0.13_dp
1587 2 : density_modes(2, 2, 1) = -0.08_dp
1588 2 : density_modes(3, 3, 1) = -0.05_dp
1589 2 : density_modes(1, 2, 2) = CMPLX(0.04_dp, -0.06_dp, KIND=dp)
1590 2 : density_modes(2, 1, 2) = CONJG(density_modes(1, 2, 2))
1591 :
1592 98 : overlap(:, :) = MATMUL(CONJG(TRANSPOSE(c0)), c1)
1593 : ! Keep the reference path independent of optimized complex MATMUL expressions.
1594 2 : CALL reference_project_operator(c0, hamiltonian_step, h0)
1595 2 : 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 2 : response_work, density_overlap, response_overlap, valid)
1599 :
1600 2 : density_step(:, :) = CMPLX(0.0_dp, 0.0_dp, KIND=dp)
1601 8 : 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 176 : SPREAD(CONJG(c0(:, i)), DIM=1, NCOPIES=nao)
1608 : END DO
1609 6 : DO i = 1, nmode
1610 4 : CALL reference_transform_density_mode(c1, density_modes(:, :, i), density_modes_ao(:, :, i))
1611 : density_overlap_ref(i) = &
1612 4 : reference_frobenius_product_real(density_step, density_modes_ao(:, :, i))/weight
1613 : response_overlap_ref(i) = &
1614 6 : reference_frobenius_product_real(hamiltonian_step, density_modes_ao(:, :, i))
1615 : END DO
1616 2 : density_norm_ref = reference_frobenius_product_real(density_step, density_step)/weight
1617 2 : response_work_ref = reference_frobenius_product_real(density_step, hamiltonian_step)
1618 2 : error = MAX(ABS(density_norm - density_norm_ref), ABS(response_work - response_work_ref))
1619 6 : error = MAX(error, MAXVAL(ABS(density_overlap - density_overlap_ref)))
1620 6 : error = MAX(error, MAXVAL(ABS(response_overlap - response_overlap_ref)))
1621 :
1622 2 : IF (.NOT. valid .OR. error > 3.0E-13_dp) nfail = nfail + 1
1623 2 : IF (mynode == 0) WRITE (io_unit, '(A,1X,ES13.6)') &
1624 1 : "moving-subspace density secant error:", error
1625 :
1626 2 : 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 4 : 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 16 : DO j = 1, SIZE(orbitals, 2)
1642 52 : DO i = 1, SIZE(orbitals, 2)
1643 36 : value = CMPLX(0.0_dp, 0.0_dp, KIND=dp)
1644 180 : DO nu = 1, SIZE(orbitals, 1)
1645 756 : DO mu = 1, SIZE(orbitals, 1)
1646 720 : value = value + CONJG(orbitals(mu, i))*hamiltonian(mu, nu)*orbitals(nu, j)
1647 : END DO
1648 : END DO
1649 48 : projected(i, j) = value
1650 : END DO
1651 : END DO
1652 :
1653 4 : 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 4 : 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 20 : DO j = 1, SIZE(orbitals, 1)
1669 84 : DO i = 1, SIZE(orbitals, 1)
1670 64 : value = CMPLX(0.0_dp, 0.0_dp, KIND=dp)
1671 256 : DO n = 1, SIZE(orbitals, 2)
1672 832 : DO m = 1, SIZE(orbitals, 2)
1673 768 : value = value + orbitals(i, m)*mode(m, n)*CONJG(orbitals(j, n))
1674 : END DO
1675 : END DO
1676 80 : transformed(i, j) = value
1677 : END DO
1678 : END DO
1679 :
1680 4 : 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 12 : 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 252 : AIMAG(matrix_a)*AIMAG(matrix_b))
1694 :
1695 12 : 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 6 : 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 6 : energy = 0.0_dp
1714 24 : DO i = 1, SIZE(occupation)
1715 : energy = energy + kpoint_weight*occupation(i)* &
1716 : REAL(DOT_PRODUCT(rotation(:, i), &
1717 330 : MATMUL(hamiltonian, rotation(:, i))), KIND=dp)
1718 : END DO
1719 :
1720 6 : 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 80 : SUBROUTINE fixed_n_fermi_occupations(energies, weights, target_electrons, temperature, &
1733 40 : 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 170 : low = MINVAL(energies) - 80.0_dp*temperature
1743 170 : high = MAXVAL(energies) + 80.0_dp*temperature
1744 6440 : DO iteration = 1, 160
1745 6400 : mu = 0.5_dp*(low + high)
1746 6400 : number = 0.0_dp
1747 27200 : DO i = 1, SIZE(energies)
1748 20800 : occupation(i) = fermi_occupation(energies(i), mu, temperature, maxocc)
1749 27200 : number = number + weights(i)*occupation(i)
1750 : END DO
1751 6440 : IF (number < target_electrons) THEN
1752 : low = mu
1753 : ELSE
1754 3530 : high = mu
1755 : END IF
1756 : END DO
1757 40 : mu = 0.5_dp*(low + high)
1758 170 : DO i = 1, SIZE(energies)
1759 170 : occupation(i) = fermi_occupation(energies(i), mu, temperature, maxocc)
1760 : END DO
1761 :
1762 40 : 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 20930 : 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 20930 : x = (energy - mu)/temperature
1779 20930 : IF (x > 40.0_dp) THEN
1780 44 : occupation = maxocc*EXP(-x)
1781 20886 : ELSE IF (x < -40.0_dp) THEN
1782 46 : occupation = maxocc*(1.0_dp - EXP(x))
1783 : ELSE
1784 20840 : occupation = maxocc/(1.0_dp + EXP(x))
1785 : END IF
1786 :
1787 20930 : 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 22 : 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 44 : REAL(KIND=dp), DIMENSION(SIZE(energies)) :: occupation
1808 :
1809 : CALL fixed_n_fermi_occupations(energies, weights, target_electrons, temperature, &
1810 22 : maxocc, occupation, mu)
1811 22 : value = 0.0_dp
1812 92 : DO i = 1, SIZE(energies)
1813 70 : 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 92 : temperature*maxocc*(p*LOG(p) + (1.0_dp - p)*LOG(1.0_dp - p)))
1816 : END DO
1817 :
1818 22 : 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 4 : SUBROUTINE fixed_n_spectral_projector(hamiltonian, weights, target_electrons, temperature, &
1830 4 : 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), &
1837 8 : SIZE(hamiltonian, 1)) :: vectors, weighted_vectors
1838 : INTEGER :: i
1839 : REAL(KIND=dp) :: mu
1840 8 : REAL(KIND=dp), DIMENSION(SIZE(hamiltonian, 1)) :: eigenvalues, occupation
1841 :
1842 4 : CALL diag_complex(hamiltonian, vectors, eigenvalues)
1843 : CALL fixed_n_fermi_occupations(eigenvalues, weights, target_electrons, temperature, &
1844 4 : maxocc, occupation, mu)
1845 52 : weighted_vectors(:, :) = vectors(:, :)
1846 16 : DO i = 1, SIZE(eigenvalues)
1847 52 : weighted_vectors(:, i) = weights(i)*occupation(i)*weighted_vectors(:, i)
1848 : END DO
1849 216 : projector(:, :) = MATMUL(weighted_vectors, CONJG(TRANSPOSE(vectors)))
1850 :
1851 4 : END SUBROUTINE fixed_n_spectral_projector
1852 :
1853 : ! **************************************************************************************************
1854 : !> \brief Check the real finite rotation pullback used by Gamma-point REF OT.
1855 : !> \param para_env parallel environment
1856 : !> \param nfail accumulated number of failures
1857 : ! **************************************************************************************************
1858 2 : SUBROUTINE test_real_rotation_frechet(para_env, nfail)
1859 : TYPE(mp_para_env_type), POINTER :: para_env
1860 : INTEGER, INTENT(INOUT) :: nfail
1861 :
1862 : INTEGER, PARAMETER :: m = 3
1863 :
1864 : INTEGER :: i
1865 4 : INTEGER, DIMENSION(:), POINTER :: col_dist, col_size, row_dist, row_size
1866 2 : INTEGER, DIMENSION(:, :), POINTER :: pgrid
1867 : LOGICAL :: found
1868 : REAL(KIND=dp) :: error, fd_slope, hstep, predicted
1869 2 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: block
1870 : REAL(KIND=dp), DIMENSION(m), TARGET :: rotation_evals
1871 : REAL(KIND=dp), DIMENSION(m, m) :: dedu, direction, generator, gradient, &
1872 : rotation_minus, rotation_plus
1873 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
1874 : TYPE(dbcsr_distribution_type) :: dist
1875 16 : TYPE(dbcsr_type), DIMENSION(7), TARGET :: matrices
1876 : TYPE(mp_para_env_type), POINTER :: local_para_env
1877 324 : TYPE(qs_ot_type) :: rotation_env
1878 :
1879 2 : NULLIFY (blacs_env, local_para_env)
1880 2 : generator(:, :) = 0.0_dp
1881 2 : generator(1, 2) = 0.23_dp
1882 2 : generator(2, 1) = -generator(1, 2)
1883 2 : generator(1, 3) = -0.11_dp
1884 2 : generator(3, 1) = -generator(1, 3)
1885 2 : generator(2, 3) = 0.17_dp
1886 2 : generator(3, 2) = -generator(2, 3)
1887 :
1888 2 : direction(:, :) = 0.0_dp
1889 2 : direction(1, 2) = -0.09_dp
1890 2 : direction(2, 1) = -direction(1, 2)
1891 2 : direction(1, 3) = 0.14_dp
1892 2 : direction(3, 1) = -direction(1, 3)
1893 2 : direction(2, 3) = 0.06_dp
1894 2 : direction(3, 2) = -direction(2, 3)
1895 :
1896 : dedu = RESHAPE([0.31_dp, -0.27_dp, 0.08_dp, &
1897 : 0.19_dp, 0.42_dp, -0.16_dp, &
1898 2 : -0.07_dp, 0.21_dp, 0.35_dp], [m, m])
1899 :
1900 0 : ALLOCATE (local_para_env)
1901 2 : CALL local_para_env%from_dup(mp_comm_self)
1902 2 : CALL cp_blacs_env_create(blacs_env=blacs_env, para_env=local_para_env)
1903 2 : ALLOCATE (pgrid(0:0, 0:0), row_dist(1), col_dist(1), row_size(1), col_size(1))
1904 6 : pgrid(:, :) = 0
1905 4 : row_dist(:) = 0
1906 4 : col_dist(:) = 0
1907 4 : row_size(:) = m
1908 4 : col_size(:) = m
1909 : CALL dbcsr_distribution_new(dist, group=mp_comm_self%get_handle(), pgrid=pgrid, &
1910 2 : row_dist=row_dist, col_dist=col_dist)
1911 16 : DO i = 1, SIZE(matrices)
1912 16 : CALL create_one_block(matrices(i), 'real_rotation', dist, row_size, col_size)
1913 : END DO
1914 2 : CALL dbcsr_put_block(matrices(2), 1, 1, generator)
1915 2 : CALL dbcsr_put_block(matrices(3), 1, 1, dedu)
1916 :
1917 2 : rotation_env%para_env => local_para_env
1918 2 : rotation_env%blacs_env => blacs_env
1919 2 : rotation_env%rot_mat_u => matrices(1)
1920 2 : rotation_env%rot_mat_x => matrices(2)
1921 2 : rotation_env%rot_mat_dedu => matrices(3)
1922 2 : rotation_env%rot_mat_gx => matrices(4)
1923 2 : rotation_env%rot_mat_evec_re => matrices(5)
1924 2 : rotation_env%rot_mat_evec_im => matrices(6)
1925 2 : rotation_env%matrix_buf1 => matrices(7)
1926 2 : rotation_env%rot_mat_evals => rotation_evals
1927 :
1928 2 : CALL qs_ot_generate_rotation(rotation_env)
1929 2 : CALL qs_ot_rot_mat_derivative(rotation_env)
1930 2 : gradient(:, :) = 0.0_dp
1931 2 : CALL dbcsr_get_block_p(matrices(4), 1, 1, block, found)
1932 26 : IF (found) gradient = block
1933 :
1934 2 : hstep = 1.0E-6_dp
1935 : rotation_plus = REAL(dense_antihermitian_exp( &
1936 50 : CMPLX(generator + hstep*direction, 0.0_dp, KIND=dp)), KIND=dp)
1937 : rotation_minus = REAL(dense_antihermitian_exp( &
1938 50 : CMPLX(generator - hstep*direction, 0.0_dp, KIND=dp)), KIND=dp)
1939 26 : fd_slope = SUM(dedu*(rotation_plus - rotation_minus))/(2.0_dp*hstep)
1940 26 : predicted = 0.5_dp*SUM(gradient*direction)
1941 2 : error = ABS(predicted - fd_slope)
1942 2 : IF (error > 5.0E-8_dp) nfail = nfail + 1
1943 2 : IF (para_env%is_source()) THEN
1944 1 : WRITE (io_unit, '(A,3(1X,ES13.6))') 'real finite rotation slope fd/predicted/error', &
1945 2 : fd_slope, predicted, error
1946 : END IF
1947 :
1948 16 : DO i = 1, SIZE(matrices)
1949 16 : CALL dbcsr_release(matrices(i))
1950 : END DO
1951 2 : CALL dbcsr_distribution_release(dist)
1952 2 : DEALLOCATE (pgrid, row_dist, col_dist, row_size, col_size)
1953 2 : CALL cp_blacs_env_release(blacs_env)
1954 2 : CALL mp_para_env_release(local_para_env)
1955 :
1956 12 : END SUBROUTINE test_real_rotation_frechet
1957 :
1958 : ! **************************************************************************************************
1959 : !> \brief Check the combined finite complex REF and rotation derivatives.
1960 : !> \param para_env parallel environment
1961 : !> \param nfail accumulated number of failures
1962 : ! **************************************************************************************************
1963 2 : SUBROUTINE test_complex_ref_rotation(para_env, nfail)
1964 : TYPE(mp_para_env_type), POINTER :: para_env
1965 : INTEGER, INTENT(INOUT) :: nfail
1966 :
1967 : COMPLEX(KIND=dp), DIMENSION(k, k) :: direction_a, generator, gradient_a, &
1968 : inverse_sqrt, rotation
1969 : COMPLEX(KIND=dp), DIMENSION(n, k) :: c_dense, c_reference, direction_x, &
1970 : gradient_x, hc_dense, q_dense, x_dense
1971 : COMPLEX(KIND=dp), DIMENSION(n, n) :: hamiltonian, overlap
1972 : INTEGER :: i, j
1973 4 : INTEGER, DIMENSION(:), POINTER :: col_dist, col_k, col_n, row_dist, row_k, &
1974 2 : row_n
1975 2 : INTEGER, DIMENSION(:, :), POINTER :: pgrid
1976 : REAL(KIND=dp) :: error_c, error_ref, error_rot, fd_ref, &
1977 : fd_rot, hstep, predicted_ref, &
1978 : predicted_rot
1979 : REAL(KIND=dp), DIMENSION(k) :: occupation
1980 : REAL(KIND=dp), DIMENSION(k), TARGET :: rotation_evals
1981 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
1982 : TYPE(dbcsr_distribution_type) :: dist_kk, dist_nk, dist_nn
1983 22 : TYPE(dbcsr_type), DIMENSION(10), TARGET :: rot_matrix
1984 6 : TYPE(dbcsr_type), DIMENSION(2), TARGET :: nk_buffer
1985 8 : TYPE(dbcsr_type), DIMENSION(3), TARGET :: nonsym_buffer
1986 10 : TYPE(dbcsr_type), DIMENSION(4), TARGET :: sym_buffer
1987 : TYPE(dbcsr_type), POINTER :: c_im_p, c_re_p, hc_im_p, hc_re_p, &
1988 : s_im_p, s_re_p
1989 : TYPE(dbcsr_type), TARGET :: c_im, c_re, f_im, f_re, g_im, g_re, &
1990 : hc_im, hc_re, s_im, s_re, sx_im, &
1991 : sx_re, x_im, x_re
1992 : TYPE(mp_para_env_type), POINTER :: local_para_env
1993 324 : TYPE(qs_ot_type) :: rotation_env
1994 :
1995 2 : NULLIFY (blacs_env, c_im_p, c_re_p, hc_im_p, hc_re_p, local_para_env, &
1996 2 : s_im_p, s_re_p)
1997 :
1998 : x_dense(:, 1) = [CMPLX(1.05_dp, 0.12_dp, KIND=dp), CMPLX(0.18_dp, -0.27_dp, KIND=dp), &
1999 8 : CMPLX(-0.14_dp, 0.16_dp, KIND=dp)]
2000 : x_dense(:, 2) = [CMPLX(0.26_dp, 0.31_dp, KIND=dp), CMPLX(0.94_dp, -0.17_dp, KIND=dp), &
2001 8 : CMPLX(0.21_dp, 0.13_dp, KIND=dp)]
2002 : direction_x(:, 1) = [CMPLX(0.11_dp, -0.07_dp, KIND=dp), CMPLX(-0.16_dp, 0.09_dp, KIND=dp), &
2003 8 : CMPLX(0.05_dp, 0.14_dp, KIND=dp)]
2004 : direction_x(:, 2) = [CMPLX(-0.08_dp, 0.15_dp, KIND=dp), CMPLX(0.13_dp, 0.04_dp, KIND=dp), &
2005 8 : CMPLX(-0.10_dp, -0.06_dp, KIND=dp)]
2006 :
2007 : generator(:, :) = CMPLX(0.0_dp, 0.0_dp, KIND=dp)
2008 2 : generator(1, 1) = CMPLX(0.0_dp, 0.19_dp, KIND=dp)
2009 2 : generator(2, 2) = CMPLX(0.0_dp, -0.08_dp, KIND=dp)
2010 2 : generator(1, 2) = CMPLX(0.23_dp, -0.11_dp, KIND=dp)
2011 2 : generator(2, 1) = -CONJG(generator(1, 2))
2012 : direction_a(:, :) = CMPLX(0.0_dp, 0.0_dp, KIND=dp)
2013 2 : direction_a(1, 1) = CMPLX(0.0_dp, -0.07_dp, KIND=dp)
2014 2 : direction_a(2, 2) = CMPLX(0.0_dp, 0.05_dp, KIND=dp)
2015 2 : direction_a(1, 2) = CMPLX(-0.09_dp, 0.06_dp, KIND=dp)
2016 2 : direction_a(2, 1) = -CONJG(direction_a(1, 2))
2017 :
2018 2 : hamiltonian(:, :) = CMPLX(0.0_dp, 0.0_dp, KIND=dp)
2019 2 : hamiltonian(1, 1) = CMPLX(0.65_dp, 0.0_dp, KIND=dp)
2020 2 : hamiltonian(2, 2) = CMPLX(-0.35_dp, 0.0_dp, KIND=dp)
2021 2 : hamiltonian(3, 3) = CMPLX(1.15_dp, 0.0_dp, KIND=dp)
2022 2 : hamiltonian(1, 2) = CMPLX(0.18_dp, 0.12_dp, KIND=dp)
2023 2 : hamiltonian(1, 3) = CMPLX(-0.09_dp, 0.07_dp, KIND=dp)
2024 2 : hamiltonian(2, 3) = CMPLX(0.14_dp, -0.17_dp, KIND=dp)
2025 8 : DO i = 1, n
2026 14 : DO j = i + 1, n
2027 12 : hamiltonian(j, i) = CONJG(hamiltonian(i, j))
2028 : END DO
2029 : END DO
2030 2 : overlap(:, :) = CMPLX(0.0_dp, 0.0_dp, KIND=dp)
2031 8 : DO i = 1, n
2032 8 : overlap(i, i) = CMPLX(1.0_dp, 0.0_dp, KIND=dp)
2033 : END DO
2034 2 : occupation = [1.0_dp, 0.37_dp]
2035 :
2036 38 : CALL dense_inverse_sqrt(MATMUL(CONJG(TRANSPOSE(x_dense)), x_dense), inverse_sqrt)
2037 38 : q_dense = MATMUL(x_dense, inverse_sqrt)
2038 2 : rotation = dense_antihermitian_exp(generator)
2039 38 : c_reference = MATMUL(q_dense, rotation)
2040 :
2041 2 : ALLOCATE (local_para_env)
2042 2 : CALL local_para_env%from_dup(mp_comm_self)
2043 2 : CALL cp_blacs_env_create(blacs_env=blacs_env, para_env=local_para_env)
2044 :
2045 : ALLOCATE (pgrid(0:0, 0:0), row_dist(1), col_dist(1), row_n(1), row_k(1), &
2046 2 : col_n(1), col_k(1))
2047 6 : pgrid(:, :) = 0
2048 4 : row_dist(:) = 0
2049 4 : col_dist(:) = 0
2050 4 : row_n(:) = n
2051 4 : col_n(:) = n
2052 4 : row_k(:) = k
2053 4 : col_k(:) = k
2054 : CALL dbcsr_distribution_new(dist_nk, group=mp_comm_self%get_handle(), pgrid=pgrid, &
2055 2 : row_dist=row_dist, col_dist=col_dist)
2056 : CALL dbcsr_distribution_new(dist_nn, group=mp_comm_self%get_handle(), pgrid=pgrid, &
2057 2 : row_dist=row_dist, col_dist=col_dist)
2058 : CALL dbcsr_distribution_new(dist_kk, group=mp_comm_self%get_handle(), pgrid=pgrid, &
2059 2 : row_dist=row_dist, col_dist=col_dist)
2060 :
2061 2 : CALL create_one_block(x_re, 'rotation_ref_x_re', dist_nk, row_n, col_k)
2062 2 : CALL create_one_block(x_im, 'rotation_ref_x_im', dist_nk, row_n, col_k)
2063 2 : CALL create_one_block(sx_re, 'rotation_ref_sx_re', dist_nk, row_n, col_k)
2064 2 : CALL create_one_block(sx_im, 'rotation_ref_sx_im', dist_nk, row_n, col_k)
2065 2 : CALL create_one_block(c_re, 'rotation_ref_c_re', dist_nk, row_n, col_k)
2066 2 : CALL create_one_block(c_im, 'rotation_ref_c_im', dist_nk, row_n, col_k)
2067 2 : CALL create_one_block(hc_re, 'rotation_ref_hc_re', dist_nk, row_n, col_k)
2068 2 : CALL create_one_block(hc_im, 'rotation_ref_hc_im', dist_nk, row_n, col_k)
2069 2 : CALL create_one_block(g_re, 'rotation_ref_g_re', dist_nk, row_n, col_k)
2070 2 : CALL create_one_block(g_im, 'rotation_ref_g_im', dist_nk, row_n, col_k)
2071 2 : CALL create_one_block(s_re, 'rotation_ref_s_re', dist_nn, row_n, col_n)
2072 2 : CALL create_one_block(s_im, 'rotation_ref_s_im', dist_nn, row_n, col_n)
2073 2 : CALL create_one_block(f_re, 'rotation_ref_f_re', dist_kk, row_k, col_k)
2074 2 : CALL create_one_block(f_im, 'rotation_ref_f_im', dist_kk, row_k, col_k)
2075 10 : DO i = 1, SIZE(sym_buffer)
2076 10 : CALL create_one_block(sym_buffer(i), 'rotation_ref_sym', dist_kk, row_k, col_k)
2077 : END DO
2078 8 : DO i = 1, SIZE(nonsym_buffer)
2079 8 : CALL create_one_block(nonsym_buffer(i), 'rotation_ref_nonsym', dist_kk, row_k, col_k)
2080 : END DO
2081 6 : DO i = 1, SIZE(nk_buffer)
2082 6 : CALL create_one_block(nk_buffer(i), 'rotation_ref_nk', dist_nk, row_n, col_k)
2083 : END DO
2084 22 : DO i = 1, SIZE(rot_matrix)
2085 22 : CALL create_one_block(rot_matrix(i), 'rotation_ref_rot', dist_kk, row_k, col_k)
2086 : END DO
2087 :
2088 2 : CALL put_complex_pair(x_re, x_im, x_dense)
2089 2 : CALL put_complex_pair(s_re, s_im, overlap)
2090 2 : CALL put_complex_pair(rot_matrix(3), rot_matrix(4), generator)
2091 :
2092 2 : rotation_env%has_complex_kpoint_state = .TRUE.
2093 2 : rotation_env%settings%do_rotation = .TRUE.
2094 2 : rotation_env%settings%eps_irac_filter_matrix = 0.0_dp
2095 2 : rotation_env%settings%ot_method = 'LBFG'
2096 2 : rotation_env%para_env => local_para_env
2097 2 : rotation_env%blacs_env => blacs_env
2098 2 : rotation_env%matrix_x => x_re
2099 2 : rotation_env%matrix_x_im => x_im
2100 2 : rotation_env%matrix_sx => sx_re
2101 2 : rotation_env%matrix_sx_im => sx_im
2102 2 : rotation_env%matrix_gx => g_re
2103 2 : rotation_env%matrix_gx_im => g_im
2104 2 : rotation_env%matrix_ref_inv_sqrt => f_re
2105 2 : rotation_env%matrix_ref_inv_sqrt_im => f_im
2106 2 : rotation_env%buf1_k_k_sym => sym_buffer(1)
2107 2 : rotation_env%buf2_k_k_sym => sym_buffer(2)
2108 2 : rotation_env%buf3_k_k_sym => sym_buffer(3)
2109 2 : rotation_env%buf4_k_k_sym => sym_buffer(4)
2110 2 : rotation_env%buf1_k_k_nosym => nonsym_buffer(1)
2111 2 : rotation_env%buf2_k_k_nosym => nonsym_buffer(2)
2112 2 : rotation_env%buf3_k_k_nosym => nonsym_buffer(3)
2113 2 : rotation_env%buf1_n_k => nk_buffer(1)
2114 2 : rotation_env%buf1_n_k_dp => nk_buffer(2)
2115 2 : rotation_env%rot_mat_u => rot_matrix(1)
2116 2 : rotation_env%rot_mat_u_im => rot_matrix(2)
2117 2 : rotation_env%rot_mat_x => rot_matrix(3)
2118 2 : rotation_env%rot_mat_x_im => rot_matrix(4)
2119 2 : rotation_env%rot_mat_dedu => rot_matrix(5)
2120 2 : rotation_env%rot_mat_dedu_im => rot_matrix(6)
2121 2 : rotation_env%rot_mat_gx => rot_matrix(7)
2122 2 : rotation_env%rot_mat_gx_im => rot_matrix(8)
2123 2 : rotation_env%rot_mat_evec_re => rot_matrix(9)
2124 2 : rotation_env%rot_mat_evec_im => rot_matrix(10)
2125 2 : rotation_env%rot_mat_evals => rotation_evals
2126 :
2127 2 : c_re_p => c_re
2128 2 : c_im_p => c_im
2129 2 : s_re_p => s_re
2130 2 : s_im_p => s_im
2131 2 : CALL qs_ot_get_orbitals_ref_complex(c_re_p, c_im_p, s_re_p, s_im_p, rotation_env)
2132 2 : CALL get_complex_pair(c_re, c_im, c_dense)
2133 18 : error_c = MAXVAL(ABS(c_dense - c_reference))
2134 :
2135 70 : hc_dense = 2.0_dp*MATMUL(hamiltonian, c_dense)
2136 6 : DO j = 1, k
2137 18 : hc_dense(:, j) = occupation(j)*hc_dense(:, j)
2138 : END DO
2139 2 : CALL put_complex_pair(hc_re, hc_im, hc_dense)
2140 2 : hc_re_p => hc_re
2141 2 : hc_im_p => hc_im
2142 2 : CALL qs_ot_get_derivative_ref_complex(hc_re_p, hc_im_p, rotation_env)
2143 2 : CALL get_complex_pair(g_re, g_im, gradient_x)
2144 2 : CALL get_complex_pair(rot_matrix(7), rot_matrix(8), gradient_a)
2145 :
2146 2 : hstep = 1.0E-6_dp
2147 : fd_ref = (rotated_polar_energy(x_dense + hstep*direction_x, generator, &
2148 : hamiltonian, occupation) - &
2149 : rotated_polar_energy(x_dense - hstep*direction_x, generator, &
2150 34 : hamiltonian, occupation))/(2.0_dp*hstep)
2151 18 : predicted_ref = REAL(SUM(CONJG(gradient_x)*direction_x), KIND=dp)
2152 2 : error_ref = ABS(predicted_ref - fd_ref)
2153 :
2154 : fd_rot = (rotated_polar_energy(x_dense, generator + hstep*direction_a, &
2155 : hamiltonian, occupation) - &
2156 : rotated_polar_energy(x_dense, generator - hstep*direction_a, &
2157 26 : hamiltonian, occupation))/(2.0_dp*hstep)
2158 : predicted_rot = 0.5_dp*(SUM(REAL(gradient_a, KIND=dp)*REAL(direction_a, KIND=dp)) + &
2159 26 : SUM(AIMAG(gradient_a)*AIMAG(direction_a)))
2160 2 : error_rot = ABS(predicted_rot - fd_rot)
2161 :
2162 2 : IF (error_c > 5.0E-11_dp .OR. error_ref > 5.0E-8_dp .OR. &
2163 0 : error_rot > 5.0E-8_dp) nfail = nfail + 1
2164 2 : IF (para_env%is_source()) THEN
2165 1 : WRITE (io_unit, '(A,1X,ES13.6)') 'complex finite REF rotation orbital error', error_c
2166 1 : WRITE (io_unit, '(A,3(1X,ES13.6))') 'complex finite REF slope fd/predicted/error', &
2167 2 : fd_ref, predicted_ref, error_ref
2168 1 : WRITE (io_unit, '(A,3(1X,ES13.6))') 'complex finite rotation slope fd/predicted/error', &
2169 2 : fd_rot, predicted_rot, error_rot
2170 : END IF
2171 :
2172 22 : DO i = 1, SIZE(rot_matrix)
2173 22 : CALL dbcsr_release(rot_matrix(i))
2174 : END DO
2175 6 : DO i = 1, SIZE(nk_buffer)
2176 6 : CALL dbcsr_release(nk_buffer(i))
2177 : END DO
2178 8 : DO i = 1, SIZE(nonsym_buffer)
2179 8 : CALL dbcsr_release(nonsym_buffer(i))
2180 : END DO
2181 10 : DO i = 1, SIZE(sym_buffer)
2182 10 : CALL dbcsr_release(sym_buffer(i))
2183 : END DO
2184 2 : CALL dbcsr_release(f_im)
2185 2 : CALL dbcsr_release(f_re)
2186 2 : CALL dbcsr_release(s_im)
2187 2 : CALL dbcsr_release(s_re)
2188 2 : CALL dbcsr_release(g_im)
2189 2 : CALL dbcsr_release(g_re)
2190 2 : CALL dbcsr_release(hc_im)
2191 2 : CALL dbcsr_release(hc_re)
2192 2 : CALL dbcsr_release(c_im)
2193 2 : CALL dbcsr_release(c_re)
2194 2 : CALL dbcsr_release(sx_im)
2195 2 : CALL dbcsr_release(sx_re)
2196 2 : CALL dbcsr_release(x_im)
2197 2 : CALL dbcsr_release(x_re)
2198 2 : CALL dbcsr_distribution_release(dist_kk)
2199 2 : CALL dbcsr_distribution_release(dist_nn)
2200 2 : CALL dbcsr_distribution_release(dist_nk)
2201 2 : DEALLOCATE (pgrid, row_dist, col_dist, row_n, row_k, col_n, col_k)
2202 2 : CALL cp_blacs_env_release(blacs_env)
2203 2 : CALL mp_para_env_release(local_para_env)
2204 :
2205 14 : END SUBROUTINE test_complex_ref_rotation
2206 :
2207 : ! **************************************************************************************************
2208 : !> \brief Check the finite complex rotation Frechet derivative against central differences.
2209 : !> \param para_env parallel environment
2210 : !> \param nfail accumulated number of failures
2211 : ! **************************************************************************************************
2212 2 : SUBROUTINE test_complex_rotation_frechet(para_env, nfail)
2213 : TYPE(mp_para_env_type), POINTER :: para_env
2214 : INTEGER, INTENT(INOUT) :: nfail
2215 :
2216 : INTEGER, PARAMETER :: m = 3
2217 :
2218 : COMPLEX(KIND=dp), DIMENSION(m, m) :: direction, fd, generator, predicted
2219 : INTEGER :: i, j
2220 : REAL(KIND=dp) :: error, hstep
2221 : REAL(KIND=dp), DIMENSION(m) :: evals
2222 :
2223 2 : evals = [-0.7_dp, 0.2000003_dp, 1.1_dp]
2224 2 : generator(:, :) = CMPLX(0.0_dp, 0.0_dp, KIND=dp)
2225 8 : DO i = 1, m
2226 8 : generator(i, i) = CMPLX(0.0_dp, -evals(i), KIND=dp)
2227 : END DO
2228 : direction(:, :) = CMPLX(0.0_dp, 0.0_dp, KIND=dp)
2229 2 : direction(1, 2) = CMPLX(0.17_dp, -0.09_dp, KIND=dp)
2230 2 : direction(2, 1) = -CONJG(direction(1, 2))
2231 2 : direction(1, 3) = CMPLX(-0.06_dp, 0.13_dp, KIND=dp)
2232 2 : direction(3, 1) = -CONJG(direction(1, 3))
2233 2 : direction(2, 3) = CMPLX(0.11_dp, 0.04_dp, KIND=dp)
2234 2 : direction(3, 2) = -CONJG(direction(2, 3))
2235 2 : direction(1, 1) = CMPLX(0.0_dp, 0.08_dp, KIND=dp)
2236 2 : direction(2, 2) = CMPLX(0.0_dp, -0.03_dp, KIND=dp)
2237 2 : direction(3, 3) = CMPLX(0.0_dp, 0.05_dp, KIND=dp)
2238 :
2239 8 : DO i = 1, m
2240 26 : DO j = 1, m
2241 : predicted(i, j) = qs_ot_complex_exp_frechet_kernel(evals(i), evals(j))* &
2242 24 : direction(i, j)
2243 : END DO
2244 : END DO
2245 2 : hstep = 1.0E-6_dp
2246 : fd = (dense_antihermitian_exp(generator + hstep*direction) - &
2247 74 : dense_antihermitian_exp(generator - hstep*direction))/(2.0_dp*hstep)
2248 26 : error = MAXVAL(ABS(fd - predicted))
2249 2 : IF (error > 5.0E-9_dp) nfail = nfail + 1
2250 2 : IF (para_env%is_source()) THEN
2251 1 : WRITE (io_unit, '(A,1X,ES13.6)') 'complex rotation Frechet finite-difference error', error
2252 : END IF
2253 :
2254 2 : END SUBROUTINE test_complex_rotation_frechet
2255 :
2256 : ! **************************************************************************************************
2257 : !> \brief Check union handling for distinct real/imaginary DBCSR block patterns.
2258 : !> \param para_env parallel environment
2259 : !> \param nfail accumulated number of failures
2260 : ! **************************************************************************************************
2261 2 : SUBROUTINE test_sparse_frechet_patterns(para_env, nfail)
2262 : TYPE(mp_para_env_type), POINTER :: para_env
2263 : INTEGER, INTENT(INOUT) :: nfail
2264 :
2265 : COMPLEX(KIND=dp) :: expected, value
2266 2 : INTEGER, DIMENSION(:), POINTER :: col_dist, col_size, row_dist, row_size
2267 : INTEGER, DIMENSION(:, :), POINTER :: pgrid
2268 : LOGICAL :: found
2269 : REAL(KIND=dp) :: error
2270 : REAL(KIND=dp), DIMENSION(2) :: evals
2271 2 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: block
2272 : TYPE(dbcsr_distribution_type) :: dist
2273 : TYPE(dbcsr_type) :: inner_im, inner_re, outer_im, outer_re
2274 :
2275 2 : ALLOCATE (pgrid(0:0, 0:0), row_dist(2), col_dist(2), row_size(2), col_size(2))
2276 6 : pgrid(:, :) = 0
2277 6 : row_dist(:) = 0
2278 6 : col_dist(:) = 0
2279 6 : row_size(:) = 1
2280 6 : col_size(:) = 1
2281 : CALL dbcsr_distribution_new(dist, group=mp_comm_self%get_handle(), pgrid=pgrid, &
2282 2 : row_dist=row_dist, col_dist=col_dist)
2283 2 : CALL dbcsr_create(inner_re, 'frechet_sparse_re', dist, dbcsr_type_no_symmetry, row_size, col_size)
2284 2 : CALL dbcsr_reserve_blocks(inner_re, rows=[1], cols=[2])
2285 2 : CALL dbcsr_finalize(inner_re)
2286 2 : CALL dbcsr_create(inner_im, 'frechet_sparse_im', dist, dbcsr_type_no_symmetry, row_size, col_size)
2287 2 : CALL dbcsr_reserve_blocks(inner_im, rows=[2], cols=[1])
2288 2 : CALL dbcsr_finalize(inner_im)
2289 2 : CALL dbcsr_put_block(inner_re, 1, 2, RESHAPE([0.37_dp], [1, 1]))
2290 2 : CALL dbcsr_put_block(inner_im, 2, 1, RESHAPE([-0.21_dp], [1, 1]))
2291 :
2292 2 : evals = [-0.4_dp, 0.9_dp]
2293 2 : CALL qs_ot_apply_complex_frechet_dbcsr(evals, inner_re, inner_im, outer_re, outer_im)
2294 : error = 0.0_dp
2295 2 : expected = qs_ot_complex_exp_frechet_kernel(evals(1), evals(2))*CMPLX(0.37_dp, 0.0_dp, KIND=dp)
2296 2 : value = CMPLX(0.0_dp, 0.0_dp, KIND=dp)
2297 2 : CALL dbcsr_get_block_p(outer_re, 1, 2, block, found)
2298 2 : IF (found) value = value + CMPLX(block(1, 1), 0.0_dp, KIND=dp)
2299 2 : CALL dbcsr_get_block_p(outer_im, 1, 2, block, found)
2300 2 : IF (found) value = value + CMPLX(0.0_dp, block(1, 1), KIND=dp)
2301 2 : error = MAX(error, ABS(value - expected))
2302 :
2303 2 : expected = qs_ot_complex_exp_frechet_kernel(evals(2), evals(1))*CMPLX(0.0_dp, -0.21_dp, KIND=dp)
2304 2 : value = CMPLX(0.0_dp, 0.0_dp, KIND=dp)
2305 2 : CALL dbcsr_get_block_p(outer_re, 2, 1, block, found)
2306 2 : IF (found) value = value + CMPLX(block(1, 1), 0.0_dp, KIND=dp)
2307 2 : CALL dbcsr_get_block_p(outer_im, 2, 1, block, found)
2308 2 : IF (found) value = value + CMPLX(0.0_dp, block(1, 1), KIND=dp)
2309 2 : error = MAX(error, ABS(value - expected))
2310 2 : IF (error > 5.0E-14_dp) nfail = nfail + 1
2311 2 : IF (para_env%is_source()) THEN
2312 1 : WRITE (io_unit, '(A,1X,ES13.6)') 'sparse complex Frechet block-union error', error
2313 : END IF
2314 :
2315 2 : CALL dbcsr_release(outer_im)
2316 2 : CALL dbcsr_release(outer_re)
2317 2 : CALL dbcsr_release(inner_im)
2318 2 : CALL dbcsr_release(inner_re)
2319 2 : CALL dbcsr_distribution_release(dist)
2320 2 : DEALLOCATE (pgrid, row_dist, col_dist, row_size, col_size)
2321 :
2322 2 : END SUBROUTINE test_sparse_frechet_patterns
2323 : ! **************************************************************************************************
2324 : !> \brief Check the finite complex STRICT chart and its analytic derivative.
2325 : !> \param dist_nk distribution for orbital-shaped matrices
2326 : !> \param dist_kk distribution for occupied-space matrices
2327 : !> \param row_n AO block sizes
2328 : !> \param row_k occupied-space row block sizes
2329 : !> \param col_k occupied-space column block sizes
2330 : !> \param para_env parallel environment
2331 : !> \param nfail accumulated number of failures
2332 : ! **************************************************************************************************
2333 2 : SUBROUTINE test_complex_strict_derivative(dist_nk, dist_kk, row_n, row_k, col_k, para_env, nfail)
2334 : TYPE(dbcsr_distribution_type), INTENT(IN) :: dist_nk, dist_kk
2335 : INTEGER, DIMENSION(:), INTENT(INOUT), POINTER :: row_n, row_k, col_k
2336 : TYPE(mp_para_env_type), POINTER :: para_env
2337 : INTEGER, INTENT(INOUT) :: nfail
2338 :
2339 : CHARACTER(LEN=32) :: name
2340 : COMPLEX(KIND=dp), DIMENSION(k, k) :: direction_a, generator, gradient_a, &
2341 : inverse_sqrt, rotation
2342 : COMPLEX(KIND=dp), DIMENSION(n, k) :: c0, direction, gradient, hc, orbitals, &
2343 : orbitals_base, sc0, sx, x
2344 : COMPLEX(KIND=dp), DIMENSION(n, n) :: hamiltonian, overlap
2345 : INTEGER :: i
2346 : LOGICAL :: found
2347 : REAL(KIND=dp) :: error, error_rot, fd_slope, &
2348 : fd_slope_rot, hstep, ortho_error, &
2349 : slope, slope_rot
2350 2 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: block
2351 : REAL(KIND=dp), DIMENSION(k) :: occupation
2352 : REAL(KIND=dp), DIMENSION(k), TARGET :: rotation_evals
2353 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
2354 22 : TYPE(dbcsr_type), DIMENSION(10), TARGET :: rot_matrix
2355 32 : TYPE(dbcsr_type), DIMENSION(15), TARGET :: nk
2356 52 : TYPE(dbcsr_type), DIMENSION(25), TARGET :: kk
2357 : TYPE(dbcsr_type), POINTER :: hc_im_p, hc_re_p, sx_im_p, sx_re_p, &
2358 : x_im_p, x_re_p
2359 : TYPE(mp_para_env_type), POINTER :: strict_para_env
2360 324 : TYPE(qs_ot_type) :: strict_env
2361 :
2362 2 : NULLIFY (blacs_env, hc_im_p, hc_re_p, strict_para_env, x_im_p, x_re_p, sx_im_p, sx_re_p)
2363 2 : ALLOCATE (strict_para_env)
2364 2 : CALL strict_para_env%from_dup(mp_comm_self)
2365 2 : CALL cp_blacs_env_create(blacs_env=blacs_env, para_env=strict_para_env)
2366 32 : DO i = 1, SIZE(nk)
2367 30 : WRITE (name, '(A,I0)') 'complex_strict_nk_', i
2368 32 : CALL create_one_block(nk(i), TRIM(name), dist_nk, row_n, col_k)
2369 : END DO
2370 52 : DO i = 1, SIZE(kk)
2371 50 : WRITE (name, '(A,I0)') 'complex_strict_kk_', i
2372 52 : CALL create_one_block(kk(i), TRIM(name), dist_kk, row_k, col_k)
2373 : END DO
2374 22 : DO i = 1, SIZE(rot_matrix)
2375 20 : WRITE (name, '(A,I0)') 'complex_strict_rot_', i
2376 22 : CALL create_one_block(rot_matrix(i), TRIM(name), dist_kk, row_k, col_k)
2377 : END DO
2378 :
2379 : overlap(:, :) = CMPLX(0.0_dp, 0.0_dp, KIND=dp)
2380 2 : overlap(1, 1) = CMPLX(1.0_dp, 0.0_dp, KIND=dp)
2381 2 : overlap(2, 2) = CMPLX(1.2_dp, 0.0_dp, KIND=dp)
2382 2 : overlap(3, 3) = CMPLX(1.4_dp, 0.0_dp, KIND=dp)
2383 2 : overlap(1, 2) = CMPLX(0.06_dp, 0.03_dp, KIND=dp)
2384 2 : overlap(2, 1) = CONJG(overlap(1, 2))
2385 2 : overlap(1, 3) = CMPLX(-0.04_dp, 0.02_dp, KIND=dp)
2386 2 : overlap(3, 1) = CONJG(overlap(1, 3))
2387 2 : overlap(2, 3) = CMPLX(0.05_dp, -0.01_dp, KIND=dp)
2388 2 : overlap(3, 2) = CONJG(overlap(2, 3))
2389 :
2390 : c0(:, 1) = [CMPLX(1.0_dp, 0.10_dp, KIND=dp), CMPLX(0.12_dp, -0.08_dp, KIND=dp), &
2391 8 : CMPLX(-0.07_dp, 0.04_dp, KIND=dp)]
2392 : c0(:, 2) = [CMPLX(0.09_dp, 0.05_dp, KIND=dp), CMPLX(0.85_dp, -0.12_dp, KIND=dp), &
2393 8 : CMPLX(0.11_dp, 0.06_dp, KIND=dp)]
2394 90 : CALL dense_inverse_sqrt(MATMUL(CONJG(TRANSPOSE(c0)), MATMUL(overlap, c0)), inverse_sqrt)
2395 38 : c0 = MATMUL(c0, inverse_sqrt)
2396 :
2397 2 : x(:, :) = CMPLX(0.0_dp, 0.0_dp, KIND=dp)
2398 2 : x(3, 1) = CMPLX(0.23_dp, 0.17_dp, KIND=dp)
2399 2 : x(3, 2) = CMPLX(-0.14_dp, 0.21_dp, KIND=dp)
2400 144 : x = x - MATMUL(c0, MATMUL(CONJG(TRANSPOSE(c0)), MATMUL(overlap, x)))
2401 2 : direction(:, :) = CMPLX(0.0_dp, 0.0_dp, KIND=dp)
2402 2 : direction(3, 1) = CMPLX(0.11_dp, -0.09_dp, KIND=dp)
2403 2 : direction(3, 2) = CMPLX(0.07_dp, 0.13_dp, KIND=dp)
2404 144 : direction = direction - MATMUL(c0, MATMUL(CONJG(TRANSPOSE(c0)), MATMUL(overlap, direction)))
2405 :
2406 : hamiltonian(:, :) = CMPLX(0.0_dp, 0.0_dp, KIND=dp)
2407 2 : hamiltonian(1, 1) = CMPLX(-0.70_dp, 0.0_dp, KIND=dp)
2408 2 : hamiltonian(2, 2) = CMPLX(-0.20_dp, 0.0_dp, KIND=dp)
2409 2 : hamiltonian(3, 3) = CMPLX(0.90_dp, 0.0_dp, KIND=dp)
2410 2 : hamiltonian(1, 2) = CMPLX(0.08_dp, 0.04_dp, KIND=dp)
2411 2 : hamiltonian(2, 1) = CONJG(hamiltonian(1, 2))
2412 2 : hamiltonian(1, 3) = CMPLX(-0.12_dp, 0.07_dp, KIND=dp)
2413 2 : hamiltonian(3, 1) = CONJG(hamiltonian(1, 3))
2414 2 : hamiltonian(2, 3) = CMPLX(0.05_dp, -0.09_dp, KIND=dp)
2415 2 : hamiltonian(3, 2) = CONJG(hamiltonian(2, 3))
2416 :
2417 : generator(:, :) = CMPLX(0.0_dp, 0.0_dp, KIND=dp)
2418 2 : generator(1, 1) = CMPLX(0.0_dp, 0.16_dp, KIND=dp)
2419 2 : generator(2, 2) = CMPLX(0.0_dp, -0.09_dp, KIND=dp)
2420 2 : generator(1, 2) = CMPLX(0.18_dp, -0.07_dp, KIND=dp)
2421 2 : generator(2, 1) = -CONJG(generator(1, 2))
2422 : direction_a(:, :) = CMPLX(0.0_dp, 0.0_dp, KIND=dp)
2423 2 : direction_a(1, 1) = CMPLX(0.0_dp, -0.05_dp, KIND=dp)
2424 2 : direction_a(2, 2) = CMPLX(0.0_dp, 0.04_dp, KIND=dp)
2425 2 : direction_a(1, 2) = CMPLX(-0.08_dp, 0.06_dp, KIND=dp)
2426 2 : direction_a(2, 1) = -CONJG(direction_a(1, 2))
2427 2 : occupation = [1.0_dp, 0.34_dp]
2428 :
2429 2 : CALL strict_orbitals(c0, x, overlap, orbitals_base)
2430 2 : rotation = dense_antihermitian_exp(generator)
2431 38 : orbitals = MATMUL(orbitals_base, rotation)
2432 70 : hc = 2.0_dp*MATMUL(hamiltonian, orbitals)
2433 6 : DO i = 1, k
2434 18 : hc(:, i) = occupation(i)*hc(:, i)
2435 : END DO
2436 54 : sx = MATMUL(overlap, x)
2437 54 : sc0 = MATMUL(overlap, c0)
2438 2 : CALL put_complex_pair(nk(1), nk(2), x)
2439 2 : CALL put_complex_pair(nk(3), nk(4), sx)
2440 2 : CALL put_complex_pair(nk(5), nk(6), hc)
2441 2 : CALL put_complex_pair(nk(9), nk(10), c0)
2442 2 : CALL put_complex_pair(nk(11), nk(12), sc0)
2443 2 : CALL put_complex_pair(rot_matrix(3), rot_matrix(4), generator)
2444 :
2445 2 : strict_env%has_complex_kpoint_state = .TRUE.
2446 2 : strict_env%settings%ot_algorithm = 'TOD'
2447 2 : strict_env%settings%do_rotation = .TRUE.
2448 2 : strict_env%para_env => strict_para_env
2449 2 : strict_env%blacs_env => blacs_env
2450 2 : strict_env%matrix_x => nk(1)
2451 2 : strict_env%matrix_x_im => nk(2)
2452 2 : strict_env%matrix_sx => nk(3)
2453 2 : strict_env%matrix_sx_im => nk(4)
2454 2 : strict_env%matrix_gx => nk(7)
2455 2 : strict_env%matrix_gx_im => nk(8)
2456 2 : strict_env%matrix_c0 => nk(9)
2457 2 : strict_env%matrix_c0_im => nk(10)
2458 2 : strict_env%matrix_sc0 => nk(11)
2459 2 : strict_env%matrix_sc0_im => nk(12)
2460 2 : strict_env%matrix_buf_nk => nk(13)
2461 2 : strict_env%matrix_buf_nk_im => nk(14)
2462 2 : strict_env%matrix_tmp_nk => nk(15)
2463 2 : strict_env%matrix_p => kk(1)
2464 2 : strict_env%matrix_p_im => kk(2)
2465 2 : strict_env%matrix_r => kk(3)
2466 2 : strict_env%matrix_r_im => kk(4)
2467 2 : strict_env%matrix_sinp => kk(5)
2468 2 : strict_env%matrix_sinp_im => kk(6)
2469 2 : strict_env%matrix_cosp => kk(7)
2470 2 : strict_env%matrix_cosp_im => kk(8)
2471 2 : strict_env%matrix_sinp_b => kk(9)
2472 2 : strict_env%matrix_cosp_b => kk(10)
2473 2 : strict_env%matrix_buf1 => kk(11)
2474 2 : strict_env%matrix_buf1_im => kk(12)
2475 2 : strict_env%matrix_buf2 => kk(13)
2476 2 : strict_env%matrix_buf2_im => kk(14)
2477 2 : strict_env%matrix_buf3 => kk(15)
2478 2 : strict_env%matrix_buf3_im => kk(16)
2479 2 : strict_env%matrix_buf4 => kk(17)
2480 2 : strict_env%matrix_buf4_im => kk(18)
2481 2 : strict_env%matrix_os => kk(19)
2482 2 : strict_env%matrix_os_im => kk(20)
2483 2 : strict_env%matrix_buf1_ortho => kk(21)
2484 2 : strict_env%matrix_buf1_ortho_im => kk(22)
2485 2 : strict_env%matrix_buf2_ortho => kk(23)
2486 2 : strict_env%matrix_buf2_ortho_im => kk(24)
2487 2 : strict_env%matrix_tmp_ortho => kk(25)
2488 2 : strict_env%rot_mat_u => rot_matrix(1)
2489 2 : strict_env%rot_mat_u_im => rot_matrix(2)
2490 2 : strict_env%rot_mat_x => rot_matrix(3)
2491 2 : strict_env%rot_mat_x_im => rot_matrix(4)
2492 2 : strict_env%rot_mat_dedu => rot_matrix(5)
2493 2 : strict_env%rot_mat_dedu_im => rot_matrix(6)
2494 2 : strict_env%rot_mat_gx => rot_matrix(7)
2495 2 : strict_env%rot_mat_gx_im => rot_matrix(8)
2496 2 : strict_env%rot_mat_evec_re => rot_matrix(9)
2497 2 : strict_env%rot_mat_evec_im => rot_matrix(10)
2498 2 : strict_env%rot_mat_evals => rotation_evals
2499 2 : ALLOCATE (strict_env%evals(k), strict_env%dum(k))
2500 :
2501 2 : x_re_p => nk(1)
2502 2 : x_im_p => nk(2)
2503 2 : sx_re_p => nk(3)
2504 2 : sx_im_p => nk(4)
2505 2 : hc_re_p => nk(5)
2506 2 : hc_im_p => nk(6)
2507 2 : CALL qs_ot_get_p_complex(x_re_p, x_im_p, sx_re_p, sx_im_p, strict_env)
2508 2 : CALL qs_ot_get_derivative_complex(hc_re_p, hc_im_p, strict_env)
2509 :
2510 2 : gradient(:, :) = CMPLX(0.0_dp, 0.0_dp, KIND=dp)
2511 2 : CALL dbcsr_get_block_p(nk(7), 1, 1, block, found)
2512 18 : IF (found) gradient = gradient + CMPLX(block, 0.0_dp, KIND=dp)
2513 2 : CALL dbcsr_get_block_p(nk(8), 1, 1, block, found)
2514 18 : IF (found) gradient = gradient + CMPLX(0.0_dp, block, KIND=dp)
2515 2 : CALL get_complex_pair(rot_matrix(7), rot_matrix(8), gradient_a)
2516 :
2517 2 : hstep = 1.0E-6_dp
2518 : fd_slope = (strict_rotated_energy(c0, x + hstep*direction, generator, overlap, &
2519 : hamiltonian, occupation) - &
2520 : strict_rotated_energy(c0, x - hstep*direction, generator, overlap, &
2521 34 : hamiltonian, occupation))/(2.0_dp*hstep)
2522 18 : slope = REAL(SUM(CONJG(gradient)*direction), KIND=dp)
2523 2 : error = ABS(slope - fd_slope)
2524 : fd_slope_rot = (strict_rotated_energy(c0, x, generator + hstep*direction_a, overlap, &
2525 : hamiltonian, occupation) - &
2526 : strict_rotated_energy(c0, x, generator - hstep*direction_a, overlap, &
2527 26 : hamiltonian, occupation))/(2.0_dp*hstep)
2528 : slope_rot = 0.5_dp*(SUM(REAL(gradient_a, KIND=dp)*REAL(direction_a, KIND=dp)) + &
2529 26 : SUM(AIMAG(gradient_a)*AIMAG(direction_a)))
2530 2 : error_rot = ABS(slope_rot - fd_slope_rot)
2531 : ortho_error = MAXVAL(ABS(MATMUL(CONJG(TRANSPOSE(orbitals)), MATMUL(overlap, orbitals)) - &
2532 : RESHAPE([CMPLX(1.0_dp, 0.0_dp, KIND=dp), &
2533 : CMPLX(0.0_dp, 0.0_dp, KIND=dp), &
2534 : CMPLX(0.0_dp, 0.0_dp, KIND=dp), &
2535 102 : CMPLX(1.0_dp, 0.0_dp, KIND=dp)], [k, k])))
2536 2 : IF (error > 5.0E-8_dp .OR. error_rot > 5.0E-8_dp .OR. &
2537 0 : ortho_error > 5.0E-13_dp) nfail = nfail + 1
2538 2 : IF (para_env%is_source()) THEN
2539 1 : WRITE (io_unit, '(A,3(1X,ES13.6))') 'complex finite-STRICT derivative: fd/predicted/error', &
2540 2 : fd_slope, slope, error
2541 1 : WRITE (io_unit, '(A,3(1X,ES13.6))') 'complex finite-STRICT rotation: fd/predicted/error', &
2542 2 : fd_slope_rot, slope_rot, error_rot
2543 1 : WRITE (io_unit, '(A,1X,ES13.6)') 'complex finite-STRICT orthonormality error', ortho_error
2544 : END IF
2545 :
2546 2 : DEALLOCATE (strict_env%evals, strict_env%dum)
2547 22 : DO i = SIZE(rot_matrix), 1, -1
2548 22 : CALL dbcsr_release(rot_matrix(i))
2549 : END DO
2550 52 : DO i = SIZE(kk), 1, -1
2551 52 : CALL dbcsr_release(kk(i))
2552 : END DO
2553 32 : DO i = SIZE(nk), 1, -1
2554 32 : CALL dbcsr_release(nk(i))
2555 : END DO
2556 2 : CALL cp_blacs_env_release(blacs_env)
2557 2 : CALL mp_para_env_release(strict_para_env)
2558 :
2559 8 : END SUBROUTINE test_complex_strict_derivative
2560 : ! **************************************************************************************************
2561 : !> \brief Check gauge invariance, Hermiticity, and positivity of the complex inverse operator.
2562 : !> \param para_env parallel environment
2563 : !> \param nfail accumulated number of failures
2564 : ! **************************************************************************************************
2565 2 : SUBROUTINE test_complex_preconditioner_gauge(para_env, nfail)
2566 : TYPE(mp_para_env_type), POINTER :: para_env
2567 : INTEGER, INTENT(INOUT) :: nfail
2568 :
2569 : COMPLEX(KIND=dp) :: value
2570 : COMPLEX(KIND=dp), DIMENSION(k) :: phases
2571 : COMPLEX(KIND=dp), DIMENSION(n, k) :: c0, c0_phase, c0_rot, gradient, &
2572 : gradient_phase, gradient_rot, output, &
2573 : output_phase, output_rot
2574 : COMPLEX(KIND=dp), DIMENSION(n, n) :: eigenvectors, full_all_vectors, h, inverse, &
2575 : inverse_rot, inverse_s, inverse_single, inverse_t, reference_s, reference_t, s, &
2576 : spectral_metric, t
2577 : INTEGER :: i, j
2578 : REAL(KIND=dp) :: covariant_descent, covariant_gauge_error, full_all_descent, &
2579 : full_all_gauge_error, full_all_orthogonality_error, full_single_hermitian_error, &
2580 : gauge_error, hermitian_error, low_rank_descent, low_rank_gauge_error, s_error, t_error
2581 : REAL(KIND=dp), DIMENSION(k) :: c0_evals
2582 : REAL(KIND=dp), DIMENSION(n) :: eigenvalues, eigenvalues_single
2583 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
2584 : TYPE(cp_cfm_type) :: matrix_c0, matrix_c0_phase, matrix_c0_rot, matrix_gradient, &
2585 : matrix_gradient_phase, matrix_gradient_rot, matrix_h, matrix_output, matrix_output_phase, &
2586 : matrix_output_rot, matrix_s, matrix_t
2587 : TYPE(cp_fm_struct_type), POINTER :: fmstruct_nk, fmstruct_nn
2588 32 : TYPE(preconditioner_type) :: covariant, covariant_rot, full_all, full_all_phase, low_rank, &
2589 28 : low_rank_rot, preconditioner, preconditioner_rot, preconditioner_s, &
2590 14 : preconditioner_single, preconditioner_t
2591 :
2592 2 : NULLIFY (blacs_env, fmstruct_nk, fmstruct_nn)
2593 2 : CALL cp_blacs_env_create(blacs_env=blacs_env, para_env=para_env)
2594 : CALL cp_fm_struct_create(fmstruct=fmstruct_nn, para_env=para_env, context=blacs_env, &
2595 2 : nrow_global=n, ncol_global=n)
2596 : CALL cp_fm_struct_create(fmstruct=fmstruct_nk, para_env=para_env, context=blacs_env, &
2597 2 : nrow_global=n, ncol_global=k)
2598 2 : CALL cp_cfm_create(matrix_h, fmstruct_nn, name='complex preconditioner H')
2599 2 : CALL cp_cfm_create(matrix_s, fmstruct_nn, name='complex preconditioner S')
2600 2 : CALL cp_cfm_create(matrix_t, fmstruct_nn, name='complex preconditioner T')
2601 2 : CALL cp_cfm_create(matrix_c0, fmstruct_nk, name='complex preconditioner C')
2602 2 : CALL cp_cfm_create(matrix_c0_phase, fmstruct_nk, name='complex preconditioner phased C')
2603 2 : CALL cp_cfm_create(matrix_c0_rot, fmstruct_nk, name='complex preconditioner C U')
2604 2 : CALL cp_cfm_create(matrix_gradient, fmstruct_nk, name='complex preconditioner gradient')
2605 : CALL cp_cfm_create(matrix_gradient_phase, fmstruct_nk, &
2606 2 : name='complex preconditioner phased gradient')
2607 : CALL cp_cfm_create(matrix_gradient_rot, fmstruct_nk, &
2608 2 : name='complex preconditioner rotated gradient')
2609 2 : CALL cp_cfm_create(matrix_output, fmstruct_nk, name='complex preconditioner output')
2610 : CALL cp_cfm_create(matrix_output_phase, fmstruct_nk, &
2611 2 : name='complex preconditioner phased output')
2612 : CALL cp_cfm_create(matrix_output_rot, fmstruct_nk, &
2613 2 : name='complex preconditioner rotated output')
2614 2 : CALL cp_cfm_set_all(matrix_h, CMPLX(0.0_dp, 0.0_dp, KIND=dp))
2615 2 : CALL cp_cfm_set_all(matrix_s, CMPLX(0.0_dp, 0.0_dp, KIND=dp))
2616 2 : CALL cp_cfm_set_all(matrix_t, CMPLX(0.0_dp, 0.0_dp, KIND=dp))
2617 2 : CALL cp_cfm_set_all(matrix_c0, CMPLX(0.0_dp, 0.0_dp, KIND=dp))
2618 2 : CALL cp_cfm_set_all(matrix_c0_phase, CMPLX(0.0_dp, 0.0_dp, KIND=dp))
2619 2 : CALL cp_cfm_set_all(matrix_c0_rot, CMPLX(0.0_dp, 0.0_dp, KIND=dp))
2620 :
2621 : h(:, :) = CMPLX(0.0_dp, 0.0_dp, KIND=dp)
2622 2 : h(1, 1) = CMPLX(-0.35_dp, 0.0_dp, KIND=dp)
2623 2 : h(2, 2) = CMPLX(0.15_dp, 0.0_dp, KIND=dp)
2624 2 : h(3, 3) = CMPLX(0.90_dp, 0.0_dp, KIND=dp)
2625 2 : h(1, 2) = CMPLX(0.08_dp, 0.04_dp, KIND=dp)
2626 2 : h(2, 1) = CONJG(h(1, 2))
2627 2 : h(1, 3) = CMPLX(-0.06_dp, 0.03_dp, KIND=dp)
2628 2 : h(3, 1) = CONJG(h(1, 3))
2629 2 : h(2, 3) = CMPLX(0.05_dp, -0.07_dp, KIND=dp)
2630 2 : h(3, 2) = CONJG(h(2, 3))
2631 2 : s(:, :) = CMPLX(0.0_dp, 0.0_dp, KIND=dp)
2632 2 : s(1, 1) = CMPLX(1.0_dp, 0.0_dp, KIND=dp)
2633 2 : s(2, 2) = CMPLX(1.0_dp, 0.0_dp, KIND=dp)
2634 2 : s(3, 3) = CMPLX(1.4_dp, 0.0_dp, KIND=dp)
2635 2 : s(1, 3) = CMPLX(0.08_dp, 0.03_dp, KIND=dp)
2636 2 : s(3, 1) = CONJG(s(1, 3))
2637 2 : s(2, 3) = CMPLX(-0.04_dp, 0.02_dp, KIND=dp)
2638 2 : s(3, 2) = CONJG(s(2, 3))
2639 2 : t = h
2640 8 : DO i = 1, n
2641 8 : t(i, i) = t(i, i) + CMPLX(1.50_dp, 0.0_dp, KIND=dp)
2642 : END DO
2643 2 : c0(:, :) = CMPLX(0.0_dp, 0.0_dp, KIND=dp)
2644 2 : c0(1, 1) = CMPLX(1.0_dp, 0.0_dp, KIND=dp)
2645 2 : c0(2, 2) = CMPLX(1.0_dp, 0.0_dp, KIND=dp)
2646 2 : c0_rot(:, :) = CMPLX(0.0_dp, 0.0_dp, KIND=dp)
2647 2 : c0_rot(1, 1) = CMPLX(SQRT(0.5_dp), 0.0_dp, KIND=dp)
2648 2 : c0_rot(1, 2) = CMPLX(0.0_dp, SQRT(0.5_dp), KIND=dp)
2649 2 : c0_rot(2, 1) = CMPLX(0.0_dp, SQRT(0.5_dp), KIND=dp)
2650 2 : c0_rot(2, 2) = CMPLX(SQRT(0.5_dp), 0.0_dp, KIND=dp)
2651 : phases = [EXP(CMPLX(0.0_dp, 0.43_dp, KIND=dp)), &
2652 2 : EXP(CMPLX(0.0_dp, -0.67_dp, KIND=dp))]
2653 6 : DO j = 1, k
2654 18 : c0_phase(:, j) = phases(j)*c0(:, j)
2655 : END DO
2656 : gradient(:, 1) = [CMPLX(0.17_dp, 0.02_dp, KIND=dp), &
2657 : CMPLX(-0.08_dp, 0.05_dp, KIND=dp), &
2658 8 : CMPLX(0.11_dp, -0.04_dp, KIND=dp)]
2659 : gradient(:, 2) = [CMPLX(-0.03_dp, 0.06_dp, KIND=dp), &
2660 : CMPLX(0.09_dp, -0.02_dp, KIND=dp), &
2661 8 : CMPLX(-0.14_dp, 0.05_dp, KIND=dp)]
2662 6 : DO j = 1, k
2663 18 : gradient_phase(:, j) = phases(j)*gradient(:, j)
2664 : END DO
2665 38 : gradient_rot = MATMUL(gradient, c0_rot(1:k, :))
2666 6 : c0_evals = [REAL(h(1, 1), KIND=dp), REAL(h(2, 2), KIND=dp)]
2667 :
2668 8 : DO j = 1, n
2669 26 : DO i = 1, n
2670 18 : CALL cp_cfm_set_element(matrix_h, i, j, h(i, j))
2671 18 : CALL cp_cfm_set_element(matrix_s, i, j, s(i, j))
2672 24 : CALL cp_cfm_set_element(matrix_t, i, j, t(i, j))
2673 : END DO
2674 : END DO
2675 6 : DO j = 1, k
2676 18 : DO i = 1, n
2677 12 : CALL cp_cfm_set_element(matrix_c0, i, j, c0(i, j))
2678 12 : CALL cp_cfm_set_element(matrix_c0_phase, i, j, c0_phase(i, j))
2679 12 : CALL cp_cfm_set_element(matrix_c0_rot, i, j, c0_rot(i, j))
2680 12 : CALL cp_cfm_set_element(matrix_gradient, i, j, gradient(i, j))
2681 12 : CALL cp_cfm_set_element(matrix_gradient_phase, i, j, gradient_phase(i, j))
2682 16 : CALL cp_cfm_set_element(matrix_gradient_rot, i, j, gradient_rot(i, j))
2683 : END DO
2684 : END DO
2685 :
2686 2 : CALL init_preconditioner(full_all, para_env, blacs_env)
2687 2 : CALL init_preconditioner(full_all_phase, para_env, blacs_env)
2688 2 : CALL init_preconditioner(low_rank, para_env, blacs_env)
2689 2 : CALL init_preconditioner(low_rank_rot, para_env, blacs_env)
2690 2 : CALL init_preconditioner(covariant, para_env, blacs_env)
2691 2 : CALL init_preconditioner(covariant_rot, para_env, blacs_env)
2692 2 : CALL init_preconditioner(preconditioner, para_env, blacs_env)
2693 2 : CALL init_preconditioner(preconditioner_rot, para_env, blacs_env)
2694 2 : CALL init_preconditioner(preconditioner_s, para_env, blacs_env)
2695 2 : CALL init_preconditioner(preconditioner_single, para_env, blacs_env)
2696 2 : CALL init_preconditioner(preconditioner_t, para_env, blacs_env)
2697 2 : CALL make_complex_full_single_inverse(preconditioner, matrix_c0, matrix_h, matrix_s, 0.20_dp)
2698 2 : CALL make_complex_full_single_inverse(preconditioner_rot, matrix_c0_rot, matrix_h, matrix_s, 0.20_dp)
2699 2 : CALL make_complex_full_s_inverse(preconditioner_s, matrix_s)
2700 2 : CALL make_complex_full_single(preconditioner_single, matrix_h, matrix_s, 0.15_dp, 0.20_dp)
2701 2 : CALL make_complex_full_kinetic(preconditioner_t, matrix_t, matrix_s, 0.20_dp)
2702 2 : CALL make_complex_full_all(full_all, matrix_c0, matrix_h, matrix_s, c0_evals, 0.08_dp)
2703 : CALL make_complex_full_all(full_all_phase, matrix_c0_phase, matrix_h, matrix_s, &
2704 2 : c0_evals, 0.08_dp)
2705 2 : CALL apply_preconditioner(full_all, matrix_gradient, matrix_output)
2706 2 : CALL apply_preconditioner(full_all_phase, matrix_gradient_phase, matrix_output_phase)
2707 8 : DO j = 1, n
2708 26 : DO i = 1, n
2709 18 : CALL cp_cfm_get_element(preconditioner%complex_fm, i, j, value)
2710 18 : inverse(i, j) = value
2711 18 : CALL cp_cfm_get_element(preconditioner_rot%complex_fm, i, j, value)
2712 18 : inverse_rot(i, j) = value
2713 18 : CALL cp_cfm_get_element(preconditioner_s%complex_fm, i, j, value)
2714 18 : inverse_s(i, j) = value
2715 18 : CALL cp_cfm_get_element(preconditioner_single%complex_fm, i, j, value)
2716 18 : inverse_single(i, j) = value
2717 18 : CALL cp_cfm_get_element(preconditioner_t%complex_fm, i, j, value)
2718 24 : inverse_t(i, j) = value
2719 : END DO
2720 : END DO
2721 26 : gauge_error = MAXVAL(ABS(inverse - inverse_rot))
2722 26 : hermitian_error = MAXVAL(ABS(inverse - CONJG(TRANSPOSE(inverse))))
2723 2 : CALL diag_complex(inverse, eigenvectors, eigenvalues)
2724 26 : full_single_hermitian_error = MAXVAL(ABS(inverse_single - CONJG(TRANSPOSE(inverse_single))))
2725 2 : CALL diag_complex(inverse_single, eigenvectors, eigenvalues_single)
2726 2 : CALL dense_hermitian_inverse(s, reference_s)
2727 26 : CALL dense_hermitian_inverse(t + 0.20_dp*s, reference_t)
2728 26 : s_error = MAXVAL(ABS(inverse_s - reference_s))
2729 26 : t_error = MAXVAL(ABS(inverse_t - reference_t))
2730 6 : DO j = 1, k
2731 18 : DO i = 1, n
2732 12 : CALL cp_cfm_get_element(matrix_output, i, j, value)
2733 12 : output(i, j) = value
2734 12 : CALL cp_cfm_get_element(matrix_output_phase, i, j, value)
2735 16 : output_phase(i, j) = value
2736 : END DO
2737 : END DO
2738 8 : DO j = 1, n
2739 26 : DO i = 1, n
2740 18 : CALL cp_cfm_get_element(full_all%complex_fm, i, j, value)
2741 24 : full_all_vectors(i, j) = value
2742 : END DO
2743 : END DO
2744 2 : full_all_gauge_error = 0.0_dp
2745 6 : DO j = 1, k
2746 : full_all_gauge_error = MAX(full_all_gauge_error, &
2747 18 : MAXVAL(ABS(output_phase(:, j) - phases(j)*output(:, j))))
2748 : END DO
2749 : spectral_metric = MATMUL(CONJG(TRANSPOSE(full_all_vectors)), &
2750 158 : MATMUL(s, full_all_vectors))
2751 : full_all_orthogonality_error = MAXVAL(ABS(spectral_metric - &
2752 : RESHAPE([CMPLX(1.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(0.0_dp, 0.0_dp, KIND=dp), &
2756 : CMPLX(1.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(0.0_dp, 0.0_dp, KIND=dp), &
2760 26 : CMPLX(1.0_dp, 0.0_dp, KIND=dp)], [n, n])))
2761 18 : full_all_descent = REAL(SUM(CONJG(gradient)*output), KIND=dp)
2762 :
2763 : CALL make_complex_fermi_low_rank(low_rank, matrix_c0, matrix_h, matrix_s, &
2764 2 : 0.08_dp, 1, 1.0_dp, 1.0E-8_dp)
2765 : CALL make_complex_fermi_low_rank(low_rank_rot, matrix_c0_rot, matrix_h, matrix_s, &
2766 2 : 0.08_dp, 1, 1.0_dp, 1.0E-8_dp)
2767 2 : CALL apply_preconditioner(low_rank, matrix_gradient, matrix_output)
2768 2 : CALL apply_preconditioner(low_rank_rot, matrix_gradient_rot, matrix_output_rot)
2769 6 : DO j = 1, k
2770 18 : DO i = 1, n
2771 12 : CALL cp_cfm_get_element(matrix_output, i, j, value)
2772 12 : output(i, j) = value
2773 12 : CALL cp_cfm_get_element(matrix_output_rot, i, j, value)
2774 16 : output_rot(i, j) = value
2775 : END DO
2776 : END DO
2777 54 : low_rank_gauge_error = MAXVAL(ABS(output_rot - MATMUL(output, c0_rot(1:k, :))))
2778 18 : low_rank_descent = REAL(SUM(CONJG(gradient)*output), KIND=dp)
2779 :
2780 2 : CALL make_complex_full_all_covariant(covariant, matrix_c0, matrix_h, matrix_s, 0.08_dp)
2781 2 : CALL make_complex_full_all_covariant(covariant_rot, matrix_c0_rot, matrix_h, matrix_s, 0.08_dp)
2782 2 : CALL apply_preconditioner(covariant, matrix_gradient, matrix_output)
2783 2 : CALL apply_preconditioner(covariant_rot, matrix_gradient_rot, matrix_output_rot)
2784 6 : DO j = 1, k
2785 18 : DO i = 1, n
2786 12 : CALL cp_cfm_get_element(matrix_output, i, j, value)
2787 12 : output(i, j) = value
2788 12 : CALL cp_cfm_get_element(matrix_output_rot, i, j, value)
2789 16 : output_rot(i, j) = value
2790 : END DO
2791 : END DO
2792 54 : covariant_gauge_error = MAXVAL(ABS(output_rot - MATMUL(output, c0_rot(1:k, :))))
2793 18 : covariant_descent = REAL(SUM(CONJG(gradient)*output), KIND=dp)
2794 8 : IF (gauge_error > 1.0E-11_dp .OR. hermitian_error > 1.0E-12_dp .OR. &
2795 0 : MINVAL(eigenvalues) <= 0.0_dp) nfail = nfail + 1
2796 2 : IF (s_error > 1.0E-11_dp .OR. t_error > 1.0E-11_dp) nfail = nfail + 1
2797 8 : IF (full_single_hermitian_error > 1.0E-12_dp .OR. &
2798 0 : MINVAL(eigenvalues_single) <= 0.0_dp) nfail = nfail + 1
2799 : IF (full_all_gauge_error > 1.0E-11_dp .OR. &
2800 2 : full_all_orthogonality_error > 1.0E-11_dp .OR. full_all_descent <= 0.0_dp) nfail = nfail + 1
2801 2 : IF (low_rank_gauge_error > 1.0E-11_dp .OR. low_rank_descent <= 0.0_dp) nfail = nfail + 1
2802 2 : IF (covariant_gauge_error > 1.0E-11_dp .OR. covariant_descent <= 0.0_dp) nfail = nfail + 1
2803 2 : IF (para_env%is_source()) THEN
2804 1 : WRITE (io_unit, '(A,3(1X,ES13.6))') 'complex FULL_SINGLE gauge/hermitian/min-eigenvalue', &
2805 5 : gauge_error, hermitian_error, MINVAL(eigenvalues)
2806 1 : WRITE (io_unit, '(A,2(1X,ES13.6))') 'complex S/KINETIC inverse errors', s_error, t_error
2807 1 : WRITE (io_unit, '(A,2(1X,ES13.6))') 'complex FULL_SINGLE hermitian/min-eigenvalue', &
2808 5 : full_single_hermitian_error, MINVAL(eigenvalues_single)
2809 1 : WRITE (io_unit, '(A,3(1X,ES13.6))') 'complex FULL_ALL phase/orthogonality/descent', &
2810 2 : full_all_gauge_error, full_all_orthogonality_error, full_all_descent
2811 1 : WRITE (io_unit, '(A,2(1X,ES13.6))') 'complex FERMI_LOW_RANK covariance/descent', &
2812 2 : low_rank_gauge_error, low_rank_descent
2813 1 : WRITE (io_unit, '(A,2(1X,ES13.6))') 'complex FULL_ALL_COVARIANT covariance/descent', &
2814 2 : covariant_gauge_error, covariant_descent
2815 : END IF
2816 :
2817 2 : CALL destroy_preconditioner(covariant_rot)
2818 2 : CALL destroy_preconditioner(covariant)
2819 2 : CALL destroy_preconditioner(low_rank_rot)
2820 2 : CALL destroy_preconditioner(low_rank)
2821 2 : CALL destroy_preconditioner(full_all_phase)
2822 2 : CALL destroy_preconditioner(full_all)
2823 2 : CALL destroy_preconditioner(preconditioner_t)
2824 2 : CALL destroy_preconditioner(preconditioner_single)
2825 2 : CALL destroy_preconditioner(preconditioner_s)
2826 2 : CALL destroy_preconditioner(preconditioner_rot)
2827 2 : CALL destroy_preconditioner(preconditioner)
2828 2 : CALL cp_cfm_release(matrix_c0_rot)
2829 2 : CALL cp_cfm_release(matrix_c0_phase)
2830 2 : CALL cp_cfm_release(matrix_c0)
2831 2 : CALL cp_cfm_release(matrix_output_phase)
2832 2 : CALL cp_cfm_release(matrix_output)
2833 2 : CALL cp_cfm_release(matrix_output_rot)
2834 2 : CALL cp_cfm_release(matrix_gradient_rot)
2835 2 : CALL cp_cfm_release(matrix_gradient_phase)
2836 2 : CALL cp_cfm_release(matrix_gradient)
2837 2 : CALL cp_cfm_release(matrix_t)
2838 2 : CALL cp_cfm_release(matrix_s)
2839 2 : CALL cp_cfm_release(matrix_h)
2840 2 : CALL cp_fm_struct_release(fmstruct_nk)
2841 2 : CALL cp_fm_struct_release(fmstruct_nn)
2842 2 : CALL cp_blacs_env_release(blacs_env)
2843 :
2844 6 : END SUBROUTINE test_complex_preconditioner_gauge
2845 :
2846 : ! **************************************************************************************************
2847 : !> \brief Check complex covariant preconditioning within equal-occupation ROKS blocks.
2848 : !> \param para_env parallel environment
2849 : !> \param nfail accumulated number of failures
2850 : ! **************************************************************************************************
2851 2 : SUBROUTINE test_complex_occupation_block_preconditioner(para_env, nfail)
2852 : TYPE(mp_para_env_type), POINTER :: para_env
2853 : INTEGER, INTENT(INOUT) :: nfail
2854 :
2855 : INTEGER, PARAMETER :: nbas = 4, norb = 3
2856 :
2857 : COMPLEX(KIND=dp) :: value
2858 : COMPLEX(KIND=dp), DIMENSION(nbas, nbas) :: hamiltonian, overlap
2859 : COMPLEX(KIND=dp), DIMENSION(nbas, norb) :: coeff, coeff_rot, gradient, &
2860 : gradient_rot, output, output_rot
2861 : COMPLEX(KIND=dp), DIMENSION(norb, norb) :: block_rotation, private_rotation
2862 : INTEGER :: i, j
2863 : REAL(KIND=dp) :: covariance_error, cross_block_error, &
2864 : descent, shifted_eigenvalue_error
2865 : REAL(KIND=dp), DIMENSION(norb, 2) :: occupation_signature
2866 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
2867 : TYPE(cp_cfm_type) :: matrix_c, matrix_c_rot, matrix_gradient, matrix_gradient_rot, matrix_h, &
2868 : matrix_output, matrix_output_rot, matrix_s
2869 : TYPE(cp_fm_struct_type), POINTER :: fmstruct_nk, fmstruct_nn
2870 14 : TYPE(preconditioner_type) :: preconditioner, preconditioner_rot
2871 :
2872 2 : NULLIFY (blacs_env, fmstruct_nk, fmstruct_nn)
2873 2 : CALL cp_blacs_env_create(blacs_env=blacs_env, para_env=para_env)
2874 : CALL cp_fm_struct_create(fmstruct=fmstruct_nn, para_env=para_env, context=blacs_env, &
2875 2 : nrow_global=nbas, ncol_global=nbas)
2876 : CALL cp_fm_struct_create(fmstruct=fmstruct_nk, para_env=para_env, context=blacs_env, &
2877 2 : nrow_global=nbas, ncol_global=norb)
2878 2 : CALL cp_cfm_create(matrix_h, fmstruct_nn, name='occupation-block H')
2879 2 : CALL cp_cfm_create(matrix_s, fmstruct_nn, name='occupation-block S')
2880 2 : CALL cp_cfm_create(matrix_c, fmstruct_nk, name='occupation-block C')
2881 2 : CALL cp_cfm_create(matrix_c_rot, fmstruct_nk, name='occupation-block C U')
2882 2 : CALL cp_cfm_create(matrix_gradient, fmstruct_nk, name='occupation-block gradient')
2883 : CALL cp_cfm_create(matrix_gradient_rot, fmstruct_nk, &
2884 2 : name='occupation-block rotated gradient')
2885 2 : CALL cp_cfm_create(matrix_output, fmstruct_nk, name='occupation-block output')
2886 2 : CALL cp_cfm_create(matrix_output_rot, fmstruct_nk, name='occupation-block rotated output')
2887 :
2888 : hamiltonian(:, :) = CMPLX(0.0_dp, 0.0_dp, KIND=dp)
2889 2 : hamiltonian(1, 1) = CMPLX(-0.42_dp, 0.0_dp, KIND=dp)
2890 2 : hamiltonian(2, 2) = CMPLX(-0.18_dp, 0.0_dp, KIND=dp)
2891 2 : hamiltonian(3, 3) = CMPLX(0.07_dp, 0.0_dp, KIND=dp)
2892 2 : hamiltonian(4, 4) = CMPLX(0.83_dp, 0.0_dp, KIND=dp)
2893 2 : hamiltonian(1, 2) = CMPLX(0.09_dp, 0.04_dp, KIND=dp)
2894 2 : hamiltonian(2, 1) = CONJG(hamiltonian(1, 2))
2895 2 : hamiltonian(1, 3) = CMPLX(-0.06_dp, 0.03_dp, KIND=dp)
2896 2 : hamiltonian(3, 1) = CONJG(hamiltonian(1, 3))
2897 2 : hamiltonian(2, 3) = CMPLX(0.05_dp, -0.02_dp, KIND=dp)
2898 2 : hamiltonian(3, 2) = CONJG(hamiltonian(2, 3))
2899 2 : hamiltonian(1, 4) = CMPLX(0.04_dp, -0.03_dp, KIND=dp)
2900 2 : hamiltonian(4, 1) = CONJG(hamiltonian(1, 4))
2901 2 : hamiltonian(2, 4) = CMPLX(-0.02_dp, 0.05_dp, KIND=dp)
2902 2 : hamiltonian(4, 2) = CONJG(hamiltonian(2, 4))
2903 2 : hamiltonian(3, 4) = CMPLX(0.07_dp, 0.01_dp, KIND=dp)
2904 2 : hamiltonian(4, 3) = CONJG(hamiltonian(3, 4))
2905 2 : overlap(:, :) = CMPLX(0.0_dp, 0.0_dp, KIND=dp)
2906 10 : DO i = 1, nbas
2907 10 : overlap(i, i) = CMPLX(1.0_dp, 0.0_dp, KIND=dp)
2908 : END DO
2909 2 : coeff(:, :) = CMPLX(0.0_dp, 0.0_dp, KIND=dp)
2910 8 : DO i = 1, norb
2911 8 : coeff(i, i) = CMPLX(1.0_dp, 0.0_dp, KIND=dp)
2912 : END DO
2913 2 : block_rotation(:, :) = CMPLX(0.0_dp, 0.0_dp, KIND=dp)
2914 2 : block_rotation(1, 1) = CMPLX(SQRT(0.7_dp), 0.0_dp, KIND=dp)
2915 2 : block_rotation(1, 3) = CMPLX(0.0_dp, SQRT(0.3_dp), KIND=dp)
2916 2 : block_rotation(2, 2) = EXP(CMPLX(0.0_dp, 0.37_dp, KIND=dp))
2917 2 : block_rotation(3, 1) = CMPLX(0.0_dp, SQRT(0.3_dp), KIND=dp)
2918 2 : block_rotation(3, 3) = CMPLX(SQRT(0.7_dp), 0.0_dp, KIND=dp)
2919 98 : coeff_rot = MATMUL(coeff, block_rotation)
2920 : gradient(:, 1) = [CMPLX(0.12_dp, 0.03_dp, KIND=dp), &
2921 : CMPLX(-0.04_dp, 0.07_dp, KIND=dp), &
2922 : CMPLX(0.08_dp, -0.02_dp, KIND=dp), &
2923 10 : CMPLX(-0.05_dp, 0.06_dp, KIND=dp)]
2924 : gradient(:, 2) = [CMPLX(-0.06_dp, 0.02_dp, KIND=dp), &
2925 : CMPLX(0.11_dp, -0.05_dp, KIND=dp), &
2926 : CMPLX(0.03_dp, 0.09_dp, KIND=dp), &
2927 10 : CMPLX(0.07_dp, -0.01_dp, KIND=dp)]
2928 : gradient(:, 3) = [CMPLX(0.04_dp, -0.08_dp, KIND=dp), &
2929 : CMPLX(-0.02_dp, 0.05_dp, KIND=dp), &
2930 : CMPLX(0.10_dp, 0.01_dp, KIND=dp), &
2931 10 : CMPLX(-0.03_dp, -0.06_dp, KIND=dp)]
2932 98 : gradient_rot = MATMUL(gradient, block_rotation)
2933 8 : occupation_signature(:, 1) = [1.0_dp, 1.0_dp, 1.0_dp]
2934 8 : occupation_signature(:, 2) = [1.0_dp, 0.0_dp, 1.0_dp]
2935 :
2936 10 : DO j = 1, nbas
2937 42 : DO i = 1, nbas
2938 32 : CALL cp_cfm_set_element(matrix_h, i, j, hamiltonian(i, j))
2939 40 : CALL cp_cfm_set_element(matrix_s, i, j, overlap(i, j))
2940 : END DO
2941 : END DO
2942 8 : DO j = 1, norb
2943 32 : DO i = 1, nbas
2944 24 : CALL cp_cfm_set_element(matrix_c, i, j, coeff(i, j))
2945 24 : CALL cp_cfm_set_element(matrix_c_rot, i, j, coeff_rot(i, j))
2946 24 : CALL cp_cfm_set_element(matrix_gradient, i, j, gradient(i, j))
2947 30 : CALL cp_cfm_set_element(matrix_gradient_rot, i, j, gradient_rot(i, j))
2948 : END DO
2949 : END DO
2950 :
2951 2 : CALL init_preconditioner(preconditioner, para_env, blacs_env)
2952 2 : CALL init_preconditioner(preconditioner_rot, para_env, blacs_env)
2953 : CALL make_complex_full_all_covariant(preconditioner, matrix_c, matrix_h, matrix_s, &
2954 2 : 0.08_dp, occupation_signature)
2955 : CALL make_complex_full_all_covariant(preconditioner_rot, matrix_c_rot, matrix_h, matrix_s, &
2956 2 : 0.08_dp, occupation_signature)
2957 2 : CALL apply_preconditioner(preconditioner, matrix_gradient, matrix_output)
2958 2 : CALL apply_preconditioner(preconditioner_rot, matrix_gradient_rot, matrix_output_rot)
2959 8 : DO j = 1, norb
2960 32 : DO i = 1, nbas
2961 24 : CALL cp_cfm_get_element(matrix_output, i, j, value)
2962 24 : output(i, j) = value
2963 24 : CALL cp_cfm_get_element(matrix_output_rot, i, j, value)
2964 30 : output_rot(i, j) = value
2965 : END DO
2966 : END DO
2967 8 : DO j = 1, norb
2968 26 : DO i = 1, norb
2969 18 : CALL cp_cfm_get_element(preconditioner%occ_rotation_complex, i, j, value)
2970 24 : private_rotation(i, j) = value
2971 : END DO
2972 : END DO
2973 :
2974 128 : covariance_error = MAXVAL(ABS(output_rot - MATMUL(output, block_rotation)))
2975 : cross_block_error = MAX(MAXVAL(ABS(private_rotation([1, 3], 3))), &
2976 10 : MAXVAL(ABS(private_rotation(2, 1:2))))
2977 8 : shifted_eigenvalue_error = MAXVAL(ABS(preconditioner%occ_evals))
2978 32 : descent = REAL(SUM(CONJG(gradient)*output), KIND=dp)
2979 : IF (covariance_error > 1.0E-11_dp .OR. cross_block_error > 1.0E-12_dp .OR. &
2980 2 : shifted_eigenvalue_error > 1.0_dp .OR. descent <= 0.0_dp) nfail = nfail + 1
2981 2 : IF (para_env%is_source()) THEN
2982 : WRITE (io_unit, '(A,4(1X,ES13.6))') &
2983 1 : 'complex occupation-block covariance/cross-block/max-eigenvalue/descent', &
2984 2 : covariance_error, cross_block_error, shifted_eigenvalue_error, descent
2985 : END IF
2986 :
2987 2 : CALL destroy_preconditioner(preconditioner_rot)
2988 2 : CALL destroy_preconditioner(preconditioner)
2989 2 : CALL cp_cfm_release(matrix_output_rot)
2990 2 : CALL cp_cfm_release(matrix_output)
2991 2 : CALL cp_cfm_release(matrix_gradient_rot)
2992 2 : CALL cp_cfm_release(matrix_gradient)
2993 2 : CALL cp_cfm_release(matrix_c_rot)
2994 2 : CALL cp_cfm_release(matrix_c)
2995 2 : CALL cp_cfm_release(matrix_s)
2996 2 : CALL cp_cfm_release(matrix_h)
2997 2 : CALL cp_fm_struct_release(fmstruct_nk)
2998 2 : CALL cp_fm_struct_release(fmstruct_nn)
2999 2 : CALL cp_blacs_env_release(blacs_env)
3000 :
3001 6 : END SUBROUTINE test_complex_occupation_block_preconditioner
3002 :
3003 : ! **************************************************************************************************
3004 : !> \brief Check the intentionally narrow complex k-point preconditioner support matrix.
3005 : !> \param nfail accumulated number of failures
3006 : ! **************************************************************************************************
3007 2 : SUBROUTINE test_preconditioner_support(nfail)
3008 : INTEGER, INTENT(INOUT) :: nfail
3009 :
3010 2 : IF (.NOT. qs_ot_kpoint_preconditioner_supported(ot_precond_none, .FALSE.)) nfail = nfail + 1
3011 2 : IF (.NOT. qs_ot_kpoint_preconditioner_supported(ot_precond_fermi_low_rank, .FALSE.)) nfail = nfail + 1
3012 2 : IF (.NOT. qs_ot_kpoint_preconditioner_supported(ot_precond_full_single_inverse, .FALSE.)) THEN
3013 0 : nfail = nfail + 1
3014 : END IF
3015 2 : IF (.NOT. qs_ot_kpoint_preconditioner_supported(ot_precond_full_all, .FALSE.)) nfail = nfail + 1
3016 2 : IF (.NOT. qs_ot_kpoint_preconditioner_supported(ot_precond_full_single, .FALSE.)) nfail = nfail + 1
3017 2 : IF (.NOT. qs_ot_kpoint_preconditioner_supported( &
3018 0 : ot_precond_full_all_covariant, .FALSE.)) nfail = nfail + 1
3019 2 : IF (.NOT. qs_ot_kpoint_preconditioner_supported(ot_precond_full_kinetic, .FALSE.)) nfail = nfail + 1
3020 2 : IF (.NOT. qs_ot_kpoint_preconditioner_supported(ot_precond_s_inverse, .FALSE.)) nfail = nfail + 1
3021 2 : IF (qs_ot_kpoint_preconditioner_supported(ot_precond_full_single_inverse, .TRUE.)) nfail = nfail + 1
3022 :
3023 2 : IF (.NOT. qs_ot_kpoint_preconditioner_solver_supported( &
3024 0 : ot_precond_full_single, ot_precond_solver_default, .FALSE.)) nfail = nfail + 1
3025 2 : IF (qs_ot_kpoint_preconditioner_solver_supported( &
3026 0 : ot_precond_full_single, ot_precond_solver_inv_chol, .FALSE.)) nfail = nfail + 1
3027 :
3028 2 : IF (.NOT. qs_ot_kpoint_preconditioner_solver_supported( &
3029 0 : ot_precond_full_single_inverse, ot_precond_solver_default, .FALSE.)) nfail = nfail + 1
3030 2 : IF (.NOT. qs_ot_kpoint_preconditioner_solver_supported( &
3031 0 : ot_precond_full_single_inverse, ot_precond_solver_inv_chol, .FALSE.)) nfail = nfail + 1
3032 2 : IF (qs_ot_kpoint_preconditioner_solver_supported( &
3033 0 : ot_precond_full_single_inverse, ot_precond_solver_direct, .FALSE.)) nfail = nfail + 1
3034 2 : IF (qs_ot_kpoint_preconditioner_solver_supported( &
3035 0 : ot_precond_full_single_inverse, ot_precond_solver_chebyshev, .FALSE.)) nfail = nfail + 1
3036 2 : IF (.NOT. qs_ot_kpoint_preconditioner_solver_supported( &
3037 0 : ot_precond_full_all, ot_precond_solver_default, .FALSE.)) nfail = nfail + 1
3038 2 : IF (qs_ot_kpoint_preconditioner_solver_supported( &
3039 0 : ot_precond_full_all, ot_precond_solver_inv_chol, .FALSE.)) nfail = nfail + 1
3040 2 : IF (.NOT. qs_ot_kpoint_preconditioner_solver_supported( &
3041 0 : ot_precond_fermi_low_rank, ot_precond_solver_default, .FALSE.)) nfail = nfail + 1
3042 2 : IF (qs_ot_kpoint_preconditioner_solver_supported( &
3043 0 : ot_precond_fermi_low_rank, ot_precond_solver_inv_chol, .FALSE.)) nfail = nfail + 1
3044 2 : IF (.NOT. qs_ot_kpoint_preconditioner_solver_supported( &
3045 0 : ot_precond_full_all_covariant, ot_precond_solver_default, .FALSE.)) nfail = nfail + 1
3046 2 : IF (qs_ot_kpoint_preconditioner_solver_supported( &
3047 0 : ot_precond_full_all_covariant, ot_precond_solver_inv_chol, .FALSE.)) nfail = nfail + 1
3048 2 : IF (.NOT. qs_ot_kpoint_preconditioner_solver_supported( &
3049 0 : ot_precond_full_kinetic, ot_precond_solver_default, .FALSE.)) nfail = nfail + 1
3050 2 : IF (.NOT. qs_ot_kpoint_preconditioner_solver_supported( &
3051 0 : ot_precond_full_kinetic, ot_precond_solver_inv_chol, .FALSE.)) nfail = nfail + 1
3052 2 : IF (qs_ot_kpoint_preconditioner_solver_supported( &
3053 0 : ot_precond_full_kinetic, ot_precond_solver_direct, .FALSE.)) nfail = nfail + 1
3054 2 : IF (.NOT. qs_ot_kpoint_preconditioner_solver_supported( &
3055 0 : ot_precond_s_inverse, ot_precond_solver_default, .FALSE.)) nfail = nfail + 1
3056 2 : IF (.NOT. qs_ot_kpoint_preconditioner_solver_supported( &
3057 0 : ot_precond_s_inverse, ot_precond_solver_inv_chol, .FALSE.)) nfail = nfail + 1
3058 2 : IF (qs_ot_kpoint_preconditioner_solver_supported( &
3059 0 : ot_precond_s_inverse, ot_precond_solver_chebyshev, .FALSE.)) nfail = nfail + 1
3060 2 : IF (ABS(qs_ot_kpoint_preconditioner_scale(0.25_dp) - 4.0_dp) > EPSILON(1.0_dp)) THEN
3061 0 : nfail = nfail + 1
3062 : END IF
3063 :
3064 2 : END SUBROUTINE test_preconditioner_support
3065 :
3066 : ! **************************************************************************************************
3067 : !> \brief Create and reserve a one-block DBCSR matrix.
3068 : !> \param matrix ...
3069 : !> \param name ...
3070 : !> \param dist ...
3071 : !> \param row_size ...
3072 : !> \param col_size ...
3073 : ! **************************************************************************************************
3074 210 : SUBROUTINE create_one_block(matrix, name, dist, row_size, col_size)
3075 : TYPE(dbcsr_type), INTENT(INOUT) :: matrix
3076 : CHARACTER(LEN=*), INTENT(IN) :: name
3077 : TYPE(dbcsr_distribution_type), INTENT(IN) :: dist
3078 : INTEGER, DIMENSION(:), INTENT(INOUT), POINTER :: row_size, col_size
3079 :
3080 210 : CALL dbcsr_create(matrix, name, dist, dbcsr_type_no_symmetry, row_size, col_size)
3081 210 : CALL dbcsr_reserve_blocks(matrix, rows=[1], cols=[1])
3082 210 : CALL dbcsr_finalize(matrix)
3083 :
3084 210 : END SUBROUTINE create_one_block
3085 :
3086 : ! **************************************************************************************************
3087 : !> \brief Store a dense complex matrix in a real/imaginary DBCSR pair.
3088 : !> \param matrix_re ...
3089 : !> \param matrix_im ...
3090 : !> \param matrix ...
3091 : ! **************************************************************************************************
3092 28 : SUBROUTINE put_complex_pair(matrix_re, matrix_im, matrix)
3093 : TYPE(dbcsr_type), INTENT(INOUT) :: matrix_re, matrix_im
3094 : COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(IN) :: matrix
3095 :
3096 248 : CALL dbcsr_put_block(matrix_re, 1, 1, REAL(matrix, KIND=dp))
3097 248 : CALL dbcsr_put_block(matrix_im, 1, 1, AIMAG(matrix))
3098 :
3099 28 : END SUBROUTINE put_complex_pair
3100 :
3101 : ! **************************************************************************************************
3102 : !> \brief Retrieve a dense complex matrix from a one-block DBCSR pair.
3103 : !> \param matrix_re ...
3104 : !> \param matrix_im ...
3105 : !> \param matrix ...
3106 : ! **************************************************************************************************
3107 8 : SUBROUTINE get_complex_pair(matrix_re, matrix_im, matrix)
3108 : TYPE(dbcsr_type), INTENT(INOUT) :: matrix_re, matrix_im
3109 : COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: matrix
3110 :
3111 : LOGICAL :: found
3112 8 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: block
3113 :
3114 64 : matrix(:, :) = CMPLX(0.0_dp, 0.0_dp, KIND=dp)
3115 8 : CALL dbcsr_get_block_p(matrix_re, 1, 1, block, found)
3116 64 : IF (found) matrix = matrix + CMPLX(block, 0.0_dp, KIND=dp)
3117 8 : CALL dbcsr_get_block_p(matrix_im, 1, 1, block, found)
3118 64 : IF (found) matrix = matrix + CMPLX(0.0_dp, block, KIND=dp)
3119 :
3120 8 : END SUBROUTINE get_complex_pair
3121 :
3122 : ! **************************************************************************************************
3123 : !> \brief Dense fixed-occupation energy after finite REF orthogonalization and rotation.
3124 : !> \param x REF coordinate
3125 : !> \param generator anti-Hermitian rotation generator
3126 : !> \param hamiltonian Hermitian one-particle Hamiltonian
3127 : !> \param occupation orbital occupations
3128 : !> \return energy
3129 : ! **************************************************************************************************
3130 8 : FUNCTION rotated_polar_energy(x, generator, hamiltonian, occupation) RESULT(energy)
3131 : COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(IN) :: x, generator, hamiltonian
3132 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: occupation
3133 : REAL(KIND=dp) :: energy
3134 :
3135 : COMPLEX(KIND=dp), &
3136 16 : DIMENSION(SIZE(x, 1), SIZE(x, 2)) :: orbitals
3137 : COMPLEX(KIND=dp), &
3138 16 : DIMENSION(SIZE(x, 2), SIZE(x, 2)) :: inverse_sqrt, rotation
3139 : INTEGER :: i
3140 :
3141 176 : CALL dense_inverse_sqrt(MATMUL(CONJG(TRANSPOSE(x)), x), inverse_sqrt)
3142 8 : rotation = dense_antihermitian_exp(generator)
3143 496 : orbitals = MATMUL(MATMUL(x, inverse_sqrt), rotation)
3144 8 : energy = 0.0_dp
3145 24 : DO i = 1, SIZE(occupation)
3146 : energy = energy + occupation(i)* &
3147 296 : REAL(DOT_PRODUCT(orbitals(:, i), MATMUL(hamiltonian, orbitals(:, i))), KIND=dp)
3148 : END DO
3149 :
3150 8 : END FUNCTION rotated_polar_energy
3151 :
3152 : ! **************************************************************************************************
3153 : !> \brief Dense reference exponential for an anti-Hermitian matrix.
3154 : !> \param generator anti-Hermitian generator
3155 : !> \return its unitary exponential
3156 : ! **************************************************************************************************
3157 52 : FUNCTION dense_antihermitian_exp(generator) RESULT(rotation)
3158 : COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(IN) :: generator
3159 : COMPLEX(KIND=dp), DIMENSION(SIZE(generator, 1), &
3160 : SIZE(generator, 2)) :: rotation
3161 :
3162 : COMPLEX(KIND=dp), DIMENSION(SIZE(generator, 1), &
3163 52 : SIZE(generator, 2)) :: hermitian, vectors
3164 : INTEGER :: i
3165 52 : REAL(KIND=dp), DIMENSION(SIZE(generator, 1)) :: eigenvalues
3166 :
3167 556 : hermitian = CMPLX(0.0_dp, 1.0_dp, KIND=dp)*generator
3168 52 : CALL diag_complex(hermitian, vectors, eigenvalues)
3169 556 : rotation(:, :) = CMPLX(0.0_dp, 0.0_dp, KIND=dp)
3170 188 : DO i = 1, SIZE(eigenvalues)
3171 : rotation = rotation + EXP(CMPLX(0.0_dp, -eigenvalues(i), KIND=dp))* &
3172 544 : MATMUL(RESHAPE(vectors(:, i), [SIZE(generator, 1), 1]), &
3173 4388 : RESHAPE(CONJG(vectors(:, i)), [1, SIZE(generator, 1)]))
3174 : END DO
3175 :
3176 52 : END FUNCTION dense_antihermitian_exp
3177 :
3178 : ! **************************************************************************************************
3179 : !> \brief Hermitian inverse for an independent dense reference.
3180 : !> \param matrix ...
3181 : !> \param inverse ...
3182 : ! **************************************************************************************************
3183 4 : SUBROUTINE dense_hermitian_inverse(matrix, inverse)
3184 : COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(IN) :: matrix
3185 : COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: inverse
3186 :
3187 : COMPLEX(KIND=dp), &
3188 8 : DIMENSION(SIZE(matrix, 1), SIZE(matrix, 2)) :: vectors
3189 : INTEGER :: i
3190 8 : REAL(KIND=dp), DIMENSION(SIZE(matrix, 1)) :: eigenvalues
3191 :
3192 4 : CALL diag_complex(matrix, vectors, eigenvalues)
3193 52 : inverse(:, :) = CMPLX(0.0_dp, 0.0_dp, KIND=dp)
3194 16 : DO i = 1, SIZE(eigenvalues)
3195 48 : inverse = inverse + MATMUL(RESHAPE(vectors(:, i), [SIZE(matrix, 1), 1]), &
3196 : RESHAPE(CONJG(vectors(:, i)), [1, SIZE(matrix, 1)]))/ &
3197 436 : eigenvalues(i)
3198 : END DO
3199 :
3200 4 : END SUBROUTINE dense_hermitian_inverse
3201 :
3202 : ! **************************************************************************************************
3203 : !> \brief Hermitian inverse square root for the independent dense reference.
3204 : !> \param matrix ...
3205 : !> \param inverse ...
3206 : ! **************************************************************************************************
3207 18 : SUBROUTINE dense_inverse_sqrt(matrix, inverse)
3208 : COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(IN) :: matrix
3209 : COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: inverse
3210 :
3211 : COMPLEX(KIND=dp), &
3212 36 : DIMENSION(SIZE(matrix, 1), SIZE(matrix, 2)) :: vectors
3213 : INTEGER :: i
3214 36 : REAL(KIND=dp), DIMENSION(SIZE(matrix, 1)) :: eigenvalues
3215 :
3216 18 : CALL diag_complex(matrix, vectors, eigenvalues)
3217 126 : inverse(:, :) = CMPLX(0.0_dp, 0.0_dp, KIND=dp)
3218 54 : DO i = 1, SIZE(eigenvalues)
3219 144 : inverse = inverse + MATMUL(RESHAPE(vectors(:, i), [SIZE(matrix, 1), 1]), &
3220 : RESHAPE(CONJG(vectors(:, i)), [1, SIZE(matrix, 1)]))/ &
3221 810 : SQRT(eigenvalues(i))
3222 : END DO
3223 :
3224 18 : END SUBROUTINE dense_inverse_sqrt
3225 :
3226 : ! **************************************************************************************************
3227 : !> \brief Energy of the polar-orthonormalized occupied subspace.
3228 : !> \param coordinate ...
3229 : !> \param hamiltonian ...
3230 : !> \return ...
3231 : ! **************************************************************************************************
3232 4 : FUNCTION polar_energy(coordinate, hamiltonian) RESULT(energy)
3233 : COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(IN) :: coordinate, hamiltonian
3234 : REAL(KIND=dp) :: energy
3235 :
3236 : COMPLEX(KIND=dp), DIMENSION(SIZE(coordinate, 1), &
3237 8 : SIZE(coordinate, 2)) :: orbitals
3238 : COMPLEX(KIND=dp), DIMENSION(SIZE(coordinate, 2), &
3239 8 : SIZE(coordinate, 2)) :: inverse
3240 :
3241 88 : CALL dense_inverse_sqrt(MATMUL(CONJG(TRANSPOSE(coordinate)), coordinate), inverse)
3242 116 : orbitals = MATMUL(coordinate, inverse)
3243 156 : energy = REAL(SUM(CONJG(orbitals)*MATMUL(hamiltonian, orbitals)), KIND=dp)
3244 :
3245 4 : END FUNCTION polar_energy
3246 :
3247 : ! **************************************************************************************************
3248 : !> \brief Independent dense finite STRICT transformation for a Hermitian overlap.
3249 : !> \param c0 reference orbitals
3250 : !> \param coordinate tangent coordinate
3251 : !> \param overlap Hermitian overlap matrix
3252 : !> \param orbitals transformed orbitals
3253 : ! **************************************************************************************************
3254 10 : SUBROUTINE strict_orbitals(c0, coordinate, overlap, orbitals)
3255 : COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(IN) :: c0, coordinate, overlap
3256 : COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: orbitals
3257 :
3258 : COMPLEX(KIND=dp), DIMENSION(SIZE(coordinate, 2), &
3259 20 : SIZE(coordinate, 2)) :: cosp, sinp, vectors
3260 : INTEGER :: i
3261 : REAL(KIND=dp) :: root
3262 20 : REAL(KIND=dp), DIMENSION(SIZE(coordinate, 2)) :: eigenvalues
3263 :
3264 80 : CALL diag_complex(MATMUL(CONJG(TRANSPOSE(coordinate)), MATMUL(overlap, coordinate)), &
3265 580 : vectors, eigenvalues)
3266 70 : cosp(:, :) = CMPLX(0.0_dp, 0.0_dp, KIND=dp)
3267 70 : sinp(:, :) = CMPLX(0.0_dp, 0.0_dp, KIND=dp)
3268 30 : DO i = 1, SIZE(eigenvalues)
3269 20 : root = SQRT(MAX(0.0_dp, eigenvalues(i)))
3270 80 : cosp = cosp + COS(root)*MATMUL(RESHAPE(vectors(:, i), [SIZE(vectors, 1), 1]), &
3271 440 : RESHAPE(CONJG(vectors(:, i)), [1, SIZE(vectors, 1)]))
3272 30 : IF (root < 1.0E-12_dp) THEN
3273 16 : sinp = sinp + MATMUL(RESHAPE(vectors(:, i), [SIZE(vectors, 1), 1]), &
3274 88 : RESHAPE(CONJG(vectors(:, i)), [1, SIZE(vectors, 1)]))
3275 : ELSE
3276 : sinp = sinp + SIN(root)/root* &
3277 64 : MATMUL(RESHAPE(vectors(:, i), [SIZE(vectors, 1), 1]), &
3278 352 : RESHAPE(CONJG(vectors(:, i)), [1, SIZE(vectors, 1)]))
3279 : END IF
3280 : END DO
3281 570 : orbitals = MATMUL(c0, cosp) + MATMUL(coordinate, sinp)
3282 :
3283 10 : END SUBROUTINE strict_orbitals
3284 :
3285 : ! **************************************************************************************************
3286 : !> \brief Energy of the independently transformed finite STRICT orbitals.
3287 : !> \param c0 reference orbitals
3288 : !> \param coordinate tangent coordinate
3289 : !> \param overlap Hermitian overlap matrix
3290 : !> \param hamiltonian Hermitian Hamiltonian
3291 : !> \return energy
3292 : ! **************************************************************************************************
3293 : FUNCTION strict_energy(c0, coordinate, overlap, hamiltonian) RESULT(energy)
3294 : COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(IN) :: c0, coordinate, overlap, hamiltonian
3295 : REAL(KIND=dp) :: energy
3296 :
3297 : COMPLEX(KIND=dp), DIMENSION(SIZE(coordinate, 1), &
3298 : SIZE(coordinate, 2)) :: orbitals
3299 :
3300 : CALL strict_orbitals(c0, coordinate, overlap, orbitals)
3301 : energy = REAL(SUM(CONJG(orbitals)*MATMUL(hamiltonian, orbitals)), KIND=dp)
3302 :
3303 : END FUNCTION strict_energy
3304 :
3305 : ! **************************************************************************************************
3306 : !> \brief Energy of finite STRICT orbitals followed by a complex occupied-space rotation.
3307 : !> \param c0 reference orbitals
3308 : !> \param coordinate tangent coordinate
3309 : !> \param generator anti-Hermitian rotation generator
3310 : !> \param overlap Hermitian overlap matrix
3311 : !> \param hamiltonian Hermitian Hamiltonian
3312 : !> \param occupation orbital occupations
3313 : !> \return energy
3314 : ! **************************************************************************************************
3315 8 : FUNCTION strict_rotated_energy(c0, coordinate, generator, overlap, hamiltonian, &
3316 8 : occupation) RESULT(energy)
3317 : COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(IN) :: c0, coordinate, generator, overlap, &
3318 : hamiltonian
3319 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: occupation
3320 : REAL(KIND=dp) :: energy
3321 :
3322 : COMPLEX(KIND=dp), DIMENSION(SIZE(coordinate, 1), &
3323 16 : SIZE(coordinate, 2)) :: orbitals
3324 : INTEGER :: i
3325 :
3326 8 : CALL strict_orbitals(c0, coordinate, overlap, orbitals)
3327 296 : orbitals = MATMUL(orbitals, dense_antihermitian_exp(generator))
3328 8 : energy = 0.0_dp
3329 24 : DO i = 1, SIZE(occupation)
3330 : energy = energy + occupation(i)* &
3331 296 : REAL(DOT_PRODUCT(orbitals(:, i), MATMUL(hamiltonian, orbitals(:, i))), KIND=dp)
3332 : END DO
3333 :
3334 8 : END FUNCTION strict_rotated_energy
3335 :
3336 : END PROGRAM qs_ot_complex_ref_unittest
|
