(git:374b731)
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-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! **************************************************************************************************
15MODULE 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
57CONTAINS
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
1571END MODULE qs_ot
preconditioner::apply_preconditioner
Definition preconditioner.F:68
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
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