(git:f2099e5)
Loading...
Searching...
No Matches
qs_ot.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 orbital transformations
10!> \par History
11!> Added Taylor expansion based computation of the matrix functions (01.2004)
12!> added additional rotation variables for non-equivalent occupied orbs (08.2004)
13!> \author Joost VandeVondele (06.2002)
14! **************************************************************************************************
15MODULE qs_ot
16 USE arnoldi_api, ONLY: arnoldi_extremal
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_get_occupation, dbcsr_init_p, &
20 dbcsr_iterator_blocks_left, dbcsr_iterator_next_block, dbcsr_iterator_start, &
21 dbcsr_iterator_stop, dbcsr_iterator_type, dbcsr_multiply, dbcsr_release, dbcsr_release_p, &
22 dbcsr_reserve_blocks, dbcsr_scale, dbcsr_set, dbcsr_transposed, dbcsr_type, &
23 dbcsr_type_no_symmetry
24 USE cp_dbcsr_cholesky, ONLY: cp_dbcsr_cholesky_decompose,&
25 cp_dbcsr_cholesky_invert,&
26 cp_dbcsr_cholesky_restore
27 USE cp_dbcsr_contrib, ONLY: dbcsr_add_on_diag,&
28 dbcsr_frobenius_norm,&
29 dbcsr_gershgorin_norm,&
30 dbcsr_hadamard_product,&
31 dbcsr_scale_by_vector
32 USE cp_dbcsr_diag, ONLY: cp_dbcsr_heevd,&
33 cp_dbcsr_syevd
34 USE kinds, ONLY: dp
35 USE mathlib, ONLY: diag_complex,&
36 diamat_all,&
37 gemm_square
38 USE message_passing, ONLY: mp_comm_type
39 USE preconditioner, ONLY: apply_preconditioner
40 USE preconditioner_types, ONLY: preconditioner_type
41 USE qs_ot_types, ONLY: qs_ot_type
42#include "./base/base_uses.f90"
43
44 IMPLICIT NONE
45 PRIVATE
46
47 PUBLIC :: qs_ot_get_p
48 PUBLIC :: qs_ot_get_p_complex
49 PUBLIC :: qs_ot_get_orbitals
50 PUBLIC :: qs_ot_get_orbitals_complex
51 PUBLIC :: qs_ot_get_derivative
52 PUBLIC :: qs_ot_get_derivative_complex
53 PUBLIC :: qs_ot_prepare_complex_tangent_metric
54 PUBLIC :: qs_ot_get_orbitals_ref
55 PUBLIC :: qs_ot_get_orbitals_ref_complex
56 PUBLIC :: qs_ot_get_derivative_ref
57 PUBLIC :: qs_ot_get_derivative_ref_complex
58 PUBLIC :: qs_ot_apply_complex_frechet_dbcsr
59 PUBLIC :: qs_ot_antihermitian_spectral_norm
60 PUBLIC :: qs_ot_complex_exp_frechet_kernel
61 PUBLIC :: qs_ot_density_secant_hessian
62 PUBLIC :: qs_ot_density_secant_orbital_overlaps
63 PUBLIC :: qs_ot_density_secant_projected_hessian
64 PUBLIC :: qs_ot_density_tangent
65 PUBLIC :: qs_ot_fixed_n_energy_gradient
66 PUBLIC :: qs_ot_fixed_n_energy_hessian
67 PUBLIC :: qs_ot_fixed_n_multigroup_schur_block
68 PUBLIC :: qs_ot_fixed_n_schur_block
69 PUBLIC :: qs_ot_fixed_n_projector_frechet
70 PUBLIC :: qs_ot_fixed_n_response_mu_shift
71 PUBLIC :: qs_ot_finite_rotation_response
72 PUBLIC :: qs_ot_projected_response_update
73 PUBLIC :: qs_ot_symmetric_sr1_update
74 PUBLIC :: qs_ot_symmetric_abs_solve
75 PUBLIC :: qs_ot_generate_rotation
76 PUBLIC :: qs_ot_generate_rotation_complex
77 PUBLIC :: qs_ot_rot_mat_derivative
78 PUBLIC :: qs_ot_rot_mat_derivative_complex
79 PUBLIC :: qs_ot_new_preconditioner
80 PRIVATE :: qs_ot_p2m_diag
81 PRIVATE :: qs_ot_p2m_diag_complex
82 PRIVATE :: qs_ot_complex_multiply
83 PRIVATE :: qs_ot_sinc
84 PRIVATE :: qs_ot_ref_poly
85 PRIVATE :: qs_ot_ref_chol
86 PRIVATE :: qs_ot_ref_lwdn
87 PRIVATE :: qs_ot_ref_decide
88 PRIVATE :: qs_ot_ref_update
89 PRIVATE :: qs_ot_refine
90 PRIVATE :: qs_ot_on_the_fly_localize
91
92 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_ot'
93
94CONTAINS
95
96! **************************************************************************************************
97!> \brief spectral norm of a dense anti-Hermitian rotation generator
98!> \param rotation_generator anti-Hermitian generator
99!> \return largest absolute eigenvalue of i times the generator
100! **************************************************************************************************
101 FUNCTION qs_ot_antihermitian_spectral_norm(rotation_generator) RESULT(norm)
102 COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(IN) :: rotation_generator
103 REAL(kind=dp) :: norm
104
105 COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: eigenvectors
106 INTEGER :: n
107 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: eigenvalues
108
109 n = SIZE(rotation_generator, 1)
110 cpassert(SIZE(rotation_generator, 2) == n)
111 norm = 0.0_dp
112 IF (n == 0) RETURN
113 ALLOCATE (eigenvectors(n, n), eigenvalues(n))
114 CALL diag_complex(cmplx(0.0_dp, 1.0_dp, kind=dp)*rotation_generator, &
115 eigenvectors, eigenvalues)
116 norm = maxval(abs(eigenvalues))
117 DEALLOCATE (eigenvalues, eigenvectors)
118
119 END FUNCTION qs_ot_antihermitian_spectral_norm
120
121! **************************************************************************************************
122!> \brief chemical-potential response for one fixed-electron-number group
123!> \param weighted_energy_response sum_i chi_i de_i over the perturbed local channels
124!> \param local_curvature_sum local sum_i chi_i, used as a serial fallback
125!> \param fixed_n_curvature_sum global sum_i chi_i for the complete fixed-N group
126!> \return first-order chemical-potential shift
127! **************************************************************************************************
128 PURE FUNCTION qs_ot_fixed_n_response_mu_shift( &
129 weighted_energy_response, local_curvature_sum, fixed_n_curvature_sum) RESULT(mu_shift)
130 REAL(kind=dp), INTENT(IN) :: weighted_energy_response, &
131 local_curvature_sum, &
132 fixed_n_curvature_sum
133 REAL(kind=dp) :: mu_shift
134
135 REAL(kind=dp) :: denominator
136
137 denominator = fixed_n_curvature_sum
138 IF (abs(denominator) <= epsilon(denominator)) denominator = local_curvature_sum
139 mu_shift = 0.0_dp
140 IF (abs(denominator) > epsilon(denominator)) mu_shift = weighted_energy_response/denominator
141
142 END FUNCTION qs_ot_fixed_n_response_mu_shift
143
144! **************************************************************************************************
145!> \brief fixed-N Mermin gradient in auxiliary-energy coordinates
146!> \param rayleigh_energy diagonal expectation values of the current Hamiltonian
147!> \param energy_coordinate auxiliary band energies controlling the occupations
148!> \param response_weight signed weighted occupation responses chi_i
149!> \param fixed_n_weight_sum global sum_i chi_i for the fixed-N group
150!> \param fixed_n_weighted_residual global sum_i chi_i (h_i-e_i)
151!> \param gradient projected fixed-N gradient
152! **************************************************************************************************
153 PURE SUBROUTINE qs_ot_fixed_n_energy_gradient( &
154 rayleigh_energy, energy_coordinate, response_weight, fixed_n_weight_sum, &
155 fixed_n_weighted_residual, gradient)
156 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: rayleigh_energy, energy_coordinate, &
157 response_weight
158 REAL(kind=dp), INTENT(IN) :: fixed_n_weight_sum, &
159 fixed_n_weighted_residual
160 REAL(kind=dp), DIMENSION(:), INTENT(OUT) :: gradient
161
162 REAL(kind=dp) :: fixed_n_mean
163
164 fixed_n_mean = 0.0_dp
165 IF (abs(fixed_n_weight_sum) > epsilon(fixed_n_weight_sum)) THEN
166 fixed_n_mean = fixed_n_weighted_residual/fixed_n_weight_sum
167 END IF
168 gradient(:) = response_weight(:)* &
169 (fixed_n_mean - (rayleigh_energy(:) - energy_coordinate(:)))
170
171 END SUBROUTINE qs_ot_fixed_n_energy_gradient
172
173! **************************************************************************************************
174!> \brief dense fixed-N occupation Hessian in auxiliary-energy coordinates
175!>
176!> H = diag(chi) - chi chi^T / sum(chi) is symmetric and has the constant-energy gauge as
177!> an exact null vector. It can be indefinite for non-monotone smearing distributions.
178!> \param response_weight signed weighted occupation responses chi_i
179!> \param fixed_n_weight_sum global sum_i chi_i for the fixed-N group
180!> \param hessian projected fixed-N Hessian
181! **************************************************************************************************
182 PURE SUBROUTINE qs_ot_fixed_n_energy_hessian(response_weight, fixed_n_weight_sum, hessian)
183 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: response_weight
184 REAL(kind=dp), INTENT(IN) :: fixed_n_weight_sum
185 REAL(kind=dp), DIMENSION(:, :), INTENT(OUT) :: hessian
186
187 INTEGER :: i, j, n
188
189 n = SIZE(response_weight)
190 hessian(:, :) = 0.0_dp
191 DO i = 1, n
192 hessian(i, i) = response_weight(i)
193 END DO
194 IF (abs(fixed_n_weight_sum) > epsilon(fixed_n_weight_sum)) THEN
195 DO j = 1, n
196 DO i = 1, n
197 hessian(i, j) = hessian(i, j) - &
198 response_weight(i)*response_weight(j)/fixed_n_weight_sum
199 END DO
200 END DO
201 END IF
202 hessian(:, :) = 0.5_dp*(hessian + transpose(hessian))
203
204 END SUBROUTINE qs_ot_fixed_n_energy_hessian
205
206! **************************************************************************************************
207!> \brief local block of the fixed-N rotation/energy Schur complement
208!>
209!> For C = D - chi chi^T / sum(chi), elimination of the auxiliary-energy block gives
210!>
211!> S = A - R^T C R
212!> = (A - R^T D R) + v v^T / sum(chi), v = R^T chi.
213!>
214!> This routine builds the channel-local terms. The final rank-one term is deliberately left
215!> separate so spin/k-point channels can be coupled without assembling a global dense
216!> rotation Hessian.
217!> \param rotation_hessian fixed-occupation rotation Hessian A
218!> \param rayleigh_response derivative R of the Rayleigh energies with respect to rotations
219!> \param response_weight local signed occupation responses chi
220!> \param rotation_gradient physical rotation gradient
221!> \param energy_gradient physical auxiliary-energy gradient
222!> \param schur_block local block A - R^T D R
223!> \param coupling_vector local part of v = R^T chi
224!> \param schur_rhs local right-hand side g_x + R^T g_e
225! **************************************************************************************************
226 SUBROUTINE qs_ot_fixed_n_schur_block( &
227 rotation_hessian, rayleigh_response, response_weight, rotation_gradient, energy_gradient, &
228 schur_block, coupling_vector, schur_rhs)
229 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: rotation_hessian, rayleigh_response
230 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: response_weight, rotation_gradient, &
231 energy_gradient
232 REAL(kind=dp), DIMENSION(:, :), INTENT(OUT) :: schur_block
233 REAL(kind=dp), DIMENSION(:), INTENT(OUT) :: coupling_vector, schur_rhs
234
235 INTEGER :: nenergy
236 INTEGER, ALLOCATABLE, DIMENSION(:) :: response_group
237 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: coupling_matrix
238
239 nenergy = SIZE(response_weight)
240 ALLOCATE (response_group(nenergy), coupling_matrix(SIZE(coupling_vector), 1))
241 response_group(:) = 1
242 CALL qs_ot_fixed_n_multigroup_schur_block( &
243 rotation_hessian, rayleigh_response, response_weight, response_group, &
244 rotation_gradient, energy_gradient, schur_block, coupling_matrix, schur_rhs)
245 coupling_vector(:) = coupling_matrix(:, 1)
246 DEALLOCATE (response_group, coupling_matrix)
247
248 END SUBROUTINE qs_ot_fixed_n_schur_block
249
250! **************************************************************************************************
251!> \brief eliminate spin-resolved auxiliary energies while retaining every fixed-N constraint
252!>
253!> Each energy variable belongs to one particle-number group. The diagonal occupation
254!> response is eliminated locally, while one coupling column per group is retained for the
255!> subsequent global chemical-potential projection. A single group reduces exactly to
256!> qs_ot_fixed_n_schur_block.
257!> \param rotation_hessian finite orbital-rotation Hessian
258!> \param rayleigh_response derivative of all spin-resolved Rayleigh energies
259!> \param response_weight positive occupation-response weights
260!> \param response_group fixed-N group of every energy variable
261!> \param rotation_gradient orbital-rotation gradient
262!> \param energy_gradient auxiliary-energy gradient
263!> \param schur_block energy-eliminated rotation block
264!> \param coupling_matrix retained chemical-potential coupling columns
265!> \param schur_rhs energy-eliminated rotation right-hand side
266! **************************************************************************************************
267 SUBROUTINE qs_ot_fixed_n_multigroup_schur_block( &
268 rotation_hessian, rayleigh_response, response_weight, response_group, &
269 rotation_gradient, energy_gradient, schur_block, coupling_matrix, schur_rhs)
270 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: rotation_hessian, rayleigh_response
271 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: response_weight
272 INTEGER, DIMENSION(:), INTENT(IN) :: response_group
273 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: rotation_gradient, energy_gradient
274 REAL(kind=dp), DIMENSION(:, :), INTENT(OUT) :: schur_block, coupling_matrix
275 REAL(kind=dp), DIMENSION(:), INTENT(OUT) :: schur_rhs
276
277 INTEGER :: group, i, nenergy, ngroups, nrotation
278 REAL(kind=dp) :: weight
279
280 nenergy = SIZE(response_weight)
281 nrotation = SIZE(rotation_gradient)
282 ngroups = SIZE(coupling_matrix, 2)
283 cpassert(ngroups > 0)
284 cpassert(all(shape(rotation_hessian) == [nrotation, nrotation]))
285 cpassert(all(shape(rayleigh_response) == [nenergy, nrotation]))
286 cpassert(SIZE(energy_gradient) == nenergy)
287 cpassert(SIZE(response_group) == nenergy)
288 cpassert(all(response_group >= 1 .AND. response_group <= ngroups))
289 cpassert(all(shape(schur_block) == [nrotation, nrotation]))
290 cpassert(SIZE(coupling_matrix, 1) == nrotation)
291 cpassert(SIZE(schur_rhs) == nrotation)
292
293 schur_block(:, :) = rotation_hessian(:, :)
294 coupling_matrix(:, :) = 0.0_dp
295 DO i = 1, nenergy
296 weight = response_weight(i)
297 schur_block(:, :) = schur_block(:, :) - weight* &
298 spread(rayleigh_response(i, :), dim=2, ncopies=nrotation)* &
299 spread(rayleigh_response(i, :), dim=1, ncopies=nrotation)
300 group = response_group(i)
301 coupling_matrix(:, group) = coupling_matrix(:, group) + &
302 weight*rayleigh_response(i, :)
303 END DO
304 schur_block(:, :) = 0.5_dp*(schur_block + transpose(schur_block))
305 schur_rhs(:) = rotation_gradient + &
306 matmul(transpose(rayleigh_response), energy_gradient)
307
308 END SUBROUTINE qs_ot_fixed_n_multigroup_schur_block
309
310! **************************************************************************************************
311!> \brief apply a positive spectral inverse of a real symmetric response matrix
312!>
313!> The magnitude of every resolved eigenmode is retained, including modes with negative
314!> physical curvature. Replacing lambda by abs(lambda) gives a descent metric without the
315!> loss of response information caused by discarding the negative subspace. Unresolved
316!> null modes are projected out instead of being amplified by an artificial eigenvalue floor.
317!> \param matrix real symmetric response matrix
318!> \param rhs one or more right-hand sides
319!> \param solution spectral-absolute inverse applied to rhs
320!> \param valid whether finite input and output were obtained
321!> \param relative_floor optional relative eigenvalue floor
322! **************************************************************************************************
323 SUBROUTINE qs_ot_symmetric_abs_solve(matrix, rhs, solution, valid, relative_floor)
324 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: matrix, rhs
325 REAL(kind=dp), DIMENSION(:, :), INTENT(OUT) :: solution
326 LOGICAL, INTENT(OUT) :: valid
327 REAL(kind=dp), INTENT(IN), OPTIONAL :: relative_floor
328
329 INTEGER :: i, n, nresolved
330 REAL(kind=dp) :: eigenvalue_floor, relative_floor_eff, &
331 scale
332 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: eigenvalues
333 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: eigenvectors, work
334
335 n = SIZE(matrix, 1)
336 cpassert(all(shape(matrix) == [n, n]))
337 cpassert(SIZE(rhs, 1) == n)
338 cpassert(all(shape(solution) == shape(rhs)))
339 valid = all(matrix == matrix) .AND. all(rhs == rhs)
340 IF (.NOT. valid) THEN
341 solution(:, :) = 0.0_dp
342 RETURN
343 END IF
344
345 relative_floor_eff = sqrt(epsilon(1.0_dp))
346 IF (PRESENT(relative_floor)) relative_floor_eff = max(relative_floor, relative_floor_eff)
347 ALLOCATE (eigenvalues(n), eigenvectors(n, n), work(n, SIZE(rhs, 2)))
348 eigenvectors(:, :) = 0.5_dp*(matrix + transpose(matrix))
349 CALL diamat_all(eigenvectors, eigenvalues)
350 scale = max(1.0_dp, maxval(abs(eigenvalues)))
351 eigenvalue_floor = relative_floor_eff*scale
352 work(:, :) = matmul(transpose(eigenvectors), rhs)
353 nresolved = 0
354 DO i = 1, n
355 IF (abs(eigenvalues(i)) > eigenvalue_floor) THEN
356 work(i, :) = work(i, :)/abs(eigenvalues(i))
357 nresolved = nresolved + 1
358 ELSE
359 work(i, :) = 0.0_dp
360 END IF
361 END DO
362 solution(:, :) = matmul(eigenvectors, work)
363 valid = nresolved > 0 .AND. all(solution == solution)
364 IF (.NOT. valid) solution(:, :) = 0.0_dp
365 DEALLOCATE (eigenvalues, eigenvectors, work)
366
367 END SUBROUTINE qs_ot_symmetric_abs_solve
368
369! **************************************************************************************************
370!> \brief update a baseline response direction in a small positive physical-response subspace
371!>
372!> The first basis mode is the baseline response direction. With B0 the projected frozen-H
373!> Hessian and K the accepted physical correction, this routine solves
374!>
375!> (B0 + K) c = g_Q.
376!>
377!> The optional projected physical gradient supplies g_Q. Without it, g_Q=B0*e1, so c=e1
378!> exactly when K=0. An indefinite or unresolved total projected Hessian is rejected instead
379!> of turning its negative modes into an unrelated active direction.
380!> \param reference_hessian projected frozen-H Hessian B0
381!> \param response_correction projected physical response K
382!> \param coefficients response coefficients c in the supplied basis
383!> \param valid whether a finite, positive, sufficiently resolved solve was obtained
384!> \param projected_gradient optional physical gradient projected onto the supplied basis
385!> \param relative_floor optional relative positive-eigenvalue floor
386! **************************************************************************************************
387 SUBROUTINE qs_ot_projected_response_update( &
388 reference_hessian, response_correction, coefficients, valid, projected_gradient, relative_floor)
389
390 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: reference_hessian, response_correction
391 REAL(kind=dp), DIMENSION(:), INTENT(OUT) :: coefficients
392 LOGICAL, INTENT(OUT) :: valid
393 REAL(kind=dp), DIMENSION(:), INTENT(IN), OPTIONAL :: projected_gradient
394 REAL(kind=dp), INTENT(IN), OPTIONAL :: relative_floor
395
396 INTEGER :: i, n
397 REAL(kind=dp) :: eigenvalue_floor, relative_floor_eff, &
398 scale
399 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: eigenvalues, rhs, work
400 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: eigenvectors
401
402 n = SIZE(reference_hessian, 1)
403 cpassert(n > 0)
404 cpassert(all(shape(reference_hessian) == [n, n]))
405 cpassert(all(shape(response_correction) == [n, n]))
406 cpassert(SIZE(coefficients) == n)
407 IF (PRESENT(projected_gradient)) THEN
408 cpassert(SIZE(projected_gradient) == n)
409 END IF
410 coefficients(:) = 0.0_dp
411 valid = all(reference_hessian == reference_hessian) .AND. &
412 all(response_correction == response_correction)
413 IF (PRESENT(projected_gradient)) THEN
414 valid = valid .AND. all(projected_gradient == projected_gradient)
415 END IF
416 IF (.NOT. valid) RETURN
417
418 relative_floor_eff = 1.0e-4_dp
419 IF (PRESENT(relative_floor)) relative_floor_eff = &
420 max(relative_floor, sqrt(epsilon(1.0_dp)))
421 ALLOCATE (eigenvalues(n), eigenvectors(n, n), rhs(n), work(n))
422 eigenvectors(:, :) = 0.5_dp* &
423 (reference_hessian + response_correction + &
424 transpose(reference_hessian + response_correction))
425 CALL diamat_all(eigenvectors, eigenvalues)
426 scale = max(maxval(abs(eigenvalues)), maxval(abs(reference_hessian)))
427 IF (scale <= tiny(scale)) THEN
428 valid = .false.
429 coefficients(:) = 0.0_dp
430 DEALLOCATE (eigenvalues, eigenvectors, rhs, work)
431 RETURN
432 END IF
433 eigenvalue_floor = relative_floor_eff*scale
434 valid = all(eigenvalues > eigenvalue_floor)
435 IF (valid) THEN
436 rhs(:) = reference_hessian(:, 1)
437 IF (PRESENT(projected_gradient)) rhs(:) = projected_gradient
438 work(:) = matmul(transpose(eigenvectors), rhs)
439 DO i = 1, n
440 work(i) = work(i)/eigenvalues(i)
441 END DO
442 coefficients(:) = matmul(eigenvectors, work)
443 valid = all(coefficients == coefficients) .AND. &
444 all(abs(coefficients) <= huge(1.0_dp))
445 END IF
446 IF (.NOT. valid) coefficients(:) = 0.0_dp
447 DEALLOCATE (eigenvalues, eigenvectors, rhs, work)
448
449 END SUBROUTINE qs_ot_projected_response_update
450
451! **************************************************************************************************
452!> \brief add one accepted symmetric response secant to a reference Hessian
453!>
454!> With r=y-B0*s, the symmetric-rank-one update B=B0+r*r^T/(r^T*s) satisfies B*s=y
455!> exactly. The signed denominator is retained because a self-consistent Hxc response can
456!> be indefinite. Nearly orthogonal residuals are rejected instead of manufacturing a
457!> large unresolved mode.
458!> \param matrix reference symmetric Hessian B0
459!> \param step accepted displacement s
460!> \param response measured gradient response y
461!> \param updated_matrix symmetric secant Hessian B
462!> \param valid whether a resolved finite update was constructed
463!> \param relative_tolerance optional SR1 denominator acceptance threshold
464! **************************************************************************************************
465 PURE SUBROUTINE qs_ot_symmetric_sr1_update( &
466 matrix, step, response, updated_matrix, valid, relative_tolerance)
467 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: matrix
468 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: step, response
469 REAL(kind=dp), DIMENSION(:, :), INTENT(OUT) :: updated_matrix
470 LOGICAL, INTENT(OUT) :: valid
471 REAL(kind=dp), INTENT(IN), OPTIONAL :: relative_tolerance
472
473 INTEGER :: n
474 REAL(kind=dp) :: denominator, residual_norm, step_norm, &
475 threshold, tolerance
476 REAL(kind=dp), DIMENSION(SIZE(step)) :: residual
477
478 n = SIZE(step)
479 updated_matrix(:, :) = 0.0_dp
480 valid = SIZE(response) == n .AND. all(shape(matrix) == [n, n]) .AND. &
481 all(shape(updated_matrix) == [n, n])
482 IF (.NOT. valid) RETURN
483 valid = all(matrix == matrix) .AND. all(step == step) .AND. all(response == response)
484 IF (.NOT. valid) RETURN
485
486 updated_matrix(:, :) = 0.5_dp*(matrix + transpose(matrix))
487 residual(:) = response - matmul(updated_matrix, step)
488 denominator = dot_product(residual, step)
489 residual_norm = sqrt(max(0.0_dp, dot_product(residual, residual)))
490 step_norm = sqrt(max(0.0_dp, dot_product(step, step)))
491 tolerance = sqrt(epsilon(1.0_dp))
492 IF (PRESENT(relative_tolerance)) tolerance = max(tolerance, relative_tolerance)
493 threshold = tolerance*residual_norm*step_norm
494 valid = residual_norm > tiny(residual_norm) .AND. step_norm > tiny(step_norm) .AND. &
495 abs(denominator) > threshold
496 IF (.NOT. valid) RETURN
497
498 updated_matrix(:, :) = updated_matrix + &
499 spread(residual, dim=2, ncopies=n)*spread(residual, dim=1, ncopies=n)/denominator
500 updated_matrix(:, :) = 0.5_dp*(updated_matrix + transpose(updated_matrix))
501 valid = all(updated_matrix == updated_matrix)
502
503 END SUBROUTINE qs_ot_symmetric_sr1_update
504
505! **************************************************************************************************
506!> \brief project a self-adjoint density/Hamiltonian secant onto density-response modes
507!>
508!> For an accepted Hermitian density change S and the corresponding self-consistent
509!> Hamiltonian change Y, the minimum-Frobenius-norm self-adjoint response satisfying
510!> K*S=Y is
511!>
512!> K = (Y<S,.> + S<Y,.>)/<S,S> - <S,Y>S<S,.>/<S,S>**2.
513!>
514!> The returned matrix is <B_q,K*B_r> for the supplied Hermitian density modes B_r. Its
515!> density-space construction is invariant under a common complex similarity transform.
516!> \param density_step accepted density-matrix change S
517!> \param hamiltonian_step accepted self-consistent Hamiltonian change Y
518!> \param density_modes density derivatives B_r of the coupled minimizer variables
519!> \param correction projected symmetric Hxc response
520!> \param valid whether a finite nonzero density secant was available
521!> \param density_norm_sq optional <S,S>
522!> \param response_work optional <S,Y>
523! **************************************************************************************************
524 SUBROUTINE qs_ot_density_secant_hessian( &
525 density_step, hamiltonian_step, density_modes, correction, valid, &
526 density_norm_sq, response_work)
527
528 COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(IN) :: density_step, hamiltonian_step
529 COMPLEX(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: density_modes
530 REAL(kind=dp), DIMENSION(:, :), INTENT(OUT) :: correction
531 LOGICAL, INTENT(OUT) :: valid
532 REAL(kind=dp), INTENT(OUT), OPTIONAL :: density_norm_sq, response_work
533
534 INTEGER :: i, j, mode, n, nmode
535 REAL(kind=dp) :: density_norm, density_response_work, &
536 scale
537 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: density_overlap, response_overlap
538
539 n = SIZE(density_step, 1)
540 nmode = SIZE(density_modes, 3)
541 cpassert(n > 0)
542 cpassert(SIZE(density_step, 2) == n)
543 cpassert(all(shape(hamiltonian_step) == [n, n]))
544 cpassert(SIZE(density_modes, 1) == n)
545 cpassert(SIZE(density_modes, 2) == n)
546 cpassert(all(shape(correction) == [nmode, nmode]))
547
548 correction(:, :) = 0.0_dp
549 valid = .false.
550 density_norm = 0.0_dp
551 density_response_work = 0.0_dp
552 DO j = 1, n
553 DO i = 1, n
554 density_norm = density_norm + &
555 REAL(conjg(density_step(i, j))*density_step(i, j), kind=dp)
556 density_response_work = density_response_work + &
557 REAL(conjg(density_step(i, j))*hamiltonian_step(i, j), kind=dp)
558 END DO
559 END DO
560 IF (PRESENT(density_norm_sq)) density_norm_sq = density_norm
561 IF (PRESENT(response_work)) response_work = density_response_work
562 scale = maxval(abs(density_step))
563 IF (scale <= tiny(1.0_dp)) RETURN
564 IF (density_norm <= 64.0_dp*epsilon(1.0_dp)*scale*scale .OR. nmode <= 0) RETURN
565
566 ALLOCATE (density_overlap(nmode), response_overlap(nmode))
567 DO mode = 1, nmode
568 density_overlap(mode) = sum(real(conjg(density_step)*density_modes(:, :, mode), kind=dp))
569 response_overlap(mode) = sum(real(conjg(hamiltonian_step)*density_modes(:, :, mode), kind=dp))
570 END DO
571 CALL qs_ot_density_secant_projected_hessian( &
572 density_norm, density_response_work, density_overlap, response_overlap, correction, valid)
573 DEALLOCATE (density_overlap, response_overlap)
574
575 END SUBROUTINE qs_ot_density_secant_hessian
576
577! **************************************************************************************************
578!> \brief form a projected self-adjoint Hxc response from distributed density-space overlaps
579!> \param density_norm_sq <S,S>
580!> \param response_work <S,Y>
581!> \param density_overlap <S,B_r>
582!> \param response_overlap <Y,B_r>
583!> \param correction projected symmetric Hxc response <B_q,K*B_r>
584!> \param valid whether finite nonzero secant data were available
585!> \param secant_mode optional mode representing the accepted full density secant divided by its
586!> line-search position
587!> \param secant_position signed line-search position of the accepted full density secant
588! **************************************************************************************************
589 SUBROUTINE qs_ot_density_secant_projected_hessian( &
590 density_norm_sq, response_work, density_overlap, response_overlap, correction, valid, &
591 secant_mode, secant_position)
592
593 REAL(kind=dp), INTENT(IN) :: density_norm_sq, response_work
594 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: density_overlap, response_overlap
595 REAL(kind=dp), DIMENSION(:, :), INTENT(OUT) :: correction
596 LOGICAL, INTENT(OUT) :: valid
597 INTEGER, INTENT(IN), OPTIONAL :: secant_mode
598 REAL(kind=dp), INTENT(IN), OPTIONAL :: secant_position
599
600 INTEGER :: i, j, nmode
601 REAL(kind=dp) :: inverse_density_norm
602 REAL(kind=dp), DIMENSION(SIZE(density_overlap)) :: projected_density_overlap, &
603 projected_response_overlap
604
605 nmode = SIZE(density_overlap)
606 cpassert(SIZE(response_overlap) == nmode)
607 cpassert(all(shape(correction) == [nmode, nmode]))
608
609 correction(:, :) = 0.0_dp
610 valid = .false.
611 IF (nmode <= 0 .OR. density_norm_sq <= tiny(1.0_dp)) RETURN
612 IF (density_norm_sq /= density_norm_sq .OR. response_work /= response_work) RETURN
613 IF (abs(density_norm_sq) > huge(1.0_dp) .OR. abs(response_work) > huge(1.0_dp)) RETURN
614 IF (any(density_overlap /= density_overlap) .OR. any(response_overlap /= response_overlap)) RETURN
615 IF (any(abs(density_overlap) > huge(1.0_dp)) .OR. &
616 any(abs(response_overlap) > huge(1.0_dp))) RETURN
617 IF (PRESENT(secant_mode) .NEQV. PRESENT(secant_position)) RETURN
618
619 projected_density_overlap(:) = density_overlap
620 projected_response_overlap(:) = response_overlap
621 IF (PRESENT(secant_mode)) THEN
622 IF (secant_mode < 1 .OR. secant_mode > nmode) RETURN
623 IF (secant_position /= secant_position .OR. &
624 abs(secant_position) > huge(1.0_dp) .OR. &
625 abs(secant_position) <= sqrt(epsilon(1.0_dp))) RETURN
626 ! The complete accepted density change can contain REF-orbital motion that is absent from
627 ! a reduced rotation/occupation tangent. S/alpha supplies its exact projected secant mode.
628 projected_density_overlap(secant_mode) = density_norm_sq/secant_position
629 projected_response_overlap(secant_mode) = response_work/secant_position
630 END IF
631
632 inverse_density_norm = 1.0_dp/density_norm_sq
633 DO j = 1, nmode
634 DO i = 1, nmode
635 correction(i, j) = inverse_density_norm* &
636 (projected_response_overlap(i)*projected_density_overlap(j) + &
637 projected_density_overlap(i)*projected_response_overlap(j) - &
638 response_work*inverse_density_norm* &
639 projected_density_overlap(i)*projected_density_overlap(j))
640 END DO
641 END DO
642 correction(:, :) = 0.5_dp*(correction + transpose(correction))
643 valid = all(correction == correction) .AND. all(abs(correction) <= huge(1.0_dp))
644 IF (.NOT. valid) correction(:, :) = 0.0_dp
645
646 END SUBROUTINE qs_ot_density_secant_projected_hessian
647
648! **************************************************************************************************
649!> \brief finite-chart density tangent for coupled complex rotations and fixed-N occupations
650!>
651!> The rotation contribution differentiates
652!>
653!> exp(X) diag(w_k f) exp(X)^H
654!>
655!> along an anti-Hermitian packed direction. The supplied weighted occupation response is
656!> added in the same chart, and the result is returned in the current physical orbital basis.
657!> \param rotation_generator current anti-Hermitian REF generator X
658!> \param occupation current occupations f
659!> \param kpoint_weight irreducible K-point weight w_k
660!> \param rotation_step interleaved real/imaginary anti-Hermitian direction
661!> \param weighted_occupation_step derivative of w_k*f, including the fixed-N mu response
662!> \param density_tangent Hermitian tangent in the current physical orbital basis
663!> \param difference_step optional finite-chart central-difference step
664! **************************************************************************************************
665 SUBROUTINE qs_ot_density_tangent( &
666 rotation_generator, occupation, kpoint_weight, rotation_step, &
667 weighted_occupation_step, density_tangent, difference_step)
668
669 COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(IN) :: rotation_generator
670 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: occupation
671 REAL(kind=dp), INTENT(IN) :: kpoint_weight
672 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: rotation_step, weighted_occupation_step
673 COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: density_tangent
674 REAL(kind=dp), INTENT(IN), OPTIONAL :: difference_step
675
676 COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: direction, generator_minus, &
677 generator_plus, mode_ref, rotation, rotation_minus, rotation_plus, weighted_rotation
678 INTEGER :: i, j, n, nrotation, r
679 REAL(kind=dp) :: step, step_scale
680
681 n = SIZE(rotation_generator, 1)
682 nrotation = n*(n - 1)
683 cpassert(n > 0)
684 cpassert(all(shape(rotation_generator) == [n, n]))
685 cpassert(SIZE(occupation) == n)
686 cpassert(SIZE(rotation_step) == nrotation)
687 cpassert(SIZE(weighted_occupation_step) == n)
688 cpassert(all(shape(density_tangent) == [n, n]))
689 cpassert(kpoint_weight > 0.0_dp)
690
691 ALLOCATE (direction(n, n), generator_minus(n, n), generator_plus(n, n), &
692 mode_ref(n, n), rotation(n, n), rotation_minus(n, n), &
693 rotation_plus(n, n), weighted_rotation(n, n))
694 direction(:, :) = cmplx(0.0_dp, 0.0_dp, kind=dp)
695 r = 0
696 DO i = 1, n - 1
697 DO j = i + 1, n
698 r = r + 1
699 direction(i, j) = cmplx(rotation_step(r), 0.0_dp, kind=dp)
700 direction(j, i) = -direction(i, j)
701 r = r + 1
702 direction(i, j) = direction(i, j) + &
703 cmplx(0.0_dp, rotation_step(r), kind=dp)
704 direction(j, i) = direction(j, i) + &
705 cmplx(0.0_dp, rotation_step(r), kind=dp)
706 END DO
707 END DO
708 cpassert(r == nrotation)
709
710 step_scale = max(1.0_dp, maxval(abs(direction)))
711 step = 1.0e-5_dp/step_scale
712 IF (PRESENT(difference_step)) step = difference_step/step_scale
713 cpassert(step > epsilon(step))
714 generator_plus(:, :) = rotation_generator + step*direction
715 generator_minus(:, :) = rotation_generator - step*direction
716 CALL qs_ot_dense_rotation_state(rotation_generator, rotation)
717 CALL qs_ot_dense_rotation_state(generator_plus, rotation_plus)
718 CALL qs_ot_dense_rotation_state(generator_minus, rotation_minus)
719
720 weighted_rotation(:, :) = rotation_plus
721 DO j = 1, n
722 weighted_rotation(:, j) = kpoint_weight*occupation(j)*weighted_rotation(:, j)
723 END DO
724 mode_ref(:, :) = matmul(weighted_rotation, conjg(transpose(rotation_plus)))
725 weighted_rotation(:, :) = rotation_minus
726 DO j = 1, n
727 weighted_rotation(:, j) = kpoint_weight*occupation(j)*weighted_rotation(:, j)
728 END DO
729 mode_ref(:, :) = (mode_ref - &
730 matmul(weighted_rotation, conjg(transpose(rotation_minus))))/(2.0_dp*step)
731
732 weighted_rotation(:, :) = rotation
733 DO j = 1, n
734 weighted_rotation(:, j) = weighted_occupation_step(j)*weighted_rotation(:, j)
735 END DO
736 mode_ref(:, :) = mode_ref + matmul(weighted_rotation, conjg(transpose(rotation)))
737 density_tangent(:, :) = matmul(conjg(transpose(rotation)), matmul(mode_ref, rotation))
738 density_tangent(:, :) = 0.5_dp*(density_tangent + conjg(transpose(density_tangent)))
739
740 DEALLOCATE (direction, generator_minus, generator_plus, mode_ref, rotation, &
741 rotation_minus, rotation_plus, weighted_rotation)
742
743 END SUBROUTINE qs_ot_density_tangent
744
745! **************************************************************************************************
746!> \brief project a physical density/Hamiltonian secant between moving orbital subspaces
747!>
748!> For separately S-orthonormal endpoint orbitals C0 and C1, O=C0^H*S*C1 retains the
749!> component of the accepted density step that leaves the old subspace. Density modes are
750!> represented in the current C1 basis and already contain the irreducible K-point weight.
751!> \param overlap_start_current cross overlap O
752!> \param occupation_start occupations at the accepted start
753!> \param occupation_current occupations at the accepted endpoint
754!> \param hamiltonian_step_start C0^H*(H1-H0)*C0
755!> \param hamiltonian_step_current C1^H*(H1-H0)*C1
756!> \param density_modes current-orbital density tangents
757!> \param kpoint_weight irreducible K-point weight
758!> \param density_norm_sq contribution to <Delta P,Delta P>
759!> \param response_work contribution to <Delta P,Delta H>
760!> \param density_overlap contributions <Delta P,B_r>
761!> \param response_overlap contributions <Delta H,B_r>
762!> \param valid whether a finite nonzero secant was available
763! **************************************************************************************************
764 SUBROUTINE qs_ot_density_secant_orbital_overlaps( &
765 overlap_start_current, occupation_start, occupation_current, hamiltonian_step_start, &
766 hamiltonian_step_current, density_modes, kpoint_weight, density_norm_sq, response_work, &
767 density_overlap, response_overlap, valid)
768
769 COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(IN) :: overlap_start_current
770 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: occupation_start, occupation_current
771 COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(IN) :: hamiltonian_step_start, &
772 hamiltonian_step_current
773 COMPLEX(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: density_modes
774 REAL(kind=dp), INTENT(IN) :: kpoint_weight
775 REAL(kind=dp), INTENT(OUT) :: density_norm_sq, response_work
776 REAL(kind=dp), DIMENSION(:), INTENT(OUT) :: density_overlap, response_overlap
777 LOGICAL, INTENT(OUT) :: valid
778
779 COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: cross_mode
780 INTEGER :: i, j, mode, n, nmode
781 REAL(kind=dp) :: cross_density, density_scale
782
783 n = SIZE(occupation_start)
784 nmode = SIZE(density_modes, 3)
785 cpassert(n > 0)
786 cpassert(SIZE(occupation_current) == n)
787 cpassert(all(shape(overlap_start_current) == [n, n]))
788 cpassert(all(shape(hamiltonian_step_start) == [n, n]))
789 cpassert(all(shape(hamiltonian_step_current) == [n, n]))
790 cpassert(SIZE(density_modes, 1) == n .AND. SIZE(density_modes, 2) == n)
791 cpassert(SIZE(density_overlap) == nmode .AND. SIZE(response_overlap) == nmode)
792
793 density_norm_sq = 0.0_dp
794 response_work = 0.0_dp
795 density_overlap(:) = 0.0_dp
796 response_overlap(:) = 0.0_dp
797 valid = .false.
798 IF (nmode <= 0 .OR. kpoint_weight <= tiny(1.0_dp)) RETURN
799
800 cross_density = 0.0_dp
801 DO j = 1, n
802 DO i = 1, n
803 cross_density = cross_density + occupation_start(i)*occupation_current(j)* &
804 abs(overlap_start_current(i, j))**2
805 END DO
806 END DO
807 density_norm_sq = kpoint_weight* &
808 (sum(occupation_start**2) + sum(occupation_current**2) - &
809 2.0_dp*cross_density)
810 response_work = 0.0_dp
811 DO i = 1, n
812 response_work = response_work + kpoint_weight* &
813 (occupation_current(i)*real(hamiltonian_step_current(i, i), kind=dp) - &
814 occupation_start(i)*real(hamiltonian_step_start(i, i), kind=dp))
815 END DO
816
817 ALLOCATE (cross_mode(n, n))
818 DO mode = 1, nmode
819 ! GCC 15 miscompiles the equivalent complex array expressions with -march=znver3.
820 CALL gemm_square(overlap_start_current, "N", density_modes(:, :, mode), "N", &
821 overlap_start_current, "C", cross_mode)
822 DO i = 1, n
823 density_overlap(mode) = density_overlap(mode) + &
824 occupation_current(i)* &
825 REAL(density_modes(i, i, mode), kind=dp) - &
826 occupation_start(i)*real(cross_mode(i, i), kind=dp)
827 END DO
828 response_overlap(mode) = &
829 sum(real(hamiltonian_step_current, kind=dp)* &
830 REAL(density_modes(:, :, mode), kind=dp) + &
831 aimag(hamiltonian_step_current)*aimag(density_modes(:, :, mode)))
832 END DO
833 DEALLOCATE (cross_mode)
834
835 density_scale = kpoint_weight* &
836 max(maxval(abs(occupation_start)), maxval(abs(occupation_current)))
837 IF (density_scale <= tiny(1.0_dp)) RETURN
838 IF (density_norm_sq <= 64.0_dp*epsilon(1.0_dp)*density_scale**2) RETURN
839 valid = density_norm_sq == density_norm_sq .AND. response_work == response_work .AND. &
840 abs(density_norm_sq) <= huge(1.0_dp) .AND. abs(response_work) <= huge(1.0_dp) .AND. &
841 all(density_overlap == density_overlap) .AND. &
842 all(response_overlap == response_overlap) .AND. &
843 all(abs(density_overlap) <= huge(1.0_dp)) .AND. &
844 all(abs(response_overlap) <= huge(1.0_dp))
845 IF (.NOT. valid) THEN
846 density_norm_sq = 0.0_dp
847 response_work = 0.0_dp
848 density_overlap(:) = 0.0_dp
849 response_overlap(:) = 0.0_dp
850 END IF
851
852 END SUBROUTINE qs_ot_density_secant_orbital_overlaps
853
854! **************************************************************************************************
855!> \brief fixed-N Frechet derivative of a smooth occupation projector
856!>
857!> The spectral divided-difference kernel is invariant under rotations inside a degenerate
858!> eigenspace. Its diagonal includes the chemical-potential response of the complete fixed-N
859!> group, while off-diagonal terms describe the physical change of the spectral projector.
860!> \param chc projected Hermitian Hamiltonian
861!> \param dchc Hermitian Hamiltonian perturbation
862!> \param occupation canonical occupations associated with the eigenvalues of chc
863!> \param kpoint_weight irreducible-k-point weight
864!> \param response_weight signed weighted occupation responses for this channel
865!> \param fixed_n_weight_sum susceptibility summed over the complete fixed-N group
866!> \param projector_derivative derivative of the weighted occupation projector
867!> \param density_factor optional representation-dependent density prefactor
868! **************************************************************************************************
869 SUBROUTINE qs_ot_fixed_n_projector_frechet( &
870 chc, dchc, occupation, kpoint_weight, response_weight, fixed_n_weight_sum, &
871 projector_derivative, density_factor)
872 COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(IN) :: chc, dchc
873 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: occupation
874 REAL(kind=dp), INTENT(IN) :: kpoint_weight
875 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: response_weight
876 REAL(kind=dp), INTENT(IN) :: fixed_n_weight_sum
877 COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: projector_derivative
878 REAL(kind=dp), INTENT(IN), OPTIONAL :: density_factor
879
880 COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: eigenvectors, kernel, work
881 INTEGER :: i, j, n
882 REAL(kind=dp) :: coefficient, denominator, factor, &
883 gap_tolerance, local_weight_sum, &
884 mu_numerator, mu_shift, scale
885 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: eigenvalues, weighted_occupation
886
887 n = SIZE(chc, 1)
888 cpassert(n > 0)
889 cpassert(all(shape(chc) == [n, n]))
890 cpassert(all(shape(dchc) == [n, n]))
891 cpassert(SIZE(occupation) == n)
892 cpassert(SIZE(response_weight) == n)
893 cpassert(all(shape(projector_derivative) == [n, n]))
894
895 factor = 1.0_dp
896 IF (PRESENT(density_factor)) factor = density_factor
897 projector_derivative(:, :) = cmplx(0.0_dp, 0.0_dp, kind=dp)
898 IF (factor == 0.0_dp .OR. kpoint_weight <= 0.0_dp) RETURN
899
900 ALLOCATE (eigenvectors(n, n), kernel(n, n), work(n, n), eigenvalues(n), &
901 weighted_occupation(n))
902 CALL diag_complex(chc, eigenvectors, eigenvalues)
903 work(:, :) = matmul(conjg(transpose(eigenvectors)), matmul(dchc, eigenvectors))
904 weighted_occupation(:) = factor*kpoint_weight*occupation(:)
905
906 local_weight_sum = sum(response_weight(:))
907 mu_numerator = 0.0_dp
908 DO i = 1, n
909 mu_numerator = mu_numerator + response_weight(i)* &
910 REAL(work(i, i), kind=dp)
911 END DO
912 mu_shift = qs_ot_fixed_n_response_mu_shift(mu_numerator, local_weight_sum, &
913 fixed_n_weight_sum)
914
915 scale = max(1.0_dp, maxval(abs(eigenvalues(:))))
916 gap_tolerance = sqrt(epsilon(1.0_dp))*scale
917 kernel(:, :) = cmplx(0.0_dp, 0.0_dp, kind=dp)
918 DO j = 1, n
919 DO i = 1, n
920 IF (i == j) THEN
921 coefficient = -factor*response_weight(i)
922 kernel(i, i) = cmplx(coefficient*(real(work(i, i), kind=dp) - mu_shift), &
923 0.0_dp, kind=dp)
924 ELSE
925 denominator = eigenvalues(i) - eigenvalues(j)
926 IF (abs(denominator) > gap_tolerance) THEN
927 coefficient = (weighted_occupation(i) - weighted_occupation(j))/denominator
928 ELSE
929 coefficient = -0.5_dp*factor*(response_weight(i) + response_weight(j))
930 END IF
931 kernel(i, j) = coefficient*work(i, j)
932 END IF
933 END DO
934 END DO
935
936 projector_derivative(:, :) = matmul(eigenvectors, &
937 matmul(kernel, conjg(transpose(eigenvectors))))
938 projector_derivative(:, :) = 0.5_dp*(projector_derivative + &
939 conjg(transpose(projector_derivative)))
940
941 DEALLOCATE (eigenvectors, kernel, work, eigenvalues, weighted_occupation)
942
943 END SUBROUTINE qs_ot_fixed_n_projector_frechet
944
945! **************************************************************************************************
946!> \brief finite complex REF rotation Hessian and Rayleigh-energy response
947!>
948!> The current projected Hamiltonian is pulled back through the finite rotation and then
949!> differentiated in the independent real-antisymmetric and imaginary-symmetric pair
950!> coordinates. This keeps the response consistent with the exponential chart used by REF
951!> OT instead of replacing it by an infinitesimal commutator away from the chart origin.
952!> \param chc current projected Hermitian Hamiltonian U^H H_ref U
953!> \param rotation_generator current anti-Hermitian REF generator
954!> \param occupation fixed occupations attached to the rotated columns
955!> \param kpoint_weight irreducible-k-point weight
956!> \param rotation_gradient gradient in interleaved real/imaginary pair coordinates
957!> \param rotation_hessian derivative of rotation_gradient in the same coordinates
958!> \param rayleigh_response derivative of diag(U^H H_ref U) with respect to the pair coordinates
959!> \param difference_step optional central finite-difference step for the Hessian action
960! **************************************************************************************************
961 SUBROUTINE qs_ot_finite_rotation_response( &
962 chc, rotation_generator, occupation, kpoint_weight, rotation_gradient, &
963 rotation_hessian, rayleigh_response, difference_step)
964 COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(IN) :: chc, rotation_generator
965 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: occupation
966 REAL(kind=dp), INTENT(IN) :: kpoint_weight
967 REAL(kind=dp), DIMENSION(:), INTENT(OUT) :: rotation_gradient
968 REAL(kind=dp), DIMENSION(:, :), INTENT(OUT) :: rotation_hessian, rayleigh_response
969 REAL(kind=dp), INTENT(IN), OPTIONAL :: difference_step
970
971 COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: base_hamiltonian, generator_minus, &
972 generator_plus, rotation
973 INTEGER :: i, j, n, nrotation, r, s
974 REAL(kind=dp) :: step
975 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: occupation_scale, rayleigh, &
976 rayleigh_minus, rayleigh_plus
977 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: gradient_imag, gradient_minus_imag, &
978 gradient_minus_real, gradient_plus_imag, gradient_plus_real, gradient_real
979
980 n = SIZE(chc, 1)
981 nrotation = n*(n - 1)
982 cpassert(n > 0)
983 cpassert(all(shape(chc) == [n, n]))
984 cpassert(all(shape(rotation_generator) == [n, n]))
985 cpassert(SIZE(occupation) == n)
986 cpassert(SIZE(rotation_gradient) == nrotation)
987 cpassert(all(shape(rotation_hessian) == [nrotation, nrotation]))
988 cpassert(all(shape(rayleigh_response) == [n, nrotation]))
989 cpassert(kpoint_weight > 0.0_dp)
990
991 step = 1.0e-4_dp
992 IF (PRESENT(difference_step)) step = difference_step
993 cpassert(step > sqrt(epsilon(step)))
994
995 ALLOCATE (base_hamiltonian(n, n), generator_minus(n, n), generator_plus(n, n), &
996 rotation(n, n), occupation_scale(n), rayleigh(n), rayleigh_minus(n), &
997 rayleigh_plus(n), gradient_imag(n, n), gradient_minus_imag(n, n), &
998 gradient_minus_real(n, n), gradient_plus_imag(n, n), &
999 gradient_plus_real(n, n), gradient_real(n, n))
1000
1001 CALL qs_ot_dense_rotation_state(rotation_generator, rotation)
1002 base_hamiltonian(:, :) = matmul(rotation, matmul(chc, conjg(transpose(rotation))))
1003 occupation_scale(:) = 2.0_dp*kpoint_weight*occupation(:)
1004 CALL qs_ot_dense_rotation_gradient(rotation_generator, base_hamiltonian, occupation_scale, &
1005 gradient_real, gradient_imag, rayleigh)
1006
1007 r = 0
1008 DO i = 1, n - 1
1009 DO j = i + 1, n
1010 r = r + 1
1011 rotation_gradient(r) = gradient_real(i, j)
1012 r = r + 1
1013 rotation_gradient(r) = gradient_imag(i, j)
1014 END DO
1015 END DO
1016 cpassert(r == nrotation)
1017
1018 s = 0
1019 DO i = 1, n - 1
1020 DO j = i + 1, n
1021 generator_plus(:, :) = rotation_generator(:, :)
1022 generator_minus(:, :) = rotation_generator(:, :)
1023 generator_plus(i, j) = generator_plus(i, j) + step
1024 generator_plus(j, i) = generator_plus(j, i) - step
1025 generator_minus(i, j) = generator_minus(i, j) - step
1026 generator_minus(j, i) = generator_minus(j, i) + step
1027 s = s + 1
1028 CALL qs_ot_dense_rotation_gradient(generator_plus, base_hamiltonian, occupation_scale, &
1029 gradient_plus_real, gradient_plus_imag, rayleigh_plus)
1030 CALL qs_ot_dense_rotation_gradient(generator_minus, base_hamiltonian, occupation_scale, &
1031 gradient_minus_real, gradient_minus_imag, rayleigh_minus)
1032 CALL qs_ot_pack_rotation_response( &
1033 gradient_plus_real, gradient_plus_imag, gradient_minus_real, gradient_minus_imag, &
1034 rayleigh_plus, rayleigh_minus, step, rotation_hessian(:, s), &
1035 rayleigh_response(:, s))
1036
1037 generator_plus(:, :) = rotation_generator(:, :)
1038 generator_minus(:, :) = rotation_generator(:, :)
1039 generator_plus(i, j) = generator_plus(i, j) + cmplx(0.0_dp, step, kind=dp)
1040 generator_plus(j, i) = generator_plus(j, i) + cmplx(0.0_dp, step, kind=dp)
1041 generator_minus(i, j) = generator_minus(i, j) - cmplx(0.0_dp, step, kind=dp)
1042 generator_minus(j, i) = generator_minus(j, i) - cmplx(0.0_dp, step, kind=dp)
1043 s = s + 1
1044 CALL qs_ot_dense_rotation_gradient(generator_plus, base_hamiltonian, occupation_scale, &
1045 gradient_plus_real, gradient_plus_imag, rayleigh_plus)
1046 CALL qs_ot_dense_rotation_gradient(generator_minus, base_hamiltonian, occupation_scale, &
1047 gradient_minus_real, gradient_minus_imag, rayleigh_minus)
1048 CALL qs_ot_pack_rotation_response( &
1049 gradient_plus_real, gradient_plus_imag, gradient_minus_real, gradient_minus_imag, &
1050 rayleigh_plus, rayleigh_minus, step, rotation_hessian(:, s), &
1051 rayleigh_response(:, s))
1052 END DO
1053 END DO
1054 cpassert(s == nrotation)
1055 rotation_hessian(:, :) = 0.5_dp*(rotation_hessian + transpose(rotation_hessian))
1056
1057 DEALLOCATE (base_hamiltonian, generator_minus, generator_plus, rotation, occupation_scale, &
1058 rayleigh, rayleigh_minus, rayleigh_plus, gradient_imag, gradient_minus_imag, &
1059 gradient_minus_real, gradient_plus_imag, gradient_plus_real, gradient_real)
1060
1061 END SUBROUTINE qs_ot_finite_rotation_response
1062
1063! **************************************************************************************************
1064!> \brief pack one finite complex rotation response column
1065!> \param gradient_plus_real real gradient at the positive endpoint
1066!> \param gradient_plus_imag imaginary gradient at the positive endpoint
1067!> \param gradient_minus_real real gradient at the negative endpoint
1068!> \param gradient_minus_imag imaginary gradient at the negative endpoint
1069!> \param rayleigh_plus Rayleigh energies at the positive endpoint
1070!> \param rayleigh_minus Rayleigh energies at the negative endpoint
1071!> \param step central finite-difference step
1072!> \param hessian_column packed rotation-Hessian column
1073!> \param rayleigh_column packed Rayleigh-response column
1074! **************************************************************************************************
1075 SUBROUTINE qs_ot_pack_rotation_response( &
1076 gradient_plus_real, gradient_plus_imag, gradient_minus_real, gradient_minus_imag, &
1077 rayleigh_plus, rayleigh_minus, step, hessian_column, rayleigh_column)
1078 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: gradient_plus_real, gradient_plus_imag, &
1079 gradient_minus_real, &
1080 gradient_minus_imag
1081 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: rayleigh_plus, rayleigh_minus
1082 REAL(kind=dp), INTENT(IN) :: step
1083 REAL(kind=dp), DIMENSION(:), INTENT(OUT) :: hessian_column, rayleigh_column
1084
1085 INTEGER :: i, j, n, r
1086
1087 n = SIZE(gradient_plus_real, 1)
1088 cpassert(all(shape(gradient_plus_real) == [n, n]))
1089 cpassert(all(shape(gradient_plus_imag) == [n, n]))
1090 cpassert(all(shape(gradient_minus_real) == [n, n]))
1091 cpassert(all(shape(gradient_minus_imag) == [n, n]))
1092 cpassert(SIZE(hessian_column) == n*(n - 1))
1093 cpassert(SIZE(rayleigh_column) == n)
1094
1095 r = 0
1096 DO i = 1, n - 1
1097 DO j = i + 1, n
1098 r = r + 1
1099 hessian_column(r) = &
1100 (gradient_plus_real(i, j) - gradient_minus_real(i, j))/(2.0_dp*step)
1101 r = r + 1
1102 hessian_column(r) = &
1103 (gradient_plus_imag(i, j) - gradient_minus_imag(i, j))/(2.0_dp*step)
1104 END DO
1105 END DO
1106 rayleigh_column(:) = (rayleigh_plus(:) - rayleigh_minus(:))/(2.0_dp*step)
1107
1108 END SUBROUTINE qs_ot_pack_rotation_response
1109
1110! **************************************************************************************************
1111!> \brief dense complex rotation and projected-Hamiltonian diagonal
1112!> \param rotation_generator anti-Hermitian REF generator
1113!> \param base_hamiltonian fixed Hamiltonian in the unrotated REF basis
1114!> \param occupation_scale twice the weighted occupations
1115!> \param gradient_real real-antisymmetric gradient component
1116!> \param gradient_imag imaginary-symmetric gradient component
1117!> \param rayleigh diagonal of the rotated Hamiltonian
1118! **************************************************************************************************
1119 SUBROUTINE qs_ot_dense_rotation_gradient( &
1120 rotation_generator, base_hamiltonian, occupation_scale, gradient_real, gradient_imag, rayleigh)
1121 COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(IN) :: rotation_generator, base_hamiltonian
1122 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: occupation_scale
1123 REAL(kind=dp), DIMENSION(:, :), INTENT(OUT) :: gradient_real, gradient_imag
1124 REAL(kind=dp), DIMENSION(:), INTENT(OUT) :: rayleigh
1125
1126 COMPLEX(KIND=dp) :: kernel
1127 COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: covector, eigenvectors, frechet, inner, &
1128 outer, rotation, work
1129 INTEGER :: i, j, n
1130 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: eigenvalues
1131
1132 n = SIZE(rotation_generator, 1)
1133 ALLOCATE (covector(n, n), eigenvectors(n, n), frechet(n, n), inner(n, n), &
1134 outer(n, n), rotation(n, n), work(n, n), eigenvalues(n))
1135 CALL qs_ot_dense_rotation_state(rotation_generator, rotation, &
1136 eigenvectors=eigenvectors, eigenvalues=eigenvalues)
1137
1138 work(:, :) = matmul(base_hamiltonian, rotation)
1139 covector(:, :) = work(:, :)
1140 DO j = 1, n
1141 covector(:, j) = occupation_scale(j)*covector(:, j)
1142 END DO
1143
1144 work(:, :) = matmul(conjg(transpose(rotation)), &
1145 matmul(base_hamiltonian, rotation))
1146 DO i = 1, n
1147 rayleigh(i) = real(work(i, i), kind=dp)
1148 END DO
1149
1150 inner(:, :) = matmul(conjg(transpose(eigenvectors)), &
1151 matmul(covector, eigenvectors))
1152 DO j = 1, n
1153 DO i = 1, n
1154 kernel = conjg(qs_ot_complex_exp_frechet_kernel(eigenvalues(i), eigenvalues(j)))
1155 outer(i, j) = inner(i, j)*kernel
1156 END DO
1157 END DO
1158 frechet(:, :) = matmul(eigenvectors, &
1159 matmul(outer, conjg(transpose(eigenvectors))))
1160 gradient_real(:, :) = real(frechet, kind=dp) - transpose(real(frechet, kind=dp))
1161 gradient_imag(:, :) = aimag(frechet) + transpose(aimag(frechet))
1162
1163 DEALLOCATE (covector, eigenvectors, frechet, inner, outer, rotation, work, eigenvalues)
1164
1165 END SUBROUTINE qs_ot_dense_rotation_gradient
1166
1167! **************************************************************************************************
1168!> \brief dense exponential of an anti-Hermitian REF generator
1169!> \param rotation_generator anti-Hermitian generator
1170!> \param rotation exp(rotation_generator)
1171!> \param eigenvectors optional eigenvectors of i*rotation_generator
1172!> \param eigenvalues optional eigenvalues of i*rotation_generator
1173! **************************************************************************************************
1174 SUBROUTINE qs_ot_dense_rotation_state(rotation_generator, rotation, eigenvectors, eigenvalues)
1175 COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(IN) :: rotation_generator
1176 COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: rotation
1177 COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(OUT), &
1178 OPTIONAL :: eigenvectors
1179 REAL(kind=dp), DIMENSION(:), INTENT(OUT), OPTIONAL :: eigenvalues
1180
1181 COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: vectors, weighted_vectors
1182 INTEGER :: j, n
1183 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: values
1184
1185 n = SIZE(rotation_generator, 1)
1186 ALLOCATE (vectors(n, n), weighted_vectors(n, n), values(n))
1187 CALL diag_complex(cmplx(0.0_dp, 1.0_dp, kind=dp)*rotation_generator, &
1188 vectors, values)
1189 weighted_vectors(:, :) = vectors(:, :)
1190 DO j = 1, n
1191 weighted_vectors(:, j) = exp(cmplx(0.0_dp, -values(j), kind=dp))* &
1192 weighted_vectors(:, j)
1193 END DO
1194 rotation(:, :) = matmul(weighted_vectors, conjg(transpose(vectors)))
1195 IF (PRESENT(eigenvectors)) eigenvectors(:, :) = vectors(:, :)
1196 IF (PRESENT(eigenvalues)) eigenvalues(:) = values(:)
1197 DEALLOCATE (vectors, weighted_vectors, values)
1198
1199 END SUBROUTINE qs_ot_dense_rotation_state
1200
1201! **************************************************************************************************
1202!> \brief Frechet divided-difference kernel for exp(-i*evals)
1203!> \param e1 ...
1204!> \param e2 ...
1205!> \return ...
1206! **************************************************************************************************
1207 PURE FUNCTION qs_ot_complex_exp_frechet_kernel(e1, e2) RESULT(kernel)
1208 REAL(kind=dp), INTENT(IN) :: e1, e2
1209 COMPLEX(KIND=dp) :: kernel
1210
1211 COMPLEX(KIND=dp) :: l1, l2, x
1212 INTEGER :: i
1213
1214 l1 = (0.0_dp, -1.0_dp)*e1
1215 l2 = (0.0_dp, -1.0_dp)*e2
1216 IF (abs(l1 - l2) > 0.5_dp) THEN
1217 kernel = (exp(l1) - exp(l2))/(l1 - l2)
1218 ELSE
1219 x = 1.0_dp
1220 kernel = 0.0_dp
1221 DO i = 1, 16
1222 kernel = kernel + x
1223 x = x*(l1 - l2)/real(i + 1, kind=dp)
1224 END DO
1225 kernel = kernel*exp(l2)
1226 END IF
1227
1228 END FUNCTION qs_ot_complex_exp_frechet_kernel
1229
1230! **************************************************************************************************
1231!> \brief apply the complex exponential Frechet kernel to sparse DBCSR Re/Im matrices
1232!> \param evals generator eigenvalues
1233!> \param inner_deriv_re real part of the matrix in the generator eigenbasis
1234!> \param inner_deriv_im imaginary part of the matrix in the generator eigenbasis
1235!> \param outer_deriv_re real part of the mapped matrix
1236!> \param outer_deriv_im imaginary part of the mapped matrix
1237!> \param adjoint use the adjoint Frechet kernel for gradients
1238! **************************************************************************************************
1239 SUBROUTINE qs_ot_apply_complex_frechet_dbcsr(evals, inner_deriv_re, inner_deriv_im, &
1240 outer_deriv_re, outer_deriv_im, adjoint)
1241 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: evals
1242 TYPE(dbcsr_type) :: inner_deriv_re, inner_deriv_im, &
1243 outer_deriv_re, outer_deriv_im
1244 LOGICAL, INTENT(IN), OPTIONAL :: adjoint
1245
1246 COMPLEX(KIND=dp) :: cval_in, kernel
1247 INTEGER :: col, i, j, max_blocks, nblocks, row
1248 INTEGER, ALLOCATABLE, DIMENSION(:) :: cols, rows
1249 INTEGER, DIMENSION(:), POINTER :: col_blk_offset, col_blk_size, &
1250 row_blk_offset, row_blk_size
1251 LOGICAL :: found_im, found_out_im, found_out_re, &
1252 found_re, use_adjoint
1253 REAL(dp), DIMENSION(:, :), POINTER :: block_in_im, block_in_re, block_out_im, &
1254 block_out_re
1255 REAL(kind=dp) :: e1, e2, im_part, re_part
1256 TYPE(dbcsr_distribution_type) :: dist
1257 TYPE(dbcsr_iterator_type) :: iter
1258
1259 use_adjoint = .false.
1260 IF (PRESENT(adjoint)) use_adjoint = adjoint
1261
1262 ! Re/Im parts can have different sparse block patterns. Build their union
1263 ! explicitly so a missing partner block is interpreted as zero.
1264 max_blocks = 0
1265 CALL dbcsr_iterator_start(iter, inner_deriv_re)
1266 DO WHILE (dbcsr_iterator_blocks_left(iter))
1267 CALL dbcsr_iterator_next_block(iter, row, col)
1268 max_blocks = max_blocks + 1
1269 END DO
1270 CALL dbcsr_iterator_stop(iter)
1271 CALL dbcsr_iterator_start(iter, inner_deriv_im)
1272 DO WHILE (dbcsr_iterator_blocks_left(iter))
1273 CALL dbcsr_iterator_next_block(iter, row, col)
1274 max_blocks = max_blocks + 1
1275 END DO
1276 CALL dbcsr_iterator_stop(iter)
1277 ALLOCATE (rows(max(max_blocks, 1)), cols(max(max_blocks, 1)))
1278 nblocks = 0
1279
1280 CALL dbcsr_iterator_start(iter, inner_deriv_re)
1281 DO WHILE (dbcsr_iterator_blocks_left(iter))
1282 CALL dbcsr_iterator_next_block(iter, row, col)
1283 CALL append_union_block(row, col)
1284 END DO
1285 CALL dbcsr_iterator_stop(iter)
1286 CALL dbcsr_iterator_start(iter, inner_deriv_im)
1287 DO WHILE (dbcsr_iterator_blocks_left(iter))
1288 CALL dbcsr_iterator_next_block(iter, row, col)
1289 CALL append_union_block(row, col)
1290 END DO
1291 CALL dbcsr_iterator_stop(iter)
1292
1293 CALL dbcsr_get_info(inner_deriv_re, distribution=dist, row_blk_size=row_blk_size, &
1294 col_blk_size=col_blk_size)
1295 CALL dbcsr_create(outer_deriv_re, "outer_deriv_re", dist, dbcsr_type_no_symmetry, &
1296 row_blk_size, col_blk_size)
1297 CALL dbcsr_create(outer_deriv_im, "outer_deriv_im", dist, dbcsr_type_no_symmetry, &
1298 row_blk_size, col_blk_size)
1299 IF (nblocks > 0) THEN
1300 CALL dbcsr_reserve_blocks(outer_deriv_re, rows=rows(1:nblocks), cols=cols(1:nblocks))
1301 CALL dbcsr_reserve_blocks(outer_deriv_im, rows=rows(1:nblocks), cols=cols(1:nblocks))
1302 END IF
1303 CALL dbcsr_finalize(outer_deriv_re)
1304 CALL dbcsr_finalize(outer_deriv_im)
1305 CALL dbcsr_set(outer_deriv_re, 0.0_dp)
1306 CALL dbcsr_set(outer_deriv_im, 0.0_dp)
1307
1308 CALL dbcsr_get_info(outer_deriv_re, row_blk_offset=row_blk_offset, &
1309 col_blk_offset=col_blk_offset)
1310 CALL dbcsr_iterator_start(iter, outer_deriv_re)
1311 DO WHILE (dbcsr_iterator_blocks_left(iter))
1312 CALL dbcsr_iterator_next_block(iter, row, col)
1313 CALL dbcsr_get_block_p(inner_deriv_re, row, col, block_in_re, found_re)
1314 CALL dbcsr_get_block_p(inner_deriv_im, row, col, block_in_im, found_im)
1315 CALL dbcsr_get_block_p(outer_deriv_re, row, col, block_out_re, found_out_re)
1316 CALL dbcsr_get_block_p(outer_deriv_im, row, col, block_out_im, found_out_im)
1317 cpassert(found_out_re .AND. found_out_im)
1318
1319 DO i = 1, SIZE(block_out_re, 1)
1320 DO j = 1, SIZE(block_out_re, 2)
1321 e1 = evals(row_blk_offset(row) + i - 1)
1322 e2 = evals(col_blk_offset(col) + j - 1)
1323 re_part = 0.0_dp
1324 im_part = 0.0_dp
1325 IF (found_re) re_part = block_in_re(i, j)
1326 IF (found_im) im_part = block_in_im(i, j)
1327 cval_in = cmplx(re_part, im_part, dp)
1328 kernel = qs_ot_complex_exp_frechet_kernel(e1, e2)
1329 IF (use_adjoint) kernel = conjg(kernel)
1330 cval_in = cval_in*kernel
1331 block_out_re(i, j) = real(cval_in, kind=dp)
1332 block_out_im(i, j) = aimag(cval_in)
1333 END DO
1334 END DO
1335 END DO
1336 CALL dbcsr_iterator_stop(iter)
1337 DEALLOCATE (rows, cols)
1338
1339 CONTAINS
1340
1341! **************************************************************************************************
1342!> \brief append a block coordinate unless it is already present
1343!> \param row_new block-row index
1344!> \param col_new block-column index
1345! **************************************************************************************************
1346 SUBROUTINE append_union_block(row_new, col_new)
1347 INTEGER, INTENT(IN) :: row_new, col_new
1348
1349 INTEGER :: iblock
1350
1351 DO iblock = 1, nblocks
1352 IF (rows(iblock) == row_new .AND. cols(iblock) == col_new) RETURN
1353 END DO
1354 nblocks = nblocks + 1
1355 rows(nblocks) = row_new
1356 cols(nblocks) = col_new
1357 END SUBROUTINE append_union_block
1358
1359 END SUBROUTINE qs_ot_apply_complex_frechet_dbcsr
1360
1361! **************************************************************************************************
1362!> \brief gets ready to use the preconditioner/ or renew the preconditioner
1363!> only keeps a pointer to the preconditioner.
1364!> If you change the preconditioner, you have to call this routine
1365!> you remain responsible of proper deallocate of your preconditioner
1366!> (or you can reuse it on the next step of the computation)
1367!> \param qs_ot_env ...
1368!> \param preconditioner ...
1369! **************************************************************************************************
1370 SUBROUTINE qs_ot_new_preconditioner(qs_ot_env, preconditioner)
1371 TYPE(qs_ot_type) :: qs_ot_env
1372 TYPE(preconditioner_type), POINTER :: preconditioner
1373
1374 INTEGER :: ncoef
1375
1376 qs_ot_env%preconditioner => preconditioner
1377 qs_ot_env%os_valid = .false.
1378 IF (.NOT. ASSOCIATED(qs_ot_env%matrix_psc0)) THEN
1379 CALL dbcsr_init_p(qs_ot_env%matrix_psc0)
1380 CALL dbcsr_copy(qs_ot_env%matrix_psc0, qs_ot_env%matrix_sc0, 'matrix_psc0')
1381 END IF
1382 IF (qs_ot_env%has_complex_kpoint_state .AND. &
1383 .NOT. ASSOCIATED(qs_ot_env%matrix_psc0_im)) THEN
1384 CALL dbcsr_init_p(qs_ot_env%matrix_psc0_im)
1385 CALL dbcsr_copy(qs_ot_env%matrix_psc0_im, qs_ot_env%matrix_sc0_im, 'matrix_psc0_im')
1386 END IF
1387
1388 IF (.NOT. qs_ot_env%use_dx) THEN
1389 qs_ot_env%use_dx = .true.
1390 CALL dbcsr_init_p(qs_ot_env%matrix_dx)
1391 CALL dbcsr_copy(qs_ot_env%matrix_dx, qs_ot_env%matrix_gx, 'matrix_dx')
1392 IF (qs_ot_env%has_complex_kpoint_state) THEN
1393 CALL dbcsr_init_p(qs_ot_env%matrix_dx_im)
1394 CALL dbcsr_copy(qs_ot_env%matrix_dx_im, qs_ot_env%matrix_gx_im, 'matrix_dx_im')
1395 END IF
1396 IF (qs_ot_env%settings%do_rotation) THEN
1397 CALL dbcsr_init_p(qs_ot_env%rot_mat_dx)
1398 CALL dbcsr_copy(qs_ot_env%rot_mat_dx, qs_ot_env%rot_mat_gx, 'rot_mat_dx')
1399 IF (qs_ot_env%has_complex_kpoint_state) THEN
1400 CALL dbcsr_init_p(qs_ot_env%rot_mat_dx_im)
1401 CALL dbcsr_copy(qs_ot_env%rot_mat_dx_im, qs_ot_env%rot_mat_gx_im, 'rot_mat_dx_im')
1402 END IF
1403 END IF
1404 IF (qs_ot_env%settings%do_ener) THEN
1405 ncoef = SIZE(qs_ot_env%ener_gx)
1406 ALLOCATE (qs_ot_env%ener_dx(ncoef))
1407 qs_ot_env%ener_dx = 0.0_dp
1408 END IF
1409 END IF
1410
1411 END SUBROUTINE qs_ot_new_preconditioner
1412
1413! **************************************************************************************************
1414!> \brief multiply paired real/imaginary DBCSR matrices, with optional conjugate transposes
1415!> \param op_a N or C
1416!> \param op_b N or C
1417!> \param a_re real part of A
1418!> \param a_im imaginary part of A
1419!> \param b_re real part of B
1420!> \param b_im imaginary part of B
1421!> \param c_re real part of A*B
1422!> \param c_im imaginary part of A*B
1423!> \param tmp real workspace shaped like C
1424! **************************************************************************************************
1425 SUBROUTINE qs_ot_complex_multiply(op_a, op_b, a_re, a_im, b_re, b_im, c_re, c_im, tmp)
1426 CHARACTER(LEN=1), INTENT(IN) :: op_a, op_b
1427 TYPE(dbcsr_type) :: a_re, a_im, b_re, b_im, c_re, c_im, tmp
1428
1429 CHARACTER(LEN=1) :: db_op_a, db_op_b
1430 REAL(kind=dp) :: sign_a, sign_b
1431
1432 SELECT CASE (op_a)
1433 CASE ('N')
1434 db_op_a = 'N'
1435 sign_a = 1.0_dp
1436 CASE ('C')
1437 db_op_a = 'T'
1438 sign_a = -1.0_dp
1439 CASE DEFAULT
1440 cpabort("Complex matrix product expects N or C for op_a")
1441 END SELECT
1442 SELECT CASE (op_b)
1443 CASE ('N')
1444 db_op_b = 'N'
1445 sign_b = 1.0_dp
1446 CASE ('C')
1447 db_op_b = 'T'
1448 sign_b = -1.0_dp
1449 CASE DEFAULT
1450 cpabort("Complex matrix product expects N or C for op_b")
1451 END SELECT
1452
1453 CALL dbcsr_multiply(db_op_a, db_op_b, 1.0_dp, a_re, b_re, 0.0_dp, c_re)
1454 CALL dbcsr_multiply(db_op_a, db_op_b, 1.0_dp, a_im, b_im, 0.0_dp, tmp)
1455 CALL dbcsr_add(c_re, tmp, alpha_scalar=1.0_dp, beta_scalar=-sign_a*sign_b)
1456
1457 CALL dbcsr_multiply(db_op_a, db_op_b, sign_b, a_re, b_im, 0.0_dp, c_im)
1458 CALL dbcsr_multiply(db_op_a, db_op_b, sign_a, a_im, b_re, 0.0_dp, tmp)
1459 CALL dbcsr_add(c_im, tmp, alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
1460
1461 END SUBROUTINE qs_ot_complex_multiply
1462
1463! **************************************************************************************************
1464!> \brief ...
1465!> \param qs_ot_env ...
1466!> \param C_NEW ...
1467!> \param SC ...
1468!> \param G_OLD ...
1469!> \param D ...
1470! **************************************************************************************************
1471 SUBROUTINE qs_ot_on_the_fly_localize(qs_ot_env, C_NEW, SC, G_OLD, D)
1472 !
1473 TYPE(qs_ot_type) :: qs_ot_env
1474 TYPE(dbcsr_type), POINTER :: c_new, sc, g_old, d
1475
1476 CHARACTER(len=*), PARAMETER :: routinen = 'qs_ot_on_the_fly_localize'
1477 INTEGER, PARAMETER :: taylor_order = 50
1478 REAL(kind=dp), PARAMETER :: alpha = 0.1_dp, f2_eps = 0.01_dp
1479
1480 INTEGER :: col, col_size, handle, i, k, n, p, row, &
1481 row_size
1482 REAL(dp), DIMENSION(:, :), POINTER :: block
1483 REAL(kind=dp) :: expfactor, f2, norm_fro, norm_gct, tmp
1484 TYPE(dbcsr_distribution_type) :: dist
1485 TYPE(dbcsr_iterator_type) :: iter
1486 TYPE(dbcsr_type), POINTER :: c, gp1, gp2, gu, u
1487 TYPE(mp_comm_type) :: group
1488
1489 CALL timeset(routinen, handle)
1490 !
1491 !
1492 CALL dbcsr_get_info(c_new, nfullrows_total=n, nfullcols_total=k)
1493 !
1494 ! C = C*expm(-G)
1495 gu => qs_ot_env%buf1_k_k_nosym ! a buffer
1496 u => qs_ot_env%buf2_k_k_nosym ! a buffer
1497 gp1 => qs_ot_env%buf3_k_k_nosym ! a buffer
1498 gp2 => qs_ot_env%buf4_k_k_nosym ! a buffer
1499 c => qs_ot_env%buf1_n_k ! a buffer
1500 !
1501 ! compute the derivative of the norm
1502 !-------------------------------------------------------------------
1503 ! (x^2+eps)^1/2
1504 f2 = 0.0_dp
1505 CALL dbcsr_copy(c, c_new)
1506 CALL dbcsr_iterator_start(iter, c)
1507 DO WHILE (dbcsr_iterator_blocks_left(iter))
1508 CALL dbcsr_iterator_next_block(iter, row, col, block, row_size=row_size, col_size=col_size)
1509 DO p = 1, col_size ! p
1510 DO i = 1, row_size ! i
1511 tmp = sqrt(block(i, p)**2 + f2_eps)
1512 f2 = f2 + tmp
1513 block(i, p) = block(i, p)/tmp
1514 END DO
1515 END DO
1516 END DO
1517 CALL dbcsr_iterator_stop(iter)
1518 CALL dbcsr_get_info(c, group=group)
1519 CALL group%sum(f2)
1520 !
1521 !
1522 CALL dbcsr_multiply('T', 'N', 1.0_dp, c, c_new, 0.0_dp, gu)
1523 !
1524 ! antisymetrize
1525 CALL dbcsr_get_info(gu, distribution=dist)
1526 CALL dbcsr_transposed(u, gu, shallow_data_copy=.false., &
1527 use_distribution=dist, &
1528 transpose_distribution=.false.)
1529 CALL dbcsr_add(gu, u, alpha_scalar=-0.5_dp, beta_scalar=0.5_dp)
1530 !-------------------------------------------------------------------
1531 !
1532 norm_fro = dbcsr_frobenius_norm(gu)
1533 norm_gct = dbcsr_gershgorin_norm(gu)
1534 !write(*,*) 'qs_ot_localize: ||P-I||_f=',norm_fro,' ||P-I||_GCT=',norm_gct
1535 !
1536 !kscale = CEILING(LOG(MIN(norm_fro,norm_gct))/LOG(2.0_dp))
1537 !scale = LOG(MIN(norm_fro,norm_gct))/LOG(2.0_dp)
1538 !write(*,*) 'qs_ot_localize: scale=',scale,' kscale=',kscale
1539 !
1540 ! rescale for steepest descent
1541 CALL dbcsr_scale(gu, -alpha)
1542 !
1543 ! compute unitary transform
1544 ! zeroth and first order
1545 expfactor = 1.0_dp
1546 CALL dbcsr_copy(u, gu)
1547 CALL dbcsr_scale(u, expfactor)
1548 CALL dbcsr_add_on_diag(u, 1.0_dp)
1549 ! other orders
1550 CALL dbcsr_copy(gp1, gu)
1551 DO i = 2, taylor_order
1552 ! new power of G
1553 CALL dbcsr_multiply('N', 'N', 1.0_dp, gu, gp1, 0.0_dp, gp2)
1554 CALL dbcsr_copy(gp1, gp2)
1555 ! add to the taylor expansion so far
1556 expfactor = expfactor/real(i, kind=dp)
1557 CALL dbcsr_add(u, gp1, alpha_scalar=1.0_dp, beta_scalar=expfactor)
1558 norm_fro = dbcsr_frobenius_norm(gp1)
1559 !write(*,*) 'Taylor expansion i=',i,' norm(X^i)/i!=',norm_fro*expfactor
1560 IF (norm_fro*expfactor < 1.0e-10_dp) EXIT
1561 END DO
1562 !
1563 ! rotate MOs
1564 CALL dbcsr_multiply('N', 'N', 1.0_dp, c_new, u, 0.0_dp, c)
1565 CALL dbcsr_copy(c_new, c)
1566 !
1567 ! rotate SC
1568 CALL dbcsr_multiply('N', 'N', 1.0_dp, sc, u, 0.0_dp, c)
1569 CALL dbcsr_copy(sc, c)
1570 !
1571 ! rotate D_i
1572 CALL dbcsr_multiply('N', 'N', 1.0_dp, d, u, 0.0_dp, c)
1573 CALL dbcsr_copy(d, c)
1574 !
1575 ! rotate G_i-1
1576 IF (ASSOCIATED(g_old)) THEN
1577 CALL dbcsr_multiply('N', 'N', 1.0_dp, g_old, u, 0.0_dp, c)
1578 CALL dbcsr_copy(g_old, c)
1579 END IF
1580 !
1581 CALL timestop(handle)
1582 END SUBROUTINE qs_ot_on_the_fly_localize
1583
1584! **************************************************************************************************
1585!> \brief ...
1586!> \param qs_ot_env ...
1587!> \param C_OLD ...
1588!> \param C_TMP ...
1589!> \param C_NEW ...
1590!> \param P ...
1591!> \param SC ...
1592!> \param update ...
1593! **************************************************************************************************
1594 SUBROUTINE qs_ot_ref_chol(qs_ot_env, C_OLD, C_TMP, C_NEW, P, SC, update)
1595 !
1596 TYPE(qs_ot_type) :: qs_ot_env
1597 TYPE(dbcsr_type) :: c_old, c_tmp, c_new, p, sc
1598 LOGICAL, INTENT(IN) :: update
1599
1600 CHARACTER(len=*), PARAMETER :: routinen = 'qs_ot_ref_chol'
1601
1602 INTEGER :: handle, k, n
1603
1604 CALL timeset(routinen, handle)
1605 !
1606 CALL dbcsr_get_info(c_new, nfullrows_total=n, nfullcols_total=k)
1607 !
1608 ! P = U'*U
1609 CALL cp_dbcsr_cholesky_decompose(p, k, qs_ot_env%para_env, qs_ot_env%blacs_env)
1610 !
1611 ! C_NEW = C_OLD*inv(U)
1612 CALL cp_dbcsr_cholesky_restore(c_old, k, p, c_new, op="SOLVE", pos="RIGHT", &
1613 transa="N", para_env=qs_ot_env%para_env, blacs_env=qs_ot_env%blacs_env)
1614 !
1615 ! Update SC if needed
1616 IF (update) THEN
1617 CALL cp_dbcsr_cholesky_restore(sc, k, p, c_tmp, op="SOLVE", pos="RIGHT", &
1618 transa="N", para_env=qs_ot_env%para_env, blacs_env=qs_ot_env%blacs_env)
1619 CALL dbcsr_copy(sc, c_tmp)
1620 END IF
1621 !
1622 CALL timestop(handle)
1623 END SUBROUTINE qs_ot_ref_chol
1624
1625! **************************************************************************************************
1626!> \brief ...
1627!> \param qs_ot_env ...
1628!> \param C_OLD ...
1629!> \param C_TMP ...
1630!> \param C_NEW ...
1631!> \param P ...
1632!> \param SC ...
1633!> \param update ...
1634! **************************************************************************************************
1635 SUBROUTINE qs_ot_ref_lwdn(qs_ot_env, C_OLD, C_TMP, C_NEW, P, SC, update)
1636 !
1637 TYPE(qs_ot_type) :: qs_ot_env
1638 TYPE(dbcsr_type) :: c_old, c_tmp, c_new, p, sc
1639 LOGICAL, INTENT(IN) :: update
1640
1641 CHARACTER(len=*), PARAMETER :: routinen = 'qs_ot_ref_lwdn'
1642
1643 INTEGER :: handle, i, k, n
1644 REAL(dp), ALLOCATABLE, DIMENSION(:) :: eig, fun
1645 TYPE(dbcsr_type), POINTER :: v, w
1646
1647 CALL timeset(routinen, handle)
1648 !
1649 CALL dbcsr_get_info(c_new, nfullrows_total=n, nfullcols_total=k)
1650 !
1651 v => qs_ot_env%buf1_k_k_nosym ! a buffer
1652 w => qs_ot_env%buf2_k_k_nosym ! a buffer
1653 ALLOCATE (eig(k), fun(k))
1654 !
1655 CALL cp_dbcsr_syevd(p, v, eig, qs_ot_env%para_env, qs_ot_env%blacs_env)
1656 !
1657 ! compute the P^(-1/2)
1658 DO i = 1, k
1659 IF (eig(i) <= 0.0_dp) THEN
1660 cpabort("P not positive definite")
1661 END IF
1662 IF (eig(i) < 1.0e-8_dp) THEN
1663 fun(i) = 0.0_dp
1664 ELSE
1665 fun(i) = 1.0_dp/sqrt(eig(i))
1666 END IF
1667 END DO
1668 CALL dbcsr_copy(w, v)
1669 CALL dbcsr_scale_by_vector(v, alpha=fun, side='right')
1670 CALL dbcsr_multiply('N', 'T', 1.0_dp, w, v, 0.0_dp, p)
1671 !
1672 ! Update C
1673 CALL dbcsr_multiply('N', 'N', 1.0_dp, c_old, p, 0.0_dp, c_new)
1674 !
1675 ! Update SC if needed
1676 IF (update) THEN
1677 CALL dbcsr_multiply('N', 'N', 1.0_dp, sc, p, 0.0_dp, c_tmp)
1678 CALL dbcsr_copy(sc, c_tmp)
1679 END IF
1680 !
1681 DEALLOCATE (eig, fun)
1682 !
1683 CALL timestop(handle)
1684 END SUBROUTINE qs_ot_ref_lwdn
1685
1686! **************************************************************************************************
1687!> \brief ...
1688!> \param qs_ot_env ...
1689!> \param C_OLD ...
1690!> \param C_TMP ...
1691!> \param C_NEW ...
1692!> \param P ...
1693!> \param SC ...
1694!> \param norm_in ...
1695!> \param update ...
1696! **************************************************************************************************
1697 SUBROUTINE qs_ot_ref_poly(qs_ot_env, C_OLD, C_TMP, C_NEW, P, SC, norm_in, update)
1698 !
1699 TYPE(qs_ot_type) :: qs_ot_env
1700 TYPE(dbcsr_type), POINTER :: c_old, c_tmp, c_new, p
1701 TYPE(dbcsr_type) :: sc
1702 REAL(dp), INTENT(IN) :: norm_in
1703 LOGICAL, INTENT(IN) :: update
1704
1705 CHARACTER(len=*), PARAMETER :: routinen = 'qs_ot_ref_poly'
1706
1707 INTEGER :: handle, irefine, k, n
1708 LOGICAL :: quick_exit
1709 REAL(dp) :: norm, norm_fro, norm_gct, occ_in, &
1710 occ_out, rescale
1711 TYPE(dbcsr_type), POINTER :: buf1, buf2, buf_nosym, ft, fy
1712
1713 CALL timeset(routinen, handle)
1714 !
1715 CALL dbcsr_get_info(c_new, nfullrows_total=n, nfullcols_total=k)
1716 !
1717 buf_nosym => qs_ot_env%buf1_k_k_nosym ! a buffer
1718 buf1 => qs_ot_env%buf1_k_k_sym ! a buffer
1719 buf2 => qs_ot_env%buf2_k_k_sym ! a buffer
1720 fy => qs_ot_env%buf3_k_k_sym ! a buffer
1721 ft => qs_ot_env%buf4_k_k_sym ! a buffer
1722 !
1723 ! initialize the norm (already computed in qs_ot_get_orbitals_ref)
1724 norm = norm_in
1725 !
1726 ! can we do a quick exit?
1727 quick_exit = .false.
1728 IF (norm < qs_ot_env%settings%eps_irac_quick_exit) quick_exit = .true.
1729 !
1730 ! lets refine
1731 rescale = 1.0_dp
1732 DO irefine = 1, qs_ot_env%settings%max_irac
1733 !
1734 ! rescaling
1735 IF (norm > 1.0_dp) THEN
1736 CALL dbcsr_scale(p, 1.0_dp/norm)
1737 rescale = rescale/sqrt(norm)
1738 END IF
1739 !
1740 ! get the refinement polynomial
1741 CALL qs_ot_refine(p, fy, buf1, buf2, qs_ot_env%settings%irac_degree, &
1742 qs_ot_env%settings%eps_irac_filter_matrix)
1743 !
1744 ! collect the transformation
1745 IF (irefine == 1) THEN
1746 CALL dbcsr_copy(ft, fy, name='FT')
1747 ELSE
1748 CALL dbcsr_multiply('N', 'N', 1.0_dp, ft, fy, 0.0_dp, buf1)
1749 IF (qs_ot_env%settings%eps_irac_filter_matrix > 0.0_dp) THEN
1750 occ_in = dbcsr_get_occupation(buf1)
1751 CALL dbcsr_filter(buf1, qs_ot_env%settings%eps_irac_filter_matrix)
1752 occ_out = dbcsr_get_occupation(buf1)
1753 END IF
1754 CALL dbcsr_copy(ft, buf1, name='FT')
1755 END IF
1756 !
1757 ! quick exit if possible
1758 IF (quick_exit) THEN
1759 EXIT
1760 END IF
1761 !
1762 ! P = FY^T * P * FY
1763 CALL dbcsr_multiply('N', 'N', 1.0_dp, p, fy, 0.0_dp, buf_nosym)
1764 IF (qs_ot_env%settings%eps_irac_filter_matrix > 0.0_dp) THEN
1765 occ_in = dbcsr_get_occupation(buf_nosym)
1766 CALL dbcsr_filter(buf_nosym, qs_ot_env%settings%eps_irac_filter_matrix)
1767 occ_out = dbcsr_get_occupation(buf_nosym)
1768 END IF
1769 CALL dbcsr_multiply('N', 'N', 1.0_dp, fy, buf_nosym, 0.0_dp, p)
1770 IF (qs_ot_env%settings%eps_irac_filter_matrix > 0.0_dp) THEN
1771 occ_in = dbcsr_get_occupation(p)
1772 CALL dbcsr_filter(p, qs_ot_env%settings%eps_irac_filter_matrix)
1773 occ_out = dbcsr_get_occupation(p)
1774 END IF
1775 !
1776 ! check ||P-1||_gct
1777 CALL dbcsr_add_on_diag(p, -1.0_dp)
1778 norm_fro = dbcsr_frobenius_norm(p)
1779 norm_gct = dbcsr_gershgorin_norm(p)
1780 CALL dbcsr_add_on_diag(p, 1.0_dp)
1781 norm = min(norm_gct, norm_fro)
1782 !
1783 ! printing
1784 !
1785 ! blows up
1786 IF (norm > 1.0e10_dp) THEN
1787 CALL cp_abort(__location__, &
1788 "Refinement blows up! "// &
1789 "We need you to improve the code, please post your input on "// &
1790 "the forum https://www.cp2k.org/")
1791 END IF
1792 !
1793 ! can we do a quick exit next step?
1794 IF (norm < qs_ot_env%settings%eps_irac_quick_exit) quick_exit = .true.
1795 !
1796 ! are we done?
1797 IF (norm < qs_ot_env%settings%eps_irac) EXIT
1798 !
1799 END DO
1800 !
1801 ! C_NEW = C_NEW * FT * rescale
1802 CALL dbcsr_multiply('N', 'N', rescale, c_old, ft, 0.0_dp, c_new)
1803 IF (qs_ot_env%settings%eps_irac_filter_matrix > 0.0_dp) THEN
1804 occ_in = dbcsr_get_occupation(c_new)
1805 CALL dbcsr_filter(c_new, qs_ot_env%settings%eps_irac_filter_matrix)
1806 occ_out = dbcsr_get_occupation(c_new)
1807 END IF
1808 !
1809 ! update SC = SC * FY * rescale
1810 IF (update) THEN
1811 CALL dbcsr_multiply('N', 'N', rescale, sc, ft, 0.0_dp, c_tmp)
1812 IF (qs_ot_env%settings%eps_irac_filter_matrix > 0.0_dp) THEN
1813 occ_in = dbcsr_get_occupation(c_tmp)
1814 CALL dbcsr_filter(c_tmp, qs_ot_env%settings%eps_irac_filter_matrix)
1815 occ_out = dbcsr_get_occupation(c_tmp)
1816 END IF
1817 CALL dbcsr_copy(sc, c_tmp)
1818 END IF
1819 !
1820 CALL timestop(handle)
1821 END SUBROUTINE qs_ot_ref_poly
1822
1823! **************************************************************************************************
1824!> \brief ...
1825!> \param qs_ot_env1 ...
1826!> \return ...
1827! **************************************************************************************************
1828 FUNCTION qs_ot_ref_update(qs_ot_env1) RESULT(update)
1829 !
1830 TYPE(qs_ot_type) :: qs_ot_env1
1831 LOGICAL :: update
1832
1833 update = .false.
1834 SELECT CASE (qs_ot_env1%settings%ot_method)
1835 CASE ("CG", "SD")
1836 SELECT CASE (qs_ot_env1%settings%line_search_method)
1837 CASE ("2PNT")
1838 IF (qs_ot_env1%line_search_count == 2) update = .true.
1839 CASE DEFAULT
1840 cpabort("NYI")
1841 END SELECT
1842 CASE ("DIIS")
1843 update = .true.
1844 CASE ("BROY", "LBFG")
1845 ! These minimizers retain positions or secants in one fixed REF chart.
1846 update = .false.
1847 CASE DEFAULT
1848 cpabort("NYI")
1849 END SELECT
1850 END FUNCTION qs_ot_ref_update
1851
1852! **************************************************************************************************
1853!> \brief ...
1854!> \param qs_ot_env1 ...
1855!> \param norm_in ...
1856!> \param ortho_irac ...
1857! **************************************************************************************************
1858 SUBROUTINE qs_ot_ref_decide(qs_ot_env1, norm_in, ortho_irac)
1859 !
1860 TYPE(qs_ot_type) :: qs_ot_env1
1861 REAL(dp), INTENT(IN) :: norm_in
1862 CHARACTER(LEN=*), INTENT(INOUT) :: ortho_irac
1863
1864 ortho_irac = qs_ot_env1%settings%ortho_irac
1865 IF (norm_in < qs_ot_env1%settings%eps_irac_switch) ortho_irac = "POLY"
1866 END SUBROUTINE qs_ot_ref_decide
1867
1868! **************************************************************************************************
1869!> \brief ...
1870!> \param matrix_c ...
1871!> \param matrix_s ...
1872!> \param matrix_x ...
1873!> \param matrix_sx ...
1874!> \param matrix_gx_old ...
1875!> \param matrix_dx ...
1876!> \param qs_ot_env ...
1877!> \param qs_ot_env1 ...
1878! **************************************************************************************************
1879 SUBROUTINE qs_ot_get_orbitals_ref(matrix_c, matrix_s, matrix_x, matrix_sx, &
1880 matrix_gx_old, matrix_dx, qs_ot_env, qs_ot_env1)
1881 !
1882 TYPE(dbcsr_type), POINTER :: matrix_c, matrix_s, matrix_x, matrix_sx, &
1883 matrix_gx_old, matrix_dx
1884 TYPE(qs_ot_type) :: qs_ot_env, qs_ot_env1
1885
1886 CHARACTER(len=*), PARAMETER :: routinen = 'qs_ot_get_orbitals_ref'
1887
1888 CHARACTER(LEN=4) :: ortho_irac
1889 INTEGER :: handle, k, n
1890 LOGICAL :: on_the_fly_loc, update
1891 REAL(dp) :: norm, norm_fro, norm_gct, occ_in, occ_out
1892 TYPE(dbcsr_type), POINTER :: c_new, c_old, c_tmp, d, g_old, p, s, sc
1893
1894 CALL timeset(routinen, handle)
1895
1896 CALL dbcsr_get_info(matrix_c, nfullrows_total=n, nfullcols_total=k)
1897 !
1898 c_new => matrix_c
1899 c_old => matrix_x ! need to be carefully updated for the gradient !
1900 sc => matrix_sx ! need to be carefully updated for the gradient !
1901 g_old => matrix_gx_old ! need to be carefully updated for localization !
1902 d => matrix_dx ! need to be carefully updated for localization !
1903 s => matrix_s
1904
1905 p => qs_ot_env%p_k_k_sym ! a buffer
1906 c_tmp => qs_ot_env%buf1_n_k ! a buffer
1907 !
1908 ! do we need to update C_OLD and SC?
1909 update = qs_ot_ref_update(qs_ot_env1)
1910 !
1911 ! do we want to on the fly localize?
1912 ! for the moment this is set from the input,
1913 ! later we might want to localize every n-step or
1914 ! when the sparsity increases...
1915 on_the_fly_loc = qs_ot_env1%settings%on_the_fly_loc
1916 !
1917 ! compute SC = S*C
1918 IF (ASSOCIATED(s)) THEN
1919 CALL dbcsr_multiply('N', 'N', 1.0_dp, s, c_old, 0.0_dp, sc)
1920 IF (qs_ot_env1%settings%eps_irac_filter_matrix > 0.0_dp) THEN
1921 occ_in = dbcsr_get_occupation(sc)
1922 CALL dbcsr_filter(sc, qs_ot_env1%settings%eps_irac_filter_matrix)
1923 occ_out = dbcsr_get_occupation(sc)
1924 END IF
1925 ELSE
1926 CALL dbcsr_copy(sc, c_old)
1927 END IF
1928 !
1929 ! compute P = C'*SC
1930 CALL dbcsr_multiply('T', 'N', 1.0_dp, c_old, sc, 0.0_dp, p)
1931 IF (qs_ot_env1%settings%eps_irac_filter_matrix > 0.0_dp) THEN
1932 occ_in = dbcsr_get_occupation(p)
1933 CALL dbcsr_filter(p, qs_ot_env1%settings%eps_irac_filter_matrix)
1934 occ_out = dbcsr_get_occupation(p)
1935 END IF
1936 !
1937 ! check ||P-1||_f and ||P-1||_gct
1938 CALL dbcsr_add_on_diag(p, -1.0_dp)
1939 norm_fro = dbcsr_frobenius_norm(p)
1940 norm_gct = dbcsr_gershgorin_norm(p)
1941 CALL dbcsr_add_on_diag(p, 1.0_dp)
1942 norm = min(norm_gct, norm_fro)
1943 CALL qs_ot_ref_decide(qs_ot_env1, norm, ortho_irac)
1944 !
1945 ! select the orthogonality method
1946 SELECT CASE (ortho_irac)
1947 CASE ("CHOL")
1948 CALL qs_ot_ref_chol(qs_ot_env, c_old, c_tmp, c_new, p, sc, update)
1949 CASE ("LWDN")
1950 CALL qs_ot_ref_lwdn(qs_ot_env, c_old, c_tmp, c_new, p, sc, update)
1951 CASE ("POLY")
1952 CALL qs_ot_ref_poly(qs_ot_env, c_old, c_tmp, c_new, p, sc, norm, update)
1953 CASE DEFAULT
1954 cpabort("Wrong argument")
1955 END SELECT
1956 !
1957 ! We update the C_i+1 and localization
1958 IF (update) THEN
1959 IF (on_the_fly_loc) THEN
1960 CALL qs_ot_on_the_fly_localize(qs_ot_env, c_new, sc, g_old, d)
1961 END IF
1962 CALL dbcsr_copy(c_old, c_new)
1963 END IF
1964
1965 IF (qs_ot_env%settings%do_rotation) THEN
1966 CALL qs_ot_generate_rotation(qs_ot_env)
1967 CALL dbcsr_multiply('N', 'N', 1.0_dp, c_new, qs_ot_env%rot_mat_u, &
1968 0.0_dp, c_tmp)
1969 CALL dbcsr_copy(c_new, c_tmp)
1970 END IF
1971 !
1972 CALL timestop(handle)
1973 END SUBROUTINE qs_ot_get_orbitals_ref
1974
1975! **************************************************************************************************
1976!> \brief update complex REF k-point orbitals and their S(k)C(k) images
1977!> \param matrix_c ...
1978!> \param matrix_c_im ...
1979!> \param matrix_s ...
1980!> \param matrix_s_im ...
1981!> \param qs_ot_env ...
1982!> \param qs_ot_env1 environment carrying the shared minimizer state
1983! **************************************************************************************************
1984 SUBROUTINE qs_ot_get_orbitals_ref_complex(matrix_c, matrix_c_im, matrix_s, matrix_s_im, &
1985 qs_ot_env, qs_ot_env1)
1986
1987 TYPE(dbcsr_type), POINTER :: matrix_c, matrix_c_im, matrix_s, &
1988 matrix_s_im
1989 TYPE(qs_ot_type) :: qs_ot_env
1990 TYPE(qs_ot_type), OPTIONAL :: qs_ot_env1
1991
1992 CHARACTER(len=*), PARAMETER :: routinen = 'qs_ot_get_orbitals_ref_complex'
1993
1994 INTEGER :: handle, i, k, n
1995 LOGICAL :: update
1996 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: eigenvalues, inverse_sqrt
1997 TYPE(dbcsr_type) :: rotated_im, rotated_re, rotation_tmp
1998 TYPE(dbcsr_type), POINTER :: c_im, c_re, f_im, f_re, p_im, p_re, &
1999 sc_im, sc_re, tmp_kk, tmp_nk, v_im, &
2000 v_re, w_im, w_re
2001
2002 CALL timeset(routinen, handle)
2003
2004 cpassert(qs_ot_env%has_complex_kpoint_state)
2005 cpassert(ASSOCIATED(matrix_s))
2006 cpassert(ASSOCIATED(matrix_s_im))
2007
2008 c_re => qs_ot_env%matrix_x
2009 c_im => qs_ot_env%matrix_x_im
2010 f_re => qs_ot_env%matrix_ref_inv_sqrt
2011 f_im => qs_ot_env%matrix_ref_inv_sqrt_im
2012 sc_re => qs_ot_env%matrix_sx
2013 sc_im => qs_ot_env%matrix_sx_im
2014 p_re => qs_ot_env%buf1_k_k_sym
2015 p_im => qs_ot_env%buf2_k_k_sym
2016 v_re => qs_ot_env%buf3_k_k_sym
2017 v_im => qs_ot_env%buf4_k_k_sym
2018 w_re => qs_ot_env%buf1_k_k_nosym
2019 w_im => qs_ot_env%buf2_k_k_nosym
2020 tmp_kk => qs_ot_env%buf3_k_k_nosym
2021 tmp_nk => qs_ot_env%buf1_n_k
2022
2023 CALL dbcsr_get_info(c_re, nfullrows_total=n, nfullcols_total=k)
2024 IF (PRESENT(qs_ot_env1)) THEN
2025 update = qs_ot_ref_update(qs_ot_env1)
2026 ELSE
2027 update = qs_ot_ref_update(qs_ot_env)
2028 END IF
2029
2030 ! SC = (S_re + i*S_im) * (C_re + i*C_im)
2031 CALL dbcsr_multiply('N', 'N', 1.0_dp, matrix_s, c_re, 0.0_dp, sc_re)
2032 CALL dbcsr_multiply('N', 'N', 1.0_dp, matrix_s_im, c_im, 0.0_dp, tmp_nk)
2033 CALL dbcsr_add(sc_re, tmp_nk, alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
2034
2035 CALL dbcsr_multiply('N', 'N', 1.0_dp, matrix_s, c_im, 0.0_dp, sc_im)
2036 CALL dbcsr_multiply('N', 'N', 1.0_dp, matrix_s_im, c_re, 0.0_dp, tmp_nk)
2037 CALL dbcsr_add(sc_im, tmp_nk, alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
2038
2039 ! P = C^H*S*C. Its imaginary component is real antisymmetric.
2040 CALL dbcsr_multiply('T', 'N', 1.0_dp, c_re, sc_re, 0.0_dp, p_re)
2041 CALL dbcsr_multiply('T', 'N', 1.0_dp, c_im, sc_im, 0.0_dp, tmp_kk)
2042 CALL dbcsr_add(p_re, tmp_kk, alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
2043
2044 CALL dbcsr_multiply('T', 'N', 1.0_dp, c_re, sc_im, 0.0_dp, p_im)
2045 CALL dbcsr_multiply('T', 'N', 1.0_dp, c_im, sc_re, 0.0_dp, tmp_kk)
2046 CALL dbcsr_add(p_im, tmp_kk, alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
2047
2048 ALLOCATE (eigenvalues(k), inverse_sqrt(k))
2049 CALL cp_dbcsr_heevd(matrix_re=p_re, matrix_im=p_im, &
2050 eigenvectors_re=v_re, eigenvectors_im=v_im, &
2051 eigenvalues=eigenvalues, para_env=qs_ot_env%para_env, &
2052 blacs_env=qs_ot_env%blacs_env)
2053 DO i = 1, k
2054 IF (eigenvalues(i) <= epsilon(1.0_dp)) THEN
2055 cpabort("Complex REF overlap is not positive definite")
2056 END IF
2057 inverse_sqrt(i) = 1.0_dp/sqrt(eigenvalues(i))
2058 END DO
2059
2060 ! P^(-1/2) = V*diag(lambda^(-1/2))*V^H.
2061 CALL dbcsr_copy(w_re, v_re)
2062 CALL dbcsr_copy(w_im, v_im)
2063 CALL dbcsr_scale_by_vector(w_re, alpha=inverse_sqrt, side='right')
2064 CALL dbcsr_scale_by_vector(w_im, alpha=inverse_sqrt, side='right')
2065
2066 CALL dbcsr_multiply('N', 'T', 1.0_dp, w_re, v_re, 0.0_dp, p_re)
2067 CALL dbcsr_multiply('N', 'T', 1.0_dp, w_im, v_im, 0.0_dp, tmp_kk)
2068 CALL dbcsr_add(p_re, tmp_kk, alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
2069
2070 CALL dbcsr_multiply('N', 'T', 1.0_dp, w_im, v_re, 0.0_dp, p_im)
2071 CALL dbcsr_multiply('N', 'T', 1.0_dp, w_re, v_im, 0.0_dp, tmp_kk)
2072 CALL dbcsr_add(p_im, tmp_kk, alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
2073
2074 CALL dbcsr_copy(f_re, p_re)
2075 CALL dbcsr_copy(f_im, p_im)
2076
2077 ! Return the physical, orthonormal orbitals without changing a rejected
2078 ! line-search coordinate.
2079 CALL dbcsr_multiply('N', 'N', 1.0_dp, c_re, p_re, 0.0_dp, matrix_c)
2080 CALL dbcsr_multiply('N', 'N', 1.0_dp, c_im, p_im, 0.0_dp, tmp_nk)
2081 CALL dbcsr_add(matrix_c, tmp_nk, alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
2082
2083 CALL dbcsr_multiply('N', 'N', 1.0_dp, c_re, p_im, 0.0_dp, matrix_c_im)
2084 CALL dbcsr_multiply('N', 'N', 1.0_dp, c_im, p_re, 0.0_dp, tmp_nk)
2085 CALL dbcsr_add(matrix_c_im, tmp_nk, alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
2086
2087 IF (update) THEN
2088 CALL dbcsr_copy(c_re, matrix_c)
2089 CALL dbcsr_copy(c_im, matrix_c_im)
2090
2091 CALL dbcsr_set(f_re, 0.0_dp)
2092 CALL dbcsr_add_on_diag(f_re, alpha=1.0_dp)
2093 CALL dbcsr_set(f_im, 0.0_dp)
2094
2095 CALL dbcsr_multiply('N', 'N', 1.0_dp, matrix_s, c_re, 0.0_dp, sc_re)
2096 CALL dbcsr_multiply('N', 'N', 1.0_dp, matrix_s_im, c_im, 0.0_dp, tmp_nk)
2097 CALL dbcsr_add(sc_re, tmp_nk, alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
2098
2099 CALL dbcsr_multiply('N', 'N', 1.0_dp, matrix_s, c_im, 0.0_dp, sc_im)
2100 CALL dbcsr_multiply('N', 'N', 1.0_dp, matrix_s_im, c_re, 0.0_dp, tmp_nk)
2101 CALL dbcsr_add(sc_im, tmp_nk, alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
2102 END IF
2103
2104 IF (qs_ot_env%settings%do_rotation) THEN
2105 CALL qs_ot_generate_rotation_complex(qs_ot_env)
2106
2107 CALL dbcsr_copy(rotated_re, matrix_c, name="rotated_re")
2108 CALL dbcsr_copy(rotated_im, matrix_c_im, name="rotated_im")
2109 CALL dbcsr_copy(rotation_tmp, matrix_c, name="rotation_tmp")
2110
2111 ! C_out = Q*U for complex Q and unitary U.
2112 CALL dbcsr_multiply('N', 'N', 1.0_dp, matrix_c, qs_ot_env%rot_mat_u, &
2113 0.0_dp, rotated_re)
2114 CALL dbcsr_multiply('N', 'N', 1.0_dp, matrix_c_im, qs_ot_env%rot_mat_u_im, &
2115 0.0_dp, rotation_tmp)
2116 CALL dbcsr_add(rotated_re, rotation_tmp, alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
2117
2118 CALL dbcsr_multiply('N', 'N', 1.0_dp, matrix_c, qs_ot_env%rot_mat_u_im, &
2119 0.0_dp, rotated_im)
2120 CALL dbcsr_multiply('N', 'N', 1.0_dp, matrix_c_im, qs_ot_env%rot_mat_u, &
2121 0.0_dp, rotation_tmp)
2122 CALL dbcsr_add(rotated_im, rotation_tmp, alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
2123
2124 CALL dbcsr_copy(matrix_c, rotated_re)
2125 CALL dbcsr_copy(matrix_c_im, rotated_im)
2126 CALL dbcsr_release(rotated_re)
2127 CALL dbcsr_release(rotated_im)
2128 CALL dbcsr_release(rotation_tmp)
2129 END IF
2130
2131 DEALLOCATE (eigenvalues, inverse_sqrt)
2132
2133 CALL timestop(handle)
2134 END SUBROUTINE qs_ot_get_orbitals_ref_complex
2135
2136! **************************************************************************************************
2137!> \brief refinement polynomial of degree 2,3 and 4 (PRB 70, 193102 (2004))
2138!> \param P ...
2139!> \param FY ...
2140!> \param P2 ...
2141!> \param T ...
2142!> \param irac_degree ...
2143!> \param eps_irac_filter_matrix ...
2144! **************************************************************************************************
2145 SUBROUTINE qs_ot_refine(P, FY, P2, T, irac_degree, eps_irac_filter_matrix)
2146 TYPE(dbcsr_type), INTENT(inout) :: p, fy, p2, t
2147 INTEGER, INTENT(in) :: irac_degree
2148 REAL(dp), INTENT(in) :: eps_irac_filter_matrix
2149
2150 CHARACTER(len=*), PARAMETER :: routinen = 'qs_ot_refine'
2151
2152 INTEGER :: handle, k
2153 REAL(dp) :: occ_in, occ_out, r
2154
2155 CALL timeset(routinen, handle)
2156
2157 CALL dbcsr_get_info(p, nfullcols_total=k)
2158 SELECT CASE (irac_degree)
2159 CASE (2)
2160 ! C_out = C_in * ( 15/8 * I - 10/8 * P + 3/8 * P^2)
2161 r = 3.0_dp/8.0_dp
2162 CALL dbcsr_multiply('N', 'N', r, p, p, 0.0_dp, fy)
2163 IF (eps_irac_filter_matrix > 0.0_dp) THEN
2164 occ_in = dbcsr_get_occupation(fy)
2165 CALL dbcsr_filter(fy, eps_irac_filter_matrix)
2166 occ_out = dbcsr_get_occupation(fy)
2167 END IF
2168 r = -10.0_dp/8.0_dp
2169 CALL dbcsr_add(fy, p, alpha_scalar=1.0_dp, beta_scalar=r)
2170 r = 15.0_dp/8.0_dp
2171 CALL dbcsr_add_on_diag(fy, alpha=r)
2172 CASE (3)
2173 ! C_out = C_in * ( 35/16 * I - 35/16 * P + 21/16 * P^2 - 5/16 P^3)
2174 CALL dbcsr_multiply('N', 'N', 1.0_dp, p, p, 0.0_dp, p2)
2175 IF (eps_irac_filter_matrix > 0.0_dp) THEN
2176 occ_in = dbcsr_get_occupation(p2)
2177 CALL dbcsr_filter(p2, eps_irac_filter_matrix)
2178 occ_out = dbcsr_get_occupation(p2)
2179 END IF
2180 r = -5.0_dp/16.0_dp
2181 CALL dbcsr_multiply('N', 'N', r, p2, p, 0.0_dp, fy)
2182 IF (eps_irac_filter_matrix > 0.0_dp) THEN
2183 occ_in = dbcsr_get_occupation(fy)
2184 CALL dbcsr_filter(fy, eps_irac_filter_matrix)
2185 occ_out = dbcsr_get_occupation(fy)
2186 END IF
2187 r = 21.0_dp/16.0_dp
2188 CALL dbcsr_add(fy, p2, alpha_scalar=1.0_dp, beta_scalar=r)
2189 r = -35.0_dp/16.0_dp
2190 CALL dbcsr_add(fy, p, alpha_scalar=1.0_dp, beta_scalar=r)
2191 r = 35.0_dp/16.0_dp
2192 CALL dbcsr_add_on_diag(fy, alpha=r)
2193 CASE (4)
2194 ! C_out = C_in * ( 315/128 * I - 420/128 * P + 378/128 * P^2 - 180/128 P^3 + 35/128 P^4 )
2195 ! = C_in * ( 315/128 * I - 420/128 * P + 378/128 * P^2 + ( - 180/128 * P + 35/128 * P^2 ) * P^2 )
2196 CALL dbcsr_multiply('N', 'N', 1.0_dp, p, p, 0.0_dp, p2) ! P^2
2197 IF (eps_irac_filter_matrix > 0.0_dp) THEN
2198 occ_in = dbcsr_get_occupation(p2)
2199 CALL dbcsr_filter(p2, eps_irac_filter_matrix)
2200 occ_out = dbcsr_get_occupation(p2)
2201 END IF
2202 r = -180.0_dp/128.0_dp
2203 CALL dbcsr_add(t, p, alpha_scalar=0.0_dp, beta_scalar=r) ! T=-180/128*P
2204 r = 35.0_dp/128.0_dp
2205 CALL dbcsr_add(t, p2, alpha_scalar=1.0_dp, beta_scalar=r) ! T=T+35/128*P^2
2206 CALL dbcsr_multiply('N', 'N', 1.0_dp, t, p2, 0.0_dp, fy) ! Y=T*P^2
2207 IF (eps_irac_filter_matrix > 0.0_dp) THEN
2208 occ_in = dbcsr_get_occupation(fy)
2209 CALL dbcsr_filter(fy, eps_irac_filter_matrix)
2210 occ_out = dbcsr_get_occupation(fy)
2211 END IF
2212 r = 378.0_dp/128.0_dp
2213 CALL dbcsr_add(fy, p2, alpha_scalar=1.0_dp, beta_scalar=r) ! Y=Y+378/128*P^2
2214 r = -420.0_dp/128.0_dp
2215 CALL dbcsr_add(fy, p, alpha_scalar=1.0_dp, beta_scalar=r) ! Y=Y-420/128*P
2216 r = 315.0_dp/128.0_dp
2217 CALL dbcsr_add_on_diag(fy, alpha=r) ! Y=Y+315/128*I
2218 CASE DEFAULT
2219 cpabort("This irac_order NYI")
2220 END SELECT
2221 CALL timestop(handle)
2222 END SUBROUTINE qs_ot_refine
2223
2224! **************************************************************************************************
2225!> \brief ...
2226!> \param matrix_hc ...
2227!> \param matrix_x ...
2228!> \param matrix_sx ...
2229!> \param matrix_gx ...
2230!> \param qs_ot_env ...
2231! **************************************************************************************************
2232 SUBROUTINE qs_ot_get_derivative_ref(matrix_hc, matrix_x, matrix_sx, matrix_gx, &
2233 qs_ot_env)
2234 TYPE(dbcsr_type), POINTER :: matrix_hc, matrix_x, matrix_sx, matrix_gx
2235 TYPE(qs_ot_type) :: qs_ot_env
2236
2237 CHARACTER(len=*), PARAMETER :: routinen = 'qs_ot_get_derivative_ref'
2238
2239 INTEGER :: handle, k, n
2240 REAL(dp) :: occ_in, occ_out
2241 TYPE(dbcsr_type), POINTER :: c, chc, g, hc, hc_work, sc
2242
2243 CALL timeset(routinen, handle)
2244
2245 CALL dbcsr_get_info(matrix_x, nfullrows_total=n, nfullcols_total=k)
2246 !
2247 c => matrix_x ! NBsf*NOcc
2248 sc => matrix_sx ! NBsf*NOcc need to be up2date
2249 hc => matrix_hc ! NBsf*NOcc
2250 g => matrix_gx ! NBsf*NOcc
2251 chc => qs_ot_env%buf1_k_k_sym ! buffer
2252
2253 IF (qs_ot_env%settings%do_rotation) THEN
2254 ! The physical orbitals are C_current=Q*U. Pull dE/dC_current
2255 ! back to the unrotated REF basis before projecting it.
2256 CALL qs_ot_rot_mat_derivative(qs_ot_env)
2257 hc_work => qs_ot_env%buf1_n_k
2258 CALL dbcsr_multiply('N', 'T', 1.0_dp, hc, qs_ot_env%rot_mat_u, &
2259 0.0_dp, hc_work)
2260 ELSE
2261 hc_work => hc
2262 END IF
2263
2264 ! C'*(H*C)
2265 CALL dbcsr_multiply('T', 'N', 1.0_dp, c, hc_work, 0.0_dp, chc)
2266 IF (qs_ot_env%settings%eps_irac_filter_matrix > 0.0_dp) THEN
2267 occ_in = dbcsr_get_occupation(chc)
2268 CALL dbcsr_filter(chc, qs_ot_env%settings%eps_irac_filter_matrix)
2269 occ_out = dbcsr_get_occupation(chc)
2270 END IF
2271 ! (S*C)*(C'*H*C)
2272 CALL dbcsr_multiply('N', 'N', 1.0_dp, sc, chc, 0.0_dp, g)
2273 IF (qs_ot_env%settings%eps_irac_filter_matrix > 0.0_dp) THEN
2274 occ_in = dbcsr_get_occupation(g)
2275 CALL dbcsr_filter(g, qs_ot_env%settings%eps_irac_filter_matrix)
2276 occ_out = dbcsr_get_occupation(g)
2277 END IF
2278 ! G = 2*(1-S*C*C')*H*C
2279 CALL dbcsr_add(g, hc_work, alpha_scalar=-1.0_dp, beta_scalar=1.0_dp)
2280 !
2281 CALL timestop(handle)
2282 END SUBROUTINE qs_ot_get_derivative_ref
2283
2284! **************************************************************************************************
2285!> \brief complex k-point REF derivative dE/dX from H(k)C(k), S(k)C(k), and C(k)
2286!> \param matrix_hc ...
2287!> \param matrix_hc_im ...
2288!> \param qs_ot_env ...
2289!> \param matrix_hc_rotation occupation-weighted H(k)C(k) for the rotation channel
2290!> \param matrix_hc_rotation_im imaginary component of matrix_hc_rotation
2291! **************************************************************************************************
2292 SUBROUTINE qs_ot_get_derivative_ref_complex(matrix_hc, matrix_hc_im, qs_ot_env, &
2293 matrix_hc_rotation, matrix_hc_rotation_im)
2294 TYPE(dbcsr_type), POINTER :: matrix_hc, matrix_hc_im
2295 TYPE(qs_ot_type) :: qs_ot_env
2296 TYPE(dbcsr_type), OPTIONAL, POINTER :: matrix_hc_rotation, matrix_hc_rotation_im
2297
2298 CHARACTER(len=*), PARAMETER :: routinen = 'qs_ot_get_derivative_ref_complex'
2299
2300 INTEGER :: handle, k, n
2301 REAL(dp) :: occ_in, occ_out
2302 TYPE(dbcsr_type) :: tmp_nk
2303 TYPE(dbcsr_type), POINTER :: b_im, b_re, f_im, f_re, g_im, g_re, hc_im, hc_re, &
2304 hc_rotation_im, hc_rotation_re, hc_work_im, hc_work_re, q_im, q_re, sc_im, sc_re, tmp_kk
2305 TYPE(dbcsr_type), TARGET :: hc_rot_im, hc_rot_re
2306
2307 CALL timeset(routinen, handle)
2308
2309 cpassert(qs_ot_env%has_complex_kpoint_state)
2310 cpassert(ASSOCIATED(matrix_hc))
2311 cpassert(ASSOCIATED(matrix_hc_im))
2312
2313 f_re => qs_ot_env%matrix_ref_inv_sqrt
2314 f_im => qs_ot_env%matrix_ref_inv_sqrt_im
2315 sc_re => qs_ot_env%matrix_sx
2316 sc_im => qs_ot_env%matrix_sx_im
2317 hc_re => matrix_hc
2318 hc_im => matrix_hc_im
2319 hc_rotation_re => hc_re
2320 hc_rotation_im => hc_im
2321 IF (PRESENT(matrix_hc_rotation) .OR. PRESENT(matrix_hc_rotation_im)) THEN
2322 cpassert(PRESENT(matrix_hc_rotation) .AND. PRESENT(matrix_hc_rotation_im))
2323 cpassert(ASSOCIATED(matrix_hc_rotation))
2324 cpassert(ASSOCIATED(matrix_hc_rotation_im))
2325 hc_rotation_re => matrix_hc_rotation
2326 hc_rotation_im => matrix_hc_rotation_im
2327 END IF
2328 hc_work_re => hc_re
2329 hc_work_im => hc_im
2330 g_re => qs_ot_env%matrix_gx
2331 g_im => qs_ot_env%matrix_gx_im
2332 b_re => qs_ot_env%buf1_k_k_sym
2333 b_im => qs_ot_env%buf2_k_k_sym
2334 tmp_kk => qs_ot_env%buf3_k_k_sym
2335 q_re => qs_ot_env%buf1_n_k
2336 q_im => qs_ot_env%buf1_n_k_dp
2337
2338 CALL dbcsr_get_info(sc_re, nfullrows_total=n, nfullcols_total=k)
2339
2340 ! Q = X*(X^H*S*X)^(-1/2), reconstructed from the current REF coordinate.
2341 CALL dbcsr_multiply('N', 'N', 1.0_dp, qs_ot_env%matrix_x, f_re, 0.0_dp, q_re)
2342 CALL dbcsr_multiply('N', 'N', 1.0_dp, qs_ot_env%matrix_x_im, f_im, 0.0_dp, g_re)
2343 CALL dbcsr_add(q_re, g_re, alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
2344
2345 CALL dbcsr_multiply('N', 'N', 1.0_dp, qs_ot_env%matrix_x, f_im, 0.0_dp, q_im)
2346 CALL dbcsr_multiply('N', 'N', 1.0_dp, qs_ot_env%matrix_x_im, f_re, 0.0_dp, g_re)
2347 CALL dbcsr_add(q_im, g_re, alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
2348
2349 IF (qs_ot_env%settings%do_rotation) THEN
2350 CALL qs_ot_generate_rotation_complex(qs_ot_env)
2351
2352 ! dF/dU = Q^H*G_C for C=Q*U.
2353 CALL dbcsr_multiply('T', 'N', 1.0_dp, q_re, hc_rotation_re, &
2354 0.0_dp, qs_ot_env%rot_mat_dedu)
2355 CALL dbcsr_multiply('T', 'N', 1.0_dp, q_im, hc_rotation_im, 0.0_dp, tmp_kk)
2356 CALL dbcsr_add(qs_ot_env%rot_mat_dedu, tmp_kk, &
2357 alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
2358
2359 CALL dbcsr_multiply('T', 'N', 1.0_dp, q_re, hc_rotation_im, &
2360 0.0_dp, qs_ot_env%rot_mat_dedu_im)
2361 CALL dbcsr_multiply('T', 'N', 1.0_dp, q_im, hc_rotation_re, 0.0_dp, tmp_kk)
2362 CALL dbcsr_add(qs_ot_env%rot_mat_dedu_im, tmp_kk, &
2363 alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
2364 CALL qs_ot_rot_mat_derivative_complex(qs_ot_env)
2365
2366 ! The REF coordinate sees G_Q=G_C*U^H.
2367 CALL dbcsr_copy(hc_rot_re, hc_re, name="hc_rot_re")
2368 CALL dbcsr_copy(hc_rot_im, hc_im, name="hc_rot_im")
2369 CALL dbcsr_copy(tmp_nk, hc_re, name="tmp_nk")
2370 CALL dbcsr_multiply('N', 'T', 1.0_dp, hc_re, qs_ot_env%rot_mat_u, &
2371 0.0_dp, hc_rot_re)
2372 CALL dbcsr_multiply('N', 'T', 1.0_dp, hc_im, qs_ot_env%rot_mat_u_im, &
2373 0.0_dp, tmp_nk)
2374 CALL dbcsr_add(hc_rot_re, tmp_nk, alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
2375
2376 CALL dbcsr_multiply('N', 'T', 1.0_dp, hc_im, qs_ot_env%rot_mat_u, &
2377 0.0_dp, hc_rot_im)
2378 CALL dbcsr_multiply('N', 'T', 1.0_dp, hc_re, qs_ot_env%rot_mat_u_im, &
2379 0.0_dp, tmp_nk)
2380 CALL dbcsr_add(hc_rot_im, tmp_nk, alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
2381 hc_work_re => hc_rot_re
2382 hc_work_im => hc_rot_im
2383 END IF
2384
2385 ! B = Q^H*G_Q. For uniform fixed occupations B is Hermitian.
2386 CALL dbcsr_multiply('T', 'N', 1.0_dp, q_re, hc_work_re, 0.0_dp, b_re)
2387 CALL dbcsr_multiply('T', 'N', 1.0_dp, q_im, hc_work_im, 0.0_dp, tmp_kk)
2388 CALL dbcsr_add(b_re, tmp_kk, alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
2389
2390 CALL dbcsr_multiply('T', 'N', 1.0_dp, q_re, hc_work_im, 0.0_dp, b_im)
2391 CALL dbcsr_multiply('T', 'N', 1.0_dp, q_im, hc_work_re, 0.0_dp, tmp_kk)
2392 CALL dbcsr_add(b_im, tmp_kk, alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
2393
2394 IF (qs_ot_env%settings%eps_irac_filter_matrix > 0.0_dp) THEN
2395 occ_in = dbcsr_get_occupation(b_re)
2396 CALL dbcsr_filter(b_re, qs_ot_env%settings%eps_irac_filter_matrix)
2397 occ_out = dbcsr_get_occupation(b_re)
2398 occ_in = dbcsr_get_occupation(b_im)
2399 CALL dbcsr_filter(b_im, qs_ot_env%settings%eps_irac_filter_matrix)
2400 occ_out = dbcsr_get_occupation(b_im)
2401 END IF
2402
2403 ! S*Q = (S*X)*F. G is temporary storage for this pair.
2404 CALL dbcsr_multiply('N', 'N', 1.0_dp, sc_re, f_re, 0.0_dp, g_re)
2405 CALL dbcsr_multiply('N', 'N', 1.0_dp, sc_im, f_im, 0.0_dp, q_re)
2406 CALL dbcsr_add(g_re, q_re, alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
2407
2408 CALL dbcsr_multiply('N', 'N', 1.0_dp, sc_re, f_im, 0.0_dp, g_im)
2409 CALL dbcsr_multiply('N', 'N', 1.0_dp, sc_im, f_re, 0.0_dp, q_re)
2410 CALL dbcsr_add(g_im, q_re, alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
2411
2412 ! Form (S*Q)*B. The Q workspaces are free after B has been formed.
2413 CALL dbcsr_multiply('N', 'N', 1.0_dp, g_re, b_re, 0.0_dp, q_re)
2414 CALL dbcsr_multiply('N', 'N', 1.0_dp, g_im, b_im, 0.0_dp, q_im)
2415 CALL dbcsr_add(q_re, q_im, alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
2416
2417 CALL dbcsr_multiply('N', 'N', 1.0_dp, g_re, b_im, 0.0_dp, q_im)
2418 CALL dbcsr_multiply('N', 'N', 1.0_dp, g_im, b_re, 0.0_dp, g_re)
2419 CALL dbcsr_add(q_im, g_re, alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
2420
2421 CALL dbcsr_add(q_re, hc_work_re, alpha_scalar=-1.0_dp, beta_scalar=1.0_dp)
2422 CALL dbcsr_add(q_im, hc_work_im, alpha_scalar=-1.0_dp, beta_scalar=1.0_dp)
2423
2424 ! Pull the physical gradient back to the finite REF coordinate: G_X = (G_C-S*Q*B)*F.
2425 CALL dbcsr_multiply('N', 'N', 1.0_dp, q_re, f_re, 0.0_dp, g_re)
2426 CALL dbcsr_multiply('N', 'N', 1.0_dp, q_im, f_im, 0.0_dp, g_im)
2427 CALL dbcsr_add(g_re, g_im, alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
2428
2429 CALL dbcsr_multiply('N', 'N', 1.0_dp, q_re, f_im, 0.0_dp, g_im)
2430 CALL dbcsr_multiply('N', 'N', 1.0_dp, q_im, f_re, 0.0_dp, q_re)
2431 CALL dbcsr_add(g_im, q_re, alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
2432
2433 IF (qs_ot_env%settings%do_rotation) THEN
2434 CALL qs_ot_add_ref_vertical_response_complex(b_re, b_im, f_re, f_im, &
2435 sc_re, sc_im, g_re, g_im, qs_ot_env)
2436 END IF
2437
2438 IF (qs_ot_env%settings%eps_irac_filter_matrix > 0.0_dp) THEN
2439 occ_in = dbcsr_get_occupation(g_re)
2440 CALL dbcsr_filter(g_re, qs_ot_env%settings%eps_irac_filter_matrix)
2441 occ_out = dbcsr_get_occupation(g_re)
2442 occ_in = dbcsr_get_occupation(g_im)
2443 CALL dbcsr_filter(g_im, qs_ot_env%settings%eps_irac_filter_matrix)
2444 occ_out = dbcsr_get_occupation(g_im)
2445 END IF
2446
2447 IF (qs_ot_env%settings%do_rotation) THEN
2448 CALL dbcsr_release(hc_rot_re)
2449 CALL dbcsr_release(hc_rot_im)
2450 CALL dbcsr_release(tmp_nk)
2451 END IF
2452
2453 CALL timestop(handle)
2454
2455 END SUBROUTINE qs_ot_get_derivative_ref_complex
2456
2457! **************************************************************************************************
2458!> \brief Transpose a square DBCSR matrix while retaining its target distribution.
2459!> \param matrix source matrix
2460!> \param transposed transposed matrix
2461!> \param identity_template square matrix template
2462! **************************************************************************************************
2463 SUBROUTINE qs_ot_square_transpose(matrix, transposed, identity_template)
2464 TYPE(dbcsr_type) :: matrix, transposed, identity_template
2465
2466 TYPE(dbcsr_type) :: identity
2467
2468 CALL dbcsr_copy(transposed, matrix, name='square_transposed')
2469 CALL dbcsr_copy(identity, identity_template, name='transpose_identity')
2470 CALL dbcsr_set(identity, 0.0_dp)
2471 CALL dbcsr_add_on_diag(identity, 1.0_dp)
2472 CALL dbcsr_multiply('T', 'N', 1.0_dp, matrix, identity, 0.0_dp, transposed)
2473 CALL dbcsr_release(identity)
2474
2475 END SUBROUTINE qs_ot_square_transpose
2476
2477! **************************************************************************************************
2478!> \brief Add the vertical response of the finite complex polar REF map.
2479!> \param a_re real part of Q^H G_Q
2480!> \param a_im imaginary part of Q^H G_Q
2481!> \param f_re real part of (X^H S X)^(-1/2)
2482!> \param f_im imaginary part of (X^H S X)^(-1/2)
2483!> \param sx_re real part of S X
2484!> \param sx_im imaginary part of S X
2485!> \param gradient_re real REF gradient, updated in place
2486!> \param gradient_im imaginary REF gradient, updated in place
2487!> \param qs_ot_env complex OT environment
2488! **************************************************************************************************
2489 SUBROUTINE qs_ot_add_ref_vertical_response_complex(a_re, a_im, f_re, f_im, &
2490 sx_re, sx_im, gradient_re, gradient_im, &
2491 qs_ot_env)
2492 TYPE(dbcsr_type) :: a_re, a_im, f_re, f_im, sx_re, sx_im, &
2493 gradient_re, gradient_im
2494 TYPE(qs_ot_type) :: qs_ot_env
2495
2496 CHARACTER(len=*), PARAMETER :: routinen = 'qs_ot_add_ref_vertical_response_complex'
2497
2498 INTEGER :: col, handle, i, j, k, row
2499 INTEGER, DIMENSION(:), POINTER :: col_blk_offset, row_blk_offset
2500 LOGICAL :: found_im, found_re
2501 REAL(kind=dp) :: denominator
2502 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: f_evals
2503 REAL(kind=dp), DIMENSION(:, :), POINTER :: block_im, block_re
2504 TYPE(dbcsr_iterator_type) :: iter
2505 TYPE(dbcsr_type) :: anti_im, anti_re, f_work_im, f_work_re, inner_im, inner_re, sq_im, &
2506 sq_re, tmp_kk, tmp_nk, v_im, v_re, vertical_im, vertical_re, work_im, work_re, z_im, z_re
2507
2508 CALL timeset(routinen, handle)
2509
2510 CALL dbcsr_get_info(a_re, nfullrows_total=k)
2511 IF (k == 0) THEN
2512 CALL timestop(handle)
2513 RETURN
2514 END IF
2515
2516 ! Only the anti-Hermitian part of A=Q^H G_Q couples to the unitary
2517 ! component of the polar differential.
2518 CALL qs_ot_square_transpose(a_re, tmp_kk, qs_ot_env%rot_mat_u)
2519 CALL dbcsr_copy(anti_re, a_re, name='anti_re')
2520 CALL dbcsr_add(anti_re, tmp_kk, alpha_scalar=0.5_dp, beta_scalar=-0.5_dp)
2521 CALL qs_ot_square_transpose(a_im, tmp_kk, qs_ot_env%rot_mat_u)
2522 CALL dbcsr_copy(anti_im, a_im, name='anti_im')
2523 CALL dbcsr_add(anti_im, tmp_kk, alpha_scalar=0.5_dp, beta_scalar=0.5_dp)
2524
2525 ! P=(X^H S X)^(1/2)=F^(-1). In the eigenbasis of F, solve
2526 ! P Z + Z P = antiherm(A) element by element.
2527 CALL dbcsr_copy(f_work_re, f_re, name='f_work_re')
2528 CALL dbcsr_copy(f_work_im, f_im, name='f_work_im')
2529 CALL dbcsr_copy(v_re, f_re, name='v_re')
2530 CALL dbcsr_copy(v_im, f_im, name='v_im')
2531 ALLOCATE (f_evals(k))
2532 CALL cp_dbcsr_heevd(matrix_re=f_work_re, matrix_im=f_work_im, &
2533 eigenvectors_re=v_re, eigenvectors_im=v_im, eigenvalues=f_evals, &
2534 para_env=qs_ot_env%para_env, blacs_env=qs_ot_env%blacs_env)
2535 IF (minval(f_evals) <= epsilon(1.0_dp)) THEN
2536 cpabort('Complex REF inverse square root is not positive definite')
2537 END IF
2538
2539 ! inner = V^H antiherm(A) V.
2540 CALL dbcsr_copy(work_re, anti_re, name='work_re')
2541 CALL dbcsr_copy(work_im, anti_im, name='work_im')
2542 CALL dbcsr_multiply('N', 'N', 1.0_dp, anti_re, v_re, 0.0_dp, work_re)
2543 CALL dbcsr_multiply('N', 'N', 1.0_dp, anti_im, v_im, 0.0_dp, tmp_kk)
2544 CALL dbcsr_add(work_re, tmp_kk, alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
2545 CALL dbcsr_multiply('N', 'N', 1.0_dp, anti_re, v_im, 0.0_dp, work_im)
2546 CALL dbcsr_multiply('N', 'N', 1.0_dp, anti_im, v_re, 0.0_dp, tmp_kk)
2547 CALL dbcsr_add(work_im, tmp_kk, alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
2548
2549 CALL dbcsr_copy(inner_re, anti_re, name='inner_re')
2550 CALL dbcsr_copy(inner_im, anti_im, name='inner_im')
2551 CALL dbcsr_multiply('T', 'N', 1.0_dp, v_re, work_re, 0.0_dp, inner_re)
2552 CALL dbcsr_multiply('T', 'N', 1.0_dp, v_im, work_im, 1.0_dp, inner_re)
2553 CALL dbcsr_multiply('T', 'N', 1.0_dp, v_re, work_im, 0.0_dp, inner_im)
2554 CALL dbcsr_multiply('T', 'N', -1.0_dp, v_im, work_re, 1.0_dp, inner_im)
2555
2556 CALL dbcsr_get_info(inner_re, row_blk_offset=row_blk_offset, &
2557 col_blk_offset=col_blk_offset)
2558 CALL dbcsr_iterator_start(iter, inner_re)
2559 DO WHILE (dbcsr_iterator_blocks_left(iter))
2560 CALL dbcsr_iterator_next_block(iter, row, col)
2561 CALL dbcsr_get_block_p(inner_re, row, col, block_re, found_re)
2562 CALL dbcsr_get_block_p(inner_im, row, col, block_im, found_im)
2563 cpassert(found_re .AND. found_im)
2564 DO i = 1, SIZE(block_re, 1)
2565 DO j = 1, SIZE(block_re, 2)
2566 denominator = 1.0_dp/f_evals(row_blk_offset(row) + i - 1) + &
2567 1.0_dp/f_evals(col_blk_offset(col) + j - 1)
2568 block_re(i, j) = block_re(i, j)/denominator
2569 block_im(i, j) = block_im(i, j)/denominator
2570 END DO
2571 END DO
2572 END DO
2573 CALL dbcsr_iterator_stop(iter)
2574
2575 ! Z = V inner V^H.
2576 CALL dbcsr_multiply('N', 'N', 1.0_dp, v_re, inner_re, 0.0_dp, work_re)
2577 CALL dbcsr_multiply('N', 'N', 1.0_dp, v_im, inner_im, 0.0_dp, tmp_kk)
2578 CALL dbcsr_add(work_re, tmp_kk, alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
2579 CALL dbcsr_multiply('N', 'N', 1.0_dp, v_re, inner_im, 0.0_dp, work_im)
2580 CALL dbcsr_multiply('N', 'N', 1.0_dp, v_im, inner_re, 0.0_dp, tmp_kk)
2581 CALL dbcsr_add(work_im, tmp_kk, alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
2582
2583 CALL dbcsr_copy(z_re, anti_re, name='z_re')
2584 CALL dbcsr_copy(z_im, anti_im, name='z_im')
2585 CALL dbcsr_multiply('N', 'T', 1.0_dp, work_re, v_re, 0.0_dp, z_re)
2586 CALL dbcsr_multiply('N', 'T', 1.0_dp, work_im, v_im, 1.0_dp, z_re)
2587 CALL dbcsr_multiply('N', 'T', 1.0_dp, work_im, v_re, 0.0_dp, z_im)
2588 CALL dbcsr_multiply('N', 'T', -1.0_dp, work_re, v_im, 1.0_dp, z_im)
2589
2590 ! The adjoint vertical contribution is 2 S Q Z, with S Q=(S X)F.
2591 CALL dbcsr_copy(sq_re, sx_re, name='sq_re')
2592 CALL dbcsr_copy(sq_im, sx_im, name='sq_im')
2593 CALL dbcsr_copy(tmp_nk, sx_re, name='tmp_nk')
2594 CALL dbcsr_multiply('N', 'N', 1.0_dp, sx_re, f_re, 0.0_dp, sq_re)
2595 CALL dbcsr_multiply('N', 'N', 1.0_dp, sx_im, f_im, 0.0_dp, tmp_nk)
2596 CALL dbcsr_add(sq_re, tmp_nk, alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
2597 CALL dbcsr_multiply('N', 'N', 1.0_dp, sx_re, f_im, 0.0_dp, sq_im)
2598 CALL dbcsr_multiply('N', 'N', 1.0_dp, sx_im, f_re, 0.0_dp, tmp_nk)
2599 CALL dbcsr_add(sq_im, tmp_nk, alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
2600
2601 CALL dbcsr_copy(vertical_re, sx_re, name='vertical_re')
2602 CALL dbcsr_copy(vertical_im, sx_im, name='vertical_im')
2603 CALL dbcsr_multiply('N', 'N', 1.0_dp, sq_re, z_re, 0.0_dp, vertical_re)
2604 CALL dbcsr_multiply('N', 'N', 1.0_dp, sq_im, z_im, 0.0_dp, tmp_nk)
2605 CALL dbcsr_add(vertical_re, tmp_nk, alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
2606 CALL dbcsr_multiply('N', 'N', 1.0_dp, sq_re, z_im, 0.0_dp, vertical_im)
2607 CALL dbcsr_multiply('N', 'N', 1.0_dp, sq_im, z_re, 0.0_dp, tmp_nk)
2608 CALL dbcsr_add(vertical_im, tmp_nk, alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
2609 CALL dbcsr_add(gradient_re, vertical_re, alpha_scalar=1.0_dp, beta_scalar=2.0_dp)
2610 CALL dbcsr_add(gradient_im, vertical_im, alpha_scalar=1.0_dp, beta_scalar=2.0_dp)
2611
2612 DEALLOCATE (f_evals)
2613 CALL dbcsr_release(anti_im)
2614 CALL dbcsr_release(anti_re)
2615 CALL dbcsr_release(f_work_im)
2616 CALL dbcsr_release(f_work_re)
2617 CALL dbcsr_release(inner_im)
2618 CALL dbcsr_release(inner_re)
2619 CALL dbcsr_release(sq_im)
2620 CALL dbcsr_release(sq_re)
2621 CALL dbcsr_release(tmp_kk)
2622 CALL dbcsr_release(tmp_nk)
2623 CALL dbcsr_release(v_im)
2624 CALL dbcsr_release(v_re)
2625 CALL dbcsr_release(vertical_im)
2626 CALL dbcsr_release(vertical_re)
2627 CALL dbcsr_release(work_im)
2628 CALL dbcsr_release(work_re)
2629 CALL dbcsr_release(z_im)
2630 CALL dbcsr_release(z_re)
2631
2632 CALL timestop(handle)
2633
2634 END SUBROUTINE qs_ot_add_ref_vertical_response_complex
2635
2636! **************************************************************************************************
2637!> \brief computes p=x*S*x and the matrix functionals related matrices
2638!> \param matrix_x ...
2639!> \param matrix_sx ...
2640!> \param qs_ot_env ...
2641! **************************************************************************************************
2642 SUBROUTINE qs_ot_get_p(matrix_x, matrix_sx, qs_ot_env)
2643
2644 TYPE(dbcsr_type), POINTER :: matrix_x, matrix_sx
2645 TYPE(qs_ot_type) :: qs_ot_env
2646
2647 CHARACTER(len=*), PARAMETER :: routinen = 'qs_ot_get_p'
2648 REAL(kind=dp), PARAMETER :: rone = 1.0_dp, rzero = 0.0_dp
2649
2650 INTEGER :: handle, k, max_iter, n
2651 LOGICAL :: converged
2652 REAL(kind=dp) :: max_ev, min_ev, threshold
2653
2654 CALL timeset(routinen, handle)
2655
2656 CALL dbcsr_get_info(matrix_x, nfullrows_total=n, nfullcols_total=k)
2657
2658 ! get the overlap
2659 CALL dbcsr_multiply('T', 'N', rone, matrix_x, matrix_sx, rzero, &
2660 qs_ot_env%matrix_p)
2661
2662 ! get an upper bound for the largest eigenvalue
2663 ! try using lancos first and fall back to gershgorin norm if it fails
2664 max_iter = 30; threshold = 1.0e-03_dp
2665 CALL arnoldi_extremal(qs_ot_env%matrix_p, max_ev, min_ev, converged, threshold, max_iter)
2666 qs_ot_env%largest_eval_upper_bound = max(max_ev, abs(min_ev))
2667
2668 IF (.NOT. converged) qs_ot_env%largest_eval_upper_bound = dbcsr_gershgorin_norm(qs_ot_env%matrix_p)
2669 CALL decide_strategy(qs_ot_env)
2670 IF (qs_ot_env%do_taylor) THEN
2671 CALL qs_ot_p2m_taylor(qs_ot_env)
2672 ELSE
2673 CALL qs_ot_p2m_diag(qs_ot_env)
2674 END IF
2675
2676 IF (qs_ot_env%settings%do_rotation) THEN
2677 CALL qs_ot_generate_rotation(qs_ot_env)
2678 END IF
2679
2680 CALL timestop(handle)
2681
2682 END SUBROUTINE qs_ot_get_p
2683
2684! **************************************************************************************************
2685!> \brief computes U=exp(A) for the complex anti-Hermitian generator
2686!> A=rot_mat_x+i*rot_mat_x_im
2687!> \param qs_ot_env a complex k-point OT environment
2688! **************************************************************************************************
2689 SUBROUTINE qs_ot_generate_rotation_complex(qs_ot_env)
2690
2691 TYPE(qs_ot_type) :: qs_ot_env
2692
2693 CHARACTER(len=*), PARAMETER :: routinen = 'qs_ot_generate_rotation_complex'
2694
2695 INTEGER :: handle, k
2696 REAL(kind=dp) :: rot_norm
2697 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: exp_evals_im, exp_evals_re
2698 TYPE(dbcsr_type) :: h_re, tmp, w_im, w_re
2699
2700 CALL timeset(routinen, handle)
2701
2702 cpassert(qs_ot_env%has_complex_kpoint_state)
2703 cpassert(ASSOCIATED(qs_ot_env%rot_mat_x_im))
2704 cpassert(ASSOCIATED(qs_ot_env%rot_mat_u_im))
2705
2706 CALL dbcsr_get_info(qs_ot_env%rot_mat_x, nfullrows_total=k)
2707 IF (k /= 0) THEN
2708 rot_norm = dbcsr_frobenius_norm(qs_ot_env%rot_mat_x) + &
2709 dbcsr_frobenius_norm(qs_ot_env%rot_mat_x_im)
2710 IF (rot_norm <= epsilon(1.0_dp)) THEN
2711 CALL dbcsr_set(qs_ot_env%rot_mat_u, 0.0_dp)
2712 CALL dbcsr_add_on_diag(qs_ot_env%rot_mat_u, 1.0_dp)
2713 CALL dbcsr_set(qs_ot_env%rot_mat_u_im, 0.0_dp)
2714 CALL timestop(handle)
2715 RETURN
2716 END IF
2717
2718 ! i*A = i*X-Y is Hermitian. Its eigenvectors give
2719 ! exp(A)=V*diag(exp(-i*lambda))*V^H.
2720 CALL dbcsr_copy(h_re, qs_ot_env%rot_mat_x_im, name="h_re")
2721 CALL dbcsr_scale(h_re, -1.0_dp)
2722 CALL cp_dbcsr_heevd(matrix_re=h_re, matrix_im=qs_ot_env%rot_mat_x, &
2723 eigenvectors_re=qs_ot_env%rot_mat_evec_re, &
2724 eigenvectors_im=qs_ot_env%rot_mat_evec_im, &
2725 eigenvalues=qs_ot_env%rot_mat_evals, &
2726 para_env=qs_ot_env%para_env, blacs_env=qs_ot_env%blacs_env)
2727
2728 ALLOCATE (exp_evals_re(k), exp_evals_im(k))
2729 exp_evals_re(:) = cos(-qs_ot_env%rot_mat_evals(:))
2730 exp_evals_im(:) = sin(-qs_ot_env%rot_mat_evals(:))
2731
2732 CALL dbcsr_copy(w_re, qs_ot_env%rot_mat_evec_re, name="w_re")
2733 CALL dbcsr_scale_by_vector(w_re, alpha=exp_evals_re, side='right')
2734 CALL dbcsr_copy(tmp, qs_ot_env%rot_mat_evec_im, name="tmp")
2735 CALL dbcsr_scale_by_vector(tmp, alpha=exp_evals_im, side='right')
2736 CALL dbcsr_add(w_re, tmp, alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
2737
2738 CALL dbcsr_copy(w_im, qs_ot_env%rot_mat_evec_im, name="w_im")
2739 CALL dbcsr_scale_by_vector(w_im, alpha=exp_evals_re, side='right')
2740 CALL dbcsr_copy(tmp, qs_ot_env%rot_mat_evec_re)
2741 CALL dbcsr_scale_by_vector(tmp, alpha=exp_evals_im, side='right')
2742 CALL dbcsr_add(w_im, tmp, alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
2743
2744 CALL dbcsr_multiply('N', 'T', 1.0_dp, w_re, qs_ot_env%rot_mat_evec_re, &
2745 0.0_dp, qs_ot_env%rot_mat_u)
2746 CALL dbcsr_multiply('N', 'T', 1.0_dp, w_im, qs_ot_env%rot_mat_evec_im, &
2747 1.0_dp, qs_ot_env%rot_mat_u)
2748 CALL dbcsr_multiply('N', 'T', 1.0_dp, w_im, qs_ot_env%rot_mat_evec_re, &
2749 0.0_dp, qs_ot_env%rot_mat_u_im)
2750 CALL dbcsr_multiply('N', 'T', -1.0_dp, w_re, qs_ot_env%rot_mat_evec_im, &
2751 1.0_dp, qs_ot_env%rot_mat_u_im)
2752
2753 CALL dbcsr_release(h_re)
2754 CALL dbcsr_release(tmp)
2755 CALL dbcsr_release(w_re)
2756 CALL dbcsr_release(w_im)
2757 DEALLOCATE (exp_evals_re, exp_evals_im)
2758 END IF
2759
2760 CALL timestop(handle)
2761
2762 END SUBROUTINE qs_ot_generate_rotation_complex
2763
2764! **************************************************************************************************
2765!> \brief pull the complex dE/dU covector back to the anti-Hermitian generator
2766!> using the adjoint Frechet derivative of exp
2767!> \param qs_ot_env a complex k-point OT environment with an up-to-date U
2768! **************************************************************************************************
2769 SUBROUTINE qs_ot_rot_mat_derivative_complex(qs_ot_env)
2770 TYPE(qs_ot_type) :: qs_ot_env
2771
2772 CHARACTER(len=*), PARAMETER :: routinen = 'qs_ot_rot_mat_derivative_complex'
2773
2774 INTEGER :: handle, k
2775 REAL(kind=dp) :: rot_norm
2776 TYPE(dbcsr_type) :: frechet_im, frechet_re, inner_deriv_im, &
2777 inner_deriv_re, outer_deriv_im, &
2778 outer_deriv_re, tmp, work_im, work_re
2779
2780 CALL timeset(routinen, handle)
2781
2782 cpassert(qs_ot_env%has_complex_kpoint_state)
2783 cpassert(ASSOCIATED(qs_ot_env%rot_mat_dedu_im))
2784 cpassert(ASSOCIATED(qs_ot_env%rot_mat_gx_im))
2785
2786 CALL dbcsr_get_info(qs_ot_env%rot_mat_u, nfullrows_total=k)
2787 IF (k /= 0) THEN
2788 rot_norm = dbcsr_frobenius_norm(qs_ot_env%rot_mat_x) + &
2789 dbcsr_frobenius_norm(qs_ot_env%rot_mat_x_im)
2790 IF (rot_norm <= epsilon(1.0_dp)) THEN
2791 CALL qs_ot_square_transpose(qs_ot_env%rot_mat_dedu, tmp, qs_ot_env%rot_mat_u)
2792 CALL dbcsr_copy(qs_ot_env%rot_mat_gx, qs_ot_env%rot_mat_dedu)
2793 CALL dbcsr_add(qs_ot_env%rot_mat_gx, tmp, &
2794 alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
2795 CALL dbcsr_release(tmp)
2796
2797 CALL qs_ot_square_transpose(qs_ot_env%rot_mat_dedu_im, tmp, qs_ot_env%rot_mat_u)
2798 CALL dbcsr_copy(qs_ot_env%rot_mat_gx_im, qs_ot_env%rot_mat_dedu_im)
2799 CALL dbcsr_add(qs_ot_env%rot_mat_gx_im, tmp, &
2800 alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
2801 CALL dbcsr_release(tmp)
2802 CALL timestop(handle)
2803 RETURN
2804 END IF
2805
2806 CALL dbcsr_copy(work_re, qs_ot_env%rot_mat_dedu, name="work_re")
2807 CALL dbcsr_copy(work_im, qs_ot_env%rot_mat_dedu, name="work_im")
2808 CALL dbcsr_copy(tmp, qs_ot_env%rot_mat_dedu, name="tmp")
2809 CALL dbcsr_copy(inner_deriv_re, qs_ot_env%rot_mat_dedu, name="inner_deriv_re")
2810 CALL dbcsr_copy(inner_deriv_im, qs_ot_env%rot_mat_dedu, name="inner_deriv_im")
2811
2812 ! V^H*(dE/dU)*V, split into real and imaginary parts.
2813 CALL dbcsr_multiply('N', 'N', 1.0_dp, qs_ot_env%rot_mat_dedu, &
2814 qs_ot_env%rot_mat_evec_re, 0.0_dp, work_re)
2815 CALL dbcsr_multiply('N', 'N', 1.0_dp, qs_ot_env%rot_mat_dedu_im, &
2816 qs_ot_env%rot_mat_evec_im, 0.0_dp, tmp)
2817 CALL dbcsr_add(work_re, tmp, alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
2818 CALL dbcsr_multiply('N', 'N', 1.0_dp, qs_ot_env%rot_mat_dedu, &
2819 qs_ot_env%rot_mat_evec_im, 0.0_dp, work_im)
2820 CALL dbcsr_multiply('N', 'N', 1.0_dp, qs_ot_env%rot_mat_dedu_im, &
2821 qs_ot_env%rot_mat_evec_re, 0.0_dp, tmp)
2822 CALL dbcsr_add(work_im, tmp, alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
2823 CALL dbcsr_multiply('T', 'N', 1.0_dp, qs_ot_env%rot_mat_evec_re, &
2824 work_re, 0.0_dp, inner_deriv_re)
2825 CALL dbcsr_multiply('T', 'N', 1.0_dp, qs_ot_env%rot_mat_evec_im, &
2826 work_im, 1.0_dp, inner_deriv_re)
2827 CALL dbcsr_multiply('T', 'N', 1.0_dp, qs_ot_env%rot_mat_evec_re, &
2828 work_im, 0.0_dp, inner_deriv_im)
2829 CALL dbcsr_multiply('T', 'N', -1.0_dp, qs_ot_env%rot_mat_evec_im, &
2830 work_re, 1.0_dp, inner_deriv_im)
2831
2832 CALL qs_ot_apply_complex_frechet_dbcsr(qs_ot_env%rot_mat_evals, &
2833 inner_deriv_re, inner_deriv_im, &
2834 outer_deriv_re, outer_deriv_im, adjoint=.true.)
2835
2836 CALL dbcsr_multiply('N', 'N', 1.0_dp, qs_ot_env%rot_mat_evec_re, &
2837 outer_deriv_re, 0.0_dp, work_re)
2838 CALL dbcsr_multiply('N', 'N', 1.0_dp, qs_ot_env%rot_mat_evec_im, &
2839 outer_deriv_im, 0.0_dp, tmp)
2840 CALL dbcsr_add(work_re, tmp, alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
2841 CALL dbcsr_multiply('N', 'N', 1.0_dp, qs_ot_env%rot_mat_evec_re, &
2842 outer_deriv_im, 0.0_dp, work_im)
2843 CALL dbcsr_multiply('N', 'N', 1.0_dp, qs_ot_env%rot_mat_evec_im, &
2844 outer_deriv_re, 0.0_dp, tmp)
2845 CALL dbcsr_add(work_im, tmp, alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
2846
2847 CALL dbcsr_copy(frechet_re, qs_ot_env%rot_mat_dedu, name="frechet_re")
2848 CALL dbcsr_copy(frechet_im, qs_ot_env%rot_mat_dedu, name="frechet_im")
2849 CALL dbcsr_multiply('N', 'T', 1.0_dp, work_re, qs_ot_env%rot_mat_evec_re, &
2850 0.0_dp, frechet_re)
2851 CALL dbcsr_multiply('N', 'T', 1.0_dp, work_im, qs_ot_env%rot_mat_evec_im, &
2852 1.0_dp, frechet_re)
2853 CALL dbcsr_multiply('N', 'T', 1.0_dp, work_im, qs_ot_env%rot_mat_evec_re, &
2854 0.0_dp, frechet_im)
2855 CALL dbcsr_multiply('N', 'T', -1.0_dp, work_re, qs_ot_env%rot_mat_evec_im, &
2856 1.0_dp, frechet_im)
2857
2858 ! Tangents satisfy X^T=-X and Y^T=Y for A=X+iY.
2859 CALL qs_ot_square_transpose(frechet_re, tmp, qs_ot_env%rot_mat_u)
2860 CALL dbcsr_copy(qs_ot_env%rot_mat_gx, frechet_re)
2861 CALL dbcsr_add(qs_ot_env%rot_mat_gx, tmp, &
2862 alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
2863 CALL qs_ot_square_transpose(frechet_im, tmp, qs_ot_env%rot_mat_u)
2864 CALL dbcsr_copy(qs_ot_env%rot_mat_gx_im, frechet_im)
2865 CALL dbcsr_add(qs_ot_env%rot_mat_gx_im, tmp, &
2866 alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
2867
2868 CALL dbcsr_release(frechet_re)
2869 CALL dbcsr_release(frechet_im)
2870 CALL dbcsr_release(inner_deriv_re)
2871 CALL dbcsr_release(inner_deriv_im)
2872 CALL dbcsr_release(outer_deriv_re)
2873 CALL dbcsr_release(outer_deriv_im)
2874 CALL dbcsr_release(tmp)
2875 CALL dbcsr_release(work_re)
2876 CALL dbcsr_release(work_im)
2877 END IF
2878
2879 CALL timestop(handle)
2880
2881 END SUBROUTINE qs_ot_rot_mat_derivative_complex
2882
2883! **************************************************************************************************
2884!> \brief compute P=X^H*S*X and the STRICT matrix functions for a complex K-point channel
2885!> \param matrix_x real part of X
2886!> \param matrix_x_im imaginary part of X
2887!> \param matrix_sx real part of S*X
2888!> \param matrix_sx_im imaginary part of S*X
2889!> \param qs_ot_env OT channel state
2890! **************************************************************************************************
2891 SUBROUTINE qs_ot_get_p_complex(matrix_x, matrix_x_im, matrix_sx, matrix_sx_im, qs_ot_env)
2892 TYPE(dbcsr_type), POINTER :: matrix_x, matrix_x_im, matrix_sx, &
2893 matrix_sx_im
2894 TYPE(qs_ot_type) :: qs_ot_env
2895
2896 CHARACTER(len=*), PARAMETER :: routinen = 'qs_ot_get_p_complex'
2897
2898 INTEGER :: handle
2899
2900 CALL timeset(routinen, handle)
2901
2902 cpassert(qs_ot_env%has_complex_kpoint_state)
2903 cpassert(qs_ot_env%settings%ot_algorithm == 'TOD')
2904
2905 CALL qs_ot_complex_multiply('C', 'N', matrix_x, matrix_x_im, matrix_sx, matrix_sx_im, &
2906 qs_ot_env%matrix_p, qs_ot_env%matrix_p_im, qs_ot_env%matrix_buf1)
2907 qs_ot_env%do_taylor = .false.
2908 CALL qs_ot_p2m_diag_complex(qs_ot_env)
2909
2910 CALL timestop(handle)
2911
2912 END SUBROUTINE qs_ot_get_p_complex
2913
2914! **************************************************************************************************
2915!> \brief computes the rotation matrix rot_mat_u that is associated to a given
2916!> rot_mat_x using rot_mat_u=exp(rot_mat_x)
2917!> \param qs_ot_env a valid qs_ot_env
2918!> \par History
2919!> 08.2004 created [Joost VandeVondele]
2920!> 12.2024 Rewrite to use only real matrices [Ole Schuett]
2921! **************************************************************************************************
2922 SUBROUTINE qs_ot_generate_rotation(qs_ot_env)
2923
2924 TYPE(qs_ot_type) :: qs_ot_env
2925
2926 CHARACTER(len=*), PARAMETER :: routinen = 'qs_ot_generate_rotation'
2927
2928 INTEGER :: handle, k
2929 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: exp_evals_im, exp_evals_re
2930 TYPE(dbcsr_type) :: buf_1, buf_2
2931
2932 CALL timeset(routinen, handle)
2933
2934 CALL dbcsr_get_info(qs_ot_env%rot_mat_x, nfullrows_total=k)
2935
2936 IF (k /= 0) THEN
2937 ! We want to compute: rot_mat_u = exp(i*rot_mat_x)
2938
2939 ! Diagonalize: matrix = i*rot_mat_x.
2940 ! Note that matrix is imaginary and hermitian because rot_mat_x is real and anti-symmetric.
2941 CALL cp_dbcsr_heevd(matrix_im=qs_ot_env%rot_mat_x, & ! matrix_re omitted because it's zero
2942 eigenvectors_re=qs_ot_env%rot_mat_evec_re, &
2943 eigenvectors_im=qs_ot_env%rot_mat_evec_im, &
2944 eigenvalues=qs_ot_env%rot_mat_evals, &
2945 para_env=qs_ot_env%para_env, &
2946 blacs_env=qs_ot_env%blacs_env)
2947
2948 ! Compute: exp_evals = EXP(-i*rot_mat_evals)
2949 ALLOCATE (exp_evals_re(k), exp_evals_im(k))
2950 exp_evals_re(:) = cos(-qs_ot_env%rot_mat_evals(:))
2951 exp_evals_im(:) = sin(-qs_ot_env%rot_mat_evals(:))
2952
2953 ! Compute: rot_mat_u = \sum_ij exp_evals_ij * |rot_mat_evec_i> <rot_mat_evec_j|
2954 ! Note that we need only two matrix multiplications because rot_mat_u is real.
2955 CALL dbcsr_copy(buf_1, qs_ot_env%rot_mat_evec_re, name="buf_1")
2956 CALL dbcsr_scale_by_vector(buf_1, alpha=exp_evals_re, side='right')
2957 CALL dbcsr_copy(buf_2, qs_ot_env%rot_mat_evec_im, name="buf_2")
2958 CALL dbcsr_scale_by_vector(buf_2, alpha=exp_evals_im, side='right')
2959 CALL dbcsr_add(buf_1, buf_2, alpha_scalar=+1.0_dp, beta_scalar=-1.0_dp)
2960 CALL dbcsr_multiply('N', 'T', 1.0_dp, buf_1, qs_ot_env%rot_mat_evec_re, 0.0_dp, qs_ot_env%rot_mat_u)
2961
2962 CALL dbcsr_copy(buf_1, qs_ot_env%rot_mat_evec_im)
2963 CALL dbcsr_scale_by_vector(buf_1, alpha=exp_evals_re, side='right')
2964 CALL dbcsr_copy(buf_2, qs_ot_env%rot_mat_evec_re)
2965 CALL dbcsr_scale_by_vector(buf_2, alpha=exp_evals_im, side='right')
2966 CALL dbcsr_add(buf_1, buf_2, alpha_scalar=+1.0_dp, beta_scalar=+1.0_dp)
2967 CALL dbcsr_multiply('N', 'T', 1.0_dp, buf_1, qs_ot_env%rot_mat_evec_im, 1.0_dp, qs_ot_env%rot_mat_u)
2968
2969 ! Clean up.
2970 CALL dbcsr_release(buf_1)
2971 CALL dbcsr_release(buf_2)
2972 DEALLOCATE (exp_evals_re, exp_evals_im)
2973 END IF
2974
2975 CALL timestop(handle)
2976
2977 END SUBROUTINE qs_ot_generate_rotation
2978
2979! **************************************************************************************************
2980!> \brief computes the derivative fields with respect to rot_mat_x
2981!> \param qs_ot_env valid qs_ot_env. In particular qs_ot_generate_rotation has to be called before
2982!> and the rot_mat_dedu matrix has to be up to date
2983!> \par History
2984!> 08.2004 created [ Joost VandeVondele ]
2985!> 12.2024 Rewrite to use only real matrices [Ole Schuett]
2986! **************************************************************************************************
2987 SUBROUTINE qs_ot_rot_mat_derivative(qs_ot_env)
2988 TYPE(qs_ot_type) :: qs_ot_env
2989
2990 CHARACTER(len=*), PARAMETER :: routinen = 'qs_ot_rot_mat_derivative'
2991
2992 INTEGER :: col, handle, i, iblock, j, k, max_blocks, nblocks, row
2993 INTEGER, ALLOCATABLE, DIMENSION(:) :: cols, rows
2994 INTEGER, DIMENSION(:), POINTER :: col_blk_offset, col_blk_size, row_blk_offset, row_blk_size
2995 REAL(kind=dp) :: e1, e2
2996 TYPE(dbcsr_type) :: outer_deriv_re, outer_deriv_im, mat_buf, &
2997 inner_deriv_re, inner_deriv_im
2998 TYPE(dbcsr_distribution_type) :: dist
2999 TYPE(dbcsr_iterator_type) :: iter
3000 REAL(dp), DIMENSION(:, :), POINTER :: block_in_re, block_in_im, block_out_re, block_out_im
3001 LOGICAL :: duplicate, found_in_im, found_in_re, found_out_im, &
3002 found_out_re
3003 REAL(kind=dp) :: im_part, re_part
3004 COMPLEX(dp) :: cval_in, cval_out
3005 CALL timeset(routinen, handle)
3006
3007 CALL dbcsr_get_info(qs_ot_env%rot_mat_u, nfullrows_total=k)
3008 IF (k /= 0) THEN
3009 CALL dbcsr_copy(qs_ot_env%matrix_buf1, qs_ot_env%rot_mat_dedu)
3010 ! now we get to the derivative wrt the antisymmetric matrix rot_mat_x
3011 CALL dbcsr_copy(mat_buf, qs_ot_env%rot_mat_dedu, "mat_buf")
3012
3013 ! inner_deriv_ij = <rot_mat_evec_i| rot_mat_dedu |rot_mat_evec_j>
3014 CALL dbcsr_copy(inner_deriv_re, qs_ot_env%rot_mat_dedu, "inner_deriv_re") ! TODO just create
3015 CALL dbcsr_copy(inner_deriv_im, qs_ot_env%rot_mat_dedu, "inner_deriv_im") ! TODO just create
3016
3017 CALL dbcsr_multiply('T', 'N', +1.0_dp, qs_ot_env%rot_mat_dedu, qs_ot_env%rot_mat_evec_im, 0.0_dp, mat_buf)
3018 CALL dbcsr_multiply('T', 'N', +1.0_dp, qs_ot_env%rot_mat_evec_im, mat_buf, 0.0_dp, inner_deriv_re)
3019 CALL dbcsr_multiply('T', 'N', +1.0_dp, qs_ot_env%rot_mat_evec_re, mat_buf, 0.0_dp, inner_deriv_im)
3020
3021 CALL dbcsr_multiply('T', 'N', +1.0_dp, qs_ot_env%rot_mat_dedu, qs_ot_env%rot_mat_evec_re, 0.0_dp, mat_buf)
3022 CALL dbcsr_multiply('T', 'N', +1.0_dp, qs_ot_env%rot_mat_evec_re, mat_buf, 1.0_dp, inner_deriv_re)
3023 CALL dbcsr_multiply('T', 'N', -1.0_dp, qs_ot_env%rot_mat_evec_im, mat_buf, 1.0_dp, inner_deriv_im)
3024
3025 ! Real and imaginary products can have different sparse block patterns.
3026 ! Form their union explicitly and treat a missing partner block as zero.
3027 max_blocks = 0
3028 CALL dbcsr_iterator_start(iter, inner_deriv_re)
3029 DO WHILE (dbcsr_iterator_blocks_left(iter))
3030 CALL dbcsr_iterator_next_block(iter, row, col)
3031 max_blocks = max_blocks + 1
3032 END DO
3033 CALL dbcsr_iterator_stop(iter)
3034 CALL dbcsr_iterator_start(iter, inner_deriv_im)
3035 DO WHILE (dbcsr_iterator_blocks_left(iter))
3036 CALL dbcsr_iterator_next_block(iter, row, col)
3037 max_blocks = max_blocks + 1
3038 END DO
3039 CALL dbcsr_iterator_stop(iter)
3040
3041 ALLOCATE (rows(max(max_blocks, 1)), cols(max(max_blocks, 1)))
3042 nblocks = 0
3043 CALL dbcsr_iterator_start(iter, inner_deriv_re)
3044 DO WHILE (dbcsr_iterator_blocks_left(iter))
3045 CALL dbcsr_iterator_next_block(iter, row, col)
3046 duplicate = .false.
3047 DO iblock = 1, nblocks
3048 duplicate = duplicate .OR. (rows(iblock) == row .AND. cols(iblock) == col)
3049 END DO
3050 IF (.NOT. duplicate) THEN
3051 nblocks = nblocks + 1
3052 rows(nblocks) = row
3053 cols(nblocks) = col
3054 END IF
3055 END DO
3056 CALL dbcsr_iterator_stop(iter)
3057 CALL dbcsr_iterator_start(iter, inner_deriv_im)
3058 DO WHILE (dbcsr_iterator_blocks_left(iter))
3059 CALL dbcsr_iterator_next_block(iter, row, col)
3060 duplicate = .false.
3061 DO iblock = 1, nblocks
3062 duplicate = duplicate .OR. (rows(iblock) == row .AND. cols(iblock) == col)
3063 END DO
3064 IF (.NOT. duplicate) THEN
3065 nblocks = nblocks + 1
3066 rows(nblocks) = row
3067 cols(nblocks) = col
3068 END IF
3069 END DO
3070 CALL dbcsr_iterator_stop(iter)
3071
3072 CALL dbcsr_get_info(inner_deriv_re, distribution=dist, row_blk_size=row_blk_size, &
3073 col_blk_size=col_blk_size)
3074 CALL dbcsr_create(outer_deriv_re, "outer_deriv_re", dist, dbcsr_type_no_symmetry, &
3075 row_blk_size, col_blk_size)
3076 CALL dbcsr_create(outer_deriv_im, "outer_deriv_im", dist, dbcsr_type_no_symmetry, &
3077 row_blk_size, col_blk_size)
3078 IF (nblocks > 0) THEN
3079 CALL dbcsr_reserve_blocks(outer_deriv_re, rows=rows(1:nblocks), cols=cols(1:nblocks))
3080 CALL dbcsr_reserve_blocks(outer_deriv_im, rows=rows(1:nblocks), cols=cols(1:nblocks))
3081 END IF
3082 CALL dbcsr_finalize(outer_deriv_re)
3083 CALL dbcsr_finalize(outer_deriv_im)
3084 CALL dbcsr_set(outer_deriv_re, 0.0_dp)
3085 CALL dbcsr_set(outer_deriv_im, 0.0_dp)
3086
3087 CALL dbcsr_get_info(outer_deriv_re, row_blk_offset=row_blk_offset, col_blk_offset=col_blk_offset)
3088 CALL dbcsr_iterator_start(iter, outer_deriv_re)
3089 DO WHILE (dbcsr_iterator_blocks_left(iter))
3090 CALL dbcsr_iterator_next_block(iter, row, col)
3091 CALL dbcsr_get_block_p(inner_deriv_re, row, col, block_in_re, found_in_re)
3092 CALL dbcsr_get_block_p(inner_deriv_im, row, col, block_in_im, found_in_im)
3093 CALL dbcsr_get_block_p(outer_deriv_re, row, col, block_out_re, found_out_re)
3094 CALL dbcsr_get_block_p(outer_deriv_im, row, col, block_out_im, found_out_im)
3095 cpassert(found_out_re .AND. found_out_im)
3096
3097 DO i = 1, SIZE(block_out_re, 1)
3098 DO j = 1, SIZE(block_out_re, 2)
3099 e1 = qs_ot_env%rot_mat_evals(row_blk_offset(row) + i - 1)
3100 e2 = qs_ot_env%rot_mat_evals(col_blk_offset(col) + j - 1)
3101 re_part = 0.0_dp
3102 im_part = 0.0_dp
3103 IF (found_in_re) re_part = block_in_re(i, j)
3104 IF (found_in_im) im_part = block_in_im(i, j)
3105 cval_in = cmplx(re_part, im_part, dp)
3106 cval_out = cval_in*cint(e1, e2)
3107 block_out_re(i, j) = real(cval_out)
3108 block_out_im(i, j) = aimag(cval_out)
3109 END DO
3110 END DO
3111 END DO
3112 CALL dbcsr_iterator_stop(iter)
3113 DEALLOCATE (rows, cols)
3114 CALL dbcsr_release(inner_deriv_re)
3115 CALL dbcsr_release(inner_deriv_im)
3116
3117 ! Compute: matrix_buf1 = \sum_i outer_deriv_ij * |rot_mat_evec_i> <rot_mat_evec_j|
3118 CALL dbcsr_multiply('N', 'N', +1.0_dp, qs_ot_env%rot_mat_evec_re, outer_deriv_re, 0.0_dp, mat_buf)
3119 CALL dbcsr_multiply('N', 'N', -1.0_dp, qs_ot_env%rot_mat_evec_im, outer_deriv_im, 1.0_dp, mat_buf)
3120 CALL dbcsr_multiply('N', 'T', +1.0_dp, mat_buf, qs_ot_env%rot_mat_evec_re, 0.0_dp, qs_ot_env%matrix_buf1)
3121
3122 CALL dbcsr_multiply('N', 'N', +1.0_dp, qs_ot_env%rot_mat_evec_re, outer_deriv_im, 0.0_dp, mat_buf)
3123 CALL dbcsr_multiply('N', 'N', +1.0_dp, qs_ot_env%rot_mat_evec_im, outer_deriv_re, 1.0_dp, mat_buf)
3124 CALL dbcsr_multiply('N', 'T', +1.0_dp, mat_buf, qs_ot_env%rot_mat_evec_im, 1.0_dp, qs_ot_env%matrix_buf1)
3125
3126 ! Account for anti-symmetry of rot_mat_x without relying on
3127 ! STRICT-only matrix buffers. REF uses the same finite chart.
3128 CALL qs_ot_square_transpose(qs_ot_env%matrix_buf1, mat_buf, qs_ot_env%rot_mat_u)
3129 CALL dbcsr_copy(qs_ot_env%rot_mat_gx, mat_buf)
3130 CALL dbcsr_add(qs_ot_env%rot_mat_gx, qs_ot_env%matrix_buf1, &
3131 alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
3132
3133 CALL dbcsr_release(mat_buf)
3134 CALL dbcsr_release(outer_deriv_re)
3135 CALL dbcsr_release(outer_deriv_im)
3136 END IF
3137 CALL timestop(handle)
3138 CONTAINS
3139
3140! **************************************************************************************************
3141!> \brief ...
3142!> \param e1 ...
3143!> \param e2 ...
3144!> \return ...
3145! **************************************************************************************************
3146 FUNCTION cint(e1, e2)
3147 REAL(kind=dp) :: e1, e2
3148 COMPLEX(KIND=dp) :: cint
3149
3150 COMPLEX(KIND=dp) :: l1, l2, x
3151 INTEGER :: i
3152
3153 l1 = (0.0_dp, -1.0_dp)*e1
3154 l2 = (0.0_dp, -1.0_dp)*e2
3155 IF (abs(l1 - l2) > 0.5_dp) THEN
3156 cint = (exp(l1) - exp(l2))/(l1 - l2)
3157 ELSE
3158 x = 1.0_dp
3159 cint = 0.0_dp
3160 DO i = 1, 16
3161 cint = cint + x
3162 x = x*(l1 - l2)/real(i + 1, kind=dp)
3163 END DO
3164 cint = cint*exp(l2)
3165 END IF
3166 END FUNCTION cint
3167 END SUBROUTINE qs_ot_rot_mat_derivative
3168
3169! **************************************************************************************************
3170!> \brief decide strategy
3171!> tries to decide if the taylor expansion of cos(sqrt(xsx)) converges rapidly enough
3172!> to make a taylor expansion of the functions cos(sqrt(xsx)) and sin(sqrt(xsx))/sqrt(xsx)
3173!> and their derivatives faster than their computation based on diagonalization since xsx can
3174!> be very small, especially during dynamics, only a few terms might indeed be needed we find
3175!> the necessary order N to have largest_eval_upper_bound**(N+1)/(2(N+1))! < eps_taylor
3176!> \param qs_ot_env ...
3177! **************************************************************************************************
3178 SUBROUTINE decide_strategy(qs_ot_env)
3179 TYPE(qs_ot_type) :: qs_ot_env
3180
3181 INTEGER :: n
3182 REAL(kind=dp) :: num_error
3183
3184 qs_ot_env%do_taylor = .false.
3185 n = 0
3186 num_error = qs_ot_env%largest_eval_upper_bound/(2.0_dp)
3187 DO WHILE (num_error > qs_ot_env%settings%eps_taylor .AND. n <= qs_ot_env%settings%max_taylor)
3188 n = n + 1
3189 num_error = num_error*qs_ot_env%largest_eval_upper_bound/real((2*n + 1)*(2*n + 2), kind=dp)
3190 END DO
3191 qs_ot_env%taylor_order = n
3192 IF (qs_ot_env%taylor_order <= qs_ot_env%settings%max_taylor) THEN
3193 qs_ot_env%do_taylor = .true.
3194 END IF
3195
3196 END SUBROUTINE decide_strategy
3197
3198! **************************************************************************************************
3199!> \brief c=(c0*cos(p^0.5)+x*sin(p^0.5)*p^(-0.5)) x rot_mat_u
3200!> this assumes that x is already ortho to S*C0, and that p is x*S*x
3201!> rot_mat_u is an optional rotation matrix
3202!> \param matrix_c ...
3203!> \param matrix_x ...
3204!> \param qs_ot_env ...
3205! **************************************************************************************************
3206 SUBROUTINE qs_ot_get_orbitals(matrix_c, matrix_x, qs_ot_env)
3207
3208 TYPE(dbcsr_type), POINTER :: matrix_c, matrix_x
3209 TYPE(qs_ot_type) :: qs_ot_env
3210
3211 CHARACTER(len=*), PARAMETER :: routinen = 'qs_ot_get_orbitals'
3212 REAL(kind=dp), PARAMETER :: rone = 1.0_dp, rzero = 0.0_dp
3213
3214 INTEGER :: handle, k, n
3215 TYPE(dbcsr_type), POINTER :: matrix_kk
3216
3217 CALL timeset(routinen, handle)
3218
3219 CALL dbcsr_get_info(matrix_x, nfullrows_total=n, nfullcols_total=k)
3220
3221 ! rotate the multiplying matrices cosp and sinp instead of the result,
3222 ! this should be cheaper for large basis sets
3223 IF (qs_ot_env%settings%do_rotation) THEN
3224 matrix_kk => qs_ot_env%matrix_buf1
3225 CALL dbcsr_multiply('N', 'N', rone, qs_ot_env%matrix_cosp, &
3226 qs_ot_env%rot_mat_u, rzero, matrix_kk)
3227 ELSE
3228 matrix_kk => qs_ot_env%matrix_cosp
3229 END IF
3230
3231 CALL dbcsr_multiply('N', 'N', rone, qs_ot_env%matrix_c0, matrix_kk, &
3232 rzero, matrix_c)
3233
3234 IF (qs_ot_env%settings%do_rotation) THEN
3235 matrix_kk => qs_ot_env%matrix_buf1
3236 CALL dbcsr_multiply('N', 'N', rone, qs_ot_env%matrix_sinp, &
3237 qs_ot_env%rot_mat_u, rzero, matrix_kk)
3238 ELSE
3239 matrix_kk => qs_ot_env%matrix_sinp
3240 END IF
3241 CALL dbcsr_multiply('N', 'N', rone, matrix_x, matrix_kk, &
3242 rone, matrix_c)
3243
3244 CALL timestop(handle)
3245
3246 END SUBROUTINE qs_ot_get_orbitals
3247
3248! **************************************************************************************************
3249!> \brief update complex K-point orbitals with the finite STRICT transformation
3250!> \param matrix_c real output orbitals
3251!> \param matrix_c_im imaginary output orbitals
3252!> \param matrix_s real overlap matrix
3253!> \param matrix_s_im imaginary overlap matrix
3254!> \param qs_ot_env OT channel state
3255! **************************************************************************************************
3256 SUBROUTINE qs_ot_get_orbitals_complex(matrix_c, matrix_c_im, matrix_s, matrix_s_im, qs_ot_env)
3257 TYPE(dbcsr_type), POINTER :: matrix_c, matrix_c_im, matrix_s, &
3258 matrix_s_im
3259 TYPE(qs_ot_type) :: qs_ot_env
3260
3261 CHARACTER(len=*), PARAMETER :: routinen = 'qs_ot_get_orbitals_complex'
3262
3263 INTEGER :: handle
3264 TYPE(dbcsr_type) :: rotated_im, rotated_re, rotation_tmp
3265
3266 CALL timeset(routinen, handle)
3267
3268 cpassert(qs_ot_env%has_complex_kpoint_state)
3269 cpassert(qs_ot_env%settings%ot_algorithm == 'TOD')
3270
3271 CALL qs_ot_complex_multiply('N', 'N', matrix_s, matrix_s_im, &
3272 qs_ot_env%matrix_x, qs_ot_env%matrix_x_im, &
3273 qs_ot_env%matrix_sx, qs_ot_env%matrix_sx_im, &
3274 qs_ot_env%matrix_tmp_nk)
3275 CALL qs_ot_get_p_complex(qs_ot_env%matrix_x, qs_ot_env%matrix_x_im, &
3276 qs_ot_env%matrix_sx, qs_ot_env%matrix_sx_im, qs_ot_env)
3277
3278 CALL qs_ot_complex_multiply('N', 'N', qs_ot_env%matrix_c0, qs_ot_env%matrix_c0_im, &
3279 qs_ot_env%matrix_cosp, qs_ot_env%matrix_cosp_im, &
3280 matrix_c, matrix_c_im, qs_ot_env%matrix_tmp_nk)
3281 CALL qs_ot_complex_multiply('N', 'N', qs_ot_env%matrix_x, qs_ot_env%matrix_x_im, &
3282 qs_ot_env%matrix_sinp, qs_ot_env%matrix_sinp_im, &
3283 qs_ot_env%matrix_buf_nk, qs_ot_env%matrix_buf_nk_im, &
3284 qs_ot_env%matrix_tmp_nk)
3285 CALL dbcsr_add(matrix_c, qs_ot_env%matrix_buf_nk, alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
3286 CALL dbcsr_add(matrix_c_im, qs_ot_env%matrix_buf_nk_im, alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
3287
3288 IF (qs_ot_env%settings%do_rotation) THEN
3289 CALL qs_ot_generate_rotation_complex(qs_ot_env)
3290
3291 CALL dbcsr_copy(rotated_re, matrix_c, name="strict_rotated_re")
3292 CALL dbcsr_copy(rotated_im, matrix_c_im, name="strict_rotated_im")
3293 CALL dbcsr_copy(rotation_tmp, matrix_c, name="strict_rotation_tmp")
3294
3295 ! C_out = Q(X)*U for the finite STRICT chart Q(X).
3296 CALL dbcsr_multiply('N', 'N', 1.0_dp, matrix_c, qs_ot_env%rot_mat_u, &
3297 0.0_dp, rotated_re)
3298 CALL dbcsr_multiply('N', 'N', 1.0_dp, matrix_c_im, qs_ot_env%rot_mat_u_im, &
3299 0.0_dp, rotation_tmp)
3300 CALL dbcsr_add(rotated_re, rotation_tmp, alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
3301
3302 CALL dbcsr_multiply('N', 'N', 1.0_dp, matrix_c, qs_ot_env%rot_mat_u_im, &
3303 0.0_dp, rotated_im)
3304 CALL dbcsr_multiply('N', 'N', 1.0_dp, matrix_c_im, qs_ot_env%rot_mat_u, &
3305 0.0_dp, rotation_tmp)
3306 CALL dbcsr_add(rotated_im, rotation_tmp, alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
3307
3308 CALL dbcsr_copy(matrix_c, rotated_re)
3309 CALL dbcsr_copy(matrix_c_im, rotated_im)
3310 CALL dbcsr_release(rotated_re)
3311 CALL dbcsr_release(rotated_im)
3312 CALL dbcsr_release(rotation_tmp)
3313 END IF
3314
3315 CALL timestop(handle)
3316
3317 END SUBROUTINE qs_ot_get_orbitals_complex
3318
3319! **************************************************************************************************
3320!> \brief this routines computes dE/dx=dx, with dx ortho to sc0
3321!> needs dE/dC=hc,C0,X,SX,p
3322!> if preconditioned it will not be the derivative, but the lagrangian multiplier
3323!> is changed so that P*dE/dx is the right derivative (i.e. in the allowed subspace)
3324!> \param matrix_hc ...
3325!> \param matrix_x ...
3326!> \param matrix_sx ...
3327!> \param matrix_gx ...
3328!> \param qs_ot_env ...
3329! **************************************************************************************************
3330 SUBROUTINE qs_ot_get_derivative(matrix_hc, matrix_x, matrix_sx, matrix_gx, &
3331 qs_ot_env)
3332 TYPE(dbcsr_type), POINTER :: matrix_hc, matrix_x, matrix_sx, matrix_gx
3333 TYPE(qs_ot_type) :: qs_ot_env
3334
3335 CHARACTER(len=*), PARAMETER :: routinen = 'qs_ot_get_derivative'
3336 REAL(kind=dp), PARAMETER :: rone = 1.0_dp, rzero = 0.0_dp
3337
3338 INTEGER :: handle, k, n, ortho_k
3339 TYPE(dbcsr_type), POINTER :: matrix_hc_local, matrix_target
3340
3341 CALL timeset(routinen, handle)
3342
3343 NULLIFY (matrix_hc_local)
3344
3345 CALL dbcsr_get_info(matrix_x, nfullrows_total=n, nfullcols_total=k)
3346
3347 ! could in principle be taken inside qs_ot_get_derivative_* for increased efficiency
3348 ! create a local rotated version of matrix_hc leaving matrix_hc untouched (needed
3349 ! for lagrangian multipliers)
3350 IF (qs_ot_env%settings%do_rotation) THEN
3351 CALL dbcsr_copy(matrix_gx, matrix_hc) ! use gx as temporary
3352 CALL dbcsr_init_p(matrix_hc_local)
3353 CALL dbcsr_copy(matrix_hc_local, matrix_hc, name='matrix_hc_local')
3354 CALL dbcsr_set(matrix_hc_local, 0.0_dp)
3355 CALL dbcsr_multiply('N', 'T', rone, matrix_gx, qs_ot_env%rot_mat_u, rzero, matrix_hc_local)
3356 ELSE
3357 matrix_hc_local => matrix_hc
3358 END IF
3359
3360 IF (qs_ot_env%do_taylor) THEN
3361 CALL qs_ot_get_derivative_taylor(matrix_hc_local, matrix_x, matrix_sx, matrix_gx, qs_ot_env)
3362 ELSE
3363 CALL qs_ot_get_derivative_diag(matrix_hc_local, matrix_x, matrix_sx, matrix_gx, qs_ot_env)
3364 END IF
3365
3366 ! and make it orthogonal
3367 CALL dbcsr_get_info(qs_ot_env%matrix_sc0, nfullcols_total=ortho_k)
3368
3369 IF (ASSOCIATED(qs_ot_env%preconditioner)) THEN
3370 matrix_target => qs_ot_env%matrix_psc0
3371 ELSE
3372 matrix_target => qs_ot_env%matrix_sc0
3373 END IF
3374 ! first make the matrix os if not yet valid
3375 IF (.NOT. qs_ot_env%os_valid) THEN
3376 ! this assumes that the preconditioner is a single matrix
3377 ! that maps sc0 onto psc0
3378
3379 IF (ASSOCIATED(qs_ot_env%preconditioner)) THEN
3380 CALL apply_preconditioner(qs_ot_env%preconditioner, qs_ot_env%matrix_sc0, &
3381 qs_ot_env%matrix_psc0)
3382 END IF
3383 CALL dbcsr_multiply('T', 'N', rone, &
3384 qs_ot_env%matrix_sc0, matrix_target, &
3385 rzero, qs_ot_env%matrix_os)
3386 CALL cp_dbcsr_cholesky_decompose(qs_ot_env%matrix_os, &
3387 para_env=qs_ot_env%para_env, blacs_env=qs_ot_env%blacs_env)
3388 CALL cp_dbcsr_cholesky_invert(qs_ot_env%matrix_os, &
3389 para_env=qs_ot_env%para_env, blacs_env=qs_ot_env%blacs_env, &
3390 uplo_to_full=.true.)
3391 qs_ot_env%os_valid = .true.
3392 END IF
3393 CALL dbcsr_multiply('T', 'N', rone, matrix_target, matrix_gx, &
3394 rzero, qs_ot_env%matrix_buf1_ortho)
3395 CALL dbcsr_multiply('N', 'N', rone, qs_ot_env%matrix_os, &
3396 qs_ot_env%matrix_buf1_ortho, rzero, qs_ot_env%matrix_buf2_ortho)
3397 CALL dbcsr_multiply('N', 'N', -rone, qs_ot_env%matrix_sc0, &
3398 qs_ot_env%matrix_buf2_ortho, rone, matrix_gx)
3399 ! also treat the rot_mat gradient here
3400 IF (qs_ot_env%settings%do_rotation) THEN
3401 CALL qs_ot_rot_mat_derivative(qs_ot_env)
3402 END IF
3403
3404 IF (qs_ot_env%settings%do_rotation) THEN
3405 CALL dbcsr_release_p(matrix_hc_local)
3406 END IF
3407
3408 CALL timestop(handle)
3409
3410 END SUBROUTINE qs_ot_get_derivative
3411
3412! **************************************************************************************************
3413!> \brief Prepare the inverse metric used to project a complex STRICT gradient.
3414!> An unusable preconditioner is detached before any minimizer history is updated.
3415!> \param qs_ot_env OT channel state
3416!> \param preconditioner_rejected true if the attached preconditioner was not positive definite
3417! **************************************************************************************************
3418 SUBROUTINE qs_ot_prepare_complex_tangent_metric(qs_ot_env, preconditioner_rejected)
3419 TYPE(qs_ot_type) :: qs_ot_env
3420 LOGICAL, INTENT(OUT), OPTIONAL :: preconditioner_rejected
3421
3422 INTEGER :: i, k
3423 REAL(kind=dp) :: eval_scale, eval_threshold
3424 TYPE(dbcsr_type), POINTER :: target_im, target_re
3425
3426 IF (PRESENT(preconditioner_rejected)) preconditioner_rejected = .false.
3427 IF (qs_ot_env%os_valid) RETURN
3428
3429 IF (ASSOCIATED(qs_ot_env%preconditioner)) THEN
3430 target_re => qs_ot_env%matrix_psc0
3431 target_im => qs_ot_env%matrix_psc0_im
3432 CALL apply_preconditioner(qs_ot_env%preconditioner, &
3433 qs_ot_env%matrix_sc0, qs_ot_env%matrix_sc0_im, &
3434 target_re, target_im)
3435 ELSE
3436 target_re => qs_ot_env%matrix_sc0
3437 target_im => qs_ot_env%matrix_sc0_im
3438 END IF
3439 CALL qs_ot_complex_multiply('C', 'N', qs_ot_env%matrix_sc0, qs_ot_env%matrix_sc0_im, &
3440 target_re, target_im, qs_ot_env%matrix_os, qs_ot_env%matrix_os_im, &
3441 qs_ot_env%matrix_buf1)
3442 CALL dbcsr_get_info(qs_ot_env%matrix_os, nfullrows_total=k)
3443 CALL cp_dbcsr_heevd(matrix_re=qs_ot_env%matrix_os, matrix_im=qs_ot_env%matrix_os_im, &
3444 eigenvectors_re=qs_ot_env%matrix_r, eigenvectors_im=qs_ot_env%matrix_r_im, &
3445 eigenvalues=qs_ot_env%evals, para_env=qs_ot_env%para_env, &
3446 blacs_env=qs_ot_env%blacs_env)
3447 eval_scale = max(1.0_dp, maxval(abs(qs_ot_env%evals(1:k))))
3448 eval_threshold = 100.0_dp*epsilon(1.0_dp)*eval_scale
3449 IF (minval(qs_ot_env%evals(1:k)) <= eval_threshold .AND. &
3450 ASSOCIATED(qs_ot_env%preconditioner)) THEN
3451 NULLIFY (qs_ot_env%preconditioner)
3452 IF (PRESENT(preconditioner_rejected)) preconditioner_rejected = .true.
3453 target_re => qs_ot_env%matrix_sc0
3454 target_im => qs_ot_env%matrix_sc0_im
3455 CALL qs_ot_complex_multiply('C', 'N', qs_ot_env%matrix_sc0, qs_ot_env%matrix_sc0_im, &
3456 target_re, target_im, qs_ot_env%matrix_os, qs_ot_env%matrix_os_im, &
3457 qs_ot_env%matrix_buf1)
3458 CALL cp_dbcsr_heevd(matrix_re=qs_ot_env%matrix_os, matrix_im=qs_ot_env%matrix_os_im, &
3459 eigenvectors_re=qs_ot_env%matrix_r, eigenvectors_im=qs_ot_env%matrix_r_im, &
3460 eigenvalues=qs_ot_env%evals, para_env=qs_ot_env%para_env, &
3461 blacs_env=qs_ot_env%blacs_env)
3462 eval_scale = max(1.0_dp, maxval(abs(qs_ot_env%evals(1:k))))
3463 eval_threshold = 100.0_dp*epsilon(1.0_dp)*eval_scale
3464 END IF
3465 cpassert(minval(qs_ot_env%evals(1:k)) > eval_threshold)
3466 DO i = 1, k
3467 qs_ot_env%dum(i) = 1.0_dp/qs_ot_env%evals(i)
3468 END DO
3469 CALL dbcsr_copy(qs_ot_env%matrix_buf1, qs_ot_env%matrix_r)
3470 CALL dbcsr_copy(qs_ot_env%matrix_buf1_im, qs_ot_env%matrix_r_im)
3471 CALL dbcsr_scale_by_vector(qs_ot_env%matrix_buf1, alpha=qs_ot_env%dum, side='right')
3472 CALL dbcsr_scale_by_vector(qs_ot_env%matrix_buf1_im, alpha=qs_ot_env%dum, side='right')
3473 CALL qs_ot_complex_multiply('N', 'C', qs_ot_env%matrix_buf1, qs_ot_env%matrix_buf1_im, &
3474 qs_ot_env%matrix_r, qs_ot_env%matrix_r_im, &
3475 qs_ot_env%matrix_os, qs_ot_env%matrix_os_im, &
3476 qs_ot_env%matrix_buf2)
3477 qs_ot_env%os_valid = .true.
3478
3479 END SUBROUTINE qs_ot_prepare_complex_tangent_metric
3480
3481! **************************************************************************************************
3482!> \brief finite complex STRICT derivative, projected onto C0^H*S*X=0
3483!> \param matrix_hc real part of H(k)*C(k)
3484!> \param matrix_hc_im imaginary part of H(k)*C(k)
3485!> \param qs_ot_env OT channel state
3486!> \param matrix_hc_rotation occupation-weighted H(k)C(k) for the rotation channel
3487!> \param matrix_hc_rotation_im imaginary component of matrix_hc_rotation
3488! **************************************************************************************************
3489 SUBROUTINE qs_ot_get_derivative_complex(matrix_hc, matrix_hc_im, qs_ot_env, &
3490 matrix_hc_rotation, matrix_hc_rotation_im)
3491 TYPE(dbcsr_type), POINTER :: matrix_hc, matrix_hc_im
3492 TYPE(qs_ot_type) :: qs_ot_env
3493 TYPE(dbcsr_type), OPTIONAL, POINTER :: matrix_hc_rotation, matrix_hc_rotation_im
3494
3495 CHARACTER(len=*), PARAMETER :: routinen = 'qs_ot_get_derivative_complex'
3496
3497 INTEGER :: handle
3498 TYPE(dbcsr_distribution_type) :: dist
3499 TYPE(dbcsr_type) :: tmp_nk
3500 TYPE(dbcsr_type), POINTER :: hc_rotation_im, hc_rotation_re, &
3501 hc_work_im, hc_work_re, target_im, &
3502 target_re
3503 TYPE(dbcsr_type), TARGET :: hc_rot_im, hc_rot_re
3504
3505 CALL timeset(routinen, handle)
3506
3507 cpassert(qs_ot_env%has_complex_kpoint_state)
3508 cpassert(qs_ot_env%settings%ot_algorithm == 'TOD')
3509
3510 hc_rotation_re => matrix_hc
3511 hc_rotation_im => matrix_hc_im
3512 IF (PRESENT(matrix_hc_rotation) .OR. PRESENT(matrix_hc_rotation_im)) THEN
3513 cpassert(PRESENT(matrix_hc_rotation) .AND. PRESENT(matrix_hc_rotation_im))
3514 cpassert(ASSOCIATED(matrix_hc_rotation))
3515 cpassert(ASSOCIATED(matrix_hc_rotation_im))
3516 hc_rotation_re => matrix_hc_rotation
3517 hc_rotation_im => matrix_hc_rotation_im
3518 END IF
3519 hc_work_re => matrix_hc
3520 hc_work_im => matrix_hc_im
3521
3522 IF (qs_ot_env%settings%do_rotation) THEN
3523 CALL qs_ot_generate_rotation_complex(qs_ot_env)
3524
3525 ! Reconstruct the unrotated finite STRICT orbitals Q(X).
3526 CALL qs_ot_complex_multiply('N', 'N', qs_ot_env%matrix_c0, qs_ot_env%matrix_c0_im, &
3527 qs_ot_env%matrix_cosp, qs_ot_env%matrix_cosp_im, &
3528 qs_ot_env%matrix_buf_nk, qs_ot_env%matrix_buf_nk_im, &
3529 qs_ot_env%matrix_tmp_nk)
3530 CALL qs_ot_complex_multiply('N', 'N', qs_ot_env%matrix_x, qs_ot_env%matrix_x_im, &
3531 qs_ot_env%matrix_sinp, qs_ot_env%matrix_sinp_im, &
3532 qs_ot_env%matrix_gx, qs_ot_env%matrix_gx_im, &
3533 qs_ot_env%matrix_tmp_nk)
3534 CALL dbcsr_add(qs_ot_env%matrix_buf_nk, qs_ot_env%matrix_gx, &
3535 alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
3536 CALL dbcsr_add(qs_ot_env%matrix_buf_nk_im, qs_ot_env%matrix_gx_im, &
3537 alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
3538
3539 ! dF/dU = Q(X)^H*G_C for C=Q(X)*U.
3540 CALL qs_ot_complex_multiply('C', 'N', &
3541 qs_ot_env%matrix_buf_nk, qs_ot_env%matrix_buf_nk_im, &
3542 hc_rotation_re, hc_rotation_im, &
3543 qs_ot_env%rot_mat_dedu, qs_ot_env%rot_mat_dedu_im, &
3544 qs_ot_env%matrix_buf1)
3545 CALL qs_ot_rot_mat_derivative_complex(qs_ot_env)
3546
3547 ! The STRICT coordinate sees G_Q=G_C*U^H.
3548 CALL dbcsr_copy(hc_rot_re, matrix_hc, name="strict_hc_rot_re")
3549 CALL dbcsr_copy(hc_rot_im, matrix_hc_im, name="strict_hc_rot_im")
3550 CALL dbcsr_copy(tmp_nk, matrix_hc, name="strict_hc_rot_tmp")
3551 CALL dbcsr_multiply('N', 'T', 1.0_dp, matrix_hc, qs_ot_env%rot_mat_u, &
3552 0.0_dp, hc_rot_re)
3553 CALL dbcsr_multiply('N', 'T', 1.0_dp, matrix_hc_im, qs_ot_env%rot_mat_u_im, &
3554 0.0_dp, tmp_nk)
3555 CALL dbcsr_add(hc_rot_re, tmp_nk, alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
3556
3557 CALL dbcsr_multiply('N', 'T', 1.0_dp, matrix_hc_im, qs_ot_env%rot_mat_u, &
3558 0.0_dp, hc_rot_im)
3559 CALL dbcsr_multiply('N', 'T', 1.0_dp, matrix_hc, qs_ot_env%rot_mat_u_im, &
3560 0.0_dp, tmp_nk)
3561 CALL dbcsr_add(hc_rot_im, tmp_nk, alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
3562 hc_work_re => hc_rot_re
3563 hc_work_im => hc_rot_im
3564 END IF
3565
3566 ! Direct X contribution, H*C sinc(sqrt(P)).
3567 CALL qs_ot_complex_multiply('N', 'N', hc_work_re, hc_work_im, &
3568 qs_ot_env%matrix_sinp, qs_ot_env%matrix_sinp_im, &
3569 qs_ot_env%matrix_gx, qs_ot_env%matrix_gx_im, &
3570 qs_ot_env%matrix_tmp_nk)
3571
3572 ! Frechet contribution from X sinc(sqrt(P)).
3573 CALL qs_ot_complex_multiply('C', 'N', hc_work_re, hc_work_im, &
3574 qs_ot_env%matrix_x, qs_ot_env%matrix_x_im, &
3575 qs_ot_env%matrix_buf2, qs_ot_env%matrix_buf2_im, &
3576 qs_ot_env%matrix_buf1)
3577 CALL qs_ot_complex_multiply('N', 'N', qs_ot_env%matrix_buf2, qs_ot_env%matrix_buf2_im, &
3578 qs_ot_env%matrix_r, qs_ot_env%matrix_r_im, &
3579 qs_ot_env%matrix_buf1, qs_ot_env%matrix_buf1_im, &
3580 qs_ot_env%matrix_buf4)
3581 CALL qs_ot_complex_multiply('C', 'N', qs_ot_env%matrix_r, qs_ot_env%matrix_r_im, &
3582 qs_ot_env%matrix_buf1, qs_ot_env%matrix_buf1_im, &
3583 qs_ot_env%matrix_buf2, qs_ot_env%matrix_buf2_im, &
3584 qs_ot_env%matrix_buf4)
3585 CALL dbcsr_hadamard_product(qs_ot_env%matrix_buf2, qs_ot_env%matrix_sinp_b, &
3586 qs_ot_env%matrix_buf3)
3587 CALL dbcsr_hadamard_product(qs_ot_env%matrix_buf2_im, qs_ot_env%matrix_sinp_b, &
3588 qs_ot_env%matrix_buf3_im)
3589
3590 ! Frechet contribution from C0 cos(sqrt(P)).
3591 CALL qs_ot_complex_multiply('C', 'N', hc_work_re, hc_work_im, &
3592 qs_ot_env%matrix_c0, qs_ot_env%matrix_c0_im, &
3593 qs_ot_env%matrix_buf2, qs_ot_env%matrix_buf2_im, &
3594 qs_ot_env%matrix_buf1)
3595 CALL qs_ot_complex_multiply('N', 'N', qs_ot_env%matrix_buf2, qs_ot_env%matrix_buf2_im, &
3596 qs_ot_env%matrix_r, qs_ot_env%matrix_r_im, &
3597 qs_ot_env%matrix_buf1, qs_ot_env%matrix_buf1_im, &
3598 qs_ot_env%matrix_buf4)
3599 CALL qs_ot_complex_multiply('C', 'N', qs_ot_env%matrix_r, qs_ot_env%matrix_r_im, &
3600 qs_ot_env%matrix_buf1, qs_ot_env%matrix_buf1_im, &
3601 qs_ot_env%matrix_buf2, qs_ot_env%matrix_buf2_im, &
3602 qs_ot_env%matrix_buf4)
3603 CALL dbcsr_hadamard_product(qs_ot_env%matrix_buf2, qs_ot_env%matrix_cosp_b, &
3604 qs_ot_env%matrix_buf4)
3605 CALL dbcsr_hadamard_product(qs_ot_env%matrix_buf2_im, qs_ot_env%matrix_cosp_b, &
3606 qs_ot_env%matrix_buf4_im)
3607 CALL dbcsr_add(qs_ot_env%matrix_buf3, qs_ot_env%matrix_buf4, &
3608 alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
3609 CALL dbcsr_add(qs_ot_env%matrix_buf3_im, qs_ot_env%matrix_buf4_im, &
3610 alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
3611
3612 ! Transform back and add the Hermitian adjoint generated by dP.
3613 CALL qs_ot_complex_multiply('N', 'N', qs_ot_env%matrix_r, qs_ot_env%matrix_r_im, &
3614 qs_ot_env%matrix_buf3, qs_ot_env%matrix_buf3_im, &
3615 qs_ot_env%matrix_buf1, qs_ot_env%matrix_buf1_im, &
3616 qs_ot_env%matrix_buf2)
3617 CALL qs_ot_complex_multiply('N', 'C', qs_ot_env%matrix_buf1, qs_ot_env%matrix_buf1_im, &
3618 qs_ot_env%matrix_r, qs_ot_env%matrix_r_im, &
3619 qs_ot_env%matrix_buf3, qs_ot_env%matrix_buf3_im, &
3620 qs_ot_env%matrix_buf2)
3621 CALL dbcsr_get_info(qs_ot_env%matrix_buf3, distribution=dist)
3622 CALL dbcsr_transposed(qs_ot_env%matrix_buf4, qs_ot_env%matrix_buf3, &
3623 shallow_data_copy=.false., use_distribution=dist, &
3624 transpose_distribution=.false.)
3625 CALL dbcsr_transposed(qs_ot_env%matrix_buf4_im, qs_ot_env%matrix_buf3_im, &
3626 shallow_data_copy=.false., use_distribution=dist, &
3627 transpose_distribution=.false.)
3628 CALL dbcsr_add(qs_ot_env%matrix_buf3, qs_ot_env%matrix_buf4, &
3629 alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
3630 CALL dbcsr_add(qs_ot_env%matrix_buf3_im, qs_ot_env%matrix_buf4_im, &
3631 alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
3632
3633 CALL qs_ot_complex_multiply('N', 'N', qs_ot_env%matrix_sx, qs_ot_env%matrix_sx_im, &
3634 qs_ot_env%matrix_buf3, qs_ot_env%matrix_buf3_im, &
3635 qs_ot_env%matrix_buf_nk, qs_ot_env%matrix_buf_nk_im, &
3636 qs_ot_env%matrix_tmp_nk)
3637 CALL dbcsr_add(qs_ot_env%matrix_gx, qs_ot_env%matrix_buf_nk, &
3638 alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
3639 CALL dbcsr_add(qs_ot_env%matrix_gx_im, qs_ot_env%matrix_buf_nk_im, &
3640 alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
3641
3642 ! Preconditioner-aware projection onto the complex STRICT tangent space.
3643 CALL qs_ot_prepare_complex_tangent_metric(qs_ot_env)
3644 IF (ASSOCIATED(qs_ot_env%preconditioner)) THEN
3645 target_re => qs_ot_env%matrix_psc0
3646 target_im => qs_ot_env%matrix_psc0_im
3647 ELSE
3648 target_re => qs_ot_env%matrix_sc0
3649 target_im => qs_ot_env%matrix_sc0_im
3650 END IF
3651
3652 CALL qs_ot_complex_multiply('C', 'N', target_re, target_im, &
3653 qs_ot_env%matrix_gx, qs_ot_env%matrix_gx_im, &
3654 qs_ot_env%matrix_buf1_ortho, qs_ot_env%matrix_buf1_ortho_im, &
3655 qs_ot_env%matrix_tmp_ortho)
3656 CALL qs_ot_complex_multiply('N', 'N', qs_ot_env%matrix_os, qs_ot_env%matrix_os_im, &
3657 qs_ot_env%matrix_buf1_ortho, qs_ot_env%matrix_buf1_ortho_im, &
3658 qs_ot_env%matrix_buf2_ortho, qs_ot_env%matrix_buf2_ortho_im, &
3659 qs_ot_env%matrix_tmp_ortho)
3660 CALL qs_ot_complex_multiply('N', 'N', qs_ot_env%matrix_sc0, qs_ot_env%matrix_sc0_im, &
3661 qs_ot_env%matrix_buf2_ortho, qs_ot_env%matrix_buf2_ortho_im, &
3662 qs_ot_env%matrix_buf_nk, qs_ot_env%matrix_buf_nk_im, &
3663 qs_ot_env%matrix_tmp_nk)
3664 CALL dbcsr_add(qs_ot_env%matrix_gx, qs_ot_env%matrix_buf_nk, &
3665 alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
3666 CALL dbcsr_add(qs_ot_env%matrix_gx_im, qs_ot_env%matrix_buf_nk_im, &
3667 alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
3668
3669 IF (qs_ot_env%settings%do_rotation) THEN
3670 CALL dbcsr_release(hc_rot_re)
3671 CALL dbcsr_release(hc_rot_im)
3672 CALL dbcsr_release(tmp_nk)
3673 END IF
3674
3675 CALL timestop(handle)
3676
3677 END SUBROUTINE qs_ot_get_derivative_complex
3678
3679! **************************************************************************************************
3680!> \brief ...
3681!> \param matrix_hc ...
3682!> \param matrix_x ...
3683!> \param matrix_sx ...
3684!> \param matrix_gx ...
3685!> \param qs_ot_env ...
3686! **************************************************************************************************
3687 SUBROUTINE qs_ot_get_derivative_diag(matrix_hc, matrix_x, matrix_sx, matrix_gx, &
3688 qs_ot_env)
3689
3690 TYPE(dbcsr_type), POINTER :: matrix_hc, matrix_x, matrix_sx, matrix_gx
3691 TYPE(qs_ot_type) :: qs_ot_env
3692
3693 CHARACTER(len=*), PARAMETER :: routinen = 'qs_ot_get_derivative_diag'
3694 REAL(kind=dp), PARAMETER :: rone = 1.0_dp, rzero = 0.0_dp
3695
3696 INTEGER :: handle, k, n
3697 TYPE(dbcsr_distribution_type) :: dist
3698
3699 CALL timeset(routinen, handle)
3700
3701 CALL dbcsr_get_info(matrix_x, nfullrows_total=n, nfullcols_total=k)
3702
3703 ! go for the derivative now
3704 ! this de/dc*(dX/dx)*sinp
3705 CALL dbcsr_multiply('N', 'N', rone, matrix_hc, qs_ot_env%matrix_sinp, rzero, matrix_gx)
3706 ! overlap hc*x
3707 CALL dbcsr_multiply('T', 'N', rone, matrix_hc, matrix_x, rzero, qs_ot_env%matrix_buf2)
3708 ! get it in the basis of the eigenvectors
3709 CALL dbcsr_multiply('N', 'N', rone, qs_ot_env%matrix_buf2, qs_ot_env%matrix_r, &
3710 rzero, qs_ot_env%matrix_buf1)
3711 CALL dbcsr_multiply('T', 'N', rone, qs_ot_env%matrix_r, qs_ot_env%matrix_buf1, &
3712 rzero, qs_ot_env%matrix_buf2)
3713
3714 ! get the schur product of O_uv*B_uv
3715 CALL dbcsr_hadamard_product(qs_ot_env%matrix_buf2, qs_ot_env%matrix_sinp_b, &
3716 qs_ot_env%matrix_buf3)
3717
3718 ! overlap hc*c0
3719 CALL dbcsr_multiply('T', 'N', rone, matrix_hc, qs_ot_env%matrix_c0, rzero, &
3720 qs_ot_env%matrix_buf2)
3721 ! get it in the basis of the eigenvectors
3722 CALL dbcsr_multiply('N', 'N', rone, qs_ot_env%matrix_buf2, qs_ot_env%matrix_r, &
3723 rzero, qs_ot_env%matrix_buf1)
3724 CALL dbcsr_multiply('T', 'N', rone, qs_ot_env%matrix_r, qs_ot_env%matrix_buf1, &
3725 rzero, qs_ot_env%matrix_buf2)
3726
3727 CALL dbcsr_hadamard_product(qs_ot_env%matrix_buf2, qs_ot_env%matrix_cosp_b, &
3728 qs_ot_env%matrix_buf4)
3729
3730 ! add the two bs and compute b+b^T
3731 CALL dbcsr_add(qs_ot_env%matrix_buf3, qs_ot_env%matrix_buf4, &
3732 alpha_scalar=rone, beta_scalar=rone)
3733
3734 ! get the b in the eigenvector basis
3735 CALL dbcsr_multiply('N', 'T', rone, qs_ot_env%matrix_buf3, qs_ot_env%matrix_r, &
3736 rzero, qs_ot_env%matrix_buf1)
3737 CALL dbcsr_multiply('N', 'N', rone, qs_ot_env%matrix_r, qs_ot_env%matrix_buf1, &
3738 rzero, qs_ot_env%matrix_buf3)
3739 CALL dbcsr_get_info(qs_ot_env%matrix_buf3, distribution=dist)
3740 CALL dbcsr_transposed(qs_ot_env%matrix_buf1, qs_ot_env%matrix_buf3, &
3741 shallow_data_copy=.false., use_distribution=dist, &
3742 transpose_distribution=.false.)
3743 CALL dbcsr_add(qs_ot_env%matrix_buf3, qs_ot_env%matrix_buf1, &
3744 alpha_scalar=rone, beta_scalar=rone)
3745
3746 ! and add to the derivative
3747 CALL dbcsr_multiply('N', 'N', rone, matrix_sx, qs_ot_env%matrix_buf3, &
3748 rone, matrix_gx)
3749 CALL timestop(handle)
3750
3751 END SUBROUTINE qs_ot_get_derivative_diag
3752
3753! **************************************************************************************************
3754!> \brief compute the derivative of the taylor expansion below
3755!> \param matrix_hc ...
3756!> \param matrix_x ...
3757!> \param matrix_sx ...
3758!> \param matrix_gx ...
3759!> \param qs_ot_env ...
3760! **************************************************************************************************
3761 SUBROUTINE qs_ot_get_derivative_taylor(matrix_hc, matrix_x, matrix_sx, matrix_gx, &
3762 qs_ot_env)
3763
3764 TYPE(dbcsr_type), POINTER :: matrix_hc, matrix_x, matrix_sx, matrix_gx
3765 TYPE(qs_ot_type) :: qs_ot_env
3766
3767 CHARACTER(len=*), PARAMETER :: routinen = 'qs_ot_get_derivative_taylor'
3768 REAL(kind=dp), PARAMETER :: rone = 1.0_dp, rzero = 0.0_dp
3769
3770 INTEGER :: handle, i, k, n
3771 REAL(kind=dp) :: cosfactor, sinfactor
3772 TYPE(dbcsr_distribution_type) :: dist
3773 TYPE(dbcsr_type), POINTER :: matrix_left, matrix_right
3774
3775 CALL timeset(routinen, handle)
3776
3777 CALL dbcsr_get_info(matrix_x, nfullrows_total=n, nfullcols_total=k)
3778
3779 ! go for the derivative now
3780 ! this de/dc*(dX/dx)*sinp i.e. zeroth order
3781 CALL dbcsr_multiply('N', 'N', rone, matrix_hc, qs_ot_env%matrix_sinp, rzero, matrix_gx)
3782
3783 IF (qs_ot_env%taylor_order <= 0) THEN
3784 CALL timestop(handle)
3785 RETURN
3786 END IF
3787
3788 ! we store the matrix that will multiply sx in matrix_r
3789 CALL dbcsr_set(qs_ot_env%matrix_r, rzero)
3790
3791 ! just better names for matrix_cosp_b and matrix_sinp_b (they are buffer space here)
3792 matrix_left => qs_ot_env%matrix_cosp_b
3793 matrix_right => qs_ot_env%matrix_sinp_b
3794
3795 ! overlap hc*x and add its transpose to matrix_left
3796 CALL dbcsr_multiply('T', 'N', rone, matrix_hc, matrix_x, rzero, matrix_left)
3797 CALL dbcsr_get_info(matrix_left, distribution=dist)
3798 CALL dbcsr_transposed(qs_ot_env%matrix_buf1, matrix_left, &
3799 shallow_data_copy=.false., use_distribution=dist, &
3800 transpose_distribution=.false.)
3801 CALL dbcsr_add(matrix_left, qs_ot_env%matrix_buf1, &
3802 alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
3803 CALL dbcsr_copy(matrix_right, matrix_left)
3804
3805 ! first order
3806 sinfactor = -1.0_dp/(2.0_dp*3.0_dp)
3807 CALL dbcsr_add(qs_ot_env%matrix_r, matrix_left, &
3808 alpha_scalar=1.0_dp, beta_scalar=sinfactor)
3809
3810 ! M
3811 ! OM+MO
3812 ! OOM+OMO+MOO
3813 ! ...
3814 DO i = 2, qs_ot_env%taylor_order
3815 sinfactor = sinfactor*(-1.0_dp)/real(2*i*(2*i + 1), kind=dp)
3816 CALL dbcsr_multiply('N', 'N', rone, qs_ot_env%matrix_p, matrix_left, rzero, qs_ot_env%matrix_buf1)
3817 CALL dbcsr_multiply('N', 'N', rone, matrix_right, qs_ot_env%matrix_p, rzero, matrix_left)
3818 CALL dbcsr_copy(matrix_right, matrix_left)
3819 CALL dbcsr_add(matrix_left, qs_ot_env%matrix_buf1, &
3820 1.0_dp, 1.0_dp)
3821 CALL dbcsr_add(qs_ot_env%matrix_r, matrix_left, &
3822 alpha_scalar=1.0_dp, beta_scalar=sinfactor)
3823 END DO
3824
3825 ! overlap hc*c0 and add its transpose to matrix_left
3826 CALL dbcsr_multiply('T', 'N', rone, matrix_hc, qs_ot_env%matrix_c0, rzero, matrix_left)
3827 CALL dbcsr_get_info(matrix_left, distribution=dist)
3828 CALL dbcsr_transposed(qs_ot_env%matrix_buf1, matrix_left, &
3829 shallow_data_copy=.false., use_distribution=dist, &
3830 transpose_distribution=.false.)
3831 CALL dbcsr_add(matrix_left, qs_ot_env%matrix_buf1, 1.0_dp, 1.0_dp)
3832 CALL dbcsr_copy(matrix_right, matrix_left)
3833
3834 ! first order
3835 cosfactor = -1.0_dp/(1.0_dp*2.0_dp)
3836 CALL dbcsr_add(qs_ot_env%matrix_r, matrix_left, &
3837 alpha_scalar=1.0_dp, beta_scalar=cosfactor)
3838
3839 ! M
3840 ! OM+MO
3841 ! OOM+OMO+MOO
3842 ! ...
3843 DO i = 2, qs_ot_env%taylor_order
3844 cosfactor = cosfactor*(-1.0_dp)/real(2*i*(2*i - 1), kind=dp)
3845 CALL dbcsr_multiply('N', 'N', rone, qs_ot_env%matrix_p, matrix_left, rzero, qs_ot_env%matrix_buf1)
3846 CALL dbcsr_multiply('N', 'N', rone, matrix_right, qs_ot_env%matrix_p, rzero, matrix_left)
3847 CALL dbcsr_copy(matrix_right, matrix_left)
3848 CALL dbcsr_add(matrix_left, qs_ot_env%matrix_buf1, 1.0_dp, 1.0_dp)
3849 CALL dbcsr_add(qs_ot_env%matrix_r, matrix_left, &
3850 alpha_scalar=1.0_dp, beta_scalar=cosfactor)
3851 END DO
3852
3853 ! and add to the derivative
3854 CALL dbcsr_multiply('N', 'N', rone, matrix_sx, qs_ot_env%matrix_r, rone, matrix_gx)
3855
3856 CALL timestop(handle)
3857
3858 END SUBROUTINE qs_ot_get_derivative_taylor
3859
3860! *************************************************************************************************
3861!> \brief computes a taylor expansion.
3862!> \param qs_ot_env ...
3863! **************************************************************************************************
3864 SUBROUTINE qs_ot_p2m_taylor(qs_ot_env)
3865 TYPE(qs_ot_type) :: qs_ot_env
3866
3867 CHARACTER(len=*), PARAMETER :: routinen = 'qs_ot_p2m_taylor'
3868 REAL(kind=dp), PARAMETER :: rone = 1.0_dp, rzero = 0.0_dp
3869
3870 INTEGER :: handle, i, k
3871 REAL(kind=dp) :: cosfactor, sinfactor
3872
3873 CALL timeset(routinen, handle)
3874
3875 ! zeroth order
3876 CALL dbcsr_set(qs_ot_env%matrix_cosp, rzero)
3877 CALL dbcsr_set(qs_ot_env%matrix_sinp, rzero)
3878 CALL dbcsr_add_on_diag(qs_ot_env%matrix_cosp, rone)
3879 CALL dbcsr_add_on_diag(qs_ot_env%matrix_sinp, rone)
3880
3881 IF (qs_ot_env%taylor_order <= 0) THEN
3882 CALL timestop(handle)
3883 RETURN
3884 END IF
3885
3886 ! first order
3887 cosfactor = -1.0_dp/(1.0_dp*2.0_dp)
3888 sinfactor = -1.0_dp/(2.0_dp*3.0_dp)
3889 CALL dbcsr_add(qs_ot_env%matrix_cosp, qs_ot_env%matrix_p, alpha_scalar=1.0_dp, beta_scalar=cosfactor)
3890 CALL dbcsr_add(qs_ot_env%matrix_sinp, qs_ot_env%matrix_p, alpha_scalar=1.0_dp, beta_scalar=sinfactor)
3891 IF (qs_ot_env%taylor_order <= 1) THEN
3892 CALL timestop(handle)
3893 RETURN
3894 END IF
3895
3896 ! other orders
3897 CALL dbcsr_get_info(qs_ot_env%matrix_p, nfullrows_total=k)
3898 CALL dbcsr_copy(qs_ot_env%matrix_r, qs_ot_env%matrix_p)
3899
3900 DO i = 2, qs_ot_env%taylor_order
3901 ! new power of p
3902 CALL dbcsr_multiply('N', 'N', rone, qs_ot_env%matrix_p, qs_ot_env%matrix_r, &
3903 rzero, qs_ot_env%matrix_buf1)
3904 CALL dbcsr_copy(qs_ot_env%matrix_r, qs_ot_env%matrix_buf1)
3905 ! add to the taylor expansion so far
3906 cosfactor = cosfactor*(-1.0_dp)/real(2*i*(2*i - 1), kind=dp)
3907 sinfactor = sinfactor*(-1.0_dp)/real(2*i*(2*i + 1), kind=dp)
3908 CALL dbcsr_add(qs_ot_env%matrix_cosp, qs_ot_env%matrix_r, &
3909 alpha_scalar=1.0_dp, beta_scalar=cosfactor)
3910 CALL dbcsr_add(qs_ot_env%matrix_sinp, qs_ot_env%matrix_r, &
3911 alpha_scalar=1.0_dp, beta_scalar=sinfactor)
3912 END DO
3913
3914 CALL timestop(handle)
3915
3916 END SUBROUTINE qs_ot_p2m_taylor
3917
3918! **************************************************************************************************
3919!> \brief given p, computes - eigenstuff (matrix_r,evals)
3920!> - cos(p^0.5),p^(-0.5)*sin(p^0.5)
3921!> - the real b matrices, needed for the derivatives of these guys
3922!> cosp_b_ij=(1/(2pii) * int(cos(z^1/2)/((z-eval(i))*(z-eval(j))))
3923!> sinp_b_ij=(1/(2pii) * int(z^(-1/2)*sin(z^1/2)/((z-eval(i))*(z-eval(j))))
3924!> \param qs_ot_env ...
3925! **************************************************************************************************
3926 SUBROUTINE qs_ot_p2m_diag(qs_ot_env)
3927
3928 TYPE(qs_ot_type) :: qs_ot_env
3929
3930 CHARACTER(len=*), PARAMETER :: routinen = 'qs_ot_p2m_diag'
3931 REAL(kind=dp), PARAMETER :: rone = 1.0_dp, rzero = 0.0_dp
3932
3933 INTEGER :: col, col_offset, col_size, handle, i, j, &
3934 k, row, row_offset, row_size
3935 REAL(dp), DIMENSION(:, :), POINTER :: block
3936 REAL(kind=dp) :: a, b
3937 TYPE(dbcsr_iterator_type) :: iter
3938
3939 CALL timeset(routinen, handle)
3940
3941 CALL dbcsr_get_info(qs_ot_env%matrix_p, nfullrows_total=k)
3942 CALL dbcsr_copy(qs_ot_env%matrix_buf1, qs_ot_env%matrix_p)
3943 CALL cp_dbcsr_syevd(qs_ot_env%matrix_buf1, qs_ot_env%matrix_r, qs_ot_env%evals, &
3944 qs_ot_env%para_env, qs_ot_env%blacs_env)
3945 DO i = 1, k
3946 qs_ot_env%evals(i) = max(0.0_dp, qs_ot_env%evals(i))
3947 END DO
3948
3949 !$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(i) SHARED(k,qs_ot_env)
3950 DO i = 1, k
3951 qs_ot_env%dum(i) = cos(sqrt(qs_ot_env%evals(i)))
3952 END DO
3953 CALL dbcsr_copy(qs_ot_env%matrix_buf1, qs_ot_env%matrix_r)
3954 CALL dbcsr_scale_by_vector(qs_ot_env%matrix_buf1, alpha=qs_ot_env%dum, side='right')
3955 CALL dbcsr_multiply('N', 'T', rone, qs_ot_env%matrix_r, qs_ot_env%matrix_buf1, &
3956 rzero, qs_ot_env%matrix_cosp)
3957
3958 !$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(i) SHARED(k,qs_ot_env)
3959 DO i = 1, k
3960 qs_ot_env%dum(i) = qs_ot_sinc(sqrt(qs_ot_env%evals(i)))
3961 END DO
3962 CALL dbcsr_copy(qs_ot_env%matrix_buf1, qs_ot_env%matrix_r)
3963 CALL dbcsr_scale_by_vector(qs_ot_env%matrix_buf1, alpha=qs_ot_env%dum, side='right')
3964 CALL dbcsr_multiply('N', 'T', rone, qs_ot_env%matrix_r, qs_ot_env%matrix_buf1, &
3965 rzero, qs_ot_env%matrix_sinp)
3966
3967 CALL dbcsr_copy(qs_ot_env%matrix_cosp_b, qs_ot_env%matrix_cosp)
3968 CALL dbcsr_iterator_start(iter, qs_ot_env%matrix_cosp_b)
3969 DO WHILE (dbcsr_iterator_blocks_left(iter))
3970 CALL dbcsr_iterator_next_block(iter, row, col, block, &
3971 row_size=row_size, col_size=col_size, &
3972 row_offset=row_offset, col_offset=col_offset)
3973 DO j = 1, col_size
3974 DO i = 1, row_size
3975 a = (sqrt(qs_ot_env%evals(row_offset + i - 1)) &
3976 - sqrt(qs_ot_env%evals(col_offset + j - 1)))/2.0_dp
3977 b = (sqrt(qs_ot_env%evals(row_offset + i - 1)) &
3978 + sqrt(qs_ot_env%evals(col_offset + j - 1)))/2.0_dp
3979 block(i, j) = -0.5_dp*qs_ot_sinc(a)*qs_ot_sinc(b)
3980 END DO
3981 END DO
3982 END DO
3983 CALL dbcsr_iterator_stop(iter)
3984
3985 CALL dbcsr_copy(qs_ot_env%matrix_sinp_b, qs_ot_env%matrix_sinp)
3986 CALL dbcsr_iterator_start(iter, qs_ot_env%matrix_sinp_b)
3987 DO WHILE (dbcsr_iterator_blocks_left(iter))
3988 CALL dbcsr_iterator_next_block(iter, row, col, block, &
3989 row_size=row_size, col_size=col_size, &
3990 row_offset=row_offset, col_offset=col_offset)
3991 DO j = 1, col_size
3992 DO i = 1, row_size
3993 a = sqrt(qs_ot_env%evals(row_offset + i - 1))
3994 b = sqrt(qs_ot_env%evals(col_offset + j - 1))
3995 block(i, j) = qs_ot_sincf(a, b)
3996 END DO
3997 END DO
3998 END DO
3999 CALL dbcsr_iterator_stop(iter)
4000
4001 CALL timestop(handle)
4002
4003 END SUBROUTINE qs_ot_p2m_diag
4004
4005! **************************************************************************************************
4006!> \brief diagonalize Hermitian P and build cos(sqrt(P)), sinc(sqrt(P)), and Frechet kernels
4007!> \param qs_ot_env complex STRICT channel state
4008! **************************************************************************************************
4009 SUBROUTINE qs_ot_p2m_diag_complex(qs_ot_env)
4010 TYPE(qs_ot_type) :: qs_ot_env
4011
4012 CHARACTER(len=*), PARAMETER :: routinen = 'qs_ot_p2m_diag_complex'
4013
4014 INTEGER :: col, col_offset, col_size, handle, i, j, &
4015 k, row, row_offset, row_size
4016 REAL(dp), DIMENSION(:, :), POINTER :: block
4017 REAL(kind=dp) :: a, b
4018 TYPE(dbcsr_iterator_type) :: iter
4019
4020 CALL timeset(routinen, handle)
4021
4022 CALL dbcsr_get_info(qs_ot_env%matrix_p, nfullrows_total=k)
4023 CALL cp_dbcsr_heevd(matrix_re=qs_ot_env%matrix_p, matrix_im=qs_ot_env%matrix_p_im, &
4024 eigenvectors_re=qs_ot_env%matrix_r, eigenvectors_im=qs_ot_env%matrix_r_im, &
4025 eigenvalues=qs_ot_env%evals, para_env=qs_ot_env%para_env, &
4026 blacs_env=qs_ot_env%blacs_env)
4027 DO i = 1, k
4028 qs_ot_env%evals(i) = max(0.0_dp, qs_ot_env%evals(i))
4029 END DO
4030
4031 DO i = 1, k
4032 qs_ot_env%dum(i) = cos(sqrt(qs_ot_env%evals(i)))
4033 END DO
4034 CALL dbcsr_copy(qs_ot_env%matrix_buf1, qs_ot_env%matrix_r)
4035 CALL dbcsr_copy(qs_ot_env%matrix_buf1_im, qs_ot_env%matrix_r_im)
4036 CALL dbcsr_scale_by_vector(qs_ot_env%matrix_buf1, alpha=qs_ot_env%dum, side='right')
4037 CALL dbcsr_scale_by_vector(qs_ot_env%matrix_buf1_im, alpha=qs_ot_env%dum, side='right')
4038 CALL qs_ot_complex_multiply('N', 'C', qs_ot_env%matrix_buf1, qs_ot_env%matrix_buf1_im, &
4039 qs_ot_env%matrix_r, qs_ot_env%matrix_r_im, &
4040 qs_ot_env%matrix_cosp, qs_ot_env%matrix_cosp_im, &
4041 qs_ot_env%matrix_buf2)
4042
4043 DO i = 1, k
4044 qs_ot_env%dum(i) = qs_ot_sinc(sqrt(qs_ot_env%evals(i)))
4045 END DO
4046 CALL dbcsr_copy(qs_ot_env%matrix_buf1, qs_ot_env%matrix_r)
4047 CALL dbcsr_copy(qs_ot_env%matrix_buf1_im, qs_ot_env%matrix_r_im)
4048 CALL dbcsr_scale_by_vector(qs_ot_env%matrix_buf1, alpha=qs_ot_env%dum, side='right')
4049 CALL dbcsr_scale_by_vector(qs_ot_env%matrix_buf1_im, alpha=qs_ot_env%dum, side='right')
4050 CALL qs_ot_complex_multiply('N', 'C', qs_ot_env%matrix_buf1, qs_ot_env%matrix_buf1_im, &
4051 qs_ot_env%matrix_r, qs_ot_env%matrix_r_im, &
4052 qs_ot_env%matrix_sinp, qs_ot_env%matrix_sinp_im, &
4053 qs_ot_env%matrix_buf2)
4054
4055 CALL dbcsr_copy(qs_ot_env%matrix_cosp_b, qs_ot_env%matrix_cosp)
4056 CALL dbcsr_iterator_start(iter, qs_ot_env%matrix_cosp_b)
4057 DO WHILE (dbcsr_iterator_blocks_left(iter))
4058 CALL dbcsr_iterator_next_block(iter, row, col, block, &
4059 row_size=row_size, col_size=col_size, &
4060 row_offset=row_offset, col_offset=col_offset)
4061 DO j = 1, col_size
4062 DO i = 1, row_size
4063 a = (sqrt(qs_ot_env%evals(row_offset + i - 1)) &
4064 - sqrt(qs_ot_env%evals(col_offset + j - 1)))/2.0_dp
4065 b = (sqrt(qs_ot_env%evals(row_offset + i - 1)) &
4066 + sqrt(qs_ot_env%evals(col_offset + j - 1)))/2.0_dp
4067 block(i, j) = -0.5_dp*qs_ot_sinc(a)*qs_ot_sinc(b)
4068 END DO
4069 END DO
4070 END DO
4071 CALL dbcsr_iterator_stop(iter)
4072
4073 CALL dbcsr_copy(qs_ot_env%matrix_sinp_b, qs_ot_env%matrix_sinp)
4074 CALL dbcsr_iterator_start(iter, qs_ot_env%matrix_sinp_b)
4075 DO WHILE (dbcsr_iterator_blocks_left(iter))
4076 CALL dbcsr_iterator_next_block(iter, row, col, block, &
4077 row_size=row_size, col_size=col_size, &
4078 row_offset=row_offset, col_offset=col_offset)
4079 DO j = 1, col_size
4080 DO i = 1, row_size
4081 a = sqrt(qs_ot_env%evals(row_offset + i - 1))
4082 b = sqrt(qs_ot_env%evals(col_offset + j - 1))
4083 block(i, j) = qs_ot_sincf(a, b)
4084 END DO
4085 END DO
4086 END DO
4087 CALL dbcsr_iterator_stop(iter)
4088
4089 CALL timestop(handle)
4090
4091 END SUBROUTINE qs_ot_p2m_diag_complex
4092
4093! **************************************************************************************************
4094!> \brief computes sin(x)/x for all values of the argument
4095!> \param x ...
4096!> \return ...
4097! **************************************************************************************************
4098 FUNCTION qs_ot_sinc(x)
4099
4100 REAL(kind=dp), INTENT(IN) :: x
4101 REAL(kind=dp) :: qs_ot_sinc
4102
4103 REAL(kind=dp), PARAMETER :: q1 = 1.0_dp, q2 = -q1/(2.0_dp*3.0_dp), q3 = -q2/(4.0_dp*5.0_dp), &
4104 q4 = -q3/(6.0_dp*7.0_dp), q5 = -q4/(8.0_dp*9.0_dp), q6 = -q5/(10.0_dp*11.0_dp), &
4105 q7 = -q6/(12.0_dp*13.0_dp), q8 = -q7/(14.0_dp*15.0_dp), q9 = -q8/(16.0_dp*17.0_dp), &
4106 q10 = -q9/(18.0_dp*19.0_dp)
4107
4108 REAL(kind=dp) :: y
4109
4110 IF (abs(x) > 0.5_dp) THEN
4111 qs_ot_sinc = sin(x)/x
4112 ELSE
4113 y = x*x
4114 qs_ot_sinc = q1 + y*(q2 + y*(q3 + y*(q4 + y*(q5 + y*(q6 + y*(q7 + y*(q8 + y*(q9 + y*(q10)))))))))
4115 END IF
4116 END FUNCTION qs_ot_sinc
4117
4118! **************************************************************************************************
4119!> \brief computes (1/(x^2-y^2))*(sinc(x)-sinc(y)) for all positive values of the arguments
4120!> \param xa ...
4121!> \param ya ...
4122!> \return ...
4123! **************************************************************************************************
4124 FUNCTION qs_ot_sincf(xa, ya)
4125
4126 REAL(kind=dp), INTENT(IN) :: xa, ya
4127 REAL(kind=dp) :: qs_ot_sincf
4128
4129 INTEGER :: i
4130 REAL(kind=dp) :: a, b, rs, sf, x, xs, y, ybx, ybxs
4131
4132 ! this is currently a limit of the routine, could be removed rather easily
4133 IF (xa < 0) cpabort("x is negative")
4134 IF (ya < 0) cpabort("y is negative")
4135
4136 IF (xa < ya) THEN
4137 x = ya
4138 y = xa
4139 ELSE
4140 x = xa
4141 y = ya
4142 END IF
4143
4144 IF (x < 0.5_dp) THEN ! use series, keeping in mind that x,y,x+y,x-y can all be zero
4145
4146 qs_ot_sincf = 0.0_dp
4147 IF (x > 0.0_dp) THEN
4148 ybx = y/x
4149 ELSE ! should be irrelevant !?
4150 ybx = 0.0_dp
4151 END IF
4152
4153 sf = -1.0_dp/((1.0_dp + ybx)*6.0_dp)
4154 rs = 1.0_dp
4155 ybxs = ybx
4156 xs = 1.0_dp
4157
4158 DO i = 1, 10
4159 qs_ot_sincf = qs_ot_sincf + sf*rs*xs*(1.0_dp + ybxs)
4160 sf = -sf/(real((2*i + 2), dp)*real((2*i + 3), dp))
4161 rs = rs + ybxs
4162 ybxs = ybxs*ybx
4163 xs = xs*x*x
4164 END DO
4165
4166 ELSE ! no series expansion
4167 IF (x - y > 0.1_dp) THEN ! safe to use the normal form
4168 qs_ot_sincf = (qs_ot_sinc(x) - qs_ot_sinc(y))/((x + y)*(x - y))
4169 ELSE
4170 a = (x + y)/2.0_dp
4171 b = (x - y)/2.0_dp ! might be close to zero
4172 ! y (=(a-b)) can not be close to zero since it is close to x>0.5
4173 qs_ot_sincf = (qs_ot_sinc(b)*cos(a) - qs_ot_sinc(a)*cos(b))/(2*x*y)
4174 END IF
4175 END IF
4176
4177 END FUNCTION qs_ot_sincf
4178
4179END MODULE qs_ot
cp_dbcsr_api::dbcsr_create
Definition cp_dbcsr_api.F:194
mathlib::gemm_square
Definition mathlib.F:90
preconditioner::apply_preconditioner
Definition preconditioner.F:115
arnoldi_api
arnoldi iteration using dbcsr
Definition arnoldi_api.F:16
arnoldi_api::arnoldi_extremal
subroutine, public arnoldi_extremal(matrix_a, max_ev, min_ev, converged, threshold, max_iter)
simple wrapper to estimate extremal eigenvalues with arnoldi, using the old lanczos interface this hi...
Definition arnoldi_api.F:254
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_release_p
subroutine, public dbcsr_release_p(matrix)
...
Definition cp_dbcsr_api.F:222
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_init_p
subroutine, public dbcsr_init_p(matrix)
...
Definition cp_dbcsr_api.F:207
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_get_occupation
real(kind=dp) function, public dbcsr_get_occupation(matrix)
...
Definition cp_dbcsr_api.F:893
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_cholesky
Interface to (sca)lapack for the Cholesky based procedures.
Definition cp_dbcsr_cholesky.F:17
cp_dbcsr_cholesky::cp_dbcsr_cholesky_decompose
subroutine, public cp_dbcsr_cholesky_decompose(matrix, n, para_env, blacs_env)
used to replace a symmetric positive def. matrix M with its cholesky decomposition U: M = U^T * U,...
Definition cp_dbcsr_cholesky.F:61
cp_dbcsr_cholesky::cp_dbcsr_cholesky_restore
subroutine, public cp_dbcsr_cholesky_restore(matrix, neig, matrixb, matrixout, op, pos, transa, para_env, blacs_env)
...
Definition cp_dbcsr_cholesky.F:176
cp_dbcsr_cholesky::cp_dbcsr_cholesky_invert
subroutine, public cp_dbcsr_cholesky_invert(matrix, n, para_env, blacs_env, uplo_to_full)
used to replace the cholesky decomposition by the inverse
Definition cp_dbcsr_cholesky.F:114
cp_dbcsr_contrib
Definition cp_dbcsr_contrib.F:8
cp_dbcsr_contrib::dbcsr_gershgorin_norm
real(dp) function, public dbcsr_gershgorin_norm(matrix)
Compute the gershgorin norm of a dbcsr matrix.
Definition cp_dbcsr_contrib.F:167
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:331
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:127
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:486
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_diag::cp_dbcsr_syevd
subroutine, public cp_dbcsr_syevd(matrix, eigenvectors, eigenvalues, para_env, blacs_env)
...
Definition cp_dbcsr_diag.F:67
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
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:1882
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
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:27
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_types
orbital transformations
Definition qs_ot_types.F:15
qs_ot
orbital transformations
Definition qs_ot.F:15
qs_ot::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.F:2294
qs_ot::qs_ot_get_p
subroutine, public qs_ot_get_p(matrix_x, matrix_sx, qs_ot_env)
computes p=x*S*x and the matrix functionals related matrices
Definition qs_ot.F:2643
qs_ot::qs_ot_density_secant_hessian
subroutine, public qs_ot_density_secant_hessian(density_step, hamiltonian_step, density_modes, correction, valid, density_norm_sq, response_work)
project a self-adjoint density/Hamiltonian secant onto density-response modes
Definition qs_ot.F:527
qs_ot::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.F:1986
qs_ot::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.F:324
qs_ot::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.F:2690
qs_ot::qs_ot_symmetric_sr1_update
pure subroutine, public qs_ot_symmetric_sr1_update(matrix, step, response, updated_matrix, valid, relative_tolerance)
add one accepted symmetric response secant to a reference Hessian
Definition qs_ot.F:467
qs_ot::qs_ot_fixed_n_projector_frechet
subroutine, public qs_ot_fixed_n_projector_frechet(chc, dchc, occupation, kpoint_weight, response_weight, fixed_n_weight_sum, projector_derivative, density_factor)
fixed-N Frechet derivative of a smooth occupation projector
Definition qs_ot.F:872
qs_ot::qs_ot_fixed_n_schur_block
subroutine, public qs_ot_fixed_n_schur_block(rotation_hessian, rayleigh_response, response_weight, rotation_gradient, energy_gradient, schur_block, coupling_vector, schur_rhs)
local block of the fixed-N rotation/energy Schur complement
Definition qs_ot.F:229
qs_ot::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.F:964
qs_ot::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.F:156
qs_ot::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.F:102
qs_ot::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.F:3491
qs_ot::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.F:2770
qs_ot::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.F:389
qs_ot::qs_ot_fixed_n_energy_hessian
pure subroutine, public qs_ot_fixed_n_energy_hessian(response_weight, fixed_n_weight_sum, hessian)
dense fixed-N occupation Hessian in auxiliary-energy coordinates
Definition qs_ot.F:183
qs_ot::qs_ot_prepare_complex_tangent_metric
subroutine, public qs_ot_prepare_complex_tangent_metric(qs_ot_env, preconditioner_rejected)
Prepare the inverse metric used to project a complex STRICT gradient. An unusable preconditioner is d...
Definition qs_ot.F:3419
qs_ot::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.F:768
qs_ot::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.F:1208
qs_ot::qs_ot_get_derivative_ref
subroutine, public qs_ot_get_derivative_ref(matrix_hc, matrix_x, matrix_sx, matrix_gx, qs_ot_env)
...
Definition qs_ot.F:2234
qs_ot::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.F:1241
qs_ot::qs_ot_get_orbitals
subroutine, public qs_ot_get_orbitals(matrix_c, matrix_x, qs_ot_env)
c=(c0*cos(p^0.5)+x*sin(p^0.5)*p^(-0.5)) x rot_mat_u this assumes that x is already ortho to S*C0,...
Definition qs_ot.F:3207
qs_ot::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.F:592
qs_ot::qs_ot_rot_mat_derivative
subroutine, public qs_ot_rot_mat_derivative(qs_ot_env)
computes the derivative fields with respect to rot_mat_x
Definition qs_ot.F:2988
qs_ot::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.F:130
qs_ot::qs_ot_get_orbitals_ref
subroutine, public qs_ot_get_orbitals_ref(matrix_c, matrix_s, matrix_x, matrix_sx, matrix_gx_old, matrix_dx, qs_ot_env, qs_ot_env1)
...
Definition qs_ot.F:1881
qs_ot::qs_ot_get_derivative
subroutine, public qs_ot_get_derivative(matrix_hc, matrix_x, matrix_sx, matrix_gx, qs_ot_env)
this routines computes dE/dx=dx, with dx ortho to sc0 needs dE/dC=hc,C0,X,SX,p if preconditioned it w...
Definition qs_ot.F:3332
qs_ot::qs_ot_new_preconditioner
subroutine, public qs_ot_new_preconditioner(qs_ot_env, preconditioner)
gets ready to use the preconditioner/ or renew the preconditioner only keeps a pointer to the precond...
Definition qs_ot.F:1371
qs_ot::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.F:270
qs_ot::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.F:668
qs_ot::qs_ot_generate_rotation
subroutine, public qs_ot_generate_rotation(qs_ot_env)
computes the rotation matrix rot_mat_u that is associated to a given rot_mat_x using rot_mat_u=exp(ro...
Definition qs_ot.F:2923
qs_ot::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)
compute P=X^H*S*X and the STRICT matrix functions for a complex K-point channel
Definition qs_ot.F:2892
qs_ot::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.F:3257
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_type
Definition cp_dbcsr_api.F:176
message_passing::mp_comm_type
Definition message_passing.F:179
preconditioner_types::preconditioner_type
Definition preconditioner_types.F:60
qs_ot_types::qs_ot_type
Definition qs_ot_types.F:139