Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief orbital transformations
10 : !> \par History
11 : !> Added Taylor expansion based computation of the matrix functions (01.2004)
12 : !> added additional rotation variables for non-equivalent occupied orbs (08.2004)
13 : !> \author Joost VandeVondele (06.2002)
14 : ! **************************************************************************************************
15 : MODULE qs_ot
16 : USE arnoldi_api, ONLY: arnoldi_extremal
17 : USE cp_dbcsr_api, ONLY: &
18 : dbcsr_add, dbcsr_copy, dbcsr_create, dbcsr_distribution_type, dbcsr_filter, &
19 : dbcsr_finalize, dbcsr_get_block_p, dbcsr_get_info, dbcsr_get_occupation, dbcsr_init_p, &
20 : dbcsr_iterator_blocks_left, dbcsr_iterator_next_block, dbcsr_iterator_start, &
21 : dbcsr_iterator_stop, dbcsr_iterator_type, dbcsr_multiply, dbcsr_release, dbcsr_release_p, &
22 : dbcsr_reserve_blocks, dbcsr_scale, dbcsr_set, dbcsr_transposed, dbcsr_type, &
23 : dbcsr_type_no_symmetry
24 : USE cp_dbcsr_cholesky, ONLY: cp_dbcsr_cholesky_decompose,&
25 : cp_dbcsr_cholesky_invert,&
26 : cp_dbcsr_cholesky_restore
27 : USE cp_dbcsr_contrib, ONLY: dbcsr_add_on_diag,&
28 : dbcsr_frobenius_norm,&
29 : dbcsr_gershgorin_norm,&
30 : dbcsr_hadamard_product,&
31 : dbcsr_scale_by_vector
32 : USE cp_dbcsr_diag, ONLY: cp_dbcsr_heevd,&
33 : cp_dbcsr_syevd
34 : USE kinds, ONLY: dp
35 : USE message_passing, ONLY: mp_comm_type
36 : USE preconditioner, ONLY: apply_preconditioner
37 : USE preconditioner_types, ONLY: preconditioner_type
38 : USE qs_ot_types, ONLY: qs_ot_type
39 : #include "./base/base_uses.f90"
40 :
41 : IMPLICIT NONE
42 : PRIVATE
43 :
44 : PUBLIC :: qs_ot_get_p
45 : PUBLIC :: qs_ot_get_orbitals
46 : PUBLIC :: qs_ot_get_derivative
47 : PUBLIC :: qs_ot_project_gradient
48 : PUBLIC :: qs_ot_get_orbitals_ref
49 : PUBLIC :: qs_ot_get_derivative_ref
50 : PUBLIC :: qs_ot_generate_rotation
51 : PUBLIC :: qs_ot_rot_mat_derivative
52 : PUBLIC :: qs_ot_new_preconditioner
53 : PUBLIC :: qs_ot_complex_multiply
54 : PUBLIC :: qs_ot_p2m_diag
55 : PUBLIC :: qs_ot_ref_update
56 : PUBLIC :: qs_ot_spectral_kernel
57 : PUBLIC :: qs_ot_spectral_values
58 : PUBLIC :: qs_ot_square_transpose
59 :
60 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_ot'
61 :
62 : CONTAINS
63 :
64 : ! **************************************************************************************************
65 : !> \brief gets ready to use the preconditioner/ or renew the preconditioner
66 : !> only keeps a pointer to the preconditioner.
67 : !> If you change the preconditioner, you have to call this routine
68 : !> you remain responsible of proper deallocate of your preconditioner
69 : !> (or you can reuse it on the next step of the computation)
70 : !> \param qs_ot_env ...
71 : !> \param preconditioner ...
72 : ! **************************************************************************************************
73 9515 : SUBROUTINE qs_ot_new_preconditioner(qs_ot_env, preconditioner)
74 : TYPE(qs_ot_type) :: qs_ot_env
75 : TYPE(preconditioner_type), POINTER :: preconditioner
76 :
77 : INTEGER :: ncoef
78 :
79 9515 : qs_ot_env%preconditioner => preconditioner
80 9515 : qs_ot_env%preconditioner_changed = .TRUE.
81 9515 : qs_ot_env%os_valid = .FALSE.
82 9515 : IF (.NOT. ASSOCIATED(qs_ot_env%matrix_psc0)) THEN
83 9515 : CALL dbcsr_init_p(qs_ot_env%matrix_psc0)
84 9515 : CALL dbcsr_copy(qs_ot_env%matrix_psc0, qs_ot_env%matrix_sc0, 'matrix_psc0')
85 : END IF
86 9515 : IF (qs_ot_env%has_complex_kpoint_state .AND. &
87 : .NOT. ASSOCIATED(qs_ot_env%matrix_psc0_im)) THEN
88 749 : CALL dbcsr_init_p(qs_ot_env%matrix_psc0_im)
89 749 : CALL dbcsr_copy(qs_ot_env%matrix_psc0_im, qs_ot_env%matrix_sc0_im, 'matrix_psc0_im')
90 : END IF
91 :
92 9515 : IF (.NOT. qs_ot_env%use_dx) THEN
93 5775 : qs_ot_env%use_dx = .TRUE.
94 5775 : CALL dbcsr_init_p(qs_ot_env%matrix_dx)
95 5775 : CALL dbcsr_copy(qs_ot_env%matrix_dx, qs_ot_env%matrix_gx, 'matrix_dx')
96 5775 : IF (qs_ot_env%has_complex_kpoint_state) THEN
97 444 : CALL dbcsr_init_p(qs_ot_env%matrix_dx_im)
98 444 : CALL dbcsr_copy(qs_ot_env%matrix_dx_im, qs_ot_env%matrix_gx_im, 'matrix_dx_im')
99 : END IF
100 5775 : IF (qs_ot_env%settings%do_rotation) THEN
101 74 : CALL dbcsr_init_p(qs_ot_env%rot_mat_dx)
102 74 : CALL dbcsr_copy(qs_ot_env%rot_mat_dx, qs_ot_env%rot_mat_gx, 'rot_mat_dx')
103 74 : IF (qs_ot_env%has_complex_kpoint_state) THEN
104 48 : CALL dbcsr_init_p(qs_ot_env%rot_mat_dx_im)
105 48 : CALL dbcsr_copy(qs_ot_env%rot_mat_dx_im, qs_ot_env%rot_mat_gx_im, 'rot_mat_dx_im')
106 : END IF
107 : END IF
108 5775 : IF (qs_ot_env%settings%do_ener) THEN
109 40 : ncoef = SIZE(qs_ot_env%ener_gx)
110 120 : ALLOCATE (qs_ot_env%ener_dx(ncoef))
111 324 : qs_ot_env%ener_dx = 0.0_dp
112 : END IF
113 : END IF
114 :
115 9515 : END SUBROUTINE qs_ot_new_preconditioner
116 :
117 : ! **************************************************************************************************
118 : !> \brief multiply paired real/imaginary DBCSR matrices, with optional conjugate transposes
119 : !> \param op_a N or C
120 : !> \param op_b N or C
121 : !> \param a_re real part of A
122 : !> \param a_im imaginary part of A
123 : !> \param b_re real part of B
124 : !> \param b_im imaginary part of B
125 : !> \param c_re real part of A*B
126 : !> \param c_im imaginary part of A*B
127 : !> \param tmp real workspace shaped like C
128 : ! **************************************************************************************************
129 42726 : SUBROUTINE qs_ot_complex_multiply(op_a, op_b, a_re, a_im, b_re, b_im, c_re, c_im, tmp)
130 : CHARACTER(LEN=1), INTENT(IN) :: op_a, op_b
131 : TYPE(dbcsr_type) :: a_re, a_im, b_re, b_im, c_re, c_im, tmp
132 :
133 : CHARACTER(LEN=1) :: db_op_a, db_op_b
134 : REAL(KIND=dp) :: sign_a, sign_b
135 :
136 42726 : IF (op_a /= 'N' .AND. op_a /= 'C') THEN
137 0 : CPABORT("Complex matrix product expects N or C for op_a")
138 : END IF
139 42726 : IF (op_b /= 'N' .AND. op_b /= 'C') THEN
140 0 : CPABORT("Complex matrix product expects N or C for op_b")
141 : END IF
142 42726 : db_op_a = MERGE('N', 'T', op_a == 'N')
143 42726 : db_op_b = MERGE('N', 'T', op_b == 'N')
144 42726 : sign_a = MERGE(1.0_dp, -1.0_dp, op_a == 'N')
145 42726 : sign_b = MERGE(1.0_dp, -1.0_dp, op_b == 'N')
146 :
147 42726 : CALL dbcsr_multiply(db_op_a, db_op_b, 1.0_dp, a_re, b_re, 0.0_dp, c_re)
148 42726 : CALL dbcsr_multiply(db_op_a, db_op_b, 1.0_dp, a_im, b_im, 0.0_dp, tmp)
149 42726 : CALL dbcsr_add(c_re, tmp, alpha_scalar=1.0_dp, beta_scalar=-sign_a*sign_b)
150 :
151 42726 : CALL dbcsr_multiply(db_op_a, db_op_b, sign_b, a_re, b_im, 0.0_dp, c_im)
152 42726 : CALL dbcsr_multiply(db_op_a, db_op_b, sign_a, a_im, b_re, 0.0_dp, tmp)
153 42726 : CALL dbcsr_add(c_im, tmp, alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
154 :
155 42726 : END SUBROUTINE qs_ot_complex_multiply
156 :
157 : ! **************************************************************************************************
158 : !> \brief ...
159 : !> \param qs_ot_env ...
160 : !> \param C_NEW ...
161 : !> \param SC ...
162 : !> \param G_OLD ...
163 : !> \param D ...
164 : ! **************************************************************************************************
165 840 : SUBROUTINE qs_ot_on_the_fly_localize(qs_ot_env, C_NEW, SC, G_OLD, D)
166 : !
167 : TYPE(qs_ot_type) :: qs_ot_env
168 : TYPE(dbcsr_type), POINTER :: C_NEW, SC, G_OLD, D
169 :
170 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_ot_on_the_fly_localize'
171 : INTEGER, PARAMETER :: taylor_order = 50
172 : REAL(KIND=dp), PARAMETER :: alpha = 0.1_dp, f2_eps = 0.01_dp
173 :
174 : INTEGER :: col, col_size, handle, i, k, n, p, row, &
175 : row_size
176 168 : REAL(dp), DIMENSION(:, :), POINTER :: block
177 : REAL(KIND=dp) :: expfactor, f2, norm_fro, norm_gct, tmp
178 : TYPE(dbcsr_distribution_type) :: dist
179 : TYPE(dbcsr_iterator_type) :: iter
180 : TYPE(dbcsr_type), POINTER :: C, Gp1, Gp2, GU, U
181 : TYPE(mp_comm_type) :: group
182 :
183 168 : CALL timeset(routineN, handle)
184 : !
185 : !
186 168 : CALL dbcsr_get_info(C_NEW, nfullrows_total=n, nfullcols_total=k)
187 : !
188 : ! C = C*expm(-G)
189 168 : GU => qs_ot_env%buf1_k_k_nosym ! a buffer
190 168 : U => qs_ot_env%buf2_k_k_nosym ! a buffer
191 168 : Gp1 => qs_ot_env%buf3_k_k_nosym ! a buffer
192 168 : Gp2 => qs_ot_env%buf4_k_k_nosym ! a buffer
193 168 : C => qs_ot_env%buf1_n_k ! a buffer
194 : !
195 : ! compute the derivative of the norm
196 : !-------------------------------------------------------------------
197 : ! (x^2+eps)^1/2
198 168 : f2 = 0.0_dp
199 168 : CALL dbcsr_copy(C, C_NEW)
200 168 : CALL dbcsr_iterator_start(iter, C)
201 392 : DO WHILE (dbcsr_iterator_blocks_left(iter))
202 224 : CALL dbcsr_iterator_next_block(iter, row, col, block, row_size=row_size, col_size=col_size)
203 1400 : DO p = 1, col_size ! p
204 9912 : DO i = 1, row_size ! i
205 8680 : tmp = SQRT(block(i, p)**2 + f2_eps)
206 8680 : f2 = f2 + tmp
207 9688 : block(i, p) = block(i, p)/tmp
208 : END DO
209 : END DO
210 : END DO
211 168 : CALL dbcsr_iterator_stop(iter)
212 168 : CALL dbcsr_get_info(C, group=group)
213 168 : CALL group%sum(f2)
214 : !
215 : !
216 168 : CALL dbcsr_multiply('T', 'N', 1.0_dp, C, C_NEW, 0.0_dp, GU)
217 : !
218 : ! antisymetrize
219 168 : CALL dbcsr_get_info(GU, distribution=dist)
220 : CALL dbcsr_transposed(U, GU, shallow_data_copy=.FALSE., &
221 : use_distribution=dist, &
222 168 : transpose_distribution=.FALSE.)
223 168 : CALL dbcsr_add(GU, U, alpha_scalar=-0.5_dp, beta_scalar=0.5_dp)
224 : !-------------------------------------------------------------------
225 : !
226 168 : norm_fro = dbcsr_frobenius_norm(GU)
227 168 : norm_gct = dbcsr_gershgorin_norm(GU)
228 : !write(*,*) 'qs_ot_localize: ||P-I||_f=',norm_fro,' ||P-I||_GCT=',norm_gct
229 : !
230 : !kscale = CEILING(LOG(MIN(norm_fro,norm_gct))/LOG(2.0_dp))
231 : !scale = LOG(MIN(norm_fro,norm_gct))/LOG(2.0_dp)
232 : !write(*,*) 'qs_ot_localize: scale=',scale,' kscale=',kscale
233 : !
234 : ! rescale for steepest descent
235 168 : CALL dbcsr_scale(GU, -alpha)
236 : !
237 : ! compute unitary transform
238 : ! zeroth and first order
239 168 : expfactor = 1.0_dp
240 168 : CALL dbcsr_copy(U, GU)
241 168 : CALL dbcsr_scale(U, expfactor)
242 168 : CALL dbcsr_add_on_diag(U, 1.0_dp)
243 : ! other orders
244 168 : CALL dbcsr_copy(Gp1, GU)
245 806 : DO i = 2, taylor_order
246 : ! new power of G
247 806 : CALL dbcsr_multiply('N', 'N', 1.0_dp, GU, Gp1, 0.0_dp, Gp2)
248 806 : CALL dbcsr_copy(Gp1, Gp2)
249 : ! add to the taylor expansion so far
250 806 : expfactor = expfactor/REAL(i, KIND=dp)
251 806 : CALL dbcsr_add(U, Gp1, alpha_scalar=1.0_dp, beta_scalar=expfactor)
252 806 : norm_fro = dbcsr_frobenius_norm(Gp1)
253 : !write(*,*) 'Taylor expansion i=',i,' norm(X^i)/i!=',norm_fro*expfactor
254 806 : IF (norm_fro*expfactor < 1.0E-10_dp) EXIT
255 : END DO
256 : !
257 : ! rotate MOs
258 168 : CALL dbcsr_multiply('N', 'N', 1.0_dp, C_NEW, U, 0.0_dp, C)
259 168 : CALL dbcsr_copy(C_NEW, C)
260 : !
261 : ! rotate SC
262 168 : CALL dbcsr_multiply('N', 'N', 1.0_dp, SC, U, 0.0_dp, C)
263 168 : CALL dbcsr_copy(SC, C)
264 : !
265 : ! rotate D_i
266 168 : CALL dbcsr_multiply('N', 'N', 1.0_dp, D, U, 0.0_dp, C)
267 168 : CALL dbcsr_copy(D, C)
268 : !
269 : ! rotate G_i-1
270 168 : IF (ASSOCIATED(G_OLD)) THEN
271 168 : CALL dbcsr_multiply('N', 'N', 1.0_dp, G_OLD, U, 0.0_dp, C)
272 168 : CALL dbcsr_copy(G_OLD, C)
273 : END IF
274 : !
275 168 : CALL timestop(handle)
276 168 : END SUBROUTINE qs_ot_on_the_fly_localize
277 :
278 : ! **************************************************************************************************
279 : !> \brief ...
280 : !> \param qs_ot_env ...
281 : !> \param C_OLD ...
282 : !> \param C_TMP ...
283 : !> \param C_NEW ...
284 : !> \param P ...
285 : !> \param SC ...
286 : !> \param update ...
287 : ! **************************************************************************************************
288 1736 : SUBROUTINE qs_ot_ref_chol(qs_ot_env, C_OLD, C_TMP, C_NEW, P, SC, update)
289 : !
290 : TYPE(qs_ot_type) :: qs_ot_env
291 : TYPE(dbcsr_type) :: C_OLD, C_TMP, C_NEW, P, SC
292 : LOGICAL, INTENT(IN) :: update
293 :
294 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_ot_ref_chol'
295 :
296 : INTEGER :: handle, k, n
297 :
298 868 : CALL timeset(routineN, handle)
299 : !
300 868 : CALL dbcsr_get_info(C_NEW, nfullrows_total=n, nfullcols_total=k)
301 : !
302 : ! P = U'*U
303 868 : CALL cp_dbcsr_cholesky_decompose(P, k, qs_ot_env%para_env, qs_ot_env%blacs_env)
304 : !
305 : ! C_NEW = C_OLD*inv(U)
306 : CALL cp_dbcsr_cholesky_restore(C_OLD, k, P, C_NEW, op="SOLVE", pos="RIGHT", &
307 868 : transa="N", para_env=qs_ot_env%para_env, blacs_env=qs_ot_env%blacs_env)
308 : !
309 : ! Update SC if needed
310 868 : IF (update) THEN
311 : CALL cp_dbcsr_cholesky_restore(SC, k, P, C_TMP, op="SOLVE", pos="RIGHT", &
312 472 : transa="N", para_env=qs_ot_env%para_env, blacs_env=qs_ot_env%blacs_env)
313 472 : CALL dbcsr_copy(SC, C_TMP)
314 : END IF
315 : !
316 868 : CALL timestop(handle)
317 868 : END SUBROUTINE qs_ot_ref_chol
318 :
319 : ! **************************************************************************************************
320 : !> \brief ...
321 : !> \param qs_ot_env ...
322 : !> \param C_OLD ...
323 : !> \param C_TMP ...
324 : !> \param C_NEW ...
325 : !> \param P ...
326 : !> \param SC ...
327 : !> \param update ...
328 : ! **************************************************************************************************
329 478 : SUBROUTINE qs_ot_ref_lwdn(qs_ot_env, C_OLD, C_TMP, C_NEW, P, SC, update)
330 : !
331 : TYPE(qs_ot_type) :: qs_ot_env
332 : TYPE(dbcsr_type) :: C_OLD, C_TMP, C_NEW, P, SC
333 : LOGICAL, INTENT(IN) :: update
334 :
335 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_ot_ref_lwdn'
336 :
337 : INTEGER :: handle, i, k, n
338 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: eig, fun
339 : TYPE(dbcsr_type), POINTER :: V, W
340 :
341 478 : CALL timeset(routineN, handle)
342 : !
343 478 : CALL dbcsr_get_info(C_NEW, nfullrows_total=n, nfullcols_total=k)
344 : !
345 478 : V => qs_ot_env%buf1_k_k_nosym ! a buffer
346 478 : W => qs_ot_env%buf2_k_k_nosym ! a buffer
347 1912 : ALLOCATE (eig(k), fun(k))
348 : !
349 478 : CALL cp_dbcsr_syevd(P, V, eig, qs_ot_env%para_env, qs_ot_env%blacs_env)
350 : !
351 : ! compute the P^(-1/2)
352 2646 : DO i = 1, k
353 2168 : IF (eig(i) <= 0.0_dp) THEN
354 0 : CPABORT("P not positive definite")
355 : END IF
356 2646 : IF (eig(i) < 1.0E-8_dp) THEN
357 0 : fun(i) = 0.0_dp
358 : ELSE
359 2168 : fun(i) = 1.0_dp/SQRT(eig(i))
360 : END IF
361 : END DO
362 478 : CALL dbcsr_copy(W, V)
363 478 : CALL dbcsr_scale_by_vector(V, alpha=fun, side='right')
364 478 : CALL dbcsr_multiply('N', 'T', 1.0_dp, W, V, 0.0_dp, P)
365 : !
366 : ! Update C
367 478 : CALL dbcsr_multiply('N', 'N', 1.0_dp, C_OLD, P, 0.0_dp, C_NEW)
368 : !
369 : ! Update SC if needed
370 478 : IF (update) THEN
371 300 : CALL dbcsr_multiply('N', 'N', 1.0_dp, SC, P, 0.0_dp, C_TMP)
372 300 : CALL dbcsr_copy(SC, C_TMP)
373 : END IF
374 : !
375 478 : DEALLOCATE (eig, fun)
376 : !
377 478 : CALL timestop(handle)
378 478 : END SUBROUTINE qs_ot_ref_lwdn
379 :
380 : ! **************************************************************************************************
381 : !> \brief ...
382 : !> \param qs_ot_env ...
383 : !> \param C_OLD ...
384 : !> \param C_TMP ...
385 : !> \param C_NEW ...
386 : !> \param P ...
387 : !> \param SC ...
388 : !> \param norm_in ...
389 : !> \param update ...
390 : ! **************************************************************************************************
391 8560 : SUBROUTINE qs_ot_ref_poly(qs_ot_env, C_OLD, C_TMP, C_NEW, P, SC, norm_in, update)
392 : !
393 : TYPE(qs_ot_type) :: qs_ot_env
394 : TYPE(dbcsr_type), POINTER :: C_OLD, C_TMP, C_NEW, P
395 : TYPE(dbcsr_type) :: SC
396 : REAL(dp), INTENT(IN) :: norm_in
397 : LOGICAL, INTENT(IN) :: update
398 :
399 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_ot_ref_poly'
400 :
401 : INTEGER :: handle, irefine, k, n
402 : LOGICAL :: quick_exit
403 : REAL(dp) :: norm, norm_fro, norm_gct, occ_in, &
404 : occ_out, rescale
405 : TYPE(dbcsr_type), POINTER :: BUF1, BUF2, BUF_NOSYM, FT, FY
406 :
407 4280 : CALL timeset(routineN, handle)
408 : !
409 4280 : CALL dbcsr_get_info(C_NEW, nfullrows_total=n, nfullcols_total=k)
410 : !
411 4280 : BUF_NOSYM => qs_ot_env%buf1_k_k_nosym ! a buffer
412 4280 : BUF1 => qs_ot_env%buf1_k_k_sym ! a buffer
413 4280 : BUF2 => qs_ot_env%buf2_k_k_sym ! a buffer
414 4280 : FY => qs_ot_env%buf3_k_k_sym ! a buffer
415 4280 : FT => qs_ot_env%buf4_k_k_sym ! a buffer
416 : !
417 : ! initialize the norm (already computed in qs_ot_get_orbitals_ref)
418 4280 : norm = norm_in
419 : !
420 : ! can we do a quick exit?
421 4280 : quick_exit = .FALSE.
422 4280 : IF (norm < qs_ot_env%settings%eps_irac_quick_exit) quick_exit = .TRUE.
423 : !
424 : ! lets refine
425 4280 : rescale = 1.0_dp
426 4714 : DO irefine = 1, qs_ot_env%settings%max_irac
427 : !
428 : ! rescaling
429 4714 : IF (norm > 1.0_dp) THEN
430 12 : CALL dbcsr_scale(P, 1.0_dp/norm)
431 12 : rescale = rescale/SQRT(norm)
432 : END IF
433 : !
434 : ! get the refinement polynomial
435 : CALL qs_ot_refine(P, FY, BUF1, BUF2, qs_ot_env%settings%irac_degree, &
436 4714 : qs_ot_env%settings%eps_irac_filter_matrix)
437 : !
438 : ! collect the transformation
439 4714 : IF (irefine == 1) THEN
440 4280 : CALL dbcsr_copy(FT, FY, name='FT')
441 : ELSE
442 434 : CALL dbcsr_multiply('N', 'N', 1.0_dp, FT, FY, 0.0_dp, BUF1)
443 434 : IF (qs_ot_env%settings%eps_irac_filter_matrix > 0.0_dp) THEN
444 4 : occ_in = dbcsr_get_occupation(buf1)
445 4 : CALL dbcsr_filter(buf1, qs_ot_env%settings%eps_irac_filter_matrix)
446 4 : occ_out = dbcsr_get_occupation(buf1)
447 : END IF
448 434 : CALL dbcsr_copy(FT, BUF1, name='FT')
449 : END IF
450 : !
451 : ! quick exit if possible
452 4714 : IF (quick_exit) THEN
453 : EXIT
454 : END IF
455 : !
456 : ! P = FY^T * P * FY
457 2036 : CALL dbcsr_multiply('N', 'N', 1.0_dp, P, FY, 0.0_dp, BUF_NOSYM)
458 2036 : IF (qs_ot_env%settings%eps_irac_filter_matrix > 0.0_dp) THEN
459 8 : occ_in = dbcsr_get_occupation(buf_nosym)
460 8 : CALL dbcsr_filter(buf_nosym, qs_ot_env%settings%eps_irac_filter_matrix)
461 8 : occ_out = dbcsr_get_occupation(buf_nosym)
462 : END IF
463 2036 : CALL dbcsr_multiply('N', 'N', 1.0_dp, FY, BUF_NOSYM, 0.0_dp, P)
464 2036 : IF (qs_ot_env%settings%eps_irac_filter_matrix > 0.0_dp) THEN
465 8 : occ_in = dbcsr_get_occupation(p)
466 8 : CALL dbcsr_filter(p, qs_ot_env%settings%eps_irac_filter_matrix)
467 8 : occ_out = dbcsr_get_occupation(p)
468 : END IF
469 : !
470 : ! check ||P-1||_gct
471 2036 : CALL dbcsr_add_on_diag(P, -1.0_dp)
472 2036 : norm_fro = dbcsr_frobenius_norm(P)
473 2036 : norm_gct = dbcsr_gershgorin_norm(P)
474 2036 : CALL dbcsr_add_on_diag(P, 1.0_dp)
475 2036 : norm = MIN(norm_gct, norm_fro)
476 : !
477 : ! printing
478 : !
479 : ! blows up
480 2036 : IF (norm > 1.0E10_dp) THEN
481 : CALL cp_abort(__LOCATION__, &
482 : "Refinement blows up! "// &
483 : "We need you to improve the code, please post your input on "// &
484 0 : "the forum https://www.cp2k.org/")
485 : END IF
486 : !
487 : ! can we do a quick exit next step?
488 2036 : IF (norm < qs_ot_env%settings%eps_irac_quick_exit) quick_exit = .TRUE.
489 : !
490 : ! are we done?
491 4714 : IF (norm < qs_ot_env%settings%eps_irac) EXIT
492 : !
493 : END DO
494 : !
495 : ! C_NEW = C_NEW * FT * rescale
496 4280 : CALL dbcsr_multiply('N', 'N', rescale, C_OLD, FT, 0.0_dp, C_NEW)
497 4280 : IF (qs_ot_env%settings%eps_irac_filter_matrix > 0.0_dp) THEN
498 4 : occ_in = dbcsr_get_occupation(c_new)
499 4 : CALL dbcsr_filter(c_new, qs_ot_env%settings%eps_irac_filter_matrix)
500 4 : occ_out = dbcsr_get_occupation(c_new)
501 : END IF
502 : !
503 : ! update SC = SC * FY * rescale
504 4280 : IF (update) THEN
505 1754 : CALL dbcsr_multiply('N', 'N', rescale, SC, FT, 0.0_dp, C_TMP)
506 1754 : IF (qs_ot_env%settings%eps_irac_filter_matrix > 0.0_dp) THEN
507 4 : occ_in = dbcsr_get_occupation(c_tmp)
508 4 : CALL dbcsr_filter(c_tmp, qs_ot_env%settings%eps_irac_filter_matrix)
509 4 : occ_out = dbcsr_get_occupation(c_tmp)
510 : END IF
511 1754 : CALL dbcsr_copy(SC, C_TMP)
512 : END IF
513 : !
514 4280 : CALL timestop(handle)
515 4280 : END SUBROUTINE qs_ot_ref_poly
516 :
517 : ! **************************************************************************************************
518 : !> \brief ...
519 : !> \param qs_ot_env1 ...
520 : !> \return ...
521 : ! **************************************************************************************************
522 8344 : FUNCTION qs_ot_ref_update(qs_ot_env1) RESULT(update)
523 : !
524 : TYPE(qs_ot_type) :: qs_ot_env1
525 : LOGICAL :: update
526 :
527 8344 : update = .FALSE.
528 6802 : SELECT CASE (qs_ot_env1%settings%ot_method)
529 : CASE ("CG", "SD")
530 6802 : SELECT CASE (qs_ot_env1%settings%line_search_method)
531 : CASE ("2PNT")
532 6802 : IF (qs_ot_env1%line_search_count == 2) update = .TRUE.
533 : CASE DEFAULT
534 6802 : CPABORT("NYI")
535 : END SELECT
536 : CASE ("DIIS")
537 628 : update = .TRUE.
538 : CASE ("BROY", "LBFG")
539 : ! These minimizers retain positions or secants in one fixed REF chart.
540 628 : update = .FALSE.
541 : CASE DEFAULT
542 8344 : CPABORT("NYI")
543 : END SELECT
544 8344 : END FUNCTION qs_ot_ref_update
545 :
546 : ! **************************************************************************************************
547 : !> \brief ...
548 : !> \param qs_ot_env1 ...
549 : !> \param norm_in ...
550 : !> \param ortho_irac ...
551 : ! **************************************************************************************************
552 5626 : SUBROUTINE qs_ot_ref_decide(qs_ot_env1, norm_in, ortho_irac)
553 : !
554 : TYPE(qs_ot_type) :: qs_ot_env1
555 : REAL(dp), INTENT(IN) :: norm_in
556 : CHARACTER(LEN=*), INTENT(INOUT) :: ortho_irac
557 :
558 5626 : ortho_irac = qs_ot_env1%settings%ortho_irac
559 5626 : IF (norm_in < qs_ot_env1%settings%eps_irac_switch) ortho_irac = "POLY"
560 5626 : END SUBROUTINE qs_ot_ref_decide
561 :
562 : ! **************************************************************************************************
563 : !> \brief ...
564 : !> \param matrix_c ...
565 : !> \param matrix_s ...
566 : !> \param matrix_x ...
567 : !> \param matrix_sx ...
568 : !> \param matrix_gx_old ...
569 : !> \param matrix_dx ...
570 : !> \param qs_ot_env ...
571 : !> \param qs_ot_env1 ...
572 : ! **************************************************************************************************
573 11252 : SUBROUTINE qs_ot_get_orbitals_ref(matrix_c, matrix_s, matrix_x, matrix_sx, &
574 : matrix_gx_old, matrix_dx, qs_ot_env, qs_ot_env1)
575 : !
576 : TYPE(dbcsr_type), POINTER :: matrix_c, matrix_s, matrix_x, matrix_sx, &
577 : matrix_gx_old, matrix_dx
578 : TYPE(qs_ot_type) :: qs_ot_env, qs_ot_env1
579 :
580 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_ot_get_orbitals_ref'
581 :
582 : CHARACTER(LEN=4) :: ortho_irac
583 : INTEGER :: handle, k, n
584 : LOGICAL :: on_the_fly_loc, update
585 : REAL(dp) :: norm, norm_fro, norm_gct, occ_in, occ_out
586 : TYPE(dbcsr_type), POINTER :: C_NEW, C_OLD, C_TMP, D, G_OLD, P, S, SC
587 :
588 5626 : CALL timeset(routineN, handle)
589 :
590 5626 : CALL dbcsr_get_info(matrix_c, nfullrows_total=n, nfullcols_total=k)
591 : !
592 5626 : C_NEW => matrix_c
593 5626 : C_OLD => matrix_x ! need to be carefully updated for the gradient !
594 5626 : SC => matrix_sx ! need to be carefully updated for the gradient !
595 5626 : G_OLD => matrix_gx_old ! need to be carefully updated for localization !
596 5626 : D => matrix_dx ! need to be carefully updated for localization !
597 5626 : S => matrix_s
598 :
599 5626 : P => qs_ot_env%p_k_k_sym ! a buffer
600 5626 : C_TMP => qs_ot_env%buf1_n_k ! a buffer
601 : !
602 : ! do we need to update C_OLD and SC?
603 5626 : update = qs_ot_ref_update(qs_ot_env1)
604 : !
605 : ! do we want to on the fly localize?
606 : ! for the moment this is set from the input,
607 : ! later we might want to localize every n-step or
608 : ! when the sparsity increases...
609 5626 : on_the_fly_loc = qs_ot_env1%settings%on_the_fly_loc
610 : !
611 : ! compute SC = S*C
612 5626 : IF (ASSOCIATED(S)) THEN
613 5626 : CALL dbcsr_multiply('N', 'N', 1.0_dp, S, C_OLD, 0.0_dp, SC)
614 5626 : IF (qs_ot_env1%settings%eps_irac_filter_matrix > 0.0_dp) THEN
615 4 : occ_in = dbcsr_get_occupation(sc)
616 4 : CALL dbcsr_filter(sc, qs_ot_env1%settings%eps_irac_filter_matrix)
617 4 : occ_out = dbcsr_get_occupation(sc)
618 : END IF
619 : ELSE
620 0 : CALL dbcsr_copy(SC, C_OLD)
621 : END IF
622 : !
623 : ! compute P = C'*SC
624 5626 : CALL dbcsr_multiply('T', 'N', 1.0_dp, C_OLD, SC, 0.0_dp, P)
625 5626 : IF (qs_ot_env1%settings%eps_irac_filter_matrix > 0.0_dp) THEN
626 4 : occ_in = dbcsr_get_occupation(p)
627 4 : CALL dbcsr_filter(p, qs_ot_env1%settings%eps_irac_filter_matrix)
628 4 : occ_out = dbcsr_get_occupation(p)
629 : END IF
630 : !
631 : ! check ||P-1||_f and ||P-1||_gct
632 5626 : CALL dbcsr_add_on_diag(P, -1.0_dp)
633 5626 : norm_fro = dbcsr_frobenius_norm(P)
634 5626 : norm_gct = dbcsr_gershgorin_norm(P)
635 5626 : CALL dbcsr_add_on_diag(P, 1.0_dp)
636 5626 : norm = MIN(norm_gct, norm_fro)
637 5626 : CALL qs_ot_ref_decide(qs_ot_env1, norm, ortho_irac)
638 : !
639 : ! select the orthogonality method
640 868 : SELECT CASE (ortho_irac)
641 : CASE ("CHOL")
642 868 : CALL qs_ot_ref_chol(qs_ot_env, C_OLD, C_TMP, C_NEW, P, SC, update)
643 : CASE ("LWDN")
644 478 : CALL qs_ot_ref_lwdn(qs_ot_env, C_OLD, C_TMP, C_NEW, P, SC, update)
645 : CASE ("POLY")
646 4280 : CALL qs_ot_ref_poly(qs_ot_env, C_OLD, C_TMP, C_NEW, P, SC, norm, update)
647 : CASE DEFAULT
648 5626 : CPABORT("Wrong argument")
649 : END SELECT
650 : !
651 : ! We update the C_i+1 and localization
652 5626 : IF (update) THEN
653 2526 : IF (on_the_fly_loc) THEN
654 168 : CALL qs_ot_on_the_fly_localize(qs_ot_env, C_NEW, SC, G_OLD, D)
655 : END IF
656 2526 : CALL dbcsr_copy(C_OLD, C_NEW)
657 : END IF
658 :
659 5626 : IF (qs_ot_env%settings%do_rotation) THEN
660 14 : CALL qs_ot_generate_rotation(qs_ot_env)
661 : CALL dbcsr_multiply('N', 'N', 1.0_dp, C_NEW, qs_ot_env%rot_mat_u, &
662 14 : 0.0_dp, C_TMP)
663 14 : CALL dbcsr_copy(C_NEW, C_TMP)
664 : END IF
665 : !
666 5626 : CALL timestop(handle)
667 5626 : END SUBROUTINE qs_ot_get_orbitals_ref
668 :
669 : ! **************************************************************************************************
670 : !> \brief refinement polynomial of degree 2,3 and 4 (PRB 70, 193102 (2004))
671 : !> \param P ...
672 : !> \param FY ...
673 : !> \param P2 ...
674 : !> \param T ...
675 : !> \param irac_degree ...
676 : !> \param eps_irac_filter_matrix ...
677 : ! **************************************************************************************************
678 9428 : SUBROUTINE qs_ot_refine(P, FY, P2, T, irac_degree, eps_irac_filter_matrix)
679 : TYPE(dbcsr_type), INTENT(inout) :: P, FY, P2, T
680 : INTEGER, INTENT(in) :: irac_degree
681 : REAL(dp), INTENT(in) :: eps_irac_filter_matrix
682 :
683 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_ot_refine'
684 :
685 : INTEGER :: handle, k
686 : REAL(dp) :: occ_in, occ_out, r
687 :
688 4714 : CALL timeset(routineN, handle)
689 :
690 4714 : CALL dbcsr_get_info(P, nfullcols_total=k)
691 4714 : SELECT CASE (irac_degree)
692 : CASE (2)
693 : ! C_out = C_in * ( 15/8 * I - 10/8 * P + 3/8 * P^2)
694 0 : r = 3.0_dp/8.0_dp
695 0 : CALL dbcsr_multiply('N', 'N', r, P, P, 0.0_dp, FY)
696 0 : IF (eps_irac_filter_matrix > 0.0_dp) THEN
697 0 : occ_in = dbcsr_get_occupation(fy)
698 0 : CALL dbcsr_filter(fy, eps_irac_filter_matrix)
699 0 : occ_out = dbcsr_get_occupation(fy)
700 : END IF
701 0 : r = -10.0_dp/8.0_dp
702 0 : CALL dbcsr_add(FY, P, alpha_scalar=1.0_dp, beta_scalar=r)
703 0 : r = 15.0_dp/8.0_dp
704 0 : CALL dbcsr_add_on_diag(FY, alpha=r)
705 : CASE (3)
706 : ! C_out = C_in * ( 35/16 * I - 35/16 * P + 21/16 * P^2 - 5/16 P^3)
707 0 : CALL dbcsr_multiply('N', 'N', 1.0_dp, P, P, 0.0_dp, P2)
708 0 : IF (eps_irac_filter_matrix > 0.0_dp) THEN
709 0 : occ_in = dbcsr_get_occupation(p2)
710 0 : CALL dbcsr_filter(p2, eps_irac_filter_matrix)
711 0 : occ_out = dbcsr_get_occupation(p2)
712 : END IF
713 0 : r = -5.0_dp/16.0_dp
714 0 : CALL dbcsr_multiply('N', 'N', r, P2, P, 0.0_dp, FY)
715 0 : IF (eps_irac_filter_matrix > 0.0_dp) THEN
716 0 : occ_in = dbcsr_get_occupation(fy)
717 0 : CALL dbcsr_filter(fy, eps_irac_filter_matrix)
718 0 : occ_out = dbcsr_get_occupation(fy)
719 : END IF
720 0 : r = 21.0_dp/16.0_dp
721 0 : CALL dbcsr_add(FY, P2, alpha_scalar=1.0_dp, beta_scalar=r)
722 0 : r = -35.0_dp/16.0_dp
723 0 : CALL dbcsr_add(FY, P, alpha_scalar=1.0_dp, beta_scalar=r)
724 0 : r = 35.0_dp/16.0_dp
725 0 : CALL dbcsr_add_on_diag(FY, alpha=r)
726 : CASE (4)
727 : ! C_out = C_in * ( 315/128 * I - 420/128 * P + 378/128 * P^2 - 180/128 P^3 + 35/128 P^4 )
728 : ! = C_in * ( 315/128 * I - 420/128 * P + 378/128 * P^2 + ( - 180/128 * P + 35/128 * P^2 ) * P^2 )
729 4714 : CALL dbcsr_multiply('N', 'N', 1.0_dp, P, P, 0.0_dp, P2) ! P^2
730 4714 : IF (eps_irac_filter_matrix > 0.0_dp) THEN
731 8 : occ_in = dbcsr_get_occupation(p2)
732 8 : CALL dbcsr_filter(p2, eps_irac_filter_matrix)
733 8 : occ_out = dbcsr_get_occupation(p2)
734 : END IF
735 4714 : r = -180.0_dp/128.0_dp
736 4714 : CALL dbcsr_add(T, P, alpha_scalar=0.0_dp, beta_scalar=r) ! T=-180/128*P
737 4714 : r = 35.0_dp/128.0_dp
738 4714 : CALL dbcsr_add(T, P2, alpha_scalar=1.0_dp, beta_scalar=r) ! T=T+35/128*P^2
739 4714 : CALL dbcsr_multiply('N', 'N', 1.0_dp, T, P2, 0.0_dp, FY) ! Y=T*P^2
740 4714 : IF (eps_irac_filter_matrix > 0.0_dp) THEN
741 8 : occ_in = dbcsr_get_occupation(fy)
742 8 : CALL dbcsr_filter(fy, eps_irac_filter_matrix)
743 8 : occ_out = dbcsr_get_occupation(fy)
744 : END IF
745 4714 : r = 378.0_dp/128.0_dp
746 4714 : CALL dbcsr_add(FY, P2, alpha_scalar=1.0_dp, beta_scalar=r) ! Y=Y+378/128*P^2
747 4714 : r = -420.0_dp/128.0_dp
748 4714 : CALL dbcsr_add(FY, P, alpha_scalar=1.0_dp, beta_scalar=r) ! Y=Y-420/128*P
749 4714 : r = 315.0_dp/128.0_dp
750 4714 : CALL dbcsr_add_on_diag(FY, alpha=r) ! Y=Y+315/128*I
751 : CASE DEFAULT
752 4714 : CPABORT("This irac_order NYI")
753 : END SELECT
754 4714 : CALL timestop(handle)
755 4714 : END SUBROUTINE qs_ot_refine
756 :
757 : ! **************************************************************************************************
758 : !> \brief ...
759 : !> \param matrix_hc ...
760 : !> \param matrix_x ...
761 : !> \param matrix_sx ...
762 : !> \param matrix_gx ...
763 : !> \param qs_ot_env ...
764 : ! **************************************************************************************************
765 7064 : SUBROUTINE qs_ot_get_derivative_ref(matrix_hc, matrix_x, matrix_sx, matrix_gx, &
766 : qs_ot_env)
767 : TYPE(dbcsr_type), POINTER :: matrix_hc, matrix_x, matrix_sx, matrix_gx
768 : TYPE(qs_ot_type) :: qs_ot_env
769 :
770 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_ot_get_derivative_ref'
771 :
772 : INTEGER :: handle, k, n
773 : REAL(dp) :: occ_in, occ_out
774 : TYPE(dbcsr_type), POINTER :: C, CHC, G, HC, HC_work, SC
775 :
776 3532 : CALL timeset(routineN, handle)
777 :
778 3532 : CALL dbcsr_get_info(matrix_x, nfullrows_total=n, nfullcols_total=k)
779 : !
780 3532 : C => matrix_x ! NBsf*NOcc
781 3532 : SC => matrix_sx ! NBsf*NOcc need to be up2date
782 3532 : HC => matrix_hc ! NBsf*NOcc
783 3532 : G => matrix_gx ! NBsf*NOcc
784 3532 : CHC => qs_ot_env%buf1_k_k_sym ! buffer
785 :
786 3532 : IF (qs_ot_env%settings%do_rotation) THEN
787 : ! The physical orbitals are C_current=Q*U. Pull dE/dC_current
788 : ! back to the unrotated REF basis before projecting it.
789 8 : CALL qs_ot_rot_mat_derivative(qs_ot_env)
790 8 : HC_work => qs_ot_env%buf1_n_k
791 : CALL dbcsr_multiply('N', 'T', 1.0_dp, HC, qs_ot_env%rot_mat_u, &
792 8 : 0.0_dp, HC_work)
793 : ELSE
794 : HC_work => HC
795 : END IF
796 :
797 : ! C'*(H*C)
798 3532 : CALL dbcsr_multiply('T', 'N', 1.0_dp, C, HC_work, 0.0_dp, CHC)
799 3532 : IF (qs_ot_env%settings%eps_irac_filter_matrix > 0.0_dp) THEN
800 4 : occ_in = dbcsr_get_occupation(chc)
801 4 : CALL dbcsr_filter(chc, qs_ot_env%settings%eps_irac_filter_matrix)
802 4 : occ_out = dbcsr_get_occupation(chc)
803 : END IF
804 : ! (S*C)*(C'*H*C)
805 3532 : CALL dbcsr_multiply('N', 'N', 1.0_dp, SC, CHC, 0.0_dp, G)
806 3532 : IF (qs_ot_env%settings%eps_irac_filter_matrix > 0.0_dp) THEN
807 4 : occ_in = dbcsr_get_occupation(g)
808 4 : CALL dbcsr_filter(g, qs_ot_env%settings%eps_irac_filter_matrix)
809 4 : occ_out = dbcsr_get_occupation(g)
810 : END IF
811 : ! G = 2*(1-S*C*C')*H*C
812 3532 : CALL dbcsr_add(G, HC_work, alpha_scalar=-1.0_dp, beta_scalar=1.0_dp)
813 : !
814 3532 : CALL timestop(handle)
815 3532 : END SUBROUTINE qs_ot_get_derivative_ref
816 :
817 : ! **************************************************************************************************
818 : !> \brief Transpose a square DBCSR matrix while retaining its target distribution.
819 : !> \param matrix source matrix
820 : !> \param transposed transposed matrix
821 : !> \param identity_template square matrix template
822 : ! **************************************************************************************************
823 6166 : SUBROUTINE qs_ot_square_transpose(matrix, transposed, identity_template)
824 : TYPE(dbcsr_type) :: matrix, transposed, identity_template
825 :
826 : TYPE(dbcsr_type) :: identity
827 :
828 6166 : CALL dbcsr_copy(transposed, matrix, name='square_transposed')
829 6166 : CALL dbcsr_copy(identity, identity_template, name='transpose_identity')
830 6166 : CALL dbcsr_set(identity, 0.0_dp)
831 6166 : CALL dbcsr_add_on_diag(identity, 1.0_dp)
832 6166 : CALL dbcsr_multiply('T', 'N', 1.0_dp, matrix, identity, 0.0_dp, transposed)
833 6166 : CALL dbcsr_release(identity)
834 :
835 6166 : END SUBROUTINE qs_ot_square_transpose
836 :
837 : ! **************************************************************************************************
838 : !> \brief computes p=x*S*x and the matrix functionals related matrices
839 : !> \param matrix_x ...
840 : !> \param matrix_sx ...
841 : !> \param qs_ot_env ...
842 : ! **************************************************************************************************
843 333519 : SUBROUTINE qs_ot_get_p(matrix_x, matrix_sx, qs_ot_env)
844 :
845 : TYPE(dbcsr_type), POINTER :: matrix_x, matrix_sx
846 : TYPE(qs_ot_type) :: qs_ot_env
847 :
848 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_ot_get_p'
849 : REAL(KIND=dp), PARAMETER :: rone = 1.0_dp, rzero = 0.0_dp
850 :
851 : INTEGER :: handle, k, max_iter, n
852 : LOGICAL :: converged
853 : REAL(KIND=dp) :: max_ev, min_ev, threshold
854 :
855 111173 : CALL timeset(routineN, handle)
856 :
857 111173 : CALL dbcsr_get_info(matrix_x, nfullrows_total=n, nfullcols_total=k)
858 :
859 : ! get the overlap
860 : CALL dbcsr_multiply('T', 'N', rone, matrix_x, matrix_sx, rzero, &
861 111173 : qs_ot_env%matrix_p)
862 :
863 : ! get an upper bound for the largest eigenvalue
864 : ! try using lancos first and fall back to gershgorin norm if it fails
865 111173 : max_iter = 30; threshold = 1.0E-03_dp
866 111173 : CALL arnoldi_extremal(qs_ot_env%matrix_p, max_ev, min_ev, converged, threshold, max_iter)
867 111173 : qs_ot_env%largest_eval_upper_bound = MAX(max_ev, ABS(min_ev))
868 :
869 111173 : IF (.NOT. converged) qs_ot_env%largest_eval_upper_bound = dbcsr_gershgorin_norm(qs_ot_env%matrix_p)
870 111173 : CALL decide_strategy(qs_ot_env)
871 111173 : IF (qs_ot_env%do_taylor) THEN
872 59864 : CALL qs_ot_p2m_taylor(qs_ot_env)
873 : ELSE
874 51309 : CALL qs_ot_p2m_diag(qs_ot_env)
875 : END IF
876 :
877 111173 : IF (qs_ot_env%settings%do_rotation) THEN
878 3142 : CALL qs_ot_generate_rotation(qs_ot_env)
879 : END IF
880 :
881 111173 : CALL timestop(handle)
882 :
883 111173 : END SUBROUTINE qs_ot_get_p
884 :
885 : ! **************************************************************************************************
886 : !> \brief computes the rotation matrix rot_mat_u that is associated to a given
887 : !> rot_mat_x using rot_mat_u=exp(rot_mat_x)
888 : !> \param qs_ot_env a valid qs_ot_env
889 : !> \par History
890 : !> 08.2004 created [Joost VandeVondele]
891 : !> 12.2024 Rewrite to use only real matrices [Ole Schuett]
892 : ! **************************************************************************************************
893 3156 : SUBROUTINE qs_ot_generate_rotation(qs_ot_env)
894 :
895 : TYPE(qs_ot_type) :: qs_ot_env
896 :
897 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_ot_generate_rotation'
898 :
899 : INTEGER :: handle, k
900 3156 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: exp_evals_im, exp_evals_re
901 : TYPE(dbcsr_type) :: buf_1, buf_2
902 :
903 3156 : CALL timeset(routineN, handle)
904 :
905 3156 : CALL dbcsr_get_info(qs_ot_env%rot_mat_x, nfullrows_total=k)
906 :
907 3156 : IF (k /= 0) THEN
908 : ! We want to compute: rot_mat_u = exp(i*rot_mat_x)
909 :
910 : ! Diagonalize: matrix = i*rot_mat_x.
911 : ! Note that matrix is imaginary and hermitian because rot_mat_x is real and anti-symmetric.
912 : CALL cp_dbcsr_heevd(matrix_im=qs_ot_env%rot_mat_x, & ! matrix_re omitted because it's zero
913 : eigenvectors_re=qs_ot_env%rot_mat_evec_re, &
914 : eigenvectors_im=qs_ot_env%rot_mat_evec_im, &
915 : eigenvalues=qs_ot_env%rot_mat_evals, &
916 : para_env=qs_ot_env%para_env, &
917 3104 : blacs_env=qs_ot_env%blacs_env)
918 :
919 : ! Compute: exp_evals = EXP(-i*rot_mat_evals)
920 12416 : ALLOCATE (exp_evals_re(k), exp_evals_im(k))
921 16758 : exp_evals_re(:) = COS(-qs_ot_env%rot_mat_evals(:))
922 16758 : exp_evals_im(:) = SIN(-qs_ot_env%rot_mat_evals(:))
923 :
924 : ! Compute: rot_mat_u = \sum_ij exp_evals_ij * |rot_mat_evec_i> <rot_mat_evec_j|
925 : ! Note that we need only two matrix multiplications because rot_mat_u is real.
926 3104 : CALL dbcsr_copy(buf_1, qs_ot_env%rot_mat_evec_re, name="buf_1")
927 3104 : CALL dbcsr_scale_by_vector(buf_1, alpha=exp_evals_re, side='right')
928 3104 : CALL dbcsr_copy(buf_2, qs_ot_env%rot_mat_evec_im, name="buf_2")
929 3104 : CALL dbcsr_scale_by_vector(buf_2, alpha=exp_evals_im, side='right')
930 3104 : CALL dbcsr_add(buf_1, buf_2, alpha_scalar=+1.0_dp, beta_scalar=-1.0_dp)
931 3104 : CALL dbcsr_multiply('N', 'T', 1.0_dp, buf_1, qs_ot_env%rot_mat_evec_re, 0.0_dp, qs_ot_env%rot_mat_u)
932 :
933 3104 : CALL dbcsr_copy(buf_1, qs_ot_env%rot_mat_evec_im)
934 3104 : CALL dbcsr_scale_by_vector(buf_1, alpha=exp_evals_re, side='right')
935 3104 : CALL dbcsr_copy(buf_2, qs_ot_env%rot_mat_evec_re)
936 3104 : CALL dbcsr_scale_by_vector(buf_2, alpha=exp_evals_im, side='right')
937 3104 : CALL dbcsr_add(buf_1, buf_2, alpha_scalar=+1.0_dp, beta_scalar=+1.0_dp)
938 3104 : CALL dbcsr_multiply('N', 'T', 1.0_dp, buf_1, qs_ot_env%rot_mat_evec_im, 1.0_dp, qs_ot_env%rot_mat_u)
939 :
940 : ! Clean up.
941 3104 : CALL dbcsr_release(buf_1)
942 3104 : CALL dbcsr_release(buf_2)
943 3104 : DEALLOCATE (exp_evals_re, exp_evals_im)
944 : END IF
945 :
946 3156 : CALL timestop(handle)
947 :
948 6312 : END SUBROUTINE qs_ot_generate_rotation
949 :
950 : ! **************************************************************************************************
951 : !> \brief computes the derivative fields with respect to rot_mat_x
952 : !> \param qs_ot_env valid qs_ot_env. In particular qs_ot_generate_rotation has to be called before
953 : !> and the rot_mat_dedu matrix has to be up to date
954 : !> \par History
955 : !> 08.2004 created [ Joost VandeVondele ]
956 : !> 12.2024 Rewrite to use only real matrices [Ole Schuett]
957 : ! **************************************************************************************************
958 1596 : SUBROUTINE qs_ot_rot_mat_derivative(qs_ot_env)
959 : TYPE(qs_ot_type) :: qs_ot_env
960 :
961 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_ot_rot_mat_derivative'
962 :
963 : INTEGER :: col, handle, i, iblock, j, k, max_blocks, nblocks, row
964 1596 : INTEGER, ALLOCATABLE, DIMENSION(:) :: cols, rows
965 1596 : INTEGER, DIMENSION(:), POINTER :: col_blk_offset, col_blk_size, row_blk_offset, row_blk_size
966 : REAL(KIND=dp) :: e1, e2
967 : TYPE(dbcsr_type) :: outer_deriv_re, outer_deriv_im, mat_buf, &
968 : inner_deriv_re, inner_deriv_im
969 : TYPE(dbcsr_distribution_type) :: dist
970 : TYPE(dbcsr_iterator_type) :: iter
971 1596 : REAL(dp), DIMENSION(:, :), POINTER :: block_in_re, block_in_im, block_out_re, block_out_im
972 : LOGICAL :: duplicate, found_in_im, found_in_re, found_out_im, &
973 : found_out_re
974 : REAL(KIND=dp) :: im_part, re_part
975 : COMPLEX(dp) :: cval_in, cval_out
976 1596 : CALL timeset(routineN, handle)
977 :
978 1596 : CALL dbcsr_get_info(qs_ot_env%rot_mat_u, nfullrows_total=k)
979 1596 : IF (k /= 0) THEN
980 1570 : CALL dbcsr_copy(qs_ot_env%matrix_buf1, qs_ot_env%rot_mat_dedu)
981 : ! now we get to the derivative wrt the antisymmetric matrix rot_mat_x
982 1570 : CALL dbcsr_copy(mat_buf, qs_ot_env%rot_mat_dedu, "mat_buf")
983 :
984 : ! inner_deriv_ij = <rot_mat_evec_i| rot_mat_dedu |rot_mat_evec_j>
985 1570 : CALL dbcsr_copy(inner_deriv_re, qs_ot_env%rot_mat_dedu, "inner_deriv_re") ! TODO just create
986 1570 : CALL dbcsr_copy(inner_deriv_im, qs_ot_env%rot_mat_dedu, "inner_deriv_im") ! TODO just create
987 :
988 1570 : CALL dbcsr_multiply('T', 'N', +1.0_dp, qs_ot_env%rot_mat_dedu, qs_ot_env%rot_mat_evec_im, 0.0_dp, mat_buf)
989 1570 : CALL dbcsr_multiply('T', 'N', +1.0_dp, qs_ot_env%rot_mat_evec_im, mat_buf, 0.0_dp, inner_deriv_re)
990 1570 : CALL dbcsr_multiply('T', 'N', +1.0_dp, qs_ot_env%rot_mat_evec_re, mat_buf, 0.0_dp, inner_deriv_im)
991 :
992 1570 : CALL dbcsr_multiply('T', 'N', +1.0_dp, qs_ot_env%rot_mat_dedu, qs_ot_env%rot_mat_evec_re, 0.0_dp, mat_buf)
993 1570 : CALL dbcsr_multiply('T', 'N', +1.0_dp, qs_ot_env%rot_mat_evec_re, mat_buf, 1.0_dp, inner_deriv_re)
994 1570 : CALL dbcsr_multiply('T', 'N', -1.0_dp, qs_ot_env%rot_mat_evec_im, mat_buf, 1.0_dp, inner_deriv_im)
995 :
996 : ! Real and imaginary products can have different sparse block patterns.
997 : ! Form their union explicitly and treat a missing partner block as zero.
998 1570 : max_blocks = 0
999 1570 : CALL dbcsr_iterator_start(iter, inner_deriv_re)
1000 2354 : DO WHILE (dbcsr_iterator_blocks_left(iter))
1001 784 : CALL dbcsr_iterator_next_block(iter, row, col)
1002 784 : max_blocks = max_blocks + 1
1003 : END DO
1004 1570 : CALL dbcsr_iterator_stop(iter)
1005 1570 : CALL dbcsr_iterator_start(iter, inner_deriv_im)
1006 2354 : DO WHILE (dbcsr_iterator_blocks_left(iter))
1007 784 : CALL dbcsr_iterator_next_block(iter, row, col)
1008 784 : max_blocks = max_blocks + 1
1009 : END DO
1010 1570 : CALL dbcsr_iterator_stop(iter)
1011 :
1012 6280 : ALLOCATE (rows(MAX(max_blocks, 1)), cols(MAX(max_blocks, 1)))
1013 1570 : nblocks = 0
1014 1570 : CALL dbcsr_iterator_start(iter, inner_deriv_re)
1015 2354 : DO WHILE (dbcsr_iterator_blocks_left(iter))
1016 784 : CALL dbcsr_iterator_next_block(iter, row, col)
1017 784 : duplicate = .FALSE.
1018 784 : DO iblock = 1, nblocks
1019 784 : duplicate = duplicate .OR. (rows(iblock) == row .AND. cols(iblock) == col)
1020 : END DO
1021 2354 : IF (.NOT. duplicate) THEN
1022 784 : nblocks = nblocks + 1
1023 784 : rows(nblocks) = row
1024 784 : cols(nblocks) = col
1025 : END IF
1026 : END DO
1027 1570 : CALL dbcsr_iterator_stop(iter)
1028 1570 : CALL dbcsr_iterator_start(iter, inner_deriv_im)
1029 2354 : DO WHILE (dbcsr_iterator_blocks_left(iter))
1030 784 : CALL dbcsr_iterator_next_block(iter, row, col)
1031 784 : duplicate = .FALSE.
1032 1568 : DO iblock = 1, nblocks
1033 1568 : duplicate = duplicate .OR. (rows(iblock) == row .AND. cols(iblock) == col)
1034 : END DO
1035 2354 : IF (.NOT. duplicate) THEN
1036 0 : nblocks = nblocks + 1
1037 0 : rows(nblocks) = row
1038 0 : cols(nblocks) = col
1039 : END IF
1040 : END DO
1041 1570 : CALL dbcsr_iterator_stop(iter)
1042 :
1043 : CALL dbcsr_get_info(inner_deriv_re, distribution=dist, row_blk_size=row_blk_size, &
1044 1570 : col_blk_size=col_blk_size)
1045 : CALL dbcsr_create(outer_deriv_re, "outer_deriv_re", dist, dbcsr_type_no_symmetry, &
1046 1570 : row_blk_size, col_blk_size)
1047 : CALL dbcsr_create(outer_deriv_im, "outer_deriv_im", dist, dbcsr_type_no_symmetry, &
1048 1570 : row_blk_size, col_blk_size)
1049 1570 : IF (nblocks > 0) THEN
1050 784 : CALL dbcsr_reserve_blocks(outer_deriv_re, rows=rows(1:nblocks), cols=cols(1:nblocks))
1051 784 : CALL dbcsr_reserve_blocks(outer_deriv_im, rows=rows(1:nblocks), cols=cols(1:nblocks))
1052 : END IF
1053 1570 : CALL dbcsr_finalize(outer_deriv_re)
1054 1570 : CALL dbcsr_finalize(outer_deriv_im)
1055 1570 : CALL dbcsr_set(outer_deriv_re, 0.0_dp)
1056 1570 : CALL dbcsr_set(outer_deriv_im, 0.0_dp)
1057 :
1058 1570 : CALL dbcsr_get_info(outer_deriv_re, row_blk_offset=row_blk_offset, col_blk_offset=col_blk_offset)
1059 1570 : CALL dbcsr_iterator_start(iter, outer_deriv_re)
1060 2354 : DO WHILE (dbcsr_iterator_blocks_left(iter))
1061 784 : CALL dbcsr_iterator_next_block(iter, row, col)
1062 784 : CALL dbcsr_get_block_p(inner_deriv_re, row, col, block_in_re, found_in_re)
1063 784 : CALL dbcsr_get_block_p(inner_deriv_im, row, col, block_in_im, found_in_im)
1064 784 : CALL dbcsr_get_block_p(outer_deriv_re, row, col, block_out_re, found_out_re)
1065 784 : CALL dbcsr_get_block_p(outer_deriv_im, row, col, block_out_im, found_out_im)
1066 784 : CPASSERT(found_out_re .AND. found_out_im)
1067 :
1068 5943 : DO i = 1, SIZE(block_out_re, 1)
1069 25726 : DO j = 1, SIZE(block_out_re, 2)
1070 21353 : e1 = qs_ot_env%rot_mat_evals(row_blk_offset(row) + i - 1)
1071 21353 : e2 = qs_ot_env%rot_mat_evals(col_blk_offset(col) + j - 1)
1072 21353 : re_part = 0.0_dp
1073 21353 : im_part = 0.0_dp
1074 21353 : IF (found_in_re) re_part = block_in_re(i, j)
1075 21353 : IF (found_in_im) im_part = block_in_im(i, j)
1076 21353 : cval_in = CMPLX(re_part, im_part, dp)
1077 21353 : cval_out = cval_in*cint(e1, e2)
1078 21353 : block_out_re(i, j) = REAL(cval_out)
1079 24942 : block_out_im(i, j) = AIMAG(cval_out)
1080 : END DO
1081 : END DO
1082 : END DO
1083 1570 : CALL dbcsr_iterator_stop(iter)
1084 1570 : DEALLOCATE (rows, cols)
1085 1570 : CALL dbcsr_release(inner_deriv_re)
1086 1570 : CALL dbcsr_release(inner_deriv_im)
1087 :
1088 : ! Compute: matrix_buf1 = \sum_i outer_deriv_ij * |rot_mat_evec_i> <rot_mat_evec_j|
1089 1570 : CALL dbcsr_multiply('N', 'N', +1.0_dp, qs_ot_env%rot_mat_evec_re, outer_deriv_re, 0.0_dp, mat_buf)
1090 1570 : CALL dbcsr_multiply('N', 'N', -1.0_dp, qs_ot_env%rot_mat_evec_im, outer_deriv_im, 1.0_dp, mat_buf)
1091 1570 : CALL dbcsr_multiply('N', 'T', +1.0_dp, mat_buf, qs_ot_env%rot_mat_evec_re, 0.0_dp, qs_ot_env%matrix_buf1)
1092 :
1093 1570 : CALL dbcsr_multiply('N', 'N', +1.0_dp, qs_ot_env%rot_mat_evec_re, outer_deriv_im, 0.0_dp, mat_buf)
1094 1570 : CALL dbcsr_multiply('N', 'N', +1.0_dp, qs_ot_env%rot_mat_evec_im, outer_deriv_re, 1.0_dp, mat_buf)
1095 1570 : CALL dbcsr_multiply('N', 'T', +1.0_dp, mat_buf, qs_ot_env%rot_mat_evec_im, 1.0_dp, qs_ot_env%matrix_buf1)
1096 :
1097 : ! Account for anti-symmetry of rot_mat_x without relying on
1098 : ! STRICT-only matrix buffers. REF uses the same finite chart.
1099 1570 : CALL qs_ot_square_transpose(qs_ot_env%matrix_buf1, mat_buf, qs_ot_env%rot_mat_u)
1100 1570 : CALL dbcsr_copy(qs_ot_env%rot_mat_gx, mat_buf)
1101 : CALL dbcsr_add(qs_ot_env%rot_mat_gx, qs_ot_env%matrix_buf1, &
1102 1570 : alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
1103 :
1104 1570 : CALL dbcsr_release(mat_buf)
1105 1570 : CALL dbcsr_release(outer_deriv_re)
1106 10990 : CALL dbcsr_release(outer_deriv_im)
1107 : END IF
1108 3192 : CALL timestop(handle)
1109 : CONTAINS
1110 :
1111 : ! **************************************************************************************************
1112 : !> \brief ...
1113 : !> \param e1 ...
1114 : !> \param e2 ...
1115 : !> \return ...
1116 : ! **************************************************************************************************
1117 21353 : FUNCTION cint(e1, e2)
1118 : REAL(KIND=dp) :: e1, e2
1119 : COMPLEX(KIND=dp) :: cint
1120 :
1121 : COMPLEX(KIND=dp) :: l1, l2, x
1122 : INTEGER :: I
1123 :
1124 21353 : l1 = (0.0_dp, -1.0_dp)*e1
1125 21353 : l2 = (0.0_dp, -1.0_dp)*e2
1126 21353 : IF (ABS(l1 - l2) > 0.5_dp) THEN
1127 1424 : cint = (EXP(l1) - EXP(l2))/(l1 - l2)
1128 : ELSE
1129 : x = 1.0_dp
1130 : cint = 0.0_dp
1131 338793 : DO I = 1, 16
1132 318864 : cint = cint + x
1133 338793 : x = x*(l1 - l2)/REAL(I + 1, KIND=dp)
1134 : END DO
1135 19929 : cint = cint*EXP(l2)
1136 : END IF
1137 21353 : END FUNCTION cint
1138 : END SUBROUTINE qs_ot_rot_mat_derivative
1139 :
1140 : ! **************************************************************************************************
1141 : !> \brief decide strategy
1142 : !> tries to decide if the taylor expansion of cos(sqrt(xsx)) converges rapidly enough
1143 : !> to make a taylor expansion of the functions cos(sqrt(xsx)) and sin(sqrt(xsx))/sqrt(xsx)
1144 : !> and their derivatives faster than their computation based on diagonalization since xsx can
1145 : !> be very small, especially during dynamics, only a few terms might indeed be needed we find
1146 : !> the necessary order N to have largest_eval_upper_bound**(N+1)/(2(N+1))! < eps_taylor
1147 : !> \param qs_ot_env ...
1148 : ! **************************************************************************************************
1149 111173 : SUBROUTINE decide_strategy(qs_ot_env)
1150 : TYPE(qs_ot_type) :: qs_ot_env
1151 :
1152 : INTEGER :: N
1153 : REAL(KIND=dp) :: num_error
1154 :
1155 111173 : qs_ot_env%do_taylor = .FALSE.
1156 111173 : N = 0
1157 111173 : num_error = qs_ot_env%largest_eval_upper_bound/(2.0_dp)
1158 472429 : DO WHILE (num_error > qs_ot_env%settings%eps_taylor .AND. N <= qs_ot_env%settings%max_taylor)
1159 361256 : N = N + 1
1160 407279 : num_error = num_error*qs_ot_env%largest_eval_upper_bound/REAL((2*N + 1)*(2*N + 2), KIND=dp)
1161 : END DO
1162 : ! An O(P) matrix-function term can be below eps_taylor while its O(X)
1163 : ! orbital derivative is still needed near P=X^T*S*X=0.
1164 111173 : qs_ot_env%taylor_order = MAX(1, N)
1165 111173 : IF (qs_ot_env%taylor_order <= qs_ot_env%settings%max_taylor) THEN
1166 59864 : qs_ot_env%do_taylor = .TRUE.
1167 : END IF
1168 :
1169 111173 : END SUBROUTINE decide_strategy
1170 :
1171 : ! **************************************************************************************************
1172 : !> \brief c=(c0*cos(p^0.5)+x*sin(p^0.5)*p^(-0.5)) x rot_mat_u
1173 : !> this assumes that x is already ortho to S*C0, and that p is x*S*x
1174 : !> rot_mat_u is an optional rotation matrix
1175 : !> \param matrix_c ...
1176 : !> \param matrix_x ...
1177 : !> \param qs_ot_env ...
1178 : ! **************************************************************************************************
1179 206450 : SUBROUTINE qs_ot_get_orbitals(matrix_c, matrix_x, qs_ot_env)
1180 :
1181 : TYPE(dbcsr_type), POINTER :: matrix_c, matrix_x
1182 : TYPE(qs_ot_type) :: qs_ot_env
1183 :
1184 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_ot_get_orbitals'
1185 : REAL(KIND=dp), PARAMETER :: rone = 1.0_dp, rzero = 0.0_dp
1186 :
1187 : INTEGER :: handle, k, n
1188 : TYPE(dbcsr_type), POINTER :: matrix_kk
1189 :
1190 103225 : CALL timeset(routineN, handle)
1191 :
1192 103225 : CALL dbcsr_get_info(matrix_x, nfullrows_total=n, nfullcols_total=k)
1193 :
1194 : ! rotate the multiplying matrices cosp and sinp instead of the result,
1195 : ! this should be cheaper for large basis sets
1196 103225 : IF (qs_ot_env%settings%do_rotation) THEN
1197 2938 : matrix_kk => qs_ot_env%matrix_buf1
1198 : CALL dbcsr_multiply('N', 'N', rone, qs_ot_env%matrix_cosp, &
1199 2938 : qs_ot_env%rot_mat_u, rzero, matrix_kk)
1200 : ELSE
1201 100287 : matrix_kk => qs_ot_env%matrix_cosp
1202 : END IF
1203 :
1204 : CALL dbcsr_multiply('N', 'N', rone, qs_ot_env%matrix_c0, matrix_kk, &
1205 103225 : rzero, matrix_c)
1206 :
1207 103225 : IF (qs_ot_env%settings%do_rotation) THEN
1208 2938 : matrix_kk => qs_ot_env%matrix_buf1
1209 : CALL dbcsr_multiply('N', 'N', rone, qs_ot_env%matrix_sinp, &
1210 2938 : qs_ot_env%rot_mat_u, rzero, matrix_kk)
1211 : ELSE
1212 100287 : matrix_kk => qs_ot_env%matrix_sinp
1213 : END IF
1214 : CALL dbcsr_multiply('N', 'N', rone, matrix_x, matrix_kk, &
1215 103225 : rone, matrix_c)
1216 :
1217 103225 : CALL timestop(handle)
1218 :
1219 103225 : END SUBROUTINE qs_ot_get_orbitals
1220 :
1221 : ! **************************************************************************************************
1222 : !> \brief this routines computes dE/dx=dx, with dx ortho to sc0
1223 : !> needs dE/dC=hc,C0,X,SX,p
1224 : !> if preconditioned it will not be the derivative, but the lagrangian multiplier
1225 : !> is changed so that P*dE/dx is the right derivative (i.e. in the allowed subspace)
1226 : !> \param matrix_hc ...
1227 : !> \param matrix_x ...
1228 : !> \param matrix_sx ...
1229 : !> \param matrix_gx ...
1230 : !> \param qs_ot_env ...
1231 : ! **************************************************************************************************
1232 79753 : SUBROUTINE qs_ot_get_derivative(matrix_hc, matrix_x, matrix_sx, matrix_gx, &
1233 : qs_ot_env)
1234 : TYPE(dbcsr_type), POINTER :: matrix_hc, matrix_x, matrix_sx, matrix_gx
1235 : TYPE(qs_ot_type) :: qs_ot_env
1236 :
1237 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_ot_get_derivative'
1238 : REAL(KIND=dp), PARAMETER :: rone = 1.0_dp, rzero = 0.0_dp
1239 :
1240 : INTEGER :: handle
1241 : TYPE(dbcsr_type), POINTER :: matrix_hc_local
1242 :
1243 79753 : CALL timeset(routineN, handle)
1244 :
1245 79753 : NULLIFY (matrix_hc_local)
1246 :
1247 : ! could in principle be taken inside qs_ot_get_derivative_* for increased efficiency
1248 : ! create a local rotated version of matrix_hc leaving matrix_hc untouched (needed
1249 : ! for lagrangian multipliers)
1250 79753 : IF (qs_ot_env%settings%do_rotation) THEN
1251 1588 : CALL dbcsr_copy(matrix_gx, matrix_hc) ! use gx as temporary
1252 1588 : CALL dbcsr_init_p(matrix_hc_local)
1253 1588 : CALL dbcsr_copy(matrix_hc_local, matrix_hc, name='matrix_hc_local')
1254 1588 : CALL dbcsr_set(matrix_hc_local, 0.0_dp)
1255 1588 : CALL dbcsr_multiply('N', 'T', rone, matrix_gx, qs_ot_env%rot_mat_u, rzero, matrix_hc_local)
1256 : ELSE
1257 78165 : matrix_hc_local => matrix_hc
1258 : END IF
1259 :
1260 79753 : IF (qs_ot_env%do_taylor) THEN
1261 44113 : CALL qs_ot_get_derivative_taylor(matrix_hc_local, matrix_x, matrix_sx, matrix_gx, qs_ot_env)
1262 : ELSE
1263 35640 : CALL qs_ot_get_derivative_diag(matrix_hc_local, matrix_x, matrix_sx, matrix_gx, qs_ot_env)
1264 : END IF
1265 :
1266 79753 : CALL qs_ot_project_gradient(qs_ot_env, matrix_gx)
1267 : ! also treat the rot_mat gradient here
1268 79753 : IF (qs_ot_env%settings%do_rotation) THEN
1269 1588 : CALL qs_ot_rot_mat_derivative(qs_ot_env)
1270 : END IF
1271 :
1272 79753 : IF (qs_ot_env%settings%do_rotation) THEN
1273 1588 : CALL dbcsr_release_p(matrix_hc_local)
1274 : END IF
1275 :
1276 79753 : CALL timestop(handle)
1277 :
1278 79753 : END SUBROUTINE qs_ot_get_derivative
1279 :
1280 : ! **************************************************************************************************
1281 : !> \brief Project a real STRICT gradient, or an identity-metric fallback direction.
1282 : !> \param qs_ot_env OT channel state
1283 : !> \param matrix_gx real component to project in place
1284 : !> \param use_preconditioner use the attached preconditioner, default true
1285 : ! **************************************************************************************************
1286 79753 : SUBROUTINE qs_ot_project_gradient(qs_ot_env, matrix_gx, use_preconditioner)
1287 : TYPE(qs_ot_type) :: qs_ot_env
1288 : TYPE(dbcsr_type), POINTER :: matrix_gx
1289 : LOGICAL, INTENT(IN), OPTIONAL :: use_preconditioner
1290 :
1291 : REAL(KIND=dp), PARAMETER :: rone = 1.0_dp, rzero = 0.0_dp
1292 :
1293 : LOGICAL :: preconditioned
1294 : TYPE(dbcsr_type), POINTER :: target_re
1295 :
1296 79753 : preconditioned = .TRUE.
1297 79753 : IF (PRESENT(use_preconditioner)) preconditioned = use_preconditioner
1298 0 : IF (.NOT. preconditioned) qs_ot_env%os_valid = .FALSE.
1299 79753 : IF (preconditioned .AND. ASSOCIATED(qs_ot_env%preconditioner)) THEN
1300 69255 : target_re => qs_ot_env%matrix_psc0
1301 : ELSE
1302 10498 : target_re => qs_ot_env%matrix_sc0
1303 : END IF
1304 : ! first make the matrix os if not yet valid
1305 79753 : IF (.NOT. qs_ot_env%os_valid) THEN
1306 : ! this assumes that the preconditioner is a single matrix
1307 : ! that maps sc0 onto psc0
1308 :
1309 8630 : IF (preconditioned .AND. ASSOCIATED(qs_ot_env%preconditioner)) THEN
1310 : CALL apply_preconditioner(qs_ot_env%preconditioner, qs_ot_env%matrix_sc0, &
1311 7632 : qs_ot_env%matrix_psc0)
1312 : END IF
1313 : CALL dbcsr_multiply('T', 'N', rone, &
1314 : qs_ot_env%matrix_sc0, target_re, &
1315 8630 : rzero, qs_ot_env%matrix_os)
1316 : CALL cp_dbcsr_cholesky_decompose(qs_ot_env%matrix_os, &
1317 8630 : para_env=qs_ot_env%para_env, blacs_env=qs_ot_env%blacs_env)
1318 : CALL cp_dbcsr_cholesky_invert(qs_ot_env%matrix_os, &
1319 : para_env=qs_ot_env%para_env, blacs_env=qs_ot_env%blacs_env, &
1320 8630 : uplo_to_full=.TRUE.)
1321 8630 : qs_ot_env%os_valid = preconditioned
1322 : END IF
1323 : CALL dbcsr_multiply('T', 'N', rone, target_re, matrix_gx, &
1324 79753 : rzero, qs_ot_env%matrix_buf1_ortho)
1325 : CALL dbcsr_multiply('N', 'N', rone, qs_ot_env%matrix_os, &
1326 79753 : qs_ot_env%matrix_buf1_ortho, rzero, qs_ot_env%matrix_buf2_ortho)
1327 : CALL dbcsr_multiply('N', 'N', -rone, qs_ot_env%matrix_sc0, &
1328 79753 : qs_ot_env%matrix_buf2_ortho, rone, matrix_gx)
1329 :
1330 79753 : END SUBROUTINE qs_ot_project_gradient
1331 :
1332 : ! **************************************************************************************************
1333 : !> \brief ...
1334 : !> \param matrix_hc ...
1335 : !> \param matrix_x ...
1336 : !> \param matrix_sx ...
1337 : !> \param matrix_gx ...
1338 : !> \param qs_ot_env ...
1339 : ! **************************************************************************************************
1340 106920 : SUBROUTINE qs_ot_get_derivative_diag(matrix_hc, matrix_x, matrix_sx, matrix_gx, &
1341 : qs_ot_env)
1342 :
1343 : TYPE(dbcsr_type), POINTER :: matrix_hc, matrix_x, matrix_sx, matrix_gx
1344 : TYPE(qs_ot_type) :: qs_ot_env
1345 :
1346 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_ot_get_derivative_diag'
1347 : REAL(KIND=dp), PARAMETER :: rone = 1.0_dp, rzero = 0.0_dp
1348 :
1349 : INTEGER :: handle, k, n
1350 : TYPE(dbcsr_distribution_type) :: dist
1351 :
1352 35640 : CALL timeset(routineN, handle)
1353 :
1354 35640 : CALL dbcsr_get_info(matrix_x, nfullrows_total=n, nfullcols_total=k)
1355 :
1356 : ! go for the derivative now
1357 : ! this de/dc*(dX/dx)*sinp
1358 35640 : CALL dbcsr_multiply('N', 'N', rone, matrix_hc, qs_ot_env%matrix_sinp, rzero, matrix_gx)
1359 : ! overlap hc*x
1360 35640 : CALL dbcsr_multiply('T', 'N', rone, matrix_hc, matrix_x, rzero, qs_ot_env%matrix_buf2)
1361 : ! get it in the basis of the eigenvectors
1362 : CALL dbcsr_multiply('N', 'N', rone, qs_ot_env%matrix_buf2, qs_ot_env%matrix_r, &
1363 35640 : rzero, qs_ot_env%matrix_buf1)
1364 : CALL dbcsr_multiply('T', 'N', rone, qs_ot_env%matrix_r, qs_ot_env%matrix_buf1, &
1365 35640 : rzero, qs_ot_env%matrix_buf2)
1366 :
1367 : ! get the schur product of O_uv*B_uv
1368 : CALL dbcsr_hadamard_product(qs_ot_env%matrix_buf2, qs_ot_env%matrix_sinp_b, &
1369 35640 : qs_ot_env%matrix_buf3)
1370 :
1371 : ! overlap hc*c0
1372 : CALL dbcsr_multiply('T', 'N', rone, matrix_hc, qs_ot_env%matrix_c0, rzero, &
1373 35640 : qs_ot_env%matrix_buf2)
1374 : ! get it in the basis of the eigenvectors
1375 : CALL dbcsr_multiply('N', 'N', rone, qs_ot_env%matrix_buf2, qs_ot_env%matrix_r, &
1376 35640 : rzero, qs_ot_env%matrix_buf1)
1377 : CALL dbcsr_multiply('T', 'N', rone, qs_ot_env%matrix_r, qs_ot_env%matrix_buf1, &
1378 35640 : rzero, qs_ot_env%matrix_buf2)
1379 :
1380 : CALL dbcsr_hadamard_product(qs_ot_env%matrix_buf2, qs_ot_env%matrix_cosp_b, &
1381 35640 : qs_ot_env%matrix_buf4)
1382 :
1383 : ! add the two bs and compute b+b^T
1384 : CALL dbcsr_add(qs_ot_env%matrix_buf3, qs_ot_env%matrix_buf4, &
1385 35640 : alpha_scalar=rone, beta_scalar=rone)
1386 :
1387 : ! get the b in the eigenvector basis
1388 : CALL dbcsr_multiply('N', 'T', rone, qs_ot_env%matrix_buf3, qs_ot_env%matrix_r, &
1389 35640 : rzero, qs_ot_env%matrix_buf1)
1390 : CALL dbcsr_multiply('N', 'N', rone, qs_ot_env%matrix_r, qs_ot_env%matrix_buf1, &
1391 35640 : rzero, qs_ot_env%matrix_buf3)
1392 35640 : CALL dbcsr_get_info(qs_ot_env%matrix_buf3, distribution=dist)
1393 : CALL dbcsr_transposed(qs_ot_env%matrix_buf1, qs_ot_env%matrix_buf3, &
1394 : shallow_data_copy=.FALSE., use_distribution=dist, &
1395 35640 : transpose_distribution=.FALSE.)
1396 : CALL dbcsr_add(qs_ot_env%matrix_buf3, qs_ot_env%matrix_buf1, &
1397 35640 : alpha_scalar=rone, beta_scalar=rone)
1398 :
1399 : ! and add to the derivative
1400 : CALL dbcsr_multiply('N', 'N', rone, matrix_sx, qs_ot_env%matrix_buf3, &
1401 35640 : rone, matrix_gx)
1402 35640 : CALL timestop(handle)
1403 :
1404 35640 : END SUBROUTINE qs_ot_get_derivative_diag
1405 :
1406 : ! **************************************************************************************************
1407 : !> \brief compute the derivative of the taylor expansion below
1408 : !> \param matrix_hc ...
1409 : !> \param matrix_x ...
1410 : !> \param matrix_sx ...
1411 : !> \param matrix_gx ...
1412 : !> \param qs_ot_env ...
1413 : ! **************************************************************************************************
1414 176452 : SUBROUTINE qs_ot_get_derivative_taylor(matrix_hc, matrix_x, matrix_sx, matrix_gx, &
1415 : qs_ot_env)
1416 :
1417 : TYPE(dbcsr_type), POINTER :: matrix_hc, matrix_x, matrix_sx, matrix_gx
1418 : TYPE(qs_ot_type) :: qs_ot_env
1419 :
1420 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_ot_get_derivative_taylor'
1421 : REAL(KIND=dp), PARAMETER :: rone = 1.0_dp, rzero = 0.0_dp
1422 :
1423 : INTEGER :: handle, i, k, n
1424 : REAL(KIND=dp) :: cosfactor, sinfactor
1425 : TYPE(dbcsr_distribution_type) :: dist
1426 : TYPE(dbcsr_type), POINTER :: matrix_left, matrix_right
1427 :
1428 44113 : CALL timeset(routineN, handle)
1429 :
1430 44113 : CALL dbcsr_get_info(matrix_x, nfullrows_total=n, nfullcols_total=k)
1431 :
1432 : ! go for the derivative now
1433 : ! this de/dc*(dX/dx)*sinp i.e. zeroth order
1434 44113 : CALL dbcsr_multiply('N', 'N', rone, matrix_hc, qs_ot_env%matrix_sinp, rzero, matrix_gx)
1435 :
1436 44113 : IF (qs_ot_env%taylor_order <= 0) THEN
1437 0 : CALL timestop(handle)
1438 0 : RETURN
1439 : END IF
1440 :
1441 : ! we store the matrix that will multiply sx in matrix_r
1442 44113 : CALL dbcsr_set(qs_ot_env%matrix_r, rzero)
1443 :
1444 : ! just better names for matrix_cosp_b and matrix_sinp_b (they are buffer space here)
1445 44113 : matrix_left => qs_ot_env%matrix_cosp_b
1446 44113 : matrix_right => qs_ot_env%matrix_sinp_b
1447 :
1448 : ! overlap hc*x and add its transpose to matrix_left
1449 44113 : CALL dbcsr_multiply('T', 'N', rone, matrix_hc, matrix_x, rzero, matrix_left)
1450 44113 : CALL dbcsr_get_info(matrix_left, distribution=dist)
1451 : CALL dbcsr_transposed(qs_ot_env%matrix_buf1, matrix_left, &
1452 : shallow_data_copy=.FALSE., use_distribution=dist, &
1453 44113 : transpose_distribution=.FALSE.)
1454 : CALL dbcsr_add(matrix_left, qs_ot_env%matrix_buf1, &
1455 44113 : alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
1456 44113 : CALL dbcsr_copy(matrix_right, matrix_left)
1457 :
1458 : ! first order
1459 44113 : sinfactor = -1.0_dp/(2.0_dp*3.0_dp)
1460 : CALL dbcsr_add(qs_ot_env%matrix_r, matrix_left, &
1461 44113 : alpha_scalar=1.0_dp, beta_scalar=sinfactor)
1462 :
1463 : ! M
1464 : ! OM+MO
1465 : ! OOM+OMO+MOO
1466 : ! ...
1467 82991 : DO i = 2, qs_ot_env%taylor_order
1468 38878 : sinfactor = sinfactor*(-1.0_dp)/REAL(2*i*(2*i + 1), KIND=dp)
1469 38878 : CALL dbcsr_multiply('N', 'N', rone, qs_ot_env%matrix_p, matrix_left, rzero, qs_ot_env%matrix_buf1)
1470 38878 : CALL dbcsr_multiply('N', 'N', rone, matrix_right, qs_ot_env%matrix_p, rzero, matrix_left)
1471 38878 : CALL dbcsr_copy(matrix_right, matrix_left)
1472 : CALL dbcsr_add(matrix_left, qs_ot_env%matrix_buf1, &
1473 38878 : 1.0_dp, 1.0_dp)
1474 : CALL dbcsr_add(qs_ot_env%matrix_r, matrix_left, &
1475 82991 : alpha_scalar=1.0_dp, beta_scalar=sinfactor)
1476 : END DO
1477 :
1478 : ! overlap hc*c0 and add its transpose to matrix_left
1479 44113 : CALL dbcsr_multiply('T', 'N', rone, matrix_hc, qs_ot_env%matrix_c0, rzero, matrix_left)
1480 44113 : CALL dbcsr_get_info(matrix_left, distribution=dist)
1481 : CALL dbcsr_transposed(qs_ot_env%matrix_buf1, matrix_left, &
1482 : shallow_data_copy=.FALSE., use_distribution=dist, &
1483 44113 : transpose_distribution=.FALSE.)
1484 44113 : CALL dbcsr_add(matrix_left, qs_ot_env%matrix_buf1, 1.0_dp, 1.0_dp)
1485 44113 : CALL dbcsr_copy(matrix_right, matrix_left)
1486 :
1487 : ! first order
1488 44113 : cosfactor = -1.0_dp/(1.0_dp*2.0_dp)
1489 : CALL dbcsr_add(qs_ot_env%matrix_r, matrix_left, &
1490 44113 : alpha_scalar=1.0_dp, beta_scalar=cosfactor)
1491 :
1492 : ! M
1493 : ! OM+MO
1494 : ! OOM+OMO+MOO
1495 : ! ...
1496 82991 : DO i = 2, qs_ot_env%taylor_order
1497 38878 : cosfactor = cosfactor*(-1.0_dp)/REAL(2*i*(2*i - 1), KIND=dp)
1498 38878 : CALL dbcsr_multiply('N', 'N', rone, qs_ot_env%matrix_p, matrix_left, rzero, qs_ot_env%matrix_buf1)
1499 38878 : CALL dbcsr_multiply('N', 'N', rone, matrix_right, qs_ot_env%matrix_p, rzero, matrix_left)
1500 38878 : CALL dbcsr_copy(matrix_right, matrix_left)
1501 38878 : CALL dbcsr_add(matrix_left, qs_ot_env%matrix_buf1, 1.0_dp, 1.0_dp)
1502 : CALL dbcsr_add(qs_ot_env%matrix_r, matrix_left, &
1503 82991 : alpha_scalar=1.0_dp, beta_scalar=cosfactor)
1504 : END DO
1505 :
1506 : ! and add to the derivative
1507 44113 : CALL dbcsr_multiply('N', 'N', rone, matrix_sx, qs_ot_env%matrix_r, rone, matrix_gx)
1508 :
1509 44113 : CALL timestop(handle)
1510 :
1511 44113 : END SUBROUTINE qs_ot_get_derivative_taylor
1512 :
1513 : ! *************************************************************************************************
1514 : !> \brief computes a taylor expansion.
1515 : !> \param qs_ot_env ...
1516 : ! **************************************************************************************************
1517 96281 : SUBROUTINE qs_ot_p2m_taylor(qs_ot_env)
1518 : TYPE(qs_ot_type) :: qs_ot_env
1519 :
1520 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_ot_p2m_taylor'
1521 : REAL(KIND=dp), PARAMETER :: rone = 1.0_dp, rzero = 0.0_dp
1522 :
1523 : INTEGER :: handle, i, k
1524 : REAL(KIND=dp) :: cosfactor, sinfactor
1525 :
1526 59864 : CALL timeset(routineN, handle)
1527 :
1528 : ! zeroth order
1529 59864 : CALL dbcsr_set(qs_ot_env%matrix_cosp, rzero)
1530 59864 : CALL dbcsr_set(qs_ot_env%matrix_sinp, rzero)
1531 59864 : CALL dbcsr_add_on_diag(qs_ot_env%matrix_cosp, rone)
1532 59864 : CALL dbcsr_add_on_diag(qs_ot_env%matrix_sinp, rone)
1533 :
1534 59864 : IF (qs_ot_env%taylor_order <= 0) THEN
1535 0 : CALL timestop(handle)
1536 23447 : RETURN
1537 : END IF
1538 :
1539 : ! first order
1540 59864 : cosfactor = -1.0_dp/(1.0_dp*2.0_dp)
1541 59864 : sinfactor = -1.0_dp/(2.0_dp*3.0_dp)
1542 59864 : CALL dbcsr_add(qs_ot_env%matrix_cosp, qs_ot_env%matrix_p, alpha_scalar=1.0_dp, beta_scalar=cosfactor)
1543 59864 : CALL dbcsr_add(qs_ot_env%matrix_sinp, qs_ot_env%matrix_p, alpha_scalar=1.0_dp, beta_scalar=sinfactor)
1544 59864 : IF (qs_ot_env%taylor_order <= 1) THEN
1545 23447 : CALL timestop(handle)
1546 23447 : RETURN
1547 : END IF
1548 :
1549 : ! other orders
1550 36417 : CALL dbcsr_get_info(qs_ot_env%matrix_p, nfullrows_total=k)
1551 36417 : CALL dbcsr_copy(qs_ot_env%matrix_r, qs_ot_env%matrix_p)
1552 :
1553 92300 : DO i = 2, qs_ot_env%taylor_order
1554 : ! new power of p
1555 : CALL dbcsr_multiply('N', 'N', rone, qs_ot_env%matrix_p, qs_ot_env%matrix_r, &
1556 55883 : rzero, qs_ot_env%matrix_buf1)
1557 55883 : CALL dbcsr_copy(qs_ot_env%matrix_r, qs_ot_env%matrix_buf1)
1558 : ! add to the taylor expansion so far
1559 55883 : cosfactor = cosfactor*(-1.0_dp)/REAL(2*i*(2*i - 1), KIND=dp)
1560 55883 : sinfactor = sinfactor*(-1.0_dp)/REAL(2*i*(2*i + 1), KIND=dp)
1561 : CALL dbcsr_add(qs_ot_env%matrix_cosp, qs_ot_env%matrix_r, &
1562 55883 : alpha_scalar=1.0_dp, beta_scalar=cosfactor)
1563 : CALL dbcsr_add(qs_ot_env%matrix_sinp, qs_ot_env%matrix_r, &
1564 92300 : alpha_scalar=1.0_dp, beta_scalar=sinfactor)
1565 : END DO
1566 :
1567 36417 : CALL timestop(handle)
1568 :
1569 : END SUBROUTINE qs_ot_p2m_taylor
1570 :
1571 : ! **************************************************************************************************
1572 : !> \brief given p, computes - eigenstuff (matrix_r,evals)
1573 : !> - cos(p^0.5),p^(-0.5)*sin(p^0.5)
1574 : !> - the real b matrices, needed for the derivatives of these guys
1575 : !> cosp_b_ij=(1/(2pii) * int(cos(z^1/2)/((z-eval(i))*(z-eval(j))))
1576 : !> sinp_b_ij=(1/(2pii) * int(z^(-1/2)*sin(z^1/2)/((z-eval(i))*(z-eval(j))))
1577 : !> \param qs_ot_env ...
1578 : ! **************************************************************************************************
1579 52883 : SUBROUTINE qs_ot_p2m_diag(qs_ot_env)
1580 : TYPE(qs_ot_type) :: qs_ot_env
1581 :
1582 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_ot_p2m_diag'
1583 :
1584 : INTEGER :: handle, phase
1585 : LOGICAL :: cosine
1586 : TYPE(dbcsr_type), POINTER :: kernel, matrix, matrix_im
1587 :
1588 52883 : CALL timeset(routineN, handle)
1589 52883 : IF (qs_ot_env%has_complex_kpoint_state) THEN
1590 : CALL cp_dbcsr_heevd(matrix_re=qs_ot_env%matrix_p, matrix_im=qs_ot_env%matrix_p_im, &
1591 : eigenvectors_re=qs_ot_env%matrix_r, eigenvectors_im=qs_ot_env%matrix_r_im, &
1592 : eigenvalues=qs_ot_env%evals, para_env=qs_ot_env%para_env, &
1593 1574 : blacs_env=qs_ot_env%blacs_env)
1594 : ELSE
1595 51309 : CALL dbcsr_copy(qs_ot_env%matrix_buf1, qs_ot_env%matrix_p)
1596 : CALL cp_dbcsr_syevd(qs_ot_env%matrix_buf1, qs_ot_env%matrix_r, qs_ot_env%evals, &
1597 51309 : qs_ot_env%para_env, qs_ot_env%blacs_env)
1598 : END IF
1599 541256 : qs_ot_env%evals(:) = MAX(0.0_dp, qs_ot_env%evals(:))
1600 158649 : DO phase = 1, 2
1601 105766 : cosine = phase == 1
1602 105766 : IF (cosine) THEN
1603 52883 : matrix => qs_ot_env%matrix_cosp
1604 52883 : matrix_im => qs_ot_env%matrix_cosp_im
1605 52883 : kernel => qs_ot_env%matrix_cosp_b
1606 : ELSE
1607 52883 : matrix => qs_ot_env%matrix_sinp
1608 52883 : matrix_im => qs_ot_env%matrix_sinp_im
1609 52883 : kernel => qs_ot_env%matrix_sinp_b
1610 : END IF
1611 105766 : CALL qs_ot_spectral_values(qs_ot_env%evals, qs_ot_env%dum, cosine)
1612 105766 : CALL dbcsr_copy(qs_ot_env%matrix_buf1, qs_ot_env%matrix_r)
1613 105766 : CALL dbcsr_scale_by_vector(qs_ot_env%matrix_buf1, alpha=qs_ot_env%dum, side='right')
1614 105766 : IF (qs_ot_env%has_complex_kpoint_state) THEN
1615 3148 : CALL dbcsr_copy(qs_ot_env%matrix_buf1_im, qs_ot_env%matrix_r_im)
1616 3148 : CALL dbcsr_scale_by_vector(qs_ot_env%matrix_buf1_im, alpha=qs_ot_env%dum, side='right')
1617 : CALL qs_ot_complex_multiply('N', 'C', qs_ot_env%matrix_buf1, qs_ot_env%matrix_buf1_im, &
1618 : qs_ot_env%matrix_r, qs_ot_env%matrix_r_im, &
1619 3148 : matrix, matrix_im, qs_ot_env%matrix_buf2)
1620 : ELSE
1621 : ! Retain the original real-path operand order.
1622 102618 : CALL dbcsr_multiply('N', 'T', 1.0_dp, qs_ot_env%matrix_r, qs_ot_env%matrix_buf1, 0.0_dp, matrix)
1623 : END IF
1624 1082512 : qs_ot_env%dum(:) = SQRT(qs_ot_env%evals(:))
1625 105766 : CALL dbcsr_copy(kernel, matrix)
1626 158649 : CALL qs_ot_spectral_kernel_dbcsr(qs_ot_env%dum, kernel, cosine)
1627 : END DO
1628 52883 : CALL timestop(handle)
1629 52883 : END SUBROUTINE qs_ot_p2m_diag
1630 :
1631 : ! **************************************************************************************************
1632 : !> \brief Evaluate cos(sqrt(P)) or sinc(sqrt(P)) on the nonnegative eigenvalues of P.
1633 : !> \param evals nonnegative eigenvalues
1634 : !> \param values function values
1635 : !> \param cosine select cosine; otherwise sinc
1636 : ! **************************************************************************************************
1637 109450 : SUBROUTINE qs_ot_spectral_values(evals, values, cosine)
1638 : REAL(dp), INTENT(IN) :: evals(:)
1639 : REAL(dp), INTENT(OUT) :: values(:)
1640 : LOGICAL, INTENT(IN) :: cosine
1641 :
1642 : INTEGER :: i
1643 :
1644 109450 : IF (cosine) THEN
1645 552846 : DO i = 1, SIZE(evals)
1646 552846 : values(i) = COS(SQRT(evals(i)))
1647 : END DO
1648 : ELSE
1649 552846 : DO i = 1, SIZE(evals)
1650 552846 : values(i) = qs_ot_sinc(SQRT(evals(i)))
1651 : END DO
1652 : END IF
1653 109450 : END SUBROUTINE qs_ot_spectral_values
1654 :
1655 : ! **************************************************************************************************
1656 : !> \brief Build a block of the real STRICT Frechet kernel, including coincident eigenvalues.
1657 : !> \param roots_i square roots of the row eigenvalues
1658 : !> \param roots_j square roots of the column eigenvalues
1659 : !> \param block dense output block
1660 : !> \param cosine select cosine; otherwise sinc
1661 : ! **************************************************************************************************
1662 78090 : SUBROUTINE qs_ot_spectral_kernel(roots_i, roots_j, block, cosine)
1663 : REAL(dp), INTENT(IN) :: roots_i(:), roots_j(:)
1664 : REAL(dp), INTENT(OUT) :: block(:, :)
1665 : LOGICAL, INTENT(IN) :: cosine
1666 :
1667 : INTEGER :: i, j
1668 : REAL(dp) :: a, b
1669 :
1670 78090 : IF (cosine) THEN
1671 516794 : DO j = 1, SIZE(roots_j)
1672 10869414 : DO i = 1, SIZE(roots_i)
1673 10352620 : a = (roots_i(i) - roots_j(j))/2.0_dp
1674 10352620 : b = (roots_i(i) + roots_j(j))/2.0_dp
1675 10830369 : block(i, j) = -0.5_dp*qs_ot_sinc(a)*qs_ot_sinc(b)
1676 : END DO
1677 : END DO
1678 : ELSE
1679 516794 : DO j = 1, SIZE(roots_j)
1680 10869414 : DO i = 1, SIZE(roots_i)
1681 10830369 : block(i, j) = qs_ot_sincf(roots_i(i), roots_j(j))
1682 : END DO
1683 : END DO
1684 : END IF
1685 78090 : END SUBROUTINE qs_ot_spectral_kernel
1686 :
1687 : ! **************************************************************************************************
1688 : !> \brief Fill the existing blocks of a STRICT Frechet kernel with the shared dense formulas.
1689 : !> \param roots square roots of the eigenvalues
1690 : !> \param matrix distributed kernel
1691 : !> \param cosine select cosine; otherwise sinc
1692 : ! **************************************************************************************************
1693 105766 : SUBROUTINE qs_ot_spectral_kernel_dbcsr(roots, matrix, cosine)
1694 : REAL(dp), INTENT(IN) :: roots(:)
1695 : TYPE(dbcsr_type), INTENT(INOUT) :: matrix
1696 : LOGICAL, INTENT(IN) :: cosine
1697 :
1698 : INTEGER :: col, nc, nr, row
1699 : LOGICAL :: transposed
1700 105766 : REAL(dp), POINTER :: block(:, :)
1701 : TYPE(dbcsr_iterator_type) :: iter
1702 :
1703 105766 : CALL dbcsr_iterator_start(iter, matrix)
1704 180172 : DO WHILE (dbcsr_iterator_blocks_left(iter))
1705 : CALL dbcsr_iterator_next_block(iter, block=block, transposed=transposed, &
1706 74406 : row_size=nr, col_size=nc, row_offset=row, col_offset=col)
1707 180172 : IF (transposed) THEN
1708 0 : CALL qs_ot_spectral_kernel(roots(col:col + nc - 1), roots(row:row + nr - 1), block, cosine)
1709 : ELSE
1710 74406 : CALL qs_ot_spectral_kernel(roots(row:row + nr - 1), roots(col:col + nc - 1), block, cosine)
1711 : END IF
1712 : END DO
1713 105766 : CALL dbcsr_iterator_stop(iter)
1714 105766 : END SUBROUTINE qs_ot_spectral_kernel_dbcsr
1715 :
1716 : ! **************************************************************************************************
1717 : !> \brief computes sin(x)/x for all values of the argument
1718 : !> \param x ...
1719 : !> \return ...
1720 : ! **************************************************************************************************
1721 28960463 : FUNCTION qs_ot_sinc(x)
1722 :
1723 : REAL(KIND=dp), INTENT(IN) :: x
1724 : REAL(KIND=dp) :: qs_ot_sinc
1725 :
1726 : REAL(KIND=dp), PARAMETER :: q1 = 1.0_dp, q2 = -q1/(2.0_dp*3.0_dp), q3 = -q2/(4.0_dp*5.0_dp), &
1727 : q4 = -q3/(6.0_dp*7.0_dp), q5 = -q4/(8.0_dp*9.0_dp), q6 = -q5/(10.0_dp*11.0_dp), &
1728 : q7 = -q6/(12.0_dp*13.0_dp), q8 = -q7/(14.0_dp*15.0_dp), q9 = -q8/(16.0_dp*17.0_dp), &
1729 : q10 = -q9/(18.0_dp*19.0_dp)
1730 :
1731 : REAL(KIND=dp) :: y
1732 :
1733 28960463 : IF (ABS(x) > 0.5_dp) THEN
1734 8718467 : qs_ot_sinc = SIN(x)/x
1735 : ELSE
1736 20241996 : y = x*x
1737 20241996 : qs_ot_sinc = q1 + y*(q2 + y*(q3 + y*(q4 + y*(q5 + y*(q6 + y*(q7 + y*(q8 + y*(q9 + y*(q10)))))))))
1738 : END IF
1739 28960463 : END FUNCTION qs_ot_sinc
1740 :
1741 : ! **************************************************************************************************
1742 : !> \brief computes (1/(x^2-y^2))*(sinc(x)-sinc(y)) for all positive values of the arguments
1743 : !> \param xa ...
1744 : !> \param ya ...
1745 : !> \return ...
1746 : ! **************************************************************************************************
1747 10352620 : FUNCTION qs_ot_sincf(xa, ya)
1748 :
1749 : REAL(KIND=dp), INTENT(IN) :: xa, ya
1750 : REAL(KIND=dp) :: qs_ot_sincf
1751 :
1752 : INTEGER :: i
1753 : REAL(KIND=dp) :: a, b, rs, sf, x, xs, y, ybx, ybxs
1754 :
1755 : ! this is currently a limit of the routine, could be removed rather easily
1756 10352620 : IF (xa < 0) CPABORT("x is negative")
1757 10352620 : IF (ya < 0) CPABORT("y is negative")
1758 :
1759 10352620 : IF (xa < ya) THEN
1760 4949888 : x = ya
1761 4949888 : y = xa
1762 : ELSE
1763 5402732 : x = xa
1764 5402732 : y = ya
1765 : END IF
1766 :
1767 10352620 : IF (x < 0.5_dp) THEN ! use series, keeping in mind that x,y,x+y,x-y can all be zero
1768 :
1769 6474069 : qs_ot_sincf = 0.0_dp
1770 6474069 : IF (x > 0.0_dp) THEN
1771 6256886 : ybx = y/x
1772 : ELSE ! should be irrelevant !?
1773 : ybx = 0.0_dp
1774 : END IF
1775 :
1776 6474069 : sf = -1.0_dp/((1.0_dp + ybx)*6.0_dp)
1777 6474069 : rs = 1.0_dp
1778 6474069 : ybxs = ybx
1779 6474069 : xs = 1.0_dp
1780 :
1781 71214759 : DO i = 1, 10
1782 64740690 : qs_ot_sincf = qs_ot_sincf + sf*rs*xs*(1.0_dp + ybxs)
1783 64740690 : sf = -sf/(REAL((2*i + 2), dp)*REAL((2*i + 3), dp))
1784 64740690 : rs = rs + ybxs
1785 64740690 : ybxs = ybxs*ybx
1786 71214759 : xs = xs*x*x
1787 : END DO
1788 :
1789 : ELSE ! no series expansion
1790 3878551 : IF (x - y > 0.1_dp) THEN ! safe to use the normal form
1791 3591678 : qs_ot_sincf = (qs_ot_sinc(x) - qs_ot_sinc(y))/((x + y)*(x - y))
1792 : ELSE
1793 286873 : a = (x + y)/2.0_dp
1794 286873 : b = (x - y)/2.0_dp ! might be close to zero
1795 : ! y (=(a-b)) can not be close to zero since it is close to x>0.5
1796 286873 : qs_ot_sincf = (qs_ot_sinc(b)*COS(a) - qs_ot_sinc(a)*COS(b))/(2*x*y)
1797 : END IF
1798 : END IF
1799 :
1800 10352620 : END FUNCTION qs_ot_sincf
1801 :
1802 : END MODULE qs_ot
|
