(git:50ddb19)
Loading...
Searching...
No Matches
qs_ot.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \brief orbital transformations
10!> \par History
11!> Added Taylor expansion based computation of the matrix functions (01.2004)
12!> added additional rotation variables for non-equivalent occupied orbs (08.2004)
13!> \author Joost VandeVondele (06.2002)
14! **************************************************************************************************
15MODULE qs_ot
16 USE arnoldi_api, ONLY: arnoldi_extremal
17 USE cp_dbcsr_api, ONLY: &
18 dbcsr_add, dbcsr_copy, dbcsr_distribution_type, dbcsr_filter, dbcsr_get_block_p, &
19 dbcsr_get_info, dbcsr_get_occupation, dbcsr_init_p, dbcsr_iterator_blocks_left, &
20 dbcsr_iterator_next_block, dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, &
21 dbcsr_multiply, dbcsr_release, dbcsr_release_p, dbcsr_scale, dbcsr_set, dbcsr_transposed, &
22 dbcsr_type
23 USE cp_dbcsr_cholesky, ONLY: cp_dbcsr_cholesky_decompose,&
24 cp_dbcsr_cholesky_invert,&
25 cp_dbcsr_cholesky_restore
26 USE cp_dbcsr_contrib, ONLY: dbcsr_add_on_diag,&
27 dbcsr_frobenius_norm,&
28 dbcsr_gershgorin_norm,&
29 dbcsr_hadamard_product,&
30 dbcsr_scale_by_vector
31 USE cp_dbcsr_diag, ONLY: cp_dbcsr_heevd,&
32 cp_dbcsr_syevd
33 USE kinds, ONLY: dp
34 USE message_passing, ONLY: mp_comm_type
35 USE preconditioner, ONLY: apply_preconditioner
36 USE preconditioner_types, ONLY: preconditioner_type
37 USE qs_ot_types, ONLY: qs_ot_type
38#include "./base/base_uses.f90"
39
40 IMPLICIT NONE
41 PRIVATE
42
43 PUBLIC :: qs_ot_get_p
44 PUBLIC :: qs_ot_get_orbitals
45 PUBLIC :: qs_ot_get_derivative
46 PUBLIC :: qs_ot_get_orbitals_ref
47 PUBLIC :: qs_ot_get_derivative_ref
48 PUBLIC :: qs_ot_new_preconditioner
49 PRIVATE :: qs_ot_p2m_diag
50 PRIVATE :: qs_ot_sinc
51 PRIVATE :: qs_ot_ref_poly
52 PRIVATE :: qs_ot_ref_chol
53 PRIVATE :: qs_ot_ref_lwdn
54 PRIVATE :: qs_ot_ref_decide
55 PRIVATE :: qs_ot_ref_update
56 PRIVATE :: qs_ot_refine
57 PRIVATE :: qs_ot_on_the_fly_localize
58
59 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_ot'
60
61CONTAINS
62
63! **************************************************************************************************
64!> \brief gets ready to use the preconditioner/ or renew the preconditioner
65!> only keeps a pointer to the preconditioner.
66!> If you change the preconditioner, you have to call this routine
67!> you remain responsible of proper deallocate of your preconditioner
68!> (or you can reuse it on the next step of the computation)
69!> \param qs_ot_env ...
70!> \param preconditioner ...
71! **************************************************************************************************
72 SUBROUTINE qs_ot_new_preconditioner(qs_ot_env, preconditioner)
73 TYPE(qs_ot_type) :: qs_ot_env
74 TYPE(preconditioner_type), POINTER :: preconditioner
75
76 INTEGER :: ncoef
77
78 qs_ot_env%preconditioner => preconditioner
79 qs_ot_env%os_valid = .false.
80 IF (.NOT. ASSOCIATED(qs_ot_env%matrix_psc0)) THEN
81 CALL dbcsr_init_p(qs_ot_env%matrix_psc0)
82 CALL dbcsr_copy(qs_ot_env%matrix_psc0, qs_ot_env%matrix_sc0, 'matrix_psc0')
83 END IF
84
85 IF (.NOT. qs_ot_env%use_dx) THEN
86 qs_ot_env%use_dx = .true.
87 CALL dbcsr_init_p(qs_ot_env%matrix_dx)
88 CALL dbcsr_copy(qs_ot_env%matrix_dx, qs_ot_env%matrix_gx, 'matrix_dx')
89 IF (qs_ot_env%settings%do_rotation) THEN
90 CALL dbcsr_init_p(qs_ot_env%rot_mat_dx)
91 CALL dbcsr_copy(qs_ot_env%rot_mat_dx, qs_ot_env%rot_mat_gx, 'rot_mat_dx')
92 END IF
93 IF (qs_ot_env%settings%do_ener) THEN
94 ncoef = SIZE(qs_ot_env%ener_gx)
95 ALLOCATE (qs_ot_env%ener_dx(ncoef))
96 qs_ot_env%ener_dx = 0.0_dp
97 END IF
98 END IF
99
100 END SUBROUTINE qs_ot_new_preconditioner
101
102! **************************************************************************************************
103!> \brief ...
104!> \param qs_ot_env ...
105!> \param C_NEW ...
106!> \param SC ...
107!> \param G_OLD ...
108!> \param D ...
109! **************************************************************************************************
110 SUBROUTINE qs_ot_on_the_fly_localize(qs_ot_env, C_NEW, SC, G_OLD, D)
111 !
112 TYPE(qs_ot_type) :: qs_ot_env
113 TYPE(dbcsr_type), POINTER :: c_new, sc, g_old, d
114
115 CHARACTER(len=*), PARAMETER :: routinen = 'qs_ot_on_the_fly_localize'
116 INTEGER, PARAMETER :: taylor_order = 50
117 REAL(kind=dp), PARAMETER :: alpha = 0.1_dp, f2_eps = 0.01_dp
118
119 INTEGER :: col, col_size, handle, i, k, n, p, row, &
120 row_size
121 REAL(dp), DIMENSION(:, :), POINTER :: block
122 REAL(kind=dp) :: expfactor, f2, norm_fro, norm_gct, tmp
123 TYPE(dbcsr_distribution_type) :: dist
124 TYPE(dbcsr_iterator_type) :: iter
125 TYPE(dbcsr_type), POINTER :: c, gp1, gp2, gu, u
126 TYPE(mp_comm_type) :: group
127
128 CALL timeset(routinen, handle)
129 !
130 !
131 CALL dbcsr_get_info(c_new, nfullrows_total=n, nfullcols_total=k)
132 !
133 ! C = C*expm(-G)
134 gu => qs_ot_env%buf1_k_k_nosym ! a buffer
135 u => qs_ot_env%buf2_k_k_nosym ! a buffer
136 gp1 => qs_ot_env%buf3_k_k_nosym ! a buffer
137 gp2 => qs_ot_env%buf4_k_k_nosym ! a buffer
138 c => qs_ot_env%buf1_n_k ! a buffer
139 !
140 ! compute the derivative of the norm
141 !-------------------------------------------------------------------
142 ! (x^2+eps)^1/2
143 f2 = 0.0_dp
144 CALL dbcsr_copy(c, c_new)
145 CALL dbcsr_iterator_start(iter, c)
146 DO WHILE (dbcsr_iterator_blocks_left(iter))
147 CALL dbcsr_iterator_next_block(iter, row, col, block, row_size=row_size, col_size=col_size)
148 DO p = 1, col_size ! p
149 DO i = 1, row_size ! i
150 tmp = sqrt(block(i, p)**2 + f2_eps)
151 f2 = f2 + tmp
152 block(i, p) = block(i, p)/tmp
153 END DO
154 END DO
155 END DO
156 CALL dbcsr_iterator_stop(iter)
157 CALL dbcsr_get_info(c, group=group)
158 CALL group%sum(f2)
159 !
160 !
161 CALL dbcsr_multiply('T', 'N', 1.0_dp, c, c_new, 0.0_dp, gu)
162 !
163 ! antisymetrize
164 CALL dbcsr_get_info(gu, distribution=dist)
165 CALL dbcsr_transposed(u, gu, shallow_data_copy=.false., &
166 use_distribution=dist, &
167 transpose_distribution=.false.)
168 CALL dbcsr_add(gu, u, alpha_scalar=-0.5_dp, beta_scalar=0.5_dp)
169 !-------------------------------------------------------------------
170 !
171 norm_fro = dbcsr_frobenius_norm(gu)
172 norm_gct = dbcsr_gershgorin_norm(gu)
173 !write(*,*) 'qs_ot_localize: ||P-I||_f=',norm_fro,' ||P-I||_GCT=',norm_gct
174 !
175 !kscale = CEILING(LOG(MIN(norm_fro,norm_gct))/LOG(2.0_dp))
176 !scale = LOG(MIN(norm_fro,norm_gct))/LOG(2.0_dp)
177 !write(*,*) 'qs_ot_localize: scale=',scale,' kscale=',kscale
178 !
179 ! rescale for steepest descent
180 CALL dbcsr_scale(gu, -alpha)
181 !
182 ! compute unitary transform
183 ! zeroth and first order
184 expfactor = 1.0_dp
185 CALL dbcsr_copy(u, gu)
186 CALL dbcsr_scale(u, expfactor)
187 CALL dbcsr_add_on_diag(u, 1.0_dp)
188 ! other orders
189 CALL dbcsr_copy(gp1, gu)
190 DO i = 2, taylor_order
191 ! new power of G
192 CALL dbcsr_multiply('N', 'N', 1.0_dp, gu, gp1, 0.0_dp, gp2)
193 CALL dbcsr_copy(gp1, gp2)
194 ! add to the taylor expansion so far
195 expfactor = expfactor/real(i, kind=dp)
196 CALL dbcsr_add(u, gp1, alpha_scalar=1.0_dp, beta_scalar=expfactor)
197 norm_fro = dbcsr_frobenius_norm(gp1)
198 !write(*,*) 'Taylor expansion i=',i,' norm(X^i)/i!=',norm_fro*expfactor
199 IF (norm_fro*expfactor < 1.0e-10_dp) EXIT
200 END DO
201 !
202 ! rotate MOs
203 CALL dbcsr_multiply('N', 'N', 1.0_dp, c_new, u, 0.0_dp, c)
204 CALL dbcsr_copy(c_new, c)
205 !
206 ! rotate SC
207 CALL dbcsr_multiply('N', 'N', 1.0_dp, sc, u, 0.0_dp, c)
208 CALL dbcsr_copy(sc, c)
209 !
210 ! rotate D_i
211 CALL dbcsr_multiply('N', 'N', 1.0_dp, d, u, 0.0_dp, c)
212 CALL dbcsr_copy(d, c)
213 !
214 ! rotate G_i-1
215 IF (ASSOCIATED(g_old)) THEN
216 CALL dbcsr_multiply('N', 'N', 1.0_dp, g_old, u, 0.0_dp, c)
217 CALL dbcsr_copy(g_old, c)
218 END IF
219 !
220 CALL timestop(handle)
221 END SUBROUTINE qs_ot_on_the_fly_localize
222
223! **************************************************************************************************
224!> \brief ...
225!> \param qs_ot_env ...
226!> \param C_OLD ...
227!> \param C_TMP ...
228!> \param C_NEW ...
229!> \param P ...
230!> \param SC ...
231!> \param update ...
232! **************************************************************************************************
233 SUBROUTINE qs_ot_ref_chol(qs_ot_env, C_OLD, C_TMP, C_NEW, P, SC, update)
234 !
235 TYPE(qs_ot_type) :: qs_ot_env
236 TYPE(dbcsr_type) :: c_old, c_tmp, c_new, p, sc
237 LOGICAL, INTENT(IN) :: update
238
239 CHARACTER(len=*), PARAMETER :: routinen = 'qs_ot_ref_chol'
240
241 INTEGER :: handle, k, n
242
243 CALL timeset(routinen, handle)
244 !
245 CALL dbcsr_get_info(c_new, nfullrows_total=n, nfullcols_total=k)
246 !
247 ! P = U'*U
248 CALL cp_dbcsr_cholesky_decompose(p, k, qs_ot_env%para_env, qs_ot_env%blacs_env)
249 !
250 ! C_NEW = C_OLD*inv(U)
251 CALL cp_dbcsr_cholesky_restore(c_old, k, p, c_new, op="SOLVE", pos="RIGHT", &
252 transa="N", para_env=qs_ot_env%para_env, blacs_env=qs_ot_env%blacs_env)
253 !
254 ! Update SC if needed
255 IF (update) THEN
256 CALL cp_dbcsr_cholesky_restore(sc, k, p, c_tmp, op="SOLVE", pos="RIGHT", &
257 transa="N", para_env=qs_ot_env%para_env, blacs_env=qs_ot_env%blacs_env)
258 CALL dbcsr_copy(sc, c_tmp)
259 END IF
260 !
261 CALL timestop(handle)
262 END SUBROUTINE qs_ot_ref_chol
263
264! **************************************************************************************************
265!> \brief ...
266!> \param qs_ot_env ...
267!> \param C_OLD ...
268!> \param C_TMP ...
269!> \param C_NEW ...
270!> \param P ...
271!> \param SC ...
272!> \param update ...
273! **************************************************************************************************
274 SUBROUTINE qs_ot_ref_lwdn(qs_ot_env, C_OLD, C_TMP, C_NEW, P, SC, update)
275 !
276 TYPE(qs_ot_type) :: qs_ot_env
277 TYPE(dbcsr_type) :: c_old, c_tmp, c_new, p, sc
278 LOGICAL, INTENT(IN) :: update
279
280 CHARACTER(len=*), PARAMETER :: routinen = 'qs_ot_ref_lwdn'
281
282 INTEGER :: handle, i, k, n
283 REAL(dp), ALLOCATABLE, DIMENSION(:) :: eig, fun
284 TYPE(dbcsr_type), POINTER :: v, w
285
286 CALL timeset(routinen, handle)
287 !
288 CALL dbcsr_get_info(c_new, nfullrows_total=n, nfullcols_total=k)
289 !
290 v => qs_ot_env%buf1_k_k_nosym ! a buffer
291 w => qs_ot_env%buf2_k_k_nosym ! a buffer
292 ALLOCATE (eig(k), fun(k))
293 !
294 CALL cp_dbcsr_syevd(p, v, eig, qs_ot_env%para_env, qs_ot_env%blacs_env)
295 !
296 ! compute the P^(-1/2)
297 DO i = 1, k
298 IF (eig(i) <= 0.0_dp) THEN
299 cpabort("P not positive definite")
300 END IF
301 IF (eig(i) < 1.0e-8_dp) THEN
302 fun(i) = 0.0_dp
303 ELSE
304 fun(i) = 1.0_dp/sqrt(eig(i))
305 END IF
306 END DO
307 CALL dbcsr_copy(w, v)
308 CALL dbcsr_scale_by_vector(v, alpha=fun, side='right')
309 CALL dbcsr_multiply('N', 'T', 1.0_dp, w, v, 0.0_dp, p)
310 !
311 ! Update C
312 CALL dbcsr_multiply('N', 'N', 1.0_dp, c_old, p, 0.0_dp, c_new)
313 !
314 ! Update SC if needed
315 IF (update) THEN
316 CALL dbcsr_multiply('N', 'N', 1.0_dp, sc, p, 0.0_dp, c_tmp)
317 CALL dbcsr_copy(sc, c_tmp)
318 END IF
319 !
320 DEALLOCATE (eig, fun)
321 !
322 CALL timestop(handle)
323 END SUBROUTINE qs_ot_ref_lwdn
324
325! **************************************************************************************************
326!> \brief ...
327!> \param qs_ot_env ...
328!> \param C_OLD ...
329!> \param C_TMP ...
330!> \param C_NEW ...
331!> \param P ...
332!> \param SC ...
333!> \param norm_in ...
334!> \param update ...
335! **************************************************************************************************
336 SUBROUTINE qs_ot_ref_poly(qs_ot_env, C_OLD, C_TMP, C_NEW, P, SC, norm_in, update)
337 !
338 TYPE(qs_ot_type) :: qs_ot_env
339 TYPE(dbcsr_type), POINTER :: c_old, c_tmp, c_new, p
340 TYPE(dbcsr_type) :: sc
341 REAL(dp), INTENT(IN) :: norm_in
342 LOGICAL, INTENT(IN) :: update
343
344 CHARACTER(len=*), PARAMETER :: routinen = 'qs_ot_ref_poly'
345
346 INTEGER :: handle, irefine, k, n
347 LOGICAL :: quick_exit
348 REAL(dp) :: norm, norm_fro, norm_gct, occ_in, &
349 occ_out, rescale
350 TYPE(dbcsr_type), POINTER :: buf1, buf2, buf_nosym, ft, fy
351
352 CALL timeset(routinen, handle)
353 !
354 CALL dbcsr_get_info(c_new, nfullrows_total=n, nfullcols_total=k)
355 !
356 buf_nosym => qs_ot_env%buf1_k_k_nosym ! a buffer
357 buf1 => qs_ot_env%buf1_k_k_sym ! a buffer
358 buf2 => qs_ot_env%buf2_k_k_sym ! a buffer
359 fy => qs_ot_env%buf3_k_k_sym ! a buffer
360 ft => qs_ot_env%buf4_k_k_sym ! a buffer
361 !
362 ! initialize the norm (already computed in qs_ot_get_orbitals_ref)
363 norm = norm_in
364 !
365 ! can we do a quick exit?
366 quick_exit = .false.
367 IF (norm < qs_ot_env%settings%eps_irac_quick_exit) quick_exit = .true.
368 !
369 ! lets refine
370 rescale = 1.0_dp
371 DO irefine = 1, qs_ot_env%settings%max_irac
372 !
373 ! rescaling
374 IF (norm > 1.0_dp) THEN
375 CALL dbcsr_scale(p, 1.0_dp/norm)
376 rescale = rescale/sqrt(norm)
377 END IF
378 !
379 ! get the refinement polynomial
380 CALL qs_ot_refine(p, fy, buf1, buf2, qs_ot_env%settings%irac_degree, &
381 qs_ot_env%settings%eps_irac_filter_matrix)
382 !
383 ! collect the transformation
384 IF (irefine == 1) THEN
385 CALL dbcsr_copy(ft, fy, name='FT')
386 ELSE
387 CALL dbcsr_multiply('N', 'N', 1.0_dp, ft, fy, 0.0_dp, buf1)
388 IF (qs_ot_env%settings%eps_irac_filter_matrix > 0.0_dp) THEN
389 occ_in = dbcsr_get_occupation(buf1)
390 CALL dbcsr_filter(buf1, qs_ot_env%settings%eps_irac_filter_matrix)
391 occ_out = dbcsr_get_occupation(buf1)
392 END IF
393 CALL dbcsr_copy(ft, buf1, name='FT')
394 END IF
395 !
396 ! quick exit if possible
397 IF (quick_exit) THEN
398 EXIT
399 END IF
400 !
401 ! P = FY^T * P * FY
402 CALL dbcsr_multiply('N', 'N', 1.0_dp, p, fy, 0.0_dp, buf_nosym)
403 IF (qs_ot_env%settings%eps_irac_filter_matrix > 0.0_dp) THEN
404 occ_in = dbcsr_get_occupation(buf_nosym)
405 CALL dbcsr_filter(buf_nosym, qs_ot_env%settings%eps_irac_filter_matrix)
406 occ_out = dbcsr_get_occupation(buf_nosym)
407 END IF
408 CALL dbcsr_multiply('N', 'N', 1.0_dp, fy, buf_nosym, 0.0_dp, p)
409 IF (qs_ot_env%settings%eps_irac_filter_matrix > 0.0_dp) THEN
410 occ_in = dbcsr_get_occupation(p)
411 CALL dbcsr_filter(p, qs_ot_env%settings%eps_irac_filter_matrix)
412 occ_out = dbcsr_get_occupation(p)
413 END IF
414 !
415 ! check ||P-1||_gct
416 CALL dbcsr_add_on_diag(p, -1.0_dp)
417 norm_fro = dbcsr_frobenius_norm(p)
418 norm_gct = dbcsr_gershgorin_norm(p)
419 CALL dbcsr_add_on_diag(p, 1.0_dp)
420 norm = min(norm_gct, norm_fro)
421 !
422 ! printing
423 !
424 ! blows up
425 IF (norm > 1.0e10_dp) THEN
426 CALL cp_abort(__location__, &
427 "Refinement blows up! "// &
428 "We need you to improve the code, please post your input on "// &
429 "the forum https://www.cp2k.org/")
430 END IF
431 !
432 ! can we do a quick exit next step?
433 IF (norm < qs_ot_env%settings%eps_irac_quick_exit) quick_exit = .true.
434 !
435 ! are we done?
436 IF (norm < qs_ot_env%settings%eps_irac) EXIT
437 !
438 END DO
439 !
440 ! C_NEW = C_NEW * FT * rescale
441 CALL dbcsr_multiply('N', 'N', rescale, c_old, ft, 0.0_dp, c_new)
442 IF (qs_ot_env%settings%eps_irac_filter_matrix > 0.0_dp) THEN
443 occ_in = dbcsr_get_occupation(c_new)
444 CALL dbcsr_filter(c_new, qs_ot_env%settings%eps_irac_filter_matrix)
445 occ_out = dbcsr_get_occupation(c_new)
446 END IF
447 !
448 ! update SC = SC * FY * rescale
449 IF (update) THEN
450 CALL dbcsr_multiply('N', 'N', rescale, sc, ft, 0.0_dp, c_tmp)
451 IF (qs_ot_env%settings%eps_irac_filter_matrix > 0.0_dp) THEN
452 occ_in = dbcsr_get_occupation(c_tmp)
453 CALL dbcsr_filter(c_tmp, qs_ot_env%settings%eps_irac_filter_matrix)
454 occ_out = dbcsr_get_occupation(c_tmp)
455 END IF
456 CALL dbcsr_copy(sc, c_tmp)
457 END IF
458 !
459 CALL timestop(handle)
460 END SUBROUTINE qs_ot_ref_poly
461
462! **************************************************************************************************
463!> \brief ...
464!> \param qs_ot_env1 ...
465!> \return ...
466! **************************************************************************************************
467 FUNCTION qs_ot_ref_update(qs_ot_env1) RESULT(update)
468 !
469 TYPE(qs_ot_type) :: qs_ot_env1
470 LOGICAL :: update
471
472 update = .false.
473 SELECT CASE (qs_ot_env1%settings%ot_method)
474 CASE ("CG")
475 SELECT CASE (qs_ot_env1%settings%line_search_method)
476 CASE ("2PNT")
477 IF (qs_ot_env1%line_search_count == 2) update = .true.
478 CASE DEFAULT
479 cpabort("NYI")
480 END SELECT
481 CASE ("DIIS")
482 update = .true.
483 CASE DEFAULT
484 cpabort("NYI")
485 END SELECT
486 END FUNCTION qs_ot_ref_update
487
488! **************************************************************************************************
489!> \brief ...
490!> \param qs_ot_env1 ...
491!> \param norm_in ...
492!> \param ortho_irac ...
493! **************************************************************************************************
494 SUBROUTINE qs_ot_ref_decide(qs_ot_env1, norm_in, ortho_irac)
495 !
496 TYPE(qs_ot_type) :: qs_ot_env1
497 REAL(dp), INTENT(IN) :: norm_in
498 CHARACTER(LEN=*), INTENT(INOUT) :: ortho_irac
499
500 ortho_irac = qs_ot_env1%settings%ortho_irac
501 IF (norm_in < qs_ot_env1%settings%eps_irac_switch) ortho_irac = "POLY"
502 END SUBROUTINE qs_ot_ref_decide
503
504! **************************************************************************************************
505!> \brief ...
506!> \param matrix_c ...
507!> \param matrix_s ...
508!> \param matrix_x ...
509!> \param matrix_sx ...
510!> \param matrix_gx_old ...
511!> \param matrix_dx ...
512!> \param qs_ot_env ...
513!> \param qs_ot_env1 ...
514! **************************************************************************************************
515 SUBROUTINE qs_ot_get_orbitals_ref(matrix_c, matrix_s, matrix_x, matrix_sx, &
516 matrix_gx_old, matrix_dx, qs_ot_env, qs_ot_env1)
517 !
518 TYPE(dbcsr_type), POINTER :: matrix_c, matrix_s, matrix_x, matrix_sx, &
519 matrix_gx_old, matrix_dx
520 TYPE(qs_ot_type) :: qs_ot_env, qs_ot_env1
521
522 CHARACTER(len=*), PARAMETER :: routinen = 'qs_ot_get_orbitals_ref'
523
524 CHARACTER(LEN=4) :: ortho_irac
525 INTEGER :: handle, k, n
526 LOGICAL :: on_the_fly_loc, update
527 REAL(dp) :: norm, norm_fro, norm_gct, occ_in, occ_out
528 TYPE(dbcsr_type), POINTER :: c_new, c_old, c_tmp, d, g_old, p, s, sc
529
530 CALL timeset(routinen, handle)
531
532 CALL dbcsr_get_info(matrix_c, nfullrows_total=n, nfullcols_total=k)
533 !
534 c_new => matrix_c
535 c_old => matrix_x ! need to be carefully updated for the gradient !
536 sc => matrix_sx ! need to be carefully updated for the gradient !
537 g_old => matrix_gx_old ! need to be carefully updated for localization !
538 d => matrix_dx ! need to be carefully updated for localization !
539 s => matrix_s
540
541 p => qs_ot_env%p_k_k_sym ! a buffer
542 c_tmp => qs_ot_env%buf1_n_k ! a buffer
543 !
544 ! do we need to update C_OLD and SC?
545 update = qs_ot_ref_update(qs_ot_env1)
546 !
547 ! do we want to on the fly localize?
548 ! for the moment this is set from the input,
549 ! later we might want to localize every n-step or
550 ! when the sparsity increases...
551 on_the_fly_loc = qs_ot_env1%settings%on_the_fly_loc
552 !
553 ! compute SC = S*C
554 IF (ASSOCIATED(s)) THEN
555 CALL dbcsr_multiply('N', 'N', 1.0_dp, s, c_old, 0.0_dp, sc)
556 IF (qs_ot_env1%settings%eps_irac_filter_matrix > 0.0_dp) THEN
557 occ_in = dbcsr_get_occupation(sc)
558 CALL dbcsr_filter(sc, qs_ot_env1%settings%eps_irac_filter_matrix)
559 occ_out = dbcsr_get_occupation(sc)
560 END IF
561 ELSE
562 CALL dbcsr_copy(sc, c_old)
563 END IF
564 !
565 ! compute P = C'*SC
566 CALL dbcsr_multiply('T', 'N', 1.0_dp, c_old, sc, 0.0_dp, p)
567 IF (qs_ot_env1%settings%eps_irac_filter_matrix > 0.0_dp) THEN
568 occ_in = dbcsr_get_occupation(p)
569 CALL dbcsr_filter(p, qs_ot_env1%settings%eps_irac_filter_matrix)
570 occ_out = dbcsr_get_occupation(p)
571 END IF
572 !
573 ! check ||P-1||_f and ||P-1||_gct
574 CALL dbcsr_add_on_diag(p, -1.0_dp)
575 norm_fro = dbcsr_frobenius_norm(p)
576 norm_gct = dbcsr_gershgorin_norm(p)
577 CALL dbcsr_add_on_diag(p, 1.0_dp)
578 norm = min(norm_gct, norm_fro)
579 CALL qs_ot_ref_decide(qs_ot_env1, norm, ortho_irac)
580 !
581 ! select the orthogonality method
582 SELECT CASE (ortho_irac)
583 CASE ("CHOL")
584 CALL qs_ot_ref_chol(qs_ot_env, c_old, c_tmp, c_new, p, sc, update)
585 CASE ("LWDN")
586 CALL qs_ot_ref_lwdn(qs_ot_env, c_old, c_tmp, c_new, p, sc, update)
587 CASE ("POLY")
588 CALL qs_ot_ref_poly(qs_ot_env, c_old, c_tmp, c_new, p, sc, norm, update)
589 CASE DEFAULT
590 cpabort("Wrong argument")
591 END SELECT
592 !
593 ! We update the C_i+1 and localization
594 IF (update) THEN
595 IF (on_the_fly_loc) THEN
596 CALL qs_ot_on_the_fly_localize(qs_ot_env, c_new, sc, g_old, d)
597 END IF
598 CALL dbcsr_copy(c_old, c_new)
599 END IF
600 !
601 CALL timestop(handle)
602 END SUBROUTINE qs_ot_get_orbitals_ref
603
604! **************************************************************************************************
605!> \brief refinement polynomial of degree 2,3 and 4 (PRB 70, 193102 (2004))
606!> \param P ...
607!> \param FY ...
608!> \param P2 ...
609!> \param T ...
610!> \param irac_degree ...
611!> \param eps_irac_filter_matrix ...
612! **************************************************************************************************
613 SUBROUTINE qs_ot_refine(P, FY, P2, T, irac_degree, eps_irac_filter_matrix)
614 TYPE(dbcsr_type), INTENT(inout) :: p, fy, p2, t
615 INTEGER, INTENT(in) :: irac_degree
616 REAL(dp), INTENT(in) :: eps_irac_filter_matrix
617
618 CHARACTER(len=*), PARAMETER :: routinen = 'qs_ot_refine'
619
620 INTEGER :: handle, k
621 REAL(dp) :: occ_in, occ_out, r
622
623 CALL timeset(routinen, handle)
624
625 CALL dbcsr_get_info(p, nfullcols_total=k)
626 SELECT CASE (irac_degree)
627 CASE (2)
628 ! C_out = C_in * ( 15/8 * I - 10/8 * P + 3/8 * P^2)
629 r = 3.0_dp/8.0_dp
630 CALL dbcsr_multiply('N', 'N', r, p, p, 0.0_dp, fy)
631 IF (eps_irac_filter_matrix > 0.0_dp) THEN
632 occ_in = dbcsr_get_occupation(fy)
633 CALL dbcsr_filter(fy, eps_irac_filter_matrix)
634 occ_out = dbcsr_get_occupation(fy)
635 END IF
636 r = -10.0_dp/8.0_dp
637 CALL dbcsr_add(fy, p, alpha_scalar=1.0_dp, beta_scalar=r)
638 r = 15.0_dp/8.0_dp
639 CALL dbcsr_add_on_diag(fy, alpha=r)
640 CASE (3)
641 ! C_out = C_in * ( 35/16 * I - 35/16 * P + 21/16 * P^2 - 5/16 P^3)
642 CALL dbcsr_multiply('N', 'N', 1.0_dp, p, p, 0.0_dp, p2)
643 IF (eps_irac_filter_matrix > 0.0_dp) THEN
644 occ_in = dbcsr_get_occupation(p2)
645 CALL dbcsr_filter(p2, eps_irac_filter_matrix)
646 occ_out = dbcsr_get_occupation(p2)
647 END IF
648 r = -5.0_dp/16.0_dp
649 CALL dbcsr_multiply('N', 'N', r, p2, p, 0.0_dp, fy)
650 IF (eps_irac_filter_matrix > 0.0_dp) THEN
651 occ_in = dbcsr_get_occupation(fy)
652 CALL dbcsr_filter(fy, eps_irac_filter_matrix)
653 occ_out = dbcsr_get_occupation(fy)
654 END IF
655 r = 21.0_dp/16.0_dp
656 CALL dbcsr_add(fy, p2, alpha_scalar=1.0_dp, beta_scalar=r)
657 r = -35.0_dp/16.0_dp
658 CALL dbcsr_add(fy, p, alpha_scalar=1.0_dp, beta_scalar=r)
659 r = 35.0_dp/16.0_dp
660 CALL dbcsr_add_on_diag(fy, alpha=r)
661 CASE (4)
662 ! C_out = C_in * ( 315/128 * I - 420/128 * P + 378/128 * P^2 - 180/128 P^3 + 35/128 P^4 )
663 ! = C_in * ( 315/128 * I - 420/128 * P + 378/128 * P^2 + ( - 180/128 * P + 35/128 * P^2 ) * P^2 )
664 CALL dbcsr_multiply('N', 'N', 1.0_dp, p, p, 0.0_dp, p2) ! P^2
665 IF (eps_irac_filter_matrix > 0.0_dp) THEN
666 occ_in = dbcsr_get_occupation(p2)
667 CALL dbcsr_filter(p2, eps_irac_filter_matrix)
668 occ_out = dbcsr_get_occupation(p2)
669 END IF
670 r = -180.0_dp/128.0_dp
671 CALL dbcsr_add(t, p, alpha_scalar=0.0_dp, beta_scalar=r) ! T=-180/128*P
672 r = 35.0_dp/128.0_dp
673 CALL dbcsr_add(t, p2, alpha_scalar=1.0_dp, beta_scalar=r) ! T=T+35/128*P^2
674 CALL dbcsr_multiply('N', 'N', 1.0_dp, t, p2, 0.0_dp, fy) ! Y=T*P^2
675 IF (eps_irac_filter_matrix > 0.0_dp) THEN
676 occ_in = dbcsr_get_occupation(fy)
677 CALL dbcsr_filter(fy, eps_irac_filter_matrix)
678 occ_out = dbcsr_get_occupation(fy)
679 END IF
680 r = 378.0_dp/128.0_dp
681 CALL dbcsr_add(fy, p2, alpha_scalar=1.0_dp, beta_scalar=r) ! Y=Y+378/128*P^2
682 r = -420.0_dp/128.0_dp
683 CALL dbcsr_add(fy, p, alpha_scalar=1.0_dp, beta_scalar=r) ! Y=Y-420/128*P
684 r = 315.0_dp/128.0_dp
685 CALL dbcsr_add_on_diag(fy, alpha=r) ! Y=Y+315/128*I
686 CASE DEFAULT
687 cpabort("This irac_order NYI")
688 END SELECT
689 CALL timestop(handle)
690 END SUBROUTINE qs_ot_refine
691
692! **************************************************************************************************
693!> \brief ...
694!> \param matrix_hc ...
695!> \param matrix_x ...
696!> \param matrix_sx ...
697!> \param matrix_gx ...
698!> \param qs_ot_env ...
699! **************************************************************************************************
700 SUBROUTINE qs_ot_get_derivative_ref(matrix_hc, matrix_x, matrix_sx, matrix_gx, &
701 qs_ot_env)
702 TYPE(dbcsr_type), POINTER :: matrix_hc, matrix_x, matrix_sx, matrix_gx
703 TYPE(qs_ot_type) :: qs_ot_env
704
705 CHARACTER(len=*), PARAMETER :: routinen = 'qs_ot_get_derivative_ref'
706
707 INTEGER :: handle, k, n
708 REAL(dp) :: occ_in, occ_out
709 TYPE(dbcsr_type), POINTER :: c, chc, g, g_dp, hc, sc
710
711 CALL timeset(routinen, handle)
712
713 CALL dbcsr_get_info(matrix_x, nfullrows_total=n, nfullcols_total=k)
714 !
715 c => matrix_x ! NBsf*NOcc
716 sc => matrix_sx ! NBsf*NOcc need to be up2date
717 hc => matrix_hc ! NBsf*NOcc
718 g => matrix_gx ! NBsf*NOcc
719 chc => qs_ot_env%buf1_k_k_sym ! buffer
720 g_dp => qs_ot_env%buf1_n_k_dp ! buffer dp
721
722 ! C'*(H*C)
723 CALL dbcsr_multiply('T', 'N', 1.0_dp, c, hc, 0.0_dp, chc)
724 IF (qs_ot_env%settings%eps_irac_filter_matrix > 0.0_dp) THEN
725 occ_in = dbcsr_get_occupation(chc)
726 CALL dbcsr_filter(chc, qs_ot_env%settings%eps_irac_filter_matrix)
727 occ_out = dbcsr_get_occupation(chc)
728 END IF
729 ! (S*C)*(C'*H*C)
730 CALL dbcsr_multiply('N', 'N', 1.0_dp, sc, chc, 0.0_dp, g)
731 IF (qs_ot_env%settings%eps_irac_filter_matrix > 0.0_dp) THEN
732 occ_in = dbcsr_get_occupation(g)
733 CALL dbcsr_filter(g, qs_ot_env%settings%eps_irac_filter_matrix)
734 occ_out = dbcsr_get_occupation(g)
735 END IF
736 ! G = 2*(1-S*C*C')*H*C
737 CALL dbcsr_add(g, hc, alpha_scalar=-1.0_dp, beta_scalar=1.0_dp)
738 !
739 CALL timestop(handle)
740 END SUBROUTINE qs_ot_get_derivative_ref
741
742! **************************************************************************************************
743!> \brief computes p=x*S*x and the matrix functionals related matrices
744!> \param matrix_x ...
745!> \param matrix_sx ...
746!> \param qs_ot_env ...
747! **************************************************************************************************
748 SUBROUTINE qs_ot_get_p(matrix_x, matrix_sx, qs_ot_env)
749
750 TYPE(dbcsr_type), POINTER :: matrix_x, matrix_sx
751 TYPE(qs_ot_type) :: qs_ot_env
752
753 CHARACTER(len=*), PARAMETER :: routinen = 'qs_ot_get_p'
754 REAL(kind=dp), PARAMETER :: rone = 1.0_dp, rzero = 0.0_dp
755
756 INTEGER :: handle, k, max_iter, n
757 LOGICAL :: converged
758 REAL(kind=dp) :: max_ev, min_ev, threshold
759
760 CALL timeset(routinen, handle)
761
762 CALL dbcsr_get_info(matrix_x, nfullrows_total=n, nfullcols_total=k)
763
764 ! get the overlap
765 CALL dbcsr_multiply('T', 'N', rone, matrix_x, matrix_sx, rzero, &
766 qs_ot_env%matrix_p)
767
768 ! get an upper bound for the largest eigenvalue
769 ! try using lancos first and fall back to gershgorin norm if it fails
770 max_iter = 30; threshold = 1.0e-03_dp
771 CALL arnoldi_extremal(qs_ot_env%matrix_p, max_ev, min_ev, converged, threshold, max_iter)
772 qs_ot_env%largest_eval_upper_bound = max(max_ev, abs(min_ev))
773
774 IF (.NOT. converged) qs_ot_env%largest_eval_upper_bound = dbcsr_gershgorin_norm(qs_ot_env%matrix_p)
775 CALL decide_strategy(qs_ot_env)
776 IF (qs_ot_env%do_taylor) THEN
777 CALL qs_ot_p2m_taylor(qs_ot_env)
778 ELSE
779 CALL qs_ot_p2m_diag(qs_ot_env)
780 END IF
781
782 IF (qs_ot_env%settings%do_rotation) THEN
783 CALL qs_ot_generate_rotation(qs_ot_env)
784 END IF
785
786 CALL timestop(handle)
787
788 END SUBROUTINE qs_ot_get_p
789
790! **************************************************************************************************
791!> \brief computes the rotation matrix rot_mat_u that is associated to a given
792!> rot_mat_x using rot_mat_u=exp(rot_mat_x)
793!> \param qs_ot_env a valid qs_ot_env
794!> \par History
795!> 08.2004 created [Joost VandeVondele]
796!> 12.2024 Rewrite to use only real matrices [Ole Schuett]
797! **************************************************************************************************
798 SUBROUTINE qs_ot_generate_rotation(qs_ot_env)
799
800 TYPE(qs_ot_type) :: qs_ot_env
801
802 CHARACTER(len=*), PARAMETER :: routinen = 'qs_ot_generate_rotation'
803
804 INTEGER :: handle, k
805 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: exp_evals_im, exp_evals_re
806 TYPE(dbcsr_type) :: buf_1, buf_2
807
808 CALL timeset(routinen, handle)
809
810 CALL dbcsr_get_info(qs_ot_env%rot_mat_x, nfullrows_total=k)
811
812 IF (k /= 0) THEN
813 ! We want to compute: rot_mat_u = exp(i*rot_mat_x)
814
815 ! Diagonalize: matrix = i*rot_mat_x.
816 ! Note that matrix is imaginary and hermitian because rot_mat_x is real and anti-symmetric.
817 CALL cp_dbcsr_heevd(matrix_im=qs_ot_env%rot_mat_x, & ! matrix_re omitted because it's zero
818 eigenvectors_re=qs_ot_env%rot_mat_evec_re, &
819 eigenvectors_im=qs_ot_env%rot_mat_evec_im, &
820 eigenvalues=qs_ot_env%rot_mat_evals, &
821 para_env=qs_ot_env%para_env, &
822 blacs_env=qs_ot_env%blacs_env)
823
824 ! Compute: exp_evals = EXP(-i*rot_mat_evals)
825 ALLOCATE (exp_evals_re(k), exp_evals_im(k))
826 exp_evals_re(:) = cos(-qs_ot_env%rot_mat_evals(:))
827 exp_evals_im(:) = sin(-qs_ot_env%rot_mat_evals(:))
828
829 ! Compute: rot_mat_u = \sum_ij exp_evals_ij * |rot_mat_evec_i> <rot_mat_evec_j|
830 ! Note that we need only two matrix multiplications because rot_mat_u is real.
831 CALL dbcsr_copy(buf_1, qs_ot_env%rot_mat_evec_re, name="buf_1")
832 CALL dbcsr_scale_by_vector(buf_1, alpha=exp_evals_re, side='right')
833 CALL dbcsr_copy(buf_2, qs_ot_env%rot_mat_evec_im, name="buf_2")
834 CALL dbcsr_scale_by_vector(buf_2, alpha=exp_evals_im, side='right')
835 CALL dbcsr_add(buf_1, buf_2, alpha_scalar=+1.0_dp, beta_scalar=-1.0_dp)
836 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)
837
838 CALL dbcsr_copy(buf_1, qs_ot_env%rot_mat_evec_im)
839 CALL dbcsr_scale_by_vector(buf_1, alpha=exp_evals_re, side='right')
840 CALL dbcsr_copy(buf_2, qs_ot_env%rot_mat_evec_re)
841 CALL dbcsr_scale_by_vector(buf_2, alpha=exp_evals_im, side='right')
842 CALL dbcsr_add(buf_1, buf_2, alpha_scalar=+1.0_dp, beta_scalar=+1.0_dp)
843 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)
844
845 ! Clean up.
846 CALL dbcsr_release(buf_1)
847 CALL dbcsr_release(buf_2)
848 DEALLOCATE (exp_evals_re, exp_evals_im)
849 END IF
850
851 CALL timestop(handle)
852
853 END SUBROUTINE qs_ot_generate_rotation
854
855! **************************************************************************************************
856!> \brief computes the derivative fields with respect to rot_mat_x
857!> \param qs_ot_env valid qs_ot_env. In particular qs_ot_generate_rotation has to be called before
858!> and the rot_mat_dedu matrix has to be up to date
859!> \par History
860!> 08.2004 created [ Joost VandeVondele ]
861!> 12.2024 Rewrite to use only real matrices [Ole Schuett]
862! **************************************************************************************************
863 SUBROUTINE qs_ot_rot_mat_derivative(qs_ot_env)
864 TYPE(qs_ot_type) :: qs_ot_env
865
866 CHARACTER(len=*), PARAMETER :: routinen = 'qs_ot_rot_mat_derivative'
867
868 INTEGER :: handle, i, j, k
869 REAL(kind=dp) :: e1, e2
870 TYPE(dbcsr_type) :: outer_deriv_re, outer_deriv_im, mat_buf, &
871 inner_deriv_re, inner_deriv_im
872 TYPE(dbcsr_iterator_type) :: iter
873 INTEGER, DIMENSION(:), POINTER :: row_blk_offset, col_blk_offset
874 REAL(dp), DIMENSION(:, :), POINTER :: block_in_re, block_in_im, block_out_re, block_out_im
875 INTEGER :: row, col
876 LOGICAL :: found
877 COMPLEX(dp) :: cval_in, cval_out
878 TYPE(dbcsr_distribution_type) :: dist
879
880 CALL timeset(routinen, handle)
881
882 CALL dbcsr_get_info(qs_ot_env%rot_mat_u, nfullrows_total=k)
883 IF (k /= 0) THEN
884 CALL dbcsr_copy(qs_ot_env%matrix_buf1, qs_ot_env%rot_mat_dedu)
885 ! now we get to the derivative wrt the antisymmetric matrix rot_mat_x
886 CALL dbcsr_copy(mat_buf, qs_ot_env%rot_mat_dedu, "mat_buf")
887
888 ! inner_deriv_ij = <rot_mat_evec_i| rot_mat_dedu |rot_mat_evec_j>
889 CALL dbcsr_copy(inner_deriv_re, qs_ot_env%rot_mat_dedu, "inner_deriv_re") ! TODO just create
890 CALL dbcsr_copy(inner_deriv_im, qs_ot_env%rot_mat_dedu, "inner_deriv_im") ! TODO just create
891
892 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)
893 CALL dbcsr_multiply('T', 'N', +1.0_dp, qs_ot_env%rot_mat_evec_im, mat_buf, 0.0_dp, inner_deriv_re)
894 CALL dbcsr_multiply('T', 'N', +1.0_dp, qs_ot_env%rot_mat_evec_re, mat_buf, 0.0_dp, inner_deriv_im)
895
896 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)
897 CALL dbcsr_multiply('T', 'N', +1.0_dp, qs_ot_env%rot_mat_evec_re, mat_buf, 1.0_dp, inner_deriv_re)
898 CALL dbcsr_multiply('T', 'N', -1.0_dp, qs_ot_env%rot_mat_evec_im, mat_buf, 1.0_dp, inner_deriv_im)
899
900 ! outer_deriv_ij = cint(eval_i, eval_j) * inner_deriv_ij
901 CALL dbcsr_copy(outer_deriv_re, qs_ot_env%rot_mat_dedu, "outer_deriv_re") ! TODO just create
902 CALL dbcsr_copy(outer_deriv_im, qs_ot_env%rot_mat_dedu, "outer_deriv_im") ! TODO just create
903
904 CALL dbcsr_get_info(qs_ot_env%rot_mat_dedu, row_blk_offset=row_blk_offset, col_blk_offset=col_blk_offset)
905 CALL dbcsr_iterator_start(iter, qs_ot_env%rot_mat_dedu)
906 DO WHILE (dbcsr_iterator_blocks_left(iter))
907 CALL dbcsr_iterator_next_block(iter, row, col)
908 CALL dbcsr_get_block_p(inner_deriv_re, row, col, block_in_re, found)
909 CALL dbcsr_get_block_p(inner_deriv_im, row, col, block_in_im, found)
910 CALL dbcsr_get_block_p(outer_deriv_re, row, col, block_out_re, found)
911 CALL dbcsr_get_block_p(outer_deriv_im, row, col, block_out_im, found)
912
913 DO i = 1, SIZE(block_in_re, 1)
914 DO j = 1, SIZE(block_in_re, 2)
915 e1 = qs_ot_env%rot_mat_evals(row_blk_offset(row) + i - 1)
916 e2 = qs_ot_env%rot_mat_evals(col_blk_offset(col) + j - 1)
917 cval_in = cmplx(block_in_re(i, j), block_in_im(i, j), dp)
918 cval_out = cval_in*cint(e1, e2)
919 block_out_re(i, j) = real(cval_out)
920 block_out_im(i, j) = aimag(cval_out)
921 END DO
922 END DO
923 END DO
924 CALL dbcsr_iterator_stop(iter)
925 CALL dbcsr_release(inner_deriv_re)
926 CALL dbcsr_release(inner_deriv_im)
927
928 ! Compute: matrix_buf1 = \sum_i outer_deriv_ij * |rot_mat_evec_i> <rot_mat_evec_j|
929 CALL dbcsr_multiply('N', 'N', +1.0_dp, qs_ot_env%rot_mat_evec_re, outer_deriv_re, 0.0_dp, mat_buf)
930 CALL dbcsr_multiply('N', 'N', -1.0_dp, qs_ot_env%rot_mat_evec_im, outer_deriv_im, 1.0_dp, mat_buf)
931 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)
932
933 CALL dbcsr_multiply('N', 'N', +1.0_dp, qs_ot_env%rot_mat_evec_re, outer_deriv_im, 0.0_dp, mat_buf)
934 CALL dbcsr_multiply('N', 'N', +1.0_dp, qs_ot_env%rot_mat_evec_im, outer_deriv_re, 1.0_dp, mat_buf)
935 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)
936
937 ! Account for anti-symmetry of rot_mat_x.
938 CALL dbcsr_get_info(qs_ot_env%matrix_buf3, distribution=dist)
939 CALL dbcsr_transposed(qs_ot_env%matrix_buf2, qs_ot_env%matrix_buf1, &
940 shallow_data_copy=.false., use_distribution=dist, &
941 transpose_distribution=.false.)
942
943 ! rot_mat_gx = matrix_buf1^T - matrix_buf1
944 CALL dbcsr_add(qs_ot_env%matrix_buf1, qs_ot_env%matrix_buf2, alpha_scalar=-1.0_dp, beta_scalar=+1.0_dp)
945 CALL dbcsr_copy(qs_ot_env%rot_mat_gx, qs_ot_env%matrix_buf1)
946
947 CALL dbcsr_release(mat_buf)
948 CALL dbcsr_release(outer_deriv_re)
949 CALL dbcsr_release(outer_deriv_im)
950 END IF
951 CALL timestop(handle)
952 CONTAINS
953
954! **************************************************************************************************
955!> \brief ...
956!> \param e1 ...
957!> \param e2 ...
958!> \return ...
959! **************************************************************************************************
960 FUNCTION cint(e1, e2)
961 REAL(kind=dp) :: e1, e2
962 COMPLEX(KIND=dp) :: cint
963
964 COMPLEX(KIND=dp) :: l1, l2, x
965 INTEGER :: i
966
967 l1 = (0.0_dp, -1.0_dp)*e1
968 l2 = (0.0_dp, -1.0_dp)*e2
969 IF (abs(l1 - l2) > 0.5_dp) THEN
970 cint = (exp(l1) - exp(l2))/(l1 - l2)
971 ELSE
972 x = 1.0_dp
973 cint = 0.0_dp
974 DO i = 1, 16
975 cint = cint + x
976 x = x*(l1 - l2)/real(i + 1, kind=dp)
977 END DO
978 cint = cint*exp(l2)
979 END IF
980 END FUNCTION cint
981 END SUBROUTINE qs_ot_rot_mat_derivative
982
983! **************************************************************************************************
984!> \brief decide strategy
985!> tries to decide if the taylor expansion of cos(sqrt(xsx)) converges rapidly enough
986!> to make a taylor expansion of the functions cos(sqrt(xsx)) and sin(sqrt(xsx))/sqrt(xsx)
987!> and their derivatives faster than their computation based on diagonalization since xsx can
988!> be very small, especially during dynamics, only a few terms might indeed be needed we find
989!> the necessary order N to have largest_eval_upper_bound**(N+1)/(2(N+1))! < eps_taylor
990!> \param qs_ot_env ...
991! **************************************************************************************************
992 SUBROUTINE decide_strategy(qs_ot_env)
993 TYPE(qs_ot_type) :: qs_ot_env
994
995 INTEGER :: n
996 REAL(kind=dp) :: num_error
997
998 qs_ot_env%do_taylor = .false.
999 n = 0
1000 num_error = qs_ot_env%largest_eval_upper_bound/(2.0_dp)
1001 DO WHILE (num_error > qs_ot_env%settings%eps_taylor .AND. n <= qs_ot_env%settings%max_taylor)
1002 n = n + 1
1003 num_error = num_error*qs_ot_env%largest_eval_upper_bound/real((2*n + 1)*(2*n + 2), kind=dp)
1004 END DO
1005 qs_ot_env%taylor_order = n
1006 IF (qs_ot_env%taylor_order <= qs_ot_env%settings%max_taylor) THEN
1007 qs_ot_env%do_taylor = .true.
1008 END IF
1009
1010 END SUBROUTINE decide_strategy
1011
1012! **************************************************************************************************
1013!> \brief c=(c0*cos(p^0.5)+x*sin(p^0.5)*p^(-0.5)) x rot_mat_u
1014!> this assumes that x is already ortho to S*C0, and that p is x*S*x
1015!> rot_mat_u is an optional rotation matrix
1016!> \param matrix_c ...
1017!> \param matrix_x ...
1018!> \param qs_ot_env ...
1019! **************************************************************************************************
1020 SUBROUTINE qs_ot_get_orbitals(matrix_c, matrix_x, qs_ot_env)
1021
1022 TYPE(dbcsr_type), POINTER :: matrix_c, matrix_x
1023 TYPE(qs_ot_type) :: qs_ot_env
1024
1025 CHARACTER(len=*), PARAMETER :: routinen = 'qs_ot_get_orbitals'
1026 REAL(kind=dp), PARAMETER :: rone = 1.0_dp, rzero = 0.0_dp
1027
1028 INTEGER :: handle, k, n
1029 TYPE(dbcsr_type), POINTER :: matrix_kk
1030
1031 CALL timeset(routinen, handle)
1032
1033 CALL dbcsr_get_info(matrix_x, nfullrows_total=n, nfullcols_total=k)
1034
1035 ! rotate the multiplying matrices cosp and sinp instead of the result,
1036 ! this should be cheaper for large basis sets
1037 IF (qs_ot_env%settings%do_rotation) THEN
1038 matrix_kk => qs_ot_env%matrix_buf1
1039 CALL dbcsr_multiply('N', 'N', rone, qs_ot_env%matrix_cosp, &
1040 qs_ot_env%rot_mat_u, rzero, matrix_kk)
1041 ELSE
1042 matrix_kk => qs_ot_env%matrix_cosp
1043 END IF
1044
1045 CALL dbcsr_multiply('N', 'N', rone, qs_ot_env%matrix_c0, matrix_kk, &
1046 rzero, matrix_c)
1047
1048 IF (qs_ot_env%settings%do_rotation) THEN
1049 matrix_kk => qs_ot_env%matrix_buf1
1050 CALL dbcsr_multiply('N', 'N', rone, qs_ot_env%matrix_sinp, &
1051 qs_ot_env%rot_mat_u, rzero, matrix_kk)
1052 ELSE
1053 matrix_kk => qs_ot_env%matrix_sinp
1054 END IF
1055 CALL dbcsr_multiply('N', 'N', rone, matrix_x, matrix_kk, &
1056 rone, matrix_c)
1057
1058 CALL timestop(handle)
1059
1060 END SUBROUTINE qs_ot_get_orbitals
1061
1062! **************************************************************************************************
1063!> \brief this routines computes dE/dx=dx, with dx ortho to sc0
1064!> needs dE/dC=hc,C0,X,SX,p
1065!> if preconditioned it will not be the derivative, but the lagrangian multiplier
1066!> is changed so that P*dE/dx is the right derivative (i.e. in the allowed subspace)
1067!> \param matrix_hc ...
1068!> \param matrix_x ...
1069!> \param matrix_sx ...
1070!> \param matrix_gx ...
1071!> \param qs_ot_env ...
1072! **************************************************************************************************
1073 SUBROUTINE qs_ot_get_derivative(matrix_hc, matrix_x, matrix_sx, matrix_gx, &
1074 qs_ot_env)
1075 TYPE(dbcsr_type), POINTER :: matrix_hc, matrix_x, matrix_sx, matrix_gx
1076 TYPE(qs_ot_type) :: qs_ot_env
1077
1078 CHARACTER(len=*), PARAMETER :: routinen = 'qs_ot_get_derivative'
1079 REAL(kind=dp), PARAMETER :: rone = 1.0_dp, rzero = 0.0_dp
1080
1081 INTEGER :: handle, k, n, ortho_k
1082 TYPE(dbcsr_type), POINTER :: matrix_hc_local, matrix_target
1083
1084 CALL timeset(routinen, handle)
1085
1086 NULLIFY (matrix_hc_local)
1087
1088 CALL dbcsr_get_info(matrix_x, nfullrows_total=n, nfullcols_total=k)
1089
1090 ! could in principle be taken inside qs_ot_get_derivative_* for increased efficiency
1091 ! create a local rotated version of matrix_hc leaving matrix_hc untouched (needed
1092 ! for lagrangian multipliers)
1093 IF (qs_ot_env%settings%do_rotation) THEN
1094 CALL dbcsr_copy(matrix_gx, matrix_hc) ! use gx as temporary
1095 CALL dbcsr_init_p(matrix_hc_local)
1096 CALL dbcsr_copy(matrix_hc_local, matrix_hc, name='matrix_hc_local')
1097 CALL dbcsr_set(matrix_hc_local, 0.0_dp)
1098 CALL dbcsr_multiply('N', 'T', rone, matrix_gx, qs_ot_env%rot_mat_u, rzero, matrix_hc_local)
1099 ELSE
1100 matrix_hc_local => matrix_hc
1101 END IF
1102
1103 IF (qs_ot_env%do_taylor) THEN
1104 CALL qs_ot_get_derivative_taylor(matrix_hc_local, matrix_x, matrix_sx, matrix_gx, qs_ot_env)
1105 ELSE
1106 CALL qs_ot_get_derivative_diag(matrix_hc_local, matrix_x, matrix_sx, matrix_gx, qs_ot_env)
1107 END IF
1108
1109 ! and make it orthogonal
1110 CALL dbcsr_get_info(qs_ot_env%matrix_sc0, nfullcols_total=ortho_k)
1111
1112 IF (ASSOCIATED(qs_ot_env%preconditioner)) THEN
1113 matrix_target => qs_ot_env%matrix_psc0
1114 ELSE
1115 matrix_target => qs_ot_env%matrix_sc0
1116 END IF
1117 ! first make the matrix os if not yet valid
1118 IF (.NOT. qs_ot_env%os_valid) THEN
1119 ! this assumes that the preconditioner is a single matrix
1120 ! that maps sc0 onto psc0
1121
1122 IF (ASSOCIATED(qs_ot_env%preconditioner)) THEN
1123 CALL apply_preconditioner(qs_ot_env%preconditioner, qs_ot_env%matrix_sc0, &
1124 qs_ot_env%matrix_psc0)
1125 END IF
1126 CALL dbcsr_multiply('T', 'N', rone, &
1127 qs_ot_env%matrix_sc0, matrix_target, &
1128 rzero, qs_ot_env%matrix_os)
1129 CALL cp_dbcsr_cholesky_decompose(qs_ot_env%matrix_os, &
1130 para_env=qs_ot_env%para_env, blacs_env=qs_ot_env%blacs_env)
1131 CALL cp_dbcsr_cholesky_invert(qs_ot_env%matrix_os, &
1132 para_env=qs_ot_env%para_env, blacs_env=qs_ot_env%blacs_env, &
1133 uplo_to_full=.true.)
1134 qs_ot_env%os_valid = .true.
1135 END IF
1136 CALL dbcsr_multiply('T', 'N', rone, matrix_target, matrix_gx, &
1137 rzero, qs_ot_env%matrix_buf1_ortho)
1138 CALL dbcsr_multiply('N', 'N', rone, qs_ot_env%matrix_os, &
1139 qs_ot_env%matrix_buf1_ortho, rzero, qs_ot_env%matrix_buf2_ortho)
1140 CALL dbcsr_multiply('N', 'N', -rone, qs_ot_env%matrix_sc0, &
1141 qs_ot_env%matrix_buf2_ortho, rone, matrix_gx)
1142 ! also treat the rot_mat gradient here
1143 IF (qs_ot_env%settings%do_rotation) THEN
1144 CALL qs_ot_rot_mat_derivative(qs_ot_env)
1145 END IF
1146
1147 IF (qs_ot_env%settings%do_rotation) THEN
1148 CALL dbcsr_release_p(matrix_hc_local)
1149 END IF
1150
1151 CALL timestop(handle)
1152
1153 END SUBROUTINE qs_ot_get_derivative
1154
1155! **************************************************************************************************
1156!> \brief ...
1157!> \param matrix_hc ...
1158!> \param matrix_x ...
1159!> \param matrix_sx ...
1160!> \param matrix_gx ...
1161!> \param qs_ot_env ...
1162! **************************************************************************************************
1163 SUBROUTINE qs_ot_get_derivative_diag(matrix_hc, matrix_x, matrix_sx, matrix_gx, &
1164 qs_ot_env)
1165
1166 TYPE(dbcsr_type), POINTER :: matrix_hc, matrix_x, matrix_sx, matrix_gx
1167 TYPE(qs_ot_type) :: qs_ot_env
1168
1169 CHARACTER(len=*), PARAMETER :: routinen = 'qs_ot_get_derivative_diag'
1170 REAL(kind=dp), PARAMETER :: rone = 1.0_dp, rzero = 0.0_dp
1171
1172 INTEGER :: handle, k, n
1173 TYPE(dbcsr_distribution_type) :: dist
1174
1175 CALL timeset(routinen, handle)
1176
1177 CALL dbcsr_get_info(matrix_x, nfullrows_total=n, nfullcols_total=k)
1178
1179 ! go for the derivative now
1180 ! this de/dc*(dX/dx)*sinp
1181 CALL dbcsr_multiply('N', 'N', rone, matrix_hc, qs_ot_env%matrix_sinp, rzero, matrix_gx)
1182 ! overlap hc*x
1183 CALL dbcsr_multiply('T', 'N', rone, matrix_hc, matrix_x, rzero, qs_ot_env%matrix_buf2)
1184 ! get it in the basis of the eigenvectors
1185 CALL dbcsr_multiply('N', 'N', rone, qs_ot_env%matrix_buf2, qs_ot_env%matrix_r, &
1186 rzero, qs_ot_env%matrix_buf1)
1187 CALL dbcsr_multiply('T', 'N', rone, qs_ot_env%matrix_r, qs_ot_env%matrix_buf1, &
1188 rzero, qs_ot_env%matrix_buf2)
1189
1190 ! get the schur product of O_uv*B_uv
1191 CALL dbcsr_hadamard_product(qs_ot_env%matrix_buf2, qs_ot_env%matrix_sinp_b, &
1192 qs_ot_env%matrix_buf3)
1193
1194 ! overlap hc*c0
1195 CALL dbcsr_multiply('T', 'N', rone, matrix_hc, qs_ot_env%matrix_c0, rzero, &
1196 qs_ot_env%matrix_buf2)
1197 ! get it in the basis of the eigenvectors
1198 CALL dbcsr_multiply('N', 'N', rone, qs_ot_env%matrix_buf2, qs_ot_env%matrix_r, &
1199 rzero, qs_ot_env%matrix_buf1)
1200 CALL dbcsr_multiply('T', 'N', rone, qs_ot_env%matrix_r, qs_ot_env%matrix_buf1, &
1201 rzero, qs_ot_env%matrix_buf2)
1202
1203 CALL dbcsr_hadamard_product(qs_ot_env%matrix_buf2, qs_ot_env%matrix_cosp_b, &
1204 qs_ot_env%matrix_buf4)
1205
1206 ! add the two bs and compute b+b^T
1207 CALL dbcsr_add(qs_ot_env%matrix_buf3, qs_ot_env%matrix_buf4, &
1208 alpha_scalar=rone, beta_scalar=rone)
1209
1210 ! get the b in the eigenvector basis
1211 CALL dbcsr_multiply('N', 'T', rone, qs_ot_env%matrix_buf3, qs_ot_env%matrix_r, &
1212 rzero, qs_ot_env%matrix_buf1)
1213 CALL dbcsr_multiply('N', 'N', rone, qs_ot_env%matrix_r, qs_ot_env%matrix_buf1, &
1214 rzero, qs_ot_env%matrix_buf3)
1215 CALL dbcsr_get_info(qs_ot_env%matrix_buf3, distribution=dist)
1216 CALL dbcsr_transposed(qs_ot_env%matrix_buf1, qs_ot_env%matrix_buf3, &
1217 shallow_data_copy=.false., use_distribution=dist, &
1218 transpose_distribution=.false.)
1219 CALL dbcsr_add(qs_ot_env%matrix_buf3, qs_ot_env%matrix_buf1, &
1220 alpha_scalar=rone, beta_scalar=rone)
1221
1222 ! and add to the derivative
1223 CALL dbcsr_multiply('N', 'N', rone, matrix_sx, qs_ot_env%matrix_buf3, &
1224 rone, matrix_gx)
1225 CALL timestop(handle)
1226
1227 END SUBROUTINE qs_ot_get_derivative_diag
1228
1229! **************************************************************************************************
1230!> \brief compute the derivative of the taylor expansion below
1231!> \param matrix_hc ...
1232!> \param matrix_x ...
1233!> \param matrix_sx ...
1234!> \param matrix_gx ...
1235!> \param qs_ot_env ...
1236! **************************************************************************************************
1237 SUBROUTINE qs_ot_get_derivative_taylor(matrix_hc, matrix_x, matrix_sx, matrix_gx, &
1238 qs_ot_env)
1239
1240 TYPE(dbcsr_type), POINTER :: matrix_hc, matrix_x, matrix_sx, matrix_gx
1241 TYPE(qs_ot_type) :: qs_ot_env
1242
1243 CHARACTER(len=*), PARAMETER :: routinen = 'qs_ot_get_derivative_taylor'
1244 REAL(kind=dp), PARAMETER :: rone = 1.0_dp, rzero = 0.0_dp
1245
1246 INTEGER :: handle, i, k, n
1247 REAL(kind=dp) :: cosfactor, sinfactor
1248 TYPE(dbcsr_distribution_type) :: dist
1249 TYPE(dbcsr_type), POINTER :: matrix_left, matrix_right
1250
1251 CALL timeset(routinen, handle)
1252
1253 CALL dbcsr_get_info(matrix_x, nfullrows_total=n, nfullcols_total=k)
1254
1255 ! go for the derivative now
1256 ! this de/dc*(dX/dx)*sinp i.e. zeroth order
1257 CALL dbcsr_multiply('N', 'N', rone, matrix_hc, qs_ot_env%matrix_sinp, rzero, matrix_gx)
1258
1259 IF (qs_ot_env%taylor_order <= 0) THEN
1260 CALL timestop(handle)
1261 RETURN
1262 END IF
1263
1264 ! we store the matrix that will multiply sx in matrix_r
1265 CALL dbcsr_set(qs_ot_env%matrix_r, rzero)
1266
1267 ! just better names for matrix_cosp_b and matrix_sinp_b (they are buffer space here)
1268 matrix_left => qs_ot_env%matrix_cosp_b
1269 matrix_right => qs_ot_env%matrix_sinp_b
1270
1271 ! overlap hc*x and add its transpose to matrix_left
1272 CALL dbcsr_multiply('T', 'N', rone, matrix_hc, matrix_x, rzero, matrix_left)
1273 CALL dbcsr_get_info(matrix_left, distribution=dist)
1274 CALL dbcsr_transposed(qs_ot_env%matrix_buf1, matrix_left, &
1275 shallow_data_copy=.false., use_distribution=dist, &
1276 transpose_distribution=.false.)
1277 CALL dbcsr_add(matrix_left, qs_ot_env%matrix_buf1, &
1278 alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
1279 CALL dbcsr_copy(matrix_right, matrix_left)
1280
1281 ! first order
1282 sinfactor = -1.0_dp/(2.0_dp*3.0_dp)
1283 CALL dbcsr_add(qs_ot_env%matrix_r, matrix_left, &
1284 alpha_scalar=1.0_dp, beta_scalar=sinfactor)
1285
1286 ! M
1287 ! OM+MO
1288 ! OOM+OMO+MOO
1289 ! ...
1290 DO i = 2, qs_ot_env%taylor_order
1291 sinfactor = sinfactor*(-1.0_dp)/real(2*i*(2*i + 1), kind=dp)
1292 CALL dbcsr_multiply('N', 'N', rone, qs_ot_env%matrix_p, matrix_left, rzero, qs_ot_env%matrix_buf1)
1293 CALL dbcsr_multiply('N', 'N', rone, matrix_right, qs_ot_env%matrix_p, rzero, matrix_left)
1294 CALL dbcsr_copy(matrix_right, matrix_left)
1295 CALL dbcsr_add(matrix_left, qs_ot_env%matrix_buf1, &
1296 1.0_dp, 1.0_dp)
1297 CALL dbcsr_add(qs_ot_env%matrix_r, matrix_left, &
1298 alpha_scalar=1.0_dp, beta_scalar=sinfactor)
1299 END DO
1300
1301 ! overlap hc*c0 and add its transpose to matrix_left
1302 CALL dbcsr_multiply('T', 'N', rone, matrix_hc, qs_ot_env%matrix_c0, rzero, matrix_left)
1303 CALL dbcsr_get_info(matrix_left, distribution=dist)
1304 CALL dbcsr_transposed(qs_ot_env%matrix_buf1, matrix_left, &
1305 shallow_data_copy=.false., use_distribution=dist, &
1306 transpose_distribution=.false.)
1307 CALL dbcsr_add(matrix_left, qs_ot_env%matrix_buf1, 1.0_dp, 1.0_dp)
1308 CALL dbcsr_copy(matrix_right, matrix_left)
1309
1310 ! first order
1311 cosfactor = -1.0_dp/(1.0_dp*2.0_dp)
1312 CALL dbcsr_add(qs_ot_env%matrix_r, matrix_left, &
1313 alpha_scalar=1.0_dp, beta_scalar=cosfactor)
1314
1315 ! M
1316 ! OM+MO
1317 ! OOM+OMO+MOO
1318 ! ...
1319 DO i = 2, qs_ot_env%taylor_order
1320 cosfactor = cosfactor*(-1.0_dp)/real(2*i*(2*i - 1), kind=dp)
1321 CALL dbcsr_multiply('N', 'N', rone, qs_ot_env%matrix_p, matrix_left, rzero, qs_ot_env%matrix_buf1)
1322 CALL dbcsr_multiply('N', 'N', rone, matrix_right, qs_ot_env%matrix_p, rzero, matrix_left)
1323 CALL dbcsr_copy(matrix_right, matrix_left)
1324 CALL dbcsr_add(matrix_left, qs_ot_env%matrix_buf1, 1.0_dp, 1.0_dp)
1325 CALL dbcsr_add(qs_ot_env%matrix_r, matrix_left, &
1326 alpha_scalar=1.0_dp, beta_scalar=cosfactor)
1327 END DO
1328
1329 ! and add to the derivative
1330 CALL dbcsr_multiply('N', 'N', rone, matrix_sx, qs_ot_env%matrix_r, rone, matrix_gx)
1331
1332 CALL timestop(handle)
1333
1334 END SUBROUTINE qs_ot_get_derivative_taylor
1335
1336! *************************************************************************************************
1337!> \brief computes a taylor expansion.
1338!> \param qs_ot_env ...
1339! **************************************************************************************************
1340 SUBROUTINE qs_ot_p2m_taylor(qs_ot_env)
1341 TYPE(qs_ot_type) :: qs_ot_env
1342
1343 CHARACTER(len=*), PARAMETER :: routinen = 'qs_ot_p2m_taylor'
1344 REAL(kind=dp), PARAMETER :: rone = 1.0_dp, rzero = 0.0_dp
1345
1346 INTEGER :: handle, i, k
1347 REAL(kind=dp) :: cosfactor, sinfactor
1348
1349 CALL timeset(routinen, handle)
1350
1351 ! zeroth order
1352 CALL dbcsr_set(qs_ot_env%matrix_cosp, rzero)
1353 CALL dbcsr_set(qs_ot_env%matrix_sinp, rzero)
1354 CALL dbcsr_add_on_diag(qs_ot_env%matrix_cosp, rone)
1355 CALL dbcsr_add_on_diag(qs_ot_env%matrix_sinp, rone)
1356
1357 IF (qs_ot_env%taylor_order <= 0) THEN
1358 CALL timestop(handle)
1359 RETURN
1360 END IF
1361
1362 ! first order
1363 cosfactor = -1.0_dp/(1.0_dp*2.0_dp)
1364 sinfactor = -1.0_dp/(2.0_dp*3.0_dp)
1365 CALL dbcsr_add(qs_ot_env%matrix_cosp, qs_ot_env%matrix_p, alpha_scalar=1.0_dp, beta_scalar=cosfactor)
1366 CALL dbcsr_add(qs_ot_env%matrix_sinp, qs_ot_env%matrix_p, alpha_scalar=1.0_dp, beta_scalar=sinfactor)
1367 IF (qs_ot_env%taylor_order <= 1) THEN
1368 CALL timestop(handle)
1369 RETURN
1370 END IF
1371
1372 ! other orders
1373 CALL dbcsr_get_info(qs_ot_env%matrix_p, nfullrows_total=k)
1374 CALL dbcsr_copy(qs_ot_env%matrix_r, qs_ot_env%matrix_p)
1375
1376 DO i = 2, qs_ot_env%taylor_order
1377 ! new power of p
1378 CALL dbcsr_multiply('N', 'N', rone, qs_ot_env%matrix_p, qs_ot_env%matrix_r, &
1379 rzero, qs_ot_env%matrix_buf1)
1380 CALL dbcsr_copy(qs_ot_env%matrix_r, qs_ot_env%matrix_buf1)
1381 ! add to the taylor expansion so far
1382 cosfactor = cosfactor*(-1.0_dp)/real(2*i*(2*i - 1), kind=dp)
1383 sinfactor = sinfactor*(-1.0_dp)/real(2*i*(2*i + 1), kind=dp)
1384 CALL dbcsr_add(qs_ot_env%matrix_cosp, qs_ot_env%matrix_r, &
1385 alpha_scalar=1.0_dp, beta_scalar=cosfactor)
1386 CALL dbcsr_add(qs_ot_env%matrix_sinp, qs_ot_env%matrix_r, &
1387 alpha_scalar=1.0_dp, beta_scalar=sinfactor)
1388 END DO
1389
1390 CALL timestop(handle)
1391
1392 END SUBROUTINE qs_ot_p2m_taylor
1393
1394! **************************************************************************************************
1395!> \brief given p, computes - eigenstuff (matrix_r,evals)
1396!> - cos(p^0.5),p^(-0.5)*sin(p^0.5)
1397!> - the real b matrices, needed for the derivatives of these guys
1398!> cosp_b_ij=(1/(2pii) * int(cos(z^1/2)/((z-eval(i))*(z-eval(j))))
1399!> sinp_b_ij=(1/(2pii) * int(z^(-1/2)*sin(z^1/2)/((z-eval(i))*(z-eval(j))))
1400!> \param qs_ot_env ...
1401! **************************************************************************************************
1402 SUBROUTINE qs_ot_p2m_diag(qs_ot_env)
1403
1404 TYPE(qs_ot_type) :: qs_ot_env
1405
1406 CHARACTER(len=*), PARAMETER :: routinen = 'qs_ot_p2m_diag'
1407 REAL(kind=dp), PARAMETER :: rone = 1.0_dp, rzero = 0.0_dp
1408
1409 INTEGER :: col, col_offset, col_size, handle, i, j, &
1410 k, row, row_offset, row_size
1411 REAL(dp), DIMENSION(:, :), POINTER :: block
1412 REAL(kind=dp) :: a, b
1413 TYPE(dbcsr_iterator_type) :: iter
1414
1415 CALL timeset(routinen, handle)
1416
1417 CALL dbcsr_get_info(qs_ot_env%matrix_p, nfullrows_total=k)
1418 CALL dbcsr_copy(qs_ot_env%matrix_buf1, qs_ot_env%matrix_p)
1419 CALL cp_dbcsr_syevd(qs_ot_env%matrix_buf1, qs_ot_env%matrix_r, qs_ot_env%evals, &
1420 qs_ot_env%para_env, qs_ot_env%blacs_env)
1421 DO i = 1, k
1422 qs_ot_env%evals(i) = max(0.0_dp, qs_ot_env%evals(i))
1423 END DO
1424
1425 !$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(i) SHARED(k,qs_ot_env)
1426 DO i = 1, k
1427 qs_ot_env%dum(i) = cos(sqrt(qs_ot_env%evals(i)))
1428 END DO
1429 CALL dbcsr_copy(qs_ot_env%matrix_buf1, qs_ot_env%matrix_r)
1430 CALL dbcsr_scale_by_vector(qs_ot_env%matrix_buf1, alpha=qs_ot_env%dum, side='right')
1431 CALL dbcsr_multiply('N', 'T', rone, qs_ot_env%matrix_r, qs_ot_env%matrix_buf1, &
1432 rzero, qs_ot_env%matrix_cosp)
1433
1434 !$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(i) SHARED(k,qs_ot_env)
1435 DO i = 1, k
1436 qs_ot_env%dum(i) = qs_ot_sinc(sqrt(qs_ot_env%evals(i)))
1437 END DO
1438 CALL dbcsr_copy(qs_ot_env%matrix_buf1, qs_ot_env%matrix_r)
1439 CALL dbcsr_scale_by_vector(qs_ot_env%matrix_buf1, alpha=qs_ot_env%dum, side='right')
1440 CALL dbcsr_multiply('N', 'T', rone, qs_ot_env%matrix_r, qs_ot_env%matrix_buf1, &
1441 rzero, qs_ot_env%matrix_sinp)
1442
1443 CALL dbcsr_copy(qs_ot_env%matrix_cosp_b, qs_ot_env%matrix_cosp)
1444 CALL dbcsr_iterator_start(iter, qs_ot_env%matrix_cosp_b)
1445 DO WHILE (dbcsr_iterator_blocks_left(iter))
1446 CALL dbcsr_iterator_next_block(iter, row, col, block, &
1447 row_size=row_size, col_size=col_size, &
1448 row_offset=row_offset, col_offset=col_offset)
1449 DO j = 1, col_size
1450 DO i = 1, row_size
1451 a = (sqrt(qs_ot_env%evals(row_offset + i - 1)) &
1452 - sqrt(qs_ot_env%evals(col_offset + j - 1)))/2.0_dp
1453 b = (sqrt(qs_ot_env%evals(row_offset + i - 1)) &
1454 + sqrt(qs_ot_env%evals(col_offset + j - 1)))/2.0_dp
1455 block(i, j) = -0.5_dp*qs_ot_sinc(a)*qs_ot_sinc(b)
1456 END DO
1457 END DO
1458 END DO
1459 CALL dbcsr_iterator_stop(iter)
1460
1461 CALL dbcsr_copy(qs_ot_env%matrix_sinp_b, qs_ot_env%matrix_sinp)
1462 CALL dbcsr_iterator_start(iter, qs_ot_env%matrix_sinp_b)
1463 DO WHILE (dbcsr_iterator_blocks_left(iter))
1464 CALL dbcsr_iterator_next_block(iter, row, col, block, &
1465 row_size=row_size, col_size=col_size, &
1466 row_offset=row_offset, col_offset=col_offset)
1467 DO j = 1, col_size
1468 DO i = 1, row_size
1469 a = sqrt(qs_ot_env%evals(row_offset + i - 1))
1470 b = sqrt(qs_ot_env%evals(col_offset + j - 1))
1471 block(i, j) = qs_ot_sincf(a, b)
1472 END DO
1473 END DO
1474 END DO
1475 CALL dbcsr_iterator_stop(iter)
1476
1477 CALL timestop(handle)
1478
1479 END SUBROUTINE qs_ot_p2m_diag
1480
1481! **************************************************************************************************
1482!> \brief computes sin(x)/x for all values of the argument
1483!> \param x ...
1484!> \return ...
1485! **************************************************************************************************
1486 FUNCTION qs_ot_sinc(x)
1487
1488 REAL(kind=dp), INTENT(IN) :: x
1489 REAL(kind=dp) :: qs_ot_sinc
1490
1491 REAL(kind=dp), PARAMETER :: q1 = 1.0_dp, q2 = -q1/(2.0_dp*3.0_dp), q3 = -q2/(4.0_dp*5.0_dp), &
1492 q4 = -q3/(6.0_dp*7.0_dp), q5 = -q4/(8.0_dp*9.0_dp), q6 = -q5/(10.0_dp*11.0_dp), &
1493 q7 = -q6/(12.0_dp*13.0_dp), q8 = -q7/(14.0_dp*15.0_dp), q9 = -q8/(16.0_dp*17.0_dp), &
1494 q10 = -q9/(18.0_dp*19.0_dp)
1495
1496 REAL(kind=dp) :: y
1497
1498 IF (abs(x) > 0.5_dp) THEN
1499 qs_ot_sinc = sin(x)/x
1500 ELSE
1501 y = x*x
1502 qs_ot_sinc = q1 + y*(q2 + y*(q3 + y*(q4 + y*(q5 + y*(q6 + y*(q7 + y*(q8 + y*(q9 + y*(q10)))))))))
1503 END IF
1504 END FUNCTION qs_ot_sinc
1505
1506! **************************************************************************************************
1507!> \brief computes (1/(x^2-y^2))*(sinc(x)-sinc(y)) for all positive values of the arguments
1508!> \param xa ...
1509!> \param ya ...
1510!> \return ...
1511! **************************************************************************************************
1512 FUNCTION qs_ot_sincf(xa, ya)
1513
1514 REAL(kind=dp), INTENT(IN) :: xa, ya
1515 REAL(kind=dp) :: qs_ot_sincf
1516
1517 INTEGER :: i
1518 REAL(kind=dp) :: a, b, rs, sf, x, xs, y, ybx, ybxs
1519
1520 ! this is currently a limit of the routine, could be removed rather easily
1521 IF (xa < 0) cpabort("x is negative")
1522 IF (ya < 0) cpabort("y is negative")
1523
1524 IF (xa < ya) THEN
1525 x = ya
1526 y = xa
1527 ELSE
1528 x = xa
1529 y = ya
1530 END IF
1531
1532 IF (x < 0.5_dp) THEN ! use series, keeping in mind that x,y,x+y,x-y can all be zero
1533
1534 qs_ot_sincf = 0.0_dp
1535 IF (x > 0.0_dp) THEN
1536 ybx = y/x
1537 ELSE ! should be irrelevant !?
1538 ybx = 0.0_dp
1539 END IF
1540
1541 sf = -1.0_dp/((1.0_dp + ybx)*6.0_dp)
1542 rs = 1.0_dp
1543 ybxs = ybx
1544 xs = 1.0_dp
1545
1546 DO i = 1, 10
1547 qs_ot_sincf = qs_ot_sincf + sf*rs*xs*(1.0_dp + ybxs)
1548 sf = -sf/(real((2*i + 2), dp)*real((2*i + 3), dp))
1549 rs = rs + ybxs
1550 ybxs = ybxs*ybx
1551 xs = xs*x*x
1552 END DO
1553
1554 ELSE ! no series expansion
1555 IF (x - y > 0.1_dp) THEN ! safe to use the normal form
1556 qs_ot_sincf = (qs_ot_sinc(x) - qs_ot_sinc(y))/((x + y)*(x - y))
1557 ELSE
1558 a = (x + y)/2.0_dp
1559 b = (x - y)/2.0_dp ! might be close to zero
1560 ! y (=(a-b)) can not be close to zero since it is close to x>0.5
1561 qs_ot_sincf = (qs_ot_sinc(b)*cos(a) - qs_ot_sinc(a)*cos(b))/(2*x*y)
1562 END IF
1563 END IF
1564
1565 END FUNCTION qs_ot_sincf
1566
1567END MODULE qs_ot
preconditioner::apply_preconditioner
Definition preconditioner.F:68
arnoldi_api
arnoldi iteration using dbcsr
Definition arnoldi_api.F:16
arnoldi_api::arnoldi_extremal
subroutine, public arnoldi_extremal(matrix_a, max_ev, min_ev, converged, threshold, max_iter)
simple wrapper to estimate extremal eigenvalues with arnoldi, using the old lanczos interface this hi...
Definition arnoldi_api.F:254
cp_dbcsr_api
Definition cp_dbcsr_api.F:8
cp_dbcsr_api::dbcsr_transposed
subroutine, public dbcsr_transposed(transposed, normal, shallow_data_copy, transpose_distribution, use_distribution)
...
Definition cp_dbcsr_api.F:1242
cp_dbcsr_api::dbcsr_release_p
subroutine, public dbcsr_release_p(matrix)
...
Definition cp_dbcsr_api.F:222
cp_dbcsr_api::dbcsr_scale
subroutine, public dbcsr_scale(matrix, alpha_scalar)
...
Definition cp_dbcsr_api.F:1188
cp_dbcsr_api::dbcsr_iterator_blocks_left
logical function, public dbcsr_iterator_blocks_left(iterator)
...
Definition cp_dbcsr_api.F:944
cp_dbcsr_api::dbcsr_iterator_stop
subroutine, public dbcsr_iterator_stop(iterator)
...
Definition cp_dbcsr_api.F:1050
cp_dbcsr_api::dbcsr_copy
subroutine, public dbcsr_copy(matrix_b, matrix_a, name, keep_sparsity, keep_imaginary)
...
Definition cp_dbcsr_api.F:371
cp_dbcsr_api::dbcsr_get_block_p
subroutine, public dbcsr_get_block_p(matrix, row, col, block, found, row_size, col_size)
...
Definition cp_dbcsr_api.F:693
cp_dbcsr_api::dbcsr_multiply
subroutine, public dbcsr_multiply(transa, transb, alpha, matrix_a, matrix_b, beta, matrix_c, first_row, last_row, first_column, last_column, first_k, last_k, retain_sparsity, filter_eps, flop)
...
Definition cp_dbcsr_api.F:1096
cp_dbcsr_api::dbcsr_get_info
subroutine, public dbcsr_get_info(matrix, nblkrows_total, nblkcols_total, nfullrows_total, nfullcols_total, nblkrows_local, nblkcols_local, nfullrows_local, nfullcols_local, my_prow, my_pcol, local_rows, local_cols, proc_row_dist, proc_col_dist, row_blk_size, col_blk_size, row_blk_offset, col_blk_offset, distribution, name, matrix_type, group)
...
Definition cp_dbcsr_api.F:808
cp_dbcsr_api::dbcsr_init_p
subroutine, public dbcsr_init_p(matrix)
...
Definition cp_dbcsr_api.F:207
cp_dbcsr_api::dbcsr_iterator_next_block
subroutine, public dbcsr_iterator_next_block(iterator, row, column, block, block_number_argument_has_been_removed, row_size, col_size, row_offset, col_offset, transposed)
...
Definition cp_dbcsr_api.F:971
cp_dbcsr_api::dbcsr_filter
subroutine, public dbcsr_filter(matrix, eps)
...
Definition cp_dbcsr_api.F:658
cp_dbcsr_api::dbcsr_get_occupation
real(kind=dp) function, public dbcsr_get_occupation(matrix)
...
Definition cp_dbcsr_api.F:893
cp_dbcsr_api::dbcsr_iterator_start
subroutine, public dbcsr_iterator_start(iterator, matrix, shared, dynamic, dynamic_byrows)
...
Definition cp_dbcsr_api.F:1012
cp_dbcsr_api::dbcsr_set
subroutine, public dbcsr_set(matrix, alpha)
...
Definition cp_dbcsr_api.F:1204
cp_dbcsr_api::dbcsr_release
subroutine, public dbcsr_release(matrix)
...
Definition cp_dbcsr_api.F:1142
cp_dbcsr_api::dbcsr_add
subroutine, public dbcsr_add(matrix_a, matrix_b, alpha_scalar, beta_scalar)
...
Definition cp_dbcsr_api.F:254
cp_dbcsr_cholesky
Interface to (sca)lapack for the Cholesky based procedures.
Definition cp_dbcsr_cholesky.F:17
cp_dbcsr_cholesky::cp_dbcsr_cholesky_decompose
subroutine, public cp_dbcsr_cholesky_decompose(matrix, n, para_env, blacs_env)
used to replace a symmetric positive def. matrix M with its cholesky decomposition U: M = U^T * U,...
Definition cp_dbcsr_cholesky.F:61
cp_dbcsr_cholesky::cp_dbcsr_cholesky_restore
subroutine, public cp_dbcsr_cholesky_restore(matrix, neig, matrixb, matrixout, op, pos, transa, para_env, blacs_env)
...
Definition cp_dbcsr_cholesky.F:176
cp_dbcsr_cholesky::cp_dbcsr_cholesky_invert
subroutine, public cp_dbcsr_cholesky_invert(matrix, n, para_env, blacs_env, uplo_to_full)
used to replace the cholesky decomposition by the inverse
Definition cp_dbcsr_cholesky.F:114
cp_dbcsr_contrib
Definition cp_dbcsr_contrib.F:8
cp_dbcsr_contrib::dbcsr_gershgorin_norm
real(dp) function, public dbcsr_gershgorin_norm(matrix)
Compute the gershgorin norm of a dbcsr matrix.
Definition cp_dbcsr_contrib.F:167
cp_dbcsr_contrib::dbcsr_add_on_diag
subroutine, public dbcsr_add_on_diag(matrix, alpha)
Adds the given scalar to the diagonal of the matrix. Reserves any missing diagonal blocks.
Definition cp_dbcsr_contrib.F:331
cp_dbcsr_contrib::dbcsr_frobenius_norm
real(dp) function, public dbcsr_frobenius_norm(matrix)
Compute the frobenius norm of a dbcsr matrix.
Definition cp_dbcsr_contrib.F:127
cp_dbcsr_contrib::dbcsr_hadamard_product
subroutine, public dbcsr_hadamard_product(matrix_a, matrix_b, matrix_c)
Hadamard product: C = A . B (C needs to be different from A and B)
Definition cp_dbcsr_contrib.F:53
cp_dbcsr_contrib::dbcsr_scale_by_vector
subroutine, public dbcsr_scale_by_vector(matrix, alpha, side)
Scales the rows/columns of given matrix.
Definition cp_dbcsr_contrib.F:486
cp_dbcsr_diag
Interface to (sca)lapack for the Cholesky based procedures.
Definition cp_dbcsr_diag.F:17
cp_dbcsr_diag::cp_dbcsr_heevd
subroutine, public cp_dbcsr_heevd(matrix_re, matrix_im, eigenvectors_re, eigenvectors_im, eigenvalues, para_env, blacs_env)
...
Definition cp_dbcsr_diag.F:186
cp_dbcsr_diag::cp_dbcsr_syevd
subroutine, public cp_dbcsr_syevd(matrix, eigenvectors, eigenvalues, para_env, blacs_env)
...
Definition cp_dbcsr_diag.F:67
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:23
preconditioner_types
types of preconditioners
Definition preconditioner_types.F:14
preconditioner
computes preconditioners, and implements methods to apply them currently used in qs_ot
Definition preconditioner.F:15
qs_ot_types
orbital transformations
Definition qs_ot_types.F:15
qs_ot
orbital transformations
Definition qs_ot.F:15
qs_ot::qs_ot_get_p
subroutine, public qs_ot_get_p(matrix_x, matrix_sx, qs_ot_env)
computes p=x*S*x and the matrix functionals related matrices
Definition qs_ot.F:749
qs_ot::qs_ot_get_derivative_ref
subroutine, public qs_ot_get_derivative_ref(matrix_hc, matrix_x, matrix_sx, matrix_gx, qs_ot_env)
...
Definition qs_ot.F:702
qs_ot::qs_ot_get_orbitals
subroutine, public qs_ot_get_orbitals(matrix_c, matrix_x, qs_ot_env)
c=(c0*cos(p^0.5)+x*sin(p^0.5)*p^(-0.5)) x rot_mat_u this assumes that x is already ortho to S*C0,...
Definition qs_ot.F:1021
qs_ot::qs_ot_get_orbitals_ref
subroutine, public qs_ot_get_orbitals_ref(matrix_c, matrix_s, matrix_x, matrix_sx, matrix_gx_old, matrix_dx, qs_ot_env, qs_ot_env1)
...
Definition qs_ot.F:517
qs_ot::qs_ot_get_derivative
subroutine, public qs_ot_get_derivative(matrix_hc, matrix_x, matrix_sx, matrix_gx, qs_ot_env)
this routines computes dE/dx=dx, with dx ortho to sc0 needs dE/dC=hc,C0,X,SX,p if preconditioned it w...
Definition qs_ot.F:1075
qs_ot::qs_ot_new_preconditioner
subroutine, public qs_ot_new_preconditioner(qs_ot_env, preconditioner)
gets ready to use the preconditioner/ or renew the preconditioner only keeps a pointer to the precond...
Definition qs_ot.F:73
cp_dbcsr_api::dbcsr_distribution_type
Definition cp_dbcsr_api.F:182
cp_dbcsr_api::dbcsr_iterator_type
Definition cp_dbcsr_api.F:188
cp_dbcsr_api::dbcsr_type
Definition cp_dbcsr_api.F:176
message_passing::mp_comm_type
Definition message_passing.F:158
preconditioner_types::preconditioner_type
Definition preconditioner_types.F:42
qs_ot_types::qs_ot_type
Definition qs_ot_types.F:101