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 2212 : PURE SUBROUTINE qs_ot_fixed_n_energy_gradient( &
154 2212 : rayleigh_energy, energy_coordinate, response_weight, fixed_n_weight_sum, &
155 2212 : 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 2212 : fixed_n_mean = 0.0_dp
165 2212 : IF (ABS(fixed_n_weight_sum) > EPSILON(fixed_n_weight_sum)) THEN
166 2168 : fixed_n_mean = fixed_n_weighted_residual/fixed_n_weight_sum
167 : END IF
168 : gradient(:) = response_weight(:)* &
169 17346 : (fixed_n_mean - (rayleigh_energy(:) - energy_coordinate(:)))
170 :
171 2212 : 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 338 : SUBROUTINE qs_ot_fixed_n_multigroup_schur_block( &
268 676 : rotation_hessian, rayleigh_response, response_weight, response_group, &
269 676 : 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 338 : nenergy = SIZE(response_weight)
281 338 : nrotation = SIZE(rotation_gradient)
282 338 : ngroups = SIZE(coupling_matrix, 2)
283 338 : CPASSERT(ngroups > 0)
284 1014 : CPASSERT(ALL(SHAPE(rotation_hessian) == [nrotation, nrotation]))
285 1014 : CPASSERT(ALL(SHAPE(rayleigh_response) == [nenergy, nrotation]))
286 338 : CPASSERT(SIZE(energy_gradient) == nenergy)
287 338 : CPASSERT(SIZE(response_group) == nenergy)
288 2304 : CPASSERT(ALL(response_group >= 1 .AND. response_group <= ngroups))
289 1014 : CPASSERT(ALL(SHAPE(schur_block) == [nrotation, nrotation]))
290 338 : CPASSERT(SIZE(coupling_matrix, 1) == nrotation)
291 338 : CPASSERT(SIZE(schur_rhs) == nrotation)
292 :
293 505370 : schur_block(:, :) = rotation_hessian(:, :)
294 11942 : coupling_matrix(:, :) = 0.0_dp
295 2304 : DO i = 1, nenergy
296 1966 : weight = response_weight(i)
297 : schur_block(:, :) = schur_block(:, :) - weight* &
298 : SPREAD(rayleigh_response(i, :), DIM=2, NCOPIES=nrotation)* &
299 5497258 : SPREAD(rayleigh_response(i, :), DIM=1, NCOPIES=nrotation)
300 1966 : group = response_group(i)
301 : coupling_matrix(:, group) = coupling_matrix(:, group) + &
302 85404 : weight*rayleigh_response(i, :)
303 : END DO
304 1010740 : schur_block(:, :) = 0.5_dp*(schur_block + TRANSPOSE(schur_block))
305 : schur_rhs(:) = rotation_gradient + &
306 10282 : MATMUL(TRANSPOSE(rayleigh_response), energy_gradient)
307 :
308 338 : 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 336 : 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 336 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eigenvalues
333 336 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: eigenvectors, work
334 :
335 336 : n = SIZE(matrix, 1)
336 1008 : CPASSERT(ALL(SHAPE(matrix) == [n, n]))
337 336 : CPASSERT(SIZE(rhs, 1) == n)
338 1008 : CPASSERT(ALL(SHAPE(solution) == SHAPE(rhs)))
339 527046 : valid = ALL(matrix == matrix) .AND. ALL(rhs == rhs)
340 336 : IF (.NOT. valid) THEN
341 0 : solution(:, :) = 0.0_dp
342 : RETURN
343 : END IF
344 :
345 336 : relative_floor_eff = SQRT(EPSILON(1.0_dp))
346 336 : IF (PRESENT(relative_floor)) relative_floor_eff = MAX(relative_floor, relative_floor_eff)
347 3024 : ALLOCATE (eigenvalues(n), eigenvectors(n, n), work(n, SIZE(rhs, 2)))
348 505224 : eigenvectors(:, :) = 0.5_dp*(matrix + TRANSPOSE(matrix))
349 336 : CALL diamat_all(eigenvectors, eigenvalues)
350 9924 : scale = MAX(1.0_dp, MAXVAL(ABS(eigenvalues)))
351 336 : eigenvalue_floor = relative_floor_eff*scale
352 2246760 : work(:, :) = MATMUL(TRANSPOSE(eigenvectors), rhs)
353 336 : nresolved = 0
354 9924 : DO i = 1, n
355 9924 : IF (ABS(eigenvalues(i)) > eigenvalue_floor) THEN
356 15711 : work(i, :) = work(i, :)/ABS(eigenvalues(i))
357 5013 : nresolved = nresolved + 1
358 : ELSE
359 14649 : work(i, :) = 0.0_dp
360 : END IF
361 : END DO
362 2259372 : solution(:, :) = MATMUL(eigenvectors, work)
363 21822 : valid = nresolved > 0 .AND. ALL(solution == solution)
364 708 : IF (.NOT. valid) solution(:, :) = 0.0_dp
365 336 : DEALLOCATE (eigenvalues, eigenvectors, work)
366 :
367 336 : 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 174 : SUBROUTINE qs_ot_density_secant_projected_hessian( &
590 174 : 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 348 : REAL(KIND=dp), DIMENSION(SIZE(density_overlap)) :: projected_density_overlap, &
603 174 : projected_response_overlap
604 :
605 174 : nmode = SIZE(density_overlap)
606 174 : CPASSERT(SIZE(response_overlap) == nmode)
607 522 : CPASSERT(ALL(SHAPE(correction) == [nmode, nmode]))
608 :
609 1218 : correction(:, :) = 0.0_dp
610 174 : valid = .FALSE.
611 174 : IF (nmode <= 0 .OR. density_norm_sq <= TINY(1.0_dp)) RETURN
612 174 : IF (density_norm_sq /= density_norm_sq .OR. response_work /= response_work) RETURN
613 174 : IF (ABS(density_norm_sq) > HUGE(1.0_dp) .OR. ABS(response_work) > HUGE(1.0_dp)) RETURN
614 1044 : IF (ANY(density_overlap /= density_overlap) .OR. ANY(response_overlap /= response_overlap)) RETURN
615 1044 : IF (ANY(ABS(density_overlap) > HUGE(1.0_dp)) .OR. &
616 : ANY(ABS(response_overlap) > HUGE(1.0_dp))) RETURN
617 174 : IF (PRESENT(secant_mode) .NEQV. PRESENT(secant_position)) RETURN
618 :
619 522 : projected_density_overlap(:) = density_overlap
620 522 : projected_response_overlap(:) = response_overlap
621 174 : 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 172 : inverse_density_norm = 1.0_dp/density_norm_sq
633 516 : DO j = 1, nmode
634 1204 : 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 1032 : projected_density_overlap(i)*projected_density_overlap(j))
640 : END DO
641 : END DO
642 2236 : correction(:, :) = 0.5_dp*(correction + TRANSPOSE(correction))
643 2408 : valid = ALL(correction == correction) .AND. ALL(ABS(correction) <= HUGE(1.0_dp))
644 172 : 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 446 : SUBROUTINE qs_ot_density_tangent( &
666 446 : rotation_generator, occupation, kpoint_weight, rotation_step, &
667 446 : 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 446 : COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: direction, generator_minus, &
677 446 : 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 446 : n = SIZE(rotation_generator, 1)
682 446 : nrotation = n*(n - 1)
683 446 : CPASSERT(n > 0)
684 1338 : CPASSERT(ALL(SHAPE(rotation_generator) == [n, n]))
685 446 : CPASSERT(SIZE(occupation) == n)
686 446 : CPASSERT(SIZE(rotation_step) == nrotation)
687 446 : CPASSERT(SIZE(weighted_occupation_step) == n)
688 1338 : CPASSERT(ALL(SHAPE(density_tangent) == [n, n]))
689 446 : 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 8028 : rotation_plus(n, n), weighted_rotation(n, n))
694 446 : direction(:, :) = CMPLX(0.0_dp, 0.0_dp, KIND=dp)
695 446 : r = 0
696 2882 : DO i = 1, n - 1
697 11768 : DO j = i + 1, n
698 8886 : r = r + 1
699 8886 : direction(i, j) = CMPLX(rotation_step(r), 0.0_dp, KIND=dp)
700 8886 : direction(j, i) = -direction(i, j)
701 8886 : r = r + 1
702 : direction(i, j) = direction(i, j) + &
703 8886 : CMPLX(0.0_dp, rotation_step(r), KIND=dp)
704 : direction(j, i) = direction(j, i) + &
705 11322 : CMPLX(0.0_dp, rotation_step(r), KIND=dp)
706 : END DO
707 : END DO
708 446 : CPASSERT(r == nrotation)
709 :
710 23982 : step_scale = MAX(1.0_dp, MAXVAL(ABS(direction)))
711 446 : step = 1.0E-5_dp/step_scale
712 446 : IF (PRESENT(difference_step)) step = difference_step/step_scale
713 446 : CPASSERT(step > EPSILON(step))
714 23982 : generator_plus(:, :) = rotation_generator + step*direction
715 23982 : generator_minus(:, :) = rotation_generator - step*direction
716 446 : CALL qs_ot_dense_rotation_state(rotation_generator, rotation)
717 446 : CALL qs_ot_dense_rotation_state(generator_plus, rotation_plus)
718 446 : CALL qs_ot_dense_rotation_state(generator_minus, rotation_minus)
719 :
720 23982 : weighted_rotation(:, :) = rotation_plus
721 3328 : DO j = 1, n
722 23982 : weighted_rotation(:, j) = kpoint_weight*occupation(j)*weighted_rotation(:, j)
723 : END DO
724 185792 : mode_ref(:, :) = MATMUL(weighted_rotation, CONJG(TRANSPOSE(rotation_plus)))
725 23982 : weighted_rotation(:, :) = rotation_minus
726 3328 : DO j = 1, n
727 23982 : weighted_rotation(:, j) = kpoint_weight*occupation(j)*weighted_rotation(:, j)
728 : END DO
729 : mode_ref(:, :) = (mode_ref - &
730 210666 : MATMUL(weighted_rotation, CONJG(TRANSPOSE(rotation_minus))))/(2.0_dp*step)
731 :
732 23982 : weighted_rotation(:, :) = rotation
733 3328 : DO j = 1, n
734 23982 : weighted_rotation(:, j) = weighted_occupation_step(j)*weighted_rotation(:, j)
735 : END DO
736 211112 : mode_ref(:, :) = mode_ref + MATMUL(weighted_rotation, CONJG(TRANSPOSE(rotation)))
737 557822 : density_tangent(:, :) = MATMUL(CONJG(TRANSPOSE(rotation)), MATMUL(mode_ref, rotation))
738 47518 : density_tangent(:, :) = 0.5_dp*(density_tangent + CONJG(TRANSPOSE(density_tangent)))
739 :
740 0 : DEALLOCATE (direction, generator_minus, generator_plus, mode_ref, rotation, &
741 446 : rotation_minus, rotation_plus, weighted_rotation)
742 :
743 446 : 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 260 : SUBROUTINE qs_ot_density_secant_orbital_overlaps( &
765 520 : overlap_start_current, occupation_start, occupation_current, hamiltonian_step_start, &
766 520 : hamiltonian_step_current, density_modes, kpoint_weight, density_norm_sq, response_work, &
767 260 : 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 260 : COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: cross_mode
780 : INTEGER :: i, j, mode, n, nmode
781 : REAL(KIND=dp) :: cross_density, density_scale
782 :
783 260 : n = SIZE(occupation_start)
784 260 : nmode = SIZE(density_modes, 3)
785 260 : CPASSERT(n > 0)
786 260 : CPASSERT(SIZE(occupation_current) == n)
787 780 : CPASSERT(ALL(SHAPE(overlap_start_current) == [n, n]))
788 780 : CPASSERT(ALL(SHAPE(hamiltonian_step_start) == [n, n]))
789 780 : CPASSERT(ALL(SHAPE(hamiltonian_step_current) == [n, n]))
790 260 : CPASSERT(SIZE(density_modes, 1) == n .AND. SIZE(density_modes, 2) == n)
791 260 : CPASSERT(SIZE(density_overlap) == nmode .AND. SIZE(response_overlap) == nmode)
792 :
793 260 : density_norm_sq = 0.0_dp
794 260 : response_work = 0.0_dp
795 780 : density_overlap(:) = 0.0_dp
796 780 : response_overlap(:) = 0.0_dp
797 260 : valid = .FALSE.
798 276 : IF (nmode <= 0 .OR. kpoint_weight <= TINY(1.0_dp)) RETURN
799 :
800 : cross_density = 0.0_dp
801 1926 : DO j = 1, n
802 13792 : DO i = 1, n
803 : cross_density = cross_density + occupation_start(i)*occupation_current(j)* &
804 13532 : 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 3592 : 2.0_dp*cross_density)
810 : response_work = 0.0_dp
811 1926 : DO i = 1, n
812 : response_work = response_work + kpoint_weight* &
813 : (occupation_current(i)*REAL(hamiltonian_step_current(i, i), KIND=dp) - &
814 1926 : occupation_start(i)*REAL(hamiltonian_step_start(i, i), KIND=dp))
815 : END DO
816 :
817 1040 : ALLOCATE (cross_mode(n, n))
818 780 : 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 520 : overlap_start_current, "C", cross_mode)
822 3852 : 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 3852 : 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 27844 : AIMAG(hamiltonian_step_current)*AIMAG(density_modes(:, :, mode)))
832 : END DO
833 260 : DEALLOCATE (cross_mode)
834 :
835 : density_scale = kpoint_weight* &
836 3592 : MAX(MAXVAL(ABS(occupation_start)), MAXVAL(ABS(occupation_current)))
837 260 : IF (density_scale <= TINY(1.0_dp)) RETURN
838 260 : 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 2928 : ALL(ABS(response_overlap) <= HUGE(1.0_dp))
845 244 : 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 354 : SUBROUTINE qs_ot_finite_rotation_response( &
962 354 : chc, rotation_generator, occupation, kpoint_weight, rotation_gradient, &
963 354 : 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 354 : COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: base_hamiltonian, generator_minus, &
972 354 : generator_plus, rotation
973 : INTEGER :: i, j, n, nrotation, r, s
974 : REAL(KIND=dp) :: step
975 354 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: occupation_scale, rayleigh, &
976 354 : rayleigh_minus, rayleigh_plus
977 354 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: gradient_imag, gradient_minus_imag, &
978 354 : gradient_minus_real, gradient_plus_imag, gradient_plus_real, gradient_real
979 :
980 354 : n = SIZE(chc, 1)
981 354 : nrotation = n*(n - 1)
982 354 : CPASSERT(n > 0)
983 1062 : CPASSERT(ALL(SHAPE(chc) == [n, n]))
984 1062 : CPASSERT(ALL(SHAPE(rotation_generator) == [n, n]))
985 354 : CPASSERT(SIZE(occupation) == n)
986 354 : CPASSERT(SIZE(rotation_gradient) == nrotation)
987 1062 : CPASSERT(ALL(SHAPE(rotation_hessian) == [nrotation, nrotation]))
988 1062 : CPASSERT(ALL(SHAPE(rayleigh_response) == [n, nrotation]))
989 354 : CPASSERT(kpoint_weight > 0.0_dp)
990 :
991 354 : step = 1.0E-4_dp
992 354 : IF (PRESENT(difference_step)) step = difference_step
993 354 : 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 9912 : gradient_plus_real(n, n), gradient_real(n, n))
1000 :
1001 354 : CALL qs_ot_dense_rotation_state(rotation_generator, rotation)
1002 318920 : base_hamiltonian(:, :) = MATMUL(rotation, MATMUL(chc, CONJG(TRANSPOSE(rotation))))
1003 2318 : occupation_scale(:) = 2.0_dp*kpoint_weight*occupation(:)
1004 : CALL qs_ot_dense_rotation_gradient(rotation_generator, base_hamiltonian, occupation_scale, &
1005 354 : gradient_real, gradient_imag, rayleigh)
1006 :
1007 354 : r = 0
1008 1964 : DO i = 1, n - 1
1009 7406 : DO j = i + 1, n
1010 5442 : r = r + 1
1011 5442 : rotation_gradient(r) = gradient_real(i, j)
1012 5442 : r = r + 1
1013 7052 : rotation_gradient(r) = gradient_imag(i, j)
1014 : END DO
1015 : END DO
1016 354 : CPASSERT(r == nrotation)
1017 :
1018 : s = 0
1019 1964 : DO i = 1, n - 1
1020 7406 : DO j = i + 1, n
1021 393030 : generator_plus(:, :) = rotation_generator(:, :)
1022 393030 : generator_minus(:, :) = rotation_generator(:, :)
1023 5442 : generator_plus(i, j) = generator_plus(i, j) + step
1024 5442 : generator_plus(j, i) = generator_plus(j, i) - step
1025 5442 : generator_minus(i, j) = generator_minus(i, j) - step
1026 5442 : generator_minus(j, i) = generator_minus(j, i) + step
1027 5442 : s = s + 1
1028 : CALL qs_ot_dense_rotation_gradient(generator_plus, base_hamiltonian, occupation_scale, &
1029 5442 : gradient_plus_real, gradient_plus_imag, rayleigh_plus)
1030 : CALL qs_ot_dense_rotation_gradient(generator_minus, base_hamiltonian, occupation_scale, &
1031 5442 : 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 5442 : rayleigh_response(:, s))
1036 :
1037 393030 : generator_plus(:, :) = rotation_generator(:, :)
1038 393030 : generator_minus(:, :) = rotation_generator(:, :)
1039 5442 : generator_plus(i, j) = generator_plus(i, j) + CMPLX(0.0_dp, step, KIND=dp)
1040 5442 : generator_plus(j, i) = generator_plus(j, i) + CMPLX(0.0_dp, step, KIND=dp)
1041 5442 : generator_minus(i, j) = generator_minus(i, j) - CMPLX(0.0_dp, step, KIND=dp)
1042 5442 : generator_minus(j, i) = generator_minus(j, i) - CMPLX(0.0_dp, step, KIND=dp)
1043 5442 : s = s + 1
1044 : CALL qs_ot_dense_rotation_gradient(generator_plus, base_hamiltonian, occupation_scale, &
1045 5442 : gradient_plus_real, gradient_plus_imag, rayleigh_plus)
1046 : CALL qs_ot_dense_rotation_gradient(generator_minus, base_hamiltonian, occupation_scale, &
1047 5442 : 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 7052 : rayleigh_response(:, s))
1052 : END DO
1053 : END DO
1054 354 : CPASSERT(s == nrotation)
1055 1240380 : 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 354 : gradient_minus_real, gradient_plus_imag, gradient_plus_real, gradient_real)
1060 :
1061 354 : 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 10884 : SUBROUTINE qs_ot_pack_rotation_response( &
1076 10884 : gradient_plus_real, gradient_plus_imag, gradient_minus_real, gradient_minus_imag, &
1077 10884 : 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 10884 : n = SIZE(gradient_plus_real, 1)
1088 32652 : CPASSERT(ALL(SHAPE(gradient_plus_real) == [n, n]))
1089 32652 : CPASSERT(ALL(SHAPE(gradient_plus_imag) == [n, n]))
1090 32652 : CPASSERT(ALL(SHAPE(gradient_minus_real) == [n, n]))
1091 32652 : CPASSERT(ALL(SHAPE(gradient_minus_imag) == [n, n]))
1092 10884 : CPASSERT(SIZE(hessian_column) == n*(n - 1))
1093 10884 : CPASSERT(SIZE(rayleigh_column) == n)
1094 :
1095 10884 : r = 0
1096 83112 : DO i = 1, n - 1
1097 387588 : DO j = i + 1, n
1098 304476 : r = r + 1
1099 : hessian_column(r) = &
1100 304476 : (gradient_plus_real(i, j) - gradient_minus_real(i, j))/(2.0_dp*step)
1101 304476 : r = r + 1
1102 : hessian_column(r) = &
1103 376704 : (gradient_plus_imag(i, j) - gradient_minus_imag(i, j))/(2.0_dp*step)
1104 : END DO
1105 : END DO
1106 93996 : rayleigh_column(:) = (rayleigh_plus(:) - rayleigh_minus(:))/(2.0_dp*step)
1107 :
1108 10884 : 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 22122 : SUBROUTINE qs_ot_dense_rotation_gradient( &
1120 22122 : 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 22122 : COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: covector, eigenvectors, frechet, inner, &
1128 22122 : outer, rotation, work
1129 : INTEGER :: i, j, n
1130 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eigenvalues
1131 :
1132 22122 : n = SIZE(rotation_generator, 1)
1133 : ALLOCATE (covector(n, n), eigenvectors(n, n), frechet(n, n), inner(n, n), &
1134 398196 : outer(n, n), rotation(n, n), work(n, n), eigenvalues(n))
1135 : CALL qs_ot_dense_rotation_state(rotation_generator, rotation, &
1136 22122 : eigenvectors=eigenvectors, eigenvalues=eigenvalues)
1137 :
1138 13892046 : work(:, :) = MATMUL(base_hamiltonian, rotation)
1139 1587286 : covector(:, :) = work(:, :)
1140 190310 : DO j = 1, n
1141 1587286 : covector(:, j) = occupation_scale(j)*covector(:, j)
1142 : END DO
1143 :
1144 22122 : work(:, :) = MATMUL(CONJG(TRANSPOSE(rotation)), &
1145 40155218 : MATMUL(base_hamiltonian, rotation))
1146 190310 : DO i = 1, n
1147 190310 : rayleigh(i) = REAL(work(i, i), KIND=dp)
1148 : END DO
1149 :
1150 22122 : inner(:, :) = MATMUL(CONJG(TRANSPOSE(eigenvectors)), &
1151 40155218 : MATMUL(covector, eigenvectors))
1152 190310 : DO j = 1, n
1153 1587286 : DO i = 1, n
1154 1396976 : kernel = CONJG(qs_ot_complex_exp_frechet_kernel(eigenvalues(i), eigenvalues(j)))
1155 1565164 : outer(i, j) = inner(i, j)*kernel
1156 : END DO
1157 : END DO
1158 22122 : frechet(:, :) = MATMUL(eigenvectors, &
1159 40155218 : MATMUL(outer, CONJG(TRANSPOSE(eigenvectors))))
1160 1587286 : gradient_real(:, :) = REAL(frechet, KIND=dp) - TRANSPOSE(REAL(frechet, KIND=dp))
1161 1587286 : gradient_imag(:, :) = AIMAG(frechet) + TRANSPOSE(AIMAG(frechet))
1162 :
1163 22122 : DEALLOCATE (covector, eigenvectors, frechet, inner, outer, rotation, work, eigenvalues)
1164 :
1165 22122 : 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 23814 : 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 23814 : n = SIZE(rotation_generator, 1)
1186 190512 : 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 1674398 : vectors, values)
1189 1674398 : weighted_vectors(:, :) = vectors(:, :)
1190 202612 : DO j = 1, n
1191 : weighted_vectors(:, j) = EXP(CMPLX(0.0_dp, -values(j), KIND=dp))* &
1192 1674398 : weighted_vectors(:, j)
1193 : END DO
1194 29097282 : rotation(:, :) = MATMUL(weighted_vectors, CONJG(TRANSPOSE(vectors)))
1195 1588978 : IF (PRESENT(eigenvectors)) eigenvectors(:, :) = vectors(:, :)
1196 192002 : IF (PRESENT(eigenvalues)) eigenvalues(:) = values(:)
1197 23814 : DEALLOCATE (vectors, weighted_vectors, values)
1198 :
1199 23814 : 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 1514546 : 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 1514546 : l1 = (0.0_dp, -1.0_dp)*e1
1215 1514546 : l2 = (0.0_dp, -1.0_dp)*e2
1216 1514546 : IF (ABS(l1 - l2) > 0.5_dp) THEN
1217 8660 : kernel = (EXP(l1) - EXP(l2))/(l1 - l2)
1218 : ELSE
1219 : x = 1.0_dp
1220 : kernel = 0.0_dp
1221 25600062 : DO i = 1, 16
1222 24094176 : kernel = kernel + x
1223 25600062 : x = x*(l1 - l2)/REAL(i + 1, KIND=dp)
1224 : END DO
1225 1505886 : kernel = kernel*EXP(l2)
1226 : END IF
1227 :
1228 1514546 : 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 2039 : 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 2039 : INTEGER, ALLOCATABLE, DIMENSION(:) :: cols, rows
1249 2039 : INTEGER, DIMENSION(:), POINTER :: col_blk_offset, col_blk_size, &
1250 2039 : row_blk_offset, row_blk_size
1251 : LOGICAL :: found_im, found_out_im, found_out_re, &
1252 : found_re, use_adjoint
1253 2039 : REAL(dp), DIMENSION(:, :), POINTER :: block_in_im, block_in_re, block_out_im, &
1254 2039 : 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 2039 : use_adjoint = .FALSE.
1260 2037 : 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 2039 : max_blocks = 0
1265 2039 : CALL dbcsr_iterator_start(iter, inner_deriv_re)
1266 3211 : DO WHILE (dbcsr_iterator_blocks_left(iter))
1267 1172 : CALL dbcsr_iterator_next_block(iter, row, col)
1268 1172 : max_blocks = max_blocks + 1
1269 : END DO
1270 2039 : CALL dbcsr_iterator_stop(iter)
1271 2039 : CALL dbcsr_iterator_start(iter, inner_deriv_im)
1272 3211 : DO WHILE (dbcsr_iterator_blocks_left(iter))
1273 1172 : CALL dbcsr_iterator_next_block(iter, row, col)
1274 1172 : max_blocks = max_blocks + 1
1275 : END DO
1276 2039 : CALL dbcsr_iterator_stop(iter)
1277 8156 : ALLOCATE (rows(MAX(max_blocks, 1)), cols(MAX(max_blocks, 1)))
1278 2039 : nblocks = 0
1279 :
1280 2039 : CALL dbcsr_iterator_start(iter, inner_deriv_re)
1281 3211 : DO WHILE (dbcsr_iterator_blocks_left(iter))
1282 1172 : CALL dbcsr_iterator_next_block(iter, row, col)
1283 1172 : CALL append_union_block(row, col)
1284 : END DO
1285 2039 : CALL dbcsr_iterator_stop(iter)
1286 2039 : CALL dbcsr_iterator_start(iter, inner_deriv_im)
1287 3211 : DO WHILE (dbcsr_iterator_blocks_left(iter))
1288 1172 : CALL dbcsr_iterator_next_block(iter, row, col)
1289 1172 : CALL append_union_block(row, col)
1290 : END DO
1291 2039 : CALL dbcsr_iterator_stop(iter)
1292 :
1293 : CALL dbcsr_get_info(inner_deriv_re, distribution=dist, row_blk_size=row_blk_size, &
1294 2039 : col_blk_size=col_blk_size)
1295 : CALL dbcsr_create(outer_deriv_re, "outer_deriv_re", dist, dbcsr_type_no_symmetry, &
1296 2039 : row_blk_size, col_blk_size)
1297 : CALL dbcsr_create(outer_deriv_im, "outer_deriv_im", dist, dbcsr_type_no_symmetry, &
1298 2039 : row_blk_size, col_blk_size)
1299 2039 : IF (nblocks > 0) THEN
1300 1172 : CALL dbcsr_reserve_blocks(outer_deriv_re, rows=rows(1:nblocks), cols=cols(1:nblocks))
1301 1172 : CALL dbcsr_reserve_blocks(outer_deriv_im, rows=rows(1:nblocks), cols=cols(1:nblocks))
1302 : END IF
1303 2039 : CALL dbcsr_finalize(outer_deriv_re)
1304 2039 : CALL dbcsr_finalize(outer_deriv_im)
1305 2039 : CALL dbcsr_set(outer_deriv_re, 0.0_dp)
1306 2039 : 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 2039 : col_blk_offset=col_blk_offset)
1310 2039 : CALL dbcsr_iterator_start(iter, outer_deriv_re)
1311 3213 : DO WHILE (dbcsr_iterator_blocks_left(iter))
1312 1174 : CALL dbcsr_iterator_next_block(iter, row, col)
1313 1174 : CALL dbcsr_get_block_p(inner_deriv_re, row, col, block_in_re, found_re)
1314 1174 : CALL dbcsr_get_block_p(inner_deriv_im, row, col, block_in_im, found_im)
1315 1174 : CALL dbcsr_get_block_p(outer_deriv_re, row, col, block_out_re, found_out_re)
1316 1174 : CALL dbcsr_get_block_p(outer_deriv_im, row, col, block_out_im, found_out_im)
1317 1174 : CPASSERT(found_out_re .AND. found_out_im)
1318 :
1319 12837 : DO i = 1, SIZE(block_out_re, 1)
1320 128346 : DO j = 1, SIZE(block_out_re, 2)
1321 117548 : e1 = evals(row_blk_offset(row) + i - 1)
1322 117548 : e2 = evals(col_blk_offset(col) + j - 1)
1323 117548 : re_part = 0.0_dp
1324 117548 : im_part = 0.0_dp
1325 117548 : IF (found_re) re_part = block_in_re(i, j)
1326 117548 : IF (found_im) im_part = block_in_im(i, j)
1327 117548 : cval_in = CMPLX(re_part, im_part, dp)
1328 117548 : kernel = qs_ot_complex_exp_frechet_kernel(e1, e2)
1329 117548 : IF (use_adjoint) kernel = CONJG(kernel)
1330 117548 : cval_in = cval_in*kernel
1331 117548 : block_out_re(i, j) = REAL(cval_in, KIND=dp)
1332 127172 : block_out_im(i, j) = AIMAG(cval_in)
1333 : END DO
1334 : END DO
1335 : END DO
1336 2039 : CALL dbcsr_iterator_stop(iter)
1337 16310 : 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 2344 : SUBROUTINE append_union_block(row_new, col_new)
1347 : INTEGER, INTENT(IN) :: row_new, col_new
1348 :
1349 : INTEGER :: iblock
1350 :
1351 2346 : DO iblock = 1, nblocks
1352 2346 : IF (rows(iblock) == row_new .AND. cols(iblock) == col_new) RETURN
1353 : END DO
1354 1174 : nblocks = nblocks + 1
1355 1174 : rows(nblocks) = row_new
1356 1174 : 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 9001 : 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 9001 : qs_ot_env%preconditioner => preconditioner
1377 9001 : qs_ot_env%os_valid = .FALSE.
1378 9001 : IF (.NOT. ASSOCIATED(qs_ot_env%matrix_psc0)) THEN
1379 9001 : CALL dbcsr_init_p(qs_ot_env%matrix_psc0)
1380 9001 : CALL dbcsr_copy(qs_ot_env%matrix_psc0, qs_ot_env%matrix_sc0, 'matrix_psc0')
1381 : END IF
1382 9001 : IF (qs_ot_env%has_complex_kpoint_state .AND. &
1383 : .NOT. ASSOCIATED(qs_ot_env%matrix_psc0_im)) THEN
1384 423 : CALL dbcsr_init_p(qs_ot_env%matrix_psc0_im)
1385 423 : CALL dbcsr_copy(qs_ot_env%matrix_psc0_im, qs_ot_env%matrix_sc0_im, 'matrix_psc0_im')
1386 : END IF
1387 :
1388 9001 : IF (.NOT. qs_ot_env%use_dx) THEN
1389 5220 : qs_ot_env%use_dx = .TRUE.
1390 5220 : CALL dbcsr_init_p(qs_ot_env%matrix_dx)
1391 5220 : CALL dbcsr_copy(qs_ot_env%matrix_dx, qs_ot_env%matrix_gx, 'matrix_dx')
1392 5220 : IF (qs_ot_env%has_complex_kpoint_state) THEN
1393 123 : CALL dbcsr_init_p(qs_ot_env%matrix_dx_im)
1394 123 : CALL dbcsr_copy(qs_ot_env%matrix_dx_im, qs_ot_env%matrix_gx_im, 'matrix_dx_im')
1395 : END IF
1396 5220 : IF (qs_ot_env%settings%do_rotation) THEN
1397 86 : CALL dbcsr_init_p(qs_ot_env%rot_mat_dx)
1398 86 : CALL dbcsr_copy(qs_ot_env%rot_mat_dx, qs_ot_env%rot_mat_gx, 'rot_mat_dx')
1399 86 : IF (qs_ot_env%has_complex_kpoint_state) THEN
1400 56 : CALL dbcsr_init_p(qs_ot_env%rot_mat_dx_im)
1401 56 : 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 5220 : IF (qs_ot_env%settings%do_ener) THEN
1405 48 : ncoef = SIZE(qs_ot_env%ener_gx)
1406 144 : ALLOCATE (qs_ot_env%ener_dx(ncoef))
1407 384 : qs_ot_env%ener_dx = 0.0_dp
1408 : END IF
1409 : END IF
1410 :
1411 9001 : 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 40184 : 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 27380 : SELECT CASE (op_a)
1433 : CASE ('N')
1434 27380 : db_op_a = 'N'
1435 27380 : sign_a = 1.0_dp
1436 : CASE ('C')
1437 12804 : db_op_a = 'T'
1438 12804 : sign_a = -1.0_dp
1439 : CASE DEFAULT
1440 40184 : CPABORT("Complex matrix product expects N or C for op_a")
1441 : END SELECT
1442 34092 : SELECT CASE (op_b)
1443 : CASE ('N')
1444 34092 : db_op_b = 'N'
1445 34092 : sign_b = 1.0_dp
1446 : CASE ('C')
1447 6092 : db_op_b = 'T'
1448 6092 : sign_b = -1.0_dp
1449 : CASE DEFAULT
1450 40184 : CPABORT("Complex matrix product expects N or C for op_b")
1451 : END SELECT
1452 :
1453 40184 : CALL dbcsr_multiply(db_op_a, db_op_b, 1.0_dp, a_re, b_re, 0.0_dp, c_re)
1454 40184 : CALL dbcsr_multiply(db_op_a, db_op_b, 1.0_dp, a_im, b_im, 0.0_dp, tmp)
1455 40184 : CALL dbcsr_add(c_re, tmp, alpha_scalar=1.0_dp, beta_scalar=-sign_a*sign_b)
1456 :
1457 40184 : CALL dbcsr_multiply(db_op_a, db_op_b, sign_b, a_re, b_im, 0.0_dp, c_im)
1458 40184 : CALL dbcsr_multiply(db_op_a, db_op_b, sign_a, a_im, b_re, 0.0_dp, tmp)
1459 40184 : CALL dbcsr_add(c_im, tmp, alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
1460 :
1461 40184 : 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 420 : 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 84 : 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 84 : CALL timeset(routineN, handle)
1490 : !
1491 : !
1492 84 : CALL dbcsr_get_info(C_NEW, nfullrows_total=n, nfullcols_total=k)
1493 : !
1494 : ! C = C*expm(-G)
1495 84 : GU => qs_ot_env%buf1_k_k_nosym ! a buffer
1496 84 : U => qs_ot_env%buf2_k_k_nosym ! a buffer
1497 84 : Gp1 => qs_ot_env%buf3_k_k_nosym ! a buffer
1498 84 : Gp2 => qs_ot_env%buf4_k_k_nosym ! a buffer
1499 84 : 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 84 : f2 = 0.0_dp
1505 84 : CALL dbcsr_copy(C, C_NEW)
1506 84 : CALL dbcsr_iterator_start(iter, C)
1507 182 : DO WHILE (dbcsr_iterator_blocks_left(iter))
1508 98 : CALL dbcsr_iterator_next_block(iter, row, col, block, row_size=row_size, col_size=col_size)
1509 686 : DO p = 1, col_size ! p
1510 6258 : DO i = 1, row_size ! i
1511 5656 : tmp = SQRT(block(i, p)**2 + f2_eps)
1512 5656 : f2 = f2 + tmp
1513 6160 : block(i, p) = block(i, p)/tmp
1514 : END DO
1515 : END DO
1516 : END DO
1517 84 : CALL dbcsr_iterator_stop(iter)
1518 84 : CALL dbcsr_get_info(C, group=group)
1519 84 : CALL group%sum(f2)
1520 : !
1521 : !
1522 84 : CALL dbcsr_multiply('T', 'N', 1.0_dp, C, C_NEW, 0.0_dp, GU)
1523 : !
1524 : ! antisymetrize
1525 84 : CALL dbcsr_get_info(GU, distribution=dist)
1526 : CALL dbcsr_transposed(U, GU, shallow_data_copy=.FALSE., &
1527 : use_distribution=dist, &
1528 84 : transpose_distribution=.FALSE.)
1529 84 : CALL dbcsr_add(GU, U, alpha_scalar=-0.5_dp, beta_scalar=0.5_dp)
1530 : !-------------------------------------------------------------------
1531 : !
1532 84 : norm_fro = dbcsr_frobenius_norm(GU)
1533 84 : 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 84 : CALL dbcsr_scale(GU, -alpha)
1542 : !
1543 : ! compute unitary transform
1544 : ! zeroth and first order
1545 84 : expfactor = 1.0_dp
1546 84 : CALL dbcsr_copy(U, GU)
1547 84 : CALL dbcsr_scale(U, expfactor)
1548 84 : CALL dbcsr_add_on_diag(U, 1.0_dp)
1549 : ! other orders
1550 84 : CALL dbcsr_copy(Gp1, GU)
1551 520 : DO i = 2, taylor_order
1552 : ! new power of G
1553 520 : CALL dbcsr_multiply('N', 'N', 1.0_dp, GU, Gp1, 0.0_dp, Gp2)
1554 520 : CALL dbcsr_copy(Gp1, Gp2)
1555 : ! add to the taylor expansion so far
1556 520 : expfactor = expfactor/REAL(i, KIND=dp)
1557 520 : CALL dbcsr_add(U, Gp1, alpha_scalar=1.0_dp, beta_scalar=expfactor)
1558 520 : norm_fro = dbcsr_frobenius_norm(Gp1)
1559 : !write(*,*) 'Taylor expansion i=',i,' norm(X^i)/i!=',norm_fro*expfactor
1560 520 : IF (norm_fro*expfactor < 1.0E-10_dp) EXIT
1561 : END DO
1562 : !
1563 : ! rotate MOs
1564 84 : CALL dbcsr_multiply('N', 'N', 1.0_dp, C_NEW, U, 0.0_dp, C)
1565 84 : CALL dbcsr_copy(C_NEW, C)
1566 : !
1567 : ! rotate SC
1568 84 : CALL dbcsr_multiply('N', 'N', 1.0_dp, SC, U, 0.0_dp, C)
1569 84 : CALL dbcsr_copy(SC, C)
1570 : !
1571 : ! rotate D_i
1572 84 : CALL dbcsr_multiply('N', 'N', 1.0_dp, D, U, 0.0_dp, C)
1573 84 : CALL dbcsr_copy(D, C)
1574 : !
1575 : ! rotate G_i-1
1576 84 : IF (ASSOCIATED(G_OLD)) THEN
1577 84 : CALL dbcsr_multiply('N', 'N', 1.0_dp, G_OLD, U, 0.0_dp, C)
1578 84 : CALL dbcsr_copy(G_OLD, C)
1579 : END IF
1580 : !
1581 84 : CALL timestop(handle)
1582 84 : 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 1660 : 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 830 : CALL timeset(routineN, handle)
1605 : !
1606 830 : CALL dbcsr_get_info(C_NEW, nfullrows_total=n, nfullcols_total=k)
1607 : !
1608 : ! P = U'*U
1609 830 : 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 830 : transa="N", para_env=qs_ot_env%para_env, blacs_env=qs_ot_env%blacs_env)
1614 : !
1615 : ! Update SC if needed
1616 830 : IF (update) THEN
1617 : CALL cp_dbcsr_cholesky_restore(SC, k, P, C_TMP, op="SOLVE", pos="RIGHT", &
1618 454 : transa="N", para_env=qs_ot_env%para_env, blacs_env=qs_ot_env%blacs_env)
1619 454 : CALL dbcsr_copy(SC, C_TMP)
1620 : END IF
1621 : !
1622 830 : CALL timestop(handle)
1623 830 : 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 308 : 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 308 : CALL timeset(routineN, handle)
1648 : !
1649 308 : CALL dbcsr_get_info(C_NEW, nfullrows_total=n, nfullcols_total=k)
1650 : !
1651 308 : V => qs_ot_env%buf1_k_k_nosym ! a buffer
1652 308 : W => qs_ot_env%buf2_k_k_nosym ! a buffer
1653 1232 : ALLOCATE (eig(k), fun(k))
1654 : !
1655 308 : 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 1796 : DO i = 1, k
1659 1488 : IF (eig(i) <= 0.0_dp) THEN
1660 0 : CPABORT("P not positive definite")
1661 : END IF
1662 1796 : IF (eig(i) < 1.0E-8_dp) THEN
1663 0 : fun(i) = 0.0_dp
1664 : ELSE
1665 1488 : fun(i) = 1.0_dp/SQRT(eig(i))
1666 : END IF
1667 : END DO
1668 308 : CALL dbcsr_copy(W, V)
1669 308 : CALL dbcsr_scale_by_vector(V, alpha=fun, side='right')
1670 308 : CALL dbcsr_multiply('N', 'T', 1.0_dp, W, V, 0.0_dp, P)
1671 : !
1672 : ! Update C
1673 308 : CALL dbcsr_multiply('N', 'N', 1.0_dp, C_OLD, P, 0.0_dp, C_NEW)
1674 : !
1675 : ! Update SC if needed
1676 308 : IF (update) THEN
1677 216 : CALL dbcsr_multiply('N', 'N', 1.0_dp, SC, P, 0.0_dp, C_TMP)
1678 216 : CALL dbcsr_copy(SC, C_TMP)
1679 : END IF
1680 : !
1681 308 : DEALLOCATE (eig, fun)
1682 : !
1683 308 : CALL timestop(handle)
1684 308 : 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 8180 : 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 4090 : CALL timeset(routineN, handle)
1714 : !
1715 4090 : CALL dbcsr_get_info(C_NEW, nfullrows_total=n, nfullcols_total=k)
1716 : !
1717 4090 : BUF_NOSYM => qs_ot_env%buf1_k_k_nosym ! a buffer
1718 4090 : BUF1 => qs_ot_env%buf1_k_k_sym ! a buffer
1719 4090 : BUF2 => qs_ot_env%buf2_k_k_sym ! a buffer
1720 4090 : FY => qs_ot_env%buf3_k_k_sym ! a buffer
1721 4090 : FT => qs_ot_env%buf4_k_k_sym ! a buffer
1722 : !
1723 : ! initialize the norm (already computed in qs_ot_get_orbitals_ref)
1724 4090 : norm = norm_in
1725 : !
1726 : ! can we do a quick exit?
1727 4090 : quick_exit = .FALSE.
1728 4090 : IF (norm < qs_ot_env%settings%eps_irac_quick_exit) quick_exit = .TRUE.
1729 : !
1730 : ! lets refine
1731 4090 : rescale = 1.0_dp
1732 4524 : DO irefine = 1, qs_ot_env%settings%max_irac
1733 : !
1734 : ! rescaling
1735 4524 : 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 4524 : qs_ot_env%settings%eps_irac_filter_matrix)
1743 : !
1744 : ! collect the transformation
1745 4524 : IF (irefine == 1) THEN
1746 4090 : 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 4524 : IF (quick_exit) THEN
1759 : EXIT
1760 : END IF
1761 : !
1762 : ! P = FY^T * P * FY
1763 1880 : CALL dbcsr_multiply('N', 'N', 1.0_dp, P, FY, 0.0_dp, BUF_NOSYM)
1764 1880 : 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 1880 : CALL dbcsr_multiply('N', 'N', 1.0_dp, FY, BUF_NOSYM, 0.0_dp, P)
1770 1880 : 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 1880 : CALL dbcsr_add_on_diag(P, -1.0_dp)
1778 1880 : norm_fro = dbcsr_frobenius_norm(P)
1779 1880 : norm_gct = dbcsr_gershgorin_norm(P)
1780 1880 : CALL dbcsr_add_on_diag(P, 1.0_dp)
1781 1880 : norm = MIN(norm_gct, norm_fro)
1782 : !
1783 : ! printing
1784 : !
1785 : ! blows up
1786 1880 : 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 1880 : IF (norm < qs_ot_env%settings%eps_irac_quick_exit) quick_exit = .TRUE.
1795 : !
1796 : ! are we done?
1797 4524 : IF (norm < qs_ot_env%settings%eps_irac) EXIT
1798 : !
1799 : END DO
1800 : !
1801 : ! C_NEW = C_NEW * FT * rescale
1802 4090 : CALL dbcsr_multiply('N', 'N', rescale, C_OLD, FT, 0.0_dp, C_NEW)
1803 4090 : 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 4090 : IF (update) THEN
1811 1670 : CALL dbcsr_multiply('N', 'N', rescale, SC, FT, 0.0_dp, C_TMP)
1812 1670 : 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 1670 : CALL dbcsr_copy(SC, C_TMP)
1818 : END IF
1819 : !
1820 4090 : CALL timestop(handle)
1821 4090 : END SUBROUTINE qs_ot_ref_poly
1822 :
1823 : ! **************************************************************************************************
1824 : !> \brief ...
1825 : !> \param qs_ot_env1 ...
1826 : !> \return ...
1827 : ! **************************************************************************************************
1828 8783 : FUNCTION qs_ot_ref_update(qs_ot_env1) RESULT(update)
1829 : !
1830 : TYPE(qs_ot_type) :: qs_ot_env1
1831 : LOGICAL :: update
1832 :
1833 8783 : update = .FALSE.
1834 6716 : SELECT CASE (qs_ot_env1%settings%ot_method)
1835 : CASE ("CG", "SD")
1836 6716 : SELECT CASE (qs_ot_env1%settings%line_search_method)
1837 : CASE ("2PNT")
1838 6716 : IF (qs_ot_env1%line_search_count == 2) update = .TRUE.
1839 : CASE DEFAULT
1840 6716 : 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 8783 : CPABORT("NYI")
1849 : END SELECT
1850 8783 : 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 5228 : 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 5228 : ortho_irac = qs_ot_env1%settings%ortho_irac
1865 5228 : IF (norm_in < qs_ot_env1%settings%eps_irac_switch) ortho_irac = "POLY"
1866 5228 : 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 10456 : 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 5228 : CALL timeset(routineN, handle)
1895 :
1896 5228 : CALL dbcsr_get_info(matrix_c, nfullrows_total=n, nfullcols_total=k)
1897 : !
1898 5228 : C_NEW => matrix_c
1899 5228 : C_OLD => matrix_x ! need to be carefully updated for the gradient !
1900 5228 : SC => matrix_sx ! need to be carefully updated for the gradient !
1901 5228 : G_OLD => matrix_gx_old ! need to be carefully updated for localization !
1902 5228 : D => matrix_dx ! need to be carefully updated for localization !
1903 5228 : S => matrix_s
1904 :
1905 5228 : P => qs_ot_env%p_k_k_sym ! a buffer
1906 5228 : C_TMP => qs_ot_env%buf1_n_k ! a buffer
1907 : !
1908 : ! do we need to update C_OLD and SC?
1909 5228 : 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 5228 : on_the_fly_loc = qs_ot_env1%settings%on_the_fly_loc
1916 : !
1917 : ! compute SC = S*C
1918 5228 : IF (ASSOCIATED(S)) THEN
1919 5228 : CALL dbcsr_multiply('N', 'N', 1.0_dp, S, C_OLD, 0.0_dp, SC)
1920 5228 : 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 5228 : CALL dbcsr_multiply('T', 'N', 1.0_dp, C_OLD, SC, 0.0_dp, P)
1931 5228 : 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 5228 : CALL dbcsr_add_on_diag(P, -1.0_dp)
1939 5228 : norm_fro = dbcsr_frobenius_norm(P)
1940 5228 : norm_gct = dbcsr_gershgorin_norm(P)
1941 5228 : CALL dbcsr_add_on_diag(P, 1.0_dp)
1942 5228 : norm = MIN(norm_gct, norm_fro)
1943 5228 : CALL qs_ot_ref_decide(qs_ot_env1, norm, ortho_irac)
1944 : !
1945 : ! select the orthogonality method
1946 830 : SELECT CASE (ortho_irac)
1947 : CASE ("CHOL")
1948 830 : CALL qs_ot_ref_chol(qs_ot_env, C_OLD, C_TMP, C_NEW, P, SC, update)
1949 : CASE ("LWDN")
1950 308 : CALL qs_ot_ref_lwdn(qs_ot_env, C_OLD, C_TMP, C_NEW, P, SC, update)
1951 : CASE ("POLY")
1952 4090 : CALL qs_ot_ref_poly(qs_ot_env, C_OLD, C_TMP, C_NEW, P, SC, norm, update)
1953 : CASE DEFAULT
1954 5228 : CPABORT("Wrong argument")
1955 : END SELECT
1956 : !
1957 : ! We update the C_i+1 and localization
1958 5228 : IF (update) THEN
1959 2340 : IF (on_the_fly_loc) THEN
1960 84 : CALL qs_ot_on_the_fly_localize(qs_ot_env, C_NEW, SC, G_OLD, D)
1961 : END IF
1962 2340 : CALL dbcsr_copy(C_OLD, C_NEW)
1963 : END IF
1964 :
1965 5228 : 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 5228 : CALL timestop(handle)
1973 5228 : 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 3555 : 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 3555 : 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 3555 : CALL timeset(routineN, handle)
2003 :
2004 3555 : CPASSERT(qs_ot_env%has_complex_kpoint_state)
2005 3555 : CPASSERT(ASSOCIATED(matrix_s))
2006 3555 : CPASSERT(ASSOCIATED(matrix_s_im))
2007 :
2008 3555 : c_re => qs_ot_env%matrix_x
2009 3555 : c_im => qs_ot_env%matrix_x_im
2010 3555 : f_re => qs_ot_env%matrix_ref_inv_sqrt
2011 3555 : f_im => qs_ot_env%matrix_ref_inv_sqrt_im
2012 3555 : sc_re => qs_ot_env%matrix_sx
2013 3555 : sc_im => qs_ot_env%matrix_sx_im
2014 3555 : p_re => qs_ot_env%buf1_k_k_sym
2015 3555 : p_im => qs_ot_env%buf2_k_k_sym
2016 3555 : v_re => qs_ot_env%buf3_k_k_sym
2017 3555 : v_im => qs_ot_env%buf4_k_k_sym
2018 3555 : w_re => qs_ot_env%buf1_k_k_nosym
2019 3555 : w_im => qs_ot_env%buf2_k_k_nosym
2020 3555 : tmp_kk => qs_ot_env%buf3_k_k_nosym
2021 3555 : tmp_nk => qs_ot_env%buf1_n_k
2022 :
2023 3555 : CALL dbcsr_get_info(c_re, nfullrows_total=n, nfullcols_total=k)
2024 3555 : IF (PRESENT(qs_ot_env1)) THEN
2025 3379 : update = qs_ot_ref_update(qs_ot_env1)
2026 : ELSE
2027 176 : 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 3555 : CALL dbcsr_multiply('N', 'N', 1.0_dp, matrix_s, c_re, 0.0_dp, sc_re)
2032 3555 : CALL dbcsr_multiply('N', 'N', 1.0_dp, matrix_s_im, c_im, 0.0_dp, tmp_nk)
2033 3555 : CALL dbcsr_add(sc_re, tmp_nk, alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
2034 :
2035 3555 : CALL dbcsr_multiply('N', 'N', 1.0_dp, matrix_s, c_im, 0.0_dp, sc_im)
2036 3555 : CALL dbcsr_multiply('N', 'N', 1.0_dp, matrix_s_im, c_re, 0.0_dp, tmp_nk)
2037 3555 : 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 3555 : CALL dbcsr_multiply('T', 'N', 1.0_dp, c_re, sc_re, 0.0_dp, p_re)
2041 3555 : CALL dbcsr_multiply('T', 'N', 1.0_dp, c_im, sc_im, 0.0_dp, tmp_kk)
2042 3555 : CALL dbcsr_add(p_re, tmp_kk, alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
2043 :
2044 3555 : CALL dbcsr_multiply('T', 'N', 1.0_dp, c_re, sc_im, 0.0_dp, p_im)
2045 3555 : CALL dbcsr_multiply('T', 'N', 1.0_dp, c_im, sc_re, 0.0_dp, tmp_kk)
2046 3555 : CALL dbcsr_add(p_im, tmp_kk, alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
2047 :
2048 14220 : 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 3555 : blacs_env=qs_ot_env%blacs_env)
2053 29733 : DO i = 1, k
2054 26178 : IF (eigenvalues(i) <= EPSILON(1.0_dp)) THEN
2055 0 : CPABORT("Complex REF overlap is not positive definite")
2056 : END IF
2057 29733 : 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 3555 : CALL dbcsr_copy(w_re, v_re)
2062 3555 : CALL dbcsr_copy(w_im, v_im)
2063 3555 : CALL dbcsr_scale_by_vector(w_re, alpha=inverse_sqrt, side='right')
2064 3555 : CALL dbcsr_scale_by_vector(w_im, alpha=inverse_sqrt, side='right')
2065 :
2066 3555 : CALL dbcsr_multiply('N', 'T', 1.0_dp, w_re, v_re, 0.0_dp, p_re)
2067 3555 : CALL dbcsr_multiply('N', 'T', 1.0_dp, w_im, v_im, 0.0_dp, tmp_kk)
2068 3555 : CALL dbcsr_add(p_re, tmp_kk, alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
2069 :
2070 3555 : CALL dbcsr_multiply('N', 'T', 1.0_dp, w_im, v_re, 0.0_dp, p_im)
2071 3555 : CALL dbcsr_multiply('N', 'T', 1.0_dp, w_re, v_im, 0.0_dp, tmp_kk)
2072 3555 : CALL dbcsr_add(p_im, tmp_kk, alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
2073 :
2074 3555 : CALL dbcsr_copy(f_re, p_re)
2075 3555 : CALL dbcsr_copy(f_im, p_im)
2076 :
2077 : ! Return the physical, orthonormal orbitals without changing a rejected
2078 : ! line-search coordinate.
2079 3555 : CALL dbcsr_multiply('N', 'N', 1.0_dp, c_re, p_re, 0.0_dp, matrix_c)
2080 3555 : CALL dbcsr_multiply('N', 'N', 1.0_dp, c_im, p_im, 0.0_dp, tmp_nk)
2081 3555 : CALL dbcsr_add(matrix_c, tmp_nk, alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
2082 :
2083 3555 : CALL dbcsr_multiply('N', 'N', 1.0_dp, c_re, p_im, 0.0_dp, matrix_c_im)
2084 3555 : CALL dbcsr_multiply('N', 'N', 1.0_dp, c_im, p_re, 0.0_dp, tmp_nk)
2085 3555 : CALL dbcsr_add(matrix_c_im, tmp_nk, alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
2086 :
2087 3555 : IF (update) THEN
2088 1485 : CALL dbcsr_copy(c_re, matrix_c)
2089 1485 : CALL dbcsr_copy(c_im, matrix_c_im)
2090 :
2091 1485 : CALL dbcsr_set(f_re, 0.0_dp)
2092 1485 : CALL dbcsr_add_on_diag(f_re, alpha=1.0_dp)
2093 1485 : CALL dbcsr_set(f_im, 0.0_dp)
2094 :
2095 1485 : CALL dbcsr_multiply('N', 'N', 1.0_dp, matrix_s, c_re, 0.0_dp, sc_re)
2096 1485 : CALL dbcsr_multiply('N', 'N', 1.0_dp, matrix_s_im, c_im, 0.0_dp, tmp_nk)
2097 1485 : CALL dbcsr_add(sc_re, tmp_nk, alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
2098 :
2099 1485 : CALL dbcsr_multiply('N', 'N', 1.0_dp, matrix_s, c_im, 0.0_dp, sc_im)
2100 1485 : CALL dbcsr_multiply('N', 'N', 1.0_dp, matrix_s_im, c_re, 0.0_dp, tmp_nk)
2101 1485 : CALL dbcsr_add(sc_im, tmp_nk, alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
2102 : END IF
2103 :
2104 3555 : IF (qs_ot_env%settings%do_rotation) THEN
2105 2150 : CALL qs_ot_generate_rotation_complex(qs_ot_env)
2106 :
2107 2150 : CALL dbcsr_copy(rotated_re, matrix_c, name="rotated_re")
2108 2150 : CALL dbcsr_copy(rotated_im, matrix_c_im, name="rotated_im")
2109 2150 : 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 2150 : 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 2150 : 0.0_dp, rotation_tmp)
2116 2150 : 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 2150 : 0.0_dp, rotated_im)
2120 : CALL dbcsr_multiply('N', 'N', 1.0_dp, matrix_c_im, qs_ot_env%rot_mat_u, &
2121 2150 : 0.0_dp, rotation_tmp)
2122 2150 : CALL dbcsr_add(rotated_im, rotation_tmp, alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
2123 :
2124 2150 : CALL dbcsr_copy(matrix_c, rotated_re)
2125 2150 : CALL dbcsr_copy(matrix_c_im, rotated_im)
2126 2150 : CALL dbcsr_release(rotated_re)
2127 2150 : CALL dbcsr_release(rotated_im)
2128 2150 : CALL dbcsr_release(rotation_tmp)
2129 : END IF
2130 :
2131 3555 : DEALLOCATE (eigenvalues, inverse_sqrt)
2132 :
2133 3555 : CALL timestop(handle)
2134 7110 : 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 9048 : 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 4524 : CALL timeset(routineN, handle)
2156 :
2157 4524 : CALL dbcsr_get_info(P, nfullcols_total=k)
2158 4524 : 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 4524 : CALL dbcsr_multiply('N', 'N', 1.0_dp, P, P, 0.0_dp, P2) ! P^2
2197 4524 : 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 4524 : r = -180.0_dp/128.0_dp
2203 4524 : CALL dbcsr_add(T, P, alpha_scalar=0.0_dp, beta_scalar=r) ! T=-180/128*P
2204 4524 : r = 35.0_dp/128.0_dp
2205 4524 : CALL dbcsr_add(T, P2, alpha_scalar=1.0_dp, beta_scalar=r) ! T=T+35/128*P^2
2206 4524 : CALL dbcsr_multiply('N', 'N', 1.0_dp, T, P2, 0.0_dp, FY) ! Y=T*P^2
2207 4524 : 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 4524 : r = 378.0_dp/128.0_dp
2213 4524 : CALL dbcsr_add(FY, P2, alpha_scalar=1.0_dp, beta_scalar=r) ! Y=Y+378/128*P^2
2214 4524 : r = -420.0_dp/128.0_dp
2215 4524 : CALL dbcsr_add(FY, P, alpha_scalar=1.0_dp, beta_scalar=r) ! Y=Y-420/128*P
2216 4524 : r = 315.0_dp/128.0_dp
2217 4524 : CALL dbcsr_add_on_diag(FY, alpha=r) ! Y=Y+315/128*I
2218 : CASE DEFAULT
2219 4524 : CPABORT("This irac_order NYI")
2220 : END SELECT
2221 4524 : CALL timestop(handle)
2222 4524 : 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 6640 : 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 3320 : CALL timeset(routineN, handle)
2244 :
2245 3320 : CALL dbcsr_get_info(matrix_x, nfullrows_total=n, nfullcols_total=k)
2246 : !
2247 3320 : C => matrix_x ! NBsf*NOcc
2248 3320 : SC => matrix_sx ! NBsf*NOcc need to be up2date
2249 3320 : HC => matrix_hc ! NBsf*NOcc
2250 3320 : G => matrix_gx ! NBsf*NOcc
2251 3320 : CHC => qs_ot_env%buf1_k_k_sym ! buffer
2252 :
2253 3320 : 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 3320 : CALL dbcsr_multiply('T', 'N', 1.0_dp, C, HC_work, 0.0_dp, CHC)
2266 3320 : 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 3320 : CALL dbcsr_multiply('N', 'N', 1.0_dp, SC, CHC, 0.0_dp, G)
2273 3320 : 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 3320 : CALL dbcsr_add(G, HC_work, alpha_scalar=-1.0_dp, beta_scalar=1.0_dp)
2280 : !
2281 3320 : CALL timestop(handle)
2282 3320 : 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 4728 : 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 2364 : CALL timeset(routineN, handle)
2308 :
2309 2364 : CPASSERT(qs_ot_env%has_complex_kpoint_state)
2310 2364 : CPASSERT(ASSOCIATED(matrix_hc))
2311 2364 : CPASSERT(ASSOCIATED(matrix_hc_im))
2312 :
2313 2364 : f_re => qs_ot_env%matrix_ref_inv_sqrt
2314 2364 : f_im => qs_ot_env%matrix_ref_inv_sqrt_im
2315 2364 : sc_re => qs_ot_env%matrix_sx
2316 2364 : sc_im => qs_ot_env%matrix_sx_im
2317 2364 : hc_re => matrix_hc
2318 2364 : hc_im => matrix_hc_im
2319 2364 : hc_rotation_re => hc_re
2320 2364 : hc_rotation_im => hc_im
2321 2364 : IF (PRESENT(matrix_hc_rotation) .OR. PRESENT(matrix_hc_rotation_im)) THEN
2322 426 : CPASSERT(PRESENT(matrix_hc_rotation) .AND. PRESENT(matrix_hc_rotation_im))
2323 426 : CPASSERT(ASSOCIATED(matrix_hc_rotation))
2324 426 : CPASSERT(ASSOCIATED(matrix_hc_rotation_im))
2325 426 : hc_rotation_re => matrix_hc_rotation
2326 426 : hc_rotation_im => matrix_hc_rotation_im
2327 : END IF
2328 2364 : hc_work_re => hc_re
2329 2364 : hc_work_im => hc_im
2330 2364 : g_re => qs_ot_env%matrix_gx
2331 2364 : g_im => qs_ot_env%matrix_gx_im
2332 2364 : b_re => qs_ot_env%buf1_k_k_sym
2333 2364 : b_im => qs_ot_env%buf2_k_k_sym
2334 2364 : tmp_kk => qs_ot_env%buf3_k_k_sym
2335 2364 : q_re => qs_ot_env%buf1_n_k
2336 2364 : q_im => qs_ot_env%buf1_n_k_dp
2337 :
2338 2364 : 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 2364 : CALL dbcsr_multiply('N', 'N', 1.0_dp, qs_ot_env%matrix_x, f_re, 0.0_dp, q_re)
2342 2364 : CALL dbcsr_multiply('N', 'N', 1.0_dp, qs_ot_env%matrix_x_im, f_im, 0.0_dp, g_re)
2343 2364 : CALL dbcsr_add(q_re, g_re, alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
2344 :
2345 2364 : CALL dbcsr_multiply('N', 'N', 1.0_dp, qs_ot_env%matrix_x, f_im, 0.0_dp, q_im)
2346 2364 : CALL dbcsr_multiply('N', 'N', 1.0_dp, qs_ot_env%matrix_x_im, f_re, 0.0_dp, g_re)
2347 2364 : CALL dbcsr_add(q_im, g_re, alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
2348 :
2349 2364 : IF (qs_ot_env%settings%do_rotation) THEN
2350 1336 : 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 1336 : 0.0_dp, qs_ot_env%rot_mat_dedu)
2355 1336 : 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 1336 : 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 1336 : 0.0_dp, qs_ot_env%rot_mat_dedu_im)
2361 1336 : 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 1336 : alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
2364 1336 : CALL qs_ot_rot_mat_derivative_complex(qs_ot_env)
2365 :
2366 : ! The REF coordinate sees G_Q=G_C*U^H.
2367 1336 : CALL dbcsr_copy(hc_rot_re, hc_re, name="hc_rot_re")
2368 1336 : CALL dbcsr_copy(hc_rot_im, hc_im, name="hc_rot_im")
2369 1336 : 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 1336 : 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 1336 : 0.0_dp, tmp_nk)
2374 1336 : 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 1336 : 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 1336 : 0.0_dp, tmp_nk)
2380 1336 : CALL dbcsr_add(hc_rot_im, tmp_nk, alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
2381 1336 : hc_work_re => hc_rot_re
2382 1336 : hc_work_im => hc_rot_im
2383 : END IF
2384 :
2385 : ! B = Q^H*G_Q. For uniform fixed occupations B is Hermitian.
2386 2364 : CALL dbcsr_multiply('T', 'N', 1.0_dp, q_re, hc_work_re, 0.0_dp, b_re)
2387 2364 : CALL dbcsr_multiply('T', 'N', 1.0_dp, q_im, hc_work_im, 0.0_dp, tmp_kk)
2388 2364 : CALL dbcsr_add(b_re, tmp_kk, alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
2389 :
2390 2364 : CALL dbcsr_multiply('T', 'N', 1.0_dp, q_re, hc_work_im, 0.0_dp, b_im)
2391 2364 : CALL dbcsr_multiply('T', 'N', 1.0_dp, q_im, hc_work_re, 0.0_dp, tmp_kk)
2392 2364 : CALL dbcsr_add(b_im, tmp_kk, alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
2393 :
2394 2364 : 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 2364 : CALL dbcsr_multiply('N', 'N', 1.0_dp, sc_re, f_re, 0.0_dp, g_re)
2405 2364 : CALL dbcsr_multiply('N', 'N', 1.0_dp, sc_im, f_im, 0.0_dp, q_re)
2406 2364 : CALL dbcsr_add(g_re, q_re, alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
2407 :
2408 2364 : CALL dbcsr_multiply('N', 'N', 1.0_dp, sc_re, f_im, 0.0_dp, g_im)
2409 2364 : CALL dbcsr_multiply('N', 'N', 1.0_dp, sc_im, f_re, 0.0_dp, q_re)
2410 2364 : 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 2364 : CALL dbcsr_multiply('N', 'N', 1.0_dp, g_re, b_re, 0.0_dp, q_re)
2414 2364 : CALL dbcsr_multiply('N', 'N', 1.0_dp, g_im, b_im, 0.0_dp, q_im)
2415 2364 : CALL dbcsr_add(q_re, q_im, alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
2416 :
2417 2364 : CALL dbcsr_multiply('N', 'N', 1.0_dp, g_re, b_im, 0.0_dp, q_im)
2418 2364 : CALL dbcsr_multiply('N', 'N', 1.0_dp, g_im, b_re, 0.0_dp, g_re)
2419 2364 : CALL dbcsr_add(q_im, g_re, alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
2420 :
2421 2364 : CALL dbcsr_add(q_re, hc_work_re, alpha_scalar=-1.0_dp, beta_scalar=1.0_dp)
2422 2364 : 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 2364 : CALL dbcsr_multiply('N', 'N', 1.0_dp, q_re, f_re, 0.0_dp, g_re)
2426 2364 : CALL dbcsr_multiply('N', 'N', 1.0_dp, q_im, f_im, 0.0_dp, g_im)
2427 2364 : CALL dbcsr_add(g_re, g_im, alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
2428 :
2429 2364 : CALL dbcsr_multiply('N', 'N', 1.0_dp, q_re, f_im, 0.0_dp, g_im)
2430 2364 : CALL dbcsr_multiply('N', 'N', 1.0_dp, q_im, f_re, 0.0_dp, q_re)
2431 2364 : CALL dbcsr_add(g_im, q_re, alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
2432 :
2433 2364 : 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 1336 : sc_re, sc_im, g_re, g_im, qs_ot_env)
2436 : END IF
2437 :
2438 2364 : 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 2364 : IF (qs_ot_env%settings%do_rotation) THEN
2448 1336 : CALL dbcsr_release(hc_rot_re)
2449 1336 : CALL dbcsr_release(hc_rot_im)
2450 1336 : CALL dbcsr_release(tmp_nk)
2451 : END IF
2452 :
2453 2364 : CALL timestop(handle)
2454 :
2455 2364 : 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 9634 : 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 9634 : CALL dbcsr_copy(transposed, matrix, name='square_transposed')
2469 9634 : CALL dbcsr_copy(identity, identity_template, name='transpose_identity')
2470 9634 : CALL dbcsr_set(identity, 0.0_dp)
2471 9634 : CALL dbcsr_add_on_diag(identity, 1.0_dp)
2472 9634 : CALL dbcsr_multiply('T', 'N', 1.0_dp, matrix, identity, 0.0_dp, transposed)
2473 9634 : CALL dbcsr_release(identity)
2474 :
2475 9634 : 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 1336 : 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 1336 : INTEGER, DIMENSION(:), POINTER :: col_blk_offset, row_blk_offset
2500 : LOGICAL :: found_im, found_re
2501 : REAL(KIND=dp) :: denominator
2502 1336 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: f_evals
2503 1336 : 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 1336 : CALL timeset(routineN, handle)
2509 :
2510 1336 : CALL dbcsr_get_info(a_re, nfullrows_total=k)
2511 1336 : 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 1336 : CALL qs_ot_square_transpose(a_re, tmp_kk, qs_ot_env%rot_mat_u)
2519 1336 : CALL dbcsr_copy(anti_re, a_re, name='anti_re')
2520 1336 : CALL dbcsr_add(anti_re, tmp_kk, alpha_scalar=0.5_dp, beta_scalar=-0.5_dp)
2521 1336 : CALL qs_ot_square_transpose(a_im, tmp_kk, qs_ot_env%rot_mat_u)
2522 1336 : CALL dbcsr_copy(anti_im, a_im, name='anti_im')
2523 1336 : 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 1336 : CALL dbcsr_copy(f_work_re, f_re, name='f_work_re')
2528 1336 : CALL dbcsr_copy(f_work_im, f_im, name='f_work_im')
2529 1336 : CALL dbcsr_copy(v_re, f_re, name='v_re')
2530 1336 : CALL dbcsr_copy(v_im, f_im, name='v_im')
2531 4008 : 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 1336 : para_env=qs_ot_env%para_env, blacs_env=qs_ot_env%blacs_env)
2535 9616 : 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 1336 : CALL dbcsr_copy(work_re, anti_re, name='work_re')
2541 1336 : CALL dbcsr_copy(work_im, anti_im, name='work_im')
2542 1336 : CALL dbcsr_multiply('N', 'N', 1.0_dp, anti_re, v_re, 0.0_dp, work_re)
2543 1336 : CALL dbcsr_multiply('N', 'N', 1.0_dp, anti_im, v_im, 0.0_dp, tmp_kk)
2544 1336 : CALL dbcsr_add(work_re, tmp_kk, alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
2545 1336 : CALL dbcsr_multiply('N', 'N', 1.0_dp, anti_re, v_im, 0.0_dp, work_im)
2546 1336 : CALL dbcsr_multiply('N', 'N', 1.0_dp, anti_im, v_re, 0.0_dp, tmp_kk)
2547 1336 : CALL dbcsr_add(work_im, tmp_kk, alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
2548 :
2549 1336 : CALL dbcsr_copy(inner_re, anti_re, name='inner_re')
2550 1336 : CALL dbcsr_copy(inner_im, anti_im, name='inner_im')
2551 1336 : CALL dbcsr_multiply('T', 'N', 1.0_dp, v_re, work_re, 0.0_dp, inner_re)
2552 1336 : CALL dbcsr_multiply('T', 'N', 1.0_dp, v_im, work_im, 1.0_dp, inner_re)
2553 1336 : CALL dbcsr_multiply('T', 'N', 1.0_dp, v_re, work_im, 0.0_dp, inner_im)
2554 1336 : 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 1336 : col_blk_offset=col_blk_offset)
2558 1336 : CALL dbcsr_iterator_start(iter, inner_re)
2559 2149 : DO WHILE (dbcsr_iterator_blocks_left(iter))
2560 813 : CALL dbcsr_iterator_next_block(iter, row, col)
2561 813 : CALL dbcsr_get_block_p(inner_re, row, col, block_re, found_re)
2562 813 : CALL dbcsr_get_block_p(inner_im, row, col, block_im, found_im)
2563 813 : CPASSERT(found_re .AND. found_im)
2564 7855 : DO i = 1, SIZE(block_re, 1)
2565 69391 : DO j = 1, SIZE(block_re, 2)
2566 : denominator = 1.0_dp/f_evals(row_blk_offset(row) + i - 1) + &
2567 62872 : 1.0_dp/f_evals(col_blk_offset(col) + j - 1)
2568 62872 : block_re(i, j) = block_re(i, j)/denominator
2569 68578 : block_im(i, j) = block_im(i, j)/denominator
2570 : END DO
2571 : END DO
2572 : END DO
2573 1336 : CALL dbcsr_iterator_stop(iter)
2574 :
2575 : ! Z = V inner V^H.
2576 1336 : CALL dbcsr_multiply('N', 'N', 1.0_dp, v_re, inner_re, 0.0_dp, work_re)
2577 1336 : CALL dbcsr_multiply('N', 'N', 1.0_dp, v_im, inner_im, 0.0_dp, tmp_kk)
2578 1336 : CALL dbcsr_add(work_re, tmp_kk, alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
2579 1336 : CALL dbcsr_multiply('N', 'N', 1.0_dp, v_re, inner_im, 0.0_dp, work_im)
2580 1336 : CALL dbcsr_multiply('N', 'N', 1.0_dp, v_im, inner_re, 0.0_dp, tmp_kk)
2581 1336 : CALL dbcsr_add(work_im, tmp_kk, alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
2582 :
2583 1336 : CALL dbcsr_copy(z_re, anti_re, name='z_re')
2584 1336 : CALL dbcsr_copy(z_im, anti_im, name='z_im')
2585 1336 : CALL dbcsr_multiply('N', 'T', 1.0_dp, work_re, v_re, 0.0_dp, z_re)
2586 1336 : CALL dbcsr_multiply('N', 'T', 1.0_dp, work_im, v_im, 1.0_dp, z_re)
2587 1336 : CALL dbcsr_multiply('N', 'T', 1.0_dp, work_im, v_re, 0.0_dp, z_im)
2588 1336 : 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 1336 : CALL dbcsr_copy(sq_re, sx_re, name='sq_re')
2592 1336 : CALL dbcsr_copy(sq_im, sx_im, name='sq_im')
2593 1336 : CALL dbcsr_copy(tmp_nk, sx_re, name='tmp_nk')
2594 1336 : CALL dbcsr_multiply('N', 'N', 1.0_dp, sx_re, f_re, 0.0_dp, sq_re)
2595 1336 : CALL dbcsr_multiply('N', 'N', 1.0_dp, sx_im, f_im, 0.0_dp, tmp_nk)
2596 1336 : CALL dbcsr_add(sq_re, tmp_nk, alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
2597 1336 : CALL dbcsr_multiply('N', 'N', 1.0_dp, sx_re, f_im, 0.0_dp, sq_im)
2598 1336 : CALL dbcsr_multiply('N', 'N', 1.0_dp, sx_im, f_re, 0.0_dp, tmp_nk)
2599 1336 : CALL dbcsr_add(sq_im, tmp_nk, alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
2600 :
2601 1336 : CALL dbcsr_copy(vertical_re, sx_re, name='vertical_re')
2602 1336 : CALL dbcsr_copy(vertical_im, sx_im, name='vertical_im')
2603 1336 : CALL dbcsr_multiply('N', 'N', 1.0_dp, sq_re, z_re, 0.0_dp, vertical_re)
2604 1336 : CALL dbcsr_multiply('N', 'N', 1.0_dp, sq_im, z_im, 0.0_dp, tmp_nk)
2605 1336 : CALL dbcsr_add(vertical_re, tmp_nk, alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
2606 1336 : CALL dbcsr_multiply('N', 'N', 1.0_dp, sq_re, z_im, 0.0_dp, vertical_im)
2607 1336 : CALL dbcsr_multiply('N', 'N', 1.0_dp, sq_im, z_re, 0.0_dp, tmp_nk)
2608 1336 : CALL dbcsr_add(vertical_im, tmp_nk, alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
2609 1336 : CALL dbcsr_add(gradient_re, vertical_re, alpha_scalar=1.0_dp, beta_scalar=2.0_dp)
2610 1336 : CALL dbcsr_add(gradient_im, vertical_im, alpha_scalar=1.0_dp, beta_scalar=2.0_dp)
2611 :
2612 1336 : DEALLOCATE (f_evals)
2613 1336 : CALL dbcsr_release(anti_im)
2614 1336 : CALL dbcsr_release(anti_re)
2615 1336 : CALL dbcsr_release(f_work_im)
2616 1336 : CALL dbcsr_release(f_work_re)
2617 1336 : CALL dbcsr_release(inner_im)
2618 1336 : CALL dbcsr_release(inner_re)
2619 1336 : CALL dbcsr_release(sq_im)
2620 1336 : CALL dbcsr_release(sq_re)
2621 1336 : CALL dbcsr_release(tmp_kk)
2622 1336 : CALL dbcsr_release(tmp_nk)
2623 1336 : CALL dbcsr_release(v_im)
2624 1336 : CALL dbcsr_release(v_re)
2625 1336 : CALL dbcsr_release(vertical_im)
2626 1336 : CALL dbcsr_release(vertical_re)
2627 1336 : CALL dbcsr_release(work_im)
2628 1336 : CALL dbcsr_release(work_re)
2629 1336 : CALL dbcsr_release(z_im)
2630 1336 : CALL dbcsr_release(z_re)
2631 :
2632 1336 : CALL timestop(handle)
2633 :
2634 2672 : 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 329283 : 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 109761 : CALL timeset(routineN, handle)
2655 :
2656 109761 : 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 109761 : 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 109761 : max_iter = 30; threshold = 1.0E-03_dp
2665 109761 : CALL arnoldi_extremal(qs_ot_env%matrix_p, max_ev, min_ev, converged, threshold, max_iter)
2666 109761 : qs_ot_env%largest_eval_upper_bound = MAX(max_ev, ABS(min_ev))
2667 :
2668 109761 : IF (.NOT. converged) qs_ot_env%largest_eval_upper_bound = dbcsr_gershgorin_norm(qs_ot_env%matrix_p)
2669 109761 : CALL decide_strategy(qs_ot_env)
2670 109761 : IF (qs_ot_env%do_taylor) THEN
2671 58064 : CALL qs_ot_p2m_taylor(qs_ot_env)
2672 : ELSE
2673 51697 : CALL qs_ot_p2m_diag(qs_ot_env)
2674 : END IF
2675 :
2676 109761 : IF (qs_ot_env%settings%do_rotation) THEN
2677 3922 : CALL qs_ot_generate_rotation(qs_ot_env)
2678 : END IF
2679 :
2680 109761 : CALL timestop(handle)
2681 :
2682 109761 : 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 5670 : 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 5670 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: exp_evals_im, exp_evals_re
2698 : TYPE(dbcsr_type) :: h_re, tmp, w_im, w_re
2699 :
2700 5670 : CALL timeset(routineN, handle)
2701 :
2702 5670 : CPASSERT(qs_ot_env%has_complex_kpoint_state)
2703 5670 : CPASSERT(ASSOCIATED(qs_ot_env%rot_mat_x_im))
2704 5670 : CPASSERT(ASSOCIATED(qs_ot_env%rot_mat_u_im))
2705 :
2706 5670 : CALL dbcsr_get_info(qs_ot_env%rot_mat_x, nfullrows_total=k)
2707 5670 : IF (k /= 0) THEN
2708 : rot_norm = dbcsr_frobenius_norm(qs_ot_env%rot_mat_x) + &
2709 5670 : dbcsr_frobenius_norm(qs_ot_env%rot_mat_x_im)
2710 5670 : IF (rot_norm <= EPSILON(1.0_dp)) THEN
2711 1020 : CALL dbcsr_set(qs_ot_env%rot_mat_u, 0.0_dp)
2712 1020 : CALL dbcsr_add_on_diag(qs_ot_env%rot_mat_u, 1.0_dp)
2713 1020 : CALL dbcsr_set(qs_ot_env%rot_mat_u_im, 0.0_dp)
2714 1020 : 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 4650 : CALL dbcsr_copy(h_re, qs_ot_env%rot_mat_x_im, name="h_re")
2721 4650 : 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 4650 : para_env=qs_ot_env%para_env, blacs_env=qs_ot_env%blacs_env)
2727 :
2728 18600 : ALLOCATE (exp_evals_re(k), exp_evals_im(k))
2729 40742 : exp_evals_re(:) = COS(-qs_ot_env%rot_mat_evals(:))
2730 40742 : exp_evals_im(:) = SIN(-qs_ot_env%rot_mat_evals(:))
2731 :
2732 4650 : CALL dbcsr_copy(w_re, qs_ot_env%rot_mat_evec_re, name="w_re")
2733 4650 : CALL dbcsr_scale_by_vector(w_re, alpha=exp_evals_re, side='right')
2734 4650 : CALL dbcsr_copy(tmp, qs_ot_env%rot_mat_evec_im, name="tmp")
2735 4650 : CALL dbcsr_scale_by_vector(tmp, alpha=exp_evals_im, side='right')
2736 4650 : CALL dbcsr_add(w_re, tmp, alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
2737 :
2738 4650 : CALL dbcsr_copy(w_im, qs_ot_env%rot_mat_evec_im, name="w_im")
2739 4650 : CALL dbcsr_scale_by_vector(w_im, alpha=exp_evals_re, side='right')
2740 4650 : CALL dbcsr_copy(tmp, qs_ot_env%rot_mat_evec_re)
2741 4650 : CALL dbcsr_scale_by_vector(tmp, alpha=exp_evals_im, side='right')
2742 4650 : 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 4650 : 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 4650 : 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 4650 : 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 4650 : 1.0_dp, qs_ot_env%rot_mat_u_im)
2752 :
2753 4650 : CALL dbcsr_release(h_re)
2754 4650 : CALL dbcsr_release(tmp)
2755 4650 : CALL dbcsr_release(w_re)
2756 4650 : CALL dbcsr_release(w_im)
2757 4650 : DEALLOCATE (exp_evals_re, exp_evals_im)
2758 : END IF
2759 :
2760 4650 : CALL timestop(handle)
2761 :
2762 10320 : 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 5104 : 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 2552 : CALL timeset(routineN, handle)
2781 :
2782 2552 : CPASSERT(qs_ot_env%has_complex_kpoint_state)
2783 2552 : CPASSERT(ASSOCIATED(qs_ot_env%rot_mat_dedu_im))
2784 2552 : CPASSERT(ASSOCIATED(qs_ot_env%rot_mat_gx_im))
2785 :
2786 2552 : CALL dbcsr_get_info(qs_ot_env%rot_mat_u, nfullrows_total=k)
2787 2552 : IF (k /= 0) THEN
2788 : rot_norm = dbcsr_frobenius_norm(qs_ot_env%rot_mat_x) + &
2789 2552 : dbcsr_frobenius_norm(qs_ot_env%rot_mat_x_im)
2790 2552 : IF (rot_norm <= EPSILON(1.0_dp)) THEN
2791 515 : CALL qs_ot_square_transpose(qs_ot_env%rot_mat_dedu, tmp, qs_ot_env%rot_mat_u)
2792 515 : 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 515 : alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
2795 515 : CALL dbcsr_release(tmp)
2796 :
2797 515 : CALL qs_ot_square_transpose(qs_ot_env%rot_mat_dedu_im, tmp, qs_ot_env%rot_mat_u)
2798 515 : 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 515 : alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
2801 515 : CALL dbcsr_release(tmp)
2802 515 : CALL timestop(handle)
2803 515 : RETURN
2804 : END IF
2805 :
2806 2037 : CALL dbcsr_copy(work_re, qs_ot_env%rot_mat_dedu, name="work_re")
2807 2037 : CALL dbcsr_copy(work_im, qs_ot_env%rot_mat_dedu, name="work_im")
2808 2037 : CALL dbcsr_copy(tmp, qs_ot_env%rot_mat_dedu, name="tmp")
2809 2037 : CALL dbcsr_copy(inner_deriv_re, qs_ot_env%rot_mat_dedu, name="inner_deriv_re")
2810 2037 : 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 2037 : 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 2037 : qs_ot_env%rot_mat_evec_im, 0.0_dp, tmp)
2817 2037 : 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 2037 : 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 2037 : qs_ot_env%rot_mat_evec_re, 0.0_dp, tmp)
2822 2037 : 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 2037 : 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 2037 : 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 2037 : 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 2037 : 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 2037 : 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 2037 : 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 2037 : outer_deriv_im, 0.0_dp, tmp)
2840 2037 : 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 2037 : 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 2037 : outer_deriv_re, 0.0_dp, tmp)
2845 2037 : CALL dbcsr_add(work_im, tmp, alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
2846 :
2847 2037 : CALL dbcsr_copy(frechet_re, qs_ot_env%rot_mat_dedu, name="frechet_re")
2848 2037 : 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 2037 : 0.0_dp, frechet_re)
2851 : CALL dbcsr_multiply('N', 'T', 1.0_dp, work_im, qs_ot_env%rot_mat_evec_im, &
2852 2037 : 1.0_dp, frechet_re)
2853 : CALL dbcsr_multiply('N', 'T', 1.0_dp, work_im, qs_ot_env%rot_mat_evec_re, &
2854 2037 : 0.0_dp, frechet_im)
2855 : CALL dbcsr_multiply('N', 'T', -1.0_dp, work_re, qs_ot_env%rot_mat_evec_im, &
2856 2037 : 1.0_dp, frechet_im)
2857 :
2858 : ! Tangents satisfy X^T=-X and Y^T=Y for A=X+iY.
2859 2037 : CALL qs_ot_square_transpose(frechet_re, tmp, qs_ot_env%rot_mat_u)
2860 2037 : CALL dbcsr_copy(qs_ot_env%rot_mat_gx, frechet_re)
2861 : CALL dbcsr_add(qs_ot_env%rot_mat_gx, tmp, &
2862 2037 : alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
2863 2037 : CALL qs_ot_square_transpose(frechet_im, tmp, qs_ot_env%rot_mat_u)
2864 2037 : 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 2037 : alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
2867 :
2868 2037 : CALL dbcsr_release(frechet_re)
2869 2037 : CALL dbcsr_release(frechet_im)
2870 2037 : CALL dbcsr_release(inner_deriv_re)
2871 2037 : CALL dbcsr_release(inner_deriv_im)
2872 2037 : CALL dbcsr_release(outer_deriv_re)
2873 2037 : CALL dbcsr_release(outer_deriv_im)
2874 2037 : CALL dbcsr_release(tmp)
2875 2037 : CALL dbcsr_release(work_re)
2876 2037 : CALL dbcsr_release(work_im)
2877 : END IF
2878 :
2879 2037 : CALL timestop(handle)
2880 :
2881 2552 : 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 2008 : 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 2008 : CALL timeset(routineN, handle)
2901 :
2902 2008 : CPASSERT(qs_ot_env%has_complex_kpoint_state)
2903 2008 : 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 2008 : qs_ot_env%matrix_p, qs_ot_env%matrix_p_im, qs_ot_env%matrix_buf1)
2907 2008 : qs_ot_env%do_taylor = .FALSE.
2908 2008 : CALL qs_ot_p2m_diag_complex(qs_ot_env)
2909 :
2910 2008 : CALL timestop(handle)
2911 :
2912 2008 : 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 3938 : 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 3938 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: exp_evals_im, exp_evals_re
2930 : TYPE(dbcsr_type) :: buf_1, buf_2
2931 :
2932 3938 : CALL timeset(routineN, handle)
2933 :
2934 3938 : CALL dbcsr_get_info(qs_ot_env%rot_mat_x, nfullrows_total=k)
2935 :
2936 3938 : 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 3886 : blacs_env=qs_ot_env%blacs_env)
2947 :
2948 : ! Compute: exp_evals = EXP(-i*rot_mat_evals)
2949 15544 : ALLOCATE (exp_evals_re(k), exp_evals_im(k))
2950 22130 : exp_evals_re(:) = COS(-qs_ot_env%rot_mat_evals(:))
2951 22130 : 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 3886 : CALL dbcsr_copy(buf_1, qs_ot_env%rot_mat_evec_re, name="buf_1")
2956 3886 : CALL dbcsr_scale_by_vector(buf_1, alpha=exp_evals_re, side='right')
2957 3886 : CALL dbcsr_copy(buf_2, qs_ot_env%rot_mat_evec_im, name="buf_2")
2958 3886 : CALL dbcsr_scale_by_vector(buf_2, alpha=exp_evals_im, side='right')
2959 3886 : CALL dbcsr_add(buf_1, buf_2, alpha_scalar=+1.0_dp, beta_scalar=-1.0_dp)
2960 3886 : 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 3886 : CALL dbcsr_copy(buf_1, qs_ot_env%rot_mat_evec_im)
2963 3886 : CALL dbcsr_scale_by_vector(buf_1, alpha=exp_evals_re, side='right')
2964 3886 : CALL dbcsr_copy(buf_2, qs_ot_env%rot_mat_evec_re)
2965 3886 : CALL dbcsr_scale_by_vector(buf_2, alpha=exp_evals_im, side='right')
2966 3886 : CALL dbcsr_add(buf_1, buf_2, alpha_scalar=+1.0_dp, beta_scalar=+1.0_dp)
2967 3886 : 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 3886 : CALL dbcsr_release(buf_1)
2971 3886 : CALL dbcsr_release(buf_2)
2972 3886 : DEALLOCATE (exp_evals_re, exp_evals_im)
2973 : END IF
2974 :
2975 3938 : CALL timestop(handle)
2976 :
2977 7876 : 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 1884 : 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 1884 : INTEGER, ALLOCATABLE, DIMENSION(:) :: cols, rows
2994 1884 : 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 1884 : 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 1884 : CALL timeset(routineN, handle)
3006 :
3007 1884 : CALL dbcsr_get_info(qs_ot_env%rot_mat_u, nfullrows_total=k)
3008 1884 : IF (k /= 0) THEN
3009 1858 : 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 1858 : 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 1858 : CALL dbcsr_copy(inner_deriv_re, qs_ot_env%rot_mat_dedu, "inner_deriv_re") ! TODO just create
3015 1858 : CALL dbcsr_copy(inner_deriv_im, qs_ot_env%rot_mat_dedu, "inner_deriv_im") ! TODO just create
3016 :
3017 1858 : 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 1858 : CALL dbcsr_multiply('T', 'N', +1.0_dp, qs_ot_env%rot_mat_evec_im, mat_buf, 0.0_dp, inner_deriv_re)
3019 1858 : 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 1858 : 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 1858 : CALL dbcsr_multiply('T', 'N', +1.0_dp, qs_ot_env%rot_mat_evec_re, mat_buf, 1.0_dp, inner_deriv_re)
3023 1858 : 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 1858 : max_blocks = 0
3028 1858 : CALL dbcsr_iterator_start(iter, inner_deriv_re)
3029 2787 : DO WHILE (dbcsr_iterator_blocks_left(iter))
3030 929 : CALL dbcsr_iterator_next_block(iter, row, col)
3031 929 : max_blocks = max_blocks + 1
3032 : END DO
3033 1858 : CALL dbcsr_iterator_stop(iter)
3034 1858 : CALL dbcsr_iterator_start(iter, inner_deriv_im)
3035 2787 : DO WHILE (dbcsr_iterator_blocks_left(iter))
3036 929 : CALL dbcsr_iterator_next_block(iter, row, col)
3037 929 : max_blocks = max_blocks + 1
3038 : END DO
3039 1858 : CALL dbcsr_iterator_stop(iter)
3040 :
3041 7432 : ALLOCATE (rows(MAX(max_blocks, 1)), cols(MAX(max_blocks, 1)))
3042 1858 : nblocks = 0
3043 1858 : CALL dbcsr_iterator_start(iter, inner_deriv_re)
3044 2787 : DO WHILE (dbcsr_iterator_blocks_left(iter))
3045 929 : CALL dbcsr_iterator_next_block(iter, row, col)
3046 929 : duplicate = .FALSE.
3047 929 : DO iblock = 1, nblocks
3048 929 : duplicate = duplicate .OR. (rows(iblock) == row .AND. cols(iblock) == col)
3049 : END DO
3050 2787 : IF (.NOT. duplicate) THEN
3051 929 : nblocks = nblocks + 1
3052 929 : rows(nblocks) = row
3053 929 : cols(nblocks) = col
3054 : END IF
3055 : END DO
3056 1858 : CALL dbcsr_iterator_stop(iter)
3057 1858 : CALL dbcsr_iterator_start(iter, inner_deriv_im)
3058 2787 : DO WHILE (dbcsr_iterator_blocks_left(iter))
3059 929 : CALL dbcsr_iterator_next_block(iter, row, col)
3060 929 : duplicate = .FALSE.
3061 1858 : DO iblock = 1, nblocks
3062 1858 : duplicate = duplicate .OR. (rows(iblock) == row .AND. cols(iblock) == col)
3063 : END DO
3064 2787 : 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 1858 : CALL dbcsr_iterator_stop(iter)
3071 :
3072 : CALL dbcsr_get_info(inner_deriv_re, distribution=dist, row_blk_size=row_blk_size, &
3073 1858 : col_blk_size=col_blk_size)
3074 : CALL dbcsr_create(outer_deriv_re, "outer_deriv_re", dist, dbcsr_type_no_symmetry, &
3075 1858 : row_blk_size, col_blk_size)
3076 : CALL dbcsr_create(outer_deriv_im, "outer_deriv_im", dist, dbcsr_type_no_symmetry, &
3077 1858 : row_blk_size, col_blk_size)
3078 1858 : IF (nblocks > 0) THEN
3079 929 : CALL dbcsr_reserve_blocks(outer_deriv_re, rows=rows(1:nblocks), cols=cols(1:nblocks))
3080 929 : CALL dbcsr_reserve_blocks(outer_deriv_im, rows=rows(1:nblocks), cols=cols(1:nblocks))
3081 : END IF
3082 1858 : CALL dbcsr_finalize(outer_deriv_re)
3083 1858 : CALL dbcsr_finalize(outer_deriv_im)
3084 1858 : CALL dbcsr_set(outer_deriv_re, 0.0_dp)
3085 1858 : CALL dbcsr_set(outer_deriv_im, 0.0_dp)
3086 :
3087 1858 : CALL dbcsr_get_info(outer_deriv_re, row_blk_offset=row_blk_offset, col_blk_offset=col_blk_offset)
3088 1858 : CALL dbcsr_iterator_start(iter, outer_deriv_re)
3089 2787 : DO WHILE (dbcsr_iterator_blocks_left(iter))
3090 929 : CALL dbcsr_iterator_next_block(iter, row, col)
3091 929 : CALL dbcsr_get_block_p(inner_deriv_re, row, col, block_in_re, found_in_re)
3092 929 : CALL dbcsr_get_block_p(inner_deriv_im, row, col, block_in_im, found_in_im)
3093 929 : CALL dbcsr_get_block_p(outer_deriv_re, row, col, block_out_re, found_out_re)
3094 929 : CALL dbcsr_get_block_p(outer_deriv_im, row, col, block_out_im, found_out_im)
3095 929 : CPASSERT(found_out_re .AND. found_out_im)
3096 :
3097 7133 : DO i = 1, SIZE(block_out_re, 1)
3098 31503 : DO j = 1, SIZE(block_out_re, 2)
3099 26228 : e1 = qs_ot_env%rot_mat_evals(row_blk_offset(row) + i - 1)
3100 26228 : e2 = qs_ot_env%rot_mat_evals(col_blk_offset(col) + j - 1)
3101 26228 : re_part = 0.0_dp
3102 26228 : im_part = 0.0_dp
3103 26228 : IF (found_in_re) re_part = block_in_re(i, j)
3104 26228 : IF (found_in_im) im_part = block_in_im(i, j)
3105 26228 : cval_in = CMPLX(re_part, im_part, dp)
3106 26228 : cval_out = cval_in*cint(e1, e2)
3107 26228 : block_out_re(i, j) = REAL(cval_out)
3108 30574 : block_out_im(i, j) = AIMAG(cval_out)
3109 : END DO
3110 : END DO
3111 : END DO
3112 1858 : CALL dbcsr_iterator_stop(iter)
3113 1858 : DEALLOCATE (rows, cols)
3114 1858 : CALL dbcsr_release(inner_deriv_re)
3115 1858 : 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 1858 : CALL dbcsr_multiply('N', 'N', +1.0_dp, qs_ot_env%rot_mat_evec_re, outer_deriv_re, 0.0_dp, mat_buf)
3119 1858 : CALL dbcsr_multiply('N', 'N', -1.0_dp, qs_ot_env%rot_mat_evec_im, outer_deriv_im, 1.0_dp, mat_buf)
3120 1858 : 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 1858 : CALL dbcsr_multiply('N', 'N', +1.0_dp, qs_ot_env%rot_mat_evec_re, outer_deriv_im, 0.0_dp, mat_buf)
3123 1858 : CALL dbcsr_multiply('N', 'N', +1.0_dp, qs_ot_env%rot_mat_evec_im, outer_deriv_re, 1.0_dp, mat_buf)
3124 1858 : 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 1858 : CALL qs_ot_square_transpose(qs_ot_env%matrix_buf1, mat_buf, qs_ot_env%rot_mat_u)
3129 1858 : 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 1858 : alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
3132 :
3133 1858 : CALL dbcsr_release(mat_buf)
3134 1858 : CALL dbcsr_release(outer_deriv_re)
3135 13006 : CALL dbcsr_release(outer_deriv_im)
3136 : END IF
3137 3768 : CALL timestop(handle)
3138 : CONTAINS
3139 :
3140 : ! **************************************************************************************************
3141 : !> \brief ...
3142 : !> \param e1 ...
3143 : !> \param e2 ...
3144 : !> \return ...
3145 : ! **************************************************************************************************
3146 26228 : 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 26228 : l1 = (0.0_dp, -1.0_dp)*e1
3154 26228 : l2 = (0.0_dp, -1.0_dp)*e2
3155 26228 : IF (ABS(l1 - l2) > 0.5_dp) THEN
3156 1020 : cint = (EXP(l1) - EXP(l2))/(l1 - l2)
3157 : ELSE
3158 : x = 1.0_dp
3159 : cint = 0.0_dp
3160 428536 : DO I = 1, 16
3161 403328 : cint = cint + x
3162 428536 : x = x*(l1 - l2)/REAL(I + 1, KIND=dp)
3163 : END DO
3164 25208 : cint = cint*EXP(l2)
3165 : END IF
3166 26228 : 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 109761 : 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 109761 : qs_ot_env%do_taylor = .FALSE.
3185 109761 : N = 0
3186 109761 : num_error = qs_ot_env%largest_eval_upper_bound/(2.0_dp)
3187 470653 : DO WHILE (num_error > qs_ot_env%settings%eps_taylor .AND. N <= qs_ot_env%settings%max_taylor)
3188 360892 : N = N + 1
3189 407367 : num_error = num_error*qs_ot_env%largest_eval_upper_bound/REAL((2*N + 1)*(2*N + 2), KIND=dp)
3190 : END DO
3191 109761 : qs_ot_env%taylor_order = N
3192 109761 : IF (qs_ot_env%taylor_order <= qs_ot_env%settings%max_taylor) THEN
3193 58064 : qs_ot_env%do_taylor = .TRUE.
3194 : END IF
3195 :
3196 109761 : END SUBROUTINE decide_strategy
3197 :
3198 : ! **************************************************************************************************
3199 : !> \brief c=(c0*cos(p^0.5)+x*sin(p^0.5)*p^(-0.5)) x rot_mat_u
3200 : !> this assumes that x is already ortho to S*C0, and that p is x*S*x
3201 : !> rot_mat_u is an optional rotation matrix
3202 : !> \param matrix_c ...
3203 : !> \param matrix_x ...
3204 : !> \param qs_ot_env ...
3205 : ! **************************************************************************************************
3206 203906 : SUBROUTINE qs_ot_get_orbitals(matrix_c, matrix_x, qs_ot_env)
3207 :
3208 : TYPE(dbcsr_type), POINTER :: matrix_c, matrix_x
3209 : TYPE(qs_ot_type) :: qs_ot_env
3210 :
3211 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_ot_get_orbitals'
3212 : REAL(KIND=dp), PARAMETER :: rone = 1.0_dp, rzero = 0.0_dp
3213 :
3214 : INTEGER :: handle, k, n
3215 : TYPE(dbcsr_type), POINTER :: matrix_kk
3216 :
3217 101953 : CALL timeset(routineN, handle)
3218 :
3219 101953 : CALL dbcsr_get_info(matrix_x, nfullrows_total=n, nfullcols_total=k)
3220 :
3221 : ! rotate the multiplying matrices cosp and sinp instead of the result,
3222 : ! this should be cheaper for large basis sets
3223 101953 : IF (qs_ot_env%settings%do_rotation) THEN
3224 3690 : matrix_kk => qs_ot_env%matrix_buf1
3225 : CALL dbcsr_multiply('N', 'N', rone, qs_ot_env%matrix_cosp, &
3226 3690 : qs_ot_env%rot_mat_u, rzero, matrix_kk)
3227 : ELSE
3228 98263 : matrix_kk => qs_ot_env%matrix_cosp
3229 : END IF
3230 :
3231 : CALL dbcsr_multiply('N', 'N', rone, qs_ot_env%matrix_c0, matrix_kk, &
3232 101953 : rzero, matrix_c)
3233 :
3234 101953 : IF (qs_ot_env%settings%do_rotation) THEN
3235 3690 : matrix_kk => qs_ot_env%matrix_buf1
3236 : CALL dbcsr_multiply('N', 'N', rone, qs_ot_env%matrix_sinp, &
3237 3690 : qs_ot_env%rot_mat_u, rzero, matrix_kk)
3238 : ELSE
3239 98263 : matrix_kk => qs_ot_env%matrix_sinp
3240 : END IF
3241 : CALL dbcsr_multiply('N', 'N', rone, matrix_x, matrix_kk, &
3242 101953 : rone, matrix_c)
3243 :
3244 101953 : CALL timestop(handle)
3245 :
3246 101953 : END SUBROUTINE qs_ot_get_orbitals
3247 :
3248 : ! **************************************************************************************************
3249 : !> \brief update complex K-point orbitals with the finite STRICT transformation
3250 : !> \param matrix_c real output orbitals
3251 : !> \param matrix_c_im imaginary output orbitals
3252 : !> \param matrix_s real overlap matrix
3253 : !> \param matrix_s_im imaginary overlap matrix
3254 : !> \param qs_ot_env OT channel state
3255 : ! **************************************************************************************************
3256 1908 : SUBROUTINE qs_ot_get_orbitals_complex(matrix_c, matrix_c_im, matrix_s, matrix_s_im, qs_ot_env)
3257 : TYPE(dbcsr_type), POINTER :: matrix_c, matrix_c_im, matrix_s, &
3258 : matrix_s_im
3259 : TYPE(qs_ot_type) :: qs_ot_env
3260 :
3261 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_ot_get_orbitals_complex'
3262 :
3263 : INTEGER :: handle
3264 : TYPE(dbcsr_type) :: rotated_im, rotated_re, rotation_tmp
3265 :
3266 1908 : CALL timeset(routineN, handle)
3267 :
3268 1908 : CPASSERT(qs_ot_env%has_complex_kpoint_state)
3269 1908 : CPASSERT(qs_ot_env%settings%ot_algorithm == 'TOD')
3270 :
3271 : CALL qs_ot_complex_multiply('N', 'N', matrix_s, matrix_s_im, &
3272 : qs_ot_env%matrix_x, qs_ot_env%matrix_x_im, &
3273 : qs_ot_env%matrix_sx, qs_ot_env%matrix_sx_im, &
3274 1908 : qs_ot_env%matrix_tmp_nk)
3275 : CALL qs_ot_get_p_complex(qs_ot_env%matrix_x, qs_ot_env%matrix_x_im, &
3276 1908 : qs_ot_env%matrix_sx, qs_ot_env%matrix_sx_im, qs_ot_env)
3277 :
3278 : CALL qs_ot_complex_multiply('N', 'N', qs_ot_env%matrix_c0, qs_ot_env%matrix_c0_im, &
3279 : qs_ot_env%matrix_cosp, qs_ot_env%matrix_cosp_im, &
3280 1908 : matrix_c, matrix_c_im, qs_ot_env%matrix_tmp_nk)
3281 : CALL qs_ot_complex_multiply('N', 'N', qs_ot_env%matrix_x, qs_ot_env%matrix_x_im, &
3282 : qs_ot_env%matrix_sinp, qs_ot_env%matrix_sinp_im, &
3283 : qs_ot_env%matrix_buf_nk, qs_ot_env%matrix_buf_nk_im, &
3284 1908 : qs_ot_env%matrix_tmp_nk)
3285 1908 : CALL dbcsr_add(matrix_c, qs_ot_env%matrix_buf_nk, alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
3286 1908 : CALL dbcsr_add(matrix_c_im, qs_ot_env%matrix_buf_nk_im, alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
3287 :
3288 1908 : IF (qs_ot_env%settings%do_rotation) THEN
3289 968 : CALL qs_ot_generate_rotation_complex(qs_ot_env)
3290 :
3291 968 : CALL dbcsr_copy(rotated_re, matrix_c, name="strict_rotated_re")
3292 968 : CALL dbcsr_copy(rotated_im, matrix_c_im, name="strict_rotated_im")
3293 968 : CALL dbcsr_copy(rotation_tmp, matrix_c, name="strict_rotation_tmp")
3294 :
3295 : ! C_out = Q(X)*U for the finite STRICT chart Q(X).
3296 : CALL dbcsr_multiply('N', 'N', 1.0_dp, matrix_c, qs_ot_env%rot_mat_u, &
3297 968 : 0.0_dp, rotated_re)
3298 : CALL dbcsr_multiply('N', 'N', 1.0_dp, matrix_c_im, qs_ot_env%rot_mat_u_im, &
3299 968 : 0.0_dp, rotation_tmp)
3300 968 : CALL dbcsr_add(rotated_re, rotation_tmp, alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
3301 :
3302 : CALL dbcsr_multiply('N', 'N', 1.0_dp, matrix_c, qs_ot_env%rot_mat_u_im, &
3303 968 : 0.0_dp, rotated_im)
3304 : CALL dbcsr_multiply('N', 'N', 1.0_dp, matrix_c_im, qs_ot_env%rot_mat_u, &
3305 968 : 0.0_dp, rotation_tmp)
3306 968 : CALL dbcsr_add(rotated_im, rotation_tmp, alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
3307 :
3308 968 : CALL dbcsr_copy(matrix_c, rotated_re)
3309 968 : CALL dbcsr_copy(matrix_c_im, rotated_im)
3310 968 : CALL dbcsr_release(rotated_re)
3311 968 : CALL dbcsr_release(rotated_im)
3312 968 : CALL dbcsr_release(rotation_tmp)
3313 : END IF
3314 :
3315 1908 : CALL timestop(handle)
3316 :
3317 1908 : END SUBROUTINE qs_ot_get_orbitals_complex
3318 :
3319 : ! **************************************************************************************************
3320 : !> \brief this routines computes dE/dx=dx, with dx ortho to sc0
3321 : !> needs dE/dC=hc,C0,X,SX,p
3322 : !> if preconditioned it will not be the derivative, but the lagrangian multiplier
3323 : !> is changed so that P*dE/dx is the right derivative (i.e. in the allowed subspace)
3324 : !> \param matrix_hc ...
3325 : !> \param matrix_x ...
3326 : !> \param matrix_sx ...
3327 : !> \param matrix_gx ...
3328 : !> \param qs_ot_env ...
3329 : ! **************************************************************************************************
3330 233319 : SUBROUTINE qs_ot_get_derivative(matrix_hc, matrix_x, matrix_sx, matrix_gx, &
3331 : qs_ot_env)
3332 : TYPE(dbcsr_type), POINTER :: matrix_hc, matrix_x, matrix_sx, matrix_gx
3333 : TYPE(qs_ot_type) :: qs_ot_env
3334 :
3335 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_ot_get_derivative'
3336 : REAL(KIND=dp), PARAMETER :: rone = 1.0_dp, rzero = 0.0_dp
3337 :
3338 : INTEGER :: handle, k, n, ortho_k
3339 : TYPE(dbcsr_type), POINTER :: matrix_hc_local, matrix_target
3340 :
3341 77773 : CALL timeset(routineN, handle)
3342 :
3343 77773 : NULLIFY (matrix_hc_local)
3344 :
3345 77773 : CALL dbcsr_get_info(matrix_x, nfullrows_total=n, nfullcols_total=k)
3346 :
3347 : ! could in principle be taken inside qs_ot_get_derivative_* for increased efficiency
3348 : ! create a local rotated version of matrix_hc leaving matrix_hc untouched (needed
3349 : ! for lagrangian multipliers)
3350 77773 : IF (qs_ot_env%settings%do_rotation) THEN
3351 1874 : CALL dbcsr_copy(matrix_gx, matrix_hc) ! use gx as temporary
3352 1874 : CALL dbcsr_init_p(matrix_hc_local)
3353 1874 : CALL dbcsr_copy(matrix_hc_local, matrix_hc, name='matrix_hc_local')
3354 1874 : CALL dbcsr_set(matrix_hc_local, 0.0_dp)
3355 1874 : CALL dbcsr_multiply('N', 'T', rone, matrix_gx, qs_ot_env%rot_mat_u, rzero, matrix_hc_local)
3356 : ELSE
3357 75899 : matrix_hc_local => matrix_hc
3358 : END IF
3359 :
3360 77773 : IF (qs_ot_env%do_taylor) THEN
3361 42613 : CALL qs_ot_get_derivative_taylor(matrix_hc_local, matrix_x, matrix_sx, matrix_gx, qs_ot_env)
3362 : ELSE
3363 35160 : CALL qs_ot_get_derivative_diag(matrix_hc_local, matrix_x, matrix_sx, matrix_gx, qs_ot_env)
3364 : END IF
3365 :
3366 : ! and make it orthogonal
3367 77773 : CALL dbcsr_get_info(qs_ot_env%matrix_sc0, nfullcols_total=ortho_k)
3368 :
3369 77773 : IF (ASSOCIATED(qs_ot_env%preconditioner)) THEN
3370 67287 : matrix_target => qs_ot_env%matrix_psc0
3371 : ELSE
3372 10486 : matrix_target => qs_ot_env%matrix_sc0
3373 : END IF
3374 : ! first make the matrix os if not yet valid
3375 77773 : IF (.NOT. qs_ot_env%os_valid) THEN
3376 : ! this assumes that the preconditioner is a single matrix
3377 : ! that maps sc0 onto psc0
3378 :
3379 8452 : IF (ASSOCIATED(qs_ot_env%preconditioner)) THEN
3380 : CALL apply_preconditioner(qs_ot_env%preconditioner, qs_ot_env%matrix_sc0, &
3381 7470 : qs_ot_env%matrix_psc0)
3382 : END IF
3383 : CALL dbcsr_multiply('T', 'N', rone, &
3384 : qs_ot_env%matrix_sc0, matrix_target, &
3385 8452 : rzero, qs_ot_env%matrix_os)
3386 : CALL cp_dbcsr_cholesky_decompose(qs_ot_env%matrix_os, &
3387 8452 : para_env=qs_ot_env%para_env, blacs_env=qs_ot_env%blacs_env)
3388 : CALL cp_dbcsr_cholesky_invert(qs_ot_env%matrix_os, &
3389 : para_env=qs_ot_env%para_env, blacs_env=qs_ot_env%blacs_env, &
3390 8452 : uplo_to_full=.TRUE.)
3391 8452 : qs_ot_env%os_valid = .TRUE.
3392 : END IF
3393 : CALL dbcsr_multiply('T', 'N', rone, matrix_target, matrix_gx, &
3394 77773 : rzero, qs_ot_env%matrix_buf1_ortho)
3395 : CALL dbcsr_multiply('N', 'N', rone, qs_ot_env%matrix_os, &
3396 77773 : qs_ot_env%matrix_buf1_ortho, rzero, qs_ot_env%matrix_buf2_ortho)
3397 : CALL dbcsr_multiply('N', 'N', -rone, qs_ot_env%matrix_sc0, &
3398 77773 : qs_ot_env%matrix_buf2_ortho, rone, matrix_gx)
3399 : ! also treat the rot_mat gradient here
3400 77773 : IF (qs_ot_env%settings%do_rotation) THEN
3401 1874 : CALL qs_ot_rot_mat_derivative(qs_ot_env)
3402 : END IF
3403 :
3404 77773 : IF (qs_ot_env%settings%do_rotation) THEN
3405 1874 : CALL dbcsr_release_p(matrix_hc_local)
3406 : END IF
3407 :
3408 77773 : CALL timestop(handle)
3409 :
3410 77773 : END SUBROUTINE qs_ot_get_derivative
3411 :
3412 : ! **************************************************************************************************
3413 : !> \brief Prepare the inverse metric used to project a complex STRICT gradient.
3414 : !> An unusable preconditioner is detached before any minimizer history is updated.
3415 : !> \param qs_ot_env OT channel state
3416 : !> \param preconditioner_rejected true if the attached preconditioner was not positive definite
3417 : ! **************************************************************************************************
3418 2176 : SUBROUTINE qs_ot_prepare_complex_tangent_metric(qs_ot_env, preconditioner_rejected)
3419 : TYPE(qs_ot_type) :: qs_ot_env
3420 : LOGICAL, INTENT(OUT), OPTIONAL :: preconditioner_rejected
3421 :
3422 : INTEGER :: i, k
3423 : REAL(KIND=dp) :: eval_scale, eval_threshold
3424 : TYPE(dbcsr_type), POINTER :: target_im, target_re
3425 :
3426 1976 : IF (PRESENT(preconditioner_rejected)) preconditioner_rejected = .FALSE.
3427 1976 : IF (qs_ot_env%os_valid) RETURN
3428 :
3429 200 : IF (ASSOCIATED(qs_ot_env%preconditioner)) THEN
3430 190 : target_re => qs_ot_env%matrix_psc0
3431 190 : target_im => qs_ot_env%matrix_psc0_im
3432 : CALL apply_preconditioner(qs_ot_env%preconditioner, &
3433 : qs_ot_env%matrix_sc0, qs_ot_env%matrix_sc0_im, &
3434 190 : target_re, target_im)
3435 : ELSE
3436 10 : target_re => qs_ot_env%matrix_sc0
3437 10 : target_im => qs_ot_env%matrix_sc0_im
3438 : END IF
3439 : CALL qs_ot_complex_multiply('C', 'N', qs_ot_env%matrix_sc0, qs_ot_env%matrix_sc0_im, &
3440 : target_re, target_im, qs_ot_env%matrix_os, qs_ot_env%matrix_os_im, &
3441 200 : qs_ot_env%matrix_buf1)
3442 200 : CALL dbcsr_get_info(qs_ot_env%matrix_os, nfullrows_total=k)
3443 : CALL cp_dbcsr_heevd(matrix_re=qs_ot_env%matrix_os, matrix_im=qs_ot_env%matrix_os_im, &
3444 : eigenvectors_re=qs_ot_env%matrix_r, eigenvectors_im=qs_ot_env%matrix_r_im, &
3445 : eigenvalues=qs_ot_env%evals, para_env=qs_ot_env%para_env, &
3446 200 : blacs_env=qs_ot_env%blacs_env)
3447 1596 : eval_scale = MAX(1.0_dp, MAXVAL(ABS(qs_ot_env%evals(1:k))))
3448 200 : eval_threshold = 100.0_dp*EPSILON(1.0_dp)*eval_scale
3449 1596 : IF (MINVAL(qs_ot_env%evals(1:k)) <= eval_threshold .AND. &
3450 : ASSOCIATED(qs_ot_env%preconditioner)) THEN
3451 0 : NULLIFY (qs_ot_env%preconditioner)
3452 0 : IF (PRESENT(preconditioner_rejected)) preconditioner_rejected = .TRUE.
3453 0 : target_re => qs_ot_env%matrix_sc0
3454 0 : target_im => qs_ot_env%matrix_sc0_im
3455 : CALL qs_ot_complex_multiply('C', 'N', qs_ot_env%matrix_sc0, qs_ot_env%matrix_sc0_im, &
3456 : target_re, target_im, qs_ot_env%matrix_os, qs_ot_env%matrix_os_im, &
3457 0 : qs_ot_env%matrix_buf1)
3458 : CALL cp_dbcsr_heevd(matrix_re=qs_ot_env%matrix_os, matrix_im=qs_ot_env%matrix_os_im, &
3459 : eigenvectors_re=qs_ot_env%matrix_r, eigenvectors_im=qs_ot_env%matrix_r_im, &
3460 : eigenvalues=qs_ot_env%evals, para_env=qs_ot_env%para_env, &
3461 0 : blacs_env=qs_ot_env%blacs_env)
3462 0 : eval_scale = MAX(1.0_dp, MAXVAL(ABS(qs_ot_env%evals(1:k))))
3463 0 : eval_threshold = 100.0_dp*EPSILON(1.0_dp)*eval_scale
3464 : END IF
3465 1596 : CPASSERT(MINVAL(qs_ot_env%evals(1:k)) > eval_threshold)
3466 1596 : DO i = 1, k
3467 1596 : qs_ot_env%dum(i) = 1.0_dp/qs_ot_env%evals(i)
3468 : END DO
3469 200 : CALL dbcsr_copy(qs_ot_env%matrix_buf1, qs_ot_env%matrix_r)
3470 200 : CALL dbcsr_copy(qs_ot_env%matrix_buf1_im, qs_ot_env%matrix_r_im)
3471 200 : CALL dbcsr_scale_by_vector(qs_ot_env%matrix_buf1, alpha=qs_ot_env%dum, side='right')
3472 200 : CALL dbcsr_scale_by_vector(qs_ot_env%matrix_buf1_im, alpha=qs_ot_env%dum, side='right')
3473 : CALL qs_ot_complex_multiply('N', 'C', qs_ot_env%matrix_buf1, qs_ot_env%matrix_buf1_im, &
3474 : qs_ot_env%matrix_r, qs_ot_env%matrix_r_im, &
3475 : qs_ot_env%matrix_os, qs_ot_env%matrix_os_im, &
3476 200 : qs_ot_env%matrix_buf2)
3477 200 : qs_ot_env%os_valid = .TRUE.
3478 :
3479 : END SUBROUTINE qs_ot_prepare_complex_tangent_metric
3480 :
3481 : ! **************************************************************************************************
3482 : !> \brief finite complex STRICT derivative, projected onto C0^H*S*X=0
3483 : !> \param matrix_hc real part of H(k)*C(k)
3484 : !> \param matrix_hc_im imaginary part of H(k)*C(k)
3485 : !> \param qs_ot_env OT channel state
3486 : !> \param matrix_hc_rotation occupation-weighted H(k)C(k) for the rotation channel
3487 : !> \param matrix_hc_rotation_im imaginary component of matrix_hc_rotation
3488 : ! **************************************************************************************************
3489 3752 : SUBROUTINE qs_ot_get_derivative_complex(matrix_hc, matrix_hc_im, qs_ot_env, &
3490 : matrix_hc_rotation, matrix_hc_rotation_im)
3491 : TYPE(dbcsr_type), POINTER :: matrix_hc, matrix_hc_im
3492 : TYPE(qs_ot_type) :: qs_ot_env
3493 : TYPE(dbcsr_type), OPTIONAL, POINTER :: matrix_hc_rotation, matrix_hc_rotation_im
3494 :
3495 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_ot_get_derivative_complex'
3496 :
3497 : INTEGER :: handle
3498 : TYPE(dbcsr_distribution_type) :: dist
3499 : TYPE(dbcsr_type) :: tmp_nk
3500 : TYPE(dbcsr_type), POINTER :: hc_rotation_im, hc_rotation_re, &
3501 : hc_work_im, hc_work_re, target_im, &
3502 : target_re
3503 : TYPE(dbcsr_type), TARGET :: hc_rot_im, hc_rot_re
3504 :
3505 1876 : CALL timeset(routineN, handle)
3506 :
3507 1876 : CPASSERT(qs_ot_env%has_complex_kpoint_state)
3508 1876 : CPASSERT(qs_ot_env%settings%ot_algorithm == 'TOD')
3509 :
3510 1876 : hc_rotation_re => matrix_hc
3511 1876 : hc_rotation_im => matrix_hc_im
3512 1876 : IF (PRESENT(matrix_hc_rotation) .OR. PRESENT(matrix_hc_rotation_im)) THEN
3513 568 : CPASSERT(PRESENT(matrix_hc_rotation) .AND. PRESENT(matrix_hc_rotation_im))
3514 568 : CPASSERT(ASSOCIATED(matrix_hc_rotation))
3515 568 : CPASSERT(ASSOCIATED(matrix_hc_rotation_im))
3516 568 : hc_rotation_re => matrix_hc_rotation
3517 568 : hc_rotation_im => matrix_hc_rotation_im
3518 : END IF
3519 1876 : hc_work_re => matrix_hc
3520 1876 : hc_work_im => matrix_hc_im
3521 :
3522 1876 : IF (qs_ot_env%settings%do_rotation) THEN
3523 1216 : CALL qs_ot_generate_rotation_complex(qs_ot_env)
3524 :
3525 : ! Reconstruct the unrotated finite STRICT orbitals Q(X).
3526 : CALL qs_ot_complex_multiply('N', 'N', qs_ot_env%matrix_c0, qs_ot_env%matrix_c0_im, &
3527 : qs_ot_env%matrix_cosp, qs_ot_env%matrix_cosp_im, &
3528 : qs_ot_env%matrix_buf_nk, qs_ot_env%matrix_buf_nk_im, &
3529 1216 : qs_ot_env%matrix_tmp_nk)
3530 : CALL qs_ot_complex_multiply('N', 'N', qs_ot_env%matrix_x, qs_ot_env%matrix_x_im, &
3531 : qs_ot_env%matrix_sinp, qs_ot_env%matrix_sinp_im, &
3532 : qs_ot_env%matrix_gx, qs_ot_env%matrix_gx_im, &
3533 1216 : qs_ot_env%matrix_tmp_nk)
3534 : CALL dbcsr_add(qs_ot_env%matrix_buf_nk, qs_ot_env%matrix_gx, &
3535 1216 : alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
3536 : CALL dbcsr_add(qs_ot_env%matrix_buf_nk_im, qs_ot_env%matrix_gx_im, &
3537 1216 : alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
3538 :
3539 : ! dF/dU = Q(X)^H*G_C for C=Q(X)*U.
3540 : CALL qs_ot_complex_multiply('C', 'N', &
3541 : qs_ot_env%matrix_buf_nk, qs_ot_env%matrix_buf_nk_im, &
3542 : hc_rotation_re, hc_rotation_im, &
3543 : qs_ot_env%rot_mat_dedu, qs_ot_env%rot_mat_dedu_im, &
3544 1216 : qs_ot_env%matrix_buf1)
3545 1216 : CALL qs_ot_rot_mat_derivative_complex(qs_ot_env)
3546 :
3547 : ! The STRICT coordinate sees G_Q=G_C*U^H.
3548 1216 : CALL dbcsr_copy(hc_rot_re, matrix_hc, name="strict_hc_rot_re")
3549 1216 : CALL dbcsr_copy(hc_rot_im, matrix_hc_im, name="strict_hc_rot_im")
3550 1216 : CALL dbcsr_copy(tmp_nk, matrix_hc, name="strict_hc_rot_tmp")
3551 : CALL dbcsr_multiply('N', 'T', 1.0_dp, matrix_hc, qs_ot_env%rot_mat_u, &
3552 1216 : 0.0_dp, hc_rot_re)
3553 : CALL dbcsr_multiply('N', 'T', 1.0_dp, matrix_hc_im, qs_ot_env%rot_mat_u_im, &
3554 1216 : 0.0_dp, tmp_nk)
3555 1216 : CALL dbcsr_add(hc_rot_re, tmp_nk, alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
3556 :
3557 : CALL dbcsr_multiply('N', 'T', 1.0_dp, matrix_hc_im, qs_ot_env%rot_mat_u, &
3558 1216 : 0.0_dp, hc_rot_im)
3559 : CALL dbcsr_multiply('N', 'T', 1.0_dp, matrix_hc, qs_ot_env%rot_mat_u_im, &
3560 1216 : 0.0_dp, tmp_nk)
3561 1216 : CALL dbcsr_add(hc_rot_im, tmp_nk, alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
3562 1216 : hc_work_re => hc_rot_re
3563 1216 : hc_work_im => hc_rot_im
3564 : END IF
3565 :
3566 : ! Direct X contribution, H*C sinc(sqrt(P)).
3567 : CALL qs_ot_complex_multiply('N', 'N', hc_work_re, hc_work_im, &
3568 : qs_ot_env%matrix_sinp, qs_ot_env%matrix_sinp_im, &
3569 : qs_ot_env%matrix_gx, qs_ot_env%matrix_gx_im, &
3570 1876 : qs_ot_env%matrix_tmp_nk)
3571 :
3572 : ! Frechet contribution from X sinc(sqrt(P)).
3573 : CALL qs_ot_complex_multiply('C', 'N', hc_work_re, hc_work_im, &
3574 : qs_ot_env%matrix_x, qs_ot_env%matrix_x_im, &
3575 : qs_ot_env%matrix_buf2, qs_ot_env%matrix_buf2_im, &
3576 1876 : qs_ot_env%matrix_buf1)
3577 : CALL qs_ot_complex_multiply('N', 'N', qs_ot_env%matrix_buf2, qs_ot_env%matrix_buf2_im, &
3578 : qs_ot_env%matrix_r, qs_ot_env%matrix_r_im, &
3579 : qs_ot_env%matrix_buf1, qs_ot_env%matrix_buf1_im, &
3580 1876 : qs_ot_env%matrix_buf4)
3581 : CALL qs_ot_complex_multiply('C', 'N', qs_ot_env%matrix_r, qs_ot_env%matrix_r_im, &
3582 : qs_ot_env%matrix_buf1, qs_ot_env%matrix_buf1_im, &
3583 : qs_ot_env%matrix_buf2, qs_ot_env%matrix_buf2_im, &
3584 1876 : qs_ot_env%matrix_buf4)
3585 : CALL dbcsr_hadamard_product(qs_ot_env%matrix_buf2, qs_ot_env%matrix_sinp_b, &
3586 1876 : qs_ot_env%matrix_buf3)
3587 : CALL dbcsr_hadamard_product(qs_ot_env%matrix_buf2_im, qs_ot_env%matrix_sinp_b, &
3588 1876 : qs_ot_env%matrix_buf3_im)
3589 :
3590 : ! Frechet contribution from C0 cos(sqrt(P)).
3591 : CALL qs_ot_complex_multiply('C', 'N', hc_work_re, hc_work_im, &
3592 : qs_ot_env%matrix_c0, qs_ot_env%matrix_c0_im, &
3593 : qs_ot_env%matrix_buf2, qs_ot_env%matrix_buf2_im, &
3594 1876 : qs_ot_env%matrix_buf1)
3595 : CALL qs_ot_complex_multiply('N', 'N', qs_ot_env%matrix_buf2, qs_ot_env%matrix_buf2_im, &
3596 : qs_ot_env%matrix_r, qs_ot_env%matrix_r_im, &
3597 : qs_ot_env%matrix_buf1, qs_ot_env%matrix_buf1_im, &
3598 1876 : qs_ot_env%matrix_buf4)
3599 : CALL qs_ot_complex_multiply('C', 'N', qs_ot_env%matrix_r, qs_ot_env%matrix_r_im, &
3600 : qs_ot_env%matrix_buf1, qs_ot_env%matrix_buf1_im, &
3601 : qs_ot_env%matrix_buf2, qs_ot_env%matrix_buf2_im, &
3602 1876 : qs_ot_env%matrix_buf4)
3603 : CALL dbcsr_hadamard_product(qs_ot_env%matrix_buf2, qs_ot_env%matrix_cosp_b, &
3604 1876 : qs_ot_env%matrix_buf4)
3605 : CALL dbcsr_hadamard_product(qs_ot_env%matrix_buf2_im, qs_ot_env%matrix_cosp_b, &
3606 1876 : qs_ot_env%matrix_buf4_im)
3607 : CALL dbcsr_add(qs_ot_env%matrix_buf3, qs_ot_env%matrix_buf4, &
3608 1876 : alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
3609 : CALL dbcsr_add(qs_ot_env%matrix_buf3_im, qs_ot_env%matrix_buf4_im, &
3610 1876 : alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
3611 :
3612 : ! Transform back and add the Hermitian adjoint generated by dP.
3613 : CALL qs_ot_complex_multiply('N', 'N', qs_ot_env%matrix_r, qs_ot_env%matrix_r_im, &
3614 : qs_ot_env%matrix_buf3, qs_ot_env%matrix_buf3_im, &
3615 : qs_ot_env%matrix_buf1, qs_ot_env%matrix_buf1_im, &
3616 1876 : qs_ot_env%matrix_buf2)
3617 : CALL qs_ot_complex_multiply('N', 'C', qs_ot_env%matrix_buf1, qs_ot_env%matrix_buf1_im, &
3618 : qs_ot_env%matrix_r, qs_ot_env%matrix_r_im, &
3619 : qs_ot_env%matrix_buf3, qs_ot_env%matrix_buf3_im, &
3620 1876 : qs_ot_env%matrix_buf2)
3621 1876 : CALL dbcsr_get_info(qs_ot_env%matrix_buf3, distribution=dist)
3622 : CALL dbcsr_transposed(qs_ot_env%matrix_buf4, qs_ot_env%matrix_buf3, &
3623 : shallow_data_copy=.FALSE., use_distribution=dist, &
3624 1876 : transpose_distribution=.FALSE.)
3625 : CALL dbcsr_transposed(qs_ot_env%matrix_buf4_im, qs_ot_env%matrix_buf3_im, &
3626 : shallow_data_copy=.FALSE., use_distribution=dist, &
3627 1876 : transpose_distribution=.FALSE.)
3628 : CALL dbcsr_add(qs_ot_env%matrix_buf3, qs_ot_env%matrix_buf4, &
3629 1876 : alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
3630 : CALL dbcsr_add(qs_ot_env%matrix_buf3_im, qs_ot_env%matrix_buf4_im, &
3631 1876 : alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
3632 :
3633 : CALL qs_ot_complex_multiply('N', 'N', qs_ot_env%matrix_sx, qs_ot_env%matrix_sx_im, &
3634 : qs_ot_env%matrix_buf3, qs_ot_env%matrix_buf3_im, &
3635 : qs_ot_env%matrix_buf_nk, qs_ot_env%matrix_buf_nk_im, &
3636 1876 : qs_ot_env%matrix_tmp_nk)
3637 : CALL dbcsr_add(qs_ot_env%matrix_gx, qs_ot_env%matrix_buf_nk, &
3638 1876 : alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
3639 : CALL dbcsr_add(qs_ot_env%matrix_gx_im, qs_ot_env%matrix_buf_nk_im, &
3640 1876 : alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
3641 :
3642 : ! Preconditioner-aware projection onto the complex STRICT tangent space.
3643 1876 : CALL qs_ot_prepare_complex_tangent_metric(qs_ot_env)
3644 1876 : IF (ASSOCIATED(qs_ot_env%preconditioner)) THEN
3645 1846 : target_re => qs_ot_env%matrix_psc0
3646 1846 : target_im => qs_ot_env%matrix_psc0_im
3647 : ELSE
3648 30 : target_re => qs_ot_env%matrix_sc0
3649 30 : target_im => qs_ot_env%matrix_sc0_im
3650 : END IF
3651 :
3652 : CALL qs_ot_complex_multiply('C', 'N', target_re, target_im, &
3653 : qs_ot_env%matrix_gx, qs_ot_env%matrix_gx_im, &
3654 : qs_ot_env%matrix_buf1_ortho, qs_ot_env%matrix_buf1_ortho_im, &
3655 1876 : qs_ot_env%matrix_tmp_ortho)
3656 : CALL qs_ot_complex_multiply('N', 'N', qs_ot_env%matrix_os, qs_ot_env%matrix_os_im, &
3657 : qs_ot_env%matrix_buf1_ortho, qs_ot_env%matrix_buf1_ortho_im, &
3658 : qs_ot_env%matrix_buf2_ortho, qs_ot_env%matrix_buf2_ortho_im, &
3659 1876 : qs_ot_env%matrix_tmp_ortho)
3660 : CALL qs_ot_complex_multiply('N', 'N', qs_ot_env%matrix_sc0, qs_ot_env%matrix_sc0_im, &
3661 : qs_ot_env%matrix_buf2_ortho, qs_ot_env%matrix_buf2_ortho_im, &
3662 : qs_ot_env%matrix_buf_nk, qs_ot_env%matrix_buf_nk_im, &
3663 1876 : qs_ot_env%matrix_tmp_nk)
3664 : CALL dbcsr_add(qs_ot_env%matrix_gx, qs_ot_env%matrix_buf_nk, &
3665 1876 : alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
3666 : CALL dbcsr_add(qs_ot_env%matrix_gx_im, qs_ot_env%matrix_buf_nk_im, &
3667 1876 : alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
3668 :
3669 1876 : IF (qs_ot_env%settings%do_rotation) THEN
3670 1216 : CALL dbcsr_release(hc_rot_re)
3671 1216 : CALL dbcsr_release(hc_rot_im)
3672 1216 : CALL dbcsr_release(tmp_nk)
3673 : END IF
3674 :
3675 1876 : CALL timestop(handle)
3676 :
3677 1876 : END SUBROUTINE qs_ot_get_derivative_complex
3678 :
3679 : ! **************************************************************************************************
3680 : !> \brief ...
3681 : !> \param matrix_hc ...
3682 : !> \param matrix_x ...
3683 : !> \param matrix_sx ...
3684 : !> \param matrix_gx ...
3685 : !> \param qs_ot_env ...
3686 : ! **************************************************************************************************
3687 105480 : SUBROUTINE qs_ot_get_derivative_diag(matrix_hc, matrix_x, matrix_sx, matrix_gx, &
3688 : qs_ot_env)
3689 :
3690 : TYPE(dbcsr_type), POINTER :: matrix_hc, matrix_x, matrix_sx, matrix_gx
3691 : TYPE(qs_ot_type) :: qs_ot_env
3692 :
3693 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_ot_get_derivative_diag'
3694 : REAL(KIND=dp), PARAMETER :: rone = 1.0_dp, rzero = 0.0_dp
3695 :
3696 : INTEGER :: handle, k, n
3697 : TYPE(dbcsr_distribution_type) :: dist
3698 :
3699 35160 : CALL timeset(routineN, handle)
3700 :
3701 35160 : CALL dbcsr_get_info(matrix_x, nfullrows_total=n, nfullcols_total=k)
3702 :
3703 : ! go for the derivative now
3704 : ! this de/dc*(dX/dx)*sinp
3705 35160 : CALL dbcsr_multiply('N', 'N', rone, matrix_hc, qs_ot_env%matrix_sinp, rzero, matrix_gx)
3706 : ! overlap hc*x
3707 35160 : CALL dbcsr_multiply('T', 'N', rone, matrix_hc, matrix_x, rzero, qs_ot_env%matrix_buf2)
3708 : ! get it in the basis of the eigenvectors
3709 : CALL dbcsr_multiply('N', 'N', rone, qs_ot_env%matrix_buf2, qs_ot_env%matrix_r, &
3710 35160 : rzero, qs_ot_env%matrix_buf1)
3711 : CALL dbcsr_multiply('T', 'N', rone, qs_ot_env%matrix_r, qs_ot_env%matrix_buf1, &
3712 35160 : rzero, qs_ot_env%matrix_buf2)
3713 :
3714 : ! get the schur product of O_uv*B_uv
3715 : CALL dbcsr_hadamard_product(qs_ot_env%matrix_buf2, qs_ot_env%matrix_sinp_b, &
3716 35160 : qs_ot_env%matrix_buf3)
3717 :
3718 : ! overlap hc*c0
3719 : CALL dbcsr_multiply('T', 'N', rone, matrix_hc, qs_ot_env%matrix_c0, rzero, &
3720 35160 : qs_ot_env%matrix_buf2)
3721 : ! get it in the basis of the eigenvectors
3722 : CALL dbcsr_multiply('N', 'N', rone, qs_ot_env%matrix_buf2, qs_ot_env%matrix_r, &
3723 35160 : rzero, qs_ot_env%matrix_buf1)
3724 : CALL dbcsr_multiply('T', 'N', rone, qs_ot_env%matrix_r, qs_ot_env%matrix_buf1, &
3725 35160 : rzero, qs_ot_env%matrix_buf2)
3726 :
3727 : CALL dbcsr_hadamard_product(qs_ot_env%matrix_buf2, qs_ot_env%matrix_cosp_b, &
3728 35160 : qs_ot_env%matrix_buf4)
3729 :
3730 : ! add the two bs and compute b+b^T
3731 : CALL dbcsr_add(qs_ot_env%matrix_buf3, qs_ot_env%matrix_buf4, &
3732 35160 : alpha_scalar=rone, beta_scalar=rone)
3733 :
3734 : ! get the b in the eigenvector basis
3735 : CALL dbcsr_multiply('N', 'T', rone, qs_ot_env%matrix_buf3, qs_ot_env%matrix_r, &
3736 35160 : rzero, qs_ot_env%matrix_buf1)
3737 : CALL dbcsr_multiply('N', 'N', rone, qs_ot_env%matrix_r, qs_ot_env%matrix_buf1, &
3738 35160 : rzero, qs_ot_env%matrix_buf3)
3739 35160 : CALL dbcsr_get_info(qs_ot_env%matrix_buf3, distribution=dist)
3740 : CALL dbcsr_transposed(qs_ot_env%matrix_buf1, qs_ot_env%matrix_buf3, &
3741 : shallow_data_copy=.FALSE., use_distribution=dist, &
3742 35160 : transpose_distribution=.FALSE.)
3743 : CALL dbcsr_add(qs_ot_env%matrix_buf3, qs_ot_env%matrix_buf1, &
3744 35160 : alpha_scalar=rone, beta_scalar=rone)
3745 :
3746 : ! and add to the derivative
3747 : CALL dbcsr_multiply('N', 'N', rone, matrix_sx, qs_ot_env%matrix_buf3, &
3748 35160 : rone, matrix_gx)
3749 35160 : CALL timestop(handle)
3750 :
3751 35160 : END SUBROUTINE qs_ot_get_derivative_diag
3752 :
3753 : ! **************************************************************************************************
3754 : !> \brief compute the derivative of the taylor expansion below
3755 : !> \param matrix_hc ...
3756 : !> \param matrix_x ...
3757 : !> \param matrix_sx ...
3758 : !> \param matrix_gx ...
3759 : !> \param qs_ot_env ...
3760 : ! **************************************************************************************************
3761 151924 : SUBROUTINE qs_ot_get_derivative_taylor(matrix_hc, matrix_x, matrix_sx, matrix_gx, &
3762 : qs_ot_env)
3763 :
3764 : TYPE(dbcsr_type), POINTER :: matrix_hc, matrix_x, matrix_sx, matrix_gx
3765 : TYPE(qs_ot_type) :: qs_ot_env
3766 :
3767 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_ot_get_derivative_taylor'
3768 : REAL(KIND=dp), PARAMETER :: rone = 1.0_dp, rzero = 0.0_dp
3769 :
3770 : INTEGER :: handle, i, k, n
3771 : REAL(KIND=dp) :: cosfactor, sinfactor
3772 : TYPE(dbcsr_distribution_type) :: dist
3773 : TYPE(dbcsr_type), POINTER :: matrix_left, matrix_right
3774 :
3775 42613 : CALL timeset(routineN, handle)
3776 :
3777 42613 : CALL dbcsr_get_info(matrix_x, nfullrows_total=n, nfullcols_total=k)
3778 :
3779 : ! go for the derivative now
3780 : ! this de/dc*(dX/dx)*sinp i.e. zeroth order
3781 42613 : CALL dbcsr_multiply('N', 'N', rone, matrix_hc, qs_ot_env%matrix_sinp, rzero, matrix_gx)
3782 :
3783 42613 : IF (qs_ot_env%taylor_order <= 0) THEN
3784 9264 : CALL timestop(handle)
3785 9264 : RETURN
3786 : END IF
3787 :
3788 : ! we store the matrix that will multiply sx in matrix_r
3789 33349 : CALL dbcsr_set(qs_ot_env%matrix_r, rzero)
3790 :
3791 : ! just better names for matrix_cosp_b and matrix_sinp_b (they are buffer space here)
3792 33349 : matrix_left => qs_ot_env%matrix_cosp_b
3793 33349 : matrix_right => qs_ot_env%matrix_sinp_b
3794 :
3795 : ! overlap hc*x and add its transpose to matrix_left
3796 33349 : CALL dbcsr_multiply('T', 'N', rone, matrix_hc, matrix_x, rzero, matrix_left)
3797 33349 : CALL dbcsr_get_info(matrix_left, distribution=dist)
3798 : CALL dbcsr_transposed(qs_ot_env%matrix_buf1, matrix_left, &
3799 : shallow_data_copy=.FALSE., use_distribution=dist, &
3800 33349 : transpose_distribution=.FALSE.)
3801 : CALL dbcsr_add(matrix_left, qs_ot_env%matrix_buf1, &
3802 33349 : alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
3803 33349 : CALL dbcsr_copy(matrix_right, matrix_left)
3804 :
3805 : ! first order
3806 33349 : sinfactor = -1.0_dp/(2.0_dp*3.0_dp)
3807 : CALL dbcsr_add(qs_ot_env%matrix_r, matrix_left, &
3808 33349 : alpha_scalar=1.0_dp, beta_scalar=sinfactor)
3809 :
3810 : ! M
3811 : ! OM+MO
3812 : ! OOM+OMO+MOO
3813 : ! ...
3814 71037 : DO i = 2, qs_ot_env%taylor_order
3815 37688 : sinfactor = sinfactor*(-1.0_dp)/REAL(2*i*(2*i + 1), KIND=dp)
3816 37688 : CALL dbcsr_multiply('N', 'N', rone, qs_ot_env%matrix_p, matrix_left, rzero, qs_ot_env%matrix_buf1)
3817 37688 : CALL dbcsr_multiply('N', 'N', rone, matrix_right, qs_ot_env%matrix_p, rzero, matrix_left)
3818 37688 : CALL dbcsr_copy(matrix_right, matrix_left)
3819 : CALL dbcsr_add(matrix_left, qs_ot_env%matrix_buf1, &
3820 37688 : 1.0_dp, 1.0_dp)
3821 : CALL dbcsr_add(qs_ot_env%matrix_r, matrix_left, &
3822 71037 : alpha_scalar=1.0_dp, beta_scalar=sinfactor)
3823 : END DO
3824 :
3825 : ! overlap hc*c0 and add its transpose to matrix_left
3826 33349 : CALL dbcsr_multiply('T', 'N', rone, matrix_hc, qs_ot_env%matrix_c0, rzero, matrix_left)
3827 33349 : CALL dbcsr_get_info(matrix_left, distribution=dist)
3828 : CALL dbcsr_transposed(qs_ot_env%matrix_buf1, matrix_left, &
3829 : shallow_data_copy=.FALSE., use_distribution=dist, &
3830 33349 : transpose_distribution=.FALSE.)
3831 33349 : CALL dbcsr_add(matrix_left, qs_ot_env%matrix_buf1, 1.0_dp, 1.0_dp)
3832 33349 : CALL dbcsr_copy(matrix_right, matrix_left)
3833 :
3834 : ! first order
3835 33349 : cosfactor = -1.0_dp/(1.0_dp*2.0_dp)
3836 : CALL dbcsr_add(qs_ot_env%matrix_r, matrix_left, &
3837 33349 : alpha_scalar=1.0_dp, beta_scalar=cosfactor)
3838 :
3839 : ! M
3840 : ! OM+MO
3841 : ! OOM+OMO+MOO
3842 : ! ...
3843 71037 : DO i = 2, qs_ot_env%taylor_order
3844 37688 : cosfactor = cosfactor*(-1.0_dp)/REAL(2*i*(2*i - 1), KIND=dp)
3845 37688 : CALL dbcsr_multiply('N', 'N', rone, qs_ot_env%matrix_p, matrix_left, rzero, qs_ot_env%matrix_buf1)
3846 37688 : CALL dbcsr_multiply('N', 'N', rone, matrix_right, qs_ot_env%matrix_p, rzero, matrix_left)
3847 37688 : CALL dbcsr_copy(matrix_right, matrix_left)
3848 37688 : CALL dbcsr_add(matrix_left, qs_ot_env%matrix_buf1, 1.0_dp, 1.0_dp)
3849 : CALL dbcsr_add(qs_ot_env%matrix_r, matrix_left, &
3850 71037 : alpha_scalar=1.0_dp, beta_scalar=cosfactor)
3851 : END DO
3852 :
3853 : ! and add to the derivative
3854 33349 : CALL dbcsr_multiply('N', 'N', rone, matrix_sx, qs_ot_env%matrix_r, rone, matrix_gx)
3855 :
3856 33349 : CALL timestop(handle)
3857 :
3858 42613 : END SUBROUTINE qs_ot_get_derivative_taylor
3859 :
3860 : ! *************************************************************************************************
3861 : !> \brief computes a taylor expansion.
3862 : !> \param qs_ot_env ...
3863 : ! **************************************************************************************************
3864 93793 : SUBROUTINE qs_ot_p2m_taylor(qs_ot_env)
3865 : TYPE(qs_ot_type) :: qs_ot_env
3866 :
3867 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_ot_p2m_taylor'
3868 : REAL(KIND=dp), PARAMETER :: rone = 1.0_dp, rzero = 0.0_dp
3869 :
3870 : INTEGER :: handle, i, k
3871 : REAL(KIND=dp) :: cosfactor, sinfactor
3872 :
3873 58064 : CALL timeset(routineN, handle)
3874 :
3875 : ! zeroth order
3876 58064 : CALL dbcsr_set(qs_ot_env%matrix_cosp, rzero)
3877 58064 : CALL dbcsr_set(qs_ot_env%matrix_sinp, rzero)
3878 58064 : CALL dbcsr_add_on_diag(qs_ot_env%matrix_cosp, rone)
3879 58064 : CALL dbcsr_add_on_diag(qs_ot_env%matrix_sinp, rone)
3880 :
3881 58064 : IF (qs_ot_env%taylor_order <= 0) THEN
3882 10008 : CALL timestop(handle)
3883 22335 : RETURN
3884 : END IF
3885 :
3886 : ! first order
3887 48056 : cosfactor = -1.0_dp/(1.0_dp*2.0_dp)
3888 48056 : sinfactor = -1.0_dp/(2.0_dp*3.0_dp)
3889 48056 : CALL dbcsr_add(qs_ot_env%matrix_cosp, qs_ot_env%matrix_p, alpha_scalar=1.0_dp, beta_scalar=cosfactor)
3890 48056 : CALL dbcsr_add(qs_ot_env%matrix_sinp, qs_ot_env%matrix_p, alpha_scalar=1.0_dp, beta_scalar=sinfactor)
3891 48056 : IF (qs_ot_env%taylor_order <= 1) THEN
3892 12327 : CALL timestop(handle)
3893 12327 : RETURN
3894 : END IF
3895 :
3896 : ! other orders
3897 35729 : CALL dbcsr_get_info(qs_ot_env%matrix_p, nfullrows_total=k)
3898 35729 : CALL dbcsr_copy(qs_ot_env%matrix_r, qs_ot_env%matrix_p)
3899 :
3900 90080 : DO i = 2, qs_ot_env%taylor_order
3901 : ! new power of p
3902 : CALL dbcsr_multiply('N', 'N', rone, qs_ot_env%matrix_p, qs_ot_env%matrix_r, &
3903 54351 : rzero, qs_ot_env%matrix_buf1)
3904 54351 : CALL dbcsr_copy(qs_ot_env%matrix_r, qs_ot_env%matrix_buf1)
3905 : ! add to the taylor expansion so far
3906 54351 : cosfactor = cosfactor*(-1.0_dp)/REAL(2*i*(2*i - 1), KIND=dp)
3907 54351 : sinfactor = sinfactor*(-1.0_dp)/REAL(2*i*(2*i + 1), KIND=dp)
3908 : CALL dbcsr_add(qs_ot_env%matrix_cosp, qs_ot_env%matrix_r, &
3909 54351 : alpha_scalar=1.0_dp, beta_scalar=cosfactor)
3910 : CALL dbcsr_add(qs_ot_env%matrix_sinp, qs_ot_env%matrix_r, &
3911 90080 : alpha_scalar=1.0_dp, beta_scalar=sinfactor)
3912 : END DO
3913 :
3914 35729 : CALL timestop(handle)
3915 :
3916 : END SUBROUTINE qs_ot_p2m_taylor
3917 :
3918 : ! **************************************************************************************************
3919 : !> \brief given p, computes - eigenstuff (matrix_r,evals)
3920 : !> - cos(p^0.5),p^(-0.5)*sin(p^0.5)
3921 : !> - the real b matrices, needed for the derivatives of these guys
3922 : !> cosp_b_ij=(1/(2pii) * int(cos(z^1/2)/((z-eval(i))*(z-eval(j))))
3923 : !> sinp_b_ij=(1/(2pii) * int(z^(-1/2)*sin(z^1/2)/((z-eval(i))*(z-eval(j))))
3924 : !> \param qs_ot_env ...
3925 : ! **************************************************************************************************
3926 206788 : SUBROUTINE qs_ot_p2m_diag(qs_ot_env)
3927 :
3928 : TYPE(qs_ot_type) :: qs_ot_env
3929 :
3930 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_ot_p2m_diag'
3931 : REAL(KIND=dp), PARAMETER :: rone = 1.0_dp, rzero = 0.0_dp
3932 :
3933 : INTEGER :: col, col_offset, col_size, handle, i, j, &
3934 : k, row, row_offset, row_size
3935 51697 : REAL(dp), DIMENSION(:, :), POINTER :: block
3936 : REAL(KIND=dp) :: a, b
3937 : TYPE(dbcsr_iterator_type) :: iter
3938 :
3939 51697 : CALL timeset(routineN, handle)
3940 :
3941 51697 : CALL dbcsr_get_info(qs_ot_env%matrix_p, nfullrows_total=k)
3942 51697 : CALL dbcsr_copy(qs_ot_env%matrix_buf1, qs_ot_env%matrix_p)
3943 : CALL cp_dbcsr_syevd(qs_ot_env%matrix_buf1, qs_ot_env%matrix_r, qs_ot_env%evals, &
3944 51697 : qs_ot_env%para_env, qs_ot_env%blacs_env)
3945 536412 : DO i = 1, k
3946 536412 : qs_ot_env%evals(i) = MAX(0.0_dp, qs_ot_env%evals(i))
3947 : END DO
3948 :
3949 51697 : !$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(i) SHARED(k,qs_ot_env)
3950 : DO i = 1, k
3951 : qs_ot_env%dum(i) = COS(SQRT(qs_ot_env%evals(i)))
3952 : END DO
3953 51697 : CALL dbcsr_copy(qs_ot_env%matrix_buf1, qs_ot_env%matrix_r)
3954 51697 : CALL dbcsr_scale_by_vector(qs_ot_env%matrix_buf1, alpha=qs_ot_env%dum, side='right')
3955 : CALL dbcsr_multiply('N', 'T', rone, qs_ot_env%matrix_r, qs_ot_env%matrix_buf1, &
3956 51697 : rzero, qs_ot_env%matrix_cosp)
3957 :
3958 51697 : !$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(i) SHARED(k,qs_ot_env)
3959 : DO i = 1, k
3960 : qs_ot_env%dum(i) = qs_ot_sinc(SQRT(qs_ot_env%evals(i)))
3961 : END DO
3962 51697 : CALL dbcsr_copy(qs_ot_env%matrix_buf1, qs_ot_env%matrix_r)
3963 51697 : CALL dbcsr_scale_by_vector(qs_ot_env%matrix_buf1, alpha=qs_ot_env%dum, side='right')
3964 : CALL dbcsr_multiply('N', 'T', rone, qs_ot_env%matrix_r, qs_ot_env%matrix_buf1, &
3965 51697 : rzero, qs_ot_env%matrix_sinp)
3966 :
3967 51697 : CALL dbcsr_copy(qs_ot_env%matrix_cosp_b, qs_ot_env%matrix_cosp)
3968 51697 : CALL dbcsr_iterator_start(iter, qs_ot_env%matrix_cosp_b)
3969 88219 : DO WHILE (dbcsr_iterator_blocks_left(iter))
3970 : CALL dbcsr_iterator_next_block(iter, row, col, block, &
3971 : row_size=row_size, col_size=col_size, &
3972 36522 : row_offset=row_offset, col_offset=col_offset)
3973 552641 : DO j = 1, col_size
3974 10706919 : DO i = 1, row_size
3975 : a = (SQRT(qs_ot_env%evals(row_offset + i - 1)) &
3976 10205975 : - SQRT(qs_ot_env%evals(col_offset + j - 1)))/2.0_dp
3977 : b = (SQRT(qs_ot_env%evals(row_offset + i - 1)) &
3978 10205975 : + SQRT(qs_ot_env%evals(col_offset + j - 1)))/2.0_dp
3979 10670397 : block(i, j) = -0.5_dp*qs_ot_sinc(a)*qs_ot_sinc(b)
3980 : END DO
3981 : END DO
3982 : END DO
3983 51697 : CALL dbcsr_iterator_stop(iter)
3984 :
3985 51697 : CALL dbcsr_copy(qs_ot_env%matrix_sinp_b, qs_ot_env%matrix_sinp)
3986 51697 : CALL dbcsr_iterator_start(iter, qs_ot_env%matrix_sinp_b)
3987 88219 : DO WHILE (dbcsr_iterator_blocks_left(iter))
3988 : CALL dbcsr_iterator_next_block(iter, row, col, block, &
3989 : row_size=row_size, col_size=col_size, &
3990 36522 : row_offset=row_offset, col_offset=col_offset)
3991 552641 : DO j = 1, col_size
3992 10706919 : DO i = 1, row_size
3993 10205975 : a = SQRT(qs_ot_env%evals(row_offset + i - 1))
3994 10205975 : b = SQRT(qs_ot_env%evals(col_offset + j - 1))
3995 10670397 : block(i, j) = qs_ot_sincf(a, b)
3996 : END DO
3997 : END DO
3998 : END DO
3999 51697 : CALL dbcsr_iterator_stop(iter)
4000 :
4001 51697 : CALL timestop(handle)
4002 :
4003 51697 : END SUBROUTINE qs_ot_p2m_diag
4004 :
4005 : ! **************************************************************************************************
4006 : !> \brief diagonalize Hermitian P and build cos(sqrt(P)), sinc(sqrt(P)), and Frechet kernels
4007 : !> \param qs_ot_env complex STRICT channel state
4008 : ! **************************************************************************************************
4009 8032 : SUBROUTINE qs_ot_p2m_diag_complex(qs_ot_env)
4010 : TYPE(qs_ot_type) :: qs_ot_env
4011 :
4012 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_ot_p2m_diag_complex'
4013 :
4014 : INTEGER :: col, col_offset, col_size, handle, i, j, &
4015 : k, row, row_offset, row_size
4016 2008 : REAL(dp), DIMENSION(:, :), POINTER :: block
4017 : REAL(KIND=dp) :: a, b
4018 : TYPE(dbcsr_iterator_type) :: iter
4019 :
4020 2008 : CALL timeset(routineN, handle)
4021 :
4022 2008 : CALL dbcsr_get_info(qs_ot_env%matrix_p, nfullrows_total=k)
4023 : CALL cp_dbcsr_heevd(matrix_re=qs_ot_env%matrix_p, matrix_im=qs_ot_env%matrix_p_im, &
4024 : eigenvectors_re=qs_ot_env%matrix_r, eigenvectors_im=qs_ot_env%matrix_r_im, &
4025 : eigenvalues=qs_ot_env%evals, para_env=qs_ot_env%para_env, &
4026 2008 : blacs_env=qs_ot_env%blacs_env)
4027 15554 : DO i = 1, k
4028 15554 : qs_ot_env%evals(i) = MAX(0.0_dp, qs_ot_env%evals(i))
4029 : END DO
4030 :
4031 15554 : DO i = 1, k
4032 15554 : qs_ot_env%dum(i) = COS(SQRT(qs_ot_env%evals(i)))
4033 : END DO
4034 2008 : CALL dbcsr_copy(qs_ot_env%matrix_buf1, qs_ot_env%matrix_r)
4035 2008 : CALL dbcsr_copy(qs_ot_env%matrix_buf1_im, qs_ot_env%matrix_r_im)
4036 2008 : CALL dbcsr_scale_by_vector(qs_ot_env%matrix_buf1, alpha=qs_ot_env%dum, side='right')
4037 2008 : CALL dbcsr_scale_by_vector(qs_ot_env%matrix_buf1_im, alpha=qs_ot_env%dum, side='right')
4038 : CALL qs_ot_complex_multiply('N', 'C', qs_ot_env%matrix_buf1, qs_ot_env%matrix_buf1_im, &
4039 : qs_ot_env%matrix_r, qs_ot_env%matrix_r_im, &
4040 : qs_ot_env%matrix_cosp, qs_ot_env%matrix_cosp_im, &
4041 2008 : qs_ot_env%matrix_buf2)
4042 :
4043 15554 : DO i = 1, k
4044 15554 : qs_ot_env%dum(i) = qs_ot_sinc(SQRT(qs_ot_env%evals(i)))
4045 : END DO
4046 2008 : CALL dbcsr_copy(qs_ot_env%matrix_buf1, qs_ot_env%matrix_r)
4047 2008 : CALL dbcsr_copy(qs_ot_env%matrix_buf1_im, qs_ot_env%matrix_r_im)
4048 2008 : CALL dbcsr_scale_by_vector(qs_ot_env%matrix_buf1, alpha=qs_ot_env%dum, side='right')
4049 2008 : CALL dbcsr_scale_by_vector(qs_ot_env%matrix_buf1_im, alpha=qs_ot_env%dum, side='right')
4050 : CALL qs_ot_complex_multiply('N', 'C', qs_ot_env%matrix_buf1, qs_ot_env%matrix_buf1_im, &
4051 : qs_ot_env%matrix_r, qs_ot_env%matrix_r_im, &
4052 : qs_ot_env%matrix_sinp, qs_ot_env%matrix_sinp_im, &
4053 2008 : qs_ot_env%matrix_buf2)
4054 :
4055 2008 : CALL dbcsr_copy(qs_ot_env%matrix_cosp_b, qs_ot_env%matrix_cosp)
4056 2008 : CALL dbcsr_iterator_start(iter, qs_ot_env%matrix_cosp_b)
4057 3091 : DO WHILE (dbcsr_iterator_blocks_left(iter))
4058 : CALL dbcsr_iterator_next_block(iter, row, col, block, &
4059 : row_size=row_size, col_size=col_size, &
4060 1083 : row_offset=row_offset, col_offset=col_offset)
4061 11582 : DO j = 1, col_size
4062 124853 : DO i = 1, row_size
4063 : a = (SQRT(qs_ot_env%evals(row_offset + i - 1)) &
4064 115279 : - SQRT(qs_ot_env%evals(col_offset + j - 1)))/2.0_dp
4065 : b = (SQRT(qs_ot_env%evals(row_offset + i - 1)) &
4066 115279 : + SQRT(qs_ot_env%evals(col_offset + j - 1)))/2.0_dp
4067 123770 : block(i, j) = -0.5_dp*qs_ot_sinc(a)*qs_ot_sinc(b)
4068 : END DO
4069 : END DO
4070 : END DO
4071 2008 : CALL dbcsr_iterator_stop(iter)
4072 :
4073 2008 : CALL dbcsr_copy(qs_ot_env%matrix_sinp_b, qs_ot_env%matrix_sinp)
4074 2008 : CALL dbcsr_iterator_start(iter, qs_ot_env%matrix_sinp_b)
4075 3091 : DO WHILE (dbcsr_iterator_blocks_left(iter))
4076 : CALL dbcsr_iterator_next_block(iter, row, col, block, &
4077 : row_size=row_size, col_size=col_size, &
4078 1083 : row_offset=row_offset, col_offset=col_offset)
4079 11582 : DO j = 1, col_size
4080 124853 : DO i = 1, row_size
4081 115279 : a = SQRT(qs_ot_env%evals(row_offset + i - 1))
4082 115279 : b = SQRT(qs_ot_env%evals(col_offset + j - 1))
4083 123770 : block(i, j) = qs_ot_sincf(a, b)
4084 : END DO
4085 : END DO
4086 : END DO
4087 2008 : CALL dbcsr_iterator_stop(iter)
4088 :
4089 2008 : CALL timestop(handle)
4090 :
4091 2008 : END SUBROUTINE qs_ot_p2m_diag_complex
4092 :
4093 : ! **************************************************************************************************
4094 : !> \brief computes sin(x)/x for all values of the argument
4095 : !> \param x ...
4096 : !> \return ...
4097 : ! **************************************************************************************************
4098 28880599 : FUNCTION qs_ot_sinc(x)
4099 :
4100 : REAL(KIND=dp), INTENT(IN) :: x
4101 : REAL(KIND=dp) :: qs_ot_sinc
4102 :
4103 : REAL(KIND=dp), PARAMETER :: q1 = 1.0_dp, q2 = -q1/(2.0_dp*3.0_dp), q3 = -q2/(4.0_dp*5.0_dp), &
4104 : q4 = -q3/(6.0_dp*7.0_dp), q5 = -q4/(8.0_dp*9.0_dp), q6 = -q5/(10.0_dp*11.0_dp), &
4105 : q7 = -q6/(12.0_dp*13.0_dp), q8 = -q7/(14.0_dp*15.0_dp), q9 = -q8/(16.0_dp*17.0_dp), &
4106 : q10 = -q9/(18.0_dp*19.0_dp)
4107 :
4108 : REAL(KIND=dp) :: y
4109 :
4110 28880599 : IF (ABS(x) > 0.5_dp) THEN
4111 8707781 : qs_ot_sinc = SIN(x)/x
4112 : ELSE
4113 20172818 : y = x*x
4114 20172818 : qs_ot_sinc = q1 + y*(q2 + y*(q3 + y*(q4 + y*(q5 + y*(q6 + y*(q7 + y*(q8 + y*(q9 + y*(q10)))))))))
4115 : END IF
4116 28880599 : END FUNCTION qs_ot_sinc
4117 :
4118 : ! **************************************************************************************************
4119 : !> \brief computes (1/(x^2-y^2))*(sinc(x)-sinc(y)) for all positive values of the arguments
4120 : !> \param xa ...
4121 : !> \param ya ...
4122 : !> \return ...
4123 : ! **************************************************************************************************
4124 10321254 : FUNCTION qs_ot_sincf(xa, ya)
4125 :
4126 : REAL(KIND=dp), INTENT(IN) :: xa, ya
4127 : REAL(KIND=dp) :: qs_ot_sincf
4128 :
4129 : INTEGER :: i
4130 : REAL(KIND=dp) :: a, b, rs, sf, x, xs, y, ybx, ybxs
4131 :
4132 : ! this is currently a limit of the routine, could be removed rather easily
4133 10321254 : IF (xa < 0) CPABORT("x is negative")
4134 10321254 : IF (ya < 0) CPABORT("y is negative")
4135 :
4136 10321254 : IF (xa < ya) THEN
4137 4935528 : x = ya
4138 4935528 : y = xa
4139 : ELSE
4140 5385726 : x = xa
4141 5385726 : y = ya
4142 : END IF
4143 :
4144 10321254 : IF (x < 0.5_dp) THEN ! use series, keeping in mind that x,y,x+y,x-y can all be zero
4145 :
4146 6451339 : qs_ot_sincf = 0.0_dp
4147 6451339 : IF (x > 0.0_dp) THEN
4148 6234744 : ybx = y/x
4149 : ELSE ! should be irrelevant !?
4150 : ybx = 0.0_dp
4151 : END IF
4152 :
4153 6451339 : sf = -1.0_dp/((1.0_dp + ybx)*6.0_dp)
4154 6451339 : rs = 1.0_dp
4155 6451339 : ybxs = ybx
4156 6451339 : xs = 1.0_dp
4157 :
4158 70964729 : DO i = 1, 10
4159 64513390 : qs_ot_sincf = qs_ot_sincf + sf*rs*xs*(1.0_dp + ybxs)
4160 64513390 : sf = -sf/(REAL((2*i + 2), dp)*REAL((2*i + 3), dp))
4161 64513390 : rs = rs + ybxs
4162 64513390 : ybxs = ybxs*ybx
4163 70964729 : xs = xs*x*x
4164 : END DO
4165 :
4166 : ELSE ! no series expansion
4167 3869915 : IF (x - y > 0.1_dp) THEN ! safe to use the normal form
4168 3582932 : qs_ot_sincf = (qs_ot_sinc(x) - qs_ot_sinc(y))/((x + y)*(x - y))
4169 : ELSE
4170 286983 : a = (x + y)/2.0_dp
4171 286983 : b = (x - y)/2.0_dp ! might be close to zero
4172 : ! y (=(a-b)) can not be close to zero since it is close to x>0.5
4173 286983 : qs_ot_sincf = (qs_ot_sinc(b)*COS(a) - qs_ot_sinc(a)*COS(b))/(2*x*y)
4174 : END IF
4175 : END IF
4176 :
4177 10321254 : END FUNCTION qs_ot_sincf
4178 :
4179 68062 : END MODULE qs_ot
|
