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