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 Complex k-point orbital transformations, response kernels and local batch execution.
10 : !> Shared real/spectral operations are provided by qs_ot.
11 : ! **************************************************************************************************
12 : MODULE qs_ot_kp
13 : USE cp_cfm_diag, ONLY: cp_cfm_diag_workspace_type,&
14 : cp_cfm_heevd_local,&
15 : cp_cfm_local_workspace_prepare
16 : USE cp_cfm_types, ONLY: cp_cfm_p_type
17 : USE cp_dbcsr_api, ONLY: &
18 : dbcsr_add, dbcsr_copy, dbcsr_create, dbcsr_distribution_type, dbcsr_filter, &
19 : dbcsr_finalize, dbcsr_get_block_p, dbcsr_get_info, dbcsr_iterator_blocks_left, &
20 : dbcsr_iterator_next_block, dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, &
21 : dbcsr_multiply, dbcsr_p_type, dbcsr_release, dbcsr_reserve_blocks, dbcsr_scale, dbcsr_set, &
22 : dbcsr_transposed, dbcsr_type, dbcsr_type_no_symmetry
23 : USE cp_dbcsr_contrib, ONLY: dbcsr_add_on_diag,&
24 : dbcsr_frobenius_norm,&
25 : dbcsr_hadamard_product,&
26 : dbcsr_reserve_all_blocks,&
27 : dbcsr_scale_by_vector
28 : USE cp_dbcsr_diag, ONLY: cp_dbcsr_heevd
29 : USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_local,&
30 : copy_local_to_dbcsr
31 : USE cp_fm_diag, ONLY: FM_DIAG_TYPE_DLAF,&
32 : diag_lib_explicit,&
33 : diag_type
34 : USE cp_fm_types, ONLY: cp_fm_p_type
35 : USE kinds, ONLY: dp
36 : USE local_gemm_api, ONLY: LOCAL_GEMM_PU_GPU,&
37 : LOCAL_GEMM_PU_HOST,&
38 : local_gemm_ctxt_type
39 : USE mathlib, ONLY: diag_complex,&
40 : diamat_all,&
41 : gemm_square
42 : USE preconditioner, ONLY: apply_preconditioner,&
43 : apply_preconditioner_complex_batch
44 : USE preconditioner_types, ONLY: preconditioner_p_type
45 : USE qs_ot, ONLY: qs_ot_complex_multiply,&
46 : qs_ot_p2m_diag,&
47 : qs_ot_project_gradient_real => qs_ot_project_gradient,&
48 : qs_ot_ref_update,&
49 : qs_ot_spectral_kernel,&
50 : qs_ot_spectral_values,&
51 : qs_ot_square_transpose
52 : USE qs_ot_types, ONLY: qs_ot_type
53 :
54 : !$ USE OMP_LIB, ONLY: omp_get_max_threads
55 : #include "./base/base_uses.f90"
56 :
57 : IMPLICIT NONE
58 : PRIVATE
59 :
60 : PUBLIC :: qs_ot_get_p_complex
61 : PUBLIC :: qs_ot_get_orbitals_complex
62 : PUBLIC :: qs_ot_get_orbitals_local_batch
63 : PUBLIC :: qs_ot_get_derivative_complex
64 : PUBLIC :: qs_ot_get_derivative_local_batch
65 : PUBLIC :: qs_ot_prepare_complex_tangent_metric
66 : PUBLIC :: qs_ot_prepare_complex_tangent_metric_batch
67 : PUBLIC :: qs_ot_get_orbitals_ref_complex
68 : PUBLIC :: qs_ot_get_derivative_ref_complex
69 : PUBLIC :: qs_ot_apply_complex_frechet_dbcsr
70 : PUBLIC :: qs_ot_antihermitian_spectral_norm
71 : PUBLIC :: qs_ot_complex_exp_frechet_kernel
72 : PUBLIC :: qs_ot_density_secant_orbital_overlaps
73 : PUBLIC :: qs_ot_density_secant_projected_hessian
74 : PUBLIC :: qs_ot_density_tangent
75 : PUBLIC :: qs_ot_fixed_n_energy_gradient
76 : PUBLIC :: qs_ot_fixed_n_multigroup_schur_block
77 : PUBLIC :: qs_ot_fixed_n_response_mu_shift
78 : PUBLIC :: qs_ot_finite_rotation_response
79 : PUBLIC :: qs_ot_projected_response_update
80 : PUBLIC :: qs_ot_symmetric_abs_solve
81 : PUBLIC :: qs_ot_project_gradient
82 :
83 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_ot_kp'
84 :
85 : INTERFACE qs_ot_complex_multiply_local
86 : MODULE PROCEDURE qs_ot_complex_multiply_pair_local, qs_ot_complex_multiply_parts_local
87 : END INTERFACE
88 :
89 : ! Private rotation spectrum and scratch, reused across response evaluations.
90 : TYPE qs_ot_rotation_work_type
91 : INTEGER :: info = 0
92 : TYPE(cp_cfm_diag_workspace_type) :: solver = cp_cfm_diag_workspace_type()
93 : COMPLEX(dp), ALLOCATABLE :: vectors(:, :), u(:, :), work(:, :), &
94 : inner(:, :)
95 : END TYPE qs_ot_rotation_work_type
96 :
97 : ! Transient per-worker storage. The last index is 1 for real and 2 for imaginary.
98 : TYPE qs_ot_local_work_type
99 : INTEGER :: n = 0, k = 0, info = 0
100 : LOGICAL :: do_rotation = .FALSE., filtered = .FALSE.
101 : TYPE(local_gemm_ctxt_type) :: gemm_ctx = local_gemm_ctxt_type()
102 : TYPE(cp_cfm_diag_workspace_type) :: solver = cp_cfm_diag_workspace_type()
103 : TYPE(qs_ot_rotation_work_type) :: rotation = qs_ot_rotation_work_type()
104 : COMPLEX(dp), ALLOCATABLE :: vectors(:, :)
105 : REAL(dp), ALLOCATABLE :: values(:)
106 : REAL(dp), ALLOCATABLE :: cosp(:, :, :), rot(:, :, :), rot_g(:, :, :)
107 : REAL(dp), ALLOCATABLE :: hc(:, :, :), x(:, :, :), sx(:, :, :), c0(:, :, :), out(:, :, :), &
108 : r(:, :, :), sinp(:, :, :), sc0(:, :, :), TARGET(:, :, :), os(:, :, :), &
109 : buf1(:, :, :), buf2(:, :, :), buf3(:, :, :), buf4(:, :, :), buf_nk(:, :, :), &
110 : buf1_ortho(:, :, :), buf2_ortho(:, :, :), sinp_b(:, :), cosp_b(:, :), tmp(:, :), tmp_ortho(:, :)
111 : END TYPE qs_ot_local_work_type
112 :
113 : CONTAINS
114 :
115 : ! **************************************************************************************************
116 : !> \brief spectral norm of a dense anti-Hermitian rotation generator
117 : !> \param rotation_generator anti-Hermitian generator
118 : !> \return largest absolute eigenvalue of i times the generator
119 : ! **************************************************************************************************
120 180 : FUNCTION qs_ot_antihermitian_spectral_norm(rotation_generator) RESULT(norm)
121 : COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(IN) :: rotation_generator
122 : REAL(KIND=dp) :: norm
123 :
124 180 : COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: eigenvectors
125 : INTEGER :: n
126 180 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eigenvalues
127 :
128 180 : n = SIZE(rotation_generator, 1)
129 180 : CPASSERT(SIZE(rotation_generator, 2) == n)
130 180 : norm = 0.0_dp
131 180 : IF (n == 0) RETURN
132 1080 : ALLOCATE (eigenvectors(n, n), eigenvalues(n))
133 : CALL diag_complex(CMPLX(0.0_dp, 1.0_dp, KIND=dp)*rotation_generator, &
134 8036 : eigenvectors, eigenvalues)
135 1228 : norm = MAXVAL(ABS(eigenvalues))
136 180 : DEALLOCATE (eigenvalues, eigenvectors)
137 :
138 180 : END FUNCTION qs_ot_antihermitian_spectral_norm
139 :
140 : ! **************************************************************************************************
141 : !> \brief chemical-potential response for one fixed-electron-number group
142 : !> \param weighted_energy_response sum_i chi_i de_i over the perturbed local channels
143 : !> \param local_curvature_sum local sum_i chi_i, used as a serial fallback
144 : !> \param fixed_n_curvature_sum global sum_i chi_i for the complete fixed-N group
145 : !> \return first-order chemical-potential shift
146 : ! **************************************************************************************************
147 0 : PURE FUNCTION qs_ot_fixed_n_response_mu_shift( &
148 : weighted_energy_response, local_curvature_sum, fixed_n_curvature_sum) RESULT(mu_shift)
149 : REAL(KIND=dp), INTENT(IN) :: weighted_energy_response, &
150 : local_curvature_sum, &
151 : fixed_n_curvature_sum
152 : REAL(KIND=dp) :: mu_shift
153 :
154 : REAL(KIND=dp) :: denominator
155 :
156 0 : denominator = fixed_n_curvature_sum
157 0 : IF (ABS(denominator) <= EPSILON(denominator)) denominator = local_curvature_sum
158 0 : mu_shift = 0.0_dp
159 0 : IF (ABS(denominator) > EPSILON(denominator)) mu_shift = weighted_energy_response/denominator
160 :
161 0 : END FUNCTION qs_ot_fixed_n_response_mu_shift
162 :
163 : ! **************************************************************************************************
164 : !> \brief fixed-N Mermin gradient in auxiliary-energy coordinates
165 : !> \param rayleigh_energy diagonal expectation values of the current Hamiltonian
166 : !> \param energy_coordinate auxiliary band energies controlling the occupations
167 : !> \param response_weight signed weighted occupation responses chi_i
168 : !> \param fixed_n_weight_sum global sum_i chi_i for the fixed-N group
169 : !> \param fixed_n_weighted_residual global sum_i chi_i (h_i-e_i)
170 : !> \param gradient projected fixed-N gradient
171 : ! **************************************************************************************************
172 2096 : PURE SUBROUTINE qs_ot_fixed_n_energy_gradient( &
173 2096 : rayleigh_energy, energy_coordinate, response_weight, fixed_n_weight_sum, &
174 2096 : fixed_n_weighted_residual, gradient)
175 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: rayleigh_energy, energy_coordinate, &
176 : response_weight
177 : REAL(KIND=dp), INTENT(IN) :: fixed_n_weight_sum, &
178 : fixed_n_weighted_residual
179 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: gradient
180 :
181 : REAL(KIND=dp) :: fixed_n_mean
182 :
183 2096 : fixed_n_mean = 0.0_dp
184 2096 : IF (ABS(fixed_n_weight_sum) > EPSILON(fixed_n_weight_sum)) THEN
185 2052 : fixed_n_mean = fixed_n_weighted_residual/fixed_n_weight_sum
186 : END IF
187 : gradient(:) = response_weight(:)* &
188 16360 : (fixed_n_mean - (rayleigh_energy(:) - energy_coordinate(:)))
189 :
190 2096 : END SUBROUTINE qs_ot_fixed_n_energy_gradient
191 :
192 : ! **************************************************************************************************
193 : !> \brief eliminate spin-resolved auxiliary energies while retaining every fixed-N constraint
194 : !>
195 : !> Each energy variable belongs to one particle-number group. The diagonal occupation
196 : !> response is eliminated locally, while one coupling column per group is retained for the
197 : !> subsequent global chemical-potential projection. A single group reduces exactly to
198 : !> the single-group Schur complement A - R^T (D - chi chi^T/sum(chi)) R.
199 : !> \param rotation_hessian finite orbital-rotation Hessian
200 : !> \param rayleigh_response derivative of all spin-resolved Rayleigh energies
201 : !> \param response_weight positive occupation-response weights
202 : !> \param response_group fixed-N group of every energy variable
203 : !> \param rotation_gradient orbital-rotation gradient
204 : !> \param energy_gradient auxiliary-energy gradient
205 : !> \param schur_block energy-eliminated rotation block
206 : !> \param coupling_matrix retained chemical-potential coupling columns
207 : !> \param schur_rhs energy-eliminated rotation right-hand side
208 : ! **************************************************************************************************
209 346 : SUBROUTINE qs_ot_fixed_n_multigroup_schur_block( &
210 692 : rotation_hessian, rayleigh_response, response_weight, response_group, &
211 692 : rotation_gradient, energy_gradient, schur_block, coupling_matrix, schur_rhs)
212 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: rotation_hessian, rayleigh_response
213 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: response_weight
214 : INTEGER, DIMENSION(:), INTENT(IN) :: response_group
215 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: rotation_gradient, energy_gradient
216 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: schur_block, coupling_matrix
217 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: schur_rhs
218 :
219 : INTEGER :: group, i, nenergy, ngroups, nrotation
220 : REAL(KIND=dp) :: weight
221 :
222 346 : nenergy = SIZE(response_weight)
223 346 : nrotation = SIZE(rotation_gradient)
224 346 : ngroups = SIZE(coupling_matrix, 2)
225 346 : CPASSERT(ngroups > 0)
226 1038 : CPASSERT(ALL(SHAPE(rotation_hessian) == [nrotation, nrotation]))
227 1038 : CPASSERT(ALL(SHAPE(rayleigh_response) == [nenergy, nrotation]))
228 346 : CPASSERT(SIZE(energy_gradient) == nenergy)
229 346 : CPASSERT(SIZE(response_group) == nenergy)
230 2350 : CPASSERT(ALL(response_group >= 1 .AND. response_group <= ngroups))
231 1038 : CPASSERT(ALL(SHAPE(schur_block) == [nrotation, nrotation]))
232 346 : CPASSERT(SIZE(coupling_matrix, 1) == nrotation)
233 346 : CPASSERT(SIZE(schur_rhs) == nrotation)
234 :
235 530038 : schur_block(:, :) = rotation_hessian(:, :)
236 12674 : coupling_matrix(:, :) = 0.0_dp
237 2350 : DO i = 1, nenergy
238 2004 : weight = response_weight(i)
239 : schur_block(:, :) = schur_block(:, :) - weight* &
240 : SPREAD(rayleigh_response(i, :), DIM=2, NCOPIES=nrotation)* &
241 5766876 : SPREAD(rayleigh_response(i, :), DIM=1, NCOPIES=nrotation)
242 2004 : group = response_group(i)
243 : coupling_matrix(:, group) = coupling_matrix(:, group) + &
244 87910 : weight*rayleigh_response(i, :)
245 : END DO
246 1060076 : schur_block(:, :) = 0.5_dp*(schur_block + TRANSPOSE(schur_block))
247 : schur_rhs(:) = rotation_gradient + &
248 10568 : MATMUL(TRANSPOSE(rayleigh_response), energy_gradient)
249 :
250 346 : END SUBROUTINE qs_ot_fixed_n_multigroup_schur_block
251 :
252 : ! **************************************************************************************************
253 : !> \brief apply a positive spectral inverse of a real symmetric response matrix
254 : !>
255 : !> The magnitude of every resolved eigenmode is retained, including modes with negative
256 : !> physical curvature. Replacing lambda by abs(lambda) gives a descent metric without the
257 : !> loss of response information caused by discarding the negative subspace. Unresolved
258 : !> null modes are projected out instead of being amplified by an artificial eigenvalue floor.
259 : !> \param matrix real symmetric response matrix
260 : !> \param rhs one or more right-hand sides
261 : !> \param solution spectral-absolute inverse applied to rhs
262 : !> \param valid whether finite input and output were obtained
263 : !> \param relative_floor optional relative eigenvalue floor
264 : ! **************************************************************************************************
265 346 : SUBROUTINE qs_ot_symmetric_abs_solve(matrix, rhs, solution, valid, relative_floor)
266 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: matrix, rhs
267 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: solution
268 : LOGICAL, INTENT(OUT) :: valid
269 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: relative_floor
270 :
271 : INTEGER :: i, n, nresolved
272 : REAL(KIND=dp) :: eigenvalue_floor, relative_floor_eff, &
273 : scale
274 346 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eigenvalues
275 346 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: eigenvectors, work
276 :
277 346 : n = SIZE(matrix, 1)
278 1038 : CPASSERT(ALL(SHAPE(matrix) == [n, n]))
279 346 : CPASSERT(SIZE(rhs, 1) == n)
280 1038 : CPASSERT(ALL(SHAPE(solution) == SHAPE(rhs)))
281 552934 : valid = ALL(matrix == matrix) .AND. ALL(rhs == rhs)
282 346 : IF (.NOT. valid) THEN
283 0 : solution(:, :) = 0.0_dp
284 : RETURN
285 : END IF
286 :
287 346 : relative_floor_eff = SQRT(EPSILON(1.0_dp))
288 346 : IF (PRESENT(relative_floor)) relative_floor_eff = MAX(relative_floor, relative_floor_eff)
289 3114 : ALLOCATE (eigenvalues(n), eigenvectors(n, n), work(n, SIZE(rhs, 2)))
290 530038 : eigenvectors(:, :) = 0.5_dp*(matrix + TRANSPOSE(matrix))
291 346 : CALL diamat_all(eigenvectors, eigenvalues)
292 10222 : scale = MAX(1.0_dp, MAXVAL(ABS(eigenvalues)))
293 346 : eigenvalue_floor = relative_floor_eff*scale
294 2356498 : work(:, :) = MATMUL(TRANSPOSE(eigenvectors), rhs)
295 346 : nresolved = 0
296 10222 : DO i = 1, n
297 10222 : IF (ABS(eigenvalues(i)) > eigenvalue_floor) THEN
298 16899 : work(i, :) = work(i, :)/ABS(eigenvalues(i))
299 5289 : nresolved = nresolved + 1
300 : ELSE
301 14757 : work(i, :) = 0.0_dp
302 : END IF
303 : END DO
304 2369942 : solution(:, :) = MATMUL(eigenvectors, work)
305 22896 : valid = nresolved > 0 .AND. ALL(solution == solution)
306 634 : IF (.NOT. valid) solution(:, :) = 0.0_dp
307 346 : DEALLOCATE (eigenvalues, eigenvectors, work)
308 :
309 346 : END SUBROUTINE qs_ot_symmetric_abs_solve
310 :
311 : ! **************************************************************************************************
312 : !> \brief update a baseline response direction in a small positive physical-response subspace
313 : !>
314 : !> The first basis mode is the baseline response direction. With B0 the projected frozen-H
315 : !> Hessian and K the accepted physical correction, this routine solves
316 : !>
317 : !> (B0 + K) c = g_Q.
318 : !>
319 : !> The optional projected physical gradient supplies g_Q. Without it, g_Q=B0*e1, so c=e1
320 : !> exactly when K=0. An indefinite or unresolved total projected Hessian is rejected instead
321 : !> of turning its negative modes into an unrelated active direction.
322 : !> \param reference_hessian projected frozen-H Hessian B0
323 : !> \param response_correction projected physical response K
324 : !> \param coefficients response coefficients c in the supplied basis
325 : !> \param valid whether a finite, positive, sufficiently resolved solve was obtained
326 : !> \param projected_gradient optional physical gradient projected onto the supplied basis
327 : !> \param relative_floor optional relative positive-eigenvalue floor
328 : ! **************************************************************************************************
329 138 : SUBROUTINE qs_ot_projected_response_update( &
330 138 : reference_hessian, response_correction, coefficients, valid, projected_gradient, relative_floor)
331 :
332 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: reference_hessian, response_correction
333 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: coefficients
334 : LOGICAL, INTENT(OUT) :: valid
335 : REAL(KIND=dp), DIMENSION(:), INTENT(IN), OPTIONAL :: projected_gradient
336 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: relative_floor
337 :
338 : INTEGER :: i, n
339 : REAL(KIND=dp) :: eigenvalue_floor, relative_floor_eff, &
340 : scale
341 138 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eigenvalues, rhs, work
342 138 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: eigenvectors
343 :
344 138 : n = SIZE(reference_hessian, 1)
345 138 : CPASSERT(n > 0)
346 414 : CPASSERT(ALL(SHAPE(reference_hessian) == [n, n]))
347 414 : CPASSERT(ALL(SHAPE(response_correction) == [n, n]))
348 138 : CPASSERT(SIZE(coefficients) == n)
349 138 : IF (PRESENT(projected_gradient)) THEN
350 138 : CPASSERT(SIZE(projected_gradient) == n)
351 : END IF
352 414 : coefficients(:) = 0.0_dp
353 : valid = ALL(reference_hessian == reference_hessian) .AND. &
354 1932 : ALL(response_correction == response_correction)
355 138 : IF (PRESENT(projected_gradient)) THEN
356 414 : valid = valid .AND. ALL(projected_gradient == projected_gradient)
357 : END IF
358 138 : IF (.NOT. valid) RETURN
359 :
360 138 : relative_floor_eff = 1.0E-4_dp
361 138 : IF (PRESENT(relative_floor)) relative_floor_eff = &
362 0 : MAX(relative_floor, SQRT(EPSILON(1.0_dp)))
363 1104 : ALLOCATE (eigenvalues(n), eigenvectors(n, n), rhs(n), work(n))
364 : eigenvectors(:, :) = 0.5_dp* &
365 : (reference_hessian + response_correction + &
366 966 : TRANSPOSE(reference_hessian + response_correction))
367 138 : CALL diamat_all(eigenvectors, eigenvalues)
368 1242 : scale = MAX(MAXVAL(ABS(eigenvalues)), MAXVAL(ABS(reference_hessian)))
369 138 : IF (scale <= TINY(scale)) THEN
370 0 : valid = .FALSE.
371 0 : coefficients(:) = 0.0_dp
372 0 : DEALLOCATE (eigenvalues, eigenvectors, rhs, work)
373 0 : RETURN
374 : END IF
375 138 : eigenvalue_floor = relative_floor_eff*scale
376 406 : valid = ALL(eigenvalues > eigenvalue_floor)
377 138 : IF (valid) THEN
378 402 : rhs(:) = reference_hessian(:, 1)
379 402 : IF (PRESENT(projected_gradient)) rhs(:) = projected_gradient
380 134 : work(:) = MATMUL(TRANSPOSE(eigenvectors), rhs)
381 402 : DO i = 1, n
382 402 : work(i) = work(i)/eigenvalues(i)
383 : END DO
384 1742 : coefficients(:) = MATMUL(eigenvectors, work)
385 : valid = ALL(coefficients == coefficients) .AND. &
386 804 : ALL(ABS(coefficients) <= HUGE(1.0_dp))
387 : END IF
388 146 : IF (.NOT. valid) coefficients(:) = 0.0_dp
389 138 : DEALLOCATE (eigenvalues, eigenvectors, rhs, work)
390 :
391 138 : END SUBROUTINE qs_ot_projected_response_update
392 :
393 : ! **************************************************************************************************
394 : !> \brief form a projected self-adjoint Hxc response from distributed density-space overlaps
395 : !> \param density_norm_sq <S,S>
396 : !> \param response_work <S,Y>
397 : !> \param density_overlap <S,B_r>
398 : !> \param response_overlap <Y,B_r>
399 : !> \param correction projected symmetric Hxc response <B_q,K*B_r>
400 : !> \param valid whether finite nonzero secant data were available
401 : !> \param secant_mode optional mode representing the accepted full density secant divided by its
402 : !> line-search position
403 : !> \param secant_position signed line-search position of the accepted full density secant
404 : ! **************************************************************************************************
405 184 : SUBROUTINE qs_ot_density_secant_projected_hessian( &
406 184 : density_norm_sq, response_work, density_overlap, response_overlap, correction, valid, &
407 : secant_mode, secant_position)
408 :
409 : REAL(KIND=dp), INTENT(IN) :: density_norm_sq, response_work
410 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: density_overlap, response_overlap
411 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: correction
412 : LOGICAL, INTENT(OUT) :: valid
413 : INTEGER, INTENT(IN), OPTIONAL :: secant_mode
414 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: secant_position
415 :
416 : INTEGER :: i, j, nmode
417 : REAL(KIND=dp) :: inverse_density_norm
418 368 : REAL(KIND=dp), DIMENSION(SIZE(density_overlap)) :: projected_density_overlap, &
419 184 : projected_response_overlap
420 :
421 184 : nmode = SIZE(density_overlap)
422 184 : CPASSERT(SIZE(response_overlap) == nmode)
423 552 : CPASSERT(ALL(SHAPE(correction) == [nmode, nmode]))
424 :
425 1288 : correction(:, :) = 0.0_dp
426 184 : valid = .FALSE.
427 184 : IF (nmode <= 0 .OR. density_norm_sq <= TINY(1.0_dp)) RETURN
428 184 : IF (density_norm_sq /= density_norm_sq .OR. response_work /= response_work) RETURN
429 184 : IF (ABS(density_norm_sq) > HUGE(1.0_dp) .OR. ABS(response_work) > HUGE(1.0_dp)) RETURN
430 1104 : IF (ANY(density_overlap /= density_overlap) .OR. ANY(response_overlap /= response_overlap)) RETURN
431 1104 : IF (ANY(ABS(density_overlap) > HUGE(1.0_dp)) .OR. &
432 : ANY(ABS(response_overlap) > HUGE(1.0_dp))) RETURN
433 184 : IF (PRESENT(secant_mode) .NEQV. PRESENT(secant_position)) RETURN
434 :
435 552 : projected_density_overlap(:) = density_overlap
436 552 : projected_response_overlap(:) = response_overlap
437 184 : IF (PRESENT(secant_mode)) THEN
438 138 : IF (secant_mode < 1 .OR. secant_mode > nmode) RETURN
439 : IF (secant_position /= secant_position .OR. &
440 138 : ABS(secant_position) > HUGE(1.0_dp) .OR. &
441 138 : ABS(secant_position) <= SQRT(EPSILON(1.0_dp))) RETURN
442 : ! The complete accepted density change can contain REF-orbital motion that is absent from
443 : ! a reduced rotation/occupation tangent. S/alpha supplies its exact projected secant mode.
444 138 : projected_density_overlap(secant_mode) = density_norm_sq/secant_position
445 138 : projected_response_overlap(secant_mode) = response_work/secant_position
446 : END IF
447 :
448 184 : inverse_density_norm = 1.0_dp/density_norm_sq
449 552 : DO j = 1, nmode
450 1288 : DO i = 1, nmode
451 : correction(i, j) = inverse_density_norm* &
452 : (projected_response_overlap(i)*projected_density_overlap(j) + &
453 : projected_density_overlap(i)*projected_response_overlap(j) - &
454 : response_work*inverse_density_norm* &
455 1104 : projected_density_overlap(i)*projected_density_overlap(j))
456 : END DO
457 : END DO
458 2392 : correction(:, :) = 0.5_dp*(correction + TRANSPOSE(correction))
459 2576 : valid = ALL(correction == correction) .AND. ALL(ABS(correction) <= HUGE(1.0_dp))
460 184 : IF (.NOT. valid) correction(:, :) = 0.0_dp
461 :
462 : END SUBROUTINE qs_ot_density_secant_projected_hessian
463 :
464 : ! **************************************************************************************************
465 : !> \brief finite-chart density tangent for coupled complex rotations and fixed-N occupations
466 : !>
467 : !> The rotation contribution differentiates
468 : !>
469 : !> exp(X) diag(w_k f) exp(X)^H
470 : !>
471 : !> along an anti-Hermitian packed direction. The supplied weighted occupation response is
472 : !> added in the same chart, and the result is returned in the current physical orbital basis.
473 : !> \param rotation_generator current anti-Hermitian REF generator X
474 : !> \param occupation current occupations f
475 : !> \param kpoint_weight irreducible K-point weight w_k
476 : !> \param rotation_step interleaved real/imaginary anti-Hermitian direction
477 : !> \param weighted_occupation_step derivative of w_k*f, including the fixed-N mu response
478 : !> \param density_tangent Hermitian tangent in the current physical orbital basis
479 : !> \param difference_step optional finite-chart central-difference step
480 : ! **************************************************************************************************
481 496 : SUBROUTINE qs_ot_density_tangent( &
482 496 : rotation_generator, occupation, kpoint_weight, rotation_step, &
483 496 : weighted_occupation_step, density_tangent, difference_step)
484 :
485 : COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(IN) :: rotation_generator
486 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: occupation
487 : REAL(KIND=dp), INTENT(IN) :: kpoint_weight
488 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: rotation_step, weighted_occupation_step
489 : COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: density_tangent
490 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: difference_step
491 :
492 496 : COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: direction, generator_minus, &
493 496 : generator_plus, mode_ref, rotation, rotation_minus, rotation_plus, weighted_rotation
494 : INTEGER :: i, j, n, nrotation, r
495 : REAL(KIND=dp) :: step, step_scale
496 : TYPE(local_gemm_ctxt_type) :: gemm_ctx
497 496 : TYPE(qs_ot_rotation_work_type) :: rotation_work
498 :
499 496 : n = SIZE(rotation_generator, 1)
500 496 : nrotation = n*(n - 1)
501 0 : CPASSERT(n > 0)
502 1488 : CPASSERT(ALL(SHAPE(rotation_generator) == [n, n]))
503 496 : CPASSERT(SIZE(occupation) == n)
504 496 : CPASSERT(SIZE(rotation_step) == nrotation)
505 496 : CPASSERT(SIZE(weighted_occupation_step) == n)
506 1488 : CPASSERT(ALL(SHAPE(density_tangent) == [n, n]))
507 496 : CPASSERT(kpoint_weight > 0.0_dp)
508 :
509 : ALLOCATE (direction(n, n), generator_minus(n, n), generator_plus(n, n), &
510 : mode_ref(n, n), rotation(n, n), rotation_minus(n, n), &
511 8928 : rotation_plus(n, n), weighted_rotation(n, n))
512 496 : direction(:, :) = CMPLX(0.0_dp, 0.0_dp, KIND=dp)
513 496 : r = 0
514 3094 : DO i = 1, n - 1
515 12430 : DO j = i + 1, n
516 9336 : r = r + 1
517 9336 : direction(i, j) = CMPLX(rotation_step(r), 0.0_dp, KIND=dp)
518 9336 : direction(j, i) = -direction(i, j)
519 9336 : r = r + 1
520 : direction(i, j) = direction(i, j) + &
521 9336 : CMPLX(0.0_dp, rotation_step(r), KIND=dp)
522 : direction(j, i) = direction(j, i) + &
523 11934 : CMPLX(0.0_dp, rotation_step(r), KIND=dp)
524 : END DO
525 : END DO
526 496 : CPASSERT(r == nrotation)
527 :
528 25356 : step_scale = MAX(1.0_dp, MAXVAL(ABS(direction)))
529 496 : step = 1.0E-5_dp/step_scale
530 496 : IF (PRESENT(difference_step)) step = difference_step/step_scale
531 496 : CPASSERT(step > EPSILON(step))
532 25356 : generator_plus(:, :) = rotation_generator + step*direction
533 25356 : generator_minus(:, :) = rotation_generator - step*direction
534 496 : CALL qs_ot_prepare_rotation_work(rotation_work, n)
535 496 : CALL gemm_ctx%create(LOCAL_GEMM_PU_HOST, timing=.FALSE.)
536 496 : CALL qs_ot_rotation_state_local(rotation_generator, rotation_work, gemm_ctx)
537 496 : CPASSERT(rotation_work%info == 0)
538 25356 : rotation(:, :) = rotation_work%u(:, :)
539 496 : CALL qs_ot_rotation_state_local(generator_plus, rotation_work, gemm_ctx)
540 496 : CPASSERT(rotation_work%info == 0)
541 25356 : rotation_plus(:, :) = rotation_work%u(:, :)
542 496 : CALL qs_ot_rotation_state_local(generator_minus, rotation_work, gemm_ctx)
543 496 : CPASSERT(rotation_work%info == 0)
544 25356 : rotation_minus(:, :) = rotation_work%u(:, :)
545 496 : CALL gemm_ctx%destroy()
546 :
547 25356 : weighted_rotation(:, :) = rotation_plus
548 3590 : DO j = 1, n
549 25356 : weighted_rotation(:, j) = kpoint_weight*occupation(j)*weighted_rotation(:, j)
550 : END DO
551 195118 : mode_ref(:, :) = MATMUL(weighted_rotation, CONJG(TRANSPOSE(rotation_plus)))
552 25356 : weighted_rotation(:, :) = rotation_minus
553 3590 : DO j = 1, n
554 25356 : weighted_rotation(:, j) = kpoint_weight*occupation(j)*weighted_rotation(:, j)
555 : END DO
556 : mode_ref(:, :) = (mode_ref - &
557 221466 : MATMUL(weighted_rotation, CONJG(TRANSPOSE(rotation_minus))))/(2.0_dp*step)
558 :
559 25356 : weighted_rotation(:, :) = rotation
560 3590 : DO j = 1, n
561 25356 : weighted_rotation(:, j) = weighted_occupation_step(j)*weighted_rotation(:, j)
562 : END DO
563 221962 : mode_ref(:, :) = mode_ref + MATMUL(weighted_rotation, CONJG(TRANSPOSE(rotation)))
564 585850 : density_tangent(:, :) = MATMUL(CONJG(TRANSPOSE(rotation)), MATMUL(mode_ref, rotation))
565 50216 : density_tangent(:, :) = 0.5_dp*(density_tangent + CONJG(TRANSPOSE(density_tangent)))
566 :
567 0 : DEALLOCATE (direction, generator_minus, generator_plus, mode_ref, rotation, &
568 496 : rotation_minus, rotation_plus, weighted_rotation)
569 :
570 992 : END SUBROUTINE qs_ot_density_tangent
571 :
572 : ! **************************************************************************************************
573 : !> \brief project a physical density/Hamiltonian secant between moving orbital subspaces
574 : !>
575 : !> For separately S-orthonormal endpoint orbitals C0 and C1, O=C0^H*S*C1 retains the
576 : !> component of the accepted density step that leaves the old subspace. Density modes are
577 : !> represented in the current C1 basis and already contain the irreducible K-point weight.
578 : !> \param overlap_start_current cross overlap O
579 : !> \param occupation_start occupations at the accepted start
580 : !> \param occupation_current occupations at the accepted endpoint
581 : !> \param hamiltonian_step_start C0^H*(H1-H0)*C0
582 : !> \param hamiltonian_step_current C1^H*(H1-H0)*C1
583 : !> \param density_modes current-orbital density tangents
584 : !> \param kpoint_weight irreducible K-point weight
585 : !> \param density_norm_sq contribution to <Delta P,Delta P>
586 : !> \param response_work contribution to <Delta P,Delta H>
587 : !> \param density_overlap contributions <Delta P,B_r>
588 : !> \param response_overlap contributions <Delta H,B_r>
589 : !> \param valid whether a finite nonzero secant was available
590 : ! **************************************************************************************************
591 284 : SUBROUTINE qs_ot_density_secant_orbital_overlaps( &
592 568 : overlap_start_current, occupation_start, occupation_current, hamiltonian_step_start, &
593 568 : hamiltonian_step_current, density_modes, kpoint_weight, density_norm_sq, response_work, &
594 284 : density_overlap, response_overlap, valid)
595 :
596 : COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(IN) :: overlap_start_current
597 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: occupation_start, occupation_current
598 : COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(IN) :: hamiltonian_step_start, &
599 : hamiltonian_step_current
600 : COMPLEX(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: density_modes
601 : REAL(KIND=dp), INTENT(IN) :: kpoint_weight
602 : REAL(KIND=dp), INTENT(OUT) :: density_norm_sq, response_work
603 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: density_overlap, response_overlap
604 : LOGICAL, INTENT(OUT) :: valid
605 :
606 284 : COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: cross_mode
607 : INTEGER :: i, j, mode, n, nmode
608 : REAL(KIND=dp) :: cross_density, density_scale
609 :
610 284 : n = SIZE(occupation_start)
611 284 : nmode = SIZE(density_modes, 3)
612 284 : CPASSERT(n > 0)
613 284 : CPASSERT(SIZE(occupation_current) == n)
614 852 : CPASSERT(ALL(SHAPE(overlap_start_current) == [n, n]))
615 852 : CPASSERT(ALL(SHAPE(hamiltonian_step_start) == [n, n]))
616 852 : CPASSERT(ALL(SHAPE(hamiltonian_step_current) == [n, n]))
617 284 : CPASSERT(SIZE(density_modes, 1) == n .AND. SIZE(density_modes, 2) == n)
618 284 : CPASSERT(SIZE(density_overlap) == nmode .AND. SIZE(response_overlap) == nmode)
619 :
620 284 : density_norm_sq = 0.0_dp
621 284 : response_work = 0.0_dp
622 852 : density_overlap(:) = 0.0_dp
623 852 : response_overlap(:) = 0.0_dp
624 284 : valid = .FALSE.
625 292 : IF (nmode <= 0 .OR. kpoint_weight <= TINY(1.0_dp)) RETURN
626 :
627 : cross_density = 0.0_dp
628 2044 : DO j = 1, n
629 14364 : DO i = 1, n
630 : cross_density = cross_density + occupation_start(i)*occupation_current(j)* &
631 14080 : ABS(overlap_start_current(i, j))**2
632 : END DO
633 : END DO
634 : density_norm_sq = kpoint_weight* &
635 : (SUM(occupation_start**2) + SUM(occupation_current**2) - &
636 3804 : 2.0_dp*cross_density)
637 : response_work = 0.0_dp
638 2044 : DO i = 1, n
639 : response_work = response_work + kpoint_weight* &
640 : (occupation_current(i)*REAL(hamiltonian_step_current(i, i), KIND=dp) - &
641 2044 : occupation_start(i)*REAL(hamiltonian_step_start(i, i), KIND=dp))
642 : END DO
643 :
644 1136 : ALLOCATE (cross_mode(n, n))
645 852 : DO mode = 1, nmode
646 : ! GCC 15 miscompiles the equivalent complex array expressions with -march=znver3.
647 : CALL gemm_square(overlap_start_current, "N", density_modes(:, :, mode), "N", &
648 568 : overlap_start_current, "C", cross_mode)
649 4088 : DO i = 1, n
650 : density_overlap(mode) = density_overlap(mode) + &
651 : occupation_current(i)* &
652 : REAL(density_modes(i, i, mode), KIND=dp) - &
653 4088 : occupation_start(i)*REAL(cross_mode(i, i), KIND=dp)
654 : END DO
655 : response_overlap(mode) = &
656 : SUM(REAL(hamiltonian_step_current, KIND=dp)* &
657 : REAL(density_modes(:, :, mode), KIND=dp) + &
658 29012 : AIMAG(hamiltonian_step_current)*AIMAG(density_modes(:, :, mode)))
659 : END DO
660 284 : DEALLOCATE (cross_mode)
661 :
662 : density_scale = kpoint_weight* &
663 3804 : MAX(MAXVAL(ABS(occupation_start)), MAXVAL(ABS(occupation_current)))
664 284 : IF (density_scale <= TINY(1.0_dp)) RETURN
665 284 : IF (density_norm_sq <= 64.0_dp*EPSILON(1.0_dp)*density_scale**2) RETURN
666 : valid = density_norm_sq == density_norm_sq .AND. response_work == response_work .AND. &
667 : ABS(density_norm_sq) <= HUGE(1.0_dp) .AND. ABS(response_work) <= HUGE(1.0_dp) .AND. &
668 : ALL(density_overlap == density_overlap) .AND. &
669 : ALL(response_overlap == response_overlap) .AND. &
670 : ALL(ABS(density_overlap) <= HUGE(1.0_dp)) .AND. &
671 3312 : ALL(ABS(response_overlap) <= HUGE(1.0_dp))
672 276 : IF (.NOT. valid) THEN
673 0 : density_norm_sq = 0.0_dp
674 0 : response_work = 0.0_dp
675 0 : density_overlap(:) = 0.0_dp
676 0 : response_overlap(:) = 0.0_dp
677 : END IF
678 :
679 8 : END SUBROUTINE qs_ot_density_secant_orbital_overlaps
680 :
681 : ! **************************************************************************************************
682 : !> \brief finite complex REF rotation Hessian and Rayleigh-energy response
683 : !>
684 : !> The current projected Hamiltonian is pulled back through the finite rotation and then
685 : !> differentiated in the independent real-antisymmetric and imaginary-symmetric pair
686 : !> coordinates. This keeps the response consistent with the exponential chart used by REF
687 : !> OT instead of replacing it by an infinitesimal commutator away from the chart origin.
688 : !> \param chc current projected Hermitian Hamiltonian U^H H_ref U
689 : !> \param rotation_generator current anti-Hermitian REF generator
690 : !> \param occupation fixed occupations attached to the rotated columns
691 : !> \param kpoint_weight irreducible-k-point weight
692 : !> \param rotation_gradient gradient in interleaved real/imaginary pair coordinates
693 : !> \param rotation_hessian derivative of rotation_gradient in the same coordinates
694 : !> \param rayleigh_response derivative of diag(U^H H_ref U) with respect to the pair coordinates
695 : !> \param difference_step optional central finite-difference step for the Hessian action
696 : ! **************************************************************************************************
697 362 : SUBROUTINE qs_ot_finite_rotation_response( &
698 362 : chc, rotation_generator, occupation, kpoint_weight, rotation_gradient, &
699 362 : rotation_hessian, rayleigh_response, difference_step)
700 : COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(IN) :: chc, rotation_generator
701 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: occupation
702 : REAL(KIND=dp), INTENT(IN) :: kpoint_weight
703 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: rotation_gradient
704 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: rotation_hessian, rayleigh_response
705 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: difference_step
706 :
707 : COMPLEX(KIND=dp) :: perturbation
708 362 : COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: base_hamiltonian, generator_minus, &
709 362 : generator_plus
710 : INTEGER :: component, i, j, n, nrotation, r, s
711 : REAL(KIND=dp) :: step
712 362 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: occupation_scale, rayleigh, &
713 362 : rayleigh_minus, rayleigh_plus
714 362 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: gradient_imag, gradient_minus_imag, &
715 362 : gradient_minus_real, gradient_plus_imag, gradient_plus_real, gradient_real
716 : TYPE(local_gemm_ctxt_type) :: gemm_ctx
717 362 : TYPE(qs_ot_rotation_work_type) :: rotation_work
718 :
719 362 : n = SIZE(chc, 1)
720 362 : nrotation = n*(n - 1)
721 0 : CPASSERT(n > 0)
722 1086 : CPASSERT(ALL(SHAPE(chc) == [n, n]))
723 1086 : CPASSERT(ALL(SHAPE(rotation_generator) == [n, n]))
724 362 : CPASSERT(SIZE(occupation) == n)
725 362 : CPASSERT(SIZE(rotation_gradient) == nrotation)
726 1086 : CPASSERT(ALL(SHAPE(rotation_hessian) == [nrotation, nrotation]))
727 1086 : CPASSERT(ALL(SHAPE(rayleigh_response) == [n, nrotation]))
728 362 : CPASSERT(kpoint_weight > 0.0_dp)
729 :
730 362 : step = 1.0E-4_dp
731 362 : IF (PRESENT(difference_step)) step = difference_step
732 362 : CPASSERT(step > SQRT(EPSILON(step)))
733 :
734 : ALLOCATE (base_hamiltonian(n, n), generator_minus(n, n), generator_plus(n, n), &
735 : occupation_scale(n), rayleigh(n), rayleigh_minus(n), &
736 : rayleigh_plus(n), gradient_imag(n, n), gradient_minus_imag(n, n), &
737 : gradient_minus_real(n, n), gradient_plus_imag(n, n), &
738 9412 : gradient_plus_real(n, n), gradient_real(n, n))
739 :
740 362 : CALL qs_ot_prepare_rotation_work(rotation_work, n)
741 362 : CALL gemm_ctx%create(LOCAL_GEMM_PU_HOST, timing=.FALSE.)
742 362 : CALL qs_ot_rotation_state_local(rotation_generator, rotation_work, gemm_ctx)
743 362 : CPASSERT(rotation_work%info == 0)
744 : CALL gemm_ctx%gemm('N', 'C', n, n, n, (1.0_dp, 0.0_dp), chc, n, rotation_work%u, n, &
745 362 : (0.0_dp, 0.0_dp), rotation_work%work, n)
746 : CALL gemm_ctx%gemm('N', 'N', n, n, n, (1.0_dp, 0.0_dp), rotation_work%u, n, rotation_work%work, n, &
747 362 : (0.0_dp, 0.0_dp), base_hamiltonian, n)
748 2370 : occupation_scale(:) = 2.0_dp*kpoint_weight*occupation(:)
749 : CALL qs_ot_dense_rotation_gradient(base_hamiltonian, occupation_scale, &
750 362 : gradient_real, gradient_imag, rayleigh, rotation_work, gemm_ctx)
751 :
752 362 : r = 0
753 2008 : DO i = 1, n - 1
754 7578 : DO j = i + 1, n
755 5570 : r = r + 1
756 5570 : rotation_gradient(r) = gradient_real(i, j)
757 5570 : r = r + 1
758 7216 : rotation_gradient(r) = gradient_imag(i, j)
759 : END DO
760 : END DO
761 362 : CPASSERT(r == nrotation)
762 :
763 : s = 0
764 2008 : DO i = 1, n - 1
765 7578 : DO j = i + 1, n
766 18356 : DO component = 1, 2
767 11140 : perturbation = CMPLX(step, 0.0_dp, KIND=dp)
768 11140 : IF (component == 2) perturbation = CMPLX(0.0_dp, step, KIND=dp)
769 815500 : generator_plus(:, :) = rotation_generator(:, :)
770 815500 : generator_minus(:, :) = rotation_generator(:, :)
771 11140 : generator_plus(i, j) = generator_plus(i, j) + perturbation
772 11140 : generator_plus(j, i) = generator_plus(j, i) - CONJG(perturbation)
773 11140 : generator_minus(i, j) = generator_minus(i, j) - perturbation
774 11140 : generator_minus(j, i) = generator_minus(j, i) + CONJG(perturbation)
775 11140 : s = s + 1
776 11140 : CALL qs_ot_rotation_state_local(generator_plus, rotation_work, gemm_ctx)
777 11140 : CPASSERT(rotation_work%info == 0)
778 : CALL qs_ot_dense_rotation_gradient(base_hamiltonian, occupation_scale, &
779 11140 : gradient_plus_real, gradient_plus_imag, rayleigh_plus, rotation_work, gemm_ctx)
780 11140 : CALL qs_ot_rotation_state_local(generator_minus, rotation_work, gemm_ctx)
781 11140 : CPASSERT(rotation_work%info == 0)
782 : CALL qs_ot_dense_rotation_gradient(base_hamiltonian, occupation_scale, &
783 11140 : gradient_minus_real, gradient_minus_imag, rayleigh_minus, rotation_work, gemm_ctx)
784 : CALL qs_ot_pack_rotation_response( &
785 : gradient_plus_real, gradient_plus_imag, gradient_minus_real, gradient_minus_imag, &
786 : rayleigh_plus, rayleigh_minus, step, rotation_hessian(:, s), &
787 16710 : rayleigh_response(:, s))
788 : END DO
789 : END DO
790 : END DO
791 362 : CPASSERT(s == nrotation)
792 1289452 : rotation_hessian(:, :) = 0.5_dp*(rotation_hessian + TRANSPOSE(rotation_hessian))
793 :
794 362 : CALL gemm_ctx%destroy()
795 0 : DEALLOCATE (base_hamiltonian, generator_minus, generator_plus, occupation_scale, &
796 0 : rayleigh, rayleigh_minus, rayleigh_plus, gradient_imag, gradient_minus_imag, &
797 362 : gradient_minus_real, gradient_plus_imag, gradient_plus_real, gradient_real)
798 :
799 724 : END SUBROUTINE qs_ot_finite_rotation_response
800 :
801 : ! **************************************************************************************************
802 : !> \brief pack one finite complex rotation response column
803 : !> \param gradient_plus_real real gradient at the positive endpoint
804 : !> \param gradient_plus_imag imaginary gradient at the positive endpoint
805 : !> \param gradient_minus_real real gradient at the negative endpoint
806 : !> \param gradient_minus_imag imaginary gradient at the negative endpoint
807 : !> \param rayleigh_plus Rayleigh energies at the positive endpoint
808 : !> \param rayleigh_minus Rayleigh energies at the negative endpoint
809 : !> \param step central finite-difference step
810 : !> \param hessian_column packed rotation-Hessian column
811 : !> \param rayleigh_column packed Rayleigh-response column
812 : ! **************************************************************************************************
813 11140 : SUBROUTINE qs_ot_pack_rotation_response( &
814 11140 : gradient_plus_real, gradient_plus_imag, gradient_minus_real, gradient_minus_imag, &
815 11140 : rayleigh_plus, rayleigh_minus, step, hessian_column, rayleigh_column)
816 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: gradient_plus_real, gradient_plus_imag, &
817 : gradient_minus_real, &
818 : gradient_minus_imag
819 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: rayleigh_plus, rayleigh_minus
820 : REAL(KIND=dp), INTENT(IN) :: step
821 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: hessian_column, rayleigh_column
822 :
823 : INTEGER :: i, j, n, r
824 :
825 11140 : n = SIZE(gradient_plus_real, 1)
826 33420 : CPASSERT(ALL(SHAPE(gradient_plus_real) == [n, n]))
827 33420 : CPASSERT(ALL(SHAPE(gradient_plus_imag) == [n, n]))
828 33420 : CPASSERT(ALL(SHAPE(gradient_minus_real) == [n, n]))
829 33420 : CPASSERT(ALL(SHAPE(gradient_minus_imag) == [n, n]))
830 11140 : CPASSERT(SIZE(hessian_column) == n*(n - 1))
831 11140 : CPASSERT(SIZE(rayleigh_column) == n)
832 :
833 11140 : r = 0
834 85568 : DO i = 1, n - 1
835 402180 : DO j = i + 1, n
836 316612 : r = r + 1
837 : hessian_column(r) = &
838 316612 : (gradient_plus_real(i, j) - gradient_minus_real(i, j))/(2.0_dp*step)
839 316612 : r = r + 1
840 : hessian_column(r) = &
841 391040 : (gradient_plus_imag(i, j) - gradient_minus_imag(i, j))/(2.0_dp*step)
842 : END DO
843 : END DO
844 96708 : rayleigh_column(:) = (rayleigh_plus(:) - rayleigh_minus(:))/(2.0_dp*step)
845 :
846 11140 : END SUBROUTINE qs_ot_pack_rotation_response
847 :
848 : ! **************************************************************************************************
849 : !> \brief Rotation gradient and projected-Hamiltonian diagonal from an already prepared rotation.
850 : !> \param base_hamiltonian fixed Hamiltonian in the unrotated REF basis
851 : !> \param occupation_scale twice the weighted occupations
852 : !> \param gradient_real real-antisymmetric gradient component
853 : !> \param gradient_imag imaginary-symmetric gradient component
854 : !> \param rayleigh diagonal of the rotated Hamiltonian
855 : !> \param w rotation state and private scratch
856 : !> \param gemm_ctx caller-owned GEMM context
857 : ! **************************************************************************************************
858 22642 : SUBROUTINE qs_ot_dense_rotation_gradient( &
859 22642 : base_hamiltonian, occupation_scale, gradient_real, gradient_imag, rayleigh, w, gemm_ctx)
860 : COMPLEX(KIND=dp), CONTIGUOUS, INTENT(IN) :: base_hamiltonian(:, :)
861 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: occupation_scale
862 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: gradient_real, gradient_imag
863 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: rayleigh
864 : TYPE(qs_ot_rotation_work_type), INTENT(INOUT) :: w
865 : TYPE(local_gemm_ctxt_type), INTENT(INOUT) :: gemm_ctx
866 :
867 : INTEGER :: i, j, n
868 :
869 22642 : n = w%solver%n
870 : CALL gemm_ctx%gemm('N', 'N', n, n, n, (1.0_dp, 0.0_dp), base_hamiltonian, n, w%u, n, &
871 22642 : (0.0_dp, 0.0_dp), w%work, n)
872 : CALL gemm_ctx%gemm('C', 'N', n, n, n, (1.0_dp, 0.0_dp), w%u, n, w%work, n, &
873 22642 : (0.0_dp, 0.0_dp), w%inner, n)
874 195786 : DO i = 1, n
875 195786 : rayleigh(i) = REAL(w%inner(i, i), KIND=dp)
876 : END DO
877 :
878 : ! Reuse H*U for the covector instead of multiplying a second time.
879 195786 : DO j = 1, n
880 1646518 : w%work(:, j) = occupation_scale(j)*w%work(:, j)
881 : END DO
882 22642 : CALL qs_ot_rotation_pullback_local(w, gemm_ctx, gradient_real, gradient_imag)
883 :
884 22642 : END SUBROUTINE qs_ot_dense_rotation_gradient
885 :
886 : ! **************************************************************************************************
887 : !> \brief Pull back dE/dU through exp(A) and project onto anti-Hermitian generator coordinates.
888 : !> \param w prepared rotation spectrum; work contains dE/dU on entry and is overwritten
889 : !> \param gemm_ctx caller-owned GEMM context
890 : !> \param gradient_real real-antisymmetric gradient component
891 : !> \param gradient_imag imaginary-symmetric gradient component
892 : ! **************************************************************************************************
893 23163 : SUBROUTINE qs_ot_rotation_pullback_local(w, gemm_ctx, gradient_real, gradient_imag)
894 : TYPE(qs_ot_rotation_work_type), INTENT(INOUT) :: w
895 : TYPE(local_gemm_ctxt_type), INTENT(INOUT) :: gemm_ctx
896 : REAL(KIND=dp), INTENT(OUT) :: gradient_real(:, :), gradient_imag(:, :)
897 :
898 : INTEGER :: i, j, n
899 :
900 23163 : n = w%solver%n
901 : CALL gemm_ctx%gemm('N', 'N', n, n, n, (1.0_dp, 0.0_dp), w%work, n, w%vectors, n, &
902 23163 : (0.0_dp, 0.0_dp), w%inner, n)
903 : CALL gemm_ctx%gemm('C', 'N', n, n, n, (1.0_dp, 0.0_dp), w%vectors, n, w%inner, n, &
904 23163 : (0.0_dp, 0.0_dp), w%work, n)
905 201605 : DO j = 1, n
906 1726581 : DO i = 1, n
907 : w%work(i, j) = w%work(i, j)* &
908 1703418 : CONJG(qs_ot_complex_exp_frechet_kernel(w%solver%evals(i), w%solver%evals(j)))
909 : END DO
910 : END DO
911 : CALL gemm_ctx%gemm('N', 'C', n, n, n, (1.0_dp, 0.0_dp), w%work, n, w%vectors, n, &
912 23163 : (0.0_dp, 0.0_dp), w%inner, n)
913 : CALL gemm_ctx%gemm('N', 'N', n, n, n, (1.0_dp, 0.0_dp), w%vectors, n, w%inner, n, &
914 23163 : (0.0_dp, 0.0_dp), w%work, n)
915 201605 : DO j = 1, n
916 1726581 : DO i = 1, n
917 1524976 : gradient_real(i, j) = REAL(w%work(i, j), KIND=dp) - REAL(w%work(j, i), KIND=dp)
918 1703418 : gradient_imag(i, j) = AIMAG(w%work(i, j)) + AIMAG(w%work(j, i))
919 : END DO
920 : END DO
921 23163 : END SUBROUTINE qs_ot_rotation_pullback_local
922 :
923 : ! **************************************************************************************************
924 : !> \brief Prepare private rotation arrays and LAPACK scratch once per matrix order.
925 : !> \param w caller-owned rotation workspace
926 : !> \param n rotation dimension
927 : ! **************************************************************************************************
928 2330 : SUBROUTINE qs_ot_prepare_rotation_work(w, n)
929 : TYPE(qs_ot_rotation_work_type), INTENT(INOUT) :: w
930 : INTEGER, INTENT(IN) :: n
931 :
932 2330 : CPASSERT(n > 0)
933 2330 : IF (ALLOCATED(w%vectors)) THEN
934 394 : IF (SIZE(w%vectors, 1) == n) RETURN
935 0 : DEALLOCATE (w%vectors, w%u, w%work, w%inner)
936 : END IF
937 19360 : ALLOCATE (w%vectors(n, n), w%u(n, n), w%work(n, n), w%inner(n, n))
938 1936 : CALL cp_cfm_local_workspace_prepare(w%solver, n, w%vectors)
939 : END SUBROUTINE qs_ot_prepare_rotation_work
940 :
941 : ! **************************************************************************************************
942 : !> \brief Evaluate exp(A) and retain its spectrum without allocating or communicating.
943 : !> \param rotation_generator anti-Hermitian generator
944 : !> \param w prepared workspace; u receives exp(A), info reports LAPACK failure
945 : !> \param gemm_ctx caller-owned GEMM context
946 : ! **************************************************************************************************
947 25602 : SUBROUTINE qs_ot_rotation_state_local(rotation_generator, w, gemm_ctx)
948 : COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(IN) :: rotation_generator
949 : TYPE(qs_ot_rotation_work_type), INTENT(INOUT) :: w
950 : TYPE(local_gemm_ctxt_type), INTENT(INOUT) :: gemm_ctx
951 :
952 : INTEGER :: j, n
953 :
954 25602 : n = w%solver%n
955 1951530 : w%vectors(:, :) = CMPLX(0.0_dp, 1.0_dp, KIND=dp)*rotation_generator
956 25602 : CALL cp_cfm_heevd_local(w%vectors, info=w%info, ws=w%solver)
957 25602 : IF (w%info /= 0) RETURN
958 222992 : DO j = 1, n
959 1951530 : w%work(:, j) = EXP(CMPLX(0.0_dp, -w%solver%evals(j), KIND=dp))*w%vectors(:, j)
960 : END DO
961 : CALL gemm_ctx%gemm('N', 'C', n, n, n, (1.0_dp, 0.0_dp), w%work, n, w%vectors, n, &
962 25602 : (0.0_dp, 0.0_dp), w%u, n)
963 : END SUBROUTINE qs_ot_rotation_state_local
964 :
965 : ! **************************************************************************************************
966 : !> \brief Frechet divided-difference kernel for exp(-i*evals)
967 : !> \param e1 ...
968 : !> \param e2 ...
969 : !> \return ...
970 : ! **************************************************************************************************
971 1569500 : PURE FUNCTION qs_ot_complex_exp_frechet_kernel(e1, e2) RESULT(kernel)
972 : REAL(KIND=dp), INTENT(IN) :: e1, e2
973 : COMPLEX(KIND=dp) :: kernel
974 :
975 : COMPLEX(KIND=dp) :: l1, l2, x
976 : INTEGER :: i
977 :
978 1569500 : l1 = (0.0_dp, -1.0_dp)*e1
979 1569500 : l2 = (0.0_dp, -1.0_dp)*e2
980 1569500 : IF (ABS(l1 - l2) > 0.5_dp) THEN
981 12696 : kernel = (EXP(l1) - EXP(l2))/(l1 - l2)
982 : ELSE
983 : x = 1.0_dp
984 : kernel = 0.0_dp
985 26465668 : DO i = 1, 16
986 24908864 : kernel = kernel + x
987 26465668 : x = x*(l1 - l2)/REAL(i + 1, KIND=dp)
988 : END DO
989 1556804 : kernel = kernel*EXP(l2)
990 : END IF
991 :
992 1569500 : END FUNCTION qs_ot_complex_exp_frechet_kernel
993 :
994 : ! **************************************************************************************************
995 : !> \brief apply the complex exponential Frechet kernel to sparse DBCSR Re/Im matrices
996 : !> \param evals generator eigenvalues
997 : !> \param inner_deriv_re real part of the matrix in the generator eigenbasis
998 : !> \param inner_deriv_im imaginary part of the matrix in the generator eigenbasis
999 : !> \param outer_deriv_re real part of the mapped matrix
1000 : !> \param outer_deriv_im imaginary part of the mapped matrix
1001 : !> \param adjoint use the adjoint Frechet kernel for gradients
1002 : ! **************************************************************************************************
1003 1098 : SUBROUTINE qs_ot_apply_complex_frechet_dbcsr(evals, inner_deriv_re, inner_deriv_im, &
1004 : outer_deriv_re, outer_deriv_im, adjoint)
1005 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: evals
1006 : TYPE(dbcsr_type) :: inner_deriv_re, inner_deriv_im, &
1007 : outer_deriv_re, outer_deriv_im
1008 : LOGICAL, INTENT(IN), OPTIONAL :: adjoint
1009 :
1010 : COMPLEX(KIND=dp) :: cval_in, kernel
1011 : INTEGER :: col, i, j, max_blocks, nblocks, row
1012 1098 : INTEGER, ALLOCATABLE, DIMENSION(:) :: cols, rows
1013 1098 : INTEGER, DIMENSION(:), POINTER :: col_blk_offset, col_blk_size, &
1014 1098 : row_blk_offset, row_blk_size
1015 : LOGICAL :: found_im, found_out_im, found_out_re, &
1016 : found_re, use_adjoint
1017 1098 : REAL(dp), DIMENSION(:, :), POINTER :: block_in_im, block_in_re, block_out_im, &
1018 1098 : block_out_re
1019 : REAL(KIND=dp) :: e1, e2, im_part, re_part
1020 : TYPE(dbcsr_distribution_type) :: dist
1021 : TYPE(dbcsr_iterator_type) :: iter
1022 :
1023 1098 : use_adjoint = .FALSE.
1024 1096 : IF (PRESENT(adjoint)) use_adjoint = adjoint
1025 :
1026 : ! Re/Im parts can have different sparse block patterns. Build their union
1027 : ! explicitly so a missing partner block is interpreted as zero.
1028 1098 : max_blocks = 0
1029 1098 : CALL dbcsr_iterator_start(iter, inner_deriv_re)
1030 1650 : DO WHILE (dbcsr_iterator_blocks_left(iter))
1031 552 : CALL dbcsr_iterator_next_block(iter, row, col)
1032 552 : max_blocks = max_blocks + 1
1033 : END DO
1034 1098 : CALL dbcsr_iterator_stop(iter)
1035 1098 : CALL dbcsr_iterator_start(iter, inner_deriv_im)
1036 1650 : DO WHILE (dbcsr_iterator_blocks_left(iter))
1037 552 : CALL dbcsr_iterator_next_block(iter, row, col)
1038 552 : max_blocks = max_blocks + 1
1039 : END DO
1040 1098 : CALL dbcsr_iterator_stop(iter)
1041 4392 : ALLOCATE (rows(MAX(max_blocks, 1)), cols(MAX(max_blocks, 1)))
1042 1098 : nblocks = 0
1043 :
1044 1098 : CALL dbcsr_iterator_start(iter, inner_deriv_re)
1045 1650 : DO WHILE (dbcsr_iterator_blocks_left(iter))
1046 552 : CALL dbcsr_iterator_next_block(iter, row, col)
1047 552 : CALL append_union_block(row, col)
1048 : END DO
1049 1098 : CALL dbcsr_iterator_stop(iter)
1050 1098 : CALL dbcsr_iterator_start(iter, inner_deriv_im)
1051 1650 : DO WHILE (dbcsr_iterator_blocks_left(iter))
1052 552 : CALL dbcsr_iterator_next_block(iter, row, col)
1053 552 : CALL append_union_block(row, col)
1054 : END DO
1055 1098 : CALL dbcsr_iterator_stop(iter)
1056 :
1057 : CALL dbcsr_get_info(inner_deriv_re, distribution=dist, row_blk_size=row_blk_size, &
1058 1098 : col_blk_size=col_blk_size)
1059 : CALL dbcsr_create(outer_deriv_re, "outer_deriv_re", dist, dbcsr_type_no_symmetry, &
1060 1098 : row_blk_size, col_blk_size)
1061 : CALL dbcsr_create(outer_deriv_im, "outer_deriv_im", dist, dbcsr_type_no_symmetry, &
1062 1098 : row_blk_size, col_blk_size)
1063 1098 : IF (nblocks > 0) THEN
1064 552 : CALL dbcsr_reserve_blocks(outer_deriv_re, rows=rows(1:nblocks), cols=cols(1:nblocks))
1065 552 : CALL dbcsr_reserve_blocks(outer_deriv_im, rows=rows(1:nblocks), cols=cols(1:nblocks))
1066 : END IF
1067 1098 : CALL dbcsr_finalize(outer_deriv_re)
1068 1098 : CALL dbcsr_finalize(outer_deriv_im)
1069 1098 : CALL dbcsr_set(outer_deriv_re, 0.0_dp)
1070 1098 : CALL dbcsr_set(outer_deriv_im, 0.0_dp)
1071 :
1072 : CALL dbcsr_get_info(outer_deriv_re, row_blk_offset=row_blk_offset, &
1073 1098 : col_blk_offset=col_blk_offset)
1074 1098 : CALL dbcsr_iterator_start(iter, outer_deriv_re)
1075 1652 : DO WHILE (dbcsr_iterator_blocks_left(iter))
1076 554 : CALL dbcsr_iterator_next_block(iter, row, col)
1077 554 : CALL dbcsr_get_block_p(inner_deriv_re, row, col, block_in_re, found_re)
1078 554 : CALL dbcsr_get_block_p(inner_deriv_im, row, col, block_in_im, found_im)
1079 554 : CALL dbcsr_get_block_p(outer_deriv_re, row, col, block_out_re, found_out_re)
1080 554 : CALL dbcsr_get_block_p(outer_deriv_im, row, col, block_out_im, found_out_im)
1081 554 : CPASSERT(found_out_re .AND. found_out_im)
1082 :
1083 5558 : DO i = 1, SIZE(block_out_re, 1)
1084 48980 : DO j = 1, SIZE(block_out_re, 2)
1085 44520 : e1 = evals(row_blk_offset(row) + i - 1)
1086 44520 : e2 = evals(col_blk_offset(col) + j - 1)
1087 44520 : re_part = 0.0_dp
1088 44520 : im_part = 0.0_dp
1089 44520 : IF (found_re) re_part = block_in_re(i, j)
1090 44520 : IF (found_im) im_part = block_in_im(i, j)
1091 44520 : cval_in = CMPLX(re_part, im_part, dp)
1092 44520 : kernel = qs_ot_complex_exp_frechet_kernel(e1, e2)
1093 44520 : IF (use_adjoint) kernel = CONJG(kernel)
1094 44520 : cval_in = cval_in*kernel
1095 44520 : block_out_re(i, j) = REAL(cval_in, KIND=dp)
1096 48426 : block_out_im(i, j) = AIMAG(cval_in)
1097 : END DO
1098 : END DO
1099 : END DO
1100 1098 : CALL dbcsr_iterator_stop(iter)
1101 8782 : DEALLOCATE (rows, cols)
1102 :
1103 : CONTAINS
1104 :
1105 : ! **************************************************************************************************
1106 : !> \brief append a block coordinate unless it is already present
1107 : !> \param row_new block-row index
1108 : !> \param col_new block-column index
1109 : ! **************************************************************************************************
1110 1104 : SUBROUTINE append_union_block(row_new, col_new)
1111 : INTEGER, INTENT(IN) :: row_new, col_new
1112 :
1113 : INTEGER :: iblock
1114 :
1115 1106 : DO iblock = 1, nblocks
1116 1106 : IF (rows(iblock) == row_new .AND. cols(iblock) == col_new) RETURN
1117 : END DO
1118 554 : nblocks = nblocks + 1
1119 554 : rows(nblocks) = row_new
1120 554 : cols(nblocks) = col_new
1121 : END SUBROUTINE append_union_block
1122 :
1123 : END SUBROUTINE qs_ot_apply_complex_frechet_dbcsr
1124 :
1125 : ! **************************************************************************************************
1126 : !> \brief update complex REF k-point orbitals and their S(k)C(k) images
1127 : !> \param matrix_c ...
1128 : !> \param matrix_c_im ...
1129 : !> \param matrix_s ...
1130 : !> \param matrix_s_im ...
1131 : !> \param qs_ot_env ...
1132 : !> \param qs_ot_env1 environment carrying the shared minimizer state
1133 : ! **************************************************************************************************
1134 1842 : SUBROUTINE qs_ot_get_orbitals_ref_complex(matrix_c, matrix_c_im, matrix_s, matrix_s_im, &
1135 : qs_ot_env, qs_ot_env1)
1136 :
1137 : TYPE(dbcsr_type), POINTER :: matrix_c, matrix_c_im, matrix_s, &
1138 : matrix_s_im
1139 : TYPE(qs_ot_type) :: qs_ot_env
1140 : TYPE(qs_ot_type), OPTIONAL :: qs_ot_env1
1141 :
1142 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_ot_get_orbitals_ref_complex'
1143 :
1144 : INTEGER :: handle, i, k, n
1145 : LOGICAL :: update
1146 1842 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eigenvalues, inverse_sqrt
1147 : TYPE(dbcsr_type) :: rotated_im, rotated_re, rotation_tmp
1148 : TYPE(dbcsr_type), POINTER :: c_im, c_re, f_im, f_re, p_im, p_re, &
1149 : sc_im, sc_re, tmp_kk, tmp_nk, v_im, &
1150 : v_re, w_im, w_re
1151 :
1152 1842 : CALL timeset(routineN, handle)
1153 :
1154 1842 : CPASSERT(qs_ot_env%has_complex_kpoint_state)
1155 1842 : CPASSERT(ASSOCIATED(matrix_s))
1156 1842 : CPASSERT(ASSOCIATED(matrix_s_im))
1157 :
1158 1842 : c_re => qs_ot_env%matrix_x
1159 1842 : c_im => qs_ot_env%matrix_x_im
1160 1842 : f_re => qs_ot_env%matrix_ref_inv_sqrt
1161 1842 : f_im => qs_ot_env%matrix_ref_inv_sqrt_im
1162 1842 : sc_re => qs_ot_env%matrix_sx
1163 1842 : sc_im => qs_ot_env%matrix_sx_im
1164 1842 : p_re => qs_ot_env%buf1_k_k_sym
1165 1842 : p_im => qs_ot_env%buf2_k_k_sym
1166 1842 : v_re => qs_ot_env%buf3_k_k_sym
1167 1842 : v_im => qs_ot_env%buf4_k_k_sym
1168 1842 : w_re => qs_ot_env%buf1_k_k_nosym
1169 1842 : w_im => qs_ot_env%buf2_k_k_nosym
1170 1842 : tmp_kk => qs_ot_env%buf3_k_k_nosym
1171 1842 : tmp_nk => qs_ot_env%buf1_n_k
1172 :
1173 1842 : CALL dbcsr_get_info(c_re, nfullrows_total=n, nfullcols_total=k)
1174 1842 : IF (PRESENT(qs_ot_env1)) THEN
1175 1724 : update = qs_ot_ref_update(qs_ot_env1)
1176 : ELSE
1177 118 : update = qs_ot_ref_update(qs_ot_env)
1178 : END IF
1179 :
1180 : ! SC = (S_re + i*S_im) * (C_re + i*C_im)
1181 1842 : CALL qs_ot_complex_multiply('N', 'N', matrix_s, matrix_s_im, c_re, c_im, sc_re, sc_im, tmp_nk)
1182 :
1183 : ! P = C^H*S*C. Its imaginary component is real antisymmetric.
1184 1842 : CALL dbcsr_multiply('T', 'N', 1.0_dp, c_re, sc_re, 0.0_dp, p_re)
1185 1842 : CALL dbcsr_multiply('T', 'N', 1.0_dp, c_im, sc_im, 0.0_dp, tmp_kk)
1186 1842 : CALL dbcsr_add(p_re, tmp_kk, alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
1187 :
1188 1842 : CALL dbcsr_multiply('T', 'N', 1.0_dp, c_re, sc_im, 0.0_dp, p_im)
1189 1842 : CALL dbcsr_multiply('T', 'N', 1.0_dp, c_im, sc_re, 0.0_dp, tmp_kk)
1190 1842 : CALL dbcsr_add(p_im, tmp_kk, alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
1191 :
1192 7368 : ALLOCATE (eigenvalues(k), inverse_sqrt(k))
1193 : CALL cp_dbcsr_heevd(matrix_re=p_re, matrix_im=p_im, &
1194 : eigenvectors_re=v_re, eigenvectors_im=v_im, &
1195 : eigenvalues=eigenvalues, para_env=qs_ot_env%para_env, &
1196 1842 : blacs_env=qs_ot_env%blacs_env)
1197 10858 : DO i = 1, k
1198 9016 : IF (eigenvalues(i) <= EPSILON(1.0_dp)) THEN
1199 0 : CPABORT("Complex REF overlap is not positive definite")
1200 : END IF
1201 10858 : inverse_sqrt(i) = 1.0_dp/SQRT(eigenvalues(i))
1202 : END DO
1203 :
1204 : ! P^(-1/2) = V*diag(lambda^(-1/2))*V^H.
1205 1842 : CALL dbcsr_copy(w_re, v_re)
1206 1842 : CALL dbcsr_copy(w_im, v_im)
1207 1842 : CALL dbcsr_scale_by_vector(w_re, alpha=inverse_sqrt, side='right')
1208 1842 : CALL dbcsr_scale_by_vector(w_im, alpha=inverse_sqrt, side='right')
1209 :
1210 1842 : CALL dbcsr_multiply('N', 'T', 1.0_dp, w_re, v_re, 0.0_dp, p_re)
1211 1842 : CALL dbcsr_multiply('N', 'T', 1.0_dp, w_im, v_im, 0.0_dp, tmp_kk)
1212 1842 : CALL dbcsr_add(p_re, tmp_kk, alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
1213 :
1214 1842 : CALL dbcsr_multiply('N', 'T', 1.0_dp, w_im, v_re, 0.0_dp, p_im)
1215 1842 : CALL dbcsr_multiply('N', 'T', 1.0_dp, w_re, v_im, 0.0_dp, tmp_kk)
1216 1842 : CALL dbcsr_add(p_im, tmp_kk, alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
1217 :
1218 1842 : CALL dbcsr_copy(f_re, p_re)
1219 1842 : CALL dbcsr_copy(f_im, p_im)
1220 :
1221 : ! Return the physical, orthonormal orbitals without changing a rejected
1222 : ! line-search coordinate.
1223 1842 : CALL qs_ot_complex_multiply('N', 'N', c_re, c_im, p_re, p_im, matrix_c, matrix_c_im, tmp_nk)
1224 :
1225 1842 : IF (update) THEN
1226 610 : CALL dbcsr_copy(c_re, matrix_c)
1227 610 : CALL dbcsr_copy(c_im, matrix_c_im)
1228 :
1229 610 : CALL dbcsr_set(f_re, 0.0_dp)
1230 610 : CALL dbcsr_add_on_diag(f_re, alpha=1.0_dp)
1231 610 : CALL dbcsr_set(f_im, 0.0_dp)
1232 :
1233 610 : CALL qs_ot_complex_multiply('N', 'N', matrix_s, matrix_s_im, c_re, c_im, sc_re, sc_im, tmp_nk)
1234 : END IF
1235 :
1236 1842 : IF (qs_ot_env%settings%do_rotation) THEN
1237 1272 : CALL qs_ot_generate_rotation_complex(qs_ot_env)
1238 :
1239 1272 : CALL dbcsr_copy(rotated_re, matrix_c, name="rotated_re")
1240 1272 : CALL dbcsr_copy(rotated_im, matrix_c_im, name="rotated_im")
1241 1272 : CALL dbcsr_copy(rotation_tmp, matrix_c, name="rotation_tmp")
1242 :
1243 : ! C_out = Q*U for complex Q and unitary U.
1244 : CALL qs_ot_complex_multiply('N', 'N', matrix_c, matrix_c_im, qs_ot_env%rot_mat_u, &
1245 1272 : qs_ot_env%rot_mat_u_im, rotated_re, rotated_im, rotation_tmp)
1246 :
1247 1272 : CALL dbcsr_copy(matrix_c, rotated_re)
1248 1272 : CALL dbcsr_copy(matrix_c_im, rotated_im)
1249 1272 : CALL dbcsr_release(rotated_re)
1250 1272 : CALL dbcsr_release(rotated_im)
1251 1272 : CALL dbcsr_release(rotation_tmp)
1252 : END IF
1253 :
1254 1842 : DEALLOCATE (eigenvalues, inverse_sqrt)
1255 :
1256 1842 : CALL timestop(handle)
1257 3684 : END SUBROUTINE qs_ot_get_orbitals_ref_complex
1258 :
1259 : ! **************************************************************************************************
1260 : !> \brief complex k-point REF derivative dE/dX from H(k)C(k), S(k)C(k), and C(k)
1261 : !> \param matrix_hc ...
1262 : !> \param matrix_hc_im ...
1263 : !> \param qs_ot_env ...
1264 : ! **************************************************************************************************
1265 2168 : SUBROUTINE qs_ot_get_derivative_ref_complex(matrix_hc, matrix_hc_im, qs_ot_env)
1266 : TYPE(dbcsr_type), POINTER :: matrix_hc, matrix_hc_im
1267 : TYPE(qs_ot_type) :: qs_ot_env
1268 :
1269 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_ot_get_derivative_ref_complex'
1270 :
1271 : INTEGER :: handle, k, n
1272 : TYPE(dbcsr_type) :: tmp_nk
1273 : TYPE(dbcsr_type), POINTER :: b_im, b_re, f_im, f_re, g_im, g_re, &
1274 : hc_im, hc_re, hc_work_im, hc_work_re, &
1275 : q_im, q_re, sc_im, sc_re, tmp_kk
1276 : TYPE(dbcsr_type), TARGET :: hc_rot_im, hc_rot_re
1277 :
1278 1084 : CALL timeset(routineN, handle)
1279 :
1280 1084 : CPASSERT(qs_ot_env%has_complex_kpoint_state)
1281 1084 : CPASSERT(ASSOCIATED(matrix_hc))
1282 1084 : CPASSERT(ASSOCIATED(matrix_hc_im))
1283 :
1284 1084 : f_re => qs_ot_env%matrix_ref_inv_sqrt
1285 1084 : f_im => qs_ot_env%matrix_ref_inv_sqrt_im
1286 1084 : sc_re => qs_ot_env%matrix_sx
1287 1084 : sc_im => qs_ot_env%matrix_sx_im
1288 1084 : hc_re => matrix_hc
1289 1084 : hc_im => matrix_hc_im
1290 1084 : hc_work_re => hc_re
1291 1084 : hc_work_im => hc_im
1292 1084 : g_re => qs_ot_env%matrix_gx
1293 1084 : g_im => qs_ot_env%matrix_gx_im
1294 1084 : b_re => qs_ot_env%buf1_k_k_sym
1295 1084 : b_im => qs_ot_env%buf2_k_k_sym
1296 1084 : tmp_kk => qs_ot_env%buf3_k_k_sym
1297 1084 : q_re => qs_ot_env%buf1_n_k
1298 1084 : q_im => qs_ot_env%buf1_n_k_dp
1299 :
1300 1084 : CALL dbcsr_get_info(sc_re, nfullrows_total=n, nfullcols_total=k)
1301 :
1302 : ! Q = X*(X^H*S*X)^(-1/2), reconstructed from the current REF coordinate.
1303 : CALL qs_ot_complex_multiply('N', 'N', qs_ot_env%matrix_x, qs_ot_env%matrix_x_im, f_re, f_im, q_re, q_im, &
1304 1084 : g_re)
1305 :
1306 1084 : IF (qs_ot_env%settings%do_rotation) THEN
1307 832 : CALL qs_ot_generate_rotation_complex(qs_ot_env)
1308 :
1309 : ! dF/dU = Q^H*G_C for C=Q*U.
1310 : CALL dbcsr_multiply('T', 'N', 1.0_dp, q_re, hc_re, &
1311 832 : 0.0_dp, qs_ot_env%rot_mat_dedu)
1312 832 : CALL dbcsr_multiply('T', 'N', 1.0_dp, q_im, hc_im, 0.0_dp, tmp_kk)
1313 : CALL dbcsr_add(qs_ot_env%rot_mat_dedu, tmp_kk, &
1314 832 : alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
1315 :
1316 : CALL dbcsr_multiply('T', 'N', 1.0_dp, q_re, hc_im, &
1317 832 : 0.0_dp, qs_ot_env%rot_mat_dedu_im)
1318 832 : CALL dbcsr_multiply('T', 'N', 1.0_dp, q_im, hc_re, 0.0_dp, tmp_kk)
1319 : CALL dbcsr_add(qs_ot_env%rot_mat_dedu_im, tmp_kk, &
1320 832 : alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
1321 832 : CALL qs_ot_rot_mat_derivative_complex(qs_ot_env)
1322 :
1323 : ! The REF coordinate sees G_Q=G_C*U^H.
1324 832 : CALL dbcsr_copy(hc_rot_re, hc_re, name="hc_rot_re")
1325 832 : CALL dbcsr_copy(hc_rot_im, hc_im, name="hc_rot_im")
1326 832 : CALL dbcsr_copy(tmp_nk, hc_re, name="tmp_nk")
1327 : CALL dbcsr_multiply('N', 'T', 1.0_dp, hc_re, qs_ot_env%rot_mat_u, &
1328 832 : 0.0_dp, hc_rot_re)
1329 : CALL dbcsr_multiply('N', 'T', 1.0_dp, hc_im, qs_ot_env%rot_mat_u_im, &
1330 832 : 0.0_dp, tmp_nk)
1331 832 : CALL dbcsr_add(hc_rot_re, tmp_nk, alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
1332 :
1333 : CALL dbcsr_multiply('N', 'T', 1.0_dp, hc_im, qs_ot_env%rot_mat_u, &
1334 832 : 0.0_dp, hc_rot_im)
1335 : CALL dbcsr_multiply('N', 'T', 1.0_dp, hc_re, qs_ot_env%rot_mat_u_im, &
1336 832 : 0.0_dp, tmp_nk)
1337 832 : CALL dbcsr_add(hc_rot_im, tmp_nk, alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
1338 832 : hc_work_re => hc_rot_re
1339 832 : hc_work_im => hc_rot_im
1340 : END IF
1341 :
1342 : ! B = Q^H*G_Q. For uniform fixed occupations B is Hermitian.
1343 1084 : CALL dbcsr_multiply('T', 'N', 1.0_dp, q_re, hc_work_re, 0.0_dp, b_re)
1344 1084 : CALL dbcsr_multiply('T', 'N', 1.0_dp, q_im, hc_work_im, 0.0_dp, tmp_kk)
1345 1084 : CALL dbcsr_add(b_re, tmp_kk, alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
1346 :
1347 1084 : CALL dbcsr_multiply('T', 'N', 1.0_dp, q_re, hc_work_im, 0.0_dp, b_im)
1348 1084 : CALL dbcsr_multiply('T', 'N', 1.0_dp, q_im, hc_work_re, 0.0_dp, tmp_kk)
1349 1084 : CALL dbcsr_add(b_im, tmp_kk, alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
1350 :
1351 1084 : IF (qs_ot_env%settings%eps_irac_filter_matrix > 0.0_dp) THEN
1352 0 : CALL dbcsr_filter(b_re, qs_ot_env%settings%eps_irac_filter_matrix)
1353 0 : CALL dbcsr_filter(b_im, qs_ot_env%settings%eps_irac_filter_matrix)
1354 : END IF
1355 :
1356 : ! S*Q = (S*X)*F. G is temporary storage for this pair.
1357 1084 : CALL qs_ot_complex_multiply('N', 'N', sc_re, sc_im, f_re, f_im, g_re, g_im, q_re)
1358 :
1359 : ! Form (S*Q)*B. The Q workspaces are free after B has been formed.
1360 1084 : CALL dbcsr_multiply('N', 'N', 1.0_dp, g_re, b_re, 0.0_dp, q_re)
1361 1084 : CALL dbcsr_multiply('N', 'N', 1.0_dp, g_im, b_im, 0.0_dp, q_im)
1362 1084 : CALL dbcsr_add(q_re, q_im, alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
1363 :
1364 1084 : CALL dbcsr_multiply('N', 'N', 1.0_dp, g_re, b_im, 0.0_dp, q_im)
1365 1084 : CALL dbcsr_multiply('N', 'N', 1.0_dp, g_im, b_re, 0.0_dp, g_re)
1366 1084 : CALL dbcsr_add(q_im, g_re, alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
1367 :
1368 1084 : CALL dbcsr_add(q_re, hc_work_re, alpha_scalar=-1.0_dp, beta_scalar=1.0_dp)
1369 1084 : CALL dbcsr_add(q_im, hc_work_im, alpha_scalar=-1.0_dp, beta_scalar=1.0_dp)
1370 :
1371 : ! Pull the physical gradient back to the finite REF coordinate: G_X = (G_C-S*Q*B)*F.
1372 1084 : CALL dbcsr_multiply('N', 'N', 1.0_dp, q_re, f_re, 0.0_dp, g_re)
1373 1084 : CALL dbcsr_multiply('N', 'N', 1.0_dp, q_im, f_im, 0.0_dp, g_im)
1374 1084 : CALL dbcsr_add(g_re, g_im, alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
1375 :
1376 1084 : CALL dbcsr_multiply('N', 'N', 1.0_dp, q_re, f_im, 0.0_dp, g_im)
1377 1084 : CALL dbcsr_multiply('N', 'N', 1.0_dp, q_im, f_re, 0.0_dp, q_re)
1378 1084 : CALL dbcsr_add(g_im, q_re, alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
1379 :
1380 1084 : IF (qs_ot_env%settings%do_rotation) THEN
1381 : CALL qs_ot_add_ref_vertical_response_complex(b_re, b_im, f_re, f_im, &
1382 832 : sc_re, sc_im, g_re, g_im, qs_ot_env)
1383 : END IF
1384 :
1385 1084 : IF (qs_ot_env%settings%eps_irac_filter_matrix > 0.0_dp) THEN
1386 0 : CALL dbcsr_filter(g_re, qs_ot_env%settings%eps_irac_filter_matrix)
1387 0 : CALL dbcsr_filter(g_im, qs_ot_env%settings%eps_irac_filter_matrix)
1388 : END IF
1389 :
1390 1084 : IF (qs_ot_env%settings%do_rotation) THEN
1391 832 : CALL dbcsr_release(hc_rot_re)
1392 832 : CALL dbcsr_release(hc_rot_im)
1393 832 : CALL dbcsr_release(tmp_nk)
1394 : END IF
1395 :
1396 1084 : CALL timestop(handle)
1397 :
1398 1084 : END SUBROUTINE qs_ot_get_derivative_ref_complex
1399 :
1400 : ! **************************************************************************************************
1401 : !> \brief Add the vertical response of the finite complex polar REF map.
1402 : !> \param a_re real part of Q^H G_Q
1403 : !> \param a_im imaginary part of Q^H G_Q
1404 : !> \param f_re real part of (X^H S X)^(-1/2)
1405 : !> \param f_im imaginary part of (X^H S X)^(-1/2)
1406 : !> \param sx_re real part of S X
1407 : !> \param sx_im imaginary part of S X
1408 : !> \param gradient_re real REF gradient, updated in place
1409 : !> \param gradient_im imaginary REF gradient, updated in place
1410 : !> \param qs_ot_env complex OT environment
1411 : ! **************************************************************************************************
1412 832 : SUBROUTINE qs_ot_add_ref_vertical_response_complex(a_re, a_im, f_re, f_im, &
1413 : sx_re, sx_im, gradient_re, gradient_im, &
1414 : qs_ot_env)
1415 : TYPE(dbcsr_type) :: a_re, a_im, f_re, f_im, sx_re, sx_im, &
1416 : gradient_re, gradient_im
1417 : TYPE(qs_ot_type) :: qs_ot_env
1418 :
1419 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_ot_add_ref_vertical_response_complex'
1420 :
1421 : INTEGER :: col, handle, i, j, k, row
1422 832 : INTEGER, DIMENSION(:), POINTER :: col_blk_offset, row_blk_offset
1423 : LOGICAL :: found_im, found_re
1424 : REAL(KIND=dp) :: denominator
1425 832 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: f_evals
1426 832 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: block_im, block_re
1427 : TYPE(dbcsr_iterator_type) :: iter
1428 : TYPE(dbcsr_type) :: anti_im, anti_re, f_work_im, f_work_re, inner_im, inner_re, sq_im, &
1429 : sq_re, tmp_kk, tmp_nk, v_im, v_re, vertical_im, vertical_re, work_im, work_re, z_im, z_re
1430 :
1431 832 : CALL timeset(routineN, handle)
1432 :
1433 832 : CALL dbcsr_get_info(a_re, nfullrows_total=k)
1434 832 : IF (k == 0) THEN
1435 0 : CALL timestop(handle)
1436 0 : RETURN
1437 : END IF
1438 :
1439 : ! Only the anti-Hermitian part of A=Q^H G_Q couples to the unitary
1440 : ! component of the polar differential.
1441 832 : CALL qs_ot_square_transpose(a_re, tmp_kk, qs_ot_env%rot_mat_u)
1442 832 : CALL dbcsr_copy(anti_re, a_re, name='anti_re')
1443 832 : CALL dbcsr_add(anti_re, tmp_kk, alpha_scalar=0.5_dp, beta_scalar=-0.5_dp)
1444 832 : CALL qs_ot_square_transpose(a_im, tmp_kk, qs_ot_env%rot_mat_u)
1445 832 : CALL dbcsr_copy(anti_im, a_im, name='anti_im')
1446 832 : CALL dbcsr_add(anti_im, tmp_kk, alpha_scalar=0.5_dp, beta_scalar=0.5_dp)
1447 :
1448 : ! P=(X^H S X)^(1/2)=F^(-1). In the eigenbasis of F, solve
1449 : ! P Z + Z P = antiherm(A) element by element.
1450 832 : CALL dbcsr_copy(f_work_re, f_re, name='f_work_re')
1451 832 : CALL dbcsr_copy(f_work_im, f_im, name='f_work_im')
1452 832 : CALL dbcsr_copy(v_re, f_re, name='v_re')
1453 832 : CALL dbcsr_copy(v_im, f_im, name='v_im')
1454 2496 : ALLOCATE (f_evals(k))
1455 : CALL cp_dbcsr_heevd(matrix_re=f_work_re, matrix_im=f_work_im, &
1456 : eigenvectors_re=v_re, eigenvectors_im=v_im, eigenvalues=f_evals, &
1457 832 : para_env=qs_ot_env%para_env, blacs_env=qs_ot_env%blacs_env)
1458 5276 : IF (MINVAL(f_evals) <= EPSILON(1.0_dp)) THEN
1459 0 : CPABORT('Complex REF inverse square root is not positive definite')
1460 : END IF
1461 :
1462 : ! inner = V^H antiherm(A) V.
1463 832 : CALL dbcsr_copy(work_re, anti_re, name='work_re')
1464 832 : CALL dbcsr_copy(work_im, anti_im, name='work_im')
1465 832 : CALL qs_ot_complex_multiply('N', 'N', anti_re, anti_im, v_re, v_im, work_re, work_im, tmp_kk)
1466 :
1467 832 : CALL dbcsr_copy(inner_re, anti_re, name='inner_re')
1468 832 : CALL dbcsr_copy(inner_im, anti_im, name='inner_im')
1469 832 : CALL dbcsr_multiply('T', 'N', 1.0_dp, v_re, work_re, 0.0_dp, inner_re)
1470 832 : CALL dbcsr_multiply('T', 'N', 1.0_dp, v_im, work_im, 1.0_dp, inner_re)
1471 832 : CALL dbcsr_multiply('T', 'N', 1.0_dp, v_re, work_im, 0.0_dp, inner_im)
1472 832 : CALL dbcsr_multiply('T', 'N', -1.0_dp, v_im, work_re, 1.0_dp, inner_im)
1473 :
1474 : CALL dbcsr_get_info(inner_re, row_blk_offset=row_blk_offset, &
1475 832 : col_blk_offset=col_blk_offset)
1476 832 : CALL dbcsr_iterator_start(iter, inner_re)
1477 1249 : DO WHILE (dbcsr_iterator_blocks_left(iter))
1478 417 : CALL dbcsr_iterator_next_block(iter, row, col)
1479 417 : CALL dbcsr_get_block_p(inner_re, row, col, block_re, found_re)
1480 417 : CALL dbcsr_get_block_p(inner_im, row, col, block_im, found_im)
1481 417 : CPASSERT(found_re .AND. found_im)
1482 3473 : DO i = 1, SIZE(block_re, 1)
1483 21793 : DO j = 1, SIZE(block_re, 2)
1484 : denominator = 1.0_dp/f_evals(row_blk_offset(row) + i - 1) + &
1485 19152 : 1.0_dp/f_evals(col_blk_offset(col) + j - 1)
1486 19152 : block_re(i, j) = block_re(i, j)/denominator
1487 21376 : block_im(i, j) = block_im(i, j)/denominator
1488 : END DO
1489 : END DO
1490 : END DO
1491 832 : CALL dbcsr_iterator_stop(iter)
1492 :
1493 : ! Z = V inner V^H.
1494 832 : CALL qs_ot_complex_multiply('N', 'N', v_re, v_im, inner_re, inner_im, work_re, work_im, tmp_kk)
1495 :
1496 832 : CALL dbcsr_copy(z_re, anti_re, name='z_re')
1497 832 : CALL dbcsr_copy(z_im, anti_im, name='z_im')
1498 832 : CALL dbcsr_multiply('N', 'T', 1.0_dp, work_re, v_re, 0.0_dp, z_re)
1499 832 : CALL dbcsr_multiply('N', 'T', 1.0_dp, work_im, v_im, 1.0_dp, z_re)
1500 832 : CALL dbcsr_multiply('N', 'T', 1.0_dp, work_im, v_re, 0.0_dp, z_im)
1501 832 : CALL dbcsr_multiply('N', 'T', -1.0_dp, work_re, v_im, 1.0_dp, z_im)
1502 :
1503 : ! The adjoint vertical contribution is 2 S Q Z, with S Q=(S X)F.
1504 832 : CALL dbcsr_copy(sq_re, sx_re, name='sq_re')
1505 832 : CALL dbcsr_copy(sq_im, sx_im, name='sq_im')
1506 832 : CALL dbcsr_copy(tmp_nk, sx_re, name='tmp_nk')
1507 832 : CALL qs_ot_complex_multiply('N', 'N', sx_re, sx_im, f_re, f_im, sq_re, sq_im, tmp_nk)
1508 :
1509 832 : CALL dbcsr_copy(vertical_re, sx_re, name='vertical_re')
1510 832 : CALL dbcsr_copy(vertical_im, sx_im, name='vertical_im')
1511 832 : CALL qs_ot_complex_multiply('N', 'N', sq_re, sq_im, z_re, z_im, vertical_re, vertical_im, tmp_nk)
1512 832 : CALL dbcsr_add(gradient_re, vertical_re, alpha_scalar=1.0_dp, beta_scalar=2.0_dp)
1513 832 : CALL dbcsr_add(gradient_im, vertical_im, alpha_scalar=1.0_dp, beta_scalar=2.0_dp)
1514 :
1515 832 : DEALLOCATE (f_evals)
1516 832 : CALL dbcsr_release(anti_im)
1517 832 : CALL dbcsr_release(anti_re)
1518 832 : CALL dbcsr_release(f_work_im)
1519 832 : CALL dbcsr_release(f_work_re)
1520 832 : CALL dbcsr_release(inner_im)
1521 832 : CALL dbcsr_release(inner_re)
1522 832 : CALL dbcsr_release(sq_im)
1523 832 : CALL dbcsr_release(sq_re)
1524 832 : CALL dbcsr_release(tmp_kk)
1525 832 : CALL dbcsr_release(tmp_nk)
1526 832 : CALL dbcsr_release(v_im)
1527 832 : CALL dbcsr_release(v_re)
1528 832 : CALL dbcsr_release(vertical_im)
1529 832 : CALL dbcsr_release(vertical_re)
1530 832 : CALL dbcsr_release(work_im)
1531 832 : CALL dbcsr_release(work_re)
1532 832 : CALL dbcsr_release(z_im)
1533 832 : CALL dbcsr_release(z_re)
1534 :
1535 832 : CALL timestop(handle)
1536 :
1537 1664 : END SUBROUTINE qs_ot_add_ref_vertical_response_complex
1538 :
1539 : ! **************************************************************************************************
1540 : !> \brief computes U=exp(A) for the complex anti-Hermitian generator
1541 : !> A=rot_mat_x+i*rot_mat_x_im
1542 : !> \param qs_ot_env a complex k-point OT environment
1543 : ! **************************************************************************************************
1544 3350 : SUBROUTINE qs_ot_generate_rotation_complex(qs_ot_env)
1545 :
1546 : TYPE(qs_ot_type) :: qs_ot_env
1547 :
1548 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_ot_generate_rotation_complex'
1549 :
1550 : INTEGER :: handle, k
1551 : REAL(KIND=dp) :: rot_norm
1552 3350 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: exp_evals_im, exp_evals_re
1553 : TYPE(dbcsr_type) :: h_re, tmp, w_im, w_re
1554 :
1555 3350 : CALL timeset(routineN, handle)
1556 :
1557 3350 : CPASSERT(qs_ot_env%has_complex_kpoint_state)
1558 3350 : CPASSERT(ASSOCIATED(qs_ot_env%rot_mat_x_im))
1559 3350 : CPASSERT(ASSOCIATED(qs_ot_env%rot_mat_u_im))
1560 :
1561 3350 : CALL dbcsr_get_info(qs_ot_env%rot_mat_x, nfullrows_total=k)
1562 3350 : IF (k /= 0) THEN
1563 : rot_norm = dbcsr_frobenius_norm(qs_ot_env%rot_mat_x) + &
1564 3350 : dbcsr_frobenius_norm(qs_ot_env%rot_mat_x_im)
1565 3350 : IF (rot_norm <= EPSILON(1.0_dp)) THEN
1566 744 : CALL dbcsr_set(qs_ot_env%rot_mat_u, 0.0_dp)
1567 744 : CALL dbcsr_add_on_diag(qs_ot_env%rot_mat_u, 1.0_dp)
1568 744 : CALL dbcsr_set(qs_ot_env%rot_mat_u_im, 0.0_dp)
1569 744 : CALL timestop(handle)
1570 : RETURN
1571 : END IF
1572 :
1573 : ! i*A = i*X-Y is Hermitian. Its eigenvectors give
1574 : ! exp(A)=V*diag(exp(-i*lambda))*V^H.
1575 2606 : CALL dbcsr_copy(h_re, qs_ot_env%rot_mat_x_im, name="h_re")
1576 2606 : CALL dbcsr_scale(h_re, -1.0_dp)
1577 : CALL cp_dbcsr_heevd(matrix_re=h_re, matrix_im=qs_ot_env%rot_mat_x, &
1578 : eigenvectors_re=qs_ot_env%rot_mat_evec_re, &
1579 : eigenvectors_im=qs_ot_env%rot_mat_evec_im, &
1580 : eigenvalues=qs_ot_env%rot_mat_evals, &
1581 2606 : para_env=qs_ot_env%para_env, blacs_env=qs_ot_env%blacs_env)
1582 :
1583 10424 : ALLOCATE (exp_evals_re(k), exp_evals_im(k))
1584 20554 : exp_evals_re(:) = COS(-qs_ot_env%rot_mat_evals(:))
1585 20554 : exp_evals_im(:) = SIN(-qs_ot_env%rot_mat_evals(:))
1586 :
1587 2606 : CALL dbcsr_copy(w_re, qs_ot_env%rot_mat_evec_re, name="w_re")
1588 2606 : CALL dbcsr_scale_by_vector(w_re, alpha=exp_evals_re, side='right')
1589 2606 : CALL dbcsr_copy(tmp, qs_ot_env%rot_mat_evec_im, name="tmp")
1590 2606 : CALL dbcsr_scale_by_vector(tmp, alpha=exp_evals_im, side='right')
1591 2606 : CALL dbcsr_add(w_re, tmp, alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
1592 :
1593 2606 : CALL dbcsr_copy(w_im, qs_ot_env%rot_mat_evec_im, name="w_im")
1594 2606 : CALL dbcsr_scale_by_vector(w_im, alpha=exp_evals_re, side='right')
1595 2606 : CALL dbcsr_copy(tmp, qs_ot_env%rot_mat_evec_re)
1596 2606 : CALL dbcsr_scale_by_vector(tmp, alpha=exp_evals_im, side='right')
1597 2606 : CALL dbcsr_add(w_im, tmp, alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
1598 :
1599 : CALL dbcsr_multiply('N', 'T', 1.0_dp, w_re, qs_ot_env%rot_mat_evec_re, &
1600 2606 : 0.0_dp, qs_ot_env%rot_mat_u)
1601 : CALL dbcsr_multiply('N', 'T', 1.0_dp, w_im, qs_ot_env%rot_mat_evec_im, &
1602 2606 : 1.0_dp, qs_ot_env%rot_mat_u)
1603 : CALL dbcsr_multiply('N', 'T', 1.0_dp, w_im, qs_ot_env%rot_mat_evec_re, &
1604 2606 : 0.0_dp, qs_ot_env%rot_mat_u_im)
1605 : CALL dbcsr_multiply('N', 'T', -1.0_dp, w_re, qs_ot_env%rot_mat_evec_im, &
1606 2606 : 1.0_dp, qs_ot_env%rot_mat_u_im)
1607 :
1608 2606 : CALL dbcsr_release(h_re)
1609 2606 : CALL dbcsr_release(tmp)
1610 2606 : CALL dbcsr_release(w_re)
1611 2606 : CALL dbcsr_release(w_im)
1612 2606 : DEALLOCATE (exp_evals_re, exp_evals_im)
1613 : END IF
1614 :
1615 2606 : CALL timestop(handle)
1616 :
1617 5956 : END SUBROUTINE qs_ot_generate_rotation_complex
1618 :
1619 : ! **************************************************************************************************
1620 : !> \brief pull the complex dE/dU covector back to the anti-Hermitian generator
1621 : !> using the adjoint Frechet derivative of exp
1622 : !> \param qs_ot_env a complex k-point OT environment with an up-to-date U
1623 : ! **************************************************************************************************
1624 2932 : SUBROUTINE qs_ot_rot_mat_derivative_complex(qs_ot_env)
1625 : TYPE(qs_ot_type) :: qs_ot_env
1626 :
1627 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_ot_rot_mat_derivative_complex'
1628 :
1629 : INTEGER :: handle, k
1630 : REAL(KIND=dp) :: rot_norm
1631 : TYPE(dbcsr_type) :: frechet_im, frechet_re, inner_deriv_im, &
1632 : inner_deriv_re, outer_deriv_im, &
1633 : outer_deriv_re, tmp, work_im, work_re
1634 :
1635 1466 : CALL timeset(routineN, handle)
1636 :
1637 1466 : CPASSERT(qs_ot_env%has_complex_kpoint_state)
1638 1466 : CPASSERT(ASSOCIATED(qs_ot_env%rot_mat_dedu_im))
1639 1466 : CPASSERT(ASSOCIATED(qs_ot_env%rot_mat_gx_im))
1640 :
1641 1466 : CALL dbcsr_get_info(qs_ot_env%rot_mat_u, nfullrows_total=k)
1642 1466 : IF (k /= 0) THEN
1643 : rot_norm = dbcsr_frobenius_norm(qs_ot_env%rot_mat_x) + &
1644 1466 : dbcsr_frobenius_norm(qs_ot_env%rot_mat_x_im)
1645 1466 : IF (rot_norm <= EPSILON(1.0_dp)) THEN
1646 370 : CALL qs_ot_square_transpose(qs_ot_env%rot_mat_dedu, tmp, qs_ot_env%rot_mat_u)
1647 370 : CALL dbcsr_copy(qs_ot_env%rot_mat_gx, qs_ot_env%rot_mat_dedu)
1648 : CALL dbcsr_add(qs_ot_env%rot_mat_gx, tmp, &
1649 370 : alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
1650 370 : CALL dbcsr_release(tmp)
1651 :
1652 370 : CALL qs_ot_square_transpose(qs_ot_env%rot_mat_dedu_im, tmp, qs_ot_env%rot_mat_u)
1653 370 : CALL dbcsr_copy(qs_ot_env%rot_mat_gx_im, qs_ot_env%rot_mat_dedu_im)
1654 : CALL dbcsr_add(qs_ot_env%rot_mat_gx_im, tmp, &
1655 370 : alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
1656 370 : CALL dbcsr_release(tmp)
1657 370 : CALL timestop(handle)
1658 370 : RETURN
1659 : END IF
1660 :
1661 1096 : CALL dbcsr_copy(work_re, qs_ot_env%rot_mat_dedu, name="work_re")
1662 1096 : CALL dbcsr_copy(work_im, qs_ot_env%rot_mat_dedu, name="work_im")
1663 1096 : CALL dbcsr_copy(tmp, qs_ot_env%rot_mat_dedu, name="tmp")
1664 1096 : CALL dbcsr_copy(inner_deriv_re, qs_ot_env%rot_mat_dedu, name="inner_deriv_re")
1665 1096 : CALL dbcsr_copy(inner_deriv_im, qs_ot_env%rot_mat_dedu, name="inner_deriv_im")
1666 :
1667 : ! V^H*(dE/dU)*V, split into real and imaginary parts.
1668 : CALL qs_ot_complex_multiply('N', 'N', qs_ot_env%rot_mat_dedu, qs_ot_env%rot_mat_dedu_im, &
1669 : qs_ot_env%rot_mat_evec_re, qs_ot_env%rot_mat_evec_im, work_re, work_im, &
1670 1096 : tmp)
1671 : CALL dbcsr_multiply('T', 'N', 1.0_dp, qs_ot_env%rot_mat_evec_re, &
1672 1096 : work_re, 0.0_dp, inner_deriv_re)
1673 : CALL dbcsr_multiply('T', 'N', 1.0_dp, qs_ot_env%rot_mat_evec_im, &
1674 1096 : work_im, 1.0_dp, inner_deriv_re)
1675 : CALL dbcsr_multiply('T', 'N', 1.0_dp, qs_ot_env%rot_mat_evec_re, &
1676 1096 : work_im, 0.0_dp, inner_deriv_im)
1677 : CALL dbcsr_multiply('T', 'N', -1.0_dp, qs_ot_env%rot_mat_evec_im, &
1678 1096 : work_re, 1.0_dp, inner_deriv_im)
1679 :
1680 : CALL qs_ot_apply_complex_frechet_dbcsr(qs_ot_env%rot_mat_evals, &
1681 : inner_deriv_re, inner_deriv_im, &
1682 1096 : outer_deriv_re, outer_deriv_im, adjoint=.TRUE.)
1683 :
1684 : CALL qs_ot_complex_multiply('N', 'N', qs_ot_env%rot_mat_evec_re, qs_ot_env%rot_mat_evec_im, &
1685 1096 : outer_deriv_re, outer_deriv_im, work_re, work_im, tmp)
1686 :
1687 1096 : CALL dbcsr_copy(frechet_re, qs_ot_env%rot_mat_dedu, name="frechet_re")
1688 1096 : CALL dbcsr_copy(frechet_im, qs_ot_env%rot_mat_dedu, name="frechet_im")
1689 : CALL dbcsr_multiply('N', 'T', 1.0_dp, work_re, qs_ot_env%rot_mat_evec_re, &
1690 1096 : 0.0_dp, frechet_re)
1691 : CALL dbcsr_multiply('N', 'T', 1.0_dp, work_im, qs_ot_env%rot_mat_evec_im, &
1692 1096 : 1.0_dp, frechet_re)
1693 : CALL dbcsr_multiply('N', 'T', 1.0_dp, work_im, qs_ot_env%rot_mat_evec_re, &
1694 1096 : 0.0_dp, frechet_im)
1695 : CALL dbcsr_multiply('N', 'T', -1.0_dp, work_re, qs_ot_env%rot_mat_evec_im, &
1696 1096 : 1.0_dp, frechet_im)
1697 :
1698 : ! Tangents satisfy X^T=-X and Y^T=Y for A=X+iY.
1699 1096 : CALL qs_ot_square_transpose(frechet_re, tmp, qs_ot_env%rot_mat_u)
1700 1096 : CALL dbcsr_copy(qs_ot_env%rot_mat_gx, frechet_re)
1701 : CALL dbcsr_add(qs_ot_env%rot_mat_gx, tmp, &
1702 1096 : alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
1703 1096 : CALL qs_ot_square_transpose(frechet_im, tmp, qs_ot_env%rot_mat_u)
1704 1096 : CALL dbcsr_copy(qs_ot_env%rot_mat_gx_im, frechet_im)
1705 : CALL dbcsr_add(qs_ot_env%rot_mat_gx_im, tmp, &
1706 1096 : alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
1707 :
1708 1096 : CALL dbcsr_release(frechet_re)
1709 1096 : CALL dbcsr_release(frechet_im)
1710 1096 : CALL dbcsr_release(inner_deriv_re)
1711 1096 : CALL dbcsr_release(inner_deriv_im)
1712 1096 : CALL dbcsr_release(outer_deriv_re)
1713 1096 : CALL dbcsr_release(outer_deriv_im)
1714 1096 : CALL dbcsr_release(tmp)
1715 1096 : CALL dbcsr_release(work_re)
1716 1096 : CALL dbcsr_release(work_im)
1717 : END IF
1718 :
1719 1096 : CALL timestop(handle)
1720 :
1721 1466 : END SUBROUTINE qs_ot_rot_mat_derivative_complex
1722 :
1723 : ! **************************************************************************************************
1724 : !> \brief compute P=X^H*S*X and the STRICT matrix functions for a complex K-point channel
1725 : !> \param matrix_x real part of X
1726 : !> \param matrix_x_im imaginary part of X
1727 : !> \param matrix_sx real part of S*X
1728 : !> \param matrix_sx_im imaginary part of S*X
1729 : !> \param qs_ot_env OT channel state
1730 : !> \param zero_coordinates initialize the known X=0 state without diagonalization
1731 : ! **************************************************************************************************
1732 1982 : SUBROUTINE qs_ot_get_p_complex(matrix_x, matrix_x_im, matrix_sx, matrix_sx_im, qs_ot_env, zero_coordinates)
1733 : TYPE(dbcsr_type), POINTER :: matrix_x, matrix_x_im, matrix_sx, &
1734 : matrix_sx_im
1735 : TYPE(qs_ot_type) :: qs_ot_env
1736 : LOGICAL, INTENT(IN), OPTIONAL :: zero_coordinates
1737 :
1738 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_ot_get_p_complex'
1739 :
1740 : INTEGER :: handle
1741 :
1742 1982 : CALL timeset(routineN, handle)
1743 :
1744 1982 : CPASSERT(qs_ot_env%has_complex_kpoint_state)
1745 1982 : CPASSERT(qs_ot_env%settings%ot_algorithm == 'TOD')
1746 :
1747 1982 : IF (PRESENT(zero_coordinates)) THEN
1748 408 : IF (zero_coordinates) THEN
1749 408 : CALL dbcsr_set(qs_ot_env%matrix_p, 0.0_dp)
1750 408 : CALL dbcsr_set(qs_ot_env%matrix_p_im, 0.0_dp)
1751 408 : CALL dbcsr_set(qs_ot_env%matrix_r, 0.0_dp)
1752 408 : CALL dbcsr_add_on_diag(qs_ot_env%matrix_r, 1.0_dp)
1753 408 : CALL dbcsr_set(qs_ot_env%matrix_r_im, 0.0_dp)
1754 408 : CALL dbcsr_copy(qs_ot_env%matrix_cosp, qs_ot_env%matrix_r)
1755 408 : CALL dbcsr_copy(qs_ot_env%matrix_sinp, qs_ot_env%matrix_r)
1756 408 : CALL dbcsr_set(qs_ot_env%matrix_cosp_im, 0.0_dp)
1757 408 : CALL dbcsr_set(qs_ot_env%matrix_sinp_im, 0.0_dp)
1758 408 : CALL dbcsr_reserve_all_blocks(qs_ot_env%matrix_cosp_b)
1759 408 : CALL dbcsr_reserve_all_blocks(qs_ot_env%matrix_sinp_b)
1760 408 : CALL dbcsr_set(qs_ot_env%matrix_cosp_b, -0.5_dp)
1761 408 : CALL dbcsr_set(qs_ot_env%matrix_sinp_b, -1.0_dp/6.0_dp)
1762 2192 : qs_ot_env%evals(:) = 0.0_dp
1763 2192 : qs_ot_env%dum(:) = 0.0_dp
1764 408 : qs_ot_env%do_taylor = .FALSE.
1765 408 : CALL timestop(handle)
1766 408 : RETURN
1767 : END IF
1768 : END IF
1769 :
1770 : CALL qs_ot_complex_multiply('C', 'N', matrix_x, matrix_x_im, matrix_sx, matrix_sx_im, &
1771 1574 : qs_ot_env%matrix_p, qs_ot_env%matrix_p_im, qs_ot_env%matrix_buf1)
1772 1574 : qs_ot_env%do_taylor = .FALSE.
1773 1574 : CALL qs_ot_p2m_diag(qs_ot_env)
1774 :
1775 1574 : CALL timestop(handle)
1776 :
1777 : END SUBROUTINE qs_ot_get_p_complex
1778 :
1779 : ! **************************************************************************************************
1780 : !> \brief update complex K-point orbitals with the finite STRICT transformation
1781 : !> \param matrix_c real output orbitals
1782 : !> \param matrix_c_im imaginary output orbitals
1783 : !> \param matrix_s real overlap matrix
1784 : !> \param matrix_s_im imaginary overlap matrix
1785 : !> \param qs_ot_env OT channel state
1786 : ! **************************************************************************************************
1787 1572 : SUBROUTINE qs_ot_get_orbitals_complex(matrix_c, matrix_c_im, matrix_s, matrix_s_im, qs_ot_env)
1788 : TYPE(dbcsr_type), POINTER :: matrix_c, matrix_c_im, matrix_s, &
1789 : matrix_s_im
1790 : TYPE(qs_ot_type) :: qs_ot_env
1791 :
1792 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_ot_get_orbitals_complex'
1793 :
1794 : INTEGER :: handle
1795 : TYPE(dbcsr_type) :: rotated_im, rotated_re, rotation_tmp
1796 :
1797 1572 : CALL timeset(routineN, handle)
1798 :
1799 1572 : CPASSERT(qs_ot_env%has_complex_kpoint_state)
1800 1572 : CPASSERT(qs_ot_env%settings%ot_algorithm == 'TOD')
1801 :
1802 : CALL qs_ot_complex_multiply('N', 'N', matrix_s, matrix_s_im, &
1803 : qs_ot_env%matrix_x, qs_ot_env%matrix_x_im, &
1804 : qs_ot_env%matrix_sx, qs_ot_env%matrix_sx_im, &
1805 1572 : qs_ot_env%matrix_tmp_nk)
1806 : CALL qs_ot_get_p_complex(qs_ot_env%matrix_x, qs_ot_env%matrix_x_im, &
1807 1572 : qs_ot_env%matrix_sx, qs_ot_env%matrix_sx_im, qs_ot_env)
1808 :
1809 : CALL qs_ot_complex_multiply('N', 'N', qs_ot_env%matrix_c0, qs_ot_env%matrix_c0_im, &
1810 : qs_ot_env%matrix_cosp, qs_ot_env%matrix_cosp_im, &
1811 1572 : matrix_c, matrix_c_im, qs_ot_env%matrix_tmp_nk)
1812 : CALL qs_ot_complex_multiply('N', 'N', qs_ot_env%matrix_x, qs_ot_env%matrix_x_im, &
1813 : qs_ot_env%matrix_sinp, qs_ot_env%matrix_sinp_im, &
1814 : qs_ot_env%matrix_buf_nk, qs_ot_env%matrix_buf_nk_im, &
1815 1572 : qs_ot_env%matrix_tmp_nk)
1816 1572 : CALL dbcsr_add(matrix_c, qs_ot_env%matrix_buf_nk, alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
1817 1572 : CALL dbcsr_add(matrix_c_im, qs_ot_env%matrix_buf_nk_im, alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
1818 :
1819 1572 : IF (qs_ot_env%settings%do_rotation) THEN
1820 612 : CALL qs_ot_generate_rotation_complex(qs_ot_env)
1821 :
1822 612 : CALL dbcsr_copy(rotated_re, matrix_c, name="strict_rotated_re")
1823 612 : CALL dbcsr_copy(rotated_im, matrix_c_im, name="strict_rotated_im")
1824 612 : CALL dbcsr_copy(rotation_tmp, matrix_c, name="strict_rotation_tmp")
1825 :
1826 : ! C_out = Q(X)*U for the finite STRICT chart Q(X).
1827 : CALL qs_ot_complex_multiply('N', 'N', matrix_c, matrix_c_im, qs_ot_env%rot_mat_u, &
1828 612 : qs_ot_env%rot_mat_u_im, rotated_re, rotated_im, rotation_tmp)
1829 :
1830 612 : CALL dbcsr_copy(matrix_c, rotated_re)
1831 612 : CALL dbcsr_copy(matrix_c_im, rotated_im)
1832 612 : CALL dbcsr_release(rotated_re)
1833 612 : CALL dbcsr_release(rotated_im)
1834 612 : CALL dbcsr_release(rotation_tmp)
1835 : END IF
1836 :
1837 1572 : CALL timestop(handle)
1838 :
1839 1572 : END SUBROUTINE qs_ot_get_orbitals_complex
1840 :
1841 : ! **************************************************************************************************
1842 : !> \brief Project a complex STRICT gradient; real minimizer fallbacks use the shared real projection.
1843 : !> \param qs_ot_env OT channel state
1844 : !> \param matrix_gx real component to project in place
1845 : !> \param matrix_gx_im imaginary component for a complex channel
1846 : !> \param use_preconditioner use the attached preconditioner, default true
1847 : ! **************************************************************************************************
1848 1306 : SUBROUTINE qs_ot_project_gradient(qs_ot_env, matrix_gx, matrix_gx_im, use_preconditioner)
1849 : TYPE(qs_ot_type) :: qs_ot_env
1850 : TYPE(dbcsr_type), POINTER :: matrix_gx
1851 : TYPE(dbcsr_type), OPTIONAL, POINTER :: matrix_gx_im
1852 : LOGICAL, INTENT(IN), OPTIONAL :: use_preconditioner
1853 :
1854 : LOGICAL :: preconditioned
1855 : TYPE(dbcsr_type), POINTER :: target_im, target_re
1856 :
1857 1306 : IF (.NOT. qs_ot_env%has_complex_kpoint_state) THEN
1858 0 : CALL qs_ot_project_gradient_real(qs_ot_env, matrix_gx, use_preconditioner)
1859 0 : RETURN
1860 : END IF
1861 1306 : preconditioned = .TRUE.
1862 1306 : IF (PRESENT(use_preconditioner)) preconditioned = use_preconditioner
1863 1306 : IF (.NOT. preconditioned) qs_ot_env%os_valid = .FALSE.
1864 1306 : CPASSERT(PRESENT(matrix_gx_im))
1865 1306 : CALL qs_ot_prepare_complex_tangent_metric(qs_ot_env, use_preconditioner=preconditioned)
1866 1306 : IF (preconditioned .AND. ASSOCIATED(qs_ot_env%preconditioner)) THEN
1867 1272 : target_re => qs_ot_env%matrix_psc0
1868 1272 : target_im => qs_ot_env%matrix_psc0_im
1869 : ELSE
1870 34 : target_re => qs_ot_env%matrix_sc0
1871 34 : target_im => qs_ot_env%matrix_sc0_im
1872 : END IF
1873 :
1874 : CALL qs_ot_complex_multiply('C', 'N', target_re, target_im, &
1875 : matrix_gx, matrix_gx_im, &
1876 : qs_ot_env%matrix_buf1_ortho, qs_ot_env%matrix_buf1_ortho_im, &
1877 1306 : qs_ot_env%matrix_tmp_ortho)
1878 : CALL qs_ot_complex_multiply('N', 'N', qs_ot_env%matrix_os, qs_ot_env%matrix_os_im, &
1879 : qs_ot_env%matrix_buf1_ortho, qs_ot_env%matrix_buf1_ortho_im, &
1880 : qs_ot_env%matrix_buf2_ortho, qs_ot_env%matrix_buf2_ortho_im, &
1881 1306 : qs_ot_env%matrix_tmp_ortho)
1882 : CALL qs_ot_complex_multiply('N', 'N', qs_ot_env%matrix_sc0, qs_ot_env%matrix_sc0_im, &
1883 : qs_ot_env%matrix_buf2_ortho, qs_ot_env%matrix_buf2_ortho_im, &
1884 : qs_ot_env%matrix_buf_nk, qs_ot_env%matrix_buf_nk_im, &
1885 1306 : qs_ot_env%matrix_tmp_nk)
1886 : CALL dbcsr_add(matrix_gx, qs_ot_env%matrix_buf_nk, &
1887 1306 : alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
1888 : CALL dbcsr_add(matrix_gx_im, qs_ot_env%matrix_buf_nk_im, &
1889 1306 : alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
1890 :
1891 : END SUBROUTINE qs_ot_project_gradient
1892 :
1893 : ! **************************************************************************************************
1894 : !> \brief Prepare the inverse metric used to project a complex STRICT gradient.
1895 : !> An unusable preconditioner is detached before any minimizer history is updated.
1896 : !> \param qs_ot_env OT channel state
1897 : !> \param preconditioner_rejected true if the attached preconditioner was not positive definite
1898 : !> \param use_preconditioner use the attached preconditioner, default true
1899 : ! **************************************************************************************************
1900 6630 : SUBROUTINE qs_ot_prepare_complex_tangent_metric(qs_ot_env, preconditioner_rejected, use_preconditioner)
1901 : TYPE(qs_ot_type) :: qs_ot_env
1902 : LOGICAL, INTENT(OUT), OPTIONAL :: preconditioner_rejected
1903 : LOGICAL, INTENT(IN), OPTIONAL :: use_preconditioner
1904 :
1905 : INTEGER :: i, k
1906 : LOGICAL :: preconditioned
1907 : REAL(KIND=dp) :: eval_scale, eval_threshold
1908 : TYPE(dbcsr_type), POINTER :: target_im, target_re
1909 :
1910 6454 : preconditioned = .TRUE.
1911 6454 : IF (PRESENT(use_preconditioner)) preconditioned = use_preconditioner
1912 1306 : IF (.NOT. preconditioned) qs_ot_env%os_valid = .FALSE.
1913 6454 : IF (PRESENT(preconditioner_rejected)) preconditioner_rejected = .FALSE.
1914 6454 : IF (qs_ot_env%os_valid) RETURN
1915 :
1916 176 : IF (preconditioned .AND. ASSOCIATED(qs_ot_env%preconditioner)) THEN
1917 166 : target_re => qs_ot_env%matrix_psc0
1918 166 : target_im => qs_ot_env%matrix_psc0_im
1919 : CALL apply_preconditioner(qs_ot_env%preconditioner, &
1920 : qs_ot_env%matrix_sc0, qs_ot_env%matrix_sc0_im, &
1921 166 : target_re, target_im)
1922 : ELSE
1923 10 : target_re => qs_ot_env%matrix_sc0
1924 10 : target_im => qs_ot_env%matrix_sc0_im
1925 : END IF
1926 : CALL qs_ot_complex_multiply('C', 'N', qs_ot_env%matrix_sc0, qs_ot_env%matrix_sc0_im, &
1927 : target_re, target_im, qs_ot_env%matrix_os, qs_ot_env%matrix_os_im, &
1928 176 : qs_ot_env%matrix_buf1)
1929 176 : CALL dbcsr_get_info(qs_ot_env%matrix_os, nfullrows_total=k)
1930 : ! Preserve the eigenvectors and eigenvalues of X^H*S*X used by the STRICT derivative.
1931 : CALL cp_dbcsr_heevd(matrix_re=qs_ot_env%matrix_os, matrix_im=qs_ot_env%matrix_os_im, &
1932 : eigenvectors_re=qs_ot_env%matrix_buf3, eigenvectors_im=qs_ot_env%matrix_buf3_im, &
1933 : eigenvalues=qs_ot_env%dum, para_env=qs_ot_env%para_env, &
1934 176 : blacs_env=qs_ot_env%blacs_env)
1935 1320 : eval_scale = MAX(1.0_dp, MAXVAL(ABS(qs_ot_env%dum(1:k))))
1936 176 : eval_threshold = 100.0_dp*EPSILON(1.0_dp)*eval_scale
1937 : IF (MINVAL(qs_ot_env%dum(1:k)) <= eval_threshold .AND. &
1938 1320 : preconditioned .AND. ASSOCIATED(qs_ot_env%preconditioner)) THEN
1939 0 : NULLIFY (qs_ot_env%preconditioner)
1940 0 : qs_ot_env%preconditioner_changed = .TRUE.
1941 0 : IF (PRESENT(preconditioner_rejected)) preconditioner_rejected = .TRUE.
1942 0 : target_re => qs_ot_env%matrix_sc0
1943 0 : target_im => qs_ot_env%matrix_sc0_im
1944 : CALL qs_ot_complex_multiply('C', 'N', qs_ot_env%matrix_sc0, qs_ot_env%matrix_sc0_im, &
1945 : target_re, target_im, qs_ot_env%matrix_os, qs_ot_env%matrix_os_im, &
1946 0 : qs_ot_env%matrix_buf1)
1947 : CALL cp_dbcsr_heevd(matrix_re=qs_ot_env%matrix_os, matrix_im=qs_ot_env%matrix_os_im, &
1948 : eigenvectors_re=qs_ot_env%matrix_buf3, eigenvectors_im=qs_ot_env%matrix_buf3_im, &
1949 : eigenvalues=qs_ot_env%dum, para_env=qs_ot_env%para_env, &
1950 0 : blacs_env=qs_ot_env%blacs_env)
1951 0 : eval_scale = MAX(1.0_dp, MAXVAL(ABS(qs_ot_env%dum(1:k))))
1952 0 : eval_threshold = 100.0_dp*EPSILON(1.0_dp)*eval_scale
1953 : END IF
1954 1320 : CPASSERT(MINVAL(qs_ot_env%dum(1:k)) > eval_threshold)
1955 1320 : DO i = 1, k
1956 1320 : qs_ot_env%dum(i) = 1.0_dp/qs_ot_env%dum(i)
1957 : END DO
1958 176 : CALL dbcsr_copy(qs_ot_env%matrix_buf1, qs_ot_env%matrix_buf3)
1959 176 : CALL dbcsr_copy(qs_ot_env%matrix_buf1_im, qs_ot_env%matrix_buf3_im)
1960 176 : CALL dbcsr_scale_by_vector(qs_ot_env%matrix_buf1, alpha=qs_ot_env%dum, side='right')
1961 176 : CALL dbcsr_scale_by_vector(qs_ot_env%matrix_buf1_im, alpha=qs_ot_env%dum, side='right')
1962 : CALL qs_ot_complex_multiply('N', 'C', qs_ot_env%matrix_buf1, qs_ot_env%matrix_buf1_im, &
1963 : qs_ot_env%matrix_buf3, qs_ot_env%matrix_buf3_im, &
1964 : qs_ot_env%matrix_os, qs_ot_env%matrix_os_im, &
1965 176 : qs_ot_env%matrix_buf2)
1966 : ! An identity-metric projection must not replace the cached preconditioned metric.
1967 176 : qs_ot_env%os_valid = preconditioned
1968 :
1969 : END SUBROUTINE qs_ot_prepare_complex_tangent_metric
1970 :
1971 : ! **************************************************************************************************
1972 : !> \brief Prepare uncached single-rank STRICT metrics in bounded OpenMP batches.
1973 : !> Leave unsuccessful channels to the scalar path so preconditioner rejection and history
1974 : !> invalidation retain their existing order. Workers use only private dense storage.
1975 : !> \param qs_ot_env local OT channels
1976 : ! **************************************************************************************************
1977 96868 : SUBROUTINE qs_ot_prepare_complex_tangent_metric_batch(qs_ot_env)
1978 : TYPE(qs_ot_type), INTENT(INOUT) :: qs_ot_env(:)
1979 :
1980 : INTEGER :: first, i, j, k, n, nbatch, nchannels, &
1981 : nworkers
1982 96868 : INTEGER, ALLOCATABLE :: channels(:)
1983 : LOGICAL :: handled
1984 : REAL(dp) :: eval_threshold
1985 96868 : TYPE(dbcsr_p_type), ALLOCATABLE :: matrix_in(:, :), matrix_out(:, :)
1986 96868 : TYPE(preconditioner_p_type), ALLOCATABLE :: preconditioners(:)
1987 96868 : TYPE(qs_ot_local_work_type), ALLOCATABLE :: work(:)
1988 :
1989 96868 : IF (.NOT. qs_ot_local_channels(qs_ot_env)) RETURN
1990 1734 : IF (.NOT. qs_ot_env(1)%settings%ks .OR. qs_ot_env(1)%settings%ot_algorithm /= 'TOD') RETURN
1991 792 : IF (diag_lib_explicit .OR. diag_type == FM_DIAG_TYPE_DLAF) RETURN
1992 2634 : nchannels = COUNT(.NOT. qs_ot_env(:)%os_valid)
1993 792 : IF (nchannels == 0) RETURN
1994 444 : ALLOCATE (channels(nchannels))
1995 148 : j = 0
1996 482 : DO i = 1, SIZE(qs_ot_env)
1997 334 : IF (qs_ot_env(i)%os_valid) CYCLE
1998 334 : j = j + 1
1999 482 : channels(j) = i
2000 : END DO
2001 148 : nworkers = 1
2002 148 : !$ nworkers = MIN(nchannels, omp_get_max_threads())
2003 0 : ALLOCATE (work(nworkers), preconditioners(nworkers), &
2004 2368 : matrix_in(2, nworkers), matrix_out(2, nworkers))
2005 296 : DO j = 1, nworkers
2006 296 : CALL work(j)%gemm_ctx%create(LOCAL_GEMM_PU_GPU, timing=.FALSE.)
2007 : END DO
2008 482 : DO first = 1, nchannels, nworkers
2009 334 : nbatch = MIN(nworkers, nchannels - first + 1)
2010 668 : DO j = 1, nbatch
2011 334 : i = channels(first + j - 1)
2012 334 : preconditioners(j)%preconditioner => qs_ot_env(i)%preconditioner
2013 334 : matrix_in(1, j)%matrix => qs_ot_env(i)%matrix_sc0
2014 334 : matrix_in(2, j)%matrix => qs_ot_env(i)%matrix_sc0_im
2015 1002 : matrix_out(:, j) = matrix_in(:, j)
2016 668 : IF (ASSOCIATED(qs_ot_env(i)%preconditioner)) THEN
2017 334 : matrix_out(1, j)%matrix => qs_ot_env(i)%matrix_psc0
2018 334 : matrix_out(2, j)%matrix => qs_ot_env(i)%matrix_psc0_im
2019 : END IF
2020 : END DO
2021 : CALL apply_preconditioner_complex_batch(preconditioners(1:nbatch), &
2022 334 : matrix_in(:, 1:nbatch), matrix_out(:, 1:nbatch), handled)
2023 668 : DO j = 1, nbatch
2024 334 : IF (.NOT. handled .AND. ASSOCIATED(preconditioners(j)%preconditioner)) THEN
2025 : CALL apply_preconditioner(preconditioners(j)%preconditioner, &
2026 : matrix_in(1, j)%matrix, matrix_in(2, j)%matrix, &
2027 0 : matrix_out(1, j)%matrix, matrix_out(2, j)%matrix)
2028 : END IF
2029 334 : CALL dbcsr_get_info(matrix_in(1, j)%matrix, nfullrows_total=n, nfullcols_total=k)
2030 334 : IF (ALLOCATED(work(j)%sc0)) THEN
2031 186 : IF (n /= work(j)%n .OR. k /= work(j)%k) THEN
2032 0 : DEALLOCATE (work(j)%sc0, work(j)%target, work(j)%os, work(j)%r, &
2033 0 : work(j)%buf1, work(j)%buf2, work(j)%vectors, work(j)%values)
2034 : END IF
2035 : END IF
2036 334 : IF (.NOT. ALLOCATED(work(j)%sc0)) THEN
2037 0 : ALLOCATE (work(j)%sc0(n, k, 2), work(j)%target(n, k, 2), work(j)%os(k, k, 2), &
2038 0 : work(j)%r(k, k, 2), work(j)%buf1(k, k, 2), work(j)%buf2(k, k, 2), &
2039 3256 : work(j)%vectors(k, k), work(j)%values(k))
2040 148 : work(j)%n = n
2041 148 : work(j)%k = k
2042 : END IF
2043 334 : CALL cp_cfm_local_workspace_prepare(work(j)%solver, k, work(j)%vectors)
2044 334 : CALL qs_ot_load_local_pair(matrix_in(1, j)%matrix, matrix_in(2, j)%matrix, work(j)%sc0)
2045 1002 : CALL qs_ot_load_local_pair(matrix_out(1, j)%matrix, matrix_out(2, j)%matrix, work(j)%target)
2046 : END DO
2047 : !$OMP PARALLEL DO DEFAULT(NONE) SCHEDULE(STATIC) NUM_THREADS(nbatch) IF(nbatch > 1) &
2048 334 : !$OMP SHARED(work,nbatch) PRIVATE(eval_threshold)
2049 : DO j = 1, nbatch
2050 : CALL qs_ot_complex_multiply_local('C', 'N', work(j)%sc0, work(j)%target, &
2051 : work(j)%os, work(j)%buf2(:, :, 1), work(j)%gemm_ctx)
2052 : work(j)%vectors(:, :) = CMPLX(work(j)%os(:, :, 1), work(j)%os(:, :, 2), KIND=dp)
2053 : CALL cp_cfm_heevd_local(work(j)%vectors, info=work(j)%info, ws=work(j)%solver)
2054 : IF (work(j)%info /= 0) CYCLE
2055 : eval_threshold = 100.0_dp*EPSILON(1.0_dp)*MAX(1.0_dp, MAXVAL(ABS(work(j)%solver%evals)))
2056 : IF (MINVAL(work(j)%solver%evals) <= eval_threshold) THEN
2057 : work(j)%info = -1
2058 : CYCLE
2059 : END IF
2060 : work(j)%values(:) = 1.0_dp/work(j)%solver%evals(:)
2061 : work(j)%r(:, :, 1) = REAL(work(j)%vectors, KIND=dp)
2062 : work(j)%r(:, :, 2) = AIMAG(work(j)%vectors)
2063 : CALL qs_ot_spectral_matrix_local(work(j), matrix=work(j)%os)
2064 : END DO
2065 : !$OMP END PARALLEL DO
2066 816 : DO j = 1, nbatch
2067 334 : IF (work(j)%info /= 0) CYCLE
2068 334 : i = channels(first + j - 1)
2069 334 : CALL qs_ot_store_local_pair(qs_ot_env(i)%matrix_os, qs_ot_env(i)%matrix_os_im, work(j)%os)
2070 668 : qs_ot_env(i)%os_valid = .TRUE.
2071 : END DO
2072 : END DO
2073 296 : DO j = 1, nworkers
2074 296 : CALL work(j)%gemm_ctx%destroy()
2075 : END DO
2076 97016 : END SUBROUTINE qs_ot_prepare_complex_tangent_metric_batch
2077 :
2078 : ! **************************************************************************************************
2079 : !> \brief finite complex STRICT derivative, projected onto C0^H*S*X=0
2080 : !> \param matrix_hc real part of H(k)*C(k)
2081 : !> \param matrix_hc_im imaginary part of H(k)*C(k)
2082 : !> \param qs_ot_env OT channel state
2083 : ! **************************************************************************************************
2084 2612 : SUBROUTINE qs_ot_get_derivative_complex(matrix_hc, matrix_hc_im, qs_ot_env)
2085 : TYPE(dbcsr_type), POINTER :: matrix_hc, matrix_hc_im
2086 : TYPE(qs_ot_type) :: qs_ot_env
2087 :
2088 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_ot_get_derivative_complex'
2089 :
2090 : INTEGER :: handle
2091 : TYPE(dbcsr_distribution_type) :: dist
2092 : TYPE(dbcsr_type) :: tmp_nk
2093 : TYPE(dbcsr_type), POINTER :: hc_work_im, hc_work_re
2094 : TYPE(dbcsr_type), TARGET :: hc_rot_im, hc_rot_re
2095 :
2096 1306 : CALL timeset(routineN, handle)
2097 :
2098 1306 : CPASSERT(qs_ot_env%has_complex_kpoint_state)
2099 1306 : CPASSERT(qs_ot_env%settings%ot_algorithm == 'TOD')
2100 :
2101 1306 : hc_work_re => matrix_hc
2102 1306 : hc_work_im => matrix_hc_im
2103 :
2104 1306 : IF (qs_ot_env%settings%do_rotation) THEN
2105 634 : CALL qs_ot_generate_rotation_complex(qs_ot_env)
2106 :
2107 : ! Reconstruct the unrotated finite STRICT orbitals Q(X).
2108 : CALL qs_ot_complex_multiply('N', 'N', qs_ot_env%matrix_c0, qs_ot_env%matrix_c0_im, &
2109 : qs_ot_env%matrix_cosp, qs_ot_env%matrix_cosp_im, &
2110 : qs_ot_env%matrix_buf_nk, qs_ot_env%matrix_buf_nk_im, &
2111 634 : qs_ot_env%matrix_tmp_nk)
2112 : CALL qs_ot_complex_multiply('N', 'N', qs_ot_env%matrix_x, qs_ot_env%matrix_x_im, &
2113 : qs_ot_env%matrix_sinp, qs_ot_env%matrix_sinp_im, &
2114 : qs_ot_env%matrix_gx, qs_ot_env%matrix_gx_im, &
2115 634 : qs_ot_env%matrix_tmp_nk)
2116 : CALL dbcsr_add(qs_ot_env%matrix_buf_nk, qs_ot_env%matrix_gx, &
2117 634 : alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
2118 : CALL dbcsr_add(qs_ot_env%matrix_buf_nk_im, qs_ot_env%matrix_gx_im, &
2119 634 : alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
2120 :
2121 : ! dF/dU = Q(X)^H*G_C for C=Q(X)*U.
2122 : CALL qs_ot_complex_multiply('C', 'N', &
2123 : qs_ot_env%matrix_buf_nk, qs_ot_env%matrix_buf_nk_im, &
2124 : matrix_hc, matrix_hc_im, &
2125 : qs_ot_env%rot_mat_dedu, qs_ot_env%rot_mat_dedu_im, &
2126 634 : qs_ot_env%matrix_buf1)
2127 634 : CALL qs_ot_rot_mat_derivative_complex(qs_ot_env)
2128 :
2129 : ! The STRICT coordinate sees G_Q=G_C*U^H.
2130 634 : CALL dbcsr_copy(hc_rot_re, matrix_hc, name="strict_hc_rot_re")
2131 634 : CALL dbcsr_copy(hc_rot_im, matrix_hc_im, name="strict_hc_rot_im")
2132 634 : CALL dbcsr_copy(tmp_nk, matrix_hc, name="strict_hc_rot_tmp")
2133 : CALL dbcsr_multiply('N', 'T', 1.0_dp, matrix_hc, qs_ot_env%rot_mat_u, &
2134 634 : 0.0_dp, hc_rot_re)
2135 : CALL dbcsr_multiply('N', 'T', 1.0_dp, matrix_hc_im, qs_ot_env%rot_mat_u_im, &
2136 634 : 0.0_dp, tmp_nk)
2137 634 : CALL dbcsr_add(hc_rot_re, tmp_nk, alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
2138 :
2139 : CALL dbcsr_multiply('N', 'T', 1.0_dp, matrix_hc_im, qs_ot_env%rot_mat_u, &
2140 634 : 0.0_dp, hc_rot_im)
2141 : CALL dbcsr_multiply('N', 'T', 1.0_dp, matrix_hc, qs_ot_env%rot_mat_u_im, &
2142 634 : 0.0_dp, tmp_nk)
2143 634 : CALL dbcsr_add(hc_rot_im, tmp_nk, alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
2144 634 : hc_work_re => hc_rot_re
2145 634 : hc_work_im => hc_rot_im
2146 : END IF
2147 :
2148 : ! Direct X contribution, H*C sinc(sqrt(P)).
2149 : CALL qs_ot_complex_multiply('N', 'N', hc_work_re, hc_work_im, &
2150 : qs_ot_env%matrix_sinp, qs_ot_env%matrix_sinp_im, &
2151 : qs_ot_env%matrix_gx, qs_ot_env%matrix_gx_im, &
2152 1306 : qs_ot_env%matrix_tmp_nk)
2153 :
2154 : CALL qs_ot_strict_frechet(hc_work_re, hc_work_im, qs_ot_env%matrix_x, qs_ot_env%matrix_x_im, &
2155 1306 : qs_ot_env%matrix_sinp_b, qs_ot_env%matrix_buf3, qs_ot_env%matrix_buf3_im, qs_ot_env)
2156 : CALL qs_ot_strict_frechet(hc_work_re, hc_work_im, qs_ot_env%matrix_c0, qs_ot_env%matrix_c0_im, &
2157 1306 : qs_ot_env%matrix_cosp_b, qs_ot_env%matrix_buf4, qs_ot_env%matrix_buf4_im, qs_ot_env)
2158 : CALL dbcsr_add(qs_ot_env%matrix_buf3, qs_ot_env%matrix_buf4, &
2159 1306 : alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
2160 : CALL dbcsr_add(qs_ot_env%matrix_buf3_im, qs_ot_env%matrix_buf4_im, &
2161 1306 : alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
2162 :
2163 : ! Transform back and add the Hermitian adjoint generated by dP.
2164 : CALL qs_ot_complex_multiply('N', 'N', qs_ot_env%matrix_r, qs_ot_env%matrix_r_im, &
2165 : qs_ot_env%matrix_buf3, qs_ot_env%matrix_buf3_im, &
2166 : qs_ot_env%matrix_buf1, qs_ot_env%matrix_buf1_im, &
2167 1306 : qs_ot_env%matrix_buf2)
2168 : CALL qs_ot_complex_multiply('N', 'C', qs_ot_env%matrix_buf1, qs_ot_env%matrix_buf1_im, &
2169 : qs_ot_env%matrix_r, qs_ot_env%matrix_r_im, &
2170 : qs_ot_env%matrix_buf3, qs_ot_env%matrix_buf3_im, &
2171 1306 : qs_ot_env%matrix_buf2)
2172 1306 : CALL dbcsr_get_info(qs_ot_env%matrix_buf3, distribution=dist)
2173 : CALL dbcsr_transposed(qs_ot_env%matrix_buf4, qs_ot_env%matrix_buf3, &
2174 : shallow_data_copy=.FALSE., use_distribution=dist, &
2175 1306 : transpose_distribution=.FALSE.)
2176 : CALL dbcsr_transposed(qs_ot_env%matrix_buf4_im, qs_ot_env%matrix_buf3_im, &
2177 : shallow_data_copy=.FALSE., use_distribution=dist, &
2178 1306 : transpose_distribution=.FALSE.)
2179 : CALL dbcsr_add(qs_ot_env%matrix_buf3, qs_ot_env%matrix_buf4, &
2180 1306 : alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
2181 : CALL dbcsr_add(qs_ot_env%matrix_buf3_im, qs_ot_env%matrix_buf4_im, &
2182 1306 : alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
2183 :
2184 : CALL qs_ot_complex_multiply('N', 'N', qs_ot_env%matrix_sx, qs_ot_env%matrix_sx_im, &
2185 : qs_ot_env%matrix_buf3, qs_ot_env%matrix_buf3_im, &
2186 : qs_ot_env%matrix_buf_nk, qs_ot_env%matrix_buf_nk_im, &
2187 1306 : qs_ot_env%matrix_tmp_nk)
2188 : CALL dbcsr_add(qs_ot_env%matrix_gx, qs_ot_env%matrix_buf_nk, &
2189 1306 : alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
2190 : CALL dbcsr_add(qs_ot_env%matrix_gx_im, qs_ot_env%matrix_buf_nk_im, &
2191 1306 : alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
2192 :
2193 1306 : CALL qs_ot_project_gradient(qs_ot_env, qs_ot_env%matrix_gx, qs_ot_env%matrix_gx_im)
2194 :
2195 1306 : IF (qs_ot_env%settings%do_rotation) THEN
2196 634 : CALL dbcsr_release(hc_rot_re)
2197 634 : CALL dbcsr_release(hc_rot_im)
2198 634 : CALL dbcsr_release(tmp_nk)
2199 : END IF
2200 :
2201 1306 : CALL timestop(handle)
2202 :
2203 1306 : END SUBROUTINE qs_ot_get_derivative_complex
2204 :
2205 : ! **************************************************************************************************
2206 : !> \brief Whether all channels can use the single-rank STRICT/IRAC kernels.
2207 : !> \param qs_ot_env local OT channels
2208 : !> \return true when every channel is eligible
2209 : ! **************************************************************************************************
2210 99940 : LOGICAL FUNCTION qs_ot_local_channels(qs_ot_env) RESULT(local)
2211 : TYPE(qs_ot_type), INTENT(IN) :: qs_ot_env(:)
2212 :
2213 : INTEGER :: i
2214 :
2215 99940 : local = .FALSE.
2216 99940 : IF (SIZE(qs_ot_env) == 0) RETURN
2217 109574 : DO i = 1, SIZE(qs_ot_env)
2218 104768 : IF (.NOT. qs_ot_env(i)%has_complex_kpoint_state) RETURN
2219 : SELECT CASE (qs_ot_env(i)%settings%ot_algorithm)
2220 : CASE ('TOD', 'REF')
2221 : CASE DEFAULT
2222 11518 : RETURN
2223 : END SELECT
2224 11518 : IF (qs_ot_env(i)%settings%ot_algorithm /= qs_ot_env(1)%settings%ot_algorithm) RETURN
2225 16324 : IF (qs_ot_env(i)%para_env%num_pe /= 1) RETURN
2226 : END DO
2227 99940 : local = .TRUE.
2228 : END FUNCTION qs_ot_local_channels
2229 :
2230 : ! **************************************************************************************************
2231 : !> \brief Rebuild single-rank STRICT/IRAC orbitals in bounded batches using cached overlaps.
2232 : !> DBCSR packing and writeback are outside OpenMP; workers only use dense data and private scratch.
2233 : !> \param qs_ot_env local OT channels
2234 : !> \param matrix_s pointers to cached real/imaginary overlaps, indexed by component and channel
2235 : !> \param matrix_c pointers to physical orbitals; only their active columns are updated
2236 : !> \param initialize build C0/SC0 from physical orbitals without modifying the physical state
2237 : ! **************************************************************************************************
2238 1816 : SUBROUTINE qs_ot_get_orbitals_local_batch(qs_ot_env, matrix_s, matrix_c, initialize)
2239 : TYPE(qs_ot_type), INTENT(INOUT) :: qs_ot_env(:)
2240 : TYPE(cp_fm_p_type), INTENT(IN) :: matrix_s(:, :)
2241 : TYPE(cp_cfm_p_type), INTENT(IN) :: matrix_c(:)
2242 : LOGICAL, INTENT(IN), OPTIONAL :: initialize
2243 :
2244 : CHARACTER(LEN=*), PARAMETER :: routineN = 'qs_ot_get_orbitals_local_batch'
2245 :
2246 : INTEGER :: first, handle, i, j, k, n, nbatch, &
2247 : nchannels, nworkers
2248 : LOGICAL :: initial, ref, update
2249 1816 : TYPE(qs_ot_local_work_type), ALLOCATABLE :: work(:)
2250 :
2251 1816 : initial = .FALSE.
2252 1816 : IF (PRESENT(initialize)) initial = initialize
2253 1816 : CPASSERT(qs_ot_local_channels(qs_ot_env))
2254 1816 : ref = qs_ot_env(1)%settings%ot_algorithm == 'REF'
2255 1816 : update = .FALSE.
2256 1816 : IF (ref) update = qs_ot_ref_update(qs_ot_env(1))
2257 1816 : CALL timeset(routineN, handle)
2258 1816 : nchannels = SIZE(qs_ot_env)
2259 1816 : nworkers = 1
2260 1816 : !$ nworkers = MIN(nchannels, omp_get_max_threads())
2261 7264 : ALLOCATE (work(nworkers))
2262 3632 : DO j = 1, nworkers
2263 3632 : CALL work(j)%gemm_ctx%create(LOCAL_GEMM_PU_GPU, timing=.FALSE.)
2264 : END DO
2265 5446 : DO first = 1, nchannels, nworkers
2266 3630 : nbatch = MIN(nworkers, nchannels - first + 1)
2267 7260 : DO j = 1, nbatch
2268 3630 : i = first + j - 1
2269 3630 : CALL dbcsr_get_info(qs_ot_env(i)%matrix_x, nfullrows_total=n, nfullcols_total=k)
2270 3630 : CALL qs_ot_allocate_local_work(work(j), n, k)
2271 3630 : work(j)%do_rotation = qs_ot_env(i)%settings%do_rotation .AND. .NOT. initial
2272 3630 : IF (.NOT. initial .OR. ref) CALL qs_ot_prepare_local_spectrum(work(j))
2273 3630 : IF (work(j)%do_rotation) THEN
2274 951 : CALL qs_ot_load_local_pair(qs_ot_env(i)%rot_mat_x, qs_ot_env(i)%rot_mat_x_im, work(j)%rot)
2275 : END IF
2276 10890 : IF (initial) THEN
2277 77555 : work(j)%x(:, :, 1) = REAL(matrix_c(i)%matrix%local_data(1:n, 1:k), KIND=dp)
2278 77555 : work(j)%x(:, :, 2) = AIMAG(matrix_c(i)%matrix%local_data(1:n, 1:k))
2279 : ELSE
2280 3231 : CALL qs_ot_load_local_pair(qs_ot_env(i)%matrix_x, qs_ot_env(i)%matrix_x_im, work(j)%x)
2281 3231 : IF (.NOT. ref) THEN
2282 1842 : CALL qs_ot_load_local_pair(qs_ot_env(i)%matrix_c0, qs_ot_env(i)%matrix_c0_im, work(j)%c0)
2283 : END IF
2284 : END IF
2285 : END DO
2286 : !$OMP PARALLEL DO DEFAULT(NONE) SCHEDULE(STATIC) NUM_THREADS(nbatch) IF(nbatch > 1) &
2287 3630 : !$OMP SHARED(work,matrix_s,matrix_c,nbatch,first,ref,UPDATE,initial) PRIVATE(i,n,k)
2288 : DO j = 1, nbatch
2289 : i = first + j - 1
2290 : IF (work(j)%do_rotation) THEN
2291 : CALL qs_ot_generate_rotation_local(work(j))
2292 : IF (work(j)%info /= 0) CYCLE
2293 : END IF
2294 : CALL qs_ot_orbitals_local(work(j), matrix_s(1, i)%matrix%local_data, &
2295 : matrix_s(2, i)%matrix%local_data, ref, update, initial)
2296 : IF (work(j)%info /= 0) CYCLE
2297 : IF (initial) CYCLE
2298 : IF (work(j)%do_rotation) THEN
2299 : ! IRAC retraction updates X to Q; rotate only the physical orbitals.
2300 : CALL qs_ot_complex_multiply_local('N', 'N', work(j)%out, work(j)%rot, &
2301 : work(j)%buf_nk, work(j)%tmp, work(j)%gemm_ctx)
2302 : work(j)%out(:, :, :) = work(j)%buf_nk(:, :, :)
2303 : END IF
2304 : n = work(j)%n
2305 : k = work(j)%k
2306 : matrix_c(i)%matrix%local_data(1:n, 1:k) = CMPLX(work(j)%out(:, :, 1), work(j)%out(:, :, 2), KIND=dp)
2307 : END DO
2308 : !$OMP END PARALLEL DO
2309 9076 : DO j = 1, nbatch
2310 3630 : IF (work(j)%info == -1) CPABORT('Complex REF overlap is not positive definite')
2311 3630 : IF (work(j)%info /= 0) CPABORT('Local complex OT diagonalization failed')
2312 3630 : i = first + j - 1
2313 3630 : IF (work(j)%do_rotation) THEN
2314 951 : CALL qs_ot_store_local_pair(qs_ot_env(i)%rot_mat_u, qs_ot_env(i)%rot_mat_u_im, work(j)%rot)
2315 : END IF
2316 3630 : IF (initial) THEN
2317 399 : CALL qs_ot_store_local_pair(qs_ot_env(i)%matrix_c0, qs_ot_env(i)%matrix_c0_im, work(j)%out)
2318 399 : CALL qs_ot_store_local_pair(qs_ot_env(i)%matrix_sc0, qs_ot_env(i)%matrix_sc0_im, work(j)%buf_nk)
2319 399 : IF (.NOT. ref) CYCLE
2320 : END IF
2321 3296 : CALL qs_ot_store_local_pair(qs_ot_env(i)%matrix_sx, qs_ot_env(i)%matrix_sx_im, work(j)%sx)
2322 3296 : IF (ref) THEN
2323 : CALL qs_ot_store_local_pair(qs_ot_env(i)%matrix_ref_inv_sqrt, &
2324 1454 : qs_ot_env(i)%matrix_ref_inv_sqrt_im, work(j)%sinp)
2325 1454 : IF (update .OR. initial) CALL qs_ot_store_local_pair(qs_ot_env(i)%matrix_x, qs_ot_env(i)%matrix_x_im, work(j)%x)
2326 : CYCLE
2327 : END IF
2328 1842 : CALL qs_ot_store_local_pair(qs_ot_env(i)%matrix_p, qs_ot_env(i)%matrix_p_im, work(j)%buf3)
2329 1842 : CALL qs_ot_store_local_pair(qs_ot_env(i)%matrix_r, qs_ot_env(i)%matrix_r_im, work(j)%r)
2330 1842 : CALL qs_ot_store_local_pair(qs_ot_env(i)%matrix_cosp, qs_ot_env(i)%matrix_cosp_im, work(j)%buf4)
2331 1842 : CALL qs_ot_store_local_pair(qs_ot_env(i)%matrix_sinp, qs_ot_env(i)%matrix_sinp_im, work(j)%sinp)
2332 1842 : CALL copy_local_to_dbcsr(qs_ot_env(i)%matrix_cosp_b, work(j)%cosp_b)
2333 1842 : CALL copy_local_to_dbcsr(qs_ot_env(i)%matrix_sinp_b, work(j)%sinp_b)
2334 11590 : qs_ot_env(i)%evals(:) = work(j)%solver%evals(:)
2335 5806 : qs_ot_env(i)%do_taylor = .FALSE.
2336 : END DO
2337 : END DO
2338 3632 : DO j = 1, nworkers
2339 3632 : CALL work(j)%gemm_ctx%destroy()
2340 : END DO
2341 1816 : CALL timestop(handle)
2342 5448 : END SUBROUTINE qs_ot_get_orbitals_local_batch
2343 :
2344 : ! **************************************************************************************************
2345 : !> \brief Evaluate S*X and reconstruct physical STRICT or IRAC orbitals.
2346 : !> \param w private storage and status; buf3 retains P and buf4 retains cos(sqrt(P)) for writeback
2347 : !> \param overlap_re cached real overlap, including any leading-dimension padding
2348 : !> \param overlap_im cached imaginary overlap
2349 : !> \param ref select IRAC instead of STRICT
2350 : !> \param update accept the IRAC coordinate at the existing minimizer update point
2351 : !> \param initialize retain SC0 and use the analytic zero-coordinate STRICT state
2352 : ! **************************************************************************************************
2353 3630 : SUBROUTINE qs_ot_orbitals_local(w, overlap_re, overlap_im, ref, update, initialize)
2354 : TYPE(qs_ot_local_work_type), INTENT(INOUT) :: w
2355 : REAL(dp), CONTIGUOUS, INTENT(IN) :: overlap_re(:, :), overlap_im(:, :)
2356 : LOGICAL, INTENT(IN) :: ref, update
2357 : LOGICAL, INTENT(IN), OPTIONAL :: initialize
2358 :
2359 : INTEGER :: i
2360 :
2361 : CALL qs_ot_complex_multiply_local('N', 'N', overlap_re, overlap_im, w%x(:, :, 1), w%x(:, :, 2), &
2362 3630 : w%sx(:, :, 1), w%sx(:, :, 2), w%tmp, w%gemm_ctx)
2363 3630 : IF (PRESENT(initialize)) THEN
2364 3630 : IF (initialize) THEN
2365 : ! Keep the incoming SC0, as in the distributed initializer. IRAC may
2366 : ! subsequently normalize C0 and accept a different X/SX pair.
2367 155509 : w%buf_nk(:, :, :) = w%sx(:, :, :)
2368 399 : IF (.NOT. ref) THEN
2369 82202 : w%out(:, :, :) = w%x(:, :, :)
2370 334 : w%info = 0
2371 334 : RETURN
2372 : END IF
2373 : END IF
2374 : END IF
2375 3296 : CALL qs_ot_complex_multiply_local('C', 'N', w%x, w%sx, w%buf3, w%buf2(:, :, 1), w%gemm_ctx)
2376 311196 : w%vectors(:, :) = CMPLX(w%buf3(:, :, 1), w%buf3(:, :, 2), KIND=dp)
2377 3296 : CALL cp_cfm_heevd_local(w%vectors, info=w%info, ws=w%solver)
2378 3296 : IF (w%info /= 0) RETURN
2379 311196 : w%r(:, :, 1) = REAL(w%vectors, KIND=dp)
2380 311196 : w%r(:, :, 2) = AIMAG(w%vectors)
2381 3296 : IF (ref) THEN
2382 15780 : IF (ANY(w%solver%evals <= EPSILON(1.0_dp))) THEN
2383 0 : w%info = -1
2384 0 : RETURN
2385 : END IF
2386 15780 : w%values(:) = 1.0_dp/SQRT(w%solver%evals(:))
2387 1454 : CALL qs_ot_spectral_matrix_local(w, matrix=w%sinp)
2388 1454 : CALL qs_ot_complex_multiply_local('N', 'N', w%x, w%sinp, w%out, w%tmp, w%gemm_ctx)
2389 1454 : IF (update) THEN
2390 792496 : w%x(:, :, :) = w%out(:, :, :)
2391 192320 : w%sinp(:, :, :) = 0.0_dp
2392 7494 : DO i = 1, w%k
2393 7494 : w%sinp(i, i, 1) = 1.0_dp
2394 : END DO
2395 : CALL qs_ot_complex_multiply_local('N', 'N', overlap_re, overlap_im, w%x(:, :, 1), w%x(:, :, 2), &
2396 712 : w%sx(:, :, 1), w%sx(:, :, 2), w%tmp, w%gemm_ctx)
2397 : END IF
2398 1454 : RETURN
2399 : END IF
2400 11590 : w%solver%evals(:) = MAX(0.0_dp, w%solver%evals(:))
2401 1842 : CALL qs_ot_spectral_matrix_local(w, .TRUE., w%buf4)
2402 1842 : CALL qs_ot_spectral_matrix_local(w, .FALSE., w%sinp)
2403 11590 : w%values(:) = SQRT(w%solver%evals(:))
2404 1842 : CALL qs_ot_spectral_kernel(w%values, w%values, w%cosp_b, .TRUE.)
2405 1842 : CALL qs_ot_spectral_kernel(w%values, w%values, w%sinp_b, .FALSE.)
2406 1842 : CALL qs_ot_complex_multiply_local('N', 'N', w%c0, w%buf4, w%out, w%tmp, w%gemm_ctx)
2407 1842 : CALL qs_ot_complex_multiply_local('N', 'N', w%x, w%sinp, w%buf_nk, w%tmp, w%gemm_ctx)
2408 622774 : w%out(:, :, :) = w%out(:, :, :) + w%buf_nk(:, :, :)
2409 : END SUBROUTINE qs_ot_orbitals_local
2410 :
2411 : ! **************************************************************************************************
2412 : !> \brief Reconstruct one spectral matrix function from the local eigenvectors.
2413 : !> \param w private eigenvectors, eigenvalues and scratch
2414 : !> \param cosine select cos/sinc(sqrt(P)); if absent, use the supplied spectral values
2415 : !> \param matrix dense real/imaginary result
2416 : ! **************************************************************************************************
2417 5472 : SUBROUTINE qs_ot_spectral_matrix_local(w, cosine, matrix)
2418 : TYPE(qs_ot_local_work_type), INTENT(INOUT) :: w
2419 : LOGICAL, INTENT(IN), OPTIONAL :: cosine
2420 : REAL(dp), CONTIGUOUS, INTENT(INOUT) :: matrix(:, :, :)
2421 :
2422 : INTEGER :: i
2423 :
2424 5472 : IF (PRESENT(cosine)) CALL qs_ot_spectral_values(w%solver%evals, w%values, cosine)
2425 40778 : DO i = 1, w%k
2426 417244 : w%buf1(:, i, 1) = w%r(:, i, 1)*w%values(i)
2427 422716 : w%buf1(:, i, 2) = w%r(:, i, 2)*w%values(i)
2428 : END DO
2429 5472 : CALL qs_ot_complex_multiply_local('N', 'C', w%buf1, w%r, matrix, w%buf2(:, :, 1), w%gemm_ctx)
2430 5472 : END SUBROUTINE qs_ot_spectral_matrix_local
2431 :
2432 : ! **************************************************************************************************
2433 : !> \brief Allocate common dense OT storage only when a worker's matrix dimensions change.
2434 : !> \param w worker storage
2435 : !> \param n AO count
2436 : !> \param k active MO count
2437 : ! **************************************************************************************************
2438 6301 : SUBROUTINE qs_ot_allocate_local_work(w, n, k)
2439 : TYPE(qs_ot_local_work_type), INTENT(INOUT) :: w
2440 : INTEGER, INTENT(IN) :: n, k
2441 :
2442 6301 : IF (ALLOCATED(w%x)) THEN
2443 3229 : IF (n == w%n .AND. k == w%k) RETURN
2444 0 : DEALLOCATE (w%x, w%sx, w%c0, w%out, w%r, w%sinp, w%buf1, w%buf2, w%buf3, w%buf4, &
2445 162 : w%buf_nk, w%sinp_b, w%cosp_b, w%tmp)
2446 : END IF
2447 0 : ALLOCATE (w%x(n, k, 2), w%sx(n, k, 2), w%c0(n, k, 2), w%out(n, k, 2), &
2448 0 : w%r(k, k, 2), w%sinp(k, k, 2), w%buf1(k, k, 2), w%buf2(k, k, 2), &
2449 0 : w%buf3(k, k, 2), w%buf4(k, k, 2), w%buf_nk(n, k, 2), &
2450 113190 : w%sinp_b(k, k), w%cosp_b(k, k), w%tmp(n, k))
2451 3234 : w%n = n
2452 3234 : w%k = k
2453 : END SUBROUTINE qs_ot_allocate_local_work
2454 :
2455 : ! **************************************************************************************************
2456 : !> \brief Prepare a worker's orbital and optional rotation eigensolvers outside OpenMP.
2457 : !> \param w worker with known matrix dimensions
2458 : ! **************************************************************************************************
2459 3817 : SUBROUTINE qs_ot_prepare_local_spectrum(w)
2460 : TYPE(qs_ot_local_work_type), INTENT(INOUT) :: w
2461 :
2462 3817 : IF (ALLOCATED(w%vectors)) THEN
2463 1791 : IF (SIZE(w%vectors, 1) /= w%k) DEALLOCATE (w%vectors, w%values)
2464 : END IF
2465 14487 : IF (.NOT. ALLOCATED(w%vectors)) ALLOCATE (w%vectors(w%k, w%k), w%values(w%k))
2466 3817 : CALL cp_cfm_local_workspace_prepare(w%solver, w%k, w%vectors)
2467 3817 : IF (w%do_rotation) THEN
2468 1472 : CALL qs_ot_prepare_rotation_work(w%rotation, w%k)
2469 1472 : IF (ALLOCATED(w%rot)) THEN
2470 394 : IF (SIZE(w%rot, 1) /= w%k) DEALLOCATE (w%rot, w%rot_g)
2471 : END IF
2472 7940 : IF (.NOT. ALLOCATED(w%rot)) ALLOCATE (w%rot(w%k, w%k, 2), w%rot_g(w%k, w%k, 2))
2473 : END IF
2474 3817 : END SUBROUTINE qs_ot_prepare_local_spectrum
2475 :
2476 : ! **************************************************************************************************
2477 : !> \brief Evaluate a worker's rotation and retain its spectrum for the gradient pullback.
2478 : !> \param w prepared storage; rot contains A on entry and U on return
2479 : ! **************************************************************************************************
2480 1472 : SUBROUTINE qs_ot_generate_rotation_local(w)
2481 : TYPE(qs_ot_local_work_type), INTENT(INOUT) :: w
2482 :
2483 228944 : w%vectors(:, :) = CMPLX(w%rot(:, :, 1), w%rot(:, :, 2), KIND=dp)
2484 1472 : CALL qs_ot_rotation_state_local(w%vectors, w%rotation, w%gemm_ctx)
2485 1472 : w%info = w%rotation%info
2486 1472 : IF (w%info /= 0) RETURN
2487 228944 : w%rot(:, :, 1) = REAL(w%rotation%u, KIND=dp)
2488 228944 : w%rot(:, :, 2) = AIMAG(w%rotation%u)
2489 : END SUBROUTINE qs_ot_generate_rotation_local
2490 :
2491 : ! **************************************************************************************************
2492 : !> \brief Write a dense complex pair back to two DBCSR matrices outside OpenMP.
2493 : !> \param matrix_re real DBCSR output
2494 : !> \param matrix_im imaginary DBCSR output
2495 : !> \param dense dense real/imaginary pair
2496 : ! **************************************************************************************************
2497 18783 : SUBROUTINE qs_ot_store_local_pair(matrix_re, matrix_im, dense)
2498 : TYPE(dbcsr_type), INTENT(INOUT) :: matrix_re, matrix_im
2499 : REAL(dp), INTENT(IN) :: dense(:, :, :)
2500 :
2501 18783 : CALL copy_local_to_dbcsr(matrix_re, dense(:, :, 1))
2502 18783 : CALL copy_local_to_dbcsr(matrix_im, dense(:, :, 2))
2503 18783 : END SUBROUTINE qs_ot_store_local_pair
2504 :
2505 : ! **************************************************************************************************
2506 : !> \brief Evaluate independent, single-rank STRICT/IRAC gradients in bounded OpenMP batches.
2507 : !> DBCSR packing, filtering and metric preparation stay outside the worker region.
2508 : !> \param qs_ot_env local OT channels
2509 : !> \param matrix_hc real orbital gradients
2510 : !> \param matrix_hc_im imaginary orbital gradients
2511 : !> \param matrix_hc_physical optional physical gradients for occupation preconditioning
2512 : !> \param matrix_hc_physical_im imaginary physical gradients
2513 : ! **************************************************************************************************
2514 2512 : SUBROUTINE qs_ot_get_derivative_local_batch(qs_ot_env, matrix_hc, matrix_hc_im, &
2515 1256 : matrix_hc_physical, matrix_hc_physical_im)
2516 : TYPE(qs_ot_type), INTENT(INOUT) :: qs_ot_env(:)
2517 : TYPE(dbcsr_p_type), INTENT(IN) :: matrix_hc(:), matrix_hc_im(:)
2518 : TYPE(dbcsr_p_type), INTENT(IN), OPTIONAL :: matrix_hc_physical(:), &
2519 : matrix_hc_physical_im(:)
2520 :
2521 : CHARACTER(LEN=*), PARAMETER :: routineN = 'qs_ot_get_derivative_local_batch'
2522 :
2523 : INTEGER :: first, handle, i, j, nbatch, nchannels, &
2524 : npasses, nworkers, pass
2525 : LOGICAL :: ref
2526 1256 : TYPE(qs_ot_local_work_type), ALLOCATABLE :: work(:)
2527 :
2528 1256 : CPASSERT(qs_ot_local_channels(qs_ot_env))
2529 1256 : ref = qs_ot_env(1)%settings%ot_algorithm == 'REF'
2530 1256 : npasses = 1
2531 1256 : IF (PRESENT(matrix_hc_physical)) THEN
2532 116 : CPASSERT(PRESENT(matrix_hc_physical_im))
2533 : npasses = 2
2534 : END IF
2535 1256 : nchannels = SIZE(qs_ot_env)
2536 :
2537 1256 : CALL timeset(routineN, handle)
2538 1256 : nworkers = 1
2539 1256 : !$ nworkers = MIN(nchannels, omp_get_max_threads())
2540 5024 : ALLOCATE (work(nworkers))
2541 2512 : DO j = 1, nworkers
2542 2512 : CALL work(j)%gemm_ctx%create(LOCAL_GEMM_PU_GPU, timing=.FALSE.)
2543 : END DO
2544 3927 : DO first = 1, nchannels, nworkers
2545 2671 : nbatch = MIN(nworkers, nchannels - first + 1)
2546 5342 : DO j = 1, nbatch
2547 2671 : i = first + j - 1
2548 2671 : IF (.NOT. ref) CALL qs_ot_prepare_complex_tangent_metric(qs_ot_env(i))
2549 2671 : CALL qs_ot_load_local_work(work(j), qs_ot_env(i), matrix_hc(i)%matrix, matrix_hc_im(i)%matrix)
2550 8013 : work(j)%filtered = ref .AND. qs_ot_env(i)%settings%eps_irac_filter_matrix > 0.0_dp
2551 : END DO
2552 6714 : DO pass = 1, npasses
2553 2787 : IF (pass == 2) THEN
2554 232 : DO j = 1, nbatch
2555 116 : i = first + j - 1
2556 232 : CALL qs_ot_load_local_pair(matrix_hc_physical(i)%matrix, matrix_hc_physical_im(i)%matrix, work(j)%hc)
2557 : END DO
2558 : END IF
2559 : !$OMP PARALLEL DO DEFAULT(NONE) SCHEDULE(STATIC) NUM_THREADS(nbatch) IF(nbatch > 1) &
2560 2787 : !$OMP SHARED(work,nbatch,ref,pass,npasses)
2561 : DO j = 1, nbatch
2562 : IF (ref) THEN
2563 : CALL qs_ot_complex_multiply_local('N', 'N', work(j)%x, work(j)%sinp, &
2564 : work(j)%c0, work(j)%tmp, work(j)%gemm_ctx)
2565 : END IF
2566 : IF (work(j)%do_rotation) THEN
2567 : CALL qs_ot_rotation_gradient_local(work(j), ref, prepare=pass == 1, physical=pass == npasses)
2568 : IF (work(j)%info /= 0) CYCLE
2569 : END IF
2570 : IF (ref) THEN
2571 : CALL qs_ot_complex_multiply_local('C', 'N', work(j)%c0, work(j)%hc, &
2572 : work(j)%buf3, work(j)%buf2(:, :, 1), work(j)%gemm_ctx)
2573 : IF (.NOT. work(j)%filtered) CALL qs_ot_ref_gradient_local(work(j))
2574 : ELSE
2575 : CALL qs_ot_strict_gradient_local(work(j))
2576 : END IF
2577 : END DO
2578 : !$OMP END PARALLEL DO
2579 5574 : IF (ANY(work(1:nbatch)%info == -1)) THEN
2580 0 : CPABORT('Complex REF inverse square root is not positive definite')
2581 : END IF
2582 5574 : IF (ANY(work(1:nbatch)%info /= 0)) THEN
2583 0 : CPABORT('Local complex OT diagonalization failed')
2584 : END IF
2585 5574 : IF (ANY(work(1:nbatch)%filtered)) THEN
2586 : ! Keep the existing block filter, including independent real/imaginary sparsity.
2587 0 : DO j = 1, nbatch
2588 0 : IF (.NOT. work(j)%filtered) CYCLE
2589 0 : i = first + j - 1
2590 0 : CALL qs_ot_store_local_pair(qs_ot_env(i)%buf1_k_k_sym, qs_ot_env(i)%buf2_k_k_sym, work(j)%buf3)
2591 0 : CALL dbcsr_filter(qs_ot_env(i)%buf1_k_k_sym, qs_ot_env(i)%settings%eps_irac_filter_matrix)
2592 0 : CALL dbcsr_filter(qs_ot_env(i)%buf2_k_k_sym, qs_ot_env(i)%settings%eps_irac_filter_matrix)
2593 0 : CALL qs_ot_load_local_pair(qs_ot_env(i)%buf1_k_k_sym, qs_ot_env(i)%buf2_k_k_sym, work(j)%buf3)
2594 : END DO
2595 0 : !$OMP PARALLEL DO DEFAULT(NONE) SCHEDULE(STATIC) NUM_THREADS(nbatch) IF(nbatch > 1) SHARED(work,nbatch)
2596 : DO j = 1, nbatch
2597 : IF (work(j)%filtered) CALL qs_ot_ref_gradient_local(work(j))
2598 : END DO
2599 : !$OMP END PARALLEL DO
2600 : END IF
2601 8245 : DO j = 1, nbatch
2602 2787 : i = first + j - 1
2603 2787 : IF (work(j)%do_rotation .AND. pass == npasses) THEN
2604 521 : CALL qs_ot_store_local_pair(qs_ot_env(i)%rot_mat_u, qs_ot_env(i)%rot_mat_u_im, work(j)%rot)
2605 521 : CALL qs_ot_store_local_pair(qs_ot_env(i)%rot_mat_gx, qs_ot_env(i)%rot_mat_gx_im, work(j)%rot_g)
2606 : END IF
2607 2787 : CALL qs_ot_store_local_pair(qs_ot_env(i)%matrix_gx, qs_ot_env(i)%matrix_gx_im, work(j)%out)
2608 2787 : IF (work(j)%filtered) THEN
2609 0 : CALL dbcsr_filter(qs_ot_env(i)%matrix_gx, qs_ot_env(i)%settings%eps_irac_filter_matrix)
2610 0 : CALL dbcsr_filter(qs_ot_env(i)%matrix_gx_im, qs_ot_env(i)%settings%eps_irac_filter_matrix)
2611 : END IF
2612 5574 : IF (pass < npasses) THEN
2613 116 : CALL dbcsr_copy(qs_ot_env(i)%matrix_preconditioned_gx, qs_ot_env(i)%matrix_gx)
2614 116 : CALL dbcsr_copy(qs_ot_env(i)%matrix_preconditioned_gx_im, qs_ot_env(i)%matrix_gx_im)
2615 : END IF
2616 : END DO
2617 : END DO
2618 : END DO
2619 2512 : DO j = 1, nworkers
2620 2512 : CALL work(j)%gemm_ctx%destroy()
2621 : END DO
2622 1256 : CALL timestop(handle)
2623 3768 : END SUBROUTINE qs_ot_get_derivative_local_batch
2624 :
2625 : ! **************************************************************************************************
2626 : !> \brief Pull back through rotation, sharing spectra across occupation-preconditioner passes.
2627 : !> \param w worker with G_C in hc; hc receives G_Q and rot_g receives the rotation gradient
2628 : !> \param ref select IRAC with Q already in c0 instead of STRICT with the reference C0 in c0
2629 : !> \param prepare compute the rotation and, for IRAC, polar spectra on the first pass
2630 : !> \param physical compute the rotation gradient only from the physical covector
2631 : ! **************************************************************************************************
2632 637 : SUBROUTINE qs_ot_rotation_gradient_local(w, ref, prepare, physical)
2633 : TYPE(qs_ot_local_work_type), INTENT(INOUT) :: w
2634 : LOGICAL, INTENT(IN) :: ref, prepare, physical
2635 :
2636 637 : IF (prepare) THEN
2637 521 : CALL qs_ot_generate_rotation_local(w)
2638 521 : IF (w%info /= 0) RETURN
2639 521 : IF (ref) THEN
2640 46497 : w%vectors(:, :) = CMPLX(w%sinp(:, :, 1), w%sinp(:, :, 2), KIND=dp)
2641 431 : CALL cp_cfm_heevd_local(w%vectors, info=w%info, ws=w%solver)
2642 431 : IF (w%info /= 0) RETURN
2643 4165 : IF (ANY(w%solver%evals <= EPSILON(1.0_dp))) THEN
2644 0 : w%info = -1
2645 0 : RETURN
2646 : END IF
2647 46497 : w%r(:, :, 1) = REAL(w%vectors, KIND=dp)
2648 46497 : w%r(:, :, 2) = AIMAG(w%vectors)
2649 : END IF
2650 : END IF
2651 637 : IF (physical) THEN
2652 521 : IF (ref) THEN
2653 431 : CALL qs_ot_complex_multiply_local('C', 'N', w%c0, w%hc, w%buf3, w%buf2(:, :, 1), w%gemm_ctx)
2654 : ELSE
2655 : ! Q(X) = C0*cos(sqrt(P)) + X*sinc(sqrt(P)), preserving C0 and the STRICT kernels.
2656 90 : CALL qs_ot_complex_multiply_local('N', 'N', w%c0, w%cosp, w%out, w%tmp, w%gemm_ctx)
2657 90 : CALL qs_ot_complex_multiply_local('N', 'N', w%x, w%sinp, w%buf_nk, w%tmp, w%gemm_ctx)
2658 153886 : w%out(:, :, :) = w%out(:, :, :) + w%buf_nk(:, :, :)
2659 90 : CALL qs_ot_complex_multiply_local('C', 'N', w%out, w%hc, w%buf3, w%buf2(:, :, 1), w%gemm_ctx)
2660 : END IF
2661 80063 : w%rotation%work(:, :) = CMPLX(w%buf3(:, :, 1), w%buf3(:, :, 2), KIND=dp)
2662 521 : CALL qs_ot_rotation_pullback_local(w%rotation, w%gemm_ctx, w%rot_g(:, :, 1), w%rot_g(:, :, 2))
2663 : END IF
2664 637 : CALL qs_ot_complex_multiply_local('N', 'C', w%hc, w%rot, w%buf_nk, w%tmp, w%gemm_ctx)
2665 558103 : w%hc(:, :, :) = w%buf_nk(:, :, :)
2666 : END SUBROUTINE qs_ot_rotation_gradient_local
2667 :
2668 : ! **************************************************************************************************
2669 : !> \brief Complete the IRAC pullback, including the vertical polar response for rotation.
2670 : !> \param w private storage with F in sinp and B in buf3
2671 : ! **************************************************************************************************
2672 1027 : SUBROUTINE qs_ot_ref_gradient_local(w)
2673 : TYPE(qs_ot_local_work_type), INTENT(INOUT) :: w
2674 :
2675 : INTEGER :: i, j
2676 : REAL(dp) :: denominator
2677 :
2678 1027 : CALL qs_ot_complex_multiply_local('N', 'N', w%sx, w%sinp, w%buf_nk, w%tmp, w%gemm_ctx)
2679 1027 : CALL qs_ot_complex_multiply_local('N', 'N', w%buf_nk, w%buf3, w%c0, w%tmp, w%gemm_ctx)
2680 1111177 : w%c0(:, :, :) = w%hc(:, :, :) - w%c0(:, :, :)
2681 1027 : CALL qs_ot_complex_multiply_local('N', 'N', w%c0, w%sinp, w%out, w%tmp, w%gemm_ctx)
2682 1027 : IF (.NOT. w%do_rotation) RETURN
2683 :
2684 : ! Solve P*Z+Z*P=antiherm(B), P=F^(-1), using the shared eigensystem of F.
2685 99253 : w%buf1(:, :, 1) = 0.5_dp*(w%buf3(:, :, 1) - TRANSPOSE(w%buf3(:, :, 1)))
2686 99253 : w%buf1(:, :, 2) = 0.5_dp*(w%buf3(:, :, 2) + TRANSPOSE(w%buf3(:, :, 2)))
2687 521 : CALL qs_ot_complex_multiply_local('N', 'N', w%buf1, w%r, w%buf2, w%buf3(:, :, 1), w%gemm_ctx)
2688 521 : CALL qs_ot_complex_multiply_local('C', 'N', w%r, w%buf2, w%buf1, w%buf3(:, :, 1), w%gemm_ctx)
2689 4747 : DO j = 1, w%k
2690 49887 : DO i = 1, w%k
2691 45140 : denominator = 1.0_dp/w%solver%evals(i) + 1.0_dp/w%solver%evals(j)
2692 139646 : w%buf1(i, j, :) = w%buf1(i, j, :)/denominator
2693 : END DO
2694 : END DO
2695 521 : CALL qs_ot_complex_multiply_local('N', 'N', w%r, w%buf1, w%buf2, w%buf3(:, :, 1), w%gemm_ctx)
2696 521 : CALL qs_ot_complex_multiply_local('N', 'C', w%buf2, w%r, w%buf1, w%buf3(:, :, 1), w%gemm_ctx)
2697 : ! buf_nk still contains S*Q. Add the adjoint vertical contribution 2*S*Q*Z.
2698 521 : CALL qs_ot_complex_multiply_local('N', 'N', w%buf_nk, w%buf1, w%c0, w%tmp, w%gemm_ctx)
2699 399771 : w%out(:, :, :) = w%out(:, :, :) + 2.0_dp*w%c0(:, :, :)
2700 : END SUBROUTINE qs_ot_ref_gradient_local
2701 :
2702 : ! **************************************************************************************************
2703 : !> \brief Contract one STRICT term with its Frechet kernel in the eigenvector basis.
2704 : !> \param matrix_hc real orbital gradient
2705 : !> \param matrix_hc_im imaginary orbital gradient
2706 : !> \param orbital real X or C0
2707 : !> \param orbital_im imaginary X or C0
2708 : !> \param kernel real sinc or cos divided differences
2709 : !> \param RESULT real contraction
2710 : !> \param result_im imaginary contraction
2711 : !> \param qs_ot_env channel eigenvectors and scratch
2712 : ! **************************************************************************************************
2713 2612 : SUBROUTINE qs_ot_strict_frechet(matrix_hc, matrix_hc_im, orbital, orbital_im, kernel, &
2714 : RESULT, result_im, qs_ot_env)
2715 : TYPE(dbcsr_type) :: matrix_hc, matrix_hc_im, orbital, &
2716 : orbital_im, kernel, RESULT, result_im
2717 : TYPE(qs_ot_type) :: qs_ot_env
2718 :
2719 : CALL qs_ot_complex_multiply('C', 'N', matrix_hc, matrix_hc_im, orbital, orbital_im, &
2720 2612 : qs_ot_env%matrix_buf2, qs_ot_env%matrix_buf2_im, qs_ot_env%matrix_buf1)
2721 : CALL qs_ot_complex_multiply('N', 'N', qs_ot_env%matrix_buf2, qs_ot_env%matrix_buf2_im, &
2722 : qs_ot_env%matrix_r, qs_ot_env%matrix_r_im, &
2723 2612 : qs_ot_env%matrix_buf1, qs_ot_env%matrix_buf1_im, qs_ot_env%matrix_buf4)
2724 : CALL qs_ot_complex_multiply('C', 'N', qs_ot_env%matrix_r, qs_ot_env%matrix_r_im, &
2725 : qs_ot_env%matrix_buf1, qs_ot_env%matrix_buf1_im, &
2726 2612 : qs_ot_env%matrix_buf2, qs_ot_env%matrix_buf2_im, qs_ot_env%matrix_buf4)
2727 2612 : CALL dbcsr_hadamard_product(qs_ot_env%matrix_buf2, kernel, RESULT)
2728 2612 : CALL dbcsr_hadamard_product(qs_ot_env%matrix_buf2_im, kernel, result_im)
2729 2612 : END SUBROUTINE qs_ot_strict_frechet
2730 :
2731 : ! **************************************************************************************************
2732 : !> \brief Reuse one worker's scratch and pack the local kernel inputs.
2733 : !> \param w worker storage
2734 : !> \param qs_ot_env channel with an already prepared tangent metric
2735 : !> \param matrix_hc real orbital gradient
2736 : !> \param matrix_hc_im imaginary orbital gradient
2737 : ! **************************************************************************************************
2738 4405 : SUBROUTINE qs_ot_load_local_work(w, qs_ot_env, matrix_hc, matrix_hc_im)
2739 : TYPE(qs_ot_local_work_type), INTENT(INOUT) :: w
2740 : TYPE(qs_ot_type), INTENT(IN) :: qs_ot_env
2741 : TYPE(dbcsr_type), INTENT(IN) :: matrix_hc, matrix_hc_im
2742 :
2743 : INTEGER :: k, n, o
2744 :
2745 2671 : CALL dbcsr_get_info(qs_ot_env%matrix_x, nfullrows_total=n, nfullcols_total=k)
2746 2671 : CALL qs_ot_allocate_local_work(w, n, k)
2747 2671 : IF (ALLOCATED(w%hc)) THEN
2748 1415 : IF (SIZE(w%hc, 1) /= n .OR. SIZE(w%hc, 2) /= k) DEALLOCATE (w%hc)
2749 : END IF
2750 7911 : IF (.NOT. ALLOCATED(w%hc)) ALLOCATE (w%hc(n, k, 2))
2751 2671 : CALL qs_ot_load_local_pair(matrix_hc, matrix_hc_im, w%hc)
2752 2671 : CALL qs_ot_load_local_pair(qs_ot_env%matrix_x, qs_ot_env%matrix_x_im, w%x)
2753 2671 : CALL qs_ot_load_local_pair(qs_ot_env%matrix_sx, qs_ot_env%matrix_sx_im, w%sx)
2754 2671 : w%do_rotation = qs_ot_env%settings%do_rotation
2755 2671 : IF (w%do_rotation) THEN
2756 521 : CALL qs_ot_prepare_local_spectrum(w)
2757 521 : CALL qs_ot_load_local_pair(qs_ot_env%rot_mat_x, qs_ot_env%rot_mat_x_im, w%rot)
2758 : END IF
2759 2671 : IF (qs_ot_env%settings%ot_algorithm == 'REF') THEN
2760 937 : CALL qs_ot_load_local_pair(qs_ot_env%matrix_ref_inv_sqrt, qs_ot_env%matrix_ref_inv_sqrt_im, w%sinp)
2761 937 : RETURN
2762 : END IF
2763 1734 : IF (w%do_rotation) THEN
2764 90 : IF (ALLOCATED(w%cosp)) THEN
2765 0 : IF (SIZE(w%cosp, 1) /= k) DEALLOCATE (w%cosp)
2766 : END IF
2767 450 : IF (.NOT. ALLOCATED(w%cosp)) ALLOCATE (w%cosp(k, k, 2))
2768 90 : CALL qs_ot_load_local_pair(qs_ot_env%matrix_cosp, qs_ot_env%matrix_cosp_im, w%cosp)
2769 : END IF
2770 1734 : CALL dbcsr_get_info(qs_ot_env%matrix_sc0, nfullcols_total=o)
2771 1734 : IF (ALLOCATED(w%sc0)) THEN
2772 1050 : IF (SIZE(w%sc0, 1) /= n .OR. SIZE(w%buf1_ortho, 2) /= k .OR. o /= SIZE(w%sc0, 2)) THEN
2773 0 : DEALLOCATE (w%sc0, w%target, w%os, w%buf1_ortho, w%buf2_ortho, w%tmp_ortho)
2774 : END IF
2775 : END IF
2776 1734 : IF (.NOT. ALLOCATED(w%sc0)) THEN
2777 0 : ALLOCATE (w%sc0(n, o, 2), w%target(n, o, 2), w%os(o, o, 2), &
2778 13680 : w%buf1_ortho(o, k, 2), w%buf2_ortho(o, k, 2), w%tmp_ortho(o, k))
2779 : END IF
2780 1734 : CALL qs_ot_load_local_pair(qs_ot_env%matrix_c0, qs_ot_env%matrix_c0_im, w%c0)
2781 1734 : CALL qs_ot_load_local_pair(qs_ot_env%matrix_r, qs_ot_env%matrix_r_im, w%r)
2782 1734 : CALL qs_ot_load_local_pair(qs_ot_env%matrix_sinp, qs_ot_env%matrix_sinp_im, w%sinp)
2783 1734 : CALL qs_ot_load_local_pair(qs_ot_env%matrix_sc0, qs_ot_env%matrix_sc0_im, w%sc0)
2784 1734 : CALL qs_ot_load_local_pair(qs_ot_env%matrix_os, qs_ot_env%matrix_os_im, w%os)
2785 1734 : CALL copy_dbcsr_to_local(qs_ot_env%matrix_sinp_b, w%sinp_b)
2786 1734 : CALL copy_dbcsr_to_local(qs_ot_env%matrix_cosp_b, w%cosp_b)
2787 1734 : IF (ASSOCIATED(qs_ot_env%preconditioner)) THEN
2788 1734 : CALL qs_ot_load_local_pair(qs_ot_env%matrix_psc0, qs_ot_env%matrix_psc0_im, w%target)
2789 : ELSE
2790 0 : w%target(:, :, :) = w%sc0(:, :, :)
2791 : END IF
2792 : END SUBROUTINE qs_ot_load_local_work
2793 :
2794 : ! **************************************************************************************************
2795 : !> \brief Pack a complex pair, allowing different sparse patterns in its two components.
2796 : !> \param matrix_re real DBCSR matrix
2797 : !> \param matrix_im imaginary DBCSR matrix
2798 : !> \param dense dense real/imaginary pair
2799 : ! **************************************************************************************************
2800 26773 : SUBROUTINE qs_ot_load_local_pair(matrix_re, matrix_im, dense)
2801 : TYPE(dbcsr_type), INTENT(IN) :: matrix_re, matrix_im
2802 : REAL(dp), INTENT(OUT) :: dense(:, :, :)
2803 :
2804 26773 : CALL copy_dbcsr_to_local(matrix_re, dense(:, :, 1))
2805 26773 : CALL copy_dbcsr_to_local(matrix_im, dense(:, :, 2))
2806 26773 : END SUBROUTINE qs_ot_load_local_pair
2807 :
2808 : ! **************************************************************************************************
2809 : !> \brief Multiply dense real/imaginary pairs in the same four-product order as DBCSR.
2810 : !> \param op_a N or C
2811 : !> \param op_b N or C
2812 : !> \param a left operand
2813 : !> \param b right operand
2814 : !> \param c result
2815 : !> \param tmp real workspace, at least as large as one component of C
2816 : !> \param gemm_ctx caller-owned GEMM context
2817 : ! **************************************************************************************************
2818 47149 : SUBROUTINE qs_ot_complex_multiply_pair_local(op_a, op_b, a, b, c, tmp, gemm_ctx)
2819 : CHARACTER(LEN=1), INTENT(IN) :: op_a, op_b
2820 : REAL(dp), CONTIGUOUS, INTENT(IN) :: a(:, :, :), b(:, :, :)
2821 : REAL(dp), CONTIGUOUS, INTENT(INOUT) :: c(:, :, :), tmp(:, :)
2822 : TYPE(local_gemm_ctxt_type) :: gemm_ctx
2823 :
2824 : CALL qs_ot_complex_multiply_parts_local(op_a, op_b, a(:, :, 1), a(:, :, 2), &
2825 47149 : b(:, :, 1), b(:, :, 2), c(:, :, 1), c(:, :, 2), tmp, gemm_ctx)
2826 47149 : END SUBROUTINE qs_ot_complex_multiply_pair_local
2827 :
2828 : ! **************************************************************************************************
2829 : !> \brief Multiply separate real/imaginary arrays, including cached overlap matrices.
2830 : !> \param op_a N or C
2831 : !> \param op_b N or C
2832 : !> \param a_re real left operand
2833 : !> \param a_im imaginary left operand
2834 : !> \param b_re real right operand
2835 : !> \param b_im imaginary right operand
2836 : !> \param c_re real result
2837 : !> \param c_im imaginary result
2838 : !> \param tmp real workspace
2839 : !> \param gemm_ctx caller-owned GEMM context
2840 : ! **************************************************************************************************
2841 51491 : SUBROUTINE qs_ot_complex_multiply_parts_local(op_a, op_b, a_re, a_im, b_re, b_im, &
2842 51491 : c_re, c_im, tmp, gemm_ctx)
2843 : CHARACTER(LEN=1), INTENT(IN) :: op_a, op_b
2844 : REAL(dp), CONTIGUOUS, INTENT(IN) :: a_re(:, :), a_im(:, :), b_re(:, :), &
2845 : b_im(:, :)
2846 : REAL(dp), CONTIGUOUS, INTENT(INOUT) :: c_re(:, :), c_im(:, :), tmp(:, :)
2847 : TYPE(local_gemm_ctxt_type) :: gemm_ctx
2848 :
2849 : CHARACTER(LEN=1) :: trans_a, trans_b
2850 : INTEGER :: k, m, n
2851 : REAL(dp) :: sign_a, sign_b
2852 :
2853 51491 : CPASSERT(op_a == 'N' .OR. op_a == 'C')
2854 51491 : CPASSERT(op_b == 'N' .OR. op_b == 'C')
2855 51491 : trans_a = MERGE('N', 'T', op_a == 'N')
2856 51491 : trans_b = MERGE('N', 'T', op_b == 'N')
2857 51491 : sign_a = MERGE(1.0_dp, -1.0_dp, op_a == 'N')
2858 51491 : sign_b = MERGE(1.0_dp, -1.0_dp, op_b == 'N')
2859 51491 : m = SIZE(c_re, 1)
2860 51491 : n = SIZE(c_re, 2)
2861 51491 : k = MERGE(SIZE(a_re, 2), SIZE(a_re, 1), op_a == 'N')
2862 : CALL gemm_ctx%gemm(trans_a, trans_b, m, n, k, 1.0_dp, a_re, SIZE(a_re, 1), &
2863 51491 : b_re, SIZE(b_re, 1), 0.0_dp, c_re, SIZE(c_re, 1))
2864 : CALL gemm_ctx%gemm(trans_a, trans_b, m, n, k, 1.0_dp, a_im, SIZE(a_im, 1), &
2865 51491 : b_im, SIZE(b_im, 1), 0.0_dp, tmp, SIZE(tmp, 1))
2866 9040231 : c_re(:, :) = c_re(:, :) - sign_a*sign_b*tmp(1:m, 1:n)
2867 : CALL gemm_ctx%gemm(trans_a, trans_b, m, n, k, sign_b, a_re, SIZE(a_re, 1), &
2868 51491 : b_im, SIZE(b_im, 1), 0.0_dp, c_im, SIZE(c_im, 1))
2869 : CALL gemm_ctx%gemm(trans_a, trans_b, m, n, k, sign_a, a_im, SIZE(a_im, 1), &
2870 51491 : b_re, SIZE(b_re, 1), 0.0_dp, tmp, SIZE(tmp, 1))
2871 9040231 : c_im(:, :) = c_im(:, :) + tmp(1:m, 1:n)
2872 51491 : END SUBROUTINE qs_ot_complex_multiply_parts_local
2873 :
2874 : ! **************************************************************************************************
2875 : !> \brief Contract one STRICT term in the local eigenvector basis.
2876 : !> \param w channel matrices and scratch
2877 : !> \param orbital dense X or C0 pair
2878 : !> \param kernel real sinc or cos divided differences
2879 : !> \param RESULT contraction in the eigenvector basis
2880 : ! **************************************************************************************************
2881 3520 : SUBROUTINE qs_ot_strict_frechet_local(w, orbital, kernel, RESULT)
2882 : TYPE(qs_ot_local_work_type) :: w
2883 : REAL(dp), CONTIGUOUS, INTENT(IN) :: orbital(:, :, :), kernel(:, :)
2884 : REAL(dp), CONTIGUOUS, INTENT(INOUT) :: result(:, :, :)
2885 :
2886 3520 : CALL qs_ot_complex_multiply_local('C', 'N', w%hc, orbital, w%buf2, w%buf1(:, :, 1), w%gemm_ctx)
2887 3520 : CALL qs_ot_complex_multiply_local('N', 'N', w%buf2, w%r, w%buf1, w%buf4(:, :, 1), w%gemm_ctx)
2888 3520 : CALL qs_ot_complex_multiply_local('C', 'N', w%r, w%buf1, w%buf2, w%buf4(:, :, 1), w%gemm_ctx)
2889 138416 : RESULT(:, :, 1) = w%buf2(:, :, 1)*kernel(:, :)
2890 138416 : RESULT(:, :, 2) = w%buf2(:, :, 2)*kernel(:, :)
2891 3520 : END SUBROUTINE qs_ot_strict_frechet_local
2892 :
2893 : ! **************************************************************************************************
2894 : !> \brief Evaluate and project one local STRICT gradient without allocating or calling DBCSR.
2895 : !> \param w channel matrices and scratch, with a prepared inverse tangent metric
2896 : ! **************************************************************************************************
2897 1760 : SUBROUTINE qs_ot_strict_gradient_local(w)
2898 : TYPE(qs_ot_local_work_type) :: w
2899 :
2900 : INTEGER :: i, j
2901 : REAL(dp) :: im, re
2902 :
2903 : ! Direct term, followed by the X*sinc and C0*cos Frechet contributions.
2904 1760 : CALL qs_ot_complex_multiply_local('N', 'N', w%hc, w%sinp, w%out, w%tmp, w%gemm_ctx)
2905 1760 : CALL qs_ot_strict_frechet_local(w, w%x, w%sinp_b, w%buf3)
2906 1760 : CALL qs_ot_strict_frechet_local(w, w%c0, w%cosp_b, w%buf4)
2907 140176 : w%buf3(:, :, :) = w%buf3(:, :, :) + w%buf4(:, :, :)
2908 1760 : CALL qs_ot_complex_multiply_local('N', 'N', w%r, w%buf3, w%buf1, w%buf2(:, :, 1), w%gemm_ctx)
2909 1760 : CALL qs_ot_complex_multiply_local('N', 'C', w%buf1, w%r, w%buf3, w%buf2(:, :, 1), w%gemm_ctx)
2910 :
2911 : ! Merge each conjugate pair in place, without a transpose temporary.
2912 10056 : DO j = 1, w%k
2913 8296 : w%buf3(j, j, 1) = w%buf3(j, j, 1) + w%buf3(j, j, 1)
2914 8296 : w%buf3(j, j, 2) = 0.0_dp
2915 35484 : DO i = j + 1, w%k
2916 25428 : re = w%buf3(i, j, 1) + w%buf3(j, i, 1)
2917 25428 : im = w%buf3(i, j, 2) - w%buf3(j, i, 2)
2918 25428 : w%buf3(i, j, 1) = re
2919 25428 : w%buf3(j, i, 1) = re
2920 25428 : w%buf3(i, j, 2) = im
2921 33724 : w%buf3(j, i, 2) = -im
2922 : END DO
2923 : END DO
2924 1760 : CALL qs_ot_complex_multiply_local('N', 'N', w%sx, w%buf3, w%buf_nk, w%tmp, w%gemm_ctx)
2925 478912 : w%out(:, :, :) = w%out(:, :, :) + w%buf_nk(:, :, :)
2926 :
2927 : ! Tangent projection with the same preconditioned reference and inverse metric.
2928 1760 : CALL qs_ot_complex_multiply_local('C', 'N', w%target, w%out, w%buf1_ortho, w%tmp_ortho, w%gemm_ctx)
2929 1760 : CALL qs_ot_complex_multiply_local('N', 'N', w%os, w%buf1_ortho, w%buf2_ortho, w%tmp_ortho, w%gemm_ctx)
2930 1760 : CALL qs_ot_complex_multiply_local('N', 'N', w%sc0, w%buf2_ortho, w%buf_nk, w%tmp, w%gemm_ctx)
2931 478912 : w%out(:, :, :) = w%out(:, :, :) - w%buf_nk(:, :, :)
2932 1760 : END SUBROUTINE qs_ot_strict_gradient_local
2933 :
2934 1488 : END MODULE qs_ot_kp
|
