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