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