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