(git:44e3845)
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 < 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) &
299 cpabort("P not positive definite")
300 IF (eig(i) < 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 < 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 > 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 == 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 > 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 > 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 > 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 < qs_ot_env%settings%eps_irac_quick_exit) quick_exit = .true.
433 !
434 ! are we done?
435 IF (norm < 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 > 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 > 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 == 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 < 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 > 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 > 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 > 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 > 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 > 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 > 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 > 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 > 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 > 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) > 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 <= 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 <= 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 <= 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 < 0) cpabort("x is negative")
1521 IF (ya < 0) cpabort("y is negative")
1522
1523 IF (xa < ya) THEN
1524 x = ya
1525 y = xa
1526 ELSE
1527 x = xa
1528 y = ya
1529 END IF
1530
1531 IF (x < 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 > 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 > 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:1232
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:1178
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:969
cp_dbcsr_api::dbcsr_iterator_blocks_left
logical function, public dbcsr_iterator_blocks_left(iterator)
...
Definition cp_dbcsr_api.F:943
cp_dbcsr_api::dbcsr_iterator_stop
subroutine, public dbcsr_iterator_stop(iterator)
...
Definition cp_dbcsr_api.F:1040
cp_dbcsr_api::dbcsr_copy
subroutine, public dbcsr_copy(matrix_b, matrix_a, name, keep_sparsity, keep_imaginary)
...
Definition cp_dbcsr_api.F:370
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:692
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:1086
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:807
cp_dbcsr_api::dbcsr_init_p
subroutine, public dbcsr_init_p(matrix)
...
Definition cp_dbcsr_api.F:207
cp_dbcsr_api::dbcsr_filter
subroutine, public dbcsr_filter(matrix, eps)
...
Definition cp_dbcsr_api.F:657
cp_dbcsr_api::dbcsr_get_occupation
real(kind=dp) function, public dbcsr_get_occupation(matrix)
...
Definition cp_dbcsr_api.F:892
cp_dbcsr_api::dbcsr_iterator_start
subroutine, public dbcsr_iterator_start(iterator, matrix, shared, dynamic, dynamic_byrows)
...
Definition cp_dbcsr_api.F:1002
cp_dbcsr_api::dbcsr_set
subroutine, public dbcsr_set(matrix, alpha)
...
Definition cp_dbcsr_api.F:1194
cp_dbcsr_api::dbcsr_release
subroutine, public dbcsr_release(matrix)
...
Definition cp_dbcsr_api.F:1132
cp_dbcsr_api::dbcsr_add
subroutine, public dbcsr_add(matrix_a, matrix_b, alpha_scalar, beta_scalar)
...
Definition cp_dbcsr_api.F:253
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: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:189
preconditioner_types::preconditioner_type
Definition preconditioner_types.F:42
qs_ot_types::qs_ot_type
Definition qs_ot_types.F:101