Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief 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 : ! **************************************************************************************************
15 : MODULE 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 :
94 : CONTAINS
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 152 : 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 152 : COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: eigenvectors
106 : INTEGER :: n
107 152 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eigenvalues
108 :
109 152 : n = SIZE(rotation_generator, 1)
110 152 : CPASSERT(SIZE(rotation_generator, 2) == n)
111 152 : norm = 0.0_dp
112 152 : IF (n == 0) RETURN
113 912 : ALLOCATE (eigenvectors(n, n), eigenvalues(n))
114 : CALL diag_complex(CMPLX(0.0_dp, 1.0_dp, KIND=dp)*rotation_generator, &
115 6584 : eigenvectors, eigenvalues)
116 1028 : norm = MAXVAL(ABS(eigenvalues))
117 152 : DEALLOCATE (eigenvalues, eigenvectors)
118 :
119 152 : 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 2 : 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 2 : denominator = fixed_n_curvature_sum
138 2 : IF (ABS(denominator) <= EPSILON(denominator)) denominator = local_curvature_sum
139 2 : mu_shift = 0.0_dp
140 2 : IF (ABS(denominator) > EPSILON(denominator)) mu_shift = weighted_energy_response/denominator
141 :
142 2 : 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 2038 : PURE SUBROUTINE qs_ot_fixed_n_energy_gradient( &
154 2038 : rayleigh_energy, energy_coordinate, response_weight, fixed_n_weight_sum, &
155 2038 : 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 2038 : fixed_n_mean = 0.0_dp
165 2038 : IF (ABS(fixed_n_weight_sum) > EPSILON(fixed_n_weight_sum)) THEN
166 1994 : fixed_n_mean = fixed_n_weighted_residual/fixed_n_weight_sum
167 : END IF
168 : gradient(:) = response_weight(:)* &
169 16136 : (fixed_n_mean - (rayleigh_energy(:) - energy_coordinate(:)))
170 :
171 2038 : 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 16 : 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 16 : n = SIZE(response_weight)
190 1008 : hessian(:, :) = 0.0_dp
191 122 : DO i = 1, n
192 122 : hessian(i, i) = response_weight(i)
193 : END DO
194 16 : IF (ABS(fixed_n_weight_sum) > EPSILON(fixed_n_weight_sum)) THEN
195 122 : DO j = 1, n
196 1008 : DO i = 1, n
197 : hessian(i, j) = hessian(i, j) - &
198 992 : response_weight(i)*response_weight(j)/fixed_n_weight_sum
199 : END DO
200 : END DO
201 : END IF
202 2016 : hessian(:, :) = 0.5_dp*(hessian + TRANSPOSE(hessian))
203 :
204 16 : 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 2 : SUBROUTINE qs_ot_fixed_n_schur_block( &
227 4 : rotation_hessian, rayleigh_response, response_weight, rotation_gradient, energy_gradient, &
228 4 : 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 2 : nenergy = SIZE(response_weight)
240 10 : ALLOCATE (response_group(nenergy), coupling_matrix(SIZE(coupling_vector), 1))
241 8 : response_group(:) = 1
242 : CALL qs_ot_fixed_n_multigroup_schur_block( &
243 : rotation_hessian, rayleigh_response, response_weight, response_group, &
244 2 : rotation_gradient, energy_gradient, schur_block, coupling_matrix, schur_rhs)
245 14 : coupling_vector(:) = coupling_matrix(:, 1)
246 2 : DEALLOCATE (response_group, coupling_matrix)
247 :
248 2 : 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 324 : SUBROUTINE qs_ot_fixed_n_multigroup_schur_block( &
268 648 : rotation_hessian, rayleigh_response, response_weight, response_group, &
269 648 : 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 324 : nenergy = SIZE(response_weight)
281 324 : nrotation = SIZE(rotation_gradient)
282 324 : ngroups = SIZE(coupling_matrix, 2)
283 324 : CPASSERT(ngroups > 0)
284 972 : CPASSERT(ALL(SHAPE(rotation_hessian) == [nrotation, nrotation]))
285 972 : CPASSERT(ALL(SHAPE(rayleigh_response) == [nenergy, nrotation]))
286 324 : CPASSERT(SIZE(energy_gradient) == nenergy)
287 324 : CPASSERT(SIZE(response_group) == nenergy)
288 2242 : CPASSERT(ALL(response_group >= 1 .AND. response_group <= ngroups))
289 972 : CPASSERT(ALL(SHAPE(schur_block) == [nrotation, nrotation]))
290 324 : CPASSERT(SIZE(coupling_matrix, 1) == nrotation)
291 324 : CPASSERT(SIZE(schur_rhs) == nrotation)
292 :
293 502992 : schur_block(:, :) = rotation_hessian(:, :)
294 11782 : coupling_matrix(:, :) = 0.0_dp
295 2242 : DO i = 1, nenergy
296 1918 : weight = response_weight(i)
297 : schur_block(:, :) = schur_block(:, :) - weight* &
298 : SPREAD(rayleigh_response(i, :), DIM=2, NCOPIES=nrotation)* &
299 5484538 : SPREAD(rayleigh_response(i, :), DIM=1, NCOPIES=nrotation)
300 1918 : group = response_group(i)
301 : coupling_matrix(:, group) = coupling_matrix(:, group) + &
302 84766 : weight*rayleigh_response(i, :)
303 : END DO
304 1005984 : schur_block(:, :) = 0.5_dp*(schur_block + TRANSPOSE(schur_block))
305 : schur_rhs(:) = rotation_gradient + &
306 10122 : MATMUL(TRANSPOSE(rayleigh_response), energy_gradient)
307 :
308 324 : 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 322 : 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 322 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eigenvalues
333 322 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: eigenvectors, work
334 :
335 322 : n = SIZE(matrix, 1)
336 966 : CPASSERT(ALL(SHAPE(matrix) == [n, n]))
337 322 : CPASSERT(SIZE(rhs, 1) == n)
338 966 : CPASSERT(ALL(SHAPE(solution) == SHAPE(rhs)))
339 524362 : valid = ALL(matrix == matrix) .AND. ALL(rhs == rhs)
340 322 : IF (.NOT. valid) THEN
341 0 : solution(:, :) = 0.0_dp
342 : RETURN
343 : END IF
344 :
345 322 : relative_floor_eff = SQRT(EPSILON(1.0_dp))
346 322 : IF (PRESENT(relative_floor)) relative_floor_eff = MAX(relative_floor, relative_floor_eff)
347 2898 : ALLOCATE (eigenvalues(n), eigenvectors(n, n), work(n, SIZE(rhs, 2)))
348 502846 : eigenvectors(:, :) = 0.5_dp*(matrix + TRANSPOSE(matrix))
349 322 : CALL diamat_all(eigenvectors, eigenvalues)
350 9778 : scale = MAX(1.0_dp, MAXVAL(ABS(eigenvalues)))
351 322 : eigenvalue_floor = relative_floor_eff*scale
352 2237422 : work(:, :) = MATMUL(TRANSPOSE(eigenvectors), rhs)
353 322 : nresolved = 0
354 9778 : DO i = 1, n
355 9778 : IF (ABS(eigenvalues(i)) > eigenvalue_floor) THEN
356 15489 : work(i, :) = work(i, :)/ABS(eigenvalues(i))
357 4939 : nresolved = nresolved + 1
358 : ELSE
359 14475 : work(i, :) = 0.0_dp
360 : END IF
361 : END DO
362 2249846 : solution(:, :) = MATMUL(eigenvectors, work)
363 21516 : valid = nresolved > 0 .AND. ALL(solution == solution)
364 610 : IF (.NOT. valid) solution(:, :) = 0.0_dp
365 322 : DEALLOCATE (eigenvalues, eigenvectors, work)
366 :
367 322 : 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 116 : SUBROUTINE qs_ot_projected_response_update( &
388 116 : 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 116 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eigenvalues, rhs, work
400 116 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: eigenvectors
401 :
402 116 : n = SIZE(reference_hessian, 1)
403 116 : CPASSERT(n > 0)
404 348 : CPASSERT(ALL(SHAPE(reference_hessian) == [n, n]))
405 348 : CPASSERT(ALL(SHAPE(response_correction) == [n, n]))
406 116 : CPASSERT(SIZE(coefficients) == n)
407 116 : IF (PRESENT(projected_gradient)) THEN
408 112 : CPASSERT(SIZE(projected_gradient) == n)
409 : END IF
410 348 : coefficients(:) = 0.0_dp
411 : valid = ALL(reference_hessian == reference_hessian) .AND. &
412 1624 : ALL(response_correction == response_correction)
413 116 : IF (PRESENT(projected_gradient)) THEN
414 336 : valid = valid .AND. ALL(projected_gradient == projected_gradient)
415 : END IF
416 116 : IF (.NOT. valid) RETURN
417 :
418 116 : relative_floor_eff = 1.0E-4_dp
419 116 : IF (PRESENT(relative_floor)) relative_floor_eff = &
420 0 : MAX(relative_floor, SQRT(EPSILON(1.0_dp)))
421 928 : ALLOCATE (eigenvalues(n), eigenvectors(n, n), rhs(n), work(n))
422 : eigenvectors(:, :) = 0.5_dp* &
423 : (reference_hessian + response_correction + &
424 812 : TRANSPOSE(reference_hessian + response_correction))
425 116 : CALL diamat_all(eigenvectors, eigenvalues)
426 1044 : scale = MAX(MAXVAL(ABS(eigenvalues)), MAXVAL(ABS(reference_hessian)))
427 116 : IF (scale <= TINY(scale)) THEN
428 0 : valid = .FALSE.
429 0 : coefficients(:) = 0.0_dp
430 0 : DEALLOCATE (eigenvalues, eigenvectors, rhs, work)
431 0 : RETURN
432 : END IF
433 116 : eigenvalue_floor = relative_floor_eff*scale
434 344 : valid = ALL(eigenvalues > eigenvalue_floor)
435 116 : IF (valid) THEN
436 342 : rhs(:) = reference_hessian(:, 1)
437 338 : IF (PRESENT(projected_gradient)) rhs(:) = projected_gradient
438 114 : work(:) = MATMUL(TRANSPOSE(eigenvectors), rhs)
439 342 : DO i = 1, n
440 342 : work(i) = work(i)/eigenvalues(i)
441 : END DO
442 1482 : coefficients(:) = MATMUL(eigenvectors, work)
443 : valid = ALL(coefficients == coefficients) .AND. &
444 684 : ALL(ABS(coefficients) <= HUGE(1.0_dp))
445 : END IF
446 120 : IF (.NOT. valid) coefficients(:) = 0.0_dp
447 116 : DEALLOCATE (eigenvalues, eigenvectors, rhs, work)
448 :
449 116 : 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 4 : PURE SUBROUTINE qs_ot_symmetric_sr1_update( &
466 4 : 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 2 : REAL(KIND=dp), DIMENSION(SIZE(step)) :: residual
477 :
478 4 : n = SIZE(step)
479 52 : updated_matrix(:, :) = 0.0_dp
480 : valid = SIZE(response) == n .AND. ALL(SHAPE(matrix) == [n, n]) .AND. &
481 20 : ALL(SHAPE(updated_matrix) == [n, n])
482 4 : IF (.NOT. valid) RETURN
483 84 : valid = ALL(matrix == matrix) .AND. ALL(step == step) .AND. ALL(response == response)
484 4 : IF (.NOT. valid) RETURN
485 :
486 52 : updated_matrix(:, :) = 0.5_dp*(matrix + TRANSPOSE(matrix))
487 76 : residual(:) = response - MATMUL(updated_matrix, step)
488 16 : denominator = DOT_PRODUCT(residual, step)
489 16 : residual_norm = SQRT(MAX(0.0_dp, DOT_PRODUCT(residual, residual)))
490 16 : step_norm = SQRT(MAX(0.0_dp, DOT_PRODUCT(step, step)))
491 4 : tolerance = SQRT(EPSILON(1.0_dp))
492 4 : IF (PRESENT(relative_tolerance)) tolerance = MAX(tolerance, relative_tolerance)
493 4 : threshold = tolerance*residual_norm*step_norm
494 : valid = residual_norm > TINY(residual_norm) .AND. step_norm > TINY(step_norm) .AND. &
495 4 : ABS(denominator) > threshold
496 4 : IF (.NOT. valid) RETURN
497 :
498 : updated_matrix(:, :) = updated_matrix + &
499 26 : SPREAD(residual, DIM=2, NCOPIES=n)*SPREAD(residual, DIM=1, NCOPIES=n)/denominator
500 50 : updated_matrix(:, :) = 0.5_dp*(updated_matrix + TRANSPOSE(updated_matrix))
501 26 : 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 10 : SUBROUTINE qs_ot_density_secant_hessian( &
525 10 : 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 10 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: density_overlap, response_overlap
538 :
539 10 : n = SIZE(density_step, 1)
540 10 : nmode = SIZE(density_modes, 3)
541 10 : CPASSERT(n > 0)
542 10 : CPASSERT(SIZE(density_step, 2) == n)
543 30 : CPASSERT(ALL(SHAPE(hamiltonian_step) == [n, n]))
544 10 : CPASSERT(SIZE(density_modes, 1) == n)
545 10 : CPASSERT(SIZE(density_modes, 2) == n)
546 30 : CPASSERT(ALL(SHAPE(correction) == [nmode, nmode]))
547 :
548 70 : correction(:, :) = 0.0_dp
549 10 : valid = .FALSE.
550 10 : density_norm = 0.0_dp
551 10 : density_response_work = 0.0_dp
552 38 : DO j = 1, n
553 122 : DO i = 1, n
554 : density_norm = density_norm + &
555 84 : REAL(CONJG(density_step(i, j))*density_step(i, j), KIND=dp)
556 : density_response_work = density_response_work + &
557 112 : REAL(CONJG(density_step(i, j))*hamiltonian_step(i, j), KIND=dp)
558 : END DO
559 : END DO
560 10 : IF (PRESENT(density_norm_sq)) density_norm_sq = density_norm
561 10 : IF (PRESENT(response_work)) response_work = density_response_work
562 122 : scale = MAXVAL(ABS(density_step))
563 10 : IF (scale <= TINY(1.0_dp)) RETURN
564 10 : IF (density_norm <= 64.0_dp*EPSILON(1.0_dp)*scale*scale .OR. nmode <= 0) RETURN
565 :
566 40 : ALLOCATE (density_overlap(nmode), response_overlap(nmode))
567 30 : DO mode = 1, nmode
568 244 : density_overlap(mode) = SUM(REAL(CONJG(density_step)*density_modes(:, :, mode), KIND=dp))
569 254 : 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 10 : density_norm, density_response_work, density_overlap, response_overlap, correction, valid)
573 10 : 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 172 : SUBROUTINE qs_ot_density_secant_projected_hessian( &
590 172 : 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 344 : REAL(KIND=dp), DIMENSION(SIZE(density_overlap)) :: projected_density_overlap, &
603 172 : projected_response_overlap
604 :
605 172 : nmode = SIZE(density_overlap)
606 172 : CPASSERT(SIZE(response_overlap) == nmode)
607 516 : CPASSERT(ALL(SHAPE(correction) == [nmode, nmode]))
608 :
609 1204 : correction(:, :) = 0.0_dp
610 172 : valid = .FALSE.
611 172 : IF (nmode <= 0 .OR. density_norm_sq <= TINY(1.0_dp)) RETURN
612 172 : IF (density_norm_sq /= density_norm_sq .OR. response_work /= response_work) RETURN
613 172 : IF (ABS(density_norm_sq) > HUGE(1.0_dp) .OR. ABS(response_work) > HUGE(1.0_dp)) RETURN
614 1032 : IF (ANY(density_overlap /= density_overlap) .OR. ANY(response_overlap /= response_overlap)) RETURN
615 1032 : IF (ANY(ABS(density_overlap) > HUGE(1.0_dp)) .OR. &
616 : ANY(ABS(response_overlap) > HUGE(1.0_dp))) RETURN
617 172 : IF (PRESENT(secant_mode) .NEQV. PRESENT(secant_position)) RETURN
618 :
619 516 : projected_density_overlap(:) = density_overlap
620 516 : projected_response_overlap(:) = response_overlap
621 172 : IF (PRESENT(secant_mode)) THEN
622 114 : IF (secant_mode < 1 .OR. secant_mode > nmode) RETURN
623 : IF (secant_position /= secant_position .OR. &
624 114 : ABS(secant_position) > HUGE(1.0_dp) .OR. &
625 114 : 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 112 : projected_density_overlap(secant_mode) = density_norm_sq/secant_position
629 112 : projected_response_overlap(secant_mode) = response_work/secant_position
630 : END IF
631 :
632 170 : inverse_density_norm = 1.0_dp/density_norm_sq
633 510 : DO j = 1, nmode
634 1190 : 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 1020 : projected_density_overlap(i)*projected_density_overlap(j))
640 : END DO
641 : END DO
642 2210 : correction(:, :) = 0.5_dp*(correction + TRANSPOSE(correction))
643 2380 : valid = ALL(correction == correction) .AND. ALL(ABS(correction) <= HUGE(1.0_dp))
644 170 : 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 444 : SUBROUTINE qs_ot_density_tangent( &
666 444 : rotation_generator, occupation, kpoint_weight, rotation_step, &
667 444 : 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 444 : COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: direction, generator_minus, &
677 444 : 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 444 : n = SIZE(rotation_generator, 1)
682 444 : nrotation = n*(n - 1)
683 444 : CPASSERT(n > 0)
684 1332 : CPASSERT(ALL(SHAPE(rotation_generator) == [n, n]))
685 444 : CPASSERT(SIZE(occupation) == n)
686 444 : CPASSERT(SIZE(rotation_step) == nrotation)
687 444 : CPASSERT(SIZE(weighted_occupation_step) == n)
688 1332 : CPASSERT(ALL(SHAPE(density_tangent) == [n, n]))
689 444 : 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 7992 : rotation_plus(n, n), weighted_rotation(n, n))
694 444 : direction(:, :) = CMPLX(0.0_dp, 0.0_dp, KIND=dp)
695 444 : r = 0
696 2870 : DO i = 1, n - 1
697 11726 : DO j = i + 1, n
698 8856 : r = r + 1
699 8856 : direction(i, j) = CMPLX(rotation_step(r), 0.0_dp, KIND=dp)
700 8856 : direction(j, i) = -direction(i, j)
701 8856 : r = r + 1
702 : direction(i, j) = direction(i, j) + &
703 8856 : CMPLX(0.0_dp, rotation_step(r), KIND=dp)
704 : direction(j, i) = direction(j, i) + &
705 11282 : CMPLX(0.0_dp, rotation_step(r), KIND=dp)
706 : END DO
707 : END DO
708 444 : CPASSERT(r == nrotation)
709 :
710 23896 : step_scale = MAX(1.0_dp, MAXVAL(ABS(direction)))
711 444 : step = 1.0E-5_dp/step_scale
712 444 : IF (PRESENT(difference_step)) step = difference_step/step_scale
713 444 : CPASSERT(step > EPSILON(step))
714 23896 : generator_plus(:, :) = rotation_generator + step*direction
715 23896 : generator_minus(:, :) = rotation_generator - step*direction
716 444 : CALL qs_ot_dense_rotation_state(rotation_generator, rotation)
717 444 : CALL qs_ot_dense_rotation_state(generator_plus, rotation_plus)
718 444 : CALL qs_ot_dense_rotation_state(generator_minus, rotation_minus)
719 :
720 23896 : weighted_rotation(:, :) = rotation_plus
721 3314 : DO j = 1, n
722 23896 : weighted_rotation(:, j) = kpoint_weight*occupation(j)*weighted_rotation(:, j)
723 : END DO
724 185274 : mode_ref(:, :) = MATMUL(weighted_rotation, CONJG(TRANSPOSE(rotation_plus)))
725 23896 : weighted_rotation(:, :) = rotation_minus
726 3314 : DO j = 1, n
727 23896 : weighted_rotation(:, j) = kpoint_weight*occupation(j)*weighted_rotation(:, j)
728 : END DO
729 : mode_ref(:, :) = (mode_ref - &
730 210058 : MATMUL(weighted_rotation, CONJG(TRANSPOSE(rotation_minus))))/(2.0_dp*step)
731 :
732 23896 : weighted_rotation(:, :) = rotation
733 3314 : DO j = 1, n
734 23896 : weighted_rotation(:, j) = weighted_occupation_step(j)*weighted_rotation(:, j)
735 : END DO
736 210502 : mode_ref(:, :) = mode_ref + MATMUL(weighted_rotation, CONJG(TRANSPOSE(rotation)))
737 556266 : density_tangent(:, :) = MATMUL(CONJG(TRANSPOSE(rotation)), MATMUL(mode_ref, rotation))
738 47348 : density_tangent(:, :) = 0.5_dp*(density_tangent + CONJG(TRANSPOSE(density_tangent)))
739 :
740 0 : DEALLOCATE (direction, generator_minus, generator_plus, mode_ref, rotation, &
741 444 : rotation_minus, rotation_plus, weighted_rotation)
742 :
743 444 : 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 258 : SUBROUTINE qs_ot_density_secant_orbital_overlaps( &
765 516 : overlap_start_current, occupation_start, occupation_current, hamiltonian_step_start, &
766 516 : hamiltonian_step_current, density_modes, kpoint_weight, density_norm_sq, response_work, &
767 258 : 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 258 : COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: cross_mode
780 : INTEGER :: i, j, mode, n, nmode
781 : REAL(KIND=dp) :: cross_density, density_scale
782 :
783 258 : n = SIZE(occupation_start)
784 258 : nmode = SIZE(density_modes, 3)
785 258 : CPASSERT(n > 0)
786 258 : CPASSERT(SIZE(occupation_current) == n)
787 774 : CPASSERT(ALL(SHAPE(overlap_start_current) == [n, n]))
788 774 : CPASSERT(ALL(SHAPE(hamiltonian_step_start) == [n, n]))
789 774 : CPASSERT(ALL(SHAPE(hamiltonian_step_current) == [n, n]))
790 258 : CPASSERT(SIZE(density_modes, 1) == n .AND. SIZE(density_modes, 2) == n)
791 258 : CPASSERT(SIZE(density_overlap) == nmode .AND. SIZE(response_overlap) == nmode)
792 :
793 258 : density_norm_sq = 0.0_dp
794 258 : response_work = 0.0_dp
795 774 : density_overlap(:) = 0.0_dp
796 774 : response_overlap(:) = 0.0_dp
797 258 : valid = .FALSE.
798 274 : IF (nmode <= 0 .OR. kpoint_weight <= TINY(1.0_dp)) RETURN
799 :
800 : cross_density = 0.0_dp
801 1912 : DO j = 1, n
802 13706 : DO i = 1, n
803 : cross_density = cross_density + occupation_start(i)*occupation_current(j)* &
804 13448 : 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 3566 : 2.0_dp*cross_density)
810 : response_work = 0.0_dp
811 1912 : DO i = 1, n
812 : response_work = response_work + kpoint_weight* &
813 : (occupation_current(i)*REAL(hamiltonian_step_current(i, i), KIND=dp) - &
814 1912 : occupation_start(i)*REAL(hamiltonian_step_start(i, i), KIND=dp))
815 : END DO
816 :
817 1032 : ALLOCATE (cross_mode(n, n))
818 774 : 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 516 : overlap_start_current, "C", cross_mode)
822 3824 : 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 3824 : 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 27670 : AIMAG(hamiltonian_step_current)*AIMAG(density_modes(:, :, mode)))
832 : END DO
833 258 : DEALLOCATE (cross_mode)
834 :
835 : density_scale = kpoint_weight* &
836 3566 : MAX(MAXVAL(ABS(occupation_start)), MAXVAL(ABS(occupation_current)))
837 258 : IF (density_scale <= TINY(1.0_dp)) RETURN
838 258 : 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 2904 : ALL(ABS(response_overlap) <= HUGE(1.0_dp))
845 242 : IF (.NOT. valid) THEN
846 0 : density_norm_sq = 0.0_dp
847 0 : response_work = 0.0_dp
848 0 : density_overlap(:) = 0.0_dp
849 0 : response_overlap(:) = 0.0_dp
850 : END IF
851 :
852 16 : 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 2 : SUBROUTINE qs_ot_fixed_n_projector_frechet( &
870 2 : chc, dchc, occupation, kpoint_weight, response_weight, fixed_n_weight_sum, &
871 2 : 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 2 : 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 2 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eigenvalues, weighted_occupation
886 :
887 2 : n = SIZE(chc, 1)
888 2 : CPASSERT(n > 0)
889 6 : CPASSERT(ALL(SHAPE(chc) == [n, n]))
890 6 : CPASSERT(ALL(SHAPE(dchc) == [n, n]))
891 2 : CPASSERT(SIZE(occupation) == n)
892 2 : CPASSERT(SIZE(response_weight) == n)
893 6 : CPASSERT(ALL(SHAPE(projector_derivative) == [n, n]))
894 :
895 2 : factor = 1.0_dp
896 2 : IF (PRESENT(density_factor)) factor = density_factor
897 26 : projector_derivative(:, :) = CMPLX(0.0_dp, 0.0_dp, KIND=dp)
898 2 : 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 22 : weighted_occupation(n))
902 2 : CALL diag_complex(chc, eigenvectors, eigenvalues)
903 380 : work(:, :) = MATMUL(CONJG(TRANSPOSE(eigenvectors)), MATMUL(dchc, eigenvectors))
904 8 : weighted_occupation(:) = factor*kpoint_weight*occupation(:)
905 :
906 8 : local_weight_sum = SUM(response_weight(:))
907 2 : mu_numerator = 0.0_dp
908 8 : DO i = 1, n
909 : mu_numerator = mu_numerator + response_weight(i)* &
910 8 : 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 2 : fixed_n_weight_sum)
914 :
915 8 : scale = MAX(1.0_dp, MAXVAL(ABS(eigenvalues(:))))
916 2 : gap_tolerance = SQRT(EPSILON(1.0_dp))*scale
917 2 : kernel(:, :) = CMPLX(0.0_dp, 0.0_dp, KIND=dp)
918 8 : DO j = 1, n
919 26 : DO i = 1, n
920 24 : IF (i == j) THEN
921 6 : coefficient = -factor*response_weight(i)
922 : kernel(i, i) = CMPLX(coefficient*(REAL(work(i, i), KIND=dp) - mu_shift), &
923 6 : 0.0_dp, KIND=dp)
924 : ELSE
925 12 : denominator = eigenvalues(i) - eigenvalues(j)
926 12 : IF (ABS(denominator) > gap_tolerance) THEN
927 12 : coefficient = (weighted_occupation(i) - weighted_occupation(j))/denominator
928 : ELSE
929 0 : coefficient = -0.5_dp*factor*(response_weight(i) + response_weight(j))
930 : END IF
931 12 : kernel(i, j) = coefficient*work(i, j)
932 : END IF
933 : END DO
934 : END DO
935 :
936 2 : projector_derivative(:, :) = MATMUL(eigenvectors, &
937 406 : MATMUL(kernel, CONJG(TRANSPOSE(eigenvectors))))
938 : projector_derivative(:, :) = 0.5_dp*(projector_derivative + &
939 50 : CONJG(TRANSPOSE(projector_derivative)))
940 :
941 2 : DEALLOCATE (eigenvectors, kernel, work, eigenvalues, weighted_occupation)
942 :
943 2 : 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 340 : SUBROUTINE qs_ot_finite_rotation_response( &
962 340 : chc, rotation_generator, occupation, kpoint_weight, rotation_gradient, &
963 340 : 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 340 : COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: base_hamiltonian, generator_minus, &
972 340 : generator_plus, rotation
973 : INTEGER :: i, j, n, nrotation, r, s
974 : REAL(KIND=dp) :: step
975 340 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: occupation_scale, rayleigh, &
976 340 : rayleigh_minus, rayleigh_plus
977 340 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: gradient_imag, gradient_minus_imag, &
978 340 : gradient_minus_real, gradient_plus_imag, gradient_plus_real, gradient_real
979 :
980 340 : n = SIZE(chc, 1)
981 340 : nrotation = n*(n - 1)
982 340 : CPASSERT(n > 0)
983 1020 : CPASSERT(ALL(SHAPE(chc) == [n, n]))
984 1020 : CPASSERT(ALL(SHAPE(rotation_generator) == [n, n]))
985 340 : CPASSERT(SIZE(occupation) == n)
986 340 : CPASSERT(SIZE(rotation_gradient) == nrotation)
987 1020 : CPASSERT(ALL(SHAPE(rotation_hessian) == [nrotation, nrotation]))
988 1020 : CPASSERT(ALL(SHAPE(rayleigh_response) == [n, nrotation]))
989 340 : CPASSERT(kpoint_weight > 0.0_dp)
990 :
991 340 : step = 1.0E-4_dp
992 340 : IF (PRESENT(difference_step)) step = difference_step
993 340 : 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 9520 : gradient_plus_real(n, n), gradient_real(n, n))
1000 :
1001 340 : CALL qs_ot_dense_rotation_state(rotation_generator, rotation)
1002 316140 : base_hamiltonian(:, :) = MATMUL(rotation, MATMUL(chc, CONJG(TRANSPOSE(rotation))))
1003 2256 : occupation_scale(:) = 2.0_dp*kpoint_weight*occupation(:)
1004 : CALL qs_ot_dense_rotation_gradient(rotation_generator, base_hamiltonian, occupation_scale, &
1005 340 : gradient_real, gradient_imag, rayleigh)
1006 :
1007 340 : r = 0
1008 1916 : DO i = 1, n - 1
1009 7292 : DO j = i + 1, n
1010 5376 : r = r + 1
1011 5376 : rotation_gradient(r) = gradient_real(i, j)
1012 5376 : r = r + 1
1013 6952 : rotation_gradient(r) = gradient_imag(i, j)
1014 : END DO
1015 : END DO
1016 340 : CPASSERT(r == nrotation)
1017 :
1018 : s = 0
1019 1916 : DO i = 1, n - 1
1020 7292 : DO j = i + 1, n
1021 391272 : generator_plus(:, :) = rotation_generator(:, :)
1022 391272 : generator_minus(:, :) = rotation_generator(:, :)
1023 5376 : generator_plus(i, j) = generator_plus(i, j) + step
1024 5376 : generator_plus(j, i) = generator_plus(j, i) - step
1025 5376 : generator_minus(i, j) = generator_minus(i, j) - step
1026 5376 : generator_minus(j, i) = generator_minus(j, i) + step
1027 5376 : s = s + 1
1028 : CALL qs_ot_dense_rotation_gradient(generator_plus, base_hamiltonian, occupation_scale, &
1029 5376 : gradient_plus_real, gradient_plus_imag, rayleigh_plus)
1030 : CALL qs_ot_dense_rotation_gradient(generator_minus, base_hamiltonian, occupation_scale, &
1031 5376 : 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 5376 : rayleigh_response(:, s))
1036 :
1037 391272 : generator_plus(:, :) = rotation_generator(:, :)
1038 391272 : generator_minus(:, :) = rotation_generator(:, :)
1039 5376 : generator_plus(i, j) = generator_plus(i, j) + CMPLX(0.0_dp, step, KIND=dp)
1040 5376 : generator_plus(j, i) = generator_plus(j, i) + CMPLX(0.0_dp, step, KIND=dp)
1041 5376 : generator_minus(i, j) = generator_minus(i, j) - CMPLX(0.0_dp, step, KIND=dp)
1042 5376 : generator_minus(j, i) = generator_minus(j, i) - CMPLX(0.0_dp, step, KIND=dp)
1043 5376 : s = s + 1
1044 : CALL qs_ot_dense_rotation_gradient(generator_plus, base_hamiltonian, occupation_scale, &
1045 5376 : gradient_plus_real, gradient_plus_imag, rayleigh_plus)
1046 : CALL qs_ot_dense_rotation_gradient(generator_minus, base_hamiltonian, occupation_scale, &
1047 5376 : 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 6952 : rayleigh_response(:, s))
1052 : END DO
1053 : END DO
1054 340 : CPASSERT(s == nrotation)
1055 1235624 : rotation_hessian(:, :) = 0.5_dp*(rotation_hessian + TRANSPOSE(rotation_hessian))
1056 :
1057 0 : DEALLOCATE (base_hamiltonian, generator_minus, generator_plus, rotation, occupation_scale, &
1058 0 : rayleigh, rayleigh_minus, rayleigh_plus, gradient_imag, gradient_minus_imag, &
1059 340 : gradient_minus_real, gradient_plus_imag, gradient_plus_real, gradient_real)
1060 :
1061 340 : 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 10752 : SUBROUTINE qs_ot_pack_rotation_response( &
1076 10752 : gradient_plus_real, gradient_plus_imag, gradient_minus_real, gradient_minus_imag, &
1077 10752 : 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 10752 : n = SIZE(gradient_plus_real, 1)
1088 32256 : CPASSERT(ALL(SHAPE(gradient_plus_real) == [n, n]))
1089 32256 : CPASSERT(ALL(SHAPE(gradient_plus_imag) == [n, n]))
1090 32256 : CPASSERT(ALL(SHAPE(gradient_minus_real) == [n, n]))
1091 32256 : CPASSERT(ALL(SHAPE(gradient_minus_imag) == [n, n]))
1092 10752 : CPASSERT(SIZE(hessian_column) == n*(n - 1))
1093 10752 : CPASSERT(SIZE(rayleigh_column) == n)
1094 :
1095 10752 : r = 0
1096 82536 : DO i = 1, n - 1
1097 385896 : DO j = i + 1, n
1098 303360 : r = r + 1
1099 : hessian_column(r) = &
1100 303360 : (gradient_plus_real(i, j) - gradient_minus_real(i, j))/(2.0_dp*step)
1101 303360 : r = r + 1
1102 : hessian_column(r) = &
1103 375144 : (gradient_plus_imag(i, j) - gradient_minus_imag(i, j))/(2.0_dp*step)
1104 : END DO
1105 : END DO
1106 93288 : rayleigh_column(:) = (rayleigh_plus(:) - rayleigh_minus(:))/(2.0_dp*step)
1107 :
1108 10752 : 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 21844 : SUBROUTINE qs_ot_dense_rotation_gradient( &
1120 21844 : 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 21844 : COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: covector, eigenvectors, frechet, inner, &
1128 21844 : outer, rotation, work
1129 : INTEGER :: i, j, n
1130 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eigenvalues
1131 :
1132 21844 : n = SIZE(rotation_generator, 1)
1133 : ALLOCATE (covector(n, n), eigenvectors(n, n), frechet(n, n), inner(n, n), &
1134 393192 : outer(n, n), rotation(n, n), work(n, n), eigenvalues(n))
1135 : CALL qs_ot_dense_rotation_state(rotation_generator, rotation, &
1136 21844 : eigenvectors=eigenvectors, eigenvalues=eigenvalues)
1137 :
1138 13854208 : work(:, :) = MATMUL(base_hamiltonian, rotation)
1139 1580012 : covector(:, :) = work(:, :)
1140 188832 : DO j = 1, n
1141 1580012 : covector(:, j) = occupation_scale(j)*covector(:, j)
1142 : END DO
1143 :
1144 21844 : work(:, :) = MATMUL(CONJG(TRANSPOSE(rotation)), &
1145 40048144 : MATMUL(base_hamiltonian, rotation))
1146 188832 : DO i = 1, n
1147 188832 : rayleigh(i) = REAL(work(i, i), KIND=dp)
1148 : END DO
1149 :
1150 21844 : inner(:, :) = MATMUL(CONJG(TRANSPOSE(eigenvectors)), &
1151 40048144 : MATMUL(covector, eigenvectors))
1152 188832 : DO j = 1, n
1153 1580012 : DO i = 1, n
1154 1391180 : kernel = CONJG(qs_ot_complex_exp_frechet_kernel(eigenvalues(i), eigenvalues(j)))
1155 1558168 : outer(i, j) = inner(i, j)*kernel
1156 : END DO
1157 : END DO
1158 21844 : frechet(:, :) = MATMUL(eigenvectors, &
1159 40048144 : MATMUL(outer, CONJG(TRANSPOSE(eigenvectors))))
1160 1580012 : gradient_real(:, :) = REAL(frechet, KIND=dp) - TRANSPOSE(REAL(frechet, KIND=dp))
1161 1580012 : gradient_imag(:, :) = AIMAG(frechet) + TRANSPOSE(AIMAG(frechet))
1162 :
1163 21844 : DEALLOCATE (covector, eigenvectors, frechet, inner, outer, rotation, work, eigenvalues)
1164 :
1165 21844 : 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 23516 : 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 23516 : n = SIZE(rotation_generator, 1)
1186 188128 : 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 1666624 : vectors, values)
1189 1666624 : weighted_vectors(:, :) = vectors(:, :)
1190 201030 : DO j = 1, n
1191 : weighted_vectors(:, j) = EXP(CMPLX(0.0_dp, -values(j), KIND=dp))* &
1192 1666624 : weighted_vectors(:, j)
1193 : END DO
1194 29016800 : rotation(:, :) = MATMUL(weighted_vectors, CONJG(TRANSPOSE(vectors)))
1195 1581684 : IF (PRESENT(eigenvectors)) eigenvectors(:, :) = vectors(:, :)
1196 190504 : IF (PRESENT(eigenvalues)) eigenvalues(:) = values(:)
1197 23516 : DEALLOCATE (vectors, weighted_vectors, values)
1198 :
1199 23516 : 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 1503434 : 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 1503434 : l1 = (0.0_dp, -1.0_dp)*e1
1215 1503434 : l2 = (0.0_dp, -1.0_dp)*e2
1216 1503434 : IF (ABS(l1 - l2) > 0.5_dp) THEN
1217 8584 : kernel = (EXP(l1) - EXP(l2))/(l1 - l2)
1218 : ELSE
1219 : x = 1.0_dp
1220 : kernel = 0.0_dp
1221 25412450 : DO i = 1, 16
1222 23917600 : kernel = kernel + x
1223 25412450 : x = x*(l1 - l2)/REAL(i + 1, KIND=dp)
1224 : END DO
1225 1494850 : kernel = kernel*EXP(l2)
1226 : END IF
1227 :
1228 1503434 : 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 1591 : 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 1591 : INTEGER, ALLOCATABLE, DIMENSION(:) :: cols, rows
1249 1591 : INTEGER, DIMENSION(:), POINTER :: col_blk_offset, col_blk_size, &
1250 1591 : row_blk_offset, row_blk_size
1251 : LOGICAL :: found_im, found_out_im, found_out_re, &
1252 : found_re, use_adjoint
1253 1591 : REAL(dp), DIMENSION(:, :), POINTER :: block_in_im, block_in_re, block_out_im, &
1254 1591 : 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 1591 : use_adjoint = .FALSE.
1260 1589 : 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 1591 : max_blocks = 0
1265 1591 : CALL dbcsr_iterator_start(iter, inner_deriv_re)
1266 2645 : DO WHILE (dbcsr_iterator_blocks_left(iter))
1267 1054 : CALL dbcsr_iterator_next_block(iter, row, col)
1268 1054 : max_blocks = max_blocks + 1
1269 : END DO
1270 1591 : CALL dbcsr_iterator_stop(iter)
1271 1591 : CALL dbcsr_iterator_start(iter, inner_deriv_im)
1272 2645 : DO WHILE (dbcsr_iterator_blocks_left(iter))
1273 1054 : CALL dbcsr_iterator_next_block(iter, row, col)
1274 1054 : max_blocks = max_blocks + 1
1275 : END DO
1276 1591 : CALL dbcsr_iterator_stop(iter)
1277 6364 : ALLOCATE (rows(MAX(max_blocks, 1)), cols(MAX(max_blocks, 1)))
1278 1591 : nblocks = 0
1279 :
1280 1591 : CALL dbcsr_iterator_start(iter, inner_deriv_re)
1281 2645 : DO WHILE (dbcsr_iterator_blocks_left(iter))
1282 1054 : CALL dbcsr_iterator_next_block(iter, row, col)
1283 1054 : CALL append_union_block(row, col)
1284 : END DO
1285 1591 : CALL dbcsr_iterator_stop(iter)
1286 1591 : CALL dbcsr_iterator_start(iter, inner_deriv_im)
1287 2645 : DO WHILE (dbcsr_iterator_blocks_left(iter))
1288 1054 : CALL dbcsr_iterator_next_block(iter, row, col)
1289 1054 : CALL append_union_block(row, col)
1290 : END DO
1291 1591 : CALL dbcsr_iterator_stop(iter)
1292 :
1293 : CALL dbcsr_get_info(inner_deriv_re, distribution=dist, row_blk_size=row_blk_size, &
1294 1591 : col_blk_size=col_blk_size)
1295 : CALL dbcsr_create(outer_deriv_re, "outer_deriv_re", dist, dbcsr_type_no_symmetry, &
1296 1591 : row_blk_size, col_blk_size)
1297 : CALL dbcsr_create(outer_deriv_im, "outer_deriv_im", dist, dbcsr_type_no_symmetry, &
1298 1591 : row_blk_size, col_blk_size)
1299 1591 : IF (nblocks > 0) THEN
1300 1054 : CALL dbcsr_reserve_blocks(outer_deriv_re, rows=rows(1:nblocks), cols=cols(1:nblocks))
1301 1054 : CALL dbcsr_reserve_blocks(outer_deriv_im, rows=rows(1:nblocks), cols=cols(1:nblocks))
1302 : END IF
1303 1591 : CALL dbcsr_finalize(outer_deriv_re)
1304 1591 : CALL dbcsr_finalize(outer_deriv_im)
1305 1591 : CALL dbcsr_set(outer_deriv_re, 0.0_dp)
1306 1591 : 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 1591 : col_blk_offset=col_blk_offset)
1310 1591 : CALL dbcsr_iterator_start(iter, outer_deriv_re)
1311 2647 : DO WHILE (dbcsr_iterator_blocks_left(iter))
1312 1056 : CALL dbcsr_iterator_next_block(iter, row, col)
1313 1056 : CALL dbcsr_get_block_p(inner_deriv_re, row, col, block_in_re, found_re)
1314 1056 : CALL dbcsr_get_block_p(inner_deriv_im, row, col, block_in_im, found_im)
1315 1056 : CALL dbcsr_get_block_p(outer_deriv_re, row, col, block_out_re, found_out_re)
1316 1056 : CALL dbcsr_get_block_p(outer_deriv_im, row, col, block_out_im, found_out_im)
1317 1056 : CPASSERT(found_out_re .AND. found_out_im)
1318 :
1319 11467 : DO i = 1, SIZE(block_out_re, 1)
1320 122108 : DO j = 1, SIZE(block_out_re, 2)
1321 112232 : e1 = evals(row_blk_offset(row) + i - 1)
1322 112232 : e2 = evals(col_blk_offset(col) + j - 1)
1323 112232 : re_part = 0.0_dp
1324 112232 : im_part = 0.0_dp
1325 112232 : IF (found_re) re_part = block_in_re(i, j)
1326 112232 : IF (found_im) im_part = block_in_im(i, j)
1327 112232 : cval_in = CMPLX(re_part, im_part, dp)
1328 112232 : kernel = qs_ot_complex_exp_frechet_kernel(e1, e2)
1329 112232 : IF (use_adjoint) kernel = CONJG(kernel)
1330 112232 : cval_in = cval_in*kernel
1331 112232 : block_out_re(i, j) = REAL(cval_in, KIND=dp)
1332 121052 : block_out_im(i, j) = AIMAG(cval_in)
1333 : END DO
1334 : END DO
1335 : END DO
1336 1591 : CALL dbcsr_iterator_stop(iter)
1337 12726 : 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 2108 : SUBROUTINE append_union_block(row_new, col_new)
1347 : INTEGER, INTENT(IN) :: row_new, col_new
1348 :
1349 : INTEGER :: iblock
1350 :
1351 2110 : DO iblock = 1, nblocks
1352 2110 : IF (rows(iblock) == row_new .AND. cols(iblock) == col_new) RETURN
1353 : END DO
1354 1056 : nblocks = nblocks + 1
1355 1056 : rows(nblocks) = row_new
1356 1056 : 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 9469 : 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 9469 : qs_ot_env%preconditioner => preconditioner
1377 9469 : qs_ot_env%os_valid = .FALSE.
1378 9469 : IF (.NOT. ASSOCIATED(qs_ot_env%matrix_psc0)) THEN
1379 9469 : CALL dbcsr_init_p(qs_ot_env%matrix_psc0)
1380 9469 : CALL dbcsr_copy(qs_ot_env%matrix_psc0, qs_ot_env%matrix_sc0, 'matrix_psc0')
1381 : END IF
1382 9469 : IF (qs_ot_env%has_complex_kpoint_state .AND. &
1383 : .NOT. ASSOCIATED(qs_ot_env%matrix_psc0_im)) THEN
1384 747 : CALL dbcsr_init_p(qs_ot_env%matrix_psc0_im)
1385 747 : CALL dbcsr_copy(qs_ot_env%matrix_psc0_im, qs_ot_env%matrix_sc0_im, 'matrix_psc0_im')
1386 : END IF
1387 :
1388 9469 : IF (.NOT. qs_ot_env%use_dx) THEN
1389 5761 : qs_ot_env%use_dx = .TRUE.
1390 5761 : CALL dbcsr_init_p(qs_ot_env%matrix_dx)
1391 5761 : CALL dbcsr_copy(qs_ot_env%matrix_dx, qs_ot_env%matrix_gx, 'matrix_dx')
1392 5761 : IF (qs_ot_env%has_complex_kpoint_state) THEN
1393 440 : CALL dbcsr_init_p(qs_ot_env%matrix_dx_im)
1394 440 : CALL dbcsr_copy(qs_ot_env%matrix_dx_im, qs_ot_env%matrix_gx_im, 'matrix_dx_im')
1395 : END IF
1396 5761 : IF (qs_ot_env%settings%do_rotation) THEN
1397 74 : CALL dbcsr_init_p(qs_ot_env%rot_mat_dx)
1398 74 : CALL dbcsr_copy(qs_ot_env%rot_mat_dx, qs_ot_env%rot_mat_gx, 'rot_mat_dx')
1399 74 : IF (qs_ot_env%has_complex_kpoint_state) THEN
1400 48 : CALL dbcsr_init_p(qs_ot_env%rot_mat_dx_im)
1401 48 : 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 5761 : IF (qs_ot_env%settings%do_ener) THEN
1405 40 : ncoef = SIZE(qs_ot_env%ener_gx)
1406 120 : ALLOCATE (qs_ot_env%ener_dx(ncoef))
1407 320 : qs_ot_env%ener_dx = 0.0_dp
1408 : END IF
1409 : END IF
1410 :
1411 9469 : 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 60036 : 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 41124 : SELECT CASE (op_a)
1433 : CASE ('N')
1434 41124 : db_op_a = 'N'
1435 41124 : sign_a = 1.0_dp
1436 : CASE ('C')
1437 18912 : db_op_a = 'T'
1438 18912 : sign_a = -1.0_dp
1439 : CASE DEFAULT
1440 60036 : CPABORT("Complex matrix product expects N or C for op_a")
1441 : END SELECT
1442 49596 : SELECT CASE (op_b)
1443 : CASE ('N')
1444 49596 : db_op_b = 'N'
1445 49596 : sign_b = 1.0_dp
1446 : CASE ('C')
1447 10440 : db_op_b = 'T'
1448 10440 : sign_b = -1.0_dp
1449 : CASE DEFAULT
1450 60036 : CPABORT("Complex matrix product expects N or C for op_b")
1451 : END SELECT
1452 :
1453 60036 : CALL dbcsr_multiply(db_op_a, db_op_b, 1.0_dp, a_re, b_re, 0.0_dp, c_re)
1454 60036 : CALL dbcsr_multiply(db_op_a, db_op_b, 1.0_dp, a_im, b_im, 0.0_dp, tmp)
1455 60036 : CALL dbcsr_add(c_re, tmp, alpha_scalar=1.0_dp, beta_scalar=-sign_a*sign_b)
1456 :
1457 60036 : CALL dbcsr_multiply(db_op_a, db_op_b, sign_b, a_re, b_im, 0.0_dp, c_im)
1458 60036 : CALL dbcsr_multiply(db_op_a, db_op_b, sign_a, a_im, b_re, 0.0_dp, tmp)
1459 60036 : CALL dbcsr_add(c_im, tmp, alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
1460 :
1461 60036 : 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 700 : 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 140 : 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 140 : CALL timeset(routineN, handle)
1490 : !
1491 : !
1492 140 : CALL dbcsr_get_info(C_NEW, nfullrows_total=n, nfullcols_total=k)
1493 : !
1494 : ! C = C*expm(-G)
1495 140 : GU => qs_ot_env%buf1_k_k_nosym ! a buffer
1496 140 : U => qs_ot_env%buf2_k_k_nosym ! a buffer
1497 140 : Gp1 => qs_ot_env%buf3_k_k_nosym ! a buffer
1498 140 : Gp2 => qs_ot_env%buf4_k_k_nosym ! a buffer
1499 140 : 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 140 : f2 = 0.0_dp
1505 140 : CALL dbcsr_copy(C, C_NEW)
1506 140 : CALL dbcsr_iterator_start(iter, C)
1507 322 : DO WHILE (dbcsr_iterator_blocks_left(iter))
1508 182 : CALL dbcsr_iterator_next_block(iter, row, col, block, row_size=row_size, col_size=col_size)
1509 1162 : DO p = 1, col_size ! p
1510 8694 : DO i = 1, row_size ! i
1511 7672 : tmp = SQRT(block(i, p)**2 + f2_eps)
1512 7672 : f2 = f2 + tmp
1513 8512 : block(i, p) = block(i, p)/tmp
1514 : END DO
1515 : END DO
1516 : END DO
1517 140 : CALL dbcsr_iterator_stop(iter)
1518 140 : CALL dbcsr_get_info(C, group=group)
1519 140 : CALL group%sum(f2)
1520 : !
1521 : !
1522 140 : CALL dbcsr_multiply('T', 'N', 1.0_dp, C, C_NEW, 0.0_dp, GU)
1523 : !
1524 : ! antisymetrize
1525 140 : CALL dbcsr_get_info(GU, distribution=dist)
1526 : CALL dbcsr_transposed(U, GU, shallow_data_copy=.FALSE., &
1527 : use_distribution=dist, &
1528 140 : transpose_distribution=.FALSE.)
1529 140 : CALL dbcsr_add(GU, U, alpha_scalar=-0.5_dp, beta_scalar=0.5_dp)
1530 : !-------------------------------------------------------------------
1531 : !
1532 140 : norm_fro = dbcsr_frobenius_norm(GU)
1533 140 : 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 140 : CALL dbcsr_scale(GU, -alpha)
1542 : !
1543 : ! compute unitary transform
1544 : ! zeroth and first order
1545 140 : expfactor = 1.0_dp
1546 140 : CALL dbcsr_copy(U, GU)
1547 140 : CALL dbcsr_scale(U, expfactor)
1548 140 : CALL dbcsr_add_on_diag(U, 1.0_dp)
1549 : ! other orders
1550 140 : CALL dbcsr_copy(Gp1, GU)
1551 750 : DO i = 2, taylor_order
1552 : ! new power of G
1553 750 : CALL dbcsr_multiply('N', 'N', 1.0_dp, GU, Gp1, 0.0_dp, Gp2)
1554 750 : CALL dbcsr_copy(Gp1, Gp2)
1555 : ! add to the taylor expansion so far
1556 750 : expfactor = expfactor/REAL(i, KIND=dp)
1557 750 : CALL dbcsr_add(U, Gp1, alpha_scalar=1.0_dp, beta_scalar=expfactor)
1558 750 : norm_fro = dbcsr_frobenius_norm(Gp1)
1559 : !write(*,*) 'Taylor expansion i=',i,' norm(X^i)/i!=',norm_fro*expfactor
1560 750 : IF (norm_fro*expfactor < 1.0E-10_dp) EXIT
1561 : END DO
1562 : !
1563 : ! rotate MOs
1564 140 : CALL dbcsr_multiply('N', 'N', 1.0_dp, C_NEW, U, 0.0_dp, C)
1565 140 : CALL dbcsr_copy(C_NEW, C)
1566 : !
1567 : ! rotate SC
1568 140 : CALL dbcsr_multiply('N', 'N', 1.0_dp, SC, U, 0.0_dp, C)
1569 140 : CALL dbcsr_copy(SC, C)
1570 : !
1571 : ! rotate D_i
1572 140 : CALL dbcsr_multiply('N', 'N', 1.0_dp, D, U, 0.0_dp, C)
1573 140 : CALL dbcsr_copy(D, C)
1574 : !
1575 : ! rotate G_i-1
1576 140 : IF (ASSOCIATED(G_OLD)) THEN
1577 140 : CALL dbcsr_multiply('N', 'N', 1.0_dp, G_OLD, U, 0.0_dp, C)
1578 140 : CALL dbcsr_copy(G_OLD, C)
1579 : END IF
1580 : !
1581 140 : CALL timestop(handle)
1582 140 : 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 1550 : 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 775 : CALL timeset(routineN, handle)
1605 : !
1606 775 : CALL dbcsr_get_info(C_NEW, nfullrows_total=n, nfullcols_total=k)
1607 : !
1608 : ! P = U'*U
1609 775 : 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 775 : transa="N", para_env=qs_ot_env%para_env, blacs_env=qs_ot_env%blacs_env)
1614 : !
1615 : ! Update SC if needed
1616 775 : IF (update) THEN
1617 : CALL cp_dbcsr_cholesky_restore(SC, k, P, C_TMP, op="SOLVE", pos="RIGHT", &
1618 429 : transa="N", para_env=qs_ot_env%para_env, blacs_env=qs_ot_env%blacs_env)
1619 429 : CALL dbcsr_copy(SC, C_TMP)
1620 : END IF
1621 : !
1622 775 : CALL timestop(handle)
1623 775 : 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 422 : 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 422 : CALL timeset(routineN, handle)
1648 : !
1649 422 : CALL dbcsr_get_info(C_NEW, nfullrows_total=n, nfullcols_total=k)
1650 : !
1651 422 : V => qs_ot_env%buf1_k_k_nosym ! a buffer
1652 422 : W => qs_ot_env%buf2_k_k_nosym ! a buffer
1653 1688 : ALLOCATE (eig(k), fun(k))
1654 : !
1655 422 : 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 2366 : DO i = 1, k
1659 1944 : IF (eig(i) <= 0.0_dp) THEN
1660 0 : CPABORT("P not positive definite")
1661 : END IF
1662 2366 : IF (eig(i) < 1.0E-8_dp) THEN
1663 0 : fun(i) = 0.0_dp
1664 : ELSE
1665 1944 : fun(i) = 1.0_dp/SQRT(eig(i))
1666 : END IF
1667 : END DO
1668 422 : CALL dbcsr_copy(W, V)
1669 422 : CALL dbcsr_scale_by_vector(V, alpha=fun, side='right')
1670 422 : CALL dbcsr_multiply('N', 'T', 1.0_dp, W, V, 0.0_dp, P)
1671 : !
1672 : ! Update C
1673 422 : CALL dbcsr_multiply('N', 'N', 1.0_dp, C_OLD, P, 0.0_dp, C_NEW)
1674 : !
1675 : ! Update SC if needed
1676 422 : IF (update) THEN
1677 272 : CALL dbcsr_multiply('N', 'N', 1.0_dp, SC, P, 0.0_dp, C_TMP)
1678 272 : CALL dbcsr_copy(SC, C_TMP)
1679 : END IF
1680 : !
1681 422 : DEALLOCATE (eig, fun)
1682 : !
1683 422 : CALL timestop(handle)
1684 422 : 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 8066 : 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 4033 : CALL timeset(routineN, handle)
1714 : !
1715 4033 : CALL dbcsr_get_info(C_NEW, nfullrows_total=n, nfullcols_total=k)
1716 : !
1717 4033 : BUF_NOSYM => qs_ot_env%buf1_k_k_nosym ! a buffer
1718 4033 : BUF1 => qs_ot_env%buf1_k_k_sym ! a buffer
1719 4033 : BUF2 => qs_ot_env%buf2_k_k_sym ! a buffer
1720 4033 : FY => qs_ot_env%buf3_k_k_sym ! a buffer
1721 4033 : FT => qs_ot_env%buf4_k_k_sym ! a buffer
1722 : !
1723 : ! initialize the norm (already computed in qs_ot_get_orbitals_ref)
1724 4033 : norm = norm_in
1725 : !
1726 : ! can we do a quick exit?
1727 4033 : quick_exit = .FALSE.
1728 4033 : IF (norm < qs_ot_env%settings%eps_irac_quick_exit) quick_exit = .TRUE.
1729 : !
1730 : ! lets refine
1731 4033 : rescale = 1.0_dp
1732 4467 : DO irefine = 1, qs_ot_env%settings%max_irac
1733 : !
1734 : ! rescaling
1735 4467 : IF (norm > 1.0_dp) THEN
1736 12 : CALL dbcsr_scale(P, 1.0_dp/norm)
1737 12 : 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 4467 : qs_ot_env%settings%eps_irac_filter_matrix)
1743 : !
1744 : ! collect the transformation
1745 4467 : IF (irefine == 1) THEN
1746 4033 : CALL dbcsr_copy(FT, FY, name='FT')
1747 : ELSE
1748 434 : CALL dbcsr_multiply('N', 'N', 1.0_dp, FT, FY, 0.0_dp, BUF1)
1749 434 : IF (qs_ot_env%settings%eps_irac_filter_matrix > 0.0_dp) THEN
1750 4 : occ_in = dbcsr_get_occupation(buf1)
1751 4 : CALL dbcsr_filter(buf1, qs_ot_env%settings%eps_irac_filter_matrix)
1752 4 : occ_out = dbcsr_get_occupation(buf1)
1753 : END IF
1754 434 : CALL dbcsr_copy(FT, BUF1, name='FT')
1755 : END IF
1756 : !
1757 : ! quick exit if possible
1758 4467 : IF (quick_exit) THEN
1759 : EXIT
1760 : END IF
1761 : !
1762 : ! P = FY^T * P * FY
1763 1881 : CALL dbcsr_multiply('N', 'N', 1.0_dp, P, FY, 0.0_dp, BUF_NOSYM)
1764 1881 : IF (qs_ot_env%settings%eps_irac_filter_matrix > 0.0_dp) THEN
1765 8 : occ_in = dbcsr_get_occupation(buf_nosym)
1766 8 : CALL dbcsr_filter(buf_nosym, qs_ot_env%settings%eps_irac_filter_matrix)
1767 8 : occ_out = dbcsr_get_occupation(buf_nosym)
1768 : END IF
1769 1881 : CALL dbcsr_multiply('N', 'N', 1.0_dp, FY, BUF_NOSYM, 0.0_dp, P)
1770 1881 : IF (qs_ot_env%settings%eps_irac_filter_matrix > 0.0_dp) THEN
1771 8 : occ_in = dbcsr_get_occupation(p)
1772 8 : CALL dbcsr_filter(p, qs_ot_env%settings%eps_irac_filter_matrix)
1773 8 : occ_out = dbcsr_get_occupation(p)
1774 : END IF
1775 : !
1776 : ! check ||P-1||_gct
1777 1881 : CALL dbcsr_add_on_diag(P, -1.0_dp)
1778 1881 : norm_fro = dbcsr_frobenius_norm(P)
1779 1881 : norm_gct = dbcsr_gershgorin_norm(P)
1780 1881 : CALL dbcsr_add_on_diag(P, 1.0_dp)
1781 1881 : norm = MIN(norm_gct, norm_fro)
1782 : !
1783 : ! printing
1784 : !
1785 : ! blows up
1786 1881 : 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 0 : "the forum https://www.cp2k.org/")
1791 : END IF
1792 : !
1793 : ! can we do a quick exit next step?
1794 1881 : IF (norm < qs_ot_env%settings%eps_irac_quick_exit) quick_exit = .TRUE.
1795 : !
1796 : ! are we done?
1797 4467 : IF (norm < qs_ot_env%settings%eps_irac) EXIT
1798 : !
1799 : END DO
1800 : !
1801 : ! C_NEW = C_NEW * FT * rescale
1802 4033 : CALL dbcsr_multiply('N', 'N', rescale, C_OLD, FT, 0.0_dp, C_NEW)
1803 4033 : IF (qs_ot_env%settings%eps_irac_filter_matrix > 0.0_dp) THEN
1804 4 : occ_in = dbcsr_get_occupation(c_new)
1805 4 : CALL dbcsr_filter(c_new, qs_ot_env%settings%eps_irac_filter_matrix)
1806 4 : occ_out = dbcsr_get_occupation(c_new)
1807 : END IF
1808 : !
1809 : ! update SC = SC * FY * rescale
1810 4033 : IF (update) THEN
1811 1647 : CALL dbcsr_multiply('N', 'N', rescale, SC, FT, 0.0_dp, C_TMP)
1812 1647 : IF (qs_ot_env%settings%eps_irac_filter_matrix > 0.0_dp) THEN
1813 4 : occ_in = dbcsr_get_occupation(c_tmp)
1814 4 : CALL dbcsr_filter(c_tmp, qs_ot_env%settings%eps_irac_filter_matrix)
1815 4 : occ_out = dbcsr_get_occupation(c_tmp)
1816 : END IF
1817 1647 : CALL dbcsr_copy(SC, C_TMP)
1818 : END IF
1819 : !
1820 4033 : CALL timestop(handle)
1821 4033 : END SUBROUTINE qs_ot_ref_poly
1822 :
1823 : ! **************************************************************************************************
1824 : !> \brief ...
1825 : !> \param qs_ot_env1 ...
1826 : !> \return ...
1827 : ! **************************************************************************************************
1828 8708 : FUNCTION qs_ot_ref_update(qs_ot_env1) RESULT(update)
1829 : !
1830 : TYPE(qs_ot_type) :: qs_ot_env1
1831 : LOGICAL :: update
1832 :
1833 8708 : update = .FALSE.
1834 6792 : SELECT CASE (qs_ot_env1%settings%ot_method)
1835 : CASE ("CG", "SD")
1836 6792 : SELECT CASE (qs_ot_env1%settings%line_search_method)
1837 : CASE ("2PNT")
1838 6792 : IF (qs_ot_env1%line_search_count == 2) update = .TRUE.
1839 : CASE DEFAULT
1840 6792 : CPABORT("NYI")
1841 : END SELECT
1842 : CASE ("DIIS")
1843 814 : update = .TRUE.
1844 : CASE ("BROY", "LBFG")
1845 : ! These minimizers retain positions or secants in one fixed REF chart.
1846 814 : update = .FALSE.
1847 : CASE DEFAULT
1848 8708 : CPABORT("NYI")
1849 : END SELECT
1850 8708 : 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 5230 : 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 5230 : ortho_irac = qs_ot_env1%settings%ortho_irac
1865 5230 : IF (norm_in < qs_ot_env1%settings%eps_irac_switch) ortho_irac = "POLY"
1866 5230 : 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 10460 : 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 5230 : CALL timeset(routineN, handle)
1895 :
1896 5230 : CALL dbcsr_get_info(matrix_c, nfullrows_total=n, nfullcols_total=k)
1897 : !
1898 5230 : C_NEW => matrix_c
1899 5230 : C_OLD => matrix_x ! need to be carefully updated for the gradient !
1900 5230 : SC => matrix_sx ! need to be carefully updated for the gradient !
1901 5230 : G_OLD => matrix_gx_old ! need to be carefully updated for localization !
1902 5230 : D => matrix_dx ! need to be carefully updated for localization !
1903 5230 : S => matrix_s
1904 :
1905 5230 : P => qs_ot_env%p_k_k_sym ! a buffer
1906 5230 : C_TMP => qs_ot_env%buf1_n_k ! a buffer
1907 : !
1908 : ! do we need to update C_OLD and SC?
1909 5230 : 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 5230 : on_the_fly_loc = qs_ot_env1%settings%on_the_fly_loc
1916 : !
1917 : ! compute SC = S*C
1918 5230 : IF (ASSOCIATED(S)) THEN
1919 5230 : CALL dbcsr_multiply('N', 'N', 1.0_dp, S, C_OLD, 0.0_dp, SC)
1920 5230 : IF (qs_ot_env1%settings%eps_irac_filter_matrix > 0.0_dp) THEN
1921 4 : occ_in = dbcsr_get_occupation(sc)
1922 4 : CALL dbcsr_filter(sc, qs_ot_env1%settings%eps_irac_filter_matrix)
1923 4 : occ_out = dbcsr_get_occupation(sc)
1924 : END IF
1925 : ELSE
1926 0 : CALL dbcsr_copy(SC, C_OLD)
1927 : END IF
1928 : !
1929 : ! compute P = C'*SC
1930 5230 : CALL dbcsr_multiply('T', 'N', 1.0_dp, C_OLD, SC, 0.0_dp, P)
1931 5230 : IF (qs_ot_env1%settings%eps_irac_filter_matrix > 0.0_dp) THEN
1932 4 : occ_in = dbcsr_get_occupation(p)
1933 4 : CALL dbcsr_filter(p, qs_ot_env1%settings%eps_irac_filter_matrix)
1934 4 : occ_out = dbcsr_get_occupation(p)
1935 : END IF
1936 : !
1937 : ! check ||P-1||_f and ||P-1||_gct
1938 5230 : CALL dbcsr_add_on_diag(P, -1.0_dp)
1939 5230 : norm_fro = dbcsr_frobenius_norm(P)
1940 5230 : norm_gct = dbcsr_gershgorin_norm(P)
1941 5230 : CALL dbcsr_add_on_diag(P, 1.0_dp)
1942 5230 : norm = MIN(norm_gct, norm_fro)
1943 5230 : CALL qs_ot_ref_decide(qs_ot_env1, norm, ortho_irac)
1944 : !
1945 : ! select the orthogonality method
1946 775 : SELECT CASE (ortho_irac)
1947 : CASE ("CHOL")
1948 775 : CALL qs_ot_ref_chol(qs_ot_env, C_OLD, C_TMP, C_NEW, P, SC, update)
1949 : CASE ("LWDN")
1950 422 : CALL qs_ot_ref_lwdn(qs_ot_env, C_OLD, C_TMP, C_NEW, P, SC, update)
1951 : CASE ("POLY")
1952 4033 : CALL qs_ot_ref_poly(qs_ot_env, C_OLD, C_TMP, C_NEW, P, SC, norm, update)
1953 : CASE DEFAULT
1954 5230 : CPABORT("Wrong argument")
1955 : END SELECT
1956 : !
1957 : ! We update the C_i+1 and localization
1958 5230 : IF (update) THEN
1959 2348 : IF (on_the_fly_loc) THEN
1960 140 : CALL qs_ot_on_the_fly_localize(qs_ot_env, C_NEW, SC, G_OLD, D)
1961 : END IF
1962 2348 : CALL dbcsr_copy(C_OLD, C_NEW)
1963 : END IF
1964 :
1965 5230 : IF (qs_ot_env%settings%do_rotation) THEN
1966 14 : 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 14 : 0.0_dp, C_TMP)
1969 14 : CALL dbcsr_copy(C_NEW, C_TMP)
1970 : END IF
1971 : !
1972 5230 : CALL timestop(handle)
1973 5230 : 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 3478 : 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 3478 : 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 3478 : CALL timeset(routineN, handle)
2003 :
2004 3478 : CPASSERT(qs_ot_env%has_complex_kpoint_state)
2005 3478 : CPASSERT(ASSOCIATED(matrix_s))
2006 3478 : CPASSERT(ASSOCIATED(matrix_s_im))
2007 :
2008 3478 : c_re => qs_ot_env%matrix_x
2009 3478 : c_im => qs_ot_env%matrix_x_im
2010 3478 : f_re => qs_ot_env%matrix_ref_inv_sqrt
2011 3478 : f_im => qs_ot_env%matrix_ref_inv_sqrt_im
2012 3478 : sc_re => qs_ot_env%matrix_sx
2013 3478 : sc_im => qs_ot_env%matrix_sx_im
2014 3478 : p_re => qs_ot_env%buf1_k_k_sym
2015 3478 : p_im => qs_ot_env%buf2_k_k_sym
2016 3478 : v_re => qs_ot_env%buf3_k_k_sym
2017 3478 : v_im => qs_ot_env%buf4_k_k_sym
2018 3478 : w_re => qs_ot_env%buf1_k_k_nosym
2019 3478 : w_im => qs_ot_env%buf2_k_k_nosym
2020 3478 : tmp_kk => qs_ot_env%buf3_k_k_nosym
2021 3478 : tmp_nk => qs_ot_env%buf1_n_k
2022 :
2023 3478 : CALL dbcsr_get_info(c_re, nfullrows_total=n, nfullcols_total=k)
2024 3478 : IF (PRESENT(qs_ot_env1)) THEN
2025 3279 : update = qs_ot_ref_update(qs_ot_env1)
2026 : ELSE
2027 199 : 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 3478 : CALL dbcsr_multiply('N', 'N', 1.0_dp, matrix_s, c_re, 0.0_dp, sc_re)
2032 3478 : CALL dbcsr_multiply('N', 'N', 1.0_dp, matrix_s_im, c_im, 0.0_dp, tmp_nk)
2033 3478 : CALL dbcsr_add(sc_re, tmp_nk, alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
2034 :
2035 3478 : CALL dbcsr_multiply('N', 'N', 1.0_dp, matrix_s, c_im, 0.0_dp, sc_im)
2036 3478 : CALL dbcsr_multiply('N', 'N', 1.0_dp, matrix_s_im, c_re, 0.0_dp, tmp_nk)
2037 3478 : 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 3478 : CALL dbcsr_multiply('T', 'N', 1.0_dp, c_re, sc_re, 0.0_dp, p_re)
2041 3478 : CALL dbcsr_multiply('T', 'N', 1.0_dp, c_im, sc_im, 0.0_dp, tmp_kk)
2042 3478 : CALL dbcsr_add(p_re, tmp_kk, alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
2043 :
2044 3478 : CALL dbcsr_multiply('T', 'N', 1.0_dp, c_re, sc_im, 0.0_dp, p_im)
2045 3478 : CALL dbcsr_multiply('T', 'N', 1.0_dp, c_im, sc_re, 0.0_dp, tmp_kk)
2046 3478 : CALL dbcsr_add(p_im, tmp_kk, alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
2047 :
2048 13912 : 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 3478 : blacs_env=qs_ot_env%blacs_env)
2053 28896 : DO i = 1, k
2054 25418 : IF (eigenvalues(i) <= EPSILON(1.0_dp)) THEN
2055 0 : CPABORT("Complex REF overlap is not positive definite")
2056 : END IF
2057 28896 : 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 3478 : CALL dbcsr_copy(w_re, v_re)
2062 3478 : CALL dbcsr_copy(w_im, v_im)
2063 3478 : CALL dbcsr_scale_by_vector(w_re, alpha=inverse_sqrt, side='right')
2064 3478 : CALL dbcsr_scale_by_vector(w_im, alpha=inverse_sqrt, side='right')
2065 :
2066 3478 : CALL dbcsr_multiply('N', 'T', 1.0_dp, w_re, v_re, 0.0_dp, p_re)
2067 3478 : CALL dbcsr_multiply('N', 'T', 1.0_dp, w_im, v_im, 0.0_dp, tmp_kk)
2068 3478 : CALL dbcsr_add(p_re, tmp_kk, alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
2069 :
2070 3478 : CALL dbcsr_multiply('N', 'T', 1.0_dp, w_im, v_re, 0.0_dp, p_im)
2071 3478 : CALL dbcsr_multiply('N', 'T', 1.0_dp, w_re, v_im, 0.0_dp, tmp_kk)
2072 3478 : CALL dbcsr_add(p_im, tmp_kk, alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
2073 :
2074 3478 : CALL dbcsr_copy(f_re, p_re)
2075 3478 : CALL dbcsr_copy(f_im, p_im)
2076 :
2077 : ! Return the physical, orthonormal orbitals without changing a rejected
2078 : ! line-search coordinate.
2079 3478 : CALL dbcsr_multiply('N', 'N', 1.0_dp, c_re, p_re, 0.0_dp, matrix_c)
2080 3478 : CALL dbcsr_multiply('N', 'N', 1.0_dp, c_im, p_im, 0.0_dp, tmp_nk)
2081 3478 : CALL dbcsr_add(matrix_c, tmp_nk, alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
2082 :
2083 3478 : CALL dbcsr_multiply('N', 'N', 1.0_dp, c_re, p_im, 0.0_dp, matrix_c_im)
2084 3478 : CALL dbcsr_multiply('N', 'N', 1.0_dp, c_im, p_re, 0.0_dp, tmp_nk)
2085 3478 : CALL dbcsr_add(matrix_c_im, tmp_nk, alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
2086 :
2087 3478 : IF (update) THEN
2088 1372 : CALL dbcsr_copy(c_re, matrix_c)
2089 1372 : CALL dbcsr_copy(c_im, matrix_c_im)
2090 :
2091 1372 : CALL dbcsr_set(f_re, 0.0_dp)
2092 1372 : CALL dbcsr_add_on_diag(f_re, alpha=1.0_dp)
2093 1372 : CALL dbcsr_set(f_im, 0.0_dp)
2094 :
2095 1372 : CALL dbcsr_multiply('N', 'N', 1.0_dp, matrix_s, c_re, 0.0_dp, sc_re)
2096 1372 : CALL dbcsr_multiply('N', 'N', 1.0_dp, matrix_s_im, c_im, 0.0_dp, tmp_nk)
2097 1372 : CALL dbcsr_add(sc_re, tmp_nk, alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
2098 :
2099 1372 : CALL dbcsr_multiply('N', 'N', 1.0_dp, matrix_s, c_im, 0.0_dp, sc_im)
2100 1372 : CALL dbcsr_multiply('N', 'N', 1.0_dp, matrix_s_im, c_re, 0.0_dp, tmp_nk)
2101 1372 : CALL dbcsr_add(sc_im, tmp_nk, alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
2102 : END IF
2103 :
2104 3478 : IF (qs_ot_env%settings%do_rotation) THEN
2105 2184 : CALL qs_ot_generate_rotation_complex(qs_ot_env)
2106 :
2107 2184 : CALL dbcsr_copy(rotated_re, matrix_c, name="rotated_re")
2108 2184 : CALL dbcsr_copy(rotated_im, matrix_c_im, name="rotated_im")
2109 2184 : 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 2184 : 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 2184 : 0.0_dp, rotation_tmp)
2116 2184 : 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 2184 : 0.0_dp, rotated_im)
2120 : CALL dbcsr_multiply('N', 'N', 1.0_dp, matrix_c_im, qs_ot_env%rot_mat_u, &
2121 2184 : 0.0_dp, rotation_tmp)
2122 2184 : CALL dbcsr_add(rotated_im, rotation_tmp, alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
2123 :
2124 2184 : CALL dbcsr_copy(matrix_c, rotated_re)
2125 2184 : CALL dbcsr_copy(matrix_c_im, rotated_im)
2126 2184 : CALL dbcsr_release(rotated_re)
2127 2184 : CALL dbcsr_release(rotated_im)
2128 2184 : CALL dbcsr_release(rotation_tmp)
2129 : END IF
2130 :
2131 3478 : DEALLOCATE (eigenvalues, inverse_sqrt)
2132 :
2133 3478 : CALL timestop(handle)
2134 6956 : 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 8934 : 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 4467 : CALL timeset(routineN, handle)
2156 :
2157 4467 : CALL dbcsr_get_info(P, nfullcols_total=k)
2158 4467 : SELECT CASE (irac_degree)
2159 : CASE (2)
2160 : ! C_out = C_in * ( 15/8 * I - 10/8 * P + 3/8 * P^2)
2161 0 : r = 3.0_dp/8.0_dp
2162 0 : CALL dbcsr_multiply('N', 'N', r, P, P, 0.0_dp, FY)
2163 0 : IF (eps_irac_filter_matrix > 0.0_dp) THEN
2164 0 : occ_in = dbcsr_get_occupation(fy)
2165 0 : CALL dbcsr_filter(fy, eps_irac_filter_matrix)
2166 0 : occ_out = dbcsr_get_occupation(fy)
2167 : END IF
2168 0 : r = -10.0_dp/8.0_dp
2169 0 : CALL dbcsr_add(FY, P, alpha_scalar=1.0_dp, beta_scalar=r)
2170 0 : r = 15.0_dp/8.0_dp
2171 0 : 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 0 : CALL dbcsr_multiply('N', 'N', 1.0_dp, P, P, 0.0_dp, P2)
2175 0 : IF (eps_irac_filter_matrix > 0.0_dp) THEN
2176 0 : occ_in = dbcsr_get_occupation(p2)
2177 0 : CALL dbcsr_filter(p2, eps_irac_filter_matrix)
2178 0 : occ_out = dbcsr_get_occupation(p2)
2179 : END IF
2180 0 : r = -5.0_dp/16.0_dp
2181 0 : CALL dbcsr_multiply('N', 'N', r, P2, P, 0.0_dp, FY)
2182 0 : IF (eps_irac_filter_matrix > 0.0_dp) THEN
2183 0 : occ_in = dbcsr_get_occupation(fy)
2184 0 : CALL dbcsr_filter(fy, eps_irac_filter_matrix)
2185 0 : occ_out = dbcsr_get_occupation(fy)
2186 : END IF
2187 0 : r = 21.0_dp/16.0_dp
2188 0 : CALL dbcsr_add(FY, P2, alpha_scalar=1.0_dp, beta_scalar=r)
2189 0 : r = -35.0_dp/16.0_dp
2190 0 : CALL dbcsr_add(FY, P, alpha_scalar=1.0_dp, beta_scalar=r)
2191 0 : r = 35.0_dp/16.0_dp
2192 0 : 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 4467 : CALL dbcsr_multiply('N', 'N', 1.0_dp, P, P, 0.0_dp, P2) ! P^2
2197 4467 : IF (eps_irac_filter_matrix > 0.0_dp) THEN
2198 8 : occ_in = dbcsr_get_occupation(p2)
2199 8 : CALL dbcsr_filter(p2, eps_irac_filter_matrix)
2200 8 : occ_out = dbcsr_get_occupation(p2)
2201 : END IF
2202 4467 : r = -180.0_dp/128.0_dp
2203 4467 : CALL dbcsr_add(T, P, alpha_scalar=0.0_dp, beta_scalar=r) ! T=-180/128*P
2204 4467 : r = 35.0_dp/128.0_dp
2205 4467 : CALL dbcsr_add(T, P2, alpha_scalar=1.0_dp, beta_scalar=r) ! T=T+35/128*P^2
2206 4467 : CALL dbcsr_multiply('N', 'N', 1.0_dp, T, P2, 0.0_dp, FY) ! Y=T*P^2
2207 4467 : IF (eps_irac_filter_matrix > 0.0_dp) THEN
2208 8 : occ_in = dbcsr_get_occupation(fy)
2209 8 : CALL dbcsr_filter(fy, eps_irac_filter_matrix)
2210 8 : occ_out = dbcsr_get_occupation(fy)
2211 : END IF
2212 4467 : r = 378.0_dp/128.0_dp
2213 4467 : CALL dbcsr_add(FY, P2, alpha_scalar=1.0_dp, beta_scalar=r) ! Y=Y+378/128*P^2
2214 4467 : r = -420.0_dp/128.0_dp
2215 4467 : CALL dbcsr_add(FY, P, alpha_scalar=1.0_dp, beta_scalar=r) ! Y=Y-420/128*P
2216 4467 : r = 315.0_dp/128.0_dp
2217 4467 : CALL dbcsr_add_on_diag(FY, alpha=r) ! Y=Y+315/128*I
2218 : CASE DEFAULT
2219 4467 : CPABORT("This irac_order NYI")
2220 : END SELECT
2221 4467 : CALL timestop(handle)
2222 4467 : 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 6628 : 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 3314 : CALL timeset(routineN, handle)
2244 :
2245 3314 : CALL dbcsr_get_info(matrix_x, nfullrows_total=n, nfullcols_total=k)
2246 : !
2247 3314 : C => matrix_x ! NBsf*NOcc
2248 3314 : SC => matrix_sx ! NBsf*NOcc need to be up2date
2249 3314 : HC => matrix_hc ! NBsf*NOcc
2250 3314 : G => matrix_gx ! NBsf*NOcc
2251 3314 : CHC => qs_ot_env%buf1_k_k_sym ! buffer
2252 :
2253 3314 : 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 8 : CALL qs_ot_rot_mat_derivative(qs_ot_env)
2257 8 : 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 8 : 0.0_dp, HC_work)
2260 : ELSE
2261 : HC_work => HC
2262 : END IF
2263 :
2264 : ! C'*(H*C)
2265 3314 : CALL dbcsr_multiply('T', 'N', 1.0_dp, C, HC_work, 0.0_dp, CHC)
2266 3314 : IF (qs_ot_env%settings%eps_irac_filter_matrix > 0.0_dp) THEN
2267 4 : occ_in = dbcsr_get_occupation(chc)
2268 4 : CALL dbcsr_filter(chc, qs_ot_env%settings%eps_irac_filter_matrix)
2269 4 : occ_out = dbcsr_get_occupation(chc)
2270 : END IF
2271 : ! (S*C)*(C'*H*C)
2272 3314 : CALL dbcsr_multiply('N', 'N', 1.0_dp, SC, CHC, 0.0_dp, G)
2273 3314 : IF (qs_ot_env%settings%eps_irac_filter_matrix > 0.0_dp) THEN
2274 4 : occ_in = dbcsr_get_occupation(g)
2275 4 : CALL dbcsr_filter(g, qs_ot_env%settings%eps_irac_filter_matrix)
2276 4 : occ_out = dbcsr_get_occupation(g)
2277 : END IF
2278 : ! G = 2*(1-S*C*C')*H*C
2279 3314 : CALL dbcsr_add(G, HC_work, alpha_scalar=-1.0_dp, beta_scalar=1.0_dp)
2280 : !
2281 3314 : CALL timestop(handle)
2282 3314 : 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 4430 : 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 2215 : CALL timeset(routineN, handle)
2308 :
2309 2215 : CPASSERT(qs_ot_env%has_complex_kpoint_state)
2310 2215 : CPASSERT(ASSOCIATED(matrix_hc))
2311 2215 : CPASSERT(ASSOCIATED(matrix_hc_im))
2312 :
2313 2215 : f_re => qs_ot_env%matrix_ref_inv_sqrt
2314 2215 : f_im => qs_ot_env%matrix_ref_inv_sqrt_im
2315 2215 : sc_re => qs_ot_env%matrix_sx
2316 2215 : sc_im => qs_ot_env%matrix_sx_im
2317 2215 : hc_re => matrix_hc
2318 2215 : hc_im => matrix_hc_im
2319 2215 : hc_rotation_re => hc_re
2320 2215 : hc_rotation_im => hc_im
2321 2215 : IF (PRESENT(matrix_hc_rotation) .OR. PRESENT(matrix_hc_rotation_im)) THEN
2322 408 : CPASSERT(PRESENT(matrix_hc_rotation) .AND. PRESENT(matrix_hc_rotation_im))
2323 408 : CPASSERT(ASSOCIATED(matrix_hc_rotation))
2324 408 : CPASSERT(ASSOCIATED(matrix_hc_rotation_im))
2325 408 : hc_rotation_re => matrix_hc_rotation
2326 408 : hc_rotation_im => matrix_hc_rotation_im
2327 : END IF
2328 2215 : hc_work_re => hc_re
2329 2215 : hc_work_im => hc_im
2330 2215 : g_re => qs_ot_env%matrix_gx
2331 2215 : g_im => qs_ot_env%matrix_gx_im
2332 2215 : b_re => qs_ot_env%buf1_k_k_sym
2333 2215 : b_im => qs_ot_env%buf2_k_k_sym
2334 2215 : tmp_kk => qs_ot_env%buf3_k_k_sym
2335 2215 : q_re => qs_ot_env%buf1_n_k
2336 2215 : q_im => qs_ot_env%buf1_n_k_dp
2337 :
2338 2215 : 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 2215 : CALL dbcsr_multiply('N', 'N', 1.0_dp, qs_ot_env%matrix_x, f_re, 0.0_dp, q_re)
2342 2215 : CALL dbcsr_multiply('N', 'N', 1.0_dp, qs_ot_env%matrix_x_im, f_im, 0.0_dp, g_re)
2343 2215 : CALL dbcsr_add(q_re, g_re, alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
2344 :
2345 2215 : CALL dbcsr_multiply('N', 'N', 1.0_dp, qs_ot_env%matrix_x, f_im, 0.0_dp, q_im)
2346 2215 : CALL dbcsr_multiply('N', 'N', 1.0_dp, qs_ot_env%matrix_x_im, f_re, 0.0_dp, g_re)
2347 2215 : CALL dbcsr_add(q_im, g_re, alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
2348 :
2349 2215 : IF (qs_ot_env%settings%do_rotation) THEN
2350 1415 : 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 1415 : 0.0_dp, qs_ot_env%rot_mat_dedu)
2355 1415 : 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 1415 : 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 1415 : 0.0_dp, qs_ot_env%rot_mat_dedu_im)
2361 1415 : 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 1415 : alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
2364 1415 : CALL qs_ot_rot_mat_derivative_complex(qs_ot_env)
2365 :
2366 : ! The REF coordinate sees G_Q=G_C*U^H.
2367 1415 : CALL dbcsr_copy(hc_rot_re, hc_re, name="hc_rot_re")
2368 1415 : CALL dbcsr_copy(hc_rot_im, hc_im, name="hc_rot_im")
2369 1415 : 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 1415 : 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 1415 : 0.0_dp, tmp_nk)
2374 1415 : 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 1415 : 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 1415 : 0.0_dp, tmp_nk)
2380 1415 : CALL dbcsr_add(hc_rot_im, tmp_nk, alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
2381 1415 : hc_work_re => hc_rot_re
2382 1415 : hc_work_im => hc_rot_im
2383 : END IF
2384 :
2385 : ! B = Q^H*G_Q. For uniform fixed occupations B is Hermitian.
2386 2215 : CALL dbcsr_multiply('T', 'N', 1.0_dp, q_re, hc_work_re, 0.0_dp, b_re)
2387 2215 : CALL dbcsr_multiply('T', 'N', 1.0_dp, q_im, hc_work_im, 0.0_dp, tmp_kk)
2388 2215 : CALL dbcsr_add(b_re, tmp_kk, alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
2389 :
2390 2215 : CALL dbcsr_multiply('T', 'N', 1.0_dp, q_re, hc_work_im, 0.0_dp, b_im)
2391 2215 : CALL dbcsr_multiply('T', 'N', 1.0_dp, q_im, hc_work_re, 0.0_dp, tmp_kk)
2392 2215 : CALL dbcsr_add(b_im, tmp_kk, alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
2393 :
2394 2215 : IF (qs_ot_env%settings%eps_irac_filter_matrix > 0.0_dp) THEN
2395 0 : occ_in = dbcsr_get_occupation(b_re)
2396 0 : CALL dbcsr_filter(b_re, qs_ot_env%settings%eps_irac_filter_matrix)
2397 0 : occ_out = dbcsr_get_occupation(b_re)
2398 0 : occ_in = dbcsr_get_occupation(b_im)
2399 0 : CALL dbcsr_filter(b_im, qs_ot_env%settings%eps_irac_filter_matrix)
2400 0 : 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 2215 : CALL dbcsr_multiply('N', 'N', 1.0_dp, sc_re, f_re, 0.0_dp, g_re)
2405 2215 : CALL dbcsr_multiply('N', 'N', 1.0_dp, sc_im, f_im, 0.0_dp, q_re)
2406 2215 : CALL dbcsr_add(g_re, q_re, alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
2407 :
2408 2215 : CALL dbcsr_multiply('N', 'N', 1.0_dp, sc_re, f_im, 0.0_dp, g_im)
2409 2215 : CALL dbcsr_multiply('N', 'N', 1.0_dp, sc_im, f_re, 0.0_dp, q_re)
2410 2215 : 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 2215 : CALL dbcsr_multiply('N', 'N', 1.0_dp, g_re, b_re, 0.0_dp, q_re)
2414 2215 : CALL dbcsr_multiply('N', 'N', 1.0_dp, g_im, b_im, 0.0_dp, q_im)
2415 2215 : CALL dbcsr_add(q_re, q_im, alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
2416 :
2417 2215 : CALL dbcsr_multiply('N', 'N', 1.0_dp, g_re, b_im, 0.0_dp, q_im)
2418 2215 : CALL dbcsr_multiply('N', 'N', 1.0_dp, g_im, b_re, 0.0_dp, g_re)
2419 2215 : CALL dbcsr_add(q_im, g_re, alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
2420 :
2421 2215 : CALL dbcsr_add(q_re, hc_work_re, alpha_scalar=-1.0_dp, beta_scalar=1.0_dp)
2422 2215 : 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 2215 : CALL dbcsr_multiply('N', 'N', 1.0_dp, q_re, f_re, 0.0_dp, g_re)
2426 2215 : CALL dbcsr_multiply('N', 'N', 1.0_dp, q_im, f_im, 0.0_dp, g_im)
2427 2215 : CALL dbcsr_add(g_re, g_im, alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
2428 :
2429 2215 : CALL dbcsr_multiply('N', 'N', 1.0_dp, q_re, f_im, 0.0_dp, g_im)
2430 2215 : CALL dbcsr_multiply('N', 'N', 1.0_dp, q_im, f_re, 0.0_dp, q_re)
2431 2215 : CALL dbcsr_add(g_im, q_re, alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
2432 :
2433 2215 : 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 1415 : sc_re, sc_im, g_re, g_im, qs_ot_env)
2436 : END IF
2437 :
2438 2215 : IF (qs_ot_env%settings%eps_irac_filter_matrix > 0.0_dp) THEN
2439 0 : occ_in = dbcsr_get_occupation(g_re)
2440 0 : CALL dbcsr_filter(g_re, qs_ot_env%settings%eps_irac_filter_matrix)
2441 0 : occ_out = dbcsr_get_occupation(g_re)
2442 0 : occ_in = dbcsr_get_occupation(g_im)
2443 0 : CALL dbcsr_filter(g_im, qs_ot_env%settings%eps_irac_filter_matrix)
2444 0 : occ_out = dbcsr_get_occupation(g_im)
2445 : END IF
2446 :
2447 2215 : IF (qs_ot_env%settings%do_rotation) THEN
2448 1415 : CALL dbcsr_release(hc_rot_re)
2449 1415 : CALL dbcsr_release(hc_rot_im)
2450 1415 : CALL dbcsr_release(tmp_nk)
2451 : END IF
2452 :
2453 2215 : CALL timestop(handle)
2454 :
2455 2215 : 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 8528 : 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 8528 : CALL dbcsr_copy(transposed, matrix, name='square_transposed')
2469 8528 : CALL dbcsr_copy(identity, identity_template, name='transpose_identity')
2470 8528 : CALL dbcsr_set(identity, 0.0_dp)
2471 8528 : CALL dbcsr_add_on_diag(identity, 1.0_dp)
2472 8528 : CALL dbcsr_multiply('T', 'N', 1.0_dp, matrix, identity, 0.0_dp, transposed)
2473 8528 : CALL dbcsr_release(identity)
2474 :
2475 8528 : 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 1415 : 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 1415 : INTEGER, DIMENSION(:), POINTER :: col_blk_offset, row_blk_offset
2500 : LOGICAL :: found_im, found_re
2501 : REAL(KIND=dp) :: denominator
2502 1415 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: f_evals
2503 1415 : 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 1415 : CALL timeset(routineN, handle)
2509 :
2510 1415 : CALL dbcsr_get_info(a_re, nfullrows_total=k)
2511 1415 : IF (k == 0) THEN
2512 0 : CALL timestop(handle)
2513 0 : 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 1415 : CALL qs_ot_square_transpose(a_re, tmp_kk, qs_ot_env%rot_mat_u)
2519 1415 : CALL dbcsr_copy(anti_re, a_re, name='anti_re')
2520 1415 : CALL dbcsr_add(anti_re, tmp_kk, alpha_scalar=0.5_dp, beta_scalar=-0.5_dp)
2521 1415 : CALL qs_ot_square_transpose(a_im, tmp_kk, qs_ot_env%rot_mat_u)
2522 1415 : CALL dbcsr_copy(anti_im, a_im, name='anti_im')
2523 1415 : 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 1415 : CALL dbcsr_copy(f_work_re, f_re, name='f_work_re')
2528 1415 : CALL dbcsr_copy(f_work_im, f_im, name='f_work_im')
2529 1415 : CALL dbcsr_copy(v_re, f_re, name='v_re')
2530 1415 : CALL dbcsr_copy(v_im, f_im, name='v_im')
2531 4245 : 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 1415 : para_env=qs_ot_env%para_env, blacs_env=qs_ot_env%blacs_env)
2535 10605 : IF (MINVAL(f_evals) <= EPSILON(1.0_dp)) THEN
2536 0 : CPABORT('Complex REF inverse square root is not positive definite')
2537 : END IF
2538 :
2539 : ! inner = V^H antiherm(A) V.
2540 1415 : CALL dbcsr_copy(work_re, anti_re, name='work_re')
2541 1415 : CALL dbcsr_copy(work_im, anti_im, name='work_im')
2542 1415 : CALL dbcsr_multiply('N', 'N', 1.0_dp, anti_re, v_re, 0.0_dp, work_re)
2543 1415 : CALL dbcsr_multiply('N', 'N', 1.0_dp, anti_im, v_im, 0.0_dp, tmp_kk)
2544 1415 : CALL dbcsr_add(work_re, tmp_kk, alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
2545 1415 : CALL dbcsr_multiply('N', 'N', 1.0_dp, anti_re, v_im, 0.0_dp, work_im)
2546 1415 : CALL dbcsr_multiply('N', 'N', 1.0_dp, anti_im, v_re, 0.0_dp, tmp_kk)
2547 1415 : CALL dbcsr_add(work_im, tmp_kk, alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
2548 :
2549 1415 : CALL dbcsr_copy(inner_re, anti_re, name='inner_re')
2550 1415 : CALL dbcsr_copy(inner_im, anti_im, name='inner_im')
2551 1415 : CALL dbcsr_multiply('T', 'N', 1.0_dp, v_re, work_re, 0.0_dp, inner_re)
2552 1415 : CALL dbcsr_multiply('T', 'N', 1.0_dp, v_im, work_im, 1.0_dp, inner_re)
2553 1415 : CALL dbcsr_multiply('T', 'N', 1.0_dp, v_re, work_im, 0.0_dp, inner_im)
2554 1415 : 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 1415 : col_blk_offset=col_blk_offset)
2558 1415 : CALL dbcsr_iterator_start(iter, inner_re)
2559 2384 : DO WHILE (dbcsr_iterator_blocks_left(iter))
2560 969 : CALL dbcsr_iterator_next_block(iter, row, col)
2561 969 : CALL dbcsr_get_block_p(inner_re, row, col, block_re, found_re)
2562 969 : CALL dbcsr_get_block_p(inner_im, row, col, block_im, found_im)
2563 969 : CPASSERT(found_re .AND. found_im)
2564 9186 : DO i = 1, SIZE(block_re, 1)
2565 78339 : DO j = 1, SIZE(block_re, 2)
2566 : denominator = 1.0_dp/f_evals(row_blk_offset(row) + i - 1) + &
2567 70568 : 1.0_dp/f_evals(col_blk_offset(col) + j - 1)
2568 70568 : block_re(i, j) = block_re(i, j)/denominator
2569 77370 : block_im(i, j) = block_im(i, j)/denominator
2570 : END DO
2571 : END DO
2572 : END DO
2573 1415 : CALL dbcsr_iterator_stop(iter)
2574 :
2575 : ! Z = V inner V^H.
2576 1415 : CALL dbcsr_multiply('N', 'N', 1.0_dp, v_re, inner_re, 0.0_dp, work_re)
2577 1415 : CALL dbcsr_multiply('N', 'N', 1.0_dp, v_im, inner_im, 0.0_dp, tmp_kk)
2578 1415 : CALL dbcsr_add(work_re, tmp_kk, alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
2579 1415 : CALL dbcsr_multiply('N', 'N', 1.0_dp, v_re, inner_im, 0.0_dp, work_im)
2580 1415 : CALL dbcsr_multiply('N', 'N', 1.0_dp, v_im, inner_re, 0.0_dp, tmp_kk)
2581 1415 : CALL dbcsr_add(work_im, tmp_kk, alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
2582 :
2583 1415 : CALL dbcsr_copy(z_re, anti_re, name='z_re')
2584 1415 : CALL dbcsr_copy(z_im, anti_im, name='z_im')
2585 1415 : CALL dbcsr_multiply('N', 'T', 1.0_dp, work_re, v_re, 0.0_dp, z_re)
2586 1415 : CALL dbcsr_multiply('N', 'T', 1.0_dp, work_im, v_im, 1.0_dp, z_re)
2587 1415 : CALL dbcsr_multiply('N', 'T', 1.0_dp, work_im, v_re, 0.0_dp, z_im)
2588 1415 : 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 1415 : CALL dbcsr_copy(sq_re, sx_re, name='sq_re')
2592 1415 : CALL dbcsr_copy(sq_im, sx_im, name='sq_im')
2593 1415 : CALL dbcsr_copy(tmp_nk, sx_re, name='tmp_nk')
2594 1415 : CALL dbcsr_multiply('N', 'N', 1.0_dp, sx_re, f_re, 0.0_dp, sq_re)
2595 1415 : CALL dbcsr_multiply('N', 'N', 1.0_dp, sx_im, f_im, 0.0_dp, tmp_nk)
2596 1415 : CALL dbcsr_add(sq_re, tmp_nk, alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
2597 1415 : CALL dbcsr_multiply('N', 'N', 1.0_dp, sx_re, f_im, 0.0_dp, sq_im)
2598 1415 : CALL dbcsr_multiply('N', 'N', 1.0_dp, sx_im, f_re, 0.0_dp, tmp_nk)
2599 1415 : CALL dbcsr_add(sq_im, tmp_nk, alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
2600 :
2601 1415 : CALL dbcsr_copy(vertical_re, sx_re, name='vertical_re')
2602 1415 : CALL dbcsr_copy(vertical_im, sx_im, name='vertical_im')
2603 1415 : CALL dbcsr_multiply('N', 'N', 1.0_dp, sq_re, z_re, 0.0_dp, vertical_re)
2604 1415 : CALL dbcsr_multiply('N', 'N', 1.0_dp, sq_im, z_im, 0.0_dp, tmp_nk)
2605 1415 : CALL dbcsr_add(vertical_re, tmp_nk, alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
2606 1415 : CALL dbcsr_multiply('N', 'N', 1.0_dp, sq_re, z_im, 0.0_dp, vertical_im)
2607 1415 : CALL dbcsr_multiply('N', 'N', 1.0_dp, sq_im, z_re, 0.0_dp, tmp_nk)
2608 1415 : CALL dbcsr_add(vertical_im, tmp_nk, alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
2609 1415 : CALL dbcsr_add(gradient_re, vertical_re, alpha_scalar=1.0_dp, beta_scalar=2.0_dp)
2610 1415 : CALL dbcsr_add(gradient_im, vertical_im, alpha_scalar=1.0_dp, beta_scalar=2.0_dp)
2611 :
2612 1415 : DEALLOCATE (f_evals)
2613 1415 : CALL dbcsr_release(anti_im)
2614 1415 : CALL dbcsr_release(anti_re)
2615 1415 : CALL dbcsr_release(f_work_im)
2616 1415 : CALL dbcsr_release(f_work_re)
2617 1415 : CALL dbcsr_release(inner_im)
2618 1415 : CALL dbcsr_release(inner_re)
2619 1415 : CALL dbcsr_release(sq_im)
2620 1415 : CALL dbcsr_release(sq_re)
2621 1415 : CALL dbcsr_release(tmp_kk)
2622 1415 : CALL dbcsr_release(tmp_nk)
2623 1415 : CALL dbcsr_release(v_im)
2624 1415 : CALL dbcsr_release(v_re)
2625 1415 : CALL dbcsr_release(vertical_im)
2626 1415 : CALL dbcsr_release(vertical_re)
2627 1415 : CALL dbcsr_release(work_im)
2628 1415 : CALL dbcsr_release(work_re)
2629 1415 : CALL dbcsr_release(z_im)
2630 1415 : CALL dbcsr_release(z_re)
2631 :
2632 1415 : CALL timestop(handle)
2633 :
2634 2830 : 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 333045 : 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 111015 : CALL timeset(routineN, handle)
2655 :
2656 111015 : 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 111015 : 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 111015 : max_iter = 30; threshold = 1.0E-03_dp
2665 111015 : CALL arnoldi_extremal(qs_ot_env%matrix_p, max_ev, min_ev, converged, threshold, max_iter)
2666 111015 : qs_ot_env%largest_eval_upper_bound = MAX(max_ev, ABS(min_ev))
2667 :
2668 111015 : IF (.NOT. converged) qs_ot_env%largest_eval_upper_bound = dbcsr_gershgorin_norm(qs_ot_env%matrix_p)
2669 111015 : CALL decide_strategy(qs_ot_env)
2670 111015 : IF (qs_ot_env%do_taylor) THEN
2671 59916 : CALL qs_ot_p2m_taylor(qs_ot_env)
2672 : ELSE
2673 51099 : CALL qs_ot_p2m_diag(qs_ot_env)
2674 : END IF
2675 :
2676 111015 : IF (qs_ot_env%settings%do_rotation) THEN
2677 3142 : CALL qs_ot_generate_rotation(qs_ot_env)
2678 : END IF
2679 :
2680 111015 : CALL timestop(handle)
2681 :
2682 111015 : 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 4931 : 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 4931 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: exp_evals_im, exp_evals_re
2698 : TYPE(dbcsr_type) :: h_re, tmp, w_im, w_re
2699 :
2700 4931 : CALL timeset(routineN, handle)
2701 :
2702 4931 : CPASSERT(qs_ot_env%has_complex_kpoint_state)
2703 4931 : CPASSERT(ASSOCIATED(qs_ot_env%rot_mat_x_im))
2704 4931 : CPASSERT(ASSOCIATED(qs_ot_env%rot_mat_u_im))
2705 :
2706 4931 : CALL dbcsr_get_info(qs_ot_env%rot_mat_x, nfullrows_total=k)
2707 4931 : IF (k /= 0) THEN
2708 : rot_norm = dbcsr_frobenius_norm(qs_ot_env%rot_mat_x) + &
2709 4931 : dbcsr_frobenius_norm(qs_ot_env%rot_mat_x_im)
2710 4931 : IF (rot_norm <= EPSILON(1.0_dp)) THEN
2711 990 : CALL dbcsr_set(qs_ot_env%rot_mat_u, 0.0_dp)
2712 990 : CALL dbcsr_add_on_diag(qs_ot_env%rot_mat_u, 1.0_dp)
2713 990 : CALL dbcsr_set(qs_ot_env%rot_mat_u_im, 0.0_dp)
2714 990 : 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 3941 : CALL dbcsr_copy(h_re, qs_ot_env%rot_mat_x_im, name="h_re")
2721 3941 : 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 3941 : para_env=qs_ot_env%para_env, blacs_env=qs_ot_env%blacs_env)
2727 :
2728 15764 : ALLOCATE (exp_evals_re(k), exp_evals_im(k))
2729 36059 : exp_evals_re(:) = COS(-qs_ot_env%rot_mat_evals(:))
2730 36059 : exp_evals_im(:) = SIN(-qs_ot_env%rot_mat_evals(:))
2731 :
2732 3941 : CALL dbcsr_copy(w_re, qs_ot_env%rot_mat_evec_re, name="w_re")
2733 3941 : CALL dbcsr_scale_by_vector(w_re, alpha=exp_evals_re, side='right')
2734 3941 : CALL dbcsr_copy(tmp, qs_ot_env%rot_mat_evec_im, name="tmp")
2735 3941 : CALL dbcsr_scale_by_vector(tmp, alpha=exp_evals_im, side='right')
2736 3941 : CALL dbcsr_add(w_re, tmp, alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
2737 :
2738 3941 : CALL dbcsr_copy(w_im, qs_ot_env%rot_mat_evec_im, name="w_im")
2739 3941 : CALL dbcsr_scale_by_vector(w_im, alpha=exp_evals_re, side='right')
2740 3941 : CALL dbcsr_copy(tmp, qs_ot_env%rot_mat_evec_re)
2741 3941 : CALL dbcsr_scale_by_vector(tmp, alpha=exp_evals_im, side='right')
2742 3941 : 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 3941 : 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 3941 : 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 3941 : 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 3941 : 1.0_dp, qs_ot_env%rot_mat_u_im)
2752 :
2753 3941 : CALL dbcsr_release(h_re)
2754 3941 : CALL dbcsr_release(tmp)
2755 3941 : CALL dbcsr_release(w_re)
2756 3941 : CALL dbcsr_release(w_im)
2757 3941 : DEALLOCATE (exp_evals_re, exp_evals_im)
2758 : END IF
2759 :
2760 3941 : CALL timestop(handle)
2761 :
2762 8872 : 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 4126 : 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 2063 : CALL timeset(routineN, handle)
2781 :
2782 2063 : CPASSERT(qs_ot_env%has_complex_kpoint_state)
2783 2063 : CPASSERT(ASSOCIATED(qs_ot_env%rot_mat_dedu_im))
2784 2063 : CPASSERT(ASSOCIATED(qs_ot_env%rot_mat_gx_im))
2785 :
2786 2063 : CALL dbcsr_get_info(qs_ot_env%rot_mat_u, nfullrows_total=k)
2787 2063 : IF (k /= 0) THEN
2788 : rot_norm = dbcsr_frobenius_norm(qs_ot_env%rot_mat_x) + &
2789 2063 : dbcsr_frobenius_norm(qs_ot_env%rot_mat_x_im)
2790 2063 : IF (rot_norm <= EPSILON(1.0_dp)) THEN
2791 474 : CALL qs_ot_square_transpose(qs_ot_env%rot_mat_dedu, tmp, qs_ot_env%rot_mat_u)
2792 474 : 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 474 : alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
2795 474 : CALL dbcsr_release(tmp)
2796 :
2797 474 : CALL qs_ot_square_transpose(qs_ot_env%rot_mat_dedu_im, tmp, qs_ot_env%rot_mat_u)
2798 474 : 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 474 : alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
2801 474 : CALL dbcsr_release(tmp)
2802 474 : CALL timestop(handle)
2803 474 : RETURN
2804 : END IF
2805 :
2806 1589 : CALL dbcsr_copy(work_re, qs_ot_env%rot_mat_dedu, name="work_re")
2807 1589 : CALL dbcsr_copy(work_im, qs_ot_env%rot_mat_dedu, name="work_im")
2808 1589 : CALL dbcsr_copy(tmp, qs_ot_env%rot_mat_dedu, name="tmp")
2809 1589 : CALL dbcsr_copy(inner_deriv_re, qs_ot_env%rot_mat_dedu, name="inner_deriv_re")
2810 1589 : 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 1589 : 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 1589 : qs_ot_env%rot_mat_evec_im, 0.0_dp, tmp)
2817 1589 : 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 1589 : 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 1589 : qs_ot_env%rot_mat_evec_re, 0.0_dp, tmp)
2822 1589 : 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 1589 : 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 1589 : 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 1589 : 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 1589 : 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 1589 : 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 1589 : 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 1589 : outer_deriv_im, 0.0_dp, tmp)
2840 1589 : 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 1589 : 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 1589 : outer_deriv_re, 0.0_dp, tmp)
2845 1589 : CALL dbcsr_add(work_im, tmp, alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
2846 :
2847 1589 : CALL dbcsr_copy(frechet_re, qs_ot_env%rot_mat_dedu, name="frechet_re")
2848 1589 : 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 1589 : 0.0_dp, frechet_re)
2851 : CALL dbcsr_multiply('N', 'T', 1.0_dp, work_im, qs_ot_env%rot_mat_evec_im, &
2852 1589 : 1.0_dp, frechet_re)
2853 : CALL dbcsr_multiply('N', 'T', 1.0_dp, work_im, qs_ot_env%rot_mat_evec_re, &
2854 1589 : 0.0_dp, frechet_im)
2855 : CALL dbcsr_multiply('N', 'T', -1.0_dp, work_re, qs_ot_env%rot_mat_evec_im, &
2856 1589 : 1.0_dp, frechet_im)
2857 :
2858 : ! Tangents satisfy X^T=-X and Y^T=Y for A=X+iY.
2859 1589 : CALL qs_ot_square_transpose(frechet_re, tmp, qs_ot_env%rot_mat_u)
2860 1589 : CALL dbcsr_copy(qs_ot_env%rot_mat_gx, frechet_re)
2861 : CALL dbcsr_add(qs_ot_env%rot_mat_gx, tmp, &
2862 1589 : alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
2863 1589 : CALL qs_ot_square_transpose(frechet_im, tmp, qs_ot_env%rot_mat_u)
2864 1589 : 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 1589 : alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
2867 :
2868 1589 : CALL dbcsr_release(frechet_re)
2869 1589 : CALL dbcsr_release(frechet_im)
2870 1589 : CALL dbcsr_release(inner_deriv_re)
2871 1589 : CALL dbcsr_release(inner_deriv_im)
2872 1589 : CALL dbcsr_release(outer_deriv_re)
2873 1589 : CALL dbcsr_release(outer_deriv_im)
2874 1589 : CALL dbcsr_release(tmp)
2875 1589 : CALL dbcsr_release(work_re)
2876 1589 : CALL dbcsr_release(work_im)
2877 : END IF
2878 :
2879 1589 : CALL timestop(handle)
2880 :
2881 2063 : 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 3552 : 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 3552 : CALL timeset(routineN, handle)
2901 :
2902 3552 : CPASSERT(qs_ot_env%has_complex_kpoint_state)
2903 3552 : 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 3552 : qs_ot_env%matrix_p, qs_ot_env%matrix_p_im, qs_ot_env%matrix_buf1)
2907 3552 : qs_ot_env%do_taylor = .FALSE.
2908 3552 : CALL qs_ot_p2m_diag_complex(qs_ot_env)
2909 :
2910 3552 : CALL timestop(handle)
2911 :
2912 3552 : 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 3158 : 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 3158 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: exp_evals_im, exp_evals_re
2930 : TYPE(dbcsr_type) :: buf_1, buf_2
2931 :
2932 3158 : CALL timeset(routineN, handle)
2933 :
2934 3158 : CALL dbcsr_get_info(qs_ot_env%rot_mat_x, nfullrows_total=k)
2935 :
2936 3158 : 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 3106 : blacs_env=qs_ot_env%blacs_env)
2947 :
2948 : ! Compute: exp_evals = EXP(-i*rot_mat_evals)
2949 12424 : ALLOCATE (exp_evals_re(k), exp_evals_im(k))
2950 16766 : exp_evals_re(:) = COS(-qs_ot_env%rot_mat_evals(:))
2951 16766 : 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 3106 : CALL dbcsr_copy(buf_1, qs_ot_env%rot_mat_evec_re, name="buf_1")
2956 3106 : CALL dbcsr_scale_by_vector(buf_1, alpha=exp_evals_re, side='right')
2957 3106 : CALL dbcsr_copy(buf_2, qs_ot_env%rot_mat_evec_im, name="buf_2")
2958 3106 : CALL dbcsr_scale_by_vector(buf_2, alpha=exp_evals_im, side='right')
2959 3106 : CALL dbcsr_add(buf_1, buf_2, alpha_scalar=+1.0_dp, beta_scalar=-1.0_dp)
2960 3106 : 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 3106 : CALL dbcsr_copy(buf_1, qs_ot_env%rot_mat_evec_im)
2963 3106 : CALL dbcsr_scale_by_vector(buf_1, alpha=exp_evals_re, side='right')
2964 3106 : CALL dbcsr_copy(buf_2, qs_ot_env%rot_mat_evec_re)
2965 3106 : CALL dbcsr_scale_by_vector(buf_2, alpha=exp_evals_im, side='right')
2966 3106 : CALL dbcsr_add(buf_1, buf_2, alpha_scalar=+1.0_dp, beta_scalar=+1.0_dp)
2967 3106 : 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 3106 : CALL dbcsr_release(buf_1)
2971 3106 : CALL dbcsr_release(buf_2)
2972 3106 : DEALLOCATE (exp_evals_re, exp_evals_im)
2973 : END IF
2974 :
2975 3158 : CALL timestop(handle)
2976 :
2977 6316 : 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 1598 : 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 1598 : INTEGER, ALLOCATABLE, DIMENSION(:) :: cols, rows
2994 1598 : 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 1598 : 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 1598 : CALL timeset(routineN, handle)
3006 :
3007 1598 : CALL dbcsr_get_info(qs_ot_env%rot_mat_u, nfullrows_total=k)
3008 1598 : IF (k /= 0) THEN
3009 1572 : 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 1572 : 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 1572 : CALL dbcsr_copy(inner_deriv_re, qs_ot_env%rot_mat_dedu, "inner_deriv_re") ! TODO just create
3015 1572 : CALL dbcsr_copy(inner_deriv_im, qs_ot_env%rot_mat_dedu, "inner_deriv_im") ! TODO just create
3016 :
3017 1572 : 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 1572 : CALL dbcsr_multiply('T', 'N', +1.0_dp, qs_ot_env%rot_mat_evec_im, mat_buf, 0.0_dp, inner_deriv_re)
3019 1572 : 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 1572 : 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 1572 : CALL dbcsr_multiply('T', 'N', +1.0_dp, qs_ot_env%rot_mat_evec_re, mat_buf, 1.0_dp, inner_deriv_re)
3023 1572 : 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 1572 : max_blocks = 0
3028 1572 : CALL dbcsr_iterator_start(iter, inner_deriv_re)
3029 2358 : DO WHILE (dbcsr_iterator_blocks_left(iter))
3030 786 : CALL dbcsr_iterator_next_block(iter, row, col)
3031 786 : max_blocks = max_blocks + 1
3032 : END DO
3033 1572 : CALL dbcsr_iterator_stop(iter)
3034 1572 : CALL dbcsr_iterator_start(iter, inner_deriv_im)
3035 2358 : DO WHILE (dbcsr_iterator_blocks_left(iter))
3036 786 : CALL dbcsr_iterator_next_block(iter, row, col)
3037 786 : max_blocks = max_blocks + 1
3038 : END DO
3039 1572 : CALL dbcsr_iterator_stop(iter)
3040 :
3041 6288 : ALLOCATE (rows(MAX(max_blocks, 1)), cols(MAX(max_blocks, 1)))
3042 1572 : nblocks = 0
3043 1572 : CALL dbcsr_iterator_start(iter, inner_deriv_re)
3044 2358 : DO WHILE (dbcsr_iterator_blocks_left(iter))
3045 786 : CALL dbcsr_iterator_next_block(iter, row, col)
3046 786 : duplicate = .FALSE.
3047 786 : DO iblock = 1, nblocks
3048 786 : duplicate = duplicate .OR. (rows(iblock) == row .AND. cols(iblock) == col)
3049 : END DO
3050 2358 : IF (.NOT. duplicate) THEN
3051 786 : nblocks = nblocks + 1
3052 786 : rows(nblocks) = row
3053 786 : cols(nblocks) = col
3054 : END IF
3055 : END DO
3056 1572 : CALL dbcsr_iterator_stop(iter)
3057 1572 : CALL dbcsr_iterator_start(iter, inner_deriv_im)
3058 2358 : DO WHILE (dbcsr_iterator_blocks_left(iter))
3059 786 : CALL dbcsr_iterator_next_block(iter, row, col)
3060 786 : duplicate = .FALSE.
3061 1572 : DO iblock = 1, nblocks
3062 1572 : duplicate = duplicate .OR. (rows(iblock) == row .AND. cols(iblock) == col)
3063 : END DO
3064 2358 : IF (.NOT. duplicate) THEN
3065 0 : nblocks = nblocks + 1
3066 0 : rows(nblocks) = row
3067 0 : cols(nblocks) = col
3068 : END IF
3069 : END DO
3070 1572 : CALL dbcsr_iterator_stop(iter)
3071 :
3072 : CALL dbcsr_get_info(inner_deriv_re, distribution=dist, row_blk_size=row_blk_size, &
3073 1572 : col_blk_size=col_blk_size)
3074 : CALL dbcsr_create(outer_deriv_re, "outer_deriv_re", dist, dbcsr_type_no_symmetry, &
3075 1572 : row_blk_size, col_blk_size)
3076 : CALL dbcsr_create(outer_deriv_im, "outer_deriv_im", dist, dbcsr_type_no_symmetry, &
3077 1572 : row_blk_size, col_blk_size)
3078 1572 : IF (nblocks > 0) THEN
3079 786 : CALL dbcsr_reserve_blocks(outer_deriv_re, rows=rows(1:nblocks), cols=cols(1:nblocks))
3080 786 : CALL dbcsr_reserve_blocks(outer_deriv_im, rows=rows(1:nblocks), cols=cols(1:nblocks))
3081 : END IF
3082 1572 : CALL dbcsr_finalize(outer_deriv_re)
3083 1572 : CALL dbcsr_finalize(outer_deriv_im)
3084 1572 : CALL dbcsr_set(outer_deriv_re, 0.0_dp)
3085 1572 : CALL dbcsr_set(outer_deriv_im, 0.0_dp)
3086 :
3087 1572 : CALL dbcsr_get_info(outer_deriv_re, row_blk_offset=row_blk_offset, col_blk_offset=col_blk_offset)
3088 1572 : CALL dbcsr_iterator_start(iter, outer_deriv_re)
3089 2358 : DO WHILE (dbcsr_iterator_blocks_left(iter))
3090 786 : CALL dbcsr_iterator_next_block(iter, row, col)
3091 786 : CALL dbcsr_get_block_p(inner_deriv_re, row, col, block_in_re, found_in_re)
3092 786 : CALL dbcsr_get_block_p(inner_deriv_im, row, col, block_in_im, found_in_im)
3093 786 : CALL dbcsr_get_block_p(outer_deriv_re, row, col, block_out_re, found_out_re)
3094 786 : CALL dbcsr_get_block_p(outer_deriv_im, row, col, block_out_im, found_out_im)
3095 786 : CPASSERT(found_out_re .AND. found_out_im)
3096 :
3097 5953 : DO i = 1, SIZE(block_out_re, 1)
3098 25752 : DO j = 1, SIZE(block_out_re, 2)
3099 21371 : e1 = qs_ot_env%rot_mat_evals(row_blk_offset(row) + i - 1)
3100 21371 : e2 = qs_ot_env%rot_mat_evals(col_blk_offset(col) + j - 1)
3101 21371 : re_part = 0.0_dp
3102 21371 : im_part = 0.0_dp
3103 21371 : IF (found_in_re) re_part = block_in_re(i, j)
3104 21371 : IF (found_in_im) im_part = block_in_im(i, j)
3105 21371 : cval_in = CMPLX(re_part, im_part, dp)
3106 21371 : cval_out = cval_in*cint(e1, e2)
3107 21371 : block_out_re(i, j) = REAL(cval_out)
3108 24966 : block_out_im(i, j) = AIMAG(cval_out)
3109 : END DO
3110 : END DO
3111 : END DO
3112 1572 : CALL dbcsr_iterator_stop(iter)
3113 1572 : DEALLOCATE (rows, cols)
3114 1572 : CALL dbcsr_release(inner_deriv_re)
3115 1572 : 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 1572 : CALL dbcsr_multiply('N', 'N', +1.0_dp, qs_ot_env%rot_mat_evec_re, outer_deriv_re, 0.0_dp, mat_buf)
3119 1572 : CALL dbcsr_multiply('N', 'N', -1.0_dp, qs_ot_env%rot_mat_evec_im, outer_deriv_im, 1.0_dp, mat_buf)
3120 1572 : 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 1572 : CALL dbcsr_multiply('N', 'N', +1.0_dp, qs_ot_env%rot_mat_evec_re, outer_deriv_im, 0.0_dp, mat_buf)
3123 1572 : CALL dbcsr_multiply('N', 'N', +1.0_dp, qs_ot_env%rot_mat_evec_im, outer_deriv_re, 1.0_dp, mat_buf)
3124 1572 : 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 1572 : CALL qs_ot_square_transpose(qs_ot_env%matrix_buf1, mat_buf, qs_ot_env%rot_mat_u)
3129 1572 : 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 1572 : alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
3132 :
3133 1572 : CALL dbcsr_release(mat_buf)
3134 1572 : CALL dbcsr_release(outer_deriv_re)
3135 11004 : CALL dbcsr_release(outer_deriv_im)
3136 : END IF
3137 3196 : CALL timestop(handle)
3138 : CONTAINS
3139 :
3140 : ! **************************************************************************************************
3141 : !> \brief ...
3142 : !> \param e1 ...
3143 : !> \param e2 ...
3144 : !> \return ...
3145 : ! **************************************************************************************************
3146 21371 : 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 21371 : l1 = (0.0_dp, -1.0_dp)*e1
3154 21371 : l2 = (0.0_dp, -1.0_dp)*e2
3155 21371 : IF (ABS(l1 - l2) > 0.5_dp) THEN
3156 1428 : cint = (EXP(l1) - EXP(l2))/(l1 - l2)
3157 : ELSE
3158 : x = 1.0_dp
3159 : cint = 0.0_dp
3160 339031 : DO I = 1, 16
3161 319088 : cint = cint + x
3162 339031 : x = x*(l1 - l2)/REAL(I + 1, KIND=dp)
3163 : END DO
3164 19943 : cint = cint*EXP(l2)
3165 : END IF
3166 21371 : 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 111015 : 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 111015 : qs_ot_env%do_taylor = .FALSE.
3185 111015 : N = 0
3186 111015 : num_error = qs_ot_env%largest_eval_upper_bound/(2.0_dp)
3187 470761 : DO WHILE (num_error > qs_ot_env%settings%eps_taylor .AND. N <= qs_ot_env%settings%max_taylor)
3188 359746 : N = N + 1
3189 405577 : num_error = num_error*qs_ot_env%largest_eval_upper_bound/REAL((2*N + 1)*(2*N + 2), KIND=dp)
3190 : END DO
3191 : ! An O(P) matrix-function term can be below eps_taylor while its O(X)
3192 : ! orbital derivative is still needed near P=X^T*S*X=0.
3193 111015 : qs_ot_env%taylor_order = MAX(1, N)
3194 111015 : IF (qs_ot_env%taylor_order <= qs_ot_env%settings%max_taylor) THEN
3195 59916 : qs_ot_env%do_taylor = .TRUE.
3196 : END IF
3197 :
3198 111015 : END SUBROUTINE decide_strategy
3199 :
3200 : ! **************************************************************************************************
3201 : !> \brief c=(c0*cos(p^0.5)+x*sin(p^0.5)*p^(-0.5)) x rot_mat_u
3202 : !> this assumes that x is already ortho to S*C0, and that p is x*S*x
3203 : !> rot_mat_u is an optional rotation matrix
3204 : !> \param matrix_c ...
3205 : !> \param matrix_x ...
3206 : !> \param qs_ot_env ...
3207 : ! **************************************************************************************************
3208 206150 : SUBROUTINE qs_ot_get_orbitals(matrix_c, matrix_x, qs_ot_env)
3209 :
3210 : TYPE(dbcsr_type), POINTER :: matrix_c, matrix_x
3211 : TYPE(qs_ot_type) :: qs_ot_env
3212 :
3213 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_ot_get_orbitals'
3214 : REAL(KIND=dp), PARAMETER :: rone = 1.0_dp, rzero = 0.0_dp
3215 :
3216 : INTEGER :: handle, k, n
3217 : TYPE(dbcsr_type), POINTER :: matrix_kk
3218 :
3219 103075 : CALL timeset(routineN, handle)
3220 :
3221 103075 : CALL dbcsr_get_info(matrix_x, nfullrows_total=n, nfullcols_total=k)
3222 :
3223 : ! rotate the multiplying matrices cosp and sinp instead of the result,
3224 : ! this should be cheaper for large basis sets
3225 103075 : IF (qs_ot_env%settings%do_rotation) THEN
3226 2938 : matrix_kk => qs_ot_env%matrix_buf1
3227 : CALL dbcsr_multiply('N', 'N', rone, qs_ot_env%matrix_cosp, &
3228 2938 : qs_ot_env%rot_mat_u, rzero, matrix_kk)
3229 : ELSE
3230 100137 : matrix_kk => qs_ot_env%matrix_cosp
3231 : END IF
3232 :
3233 : CALL dbcsr_multiply('N', 'N', rone, qs_ot_env%matrix_c0, matrix_kk, &
3234 103075 : rzero, matrix_c)
3235 :
3236 103075 : IF (qs_ot_env%settings%do_rotation) THEN
3237 2938 : matrix_kk => qs_ot_env%matrix_buf1
3238 : CALL dbcsr_multiply('N', 'N', rone, qs_ot_env%matrix_sinp, &
3239 2938 : qs_ot_env%rot_mat_u, rzero, matrix_kk)
3240 : ELSE
3241 100137 : matrix_kk => qs_ot_env%matrix_sinp
3242 : END IF
3243 : CALL dbcsr_multiply('N', 'N', rone, matrix_x, matrix_kk, &
3244 103075 : rone, matrix_c)
3245 :
3246 103075 : CALL timestop(handle)
3247 :
3248 103075 : END SUBROUTINE qs_ot_get_orbitals
3249 :
3250 : ! **************************************************************************************************
3251 : !> \brief update complex K-point orbitals with the finite STRICT transformation
3252 : !> \param matrix_c real output orbitals
3253 : !> \param matrix_c_im imaginary output orbitals
3254 : !> \param matrix_s real overlap matrix
3255 : !> \param matrix_s_im imaginary overlap matrix
3256 : !> \param qs_ot_env OT channel state
3257 : ! **************************************************************************************************
3258 3160 : SUBROUTINE qs_ot_get_orbitals_complex(matrix_c, matrix_c_im, matrix_s, matrix_s_im, qs_ot_env)
3259 : TYPE(dbcsr_type), POINTER :: matrix_c, matrix_c_im, matrix_s, &
3260 : matrix_s_im
3261 : TYPE(qs_ot_type) :: qs_ot_env
3262 :
3263 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_ot_get_orbitals_complex'
3264 :
3265 : INTEGER :: handle
3266 : TYPE(dbcsr_type) :: rotated_im, rotated_re, rotation_tmp
3267 :
3268 3160 : CALL timeset(routineN, handle)
3269 :
3270 3160 : CPASSERT(qs_ot_env%has_complex_kpoint_state)
3271 3160 : CPASSERT(qs_ot_env%settings%ot_algorithm == 'TOD')
3272 :
3273 : CALL qs_ot_complex_multiply('N', 'N', matrix_s, matrix_s_im, &
3274 : qs_ot_env%matrix_x, qs_ot_env%matrix_x_im, &
3275 : qs_ot_env%matrix_sx, qs_ot_env%matrix_sx_im, &
3276 3160 : qs_ot_env%matrix_tmp_nk)
3277 : CALL qs_ot_get_p_complex(qs_ot_env%matrix_x, qs_ot_env%matrix_x_im, &
3278 3160 : qs_ot_env%matrix_sx, qs_ot_env%matrix_sx_im, qs_ot_env)
3279 :
3280 : CALL qs_ot_complex_multiply('N', 'N', qs_ot_env%matrix_c0, qs_ot_env%matrix_c0_im, &
3281 : qs_ot_env%matrix_cosp, qs_ot_env%matrix_cosp_im, &
3282 3160 : matrix_c, matrix_c_im, qs_ot_env%matrix_tmp_nk)
3283 : CALL qs_ot_complex_multiply('N', 'N', qs_ot_env%matrix_x, qs_ot_env%matrix_x_im, &
3284 : qs_ot_env%matrix_sinp, qs_ot_env%matrix_sinp_im, &
3285 : qs_ot_env%matrix_buf_nk, qs_ot_env%matrix_buf_nk_im, &
3286 3160 : qs_ot_env%matrix_tmp_nk)
3287 3160 : CALL dbcsr_add(matrix_c, qs_ot_env%matrix_buf_nk, alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
3288 3160 : CALL dbcsr_add(matrix_c_im, qs_ot_env%matrix_buf_nk_im, alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
3289 :
3290 3160 : IF (qs_ot_env%settings%do_rotation) THEN
3291 684 : CALL qs_ot_generate_rotation_complex(qs_ot_env)
3292 :
3293 684 : CALL dbcsr_copy(rotated_re, matrix_c, name="strict_rotated_re")
3294 684 : CALL dbcsr_copy(rotated_im, matrix_c_im, name="strict_rotated_im")
3295 684 : CALL dbcsr_copy(rotation_tmp, matrix_c, name="strict_rotation_tmp")
3296 :
3297 : ! C_out = Q(X)*U for the finite STRICT chart Q(X).
3298 : CALL dbcsr_multiply('N', 'N', 1.0_dp, matrix_c, qs_ot_env%rot_mat_u, &
3299 684 : 0.0_dp, rotated_re)
3300 : CALL dbcsr_multiply('N', 'N', 1.0_dp, matrix_c_im, qs_ot_env%rot_mat_u_im, &
3301 684 : 0.0_dp, rotation_tmp)
3302 684 : CALL dbcsr_add(rotated_re, rotation_tmp, alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
3303 :
3304 : CALL dbcsr_multiply('N', 'N', 1.0_dp, matrix_c, qs_ot_env%rot_mat_u_im, &
3305 684 : 0.0_dp, rotated_im)
3306 : CALL dbcsr_multiply('N', 'N', 1.0_dp, matrix_c_im, qs_ot_env%rot_mat_u, &
3307 684 : 0.0_dp, rotation_tmp)
3308 684 : CALL dbcsr_add(rotated_im, rotation_tmp, alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
3309 :
3310 684 : CALL dbcsr_copy(matrix_c, rotated_re)
3311 684 : CALL dbcsr_copy(matrix_c_im, rotated_im)
3312 684 : CALL dbcsr_release(rotated_re)
3313 684 : CALL dbcsr_release(rotated_im)
3314 684 : CALL dbcsr_release(rotation_tmp)
3315 : END IF
3316 :
3317 3160 : CALL timestop(handle)
3318 :
3319 3160 : END SUBROUTINE qs_ot_get_orbitals_complex
3320 :
3321 : ! **************************************************************************************************
3322 : !> \brief this routines computes dE/dx=dx, with dx ortho to sc0
3323 : !> needs dE/dC=hc,C0,X,SX,p
3324 : !> if preconditioned it will not be the derivative, but the lagrangian multiplier
3325 : !> is changed so that P*dE/dx is the right derivative (i.e. in the allowed subspace)
3326 : !> \param matrix_hc ...
3327 : !> \param matrix_x ...
3328 : !> \param matrix_sx ...
3329 : !> \param matrix_gx ...
3330 : !> \param qs_ot_env ...
3331 : ! **************************************************************************************************
3332 238683 : SUBROUTINE qs_ot_get_derivative(matrix_hc, matrix_x, matrix_sx, matrix_gx, &
3333 : qs_ot_env)
3334 : TYPE(dbcsr_type), POINTER :: matrix_hc, matrix_x, matrix_sx, matrix_gx
3335 : TYPE(qs_ot_type) :: qs_ot_env
3336 :
3337 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_ot_get_derivative'
3338 : REAL(KIND=dp), PARAMETER :: rone = 1.0_dp, rzero = 0.0_dp
3339 :
3340 : INTEGER :: handle, k, n, ortho_k
3341 : TYPE(dbcsr_type), POINTER :: matrix_hc_local, matrix_target
3342 :
3343 79561 : CALL timeset(routineN, handle)
3344 :
3345 79561 : NULLIFY (matrix_hc_local)
3346 :
3347 79561 : CALL dbcsr_get_info(matrix_x, nfullrows_total=n, nfullcols_total=k)
3348 :
3349 : ! could in principle be taken inside qs_ot_get_derivative_* for increased efficiency
3350 : ! create a local rotated version of matrix_hc leaving matrix_hc untouched (needed
3351 : ! for lagrangian multipliers)
3352 79561 : IF (qs_ot_env%settings%do_rotation) THEN
3353 1588 : CALL dbcsr_copy(matrix_gx, matrix_hc) ! use gx as temporary
3354 1588 : CALL dbcsr_init_p(matrix_hc_local)
3355 1588 : CALL dbcsr_copy(matrix_hc_local, matrix_hc, name='matrix_hc_local')
3356 1588 : CALL dbcsr_set(matrix_hc_local, 0.0_dp)
3357 1588 : CALL dbcsr_multiply('N', 'T', rone, matrix_gx, qs_ot_env%rot_mat_u, rzero, matrix_hc_local)
3358 : ELSE
3359 77973 : matrix_hc_local => matrix_hc
3360 : END IF
3361 :
3362 79561 : IF (qs_ot_env%do_taylor) THEN
3363 44115 : CALL qs_ot_get_derivative_taylor(matrix_hc_local, matrix_x, matrix_sx, matrix_gx, qs_ot_env)
3364 : ELSE
3365 35446 : CALL qs_ot_get_derivative_diag(matrix_hc_local, matrix_x, matrix_sx, matrix_gx, qs_ot_env)
3366 : END IF
3367 :
3368 : ! and make it orthogonal
3369 79561 : CALL dbcsr_get_info(qs_ot_env%matrix_sc0, nfullcols_total=ortho_k)
3370 :
3371 79561 : IF (ASSOCIATED(qs_ot_env%preconditioner)) THEN
3372 69063 : matrix_target => qs_ot_env%matrix_psc0
3373 : ELSE
3374 10498 : matrix_target => qs_ot_env%matrix_sc0
3375 : END IF
3376 : ! first make the matrix os if not yet valid
3377 79561 : IF (.NOT. qs_ot_env%os_valid) THEN
3378 : ! this assumes that the preconditioner is a single matrix
3379 : ! that maps sc0 onto psc0
3380 :
3381 8626 : IF (ASSOCIATED(qs_ot_env%preconditioner)) THEN
3382 : CALL apply_preconditioner(qs_ot_env%preconditioner, qs_ot_env%matrix_sc0, &
3383 7628 : qs_ot_env%matrix_psc0)
3384 : END IF
3385 : CALL dbcsr_multiply('T', 'N', rone, &
3386 : qs_ot_env%matrix_sc0, matrix_target, &
3387 8626 : rzero, qs_ot_env%matrix_os)
3388 : CALL cp_dbcsr_cholesky_decompose(qs_ot_env%matrix_os, &
3389 8626 : para_env=qs_ot_env%para_env, blacs_env=qs_ot_env%blacs_env)
3390 : CALL cp_dbcsr_cholesky_invert(qs_ot_env%matrix_os, &
3391 : para_env=qs_ot_env%para_env, blacs_env=qs_ot_env%blacs_env, &
3392 8626 : uplo_to_full=.TRUE.)
3393 8626 : qs_ot_env%os_valid = .TRUE.
3394 : END IF
3395 : CALL dbcsr_multiply('T', 'N', rone, matrix_target, matrix_gx, &
3396 79561 : rzero, qs_ot_env%matrix_buf1_ortho)
3397 : CALL dbcsr_multiply('N', 'N', rone, qs_ot_env%matrix_os, &
3398 79561 : qs_ot_env%matrix_buf1_ortho, rzero, qs_ot_env%matrix_buf2_ortho)
3399 : CALL dbcsr_multiply('N', 'N', -rone, qs_ot_env%matrix_sc0, &
3400 79561 : qs_ot_env%matrix_buf2_ortho, rone, matrix_gx)
3401 : ! also treat the rot_mat gradient here
3402 79561 : IF (qs_ot_env%settings%do_rotation) THEN
3403 1588 : CALL qs_ot_rot_mat_derivative(qs_ot_env)
3404 : END IF
3405 :
3406 79561 : IF (qs_ot_env%settings%do_rotation) THEN
3407 1588 : CALL dbcsr_release_p(matrix_hc_local)
3408 : END IF
3409 :
3410 79561 : CALL timestop(handle)
3411 :
3412 79561 : END SUBROUTINE qs_ot_get_derivative
3413 :
3414 : ! **************************************************************************************************
3415 : !> \brief Prepare the inverse metric used to project a complex STRICT gradient.
3416 : !> An unusable preconditioner is detached before any minimizer history is updated.
3417 : !> \param qs_ot_env OT channel state
3418 : !> \param preconditioner_rejected true if the attached preconditioner was not positive definite
3419 : ! **************************************************************************************************
3420 3436 : SUBROUTINE qs_ot_prepare_complex_tangent_metric(qs_ot_env, preconditioner_rejected)
3421 : TYPE(qs_ot_type) :: qs_ot_env
3422 : LOGICAL, INTENT(OUT), OPTIONAL :: preconditioner_rejected
3423 :
3424 : INTEGER :: i, k
3425 : REAL(KIND=dp) :: eval_scale, eval_threshold
3426 : TYPE(dbcsr_type), POINTER :: target_im, target_re
3427 :
3428 2944 : IF (PRESENT(preconditioner_rejected)) preconditioner_rejected = .FALSE.
3429 2944 : IF (qs_ot_env%os_valid) RETURN
3430 :
3431 492 : IF (ASSOCIATED(qs_ot_env%preconditioner)) THEN
3432 482 : target_re => qs_ot_env%matrix_psc0
3433 482 : target_im => qs_ot_env%matrix_psc0_im
3434 : CALL apply_preconditioner(qs_ot_env%preconditioner, &
3435 : qs_ot_env%matrix_sc0, qs_ot_env%matrix_sc0_im, &
3436 482 : target_re, target_im)
3437 : ELSE
3438 10 : target_re => qs_ot_env%matrix_sc0
3439 10 : target_im => qs_ot_env%matrix_sc0_im
3440 : END IF
3441 : CALL qs_ot_complex_multiply('C', 'N', qs_ot_env%matrix_sc0, qs_ot_env%matrix_sc0_im, &
3442 : target_re, target_im, qs_ot_env%matrix_os, qs_ot_env%matrix_os_im, &
3443 492 : qs_ot_env%matrix_buf1)
3444 492 : CALL dbcsr_get_info(qs_ot_env%matrix_os, nfullrows_total=k)
3445 : CALL cp_dbcsr_heevd(matrix_re=qs_ot_env%matrix_os, matrix_im=qs_ot_env%matrix_os_im, &
3446 : eigenvectors_re=qs_ot_env%matrix_r, eigenvectors_im=qs_ot_env%matrix_r_im, &
3447 : eigenvalues=qs_ot_env%evals, para_env=qs_ot_env%para_env, &
3448 492 : blacs_env=qs_ot_env%blacs_env)
3449 2984 : eval_scale = MAX(1.0_dp, MAXVAL(ABS(qs_ot_env%evals(1:k))))
3450 492 : eval_threshold = 100.0_dp*EPSILON(1.0_dp)*eval_scale
3451 2984 : IF (MINVAL(qs_ot_env%evals(1:k)) <= eval_threshold .AND. &
3452 : ASSOCIATED(qs_ot_env%preconditioner)) THEN
3453 0 : NULLIFY (qs_ot_env%preconditioner)
3454 0 : IF (PRESENT(preconditioner_rejected)) preconditioner_rejected = .TRUE.
3455 0 : target_re => qs_ot_env%matrix_sc0
3456 0 : target_im => qs_ot_env%matrix_sc0_im
3457 : CALL qs_ot_complex_multiply('C', 'N', qs_ot_env%matrix_sc0, qs_ot_env%matrix_sc0_im, &
3458 : target_re, target_im, qs_ot_env%matrix_os, qs_ot_env%matrix_os_im, &
3459 0 : qs_ot_env%matrix_buf1)
3460 : CALL cp_dbcsr_heevd(matrix_re=qs_ot_env%matrix_os, matrix_im=qs_ot_env%matrix_os_im, &
3461 : eigenvectors_re=qs_ot_env%matrix_r, eigenvectors_im=qs_ot_env%matrix_r_im, &
3462 : eigenvalues=qs_ot_env%evals, para_env=qs_ot_env%para_env, &
3463 0 : blacs_env=qs_ot_env%blacs_env)
3464 0 : eval_scale = MAX(1.0_dp, MAXVAL(ABS(qs_ot_env%evals(1:k))))
3465 0 : eval_threshold = 100.0_dp*EPSILON(1.0_dp)*eval_scale
3466 : END IF
3467 2984 : CPASSERT(MINVAL(qs_ot_env%evals(1:k)) > eval_threshold)
3468 2984 : DO i = 1, k
3469 2984 : qs_ot_env%dum(i) = 1.0_dp/qs_ot_env%evals(i)
3470 : END DO
3471 492 : CALL dbcsr_copy(qs_ot_env%matrix_buf1, qs_ot_env%matrix_r)
3472 492 : CALL dbcsr_copy(qs_ot_env%matrix_buf1_im, qs_ot_env%matrix_r_im)
3473 492 : CALL dbcsr_scale_by_vector(qs_ot_env%matrix_buf1, alpha=qs_ot_env%dum, side='right')
3474 492 : CALL dbcsr_scale_by_vector(qs_ot_env%matrix_buf1_im, alpha=qs_ot_env%dum, side='right')
3475 : CALL qs_ot_complex_multiply('N', 'C', qs_ot_env%matrix_buf1, qs_ot_env%matrix_buf1_im, &
3476 : qs_ot_env%matrix_r, qs_ot_env%matrix_r_im, &
3477 : qs_ot_env%matrix_os, qs_ot_env%matrix_os_im, &
3478 492 : qs_ot_env%matrix_buf2)
3479 492 : qs_ot_env%os_valid = .TRUE.
3480 :
3481 : END SUBROUTINE qs_ot_prepare_complex_tangent_metric
3482 :
3483 : ! **************************************************************************************************
3484 : !> \brief finite complex STRICT derivative, projected onto C0^H*S*X=0
3485 : !> \param matrix_hc real part of H(k)*C(k)
3486 : !> \param matrix_hc_im imaginary part of H(k)*C(k)
3487 : !> \param qs_ot_env OT channel state
3488 : !> \param matrix_hc_rotation occupation-weighted H(k)C(k) for the rotation channel
3489 : !> \param matrix_hc_rotation_im imaginary component of matrix_hc_rotation
3490 : ! **************************************************************************************************
3491 5688 : SUBROUTINE qs_ot_get_derivative_complex(matrix_hc, matrix_hc_im, qs_ot_env, &
3492 : matrix_hc_rotation, matrix_hc_rotation_im)
3493 : TYPE(dbcsr_type), POINTER :: matrix_hc, matrix_hc_im
3494 : TYPE(qs_ot_type) :: qs_ot_env
3495 : TYPE(dbcsr_type), OPTIONAL, POINTER :: matrix_hc_rotation, matrix_hc_rotation_im
3496 :
3497 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_ot_get_derivative_complex'
3498 :
3499 : INTEGER :: handle
3500 : TYPE(dbcsr_distribution_type) :: dist
3501 : TYPE(dbcsr_type) :: tmp_nk
3502 : TYPE(dbcsr_type), POINTER :: hc_rotation_im, hc_rotation_re, &
3503 : hc_work_im, hc_work_re, target_im, &
3504 : target_re
3505 : TYPE(dbcsr_type), TARGET :: hc_rot_im, hc_rot_re
3506 :
3507 2844 : CALL timeset(routineN, handle)
3508 :
3509 2844 : CPASSERT(qs_ot_env%has_complex_kpoint_state)
3510 2844 : CPASSERT(qs_ot_env%settings%ot_algorithm == 'TOD')
3511 :
3512 2844 : hc_rotation_re => matrix_hc
3513 2844 : hc_rotation_im => matrix_hc_im
3514 2844 : IF (PRESENT(matrix_hc_rotation) .OR. PRESENT(matrix_hc_rotation_im)) THEN
3515 284 : CPASSERT(PRESENT(matrix_hc_rotation) .AND. PRESENT(matrix_hc_rotation_im))
3516 284 : CPASSERT(ASSOCIATED(matrix_hc_rotation))
3517 284 : CPASSERT(ASSOCIATED(matrix_hc_rotation_im))
3518 284 : hc_rotation_re => matrix_hc_rotation
3519 284 : hc_rotation_im => matrix_hc_rotation_im
3520 : END IF
3521 2844 : hc_work_re => matrix_hc
3522 2844 : hc_work_im => matrix_hc_im
3523 :
3524 2844 : IF (qs_ot_env%settings%do_rotation) THEN
3525 648 : CALL qs_ot_generate_rotation_complex(qs_ot_env)
3526 :
3527 : ! Reconstruct the unrotated finite STRICT orbitals Q(X).
3528 : CALL qs_ot_complex_multiply('N', 'N', qs_ot_env%matrix_c0, qs_ot_env%matrix_c0_im, &
3529 : qs_ot_env%matrix_cosp, qs_ot_env%matrix_cosp_im, &
3530 : qs_ot_env%matrix_buf_nk, qs_ot_env%matrix_buf_nk_im, &
3531 648 : qs_ot_env%matrix_tmp_nk)
3532 : CALL qs_ot_complex_multiply('N', 'N', qs_ot_env%matrix_x, qs_ot_env%matrix_x_im, &
3533 : qs_ot_env%matrix_sinp, qs_ot_env%matrix_sinp_im, &
3534 : qs_ot_env%matrix_gx, qs_ot_env%matrix_gx_im, &
3535 648 : qs_ot_env%matrix_tmp_nk)
3536 : CALL dbcsr_add(qs_ot_env%matrix_buf_nk, qs_ot_env%matrix_gx, &
3537 648 : alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
3538 : CALL dbcsr_add(qs_ot_env%matrix_buf_nk_im, qs_ot_env%matrix_gx_im, &
3539 648 : alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
3540 :
3541 : ! dF/dU = Q(X)^H*G_C for C=Q(X)*U.
3542 : CALL qs_ot_complex_multiply('C', 'N', &
3543 : qs_ot_env%matrix_buf_nk, qs_ot_env%matrix_buf_nk_im, &
3544 : hc_rotation_re, hc_rotation_im, &
3545 : qs_ot_env%rot_mat_dedu, qs_ot_env%rot_mat_dedu_im, &
3546 648 : qs_ot_env%matrix_buf1)
3547 648 : CALL qs_ot_rot_mat_derivative_complex(qs_ot_env)
3548 :
3549 : ! The STRICT coordinate sees G_Q=G_C*U^H.
3550 648 : CALL dbcsr_copy(hc_rot_re, matrix_hc, name="strict_hc_rot_re")
3551 648 : CALL dbcsr_copy(hc_rot_im, matrix_hc_im, name="strict_hc_rot_im")
3552 648 : CALL dbcsr_copy(tmp_nk, matrix_hc, name="strict_hc_rot_tmp")
3553 : CALL dbcsr_multiply('N', 'T', 1.0_dp, matrix_hc, qs_ot_env%rot_mat_u, &
3554 648 : 0.0_dp, hc_rot_re)
3555 : CALL dbcsr_multiply('N', 'T', 1.0_dp, matrix_hc_im, qs_ot_env%rot_mat_u_im, &
3556 648 : 0.0_dp, tmp_nk)
3557 648 : CALL dbcsr_add(hc_rot_re, tmp_nk, alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
3558 :
3559 : CALL dbcsr_multiply('N', 'T', 1.0_dp, matrix_hc_im, qs_ot_env%rot_mat_u, &
3560 648 : 0.0_dp, hc_rot_im)
3561 : CALL dbcsr_multiply('N', 'T', 1.0_dp, matrix_hc, qs_ot_env%rot_mat_u_im, &
3562 648 : 0.0_dp, tmp_nk)
3563 648 : CALL dbcsr_add(hc_rot_im, tmp_nk, alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
3564 648 : hc_work_re => hc_rot_re
3565 648 : hc_work_im => hc_rot_im
3566 : END IF
3567 :
3568 : ! Direct X contribution, H*C sinc(sqrt(P)).
3569 : CALL qs_ot_complex_multiply('N', 'N', hc_work_re, hc_work_im, &
3570 : qs_ot_env%matrix_sinp, qs_ot_env%matrix_sinp_im, &
3571 : qs_ot_env%matrix_gx, qs_ot_env%matrix_gx_im, &
3572 2844 : qs_ot_env%matrix_tmp_nk)
3573 :
3574 : ! Frechet contribution from X sinc(sqrt(P)).
3575 : CALL qs_ot_complex_multiply('C', 'N', hc_work_re, hc_work_im, &
3576 : qs_ot_env%matrix_x, qs_ot_env%matrix_x_im, &
3577 : qs_ot_env%matrix_buf2, qs_ot_env%matrix_buf2_im, &
3578 2844 : qs_ot_env%matrix_buf1)
3579 : CALL qs_ot_complex_multiply('N', 'N', qs_ot_env%matrix_buf2, qs_ot_env%matrix_buf2_im, &
3580 : qs_ot_env%matrix_r, qs_ot_env%matrix_r_im, &
3581 : qs_ot_env%matrix_buf1, qs_ot_env%matrix_buf1_im, &
3582 2844 : qs_ot_env%matrix_buf4)
3583 : CALL qs_ot_complex_multiply('C', 'N', qs_ot_env%matrix_r, qs_ot_env%matrix_r_im, &
3584 : qs_ot_env%matrix_buf1, qs_ot_env%matrix_buf1_im, &
3585 : qs_ot_env%matrix_buf2, qs_ot_env%matrix_buf2_im, &
3586 2844 : qs_ot_env%matrix_buf4)
3587 : CALL dbcsr_hadamard_product(qs_ot_env%matrix_buf2, qs_ot_env%matrix_sinp_b, &
3588 2844 : qs_ot_env%matrix_buf3)
3589 : CALL dbcsr_hadamard_product(qs_ot_env%matrix_buf2_im, qs_ot_env%matrix_sinp_b, &
3590 2844 : qs_ot_env%matrix_buf3_im)
3591 :
3592 : ! Frechet contribution from C0 cos(sqrt(P)).
3593 : CALL qs_ot_complex_multiply('C', 'N', hc_work_re, hc_work_im, &
3594 : qs_ot_env%matrix_c0, qs_ot_env%matrix_c0_im, &
3595 : qs_ot_env%matrix_buf2, qs_ot_env%matrix_buf2_im, &
3596 2844 : qs_ot_env%matrix_buf1)
3597 : CALL qs_ot_complex_multiply('N', 'N', qs_ot_env%matrix_buf2, qs_ot_env%matrix_buf2_im, &
3598 : qs_ot_env%matrix_r, qs_ot_env%matrix_r_im, &
3599 : qs_ot_env%matrix_buf1, qs_ot_env%matrix_buf1_im, &
3600 2844 : qs_ot_env%matrix_buf4)
3601 : CALL qs_ot_complex_multiply('C', 'N', qs_ot_env%matrix_r, qs_ot_env%matrix_r_im, &
3602 : qs_ot_env%matrix_buf1, qs_ot_env%matrix_buf1_im, &
3603 : qs_ot_env%matrix_buf2, qs_ot_env%matrix_buf2_im, &
3604 2844 : qs_ot_env%matrix_buf4)
3605 : CALL dbcsr_hadamard_product(qs_ot_env%matrix_buf2, qs_ot_env%matrix_cosp_b, &
3606 2844 : qs_ot_env%matrix_buf4)
3607 : CALL dbcsr_hadamard_product(qs_ot_env%matrix_buf2_im, qs_ot_env%matrix_cosp_b, &
3608 2844 : qs_ot_env%matrix_buf4_im)
3609 : CALL dbcsr_add(qs_ot_env%matrix_buf3, qs_ot_env%matrix_buf4, &
3610 2844 : alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
3611 : CALL dbcsr_add(qs_ot_env%matrix_buf3_im, qs_ot_env%matrix_buf4_im, &
3612 2844 : alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
3613 :
3614 : ! Transform back and add the Hermitian adjoint generated by dP.
3615 : CALL qs_ot_complex_multiply('N', 'N', qs_ot_env%matrix_r, qs_ot_env%matrix_r_im, &
3616 : qs_ot_env%matrix_buf3, qs_ot_env%matrix_buf3_im, &
3617 : qs_ot_env%matrix_buf1, qs_ot_env%matrix_buf1_im, &
3618 2844 : qs_ot_env%matrix_buf2)
3619 : CALL qs_ot_complex_multiply('N', 'C', qs_ot_env%matrix_buf1, qs_ot_env%matrix_buf1_im, &
3620 : qs_ot_env%matrix_r, qs_ot_env%matrix_r_im, &
3621 : qs_ot_env%matrix_buf3, qs_ot_env%matrix_buf3_im, &
3622 2844 : qs_ot_env%matrix_buf2)
3623 2844 : CALL dbcsr_get_info(qs_ot_env%matrix_buf3, distribution=dist)
3624 : CALL dbcsr_transposed(qs_ot_env%matrix_buf4, qs_ot_env%matrix_buf3, &
3625 : shallow_data_copy=.FALSE., use_distribution=dist, &
3626 2844 : transpose_distribution=.FALSE.)
3627 : CALL dbcsr_transposed(qs_ot_env%matrix_buf4_im, qs_ot_env%matrix_buf3_im, &
3628 : shallow_data_copy=.FALSE., use_distribution=dist, &
3629 2844 : transpose_distribution=.FALSE.)
3630 : CALL dbcsr_add(qs_ot_env%matrix_buf3, qs_ot_env%matrix_buf4, &
3631 2844 : alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
3632 : CALL dbcsr_add(qs_ot_env%matrix_buf3_im, qs_ot_env%matrix_buf4_im, &
3633 2844 : alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
3634 :
3635 : CALL qs_ot_complex_multiply('N', 'N', qs_ot_env%matrix_sx, qs_ot_env%matrix_sx_im, &
3636 : qs_ot_env%matrix_buf3, qs_ot_env%matrix_buf3_im, &
3637 : qs_ot_env%matrix_buf_nk, qs_ot_env%matrix_buf_nk_im, &
3638 2844 : qs_ot_env%matrix_tmp_nk)
3639 : CALL dbcsr_add(qs_ot_env%matrix_gx, qs_ot_env%matrix_buf_nk, &
3640 2844 : alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
3641 : CALL dbcsr_add(qs_ot_env%matrix_gx_im, qs_ot_env%matrix_buf_nk_im, &
3642 2844 : alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
3643 :
3644 : ! Preconditioner-aware projection onto the complex STRICT tangent space.
3645 2844 : CALL qs_ot_prepare_complex_tangent_metric(qs_ot_env)
3646 2844 : IF (ASSOCIATED(qs_ot_env%preconditioner)) THEN
3647 2818 : target_re => qs_ot_env%matrix_psc0
3648 2818 : target_im => qs_ot_env%matrix_psc0_im
3649 : ELSE
3650 26 : target_re => qs_ot_env%matrix_sc0
3651 26 : target_im => qs_ot_env%matrix_sc0_im
3652 : END IF
3653 :
3654 : CALL qs_ot_complex_multiply('C', 'N', target_re, target_im, &
3655 : qs_ot_env%matrix_gx, qs_ot_env%matrix_gx_im, &
3656 : qs_ot_env%matrix_buf1_ortho, qs_ot_env%matrix_buf1_ortho_im, &
3657 2844 : qs_ot_env%matrix_tmp_ortho)
3658 : CALL qs_ot_complex_multiply('N', 'N', qs_ot_env%matrix_os, qs_ot_env%matrix_os_im, &
3659 : qs_ot_env%matrix_buf1_ortho, qs_ot_env%matrix_buf1_ortho_im, &
3660 : qs_ot_env%matrix_buf2_ortho, qs_ot_env%matrix_buf2_ortho_im, &
3661 2844 : qs_ot_env%matrix_tmp_ortho)
3662 : CALL qs_ot_complex_multiply('N', 'N', qs_ot_env%matrix_sc0, qs_ot_env%matrix_sc0_im, &
3663 : qs_ot_env%matrix_buf2_ortho, qs_ot_env%matrix_buf2_ortho_im, &
3664 : qs_ot_env%matrix_buf_nk, qs_ot_env%matrix_buf_nk_im, &
3665 2844 : qs_ot_env%matrix_tmp_nk)
3666 : CALL dbcsr_add(qs_ot_env%matrix_gx, qs_ot_env%matrix_buf_nk, &
3667 2844 : alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
3668 : CALL dbcsr_add(qs_ot_env%matrix_gx_im, qs_ot_env%matrix_buf_nk_im, &
3669 2844 : alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
3670 :
3671 2844 : IF (qs_ot_env%settings%do_rotation) THEN
3672 648 : CALL dbcsr_release(hc_rot_re)
3673 648 : CALL dbcsr_release(hc_rot_im)
3674 648 : CALL dbcsr_release(tmp_nk)
3675 : END IF
3676 :
3677 2844 : CALL timestop(handle)
3678 :
3679 2844 : END SUBROUTINE qs_ot_get_derivative_complex
3680 :
3681 : ! **************************************************************************************************
3682 : !> \brief ...
3683 : !> \param matrix_hc ...
3684 : !> \param matrix_x ...
3685 : !> \param matrix_sx ...
3686 : !> \param matrix_gx ...
3687 : !> \param qs_ot_env ...
3688 : ! **************************************************************************************************
3689 106338 : SUBROUTINE qs_ot_get_derivative_diag(matrix_hc, matrix_x, matrix_sx, matrix_gx, &
3690 : qs_ot_env)
3691 :
3692 : TYPE(dbcsr_type), POINTER :: matrix_hc, matrix_x, matrix_sx, matrix_gx
3693 : TYPE(qs_ot_type) :: qs_ot_env
3694 :
3695 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_ot_get_derivative_diag'
3696 : REAL(KIND=dp), PARAMETER :: rone = 1.0_dp, rzero = 0.0_dp
3697 :
3698 : INTEGER :: handle, k, n
3699 : TYPE(dbcsr_distribution_type) :: dist
3700 :
3701 35446 : CALL timeset(routineN, handle)
3702 :
3703 35446 : CALL dbcsr_get_info(matrix_x, nfullrows_total=n, nfullcols_total=k)
3704 :
3705 : ! go for the derivative now
3706 : ! this de/dc*(dX/dx)*sinp
3707 35446 : CALL dbcsr_multiply('N', 'N', rone, matrix_hc, qs_ot_env%matrix_sinp, rzero, matrix_gx)
3708 : ! overlap hc*x
3709 35446 : CALL dbcsr_multiply('T', 'N', rone, matrix_hc, matrix_x, rzero, qs_ot_env%matrix_buf2)
3710 : ! get it in the basis of the eigenvectors
3711 : CALL dbcsr_multiply('N', 'N', rone, qs_ot_env%matrix_buf2, qs_ot_env%matrix_r, &
3712 35446 : rzero, qs_ot_env%matrix_buf1)
3713 : CALL dbcsr_multiply('T', 'N', rone, qs_ot_env%matrix_r, qs_ot_env%matrix_buf1, &
3714 35446 : rzero, qs_ot_env%matrix_buf2)
3715 :
3716 : ! get the schur product of O_uv*B_uv
3717 : CALL dbcsr_hadamard_product(qs_ot_env%matrix_buf2, qs_ot_env%matrix_sinp_b, &
3718 35446 : qs_ot_env%matrix_buf3)
3719 :
3720 : ! overlap hc*c0
3721 : CALL dbcsr_multiply('T', 'N', rone, matrix_hc, qs_ot_env%matrix_c0, rzero, &
3722 35446 : qs_ot_env%matrix_buf2)
3723 : ! get it in the basis of the eigenvectors
3724 : CALL dbcsr_multiply('N', 'N', rone, qs_ot_env%matrix_buf2, qs_ot_env%matrix_r, &
3725 35446 : rzero, qs_ot_env%matrix_buf1)
3726 : CALL dbcsr_multiply('T', 'N', rone, qs_ot_env%matrix_r, qs_ot_env%matrix_buf1, &
3727 35446 : rzero, qs_ot_env%matrix_buf2)
3728 :
3729 : CALL dbcsr_hadamard_product(qs_ot_env%matrix_buf2, qs_ot_env%matrix_cosp_b, &
3730 35446 : qs_ot_env%matrix_buf4)
3731 :
3732 : ! add the two bs and compute b+b^T
3733 : CALL dbcsr_add(qs_ot_env%matrix_buf3, qs_ot_env%matrix_buf4, &
3734 35446 : alpha_scalar=rone, beta_scalar=rone)
3735 :
3736 : ! get the b in the eigenvector basis
3737 : CALL dbcsr_multiply('N', 'T', rone, qs_ot_env%matrix_buf3, qs_ot_env%matrix_r, &
3738 35446 : rzero, qs_ot_env%matrix_buf1)
3739 : CALL dbcsr_multiply('N', 'N', rone, qs_ot_env%matrix_r, qs_ot_env%matrix_buf1, &
3740 35446 : rzero, qs_ot_env%matrix_buf3)
3741 35446 : CALL dbcsr_get_info(qs_ot_env%matrix_buf3, distribution=dist)
3742 : CALL dbcsr_transposed(qs_ot_env%matrix_buf1, qs_ot_env%matrix_buf3, &
3743 : shallow_data_copy=.FALSE., use_distribution=dist, &
3744 35446 : transpose_distribution=.FALSE.)
3745 : CALL dbcsr_add(qs_ot_env%matrix_buf3, qs_ot_env%matrix_buf1, &
3746 35446 : alpha_scalar=rone, beta_scalar=rone)
3747 :
3748 : ! and add to the derivative
3749 : CALL dbcsr_multiply('N', 'N', rone, matrix_sx, qs_ot_env%matrix_buf3, &
3750 35446 : rone, matrix_gx)
3751 35446 : CALL timestop(handle)
3752 :
3753 35446 : END SUBROUTINE qs_ot_get_derivative_diag
3754 :
3755 : ! **************************************************************************************************
3756 : !> \brief compute the derivative of the taylor expansion below
3757 : !> \param matrix_hc ...
3758 : !> \param matrix_x ...
3759 : !> \param matrix_sx ...
3760 : !> \param matrix_gx ...
3761 : !> \param qs_ot_env ...
3762 : ! **************************************************************************************************
3763 176460 : SUBROUTINE qs_ot_get_derivative_taylor(matrix_hc, matrix_x, matrix_sx, matrix_gx, &
3764 : qs_ot_env)
3765 :
3766 : TYPE(dbcsr_type), POINTER :: matrix_hc, matrix_x, matrix_sx, matrix_gx
3767 : TYPE(qs_ot_type) :: qs_ot_env
3768 :
3769 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_ot_get_derivative_taylor'
3770 : REAL(KIND=dp), PARAMETER :: rone = 1.0_dp, rzero = 0.0_dp
3771 :
3772 : INTEGER :: handle, i, k, n
3773 : REAL(KIND=dp) :: cosfactor, sinfactor
3774 : TYPE(dbcsr_distribution_type) :: dist
3775 : TYPE(dbcsr_type), POINTER :: matrix_left, matrix_right
3776 :
3777 44115 : CALL timeset(routineN, handle)
3778 :
3779 44115 : CALL dbcsr_get_info(matrix_x, nfullrows_total=n, nfullcols_total=k)
3780 :
3781 : ! go for the derivative now
3782 : ! this de/dc*(dX/dx)*sinp i.e. zeroth order
3783 44115 : CALL dbcsr_multiply('N', 'N', rone, matrix_hc, qs_ot_env%matrix_sinp, rzero, matrix_gx)
3784 :
3785 44115 : IF (qs_ot_env%taylor_order <= 0) THEN
3786 0 : CALL timestop(handle)
3787 0 : RETURN
3788 : END IF
3789 :
3790 : ! we store the matrix that will multiply sx in matrix_r
3791 44115 : CALL dbcsr_set(qs_ot_env%matrix_r, rzero)
3792 :
3793 : ! just better names for matrix_cosp_b and matrix_sinp_b (they are buffer space here)
3794 44115 : matrix_left => qs_ot_env%matrix_cosp_b
3795 44115 : matrix_right => qs_ot_env%matrix_sinp_b
3796 :
3797 : ! overlap hc*x and add its transpose to matrix_left
3798 44115 : CALL dbcsr_multiply('T', 'N', rone, matrix_hc, matrix_x, rzero, matrix_left)
3799 44115 : CALL dbcsr_get_info(matrix_left, distribution=dist)
3800 : CALL dbcsr_transposed(qs_ot_env%matrix_buf1, matrix_left, &
3801 : shallow_data_copy=.FALSE., use_distribution=dist, &
3802 44115 : transpose_distribution=.FALSE.)
3803 : CALL dbcsr_add(matrix_left, qs_ot_env%matrix_buf1, &
3804 44115 : alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
3805 44115 : CALL dbcsr_copy(matrix_right, matrix_left)
3806 :
3807 : ! first order
3808 44115 : sinfactor = -1.0_dp/(2.0_dp*3.0_dp)
3809 : CALL dbcsr_add(qs_ot_env%matrix_r, matrix_left, &
3810 44115 : alpha_scalar=1.0_dp, beta_scalar=sinfactor)
3811 :
3812 : ! M
3813 : ! OM+MO
3814 : ! OOM+OMO+MOO
3815 : ! ...
3816 82655 : DO i = 2, qs_ot_env%taylor_order
3817 38540 : sinfactor = sinfactor*(-1.0_dp)/REAL(2*i*(2*i + 1), KIND=dp)
3818 38540 : CALL dbcsr_multiply('N', 'N', rone, qs_ot_env%matrix_p, matrix_left, rzero, qs_ot_env%matrix_buf1)
3819 38540 : CALL dbcsr_multiply('N', 'N', rone, matrix_right, qs_ot_env%matrix_p, rzero, matrix_left)
3820 38540 : CALL dbcsr_copy(matrix_right, matrix_left)
3821 : CALL dbcsr_add(matrix_left, qs_ot_env%matrix_buf1, &
3822 38540 : 1.0_dp, 1.0_dp)
3823 : CALL dbcsr_add(qs_ot_env%matrix_r, matrix_left, &
3824 82655 : alpha_scalar=1.0_dp, beta_scalar=sinfactor)
3825 : END DO
3826 :
3827 : ! overlap hc*c0 and add its transpose to matrix_left
3828 44115 : CALL dbcsr_multiply('T', 'N', rone, matrix_hc, qs_ot_env%matrix_c0, rzero, matrix_left)
3829 44115 : CALL dbcsr_get_info(matrix_left, distribution=dist)
3830 : CALL dbcsr_transposed(qs_ot_env%matrix_buf1, matrix_left, &
3831 : shallow_data_copy=.FALSE., use_distribution=dist, &
3832 44115 : transpose_distribution=.FALSE.)
3833 44115 : CALL dbcsr_add(matrix_left, qs_ot_env%matrix_buf1, 1.0_dp, 1.0_dp)
3834 44115 : CALL dbcsr_copy(matrix_right, matrix_left)
3835 :
3836 : ! first order
3837 44115 : cosfactor = -1.0_dp/(1.0_dp*2.0_dp)
3838 : CALL dbcsr_add(qs_ot_env%matrix_r, matrix_left, &
3839 44115 : alpha_scalar=1.0_dp, beta_scalar=cosfactor)
3840 :
3841 : ! M
3842 : ! OM+MO
3843 : ! OOM+OMO+MOO
3844 : ! ...
3845 82655 : DO i = 2, qs_ot_env%taylor_order
3846 38540 : cosfactor = cosfactor*(-1.0_dp)/REAL(2*i*(2*i - 1), KIND=dp)
3847 38540 : CALL dbcsr_multiply('N', 'N', rone, qs_ot_env%matrix_p, matrix_left, rzero, qs_ot_env%matrix_buf1)
3848 38540 : CALL dbcsr_multiply('N', 'N', rone, matrix_right, qs_ot_env%matrix_p, rzero, matrix_left)
3849 38540 : CALL dbcsr_copy(matrix_right, matrix_left)
3850 38540 : CALL dbcsr_add(matrix_left, qs_ot_env%matrix_buf1, 1.0_dp, 1.0_dp)
3851 : CALL dbcsr_add(qs_ot_env%matrix_r, matrix_left, &
3852 82655 : alpha_scalar=1.0_dp, beta_scalar=cosfactor)
3853 : END DO
3854 :
3855 : ! and add to the derivative
3856 44115 : CALL dbcsr_multiply('N', 'N', rone, matrix_sx, qs_ot_env%matrix_r, rone, matrix_gx)
3857 :
3858 44115 : CALL timestop(handle)
3859 :
3860 44115 : END SUBROUTINE qs_ot_get_derivative_taylor
3861 :
3862 : ! *************************************************************************************************
3863 : !> \brief computes a taylor expansion.
3864 : !> \param qs_ot_env ...
3865 : ! **************************************************************************************************
3866 96287 : SUBROUTINE qs_ot_p2m_taylor(qs_ot_env)
3867 : TYPE(qs_ot_type) :: qs_ot_env
3868 :
3869 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_ot_p2m_taylor'
3870 : REAL(KIND=dp), PARAMETER :: rone = 1.0_dp, rzero = 0.0_dp
3871 :
3872 : INTEGER :: handle, i, k
3873 : REAL(KIND=dp) :: cosfactor, sinfactor
3874 :
3875 59916 : CALL timeset(routineN, handle)
3876 :
3877 : ! zeroth order
3878 59916 : CALL dbcsr_set(qs_ot_env%matrix_cosp, rzero)
3879 59916 : CALL dbcsr_set(qs_ot_env%matrix_sinp, rzero)
3880 59916 : CALL dbcsr_add_on_diag(qs_ot_env%matrix_cosp, rone)
3881 59916 : CALL dbcsr_add_on_diag(qs_ot_env%matrix_sinp, rone)
3882 :
3883 59916 : IF (qs_ot_env%taylor_order <= 0) THEN
3884 0 : CALL timestop(handle)
3885 23545 : RETURN
3886 : END IF
3887 :
3888 : ! first order
3889 59916 : cosfactor = -1.0_dp/(1.0_dp*2.0_dp)
3890 59916 : sinfactor = -1.0_dp/(2.0_dp*3.0_dp)
3891 59916 : CALL dbcsr_add(qs_ot_env%matrix_cosp, qs_ot_env%matrix_p, alpha_scalar=1.0_dp, beta_scalar=cosfactor)
3892 59916 : CALL dbcsr_add(qs_ot_env%matrix_sinp, qs_ot_env%matrix_p, alpha_scalar=1.0_dp, beta_scalar=sinfactor)
3893 59916 : IF (qs_ot_env%taylor_order <= 1) THEN
3894 23545 : CALL timestop(handle)
3895 23545 : RETURN
3896 : END IF
3897 :
3898 : ! other orders
3899 36371 : CALL dbcsr_get_info(qs_ot_env%matrix_p, nfullrows_total=k)
3900 36371 : CALL dbcsr_copy(qs_ot_env%matrix_r, qs_ot_env%matrix_p)
3901 :
3902 91750 : DO i = 2, qs_ot_env%taylor_order
3903 : ! new power of p
3904 : CALL dbcsr_multiply('N', 'N', rone, qs_ot_env%matrix_p, qs_ot_env%matrix_r, &
3905 55379 : rzero, qs_ot_env%matrix_buf1)
3906 55379 : CALL dbcsr_copy(qs_ot_env%matrix_r, qs_ot_env%matrix_buf1)
3907 : ! add to the taylor expansion so far
3908 55379 : cosfactor = cosfactor*(-1.0_dp)/REAL(2*i*(2*i - 1), KIND=dp)
3909 55379 : sinfactor = sinfactor*(-1.0_dp)/REAL(2*i*(2*i + 1), KIND=dp)
3910 : CALL dbcsr_add(qs_ot_env%matrix_cosp, qs_ot_env%matrix_r, &
3911 55379 : alpha_scalar=1.0_dp, beta_scalar=cosfactor)
3912 : CALL dbcsr_add(qs_ot_env%matrix_sinp, qs_ot_env%matrix_r, &
3913 91750 : alpha_scalar=1.0_dp, beta_scalar=sinfactor)
3914 : END DO
3915 :
3916 36371 : CALL timestop(handle)
3917 :
3918 : END SUBROUTINE qs_ot_p2m_taylor
3919 :
3920 : ! **************************************************************************************************
3921 : !> \brief given p, computes - eigenstuff (matrix_r,evals)
3922 : !> - cos(p^0.5),p^(-0.5)*sin(p^0.5)
3923 : !> - the real b matrices, needed for the derivatives of these guys
3924 : !> cosp_b_ij=(1/(2pii) * int(cos(z^1/2)/((z-eval(i))*(z-eval(j))))
3925 : !> sinp_b_ij=(1/(2pii) * int(z^(-1/2)*sin(z^1/2)/((z-eval(i))*(z-eval(j))))
3926 : !> \param qs_ot_env ...
3927 : ! **************************************************************************************************
3928 204396 : SUBROUTINE qs_ot_p2m_diag(qs_ot_env)
3929 :
3930 : TYPE(qs_ot_type) :: qs_ot_env
3931 :
3932 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_ot_p2m_diag'
3933 : REAL(KIND=dp), PARAMETER :: rone = 1.0_dp, rzero = 0.0_dp
3934 :
3935 : INTEGER :: col, col_offset, col_size, handle, i, j, &
3936 : k, row, row_offset, row_size
3937 51099 : REAL(dp), DIMENSION(:, :), POINTER :: block
3938 : REAL(KIND=dp) :: a, b
3939 : TYPE(dbcsr_iterator_type) :: iter
3940 :
3941 51099 : CALL timeset(routineN, handle)
3942 :
3943 51099 : CALL dbcsr_get_info(qs_ot_env%matrix_p, nfullrows_total=k)
3944 51099 : CALL dbcsr_copy(qs_ot_env%matrix_buf1, qs_ot_env%matrix_p)
3945 : CALL cp_dbcsr_syevd(qs_ot_env%matrix_buf1, qs_ot_env%matrix_r, qs_ot_env%evals, &
3946 51099 : qs_ot_env%para_env, qs_ot_env%blacs_env)
3947 529984 : DO i = 1, k
3948 529984 : qs_ot_env%evals(i) = MAX(0.0_dp, qs_ot_env%evals(i))
3949 : END DO
3950 :
3951 51099 : !$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(i) SHARED(k,qs_ot_env)
3952 : DO i = 1, k
3953 : qs_ot_env%dum(i) = COS(SQRT(qs_ot_env%evals(i)))
3954 : END DO
3955 51099 : CALL dbcsr_copy(qs_ot_env%matrix_buf1, qs_ot_env%matrix_r)
3956 51099 : CALL dbcsr_scale_by_vector(qs_ot_env%matrix_buf1, alpha=qs_ot_env%dum, side='right')
3957 : CALL dbcsr_multiply('N', 'T', rone, qs_ot_env%matrix_r, qs_ot_env%matrix_buf1, &
3958 51099 : rzero, qs_ot_env%matrix_cosp)
3959 :
3960 51099 : !$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(i) SHARED(k,qs_ot_env)
3961 : DO i = 1, k
3962 : qs_ot_env%dum(i) = qs_ot_sinc(SQRT(qs_ot_env%evals(i)))
3963 : END DO
3964 51099 : CALL dbcsr_copy(qs_ot_env%matrix_buf1, qs_ot_env%matrix_r)
3965 51099 : CALL dbcsr_scale_by_vector(qs_ot_env%matrix_buf1, alpha=qs_ot_env%dum, side='right')
3966 : CALL dbcsr_multiply('N', 'T', rone, qs_ot_env%matrix_r, qs_ot_env%matrix_buf1, &
3967 51099 : rzero, qs_ot_env%matrix_sinp)
3968 :
3969 51099 : CALL dbcsr_copy(qs_ot_env%matrix_cosp_b, qs_ot_env%matrix_cosp)
3970 51099 : CALL dbcsr_iterator_start(iter, qs_ot_env%matrix_cosp_b)
3971 87340 : DO WHILE (dbcsr_iterator_blocks_left(iter))
3972 : CALL dbcsr_iterator_next_block(iter, row, col, block, &
3973 : row_size=row_size, col_size=col_size, &
3974 36241 : row_offset=row_offset, col_offset=col_offset)
3975 549087 : DO j = 1, col_size
3976 10678828 : DO i = 1, row_size
3977 : a = (SQRT(qs_ot_env%evals(row_offset + i - 1)) &
3978 10180840 : - SQRT(qs_ot_env%evals(col_offset + j - 1)))/2.0_dp
3979 : b = (SQRT(qs_ot_env%evals(row_offset + i - 1)) &
3980 10180840 : + SQRT(qs_ot_env%evals(col_offset + j - 1)))/2.0_dp
3981 10642587 : block(i, j) = -0.5_dp*qs_ot_sinc(a)*qs_ot_sinc(b)
3982 : END DO
3983 : END DO
3984 : END DO
3985 51099 : CALL dbcsr_iterator_stop(iter)
3986 :
3987 51099 : CALL dbcsr_copy(qs_ot_env%matrix_sinp_b, qs_ot_env%matrix_sinp)
3988 51099 : CALL dbcsr_iterator_start(iter, qs_ot_env%matrix_sinp_b)
3989 87340 : DO WHILE (dbcsr_iterator_blocks_left(iter))
3990 : CALL dbcsr_iterator_next_block(iter, row, col, block, &
3991 : row_size=row_size, col_size=col_size, &
3992 36241 : row_offset=row_offset, col_offset=col_offset)
3993 549087 : DO j = 1, col_size
3994 10678828 : DO i = 1, row_size
3995 10180840 : a = SQRT(qs_ot_env%evals(row_offset + i - 1))
3996 10180840 : b = SQRT(qs_ot_env%evals(col_offset + j - 1))
3997 10642587 : block(i, j) = qs_ot_sincf(a, b)
3998 : END DO
3999 : END DO
4000 : END DO
4001 51099 : CALL dbcsr_iterator_stop(iter)
4002 :
4003 51099 : CALL timestop(handle)
4004 :
4005 51099 : END SUBROUTINE qs_ot_p2m_diag
4006 :
4007 : ! **************************************************************************************************
4008 : !> \brief diagonalize Hermitian P and build cos(sqrt(P)), sinc(sqrt(P)), and Frechet kernels
4009 : !> \param qs_ot_env complex STRICT channel state
4010 : ! **************************************************************************************************
4011 14208 : SUBROUTINE qs_ot_p2m_diag_complex(qs_ot_env)
4012 : TYPE(qs_ot_type) :: qs_ot_env
4013 :
4014 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_ot_p2m_diag_complex'
4015 :
4016 : INTEGER :: col, col_offset, col_size, handle, i, j, &
4017 : k, row, row_offset, row_size
4018 3552 : REAL(dp), DIMENSION(:, :), POINTER :: block
4019 : REAL(KIND=dp) :: a, b
4020 : TYPE(dbcsr_iterator_type) :: iter
4021 :
4022 3552 : CALL timeset(routineN, handle)
4023 :
4024 3552 : CALL dbcsr_get_info(qs_ot_env%matrix_p, nfullrows_total=k)
4025 : CALL cp_dbcsr_heevd(matrix_re=qs_ot_env%matrix_p, matrix_im=qs_ot_env%matrix_p_im, &
4026 : eigenvectors_re=qs_ot_env%matrix_r, eigenvectors_im=qs_ot_env%matrix_r_im, &
4027 : eigenvalues=qs_ot_env%evals, para_env=qs_ot_env%para_env, &
4028 3552 : blacs_env=qs_ot_env%blacs_env)
4029 22360 : DO i = 1, k
4030 22360 : qs_ot_env%evals(i) = MAX(0.0_dp, qs_ot_env%evals(i))
4031 : END DO
4032 :
4033 22360 : DO i = 1, k
4034 22360 : qs_ot_env%dum(i) = COS(SQRT(qs_ot_env%evals(i)))
4035 : END DO
4036 3552 : CALL dbcsr_copy(qs_ot_env%matrix_buf1, qs_ot_env%matrix_r)
4037 3552 : CALL dbcsr_copy(qs_ot_env%matrix_buf1_im, qs_ot_env%matrix_r_im)
4038 3552 : CALL dbcsr_scale_by_vector(qs_ot_env%matrix_buf1, alpha=qs_ot_env%dum, side='right')
4039 3552 : CALL dbcsr_scale_by_vector(qs_ot_env%matrix_buf1_im, alpha=qs_ot_env%dum, side='right')
4040 : CALL qs_ot_complex_multiply('N', 'C', qs_ot_env%matrix_buf1, qs_ot_env%matrix_buf1_im, &
4041 : qs_ot_env%matrix_r, qs_ot_env%matrix_r_im, &
4042 : qs_ot_env%matrix_cosp, qs_ot_env%matrix_cosp_im, &
4043 3552 : qs_ot_env%matrix_buf2)
4044 :
4045 22360 : DO i = 1, k
4046 22360 : qs_ot_env%dum(i) = qs_ot_sinc(SQRT(qs_ot_env%evals(i)))
4047 : END DO
4048 3552 : CALL dbcsr_copy(qs_ot_env%matrix_buf1, qs_ot_env%matrix_r)
4049 3552 : CALL dbcsr_copy(qs_ot_env%matrix_buf1_im, qs_ot_env%matrix_r_im)
4050 3552 : CALL dbcsr_scale_by_vector(qs_ot_env%matrix_buf1, alpha=qs_ot_env%dum, side='right')
4051 3552 : CALL dbcsr_scale_by_vector(qs_ot_env%matrix_buf1_im, alpha=qs_ot_env%dum, side='right')
4052 : CALL qs_ot_complex_multiply('N', 'C', qs_ot_env%matrix_buf1, qs_ot_env%matrix_buf1_im, &
4053 : qs_ot_env%matrix_r, qs_ot_env%matrix_r_im, &
4054 : qs_ot_env%matrix_sinp, qs_ot_env%matrix_sinp_im, &
4055 3552 : qs_ot_env%matrix_buf2)
4056 :
4057 3552 : CALL dbcsr_copy(qs_ot_env%matrix_cosp_b, qs_ot_env%matrix_cosp)
4058 3552 : CALL dbcsr_iterator_start(iter, qs_ot_env%matrix_cosp_b)
4059 6335 : DO WHILE (dbcsr_iterator_blocks_left(iter))
4060 : CALL dbcsr_iterator_next_block(iter, row, col, block, &
4061 : row_size=row_size, col_size=col_size, &
4062 2783 : row_offset=row_offset, col_offset=col_offset)
4063 21169 : DO j = 1, col_size
4064 155041 : DO i = 1, row_size
4065 : a = (SQRT(qs_ot_env%evals(row_offset + i - 1)) &
4066 137424 : - SQRT(qs_ot_env%evals(col_offset + j - 1)))/2.0_dp
4067 : b = (SQRT(qs_ot_env%evals(row_offset + i - 1)) &
4068 137424 : + SQRT(qs_ot_env%evals(col_offset + j - 1)))/2.0_dp
4069 152258 : block(i, j) = -0.5_dp*qs_ot_sinc(a)*qs_ot_sinc(b)
4070 : END DO
4071 : END DO
4072 : END DO
4073 3552 : CALL dbcsr_iterator_stop(iter)
4074 :
4075 3552 : CALL dbcsr_copy(qs_ot_env%matrix_sinp_b, qs_ot_env%matrix_sinp)
4076 3552 : CALL dbcsr_iterator_start(iter, qs_ot_env%matrix_sinp_b)
4077 6335 : DO WHILE (dbcsr_iterator_blocks_left(iter))
4078 : CALL dbcsr_iterator_next_block(iter, row, col, block, &
4079 : row_size=row_size, col_size=col_size, &
4080 2783 : row_offset=row_offset, col_offset=col_offset)
4081 21169 : DO j = 1, col_size
4082 155041 : DO i = 1, row_size
4083 137424 : a = SQRT(qs_ot_env%evals(row_offset + i - 1))
4084 137424 : b = SQRT(qs_ot_env%evals(col_offset + j - 1))
4085 152258 : block(i, j) = qs_ot_sincf(a, b)
4086 : END DO
4087 : END DO
4088 : END DO
4089 3552 : CALL dbcsr_iterator_stop(iter)
4090 :
4091 3552 : CALL timestop(handle)
4092 :
4093 3552 : END SUBROUTINE qs_ot_p2m_diag_complex
4094 :
4095 : ! **************************************************************************************************
4096 : !> \brief computes sin(x)/x for all values of the argument
4097 : !> \param x ...
4098 : !> \return ...
4099 : ! **************************************************************************************************
4100 28856597 : FUNCTION qs_ot_sinc(x)
4101 :
4102 : REAL(KIND=dp), INTENT(IN) :: x
4103 : REAL(KIND=dp) :: qs_ot_sinc
4104 :
4105 : REAL(KIND=dp), PARAMETER :: q1 = 1.0_dp, q2 = -q1/(2.0_dp*3.0_dp), q3 = -q2/(4.0_dp*5.0_dp), &
4106 : q4 = -q3/(6.0_dp*7.0_dp), q5 = -q4/(8.0_dp*9.0_dp), q6 = -q5/(10.0_dp*11.0_dp), &
4107 : q7 = -q6/(12.0_dp*13.0_dp), q8 = -q7/(14.0_dp*15.0_dp), q9 = -q8/(16.0_dp*17.0_dp), &
4108 : q10 = -q9/(18.0_dp*19.0_dp)
4109 :
4110 : REAL(KIND=dp) :: y
4111 :
4112 28856597 : IF (ABS(x) > 0.5_dp) THEN
4113 8686659 : qs_ot_sinc = SIN(x)/x
4114 : ELSE
4115 20169938 : y = x*x
4116 20169938 : qs_ot_sinc = q1 + y*(q2 + y*(q3 + y*(q4 + y*(q5 + y*(q6 + y*(q7 + y*(q8 + y*(q9 + y*(q10)))))))))
4117 : END IF
4118 28856597 : END FUNCTION qs_ot_sinc
4119 :
4120 : ! **************************************************************************************************
4121 : !> \brief computes (1/(x^2-y^2))*(sinc(x)-sinc(y)) for all positive values of the arguments
4122 : !> \param xa ...
4123 : !> \param ya ...
4124 : !> \return ...
4125 : ! **************************************************************************************************
4126 10318264 : FUNCTION qs_ot_sincf(xa, ya)
4127 :
4128 : REAL(KIND=dp), INTENT(IN) :: xa, ya
4129 : REAL(KIND=dp) :: qs_ot_sincf
4130 :
4131 : INTEGER :: i
4132 : REAL(KIND=dp) :: a, b, rs, sf, x, xs, y, ybx, ybxs
4133 :
4134 : ! this is currently a limit of the routine, could be removed rather easily
4135 10318264 : IF (xa < 0) CPABORT("x is negative")
4136 10318264 : IF (ya < 0) CPABORT("y is negative")
4137 :
4138 10318264 : IF (xa < ya) THEN
4139 4929902 : x = ya
4140 4929902 : y = xa
4141 : ELSE
4142 5388362 : x = xa
4143 5388362 : y = ya
4144 : END IF
4145 :
4146 10318264 : IF (x < 0.5_dp) THEN ! use series, keeping in mind that x,y,x+y,x-y can all be zero
4147 :
4148 6457076 : qs_ot_sincf = 0.0_dp
4149 6457076 : IF (x > 0.0_dp) THEN
4150 6234416 : ybx = y/x
4151 : ELSE ! should be irrelevant !?
4152 : ybx = 0.0_dp
4153 : END IF
4154 :
4155 6457076 : sf = -1.0_dp/((1.0_dp + ybx)*6.0_dp)
4156 6457076 : rs = 1.0_dp
4157 6457076 : ybxs = ybx
4158 6457076 : xs = 1.0_dp
4159 :
4160 71027836 : DO i = 1, 10
4161 64570760 : qs_ot_sincf = qs_ot_sincf + sf*rs*xs*(1.0_dp + ybxs)
4162 64570760 : sf = -sf/(REAL((2*i + 2), dp)*REAL((2*i + 3), dp))
4163 64570760 : rs = rs + ybxs
4164 64570760 : ybxs = ybxs*ybx
4165 71027836 : xs = xs*x*x
4166 : END DO
4167 :
4168 : ELSE ! no series expansion
4169 3861188 : IF (x - y > 0.1_dp) THEN ! safe to use the normal form
4170 3574872 : qs_ot_sincf = (qs_ot_sinc(x) - qs_ot_sinc(y))/((x + y)*(x - y))
4171 : ELSE
4172 286316 : a = (x + y)/2.0_dp
4173 286316 : b = (x - y)/2.0_dp ! might be close to zero
4174 : ! y (=(a-b)) can not be close to zero since it is close to x>0.5
4175 286316 : qs_ot_sincf = (qs_ot_sinc(b)*COS(a) - qs_ot_sinc(a)*COS(b))/(2*x*y)
4176 : END IF
4177 : END IF
4178 :
4179 10318264 : END FUNCTION qs_ot_sincf
4180 :
4181 67208 : END MODULE qs_ot
|
