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