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