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_distribution_type, dbcsr_filter, dbcsr_get_block_p, &
19 : dbcsr_get_info, dbcsr_get_occupation, dbcsr_init_p, dbcsr_iterator_blocks_left, &
20 : dbcsr_iterator_next_block, dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, &
21 : dbcsr_multiply, dbcsr_release, dbcsr_release_p, dbcsr_scale, dbcsr_set, dbcsr_transposed, &
22 : dbcsr_type
23 : USE cp_dbcsr_cholesky, ONLY: cp_dbcsr_cholesky_decompose,&
24 : cp_dbcsr_cholesky_invert,&
25 : cp_dbcsr_cholesky_restore
26 : USE cp_dbcsr_contrib, ONLY: dbcsr_add_on_diag,&
27 : dbcsr_frobenius_norm,&
28 : dbcsr_gershgorin_norm,&
29 : dbcsr_hadamard_product,&
30 : dbcsr_scale_by_vector
31 : USE cp_dbcsr_diag, ONLY: cp_dbcsr_heevd,&
32 : cp_dbcsr_syevd
33 : USE kinds, ONLY: dp
34 : USE message_passing, ONLY: mp_comm_type
35 : USE preconditioner, ONLY: apply_preconditioner
36 : USE preconditioner_types, ONLY: preconditioner_type
37 : USE qs_ot_types, ONLY: qs_ot_type
38 : #include "./base/base_uses.f90"
39 :
40 : IMPLICIT NONE
41 : PRIVATE
42 :
43 : PUBLIC :: qs_ot_get_p
44 : PUBLIC :: qs_ot_get_orbitals
45 : PUBLIC :: qs_ot_get_derivative
46 : PUBLIC :: qs_ot_get_orbitals_ref
47 : PUBLIC :: qs_ot_get_derivative_ref
48 : PUBLIC :: qs_ot_new_preconditioner
49 : PRIVATE :: qs_ot_p2m_diag
50 : PRIVATE :: qs_ot_sinc
51 : PRIVATE :: qs_ot_ref_poly
52 : PRIVATE :: qs_ot_ref_chol
53 : PRIVATE :: qs_ot_ref_lwdn
54 : PRIVATE :: qs_ot_ref_decide
55 : PRIVATE :: qs_ot_ref_update
56 : PRIVATE :: qs_ot_refine
57 : PRIVATE :: qs_ot_on_the_fly_localize
58 :
59 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_ot'
60 :
61 : CONTAINS
62 :
63 : ! **************************************************************************************************
64 : !> \brief gets ready to use the preconditioner/ or renew the preconditioner
65 : !> only keeps a pointer to the preconditioner.
66 : !> If you change the preconditioner, you have to call this routine
67 : !> you remain responsible of proper deallocate of your preconditioner
68 : !> (or you can reuse it on the next step of the computation)
69 : !> \param qs_ot_env ...
70 : !> \param preconditioner ...
71 : ! **************************************************************************************************
72 8244 : SUBROUTINE qs_ot_new_preconditioner(qs_ot_env, preconditioner)
73 : TYPE(qs_ot_type) :: qs_ot_env
74 : TYPE(preconditioner_type), POINTER :: preconditioner
75 :
76 : INTEGER :: ncoef
77 :
78 8244 : qs_ot_env%preconditioner => preconditioner
79 8244 : qs_ot_env%os_valid = .FALSE.
80 8244 : IF (.NOT. ASSOCIATED(qs_ot_env%matrix_psc0)) THEN
81 8244 : CALL dbcsr_init_p(qs_ot_env%matrix_psc0)
82 8244 : CALL dbcsr_copy(qs_ot_env%matrix_psc0, qs_ot_env%matrix_sc0, 'matrix_psc0')
83 : END IF
84 :
85 8244 : IF (.NOT. qs_ot_env%use_dx) THEN
86 4911 : qs_ot_env%use_dx = .TRUE.
87 4911 : CALL dbcsr_init_p(qs_ot_env%matrix_dx)
88 4911 : CALL dbcsr_copy(qs_ot_env%matrix_dx, qs_ot_env%matrix_gx, 'matrix_dx')
89 4911 : IF (qs_ot_env%settings%do_rotation) THEN
90 30 : CALL dbcsr_init_p(qs_ot_env%rot_mat_dx)
91 30 : CALL dbcsr_copy(qs_ot_env%rot_mat_dx, qs_ot_env%rot_mat_gx, 'rot_mat_dx')
92 : END IF
93 4911 : IF (qs_ot_env%settings%do_ener) THEN
94 0 : ncoef = SIZE(qs_ot_env%ener_gx)
95 0 : ALLOCATE (qs_ot_env%ener_dx(ncoef))
96 0 : qs_ot_env%ener_dx = 0.0_dp
97 : END IF
98 : END IF
99 :
100 8244 : END SUBROUTINE qs_ot_new_preconditioner
101 :
102 : ! **************************************************************************************************
103 : !> \brief ...
104 : !> \param qs_ot_env ...
105 : !> \param C_NEW ...
106 : !> \param SC ...
107 : !> \param G_OLD ...
108 : !> \param D ...
109 : ! **************************************************************************************************
110 420 : SUBROUTINE qs_ot_on_the_fly_localize(qs_ot_env, C_NEW, SC, G_OLD, D)
111 : !
112 : TYPE(qs_ot_type) :: qs_ot_env
113 : TYPE(dbcsr_type), POINTER :: C_NEW, SC, G_OLD, D
114 :
115 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_ot_on_the_fly_localize'
116 : INTEGER, PARAMETER :: taylor_order = 50
117 : REAL(KIND=dp), PARAMETER :: alpha = 0.1_dp, f2_eps = 0.01_dp
118 :
119 : INTEGER :: col, col_size, handle, i, k, n, p, row, &
120 : row_size
121 84 : REAL(dp), DIMENSION(:, :), POINTER :: block
122 : REAL(KIND=dp) :: expfactor, f2, norm_fro, norm_gct, tmp
123 : TYPE(dbcsr_distribution_type) :: dist
124 : TYPE(dbcsr_iterator_type) :: iter
125 : TYPE(dbcsr_type), POINTER :: C, Gp1, Gp2, GU, U
126 : TYPE(mp_comm_type) :: group
127 :
128 84 : CALL timeset(routineN, handle)
129 : !
130 : !
131 84 : CALL dbcsr_get_info(C_NEW, nfullrows_total=n, nfullcols_total=k)
132 : !
133 : ! C = C*expm(-G)
134 84 : GU => qs_ot_env%buf1_k_k_nosym ! a buffer
135 84 : U => qs_ot_env%buf2_k_k_nosym ! a buffer
136 84 : Gp1 => qs_ot_env%buf3_k_k_nosym ! a buffer
137 84 : Gp2 => qs_ot_env%buf4_k_k_nosym ! a buffer
138 84 : C => qs_ot_env%buf1_n_k ! a buffer
139 : !
140 : ! compute the derivative of the norm
141 : !-------------------------------------------------------------------
142 : ! (x^2+eps)^1/2
143 84 : f2 = 0.0_dp
144 84 : CALL dbcsr_copy(C, C_NEW)
145 84 : CALL dbcsr_iterator_start(iter, C)
146 182 : DO WHILE (dbcsr_iterator_blocks_left(iter))
147 98 : CALL dbcsr_iterator_next_block(iter, row, col, block, row_size=row_size, col_size=col_size)
148 686 : DO p = 1, col_size ! p
149 6258 : DO i = 1, row_size ! i
150 5656 : tmp = SQRT(block(i, p)**2 + f2_eps)
151 5656 : f2 = f2 + tmp
152 6160 : block(i, p) = block(i, p)/tmp
153 : END DO
154 : END DO
155 : END DO
156 84 : CALL dbcsr_iterator_stop(iter)
157 84 : CALL dbcsr_get_info(C, group=group)
158 84 : CALL group%sum(f2)
159 : !
160 : !
161 84 : CALL dbcsr_multiply('T', 'N', 1.0_dp, C, C_NEW, 0.0_dp, GU)
162 : !
163 : ! antisymetrize
164 84 : CALL dbcsr_get_info(GU, distribution=dist)
165 : CALL dbcsr_transposed(U, GU, shallow_data_copy=.FALSE., &
166 : use_distribution=dist, &
167 84 : transpose_distribution=.FALSE.)
168 84 : CALL dbcsr_add(GU, U, alpha_scalar=-0.5_dp, beta_scalar=0.5_dp)
169 : !-------------------------------------------------------------------
170 : !
171 84 : norm_fro = dbcsr_frobenius_norm(GU)
172 84 : norm_gct = dbcsr_gershgorin_norm(GU)
173 : !write(*,*) 'qs_ot_localize: ||P-I||_f=',norm_fro,' ||P-I||_GCT=',norm_gct
174 : !
175 : !kscale = CEILING(LOG(MIN(norm_fro,norm_gct))/LOG(2.0_dp))
176 : !scale = LOG(MIN(norm_fro,norm_gct))/LOG(2.0_dp)
177 : !write(*,*) 'qs_ot_localize: scale=',scale,' kscale=',kscale
178 : !
179 : ! rescale for steepest descent
180 84 : CALL dbcsr_scale(GU, -alpha)
181 : !
182 : ! compute unitary transform
183 : ! zeroth and first order
184 84 : expfactor = 1.0_dp
185 84 : CALL dbcsr_copy(U, GU)
186 84 : CALL dbcsr_scale(U, expfactor)
187 84 : CALL dbcsr_add_on_diag(U, 1.0_dp)
188 : ! other orders
189 84 : CALL dbcsr_copy(Gp1, GU)
190 520 : DO i = 2, taylor_order
191 : ! new power of G
192 520 : CALL dbcsr_multiply('N', 'N', 1.0_dp, GU, Gp1, 0.0_dp, Gp2)
193 520 : CALL dbcsr_copy(Gp1, Gp2)
194 : ! add to the taylor expansion so far
195 520 : expfactor = expfactor/REAL(i, KIND=dp)
196 520 : CALL dbcsr_add(U, Gp1, alpha_scalar=1.0_dp, beta_scalar=expfactor)
197 520 : norm_fro = dbcsr_frobenius_norm(Gp1)
198 : !write(*,*) 'Taylor expansion i=',i,' norm(X^i)/i!=',norm_fro*expfactor
199 520 : IF (norm_fro*expfactor < 1.0E-10_dp) EXIT
200 : END DO
201 : !
202 : ! rotate MOs
203 84 : CALL dbcsr_multiply('N', 'N', 1.0_dp, C_NEW, U, 0.0_dp, C)
204 84 : CALL dbcsr_copy(C_NEW, C)
205 : !
206 : ! rotate SC
207 84 : CALL dbcsr_multiply('N', 'N', 1.0_dp, SC, U, 0.0_dp, C)
208 84 : CALL dbcsr_copy(SC, C)
209 : !
210 : ! rotate D_i
211 84 : CALL dbcsr_multiply('N', 'N', 1.0_dp, D, U, 0.0_dp, C)
212 84 : CALL dbcsr_copy(D, C)
213 : !
214 : ! rotate G_i-1
215 84 : IF (ASSOCIATED(G_OLD)) THEN
216 84 : CALL dbcsr_multiply('N', 'N', 1.0_dp, G_OLD, U, 0.0_dp, C)
217 84 : CALL dbcsr_copy(G_OLD, C)
218 : END IF
219 : !
220 84 : CALL timestop(handle)
221 84 : END SUBROUTINE qs_ot_on_the_fly_localize
222 :
223 : ! **************************************************************************************************
224 : !> \brief ...
225 : !> \param qs_ot_env ...
226 : !> \param C_OLD ...
227 : !> \param C_TMP ...
228 : !> \param C_NEW ...
229 : !> \param P ...
230 : !> \param SC ...
231 : !> \param update ...
232 : ! **************************************************************************************************
233 1492 : SUBROUTINE qs_ot_ref_chol(qs_ot_env, C_OLD, C_TMP, C_NEW, P, SC, update)
234 : !
235 : TYPE(qs_ot_type) :: qs_ot_env
236 : TYPE(dbcsr_type) :: C_OLD, C_TMP, C_NEW, P, SC
237 : LOGICAL, INTENT(IN) :: update
238 :
239 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_ot_ref_chol'
240 :
241 : INTEGER :: handle, k, n
242 :
243 746 : CALL timeset(routineN, handle)
244 : !
245 746 : CALL dbcsr_get_info(C_NEW, nfullrows_total=n, nfullcols_total=k)
246 : !
247 : ! P = U'*U
248 746 : CALL cp_dbcsr_cholesky_decompose(P, k, qs_ot_env%para_env, qs_ot_env%blacs_env)
249 : !
250 : ! C_NEW = C_OLD*inv(U)
251 : CALL cp_dbcsr_cholesky_restore(C_OLD, k, P, C_NEW, op="SOLVE", pos="RIGHT", &
252 746 : transa="N", para_env=qs_ot_env%para_env, blacs_env=qs_ot_env%blacs_env)
253 : !
254 : ! Update SC if needed
255 746 : IF (update) THEN
256 : CALL cp_dbcsr_cholesky_restore(SC, k, P, C_TMP, op="SOLVE", pos="RIGHT", &
257 414 : transa="N", para_env=qs_ot_env%para_env, blacs_env=qs_ot_env%blacs_env)
258 414 : CALL dbcsr_copy(SC, C_TMP)
259 : END IF
260 : !
261 746 : CALL timestop(handle)
262 746 : END SUBROUTINE qs_ot_ref_chol
263 :
264 : ! **************************************************************************************************
265 : !> \brief ...
266 : !> \param qs_ot_env ...
267 : !> \param C_OLD ...
268 : !> \param C_TMP ...
269 : !> \param C_NEW ...
270 : !> \param P ...
271 : !> \param SC ...
272 : !> \param update ...
273 : ! **************************************************************************************************
274 308 : SUBROUTINE qs_ot_ref_lwdn(qs_ot_env, C_OLD, C_TMP, C_NEW, P, SC, update)
275 : !
276 : TYPE(qs_ot_type) :: qs_ot_env
277 : TYPE(dbcsr_type) :: C_OLD, C_TMP, C_NEW, P, SC
278 : LOGICAL, INTENT(IN) :: update
279 :
280 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_ot_ref_lwdn'
281 :
282 : INTEGER :: handle, i, k, n
283 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: eig, fun
284 : TYPE(dbcsr_type), POINTER :: V, W
285 :
286 308 : CALL timeset(routineN, handle)
287 : !
288 308 : CALL dbcsr_get_info(C_NEW, nfullrows_total=n, nfullcols_total=k)
289 : !
290 308 : V => qs_ot_env%buf1_k_k_nosym ! a buffer
291 308 : W => qs_ot_env%buf2_k_k_nosym ! a buffer
292 1232 : ALLOCATE (eig(k), fun(k))
293 : !
294 308 : CALL cp_dbcsr_syevd(P, V, eig, qs_ot_env%para_env, qs_ot_env%blacs_env)
295 : !
296 : ! compute the P^(-1/2)
297 1796 : DO i = 1, k
298 1488 : IF (eig(i) <= 0.0_dp) THEN
299 0 : CPABORT("P not positive definite")
300 : END IF
301 1796 : IF (eig(i) < 1.0E-8_dp) THEN
302 0 : fun(i) = 0.0_dp
303 : ELSE
304 1488 : fun(i) = 1.0_dp/SQRT(eig(i))
305 : END IF
306 : END DO
307 308 : CALL dbcsr_copy(W, V)
308 308 : CALL dbcsr_scale_by_vector(V, alpha=fun, side='right')
309 308 : CALL dbcsr_multiply('N', 'T', 1.0_dp, W, V, 0.0_dp, P)
310 : !
311 : ! Update C
312 308 : CALL dbcsr_multiply('N', 'N', 1.0_dp, C_OLD, P, 0.0_dp, C_NEW)
313 : !
314 : ! Update SC if needed
315 308 : IF (update) THEN
316 216 : CALL dbcsr_multiply('N', 'N', 1.0_dp, SC, P, 0.0_dp, C_TMP)
317 216 : CALL dbcsr_copy(SC, C_TMP)
318 : END IF
319 : !
320 308 : DEALLOCATE (eig, fun)
321 : !
322 308 : CALL timestop(handle)
323 308 : END SUBROUTINE qs_ot_ref_lwdn
324 :
325 : ! **************************************************************************************************
326 : !> \brief ...
327 : !> \param qs_ot_env ...
328 : !> \param C_OLD ...
329 : !> \param C_TMP ...
330 : !> \param C_NEW ...
331 : !> \param P ...
332 : !> \param SC ...
333 : !> \param norm_in ...
334 : !> \param update ...
335 : ! **************************************************************************************************
336 7104 : SUBROUTINE qs_ot_ref_poly(qs_ot_env, C_OLD, C_TMP, C_NEW, P, SC, norm_in, update)
337 : !
338 : TYPE(qs_ot_type) :: qs_ot_env
339 : TYPE(dbcsr_type), POINTER :: C_OLD, C_TMP, C_NEW, P
340 : TYPE(dbcsr_type) :: SC
341 : REAL(dp), INTENT(IN) :: norm_in
342 : LOGICAL, INTENT(IN) :: update
343 :
344 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_ot_ref_poly'
345 :
346 : INTEGER :: handle, irefine, k, n
347 : LOGICAL :: quick_exit
348 : REAL(dp) :: norm, norm_fro, norm_gct, occ_in, &
349 : occ_out, rescale
350 : TYPE(dbcsr_type), POINTER :: BUF1, BUF2, BUF_NOSYM, FT, FY
351 :
352 3552 : CALL timeset(routineN, handle)
353 : !
354 3552 : CALL dbcsr_get_info(C_NEW, nfullrows_total=n, nfullcols_total=k)
355 : !
356 3552 : BUF_NOSYM => qs_ot_env%buf1_k_k_nosym ! a buffer
357 3552 : BUF1 => qs_ot_env%buf1_k_k_sym ! a buffer
358 3552 : BUF2 => qs_ot_env%buf2_k_k_sym ! a buffer
359 3552 : FY => qs_ot_env%buf3_k_k_sym ! a buffer
360 3552 : FT => qs_ot_env%buf4_k_k_sym ! a buffer
361 : !
362 : ! initialize the norm (already computed in qs_ot_get_orbitals_ref)
363 3552 : norm = norm_in
364 : !
365 : ! can we do a quick exit?
366 3552 : quick_exit = .FALSE.
367 3552 : IF (norm < qs_ot_env%settings%eps_irac_quick_exit) quick_exit = .TRUE.
368 : !
369 : ! lets refine
370 3552 : rescale = 1.0_dp
371 3986 : DO irefine = 1, qs_ot_env%settings%max_irac
372 : !
373 : ! rescaling
374 3986 : IF (norm > 1.0_dp) THEN
375 12 : CALL dbcsr_scale(P, 1.0_dp/norm)
376 12 : rescale = rescale/SQRT(norm)
377 : END IF
378 : !
379 : ! get the refinement polynomial
380 : CALL qs_ot_refine(P, FY, BUF1, BUF2, qs_ot_env%settings%irac_degree, &
381 3986 : qs_ot_env%settings%eps_irac_filter_matrix)
382 : !
383 : ! collect the transformation
384 3986 : IF (irefine == 1) THEN
385 3552 : CALL dbcsr_copy(FT, FY, name='FT')
386 : ELSE
387 434 : CALL dbcsr_multiply('N', 'N', 1.0_dp, FT, FY, 0.0_dp, BUF1)
388 434 : IF (qs_ot_env%settings%eps_irac_filter_matrix > 0.0_dp) THEN
389 4 : occ_in = dbcsr_get_occupation(buf1)
390 4 : CALL dbcsr_filter(buf1, qs_ot_env%settings%eps_irac_filter_matrix)
391 4 : occ_out = dbcsr_get_occupation(buf1)
392 : END IF
393 434 : CALL dbcsr_copy(FT, BUF1, name='FT')
394 : END IF
395 : !
396 : ! quick exit if possible
397 3986 : IF (quick_exit) THEN
398 : EXIT
399 : END IF
400 : !
401 : ! P = FY^T * P * FY
402 1712 : CALL dbcsr_multiply('N', 'N', 1.0_dp, P, FY, 0.0_dp, BUF_NOSYM)
403 1712 : IF (qs_ot_env%settings%eps_irac_filter_matrix > 0.0_dp) THEN
404 8 : occ_in = dbcsr_get_occupation(buf_nosym)
405 8 : CALL dbcsr_filter(buf_nosym, qs_ot_env%settings%eps_irac_filter_matrix)
406 8 : occ_out = dbcsr_get_occupation(buf_nosym)
407 : END IF
408 1712 : CALL dbcsr_multiply('N', 'N', 1.0_dp, FY, BUF_NOSYM, 0.0_dp, P)
409 1712 : IF (qs_ot_env%settings%eps_irac_filter_matrix > 0.0_dp) THEN
410 8 : occ_in = dbcsr_get_occupation(p)
411 8 : CALL dbcsr_filter(p, qs_ot_env%settings%eps_irac_filter_matrix)
412 8 : occ_out = dbcsr_get_occupation(p)
413 : END IF
414 : !
415 : ! check ||P-1||_gct
416 1712 : CALL dbcsr_add_on_diag(P, -1.0_dp)
417 1712 : norm_fro = dbcsr_frobenius_norm(P)
418 1712 : norm_gct = dbcsr_gershgorin_norm(P)
419 1712 : CALL dbcsr_add_on_diag(P, 1.0_dp)
420 1712 : norm = MIN(norm_gct, norm_fro)
421 : !
422 : ! printing
423 : !
424 : ! blows up
425 1712 : IF (norm > 1.0E10_dp) THEN
426 : CALL cp_abort(__LOCATION__, &
427 : "Refinement blows up! "// &
428 : "We need you to improve the code, please post your input on "// &
429 0 : "the forum https://www.cp2k.org/")
430 : END IF
431 : !
432 : ! can we do a quick exit next step?
433 1712 : IF (norm < qs_ot_env%settings%eps_irac_quick_exit) quick_exit = .TRUE.
434 : !
435 : ! are we done?
436 3986 : IF (norm < qs_ot_env%settings%eps_irac) EXIT
437 : !
438 : END DO
439 : !
440 : ! C_NEW = C_NEW * FT * rescale
441 3552 : CALL dbcsr_multiply('N', 'N', rescale, C_OLD, FT, 0.0_dp, C_NEW)
442 3552 : IF (qs_ot_env%settings%eps_irac_filter_matrix > 0.0_dp) THEN
443 4 : occ_in = dbcsr_get_occupation(c_new)
444 4 : CALL dbcsr_filter(c_new, qs_ot_env%settings%eps_irac_filter_matrix)
445 4 : occ_out = dbcsr_get_occupation(c_new)
446 : END IF
447 : !
448 : ! update SC = SC * FY * rescale
449 3552 : IF (update) THEN
450 1412 : CALL dbcsr_multiply('N', 'N', rescale, SC, FT, 0.0_dp, C_TMP)
451 1412 : IF (qs_ot_env%settings%eps_irac_filter_matrix > 0.0_dp) THEN
452 4 : occ_in = dbcsr_get_occupation(c_tmp)
453 4 : CALL dbcsr_filter(c_tmp, qs_ot_env%settings%eps_irac_filter_matrix)
454 4 : occ_out = dbcsr_get_occupation(c_tmp)
455 : END IF
456 1412 : CALL dbcsr_copy(SC, C_TMP)
457 : END IF
458 : !
459 3552 : CALL timestop(handle)
460 3552 : END SUBROUTINE qs_ot_ref_poly
461 :
462 : ! **************************************************************************************************
463 : !> \brief ...
464 : !> \param qs_ot_env1 ...
465 : !> \return ...
466 : ! **************************************************************************************************
467 4606 : FUNCTION qs_ot_ref_update(qs_ot_env1) RESULT(update)
468 : !
469 : TYPE(qs_ot_type) :: qs_ot_env1
470 : LOGICAL :: update
471 :
472 4606 : update = .FALSE.
473 4174 : SELECT CASE (qs_ot_env1%settings%ot_method)
474 : CASE ("CG")
475 4174 : SELECT CASE (qs_ot_env1%settings%line_search_method)
476 : CASE ("2PNT")
477 4174 : IF (qs_ot_env1%line_search_count == 2) update = .TRUE.
478 : CASE DEFAULT
479 4174 : CPABORT("NYI")
480 : END SELECT
481 : CASE ("DIIS")
482 0 : update = .TRUE.
483 : CASE DEFAULT
484 4606 : CPABORT("NYI")
485 : END SELECT
486 4606 : END FUNCTION qs_ot_ref_update
487 :
488 : ! **************************************************************************************************
489 : !> \brief ...
490 : !> \param qs_ot_env1 ...
491 : !> \param norm_in ...
492 : !> \param ortho_irac ...
493 : ! **************************************************************************************************
494 4606 : SUBROUTINE qs_ot_ref_decide(qs_ot_env1, norm_in, ortho_irac)
495 : !
496 : TYPE(qs_ot_type) :: qs_ot_env1
497 : REAL(dp), INTENT(IN) :: norm_in
498 : CHARACTER(LEN=*), INTENT(INOUT) :: ortho_irac
499 :
500 4606 : ortho_irac = qs_ot_env1%settings%ortho_irac
501 4606 : IF (norm_in < qs_ot_env1%settings%eps_irac_switch) ortho_irac = "POLY"
502 4606 : END SUBROUTINE qs_ot_ref_decide
503 :
504 : ! **************************************************************************************************
505 : !> \brief ...
506 : !> \param matrix_c ...
507 : !> \param matrix_s ...
508 : !> \param matrix_x ...
509 : !> \param matrix_sx ...
510 : !> \param matrix_gx_old ...
511 : !> \param matrix_dx ...
512 : !> \param qs_ot_env ...
513 : !> \param qs_ot_env1 ...
514 : ! **************************************************************************************************
515 9212 : SUBROUTINE qs_ot_get_orbitals_ref(matrix_c, matrix_s, matrix_x, matrix_sx, &
516 : matrix_gx_old, matrix_dx, qs_ot_env, qs_ot_env1)
517 : !
518 : TYPE(dbcsr_type), POINTER :: matrix_c, matrix_s, matrix_x, matrix_sx, &
519 : matrix_gx_old, matrix_dx
520 : TYPE(qs_ot_type) :: qs_ot_env, qs_ot_env1
521 :
522 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_ot_get_orbitals_ref'
523 :
524 : CHARACTER(LEN=4) :: ortho_irac
525 : INTEGER :: handle, k, n
526 : LOGICAL :: on_the_fly_loc, update
527 : REAL(dp) :: norm, norm_fro, norm_gct, occ_in, occ_out
528 : TYPE(dbcsr_type), POINTER :: C_NEW, C_OLD, C_TMP, D, G_OLD, P, S, SC
529 :
530 4606 : CALL timeset(routineN, handle)
531 :
532 4606 : CALL dbcsr_get_info(matrix_c, nfullrows_total=n, nfullcols_total=k)
533 : !
534 4606 : C_NEW => matrix_c
535 4606 : C_OLD => matrix_x ! need to be carefully updated for the gradient !
536 4606 : SC => matrix_sx ! need to be carefully updated for the gradient !
537 4606 : G_OLD => matrix_gx_old ! need to be carefully updated for localization !
538 4606 : D => matrix_dx ! need to be carefully updated for localization !
539 4606 : S => matrix_s
540 :
541 4606 : P => qs_ot_env%p_k_k_sym ! a buffer
542 4606 : C_TMP => qs_ot_env%buf1_n_k ! a buffer
543 : !
544 : ! do we need to update C_OLD and SC?
545 4606 : update = qs_ot_ref_update(qs_ot_env1)
546 : !
547 : ! do we want to on the fly localize?
548 : ! for the moment this is set from the input,
549 : ! later we might want to localize every n-step or
550 : ! when the sparsity increases...
551 4606 : on_the_fly_loc = qs_ot_env1%settings%on_the_fly_loc
552 : !
553 : ! compute SC = S*C
554 4606 : IF (ASSOCIATED(S)) THEN
555 4606 : CALL dbcsr_multiply('N', 'N', 1.0_dp, S, C_OLD, 0.0_dp, SC)
556 4606 : IF (qs_ot_env1%settings%eps_irac_filter_matrix > 0.0_dp) THEN
557 4 : occ_in = dbcsr_get_occupation(sc)
558 4 : CALL dbcsr_filter(sc, qs_ot_env1%settings%eps_irac_filter_matrix)
559 4 : occ_out = dbcsr_get_occupation(sc)
560 : END IF
561 : ELSE
562 0 : CALL dbcsr_copy(SC, C_OLD)
563 : END IF
564 : !
565 : ! compute P = C'*SC
566 4606 : CALL dbcsr_multiply('T', 'N', 1.0_dp, C_OLD, SC, 0.0_dp, P)
567 4606 : IF (qs_ot_env1%settings%eps_irac_filter_matrix > 0.0_dp) THEN
568 4 : occ_in = dbcsr_get_occupation(p)
569 4 : CALL dbcsr_filter(p, qs_ot_env1%settings%eps_irac_filter_matrix)
570 4 : occ_out = dbcsr_get_occupation(p)
571 : END IF
572 : !
573 : ! check ||P-1||_f and ||P-1||_gct
574 4606 : CALL dbcsr_add_on_diag(P, -1.0_dp)
575 4606 : norm_fro = dbcsr_frobenius_norm(P)
576 4606 : norm_gct = dbcsr_gershgorin_norm(P)
577 4606 : CALL dbcsr_add_on_diag(P, 1.0_dp)
578 4606 : norm = MIN(norm_gct, norm_fro)
579 4606 : CALL qs_ot_ref_decide(qs_ot_env1, norm, ortho_irac)
580 : !
581 : ! select the orthogonality method
582 746 : SELECT CASE (ortho_irac)
583 : CASE ("CHOL")
584 746 : CALL qs_ot_ref_chol(qs_ot_env, C_OLD, C_TMP, C_NEW, P, SC, update)
585 : CASE ("LWDN")
586 308 : CALL qs_ot_ref_lwdn(qs_ot_env, C_OLD, C_TMP, C_NEW, P, SC, update)
587 : CASE ("POLY")
588 3552 : CALL qs_ot_ref_poly(qs_ot_env, C_OLD, C_TMP, C_NEW, P, SC, norm, update)
589 : CASE DEFAULT
590 4606 : CPABORT("Wrong argument")
591 : END SELECT
592 : !
593 : ! We update the C_i+1 and localization
594 4606 : IF (update) THEN
595 2042 : IF (on_the_fly_loc) THEN
596 84 : CALL qs_ot_on_the_fly_localize(qs_ot_env, C_NEW, SC, G_OLD, D)
597 : END IF
598 2042 : CALL dbcsr_copy(C_OLD, C_NEW)
599 : END IF
600 : !
601 4606 : CALL timestop(handle)
602 4606 : END SUBROUTINE qs_ot_get_orbitals_ref
603 :
604 : ! **************************************************************************************************
605 : !> \brief refinement polynomial of degree 2,3 and 4 (PRB 70, 193102 (2004))
606 : !> \param P ...
607 : !> \param FY ...
608 : !> \param P2 ...
609 : !> \param T ...
610 : !> \param irac_degree ...
611 : !> \param eps_irac_filter_matrix ...
612 : ! **************************************************************************************************
613 7972 : SUBROUTINE qs_ot_refine(P, FY, P2, T, irac_degree, eps_irac_filter_matrix)
614 : TYPE(dbcsr_type), INTENT(inout) :: P, FY, P2, T
615 : INTEGER, INTENT(in) :: irac_degree
616 : REAL(dp), INTENT(in) :: eps_irac_filter_matrix
617 :
618 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_ot_refine'
619 :
620 : INTEGER :: handle, k
621 : REAL(dp) :: occ_in, occ_out, r
622 :
623 3986 : CALL timeset(routineN, handle)
624 :
625 3986 : CALL dbcsr_get_info(P, nfullcols_total=k)
626 3986 : SELECT CASE (irac_degree)
627 : CASE (2)
628 : ! C_out = C_in * ( 15/8 * I - 10/8 * P + 3/8 * P^2)
629 0 : r = 3.0_dp/8.0_dp
630 0 : CALL dbcsr_multiply('N', 'N', r, P, P, 0.0_dp, FY)
631 0 : IF (eps_irac_filter_matrix > 0.0_dp) THEN
632 0 : occ_in = dbcsr_get_occupation(fy)
633 0 : CALL dbcsr_filter(fy, eps_irac_filter_matrix)
634 0 : occ_out = dbcsr_get_occupation(fy)
635 : END IF
636 0 : r = -10.0_dp/8.0_dp
637 0 : CALL dbcsr_add(FY, P, alpha_scalar=1.0_dp, beta_scalar=r)
638 0 : r = 15.0_dp/8.0_dp
639 0 : CALL dbcsr_add_on_diag(FY, alpha=r)
640 : CASE (3)
641 : ! C_out = C_in * ( 35/16 * I - 35/16 * P + 21/16 * P^2 - 5/16 P^3)
642 0 : CALL dbcsr_multiply('N', 'N', 1.0_dp, P, P, 0.0_dp, P2)
643 0 : IF (eps_irac_filter_matrix > 0.0_dp) THEN
644 0 : occ_in = dbcsr_get_occupation(p2)
645 0 : CALL dbcsr_filter(p2, eps_irac_filter_matrix)
646 0 : occ_out = dbcsr_get_occupation(p2)
647 : END IF
648 0 : r = -5.0_dp/16.0_dp
649 0 : CALL dbcsr_multiply('N', 'N', r, P2, P, 0.0_dp, FY)
650 0 : IF (eps_irac_filter_matrix > 0.0_dp) THEN
651 0 : occ_in = dbcsr_get_occupation(fy)
652 0 : CALL dbcsr_filter(fy, eps_irac_filter_matrix)
653 0 : occ_out = dbcsr_get_occupation(fy)
654 : END IF
655 0 : r = 21.0_dp/16.0_dp
656 0 : CALL dbcsr_add(FY, P2, alpha_scalar=1.0_dp, beta_scalar=r)
657 0 : r = -35.0_dp/16.0_dp
658 0 : CALL dbcsr_add(FY, P, alpha_scalar=1.0_dp, beta_scalar=r)
659 0 : r = 35.0_dp/16.0_dp
660 0 : CALL dbcsr_add_on_diag(FY, alpha=r)
661 : CASE (4)
662 : ! C_out = C_in * ( 315/128 * I - 420/128 * P + 378/128 * P^2 - 180/128 P^3 + 35/128 P^4 )
663 : ! = C_in * ( 315/128 * I - 420/128 * P + 378/128 * P^2 + ( - 180/128 * P + 35/128 * P^2 ) * P^2 )
664 3986 : CALL dbcsr_multiply('N', 'N', 1.0_dp, P, P, 0.0_dp, P2) ! P^2
665 3986 : IF (eps_irac_filter_matrix > 0.0_dp) THEN
666 8 : occ_in = dbcsr_get_occupation(p2)
667 8 : CALL dbcsr_filter(p2, eps_irac_filter_matrix)
668 8 : occ_out = dbcsr_get_occupation(p2)
669 : END IF
670 3986 : r = -180.0_dp/128.0_dp
671 3986 : CALL dbcsr_add(T, P, alpha_scalar=0.0_dp, beta_scalar=r) ! T=-180/128*P
672 3986 : r = 35.0_dp/128.0_dp
673 3986 : CALL dbcsr_add(T, P2, alpha_scalar=1.0_dp, beta_scalar=r) ! T=T+35/128*P^2
674 3986 : CALL dbcsr_multiply('N', 'N', 1.0_dp, T, P2, 0.0_dp, FY) ! Y=T*P^2
675 3986 : IF (eps_irac_filter_matrix > 0.0_dp) THEN
676 8 : occ_in = dbcsr_get_occupation(fy)
677 8 : CALL dbcsr_filter(fy, eps_irac_filter_matrix)
678 8 : occ_out = dbcsr_get_occupation(fy)
679 : END IF
680 3986 : r = 378.0_dp/128.0_dp
681 3986 : CALL dbcsr_add(FY, P2, alpha_scalar=1.0_dp, beta_scalar=r) ! Y=Y+378/128*P^2
682 3986 : r = -420.0_dp/128.0_dp
683 3986 : CALL dbcsr_add(FY, P, alpha_scalar=1.0_dp, beta_scalar=r) ! Y=Y-420/128*P
684 3986 : r = 315.0_dp/128.0_dp
685 3986 : CALL dbcsr_add_on_diag(FY, alpha=r) ! Y=Y+315/128*I
686 : CASE DEFAULT
687 3986 : CPABORT("This irac_order NYI")
688 : END SELECT
689 3986 : CALL timestop(handle)
690 3986 : END SUBROUTINE qs_ot_refine
691 :
692 : ! **************************************************************************************************
693 : !> \brief ...
694 : !> \param matrix_hc ...
695 : !> \param matrix_x ...
696 : !> \param matrix_sx ...
697 : !> \param matrix_gx ...
698 : !> \param qs_ot_env ...
699 : ! **************************************************************************************************
700 5992 : SUBROUTINE qs_ot_get_derivative_ref(matrix_hc, matrix_x, matrix_sx, matrix_gx, &
701 : qs_ot_env)
702 : TYPE(dbcsr_type), POINTER :: matrix_hc, matrix_x, matrix_sx, matrix_gx
703 : TYPE(qs_ot_type) :: qs_ot_env
704 :
705 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_ot_get_derivative_ref'
706 :
707 : INTEGER :: handle, k, n
708 : REAL(dp) :: occ_in, occ_out
709 : TYPE(dbcsr_type), POINTER :: C, CHC, G, G_dp, HC, SC
710 :
711 2996 : CALL timeset(routineN, handle)
712 :
713 2996 : CALL dbcsr_get_info(matrix_x, nfullrows_total=n, nfullcols_total=k)
714 : !
715 2996 : C => matrix_x ! NBsf*NOcc
716 2996 : SC => matrix_sx ! NBsf*NOcc need to be up2date
717 2996 : HC => matrix_hc ! NBsf*NOcc
718 2996 : G => matrix_gx ! NBsf*NOcc
719 2996 : CHC => qs_ot_env%buf1_k_k_sym ! buffer
720 2996 : G_dp => qs_ot_env%buf1_n_k_dp ! buffer dp
721 :
722 : ! C'*(H*C)
723 2996 : CALL dbcsr_multiply('T', 'N', 1.0_dp, C, HC, 0.0_dp, CHC)
724 2996 : IF (qs_ot_env%settings%eps_irac_filter_matrix > 0.0_dp) THEN
725 4 : occ_in = dbcsr_get_occupation(chc)
726 4 : CALL dbcsr_filter(chc, qs_ot_env%settings%eps_irac_filter_matrix)
727 4 : occ_out = dbcsr_get_occupation(chc)
728 : END IF
729 : ! (S*C)*(C'*H*C)
730 2996 : CALL dbcsr_multiply('N', 'N', 1.0_dp, SC, CHC, 0.0_dp, G)
731 2996 : IF (qs_ot_env%settings%eps_irac_filter_matrix > 0.0_dp) THEN
732 4 : occ_in = dbcsr_get_occupation(g)
733 4 : CALL dbcsr_filter(g, qs_ot_env%settings%eps_irac_filter_matrix)
734 4 : occ_out = dbcsr_get_occupation(g)
735 : END IF
736 : ! G = 2*(1-S*C*C')*H*C
737 2996 : CALL dbcsr_add(G, HC, alpha_scalar=-1.0_dp, beta_scalar=1.0_dp)
738 : !
739 2996 : CALL timestop(handle)
740 2996 : END SUBROUTINE qs_ot_get_derivative_ref
741 :
742 : ! **************************************************************************************************
743 : !> \brief computes p=x*S*x and the matrix functionals related matrices
744 : !> \param matrix_x ...
745 : !> \param matrix_sx ...
746 : !> \param qs_ot_env ...
747 : ! **************************************************************************************************
748 314679 : SUBROUTINE qs_ot_get_p(matrix_x, matrix_sx, qs_ot_env)
749 :
750 : TYPE(dbcsr_type), POINTER :: matrix_x, matrix_sx
751 : TYPE(qs_ot_type) :: qs_ot_env
752 :
753 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_ot_get_p'
754 : REAL(KIND=dp), PARAMETER :: rone = 1.0_dp, rzero = 0.0_dp
755 :
756 : INTEGER :: handle, k, max_iter, n
757 : LOGICAL :: converged
758 : REAL(KIND=dp) :: max_ev, min_ev, threshold
759 :
760 104893 : CALL timeset(routineN, handle)
761 :
762 104893 : CALL dbcsr_get_info(matrix_x, nfullrows_total=n, nfullcols_total=k)
763 :
764 : ! get the overlap
765 : CALL dbcsr_multiply('T', 'N', rone, matrix_x, matrix_sx, rzero, &
766 104893 : qs_ot_env%matrix_p)
767 :
768 : ! get an upper bound for the largest eigenvalue
769 : ! try using lancos first and fall back to gershgorin norm if it fails
770 104893 : max_iter = 30; threshold = 1.0E-03_dp
771 104893 : CALL arnoldi_extremal(qs_ot_env%matrix_p, max_ev, min_ev, converged, threshold, max_iter)
772 104893 : qs_ot_env%largest_eval_upper_bound = MAX(max_ev, ABS(min_ev))
773 :
774 104893 : IF (.NOT. converged) qs_ot_env%largest_eval_upper_bound = dbcsr_gershgorin_norm(qs_ot_env%matrix_p)
775 104893 : CALL decide_strategy(qs_ot_env)
776 104893 : IF (qs_ot_env%do_taylor) THEN
777 56314 : CALL qs_ot_p2m_taylor(qs_ot_env)
778 : ELSE
779 48579 : CALL qs_ot_p2m_diag(qs_ot_env)
780 : END IF
781 :
782 104893 : IF (qs_ot_env%settings%do_rotation) THEN
783 3246 : CALL qs_ot_generate_rotation(qs_ot_env)
784 : END IF
785 :
786 104893 : CALL timestop(handle)
787 :
788 104893 : END SUBROUTINE qs_ot_get_p
789 :
790 : ! **************************************************************************************************
791 : !> \brief computes the rotation matrix rot_mat_u that is associated to a given
792 : !> rot_mat_x using rot_mat_u=exp(rot_mat_x)
793 : !> \param qs_ot_env a valid qs_ot_env
794 : !> \par History
795 : !> 08.2004 created [Joost VandeVondele]
796 : !> 12.2024 Rewrite to use only real matrices [Ole Schuett]
797 : ! **************************************************************************************************
798 3246 : SUBROUTINE qs_ot_generate_rotation(qs_ot_env)
799 :
800 : TYPE(qs_ot_type) :: qs_ot_env
801 :
802 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_ot_generate_rotation'
803 :
804 : INTEGER :: handle, k
805 3246 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: exp_evals_im, exp_evals_re
806 : TYPE(dbcsr_type) :: buf_1, buf_2
807 :
808 3246 : CALL timeset(routineN, handle)
809 :
810 3246 : CALL dbcsr_get_info(qs_ot_env%rot_mat_x, nfullrows_total=k)
811 :
812 3246 : IF (k /= 0) THEN
813 : ! We want to compute: rot_mat_u = exp(i*rot_mat_x)
814 :
815 : ! Diagonalize: matrix = i*rot_mat_x.
816 : ! Note that matrix is imaginary and hermitian because rot_mat_x is real and anti-symmetric.
817 : CALL cp_dbcsr_heevd(matrix_im=qs_ot_env%rot_mat_x, & ! matrix_re omitted because it's zero
818 : eigenvectors_re=qs_ot_env%rot_mat_evec_re, &
819 : eigenvectors_im=qs_ot_env%rot_mat_evec_im, &
820 : eigenvalues=qs_ot_env%rot_mat_evals, &
821 : para_env=qs_ot_env%para_env, &
822 3194 : blacs_env=qs_ot_env%blacs_env)
823 :
824 : ! Compute: exp_evals = EXP(-i*rot_mat_evals)
825 12776 : ALLOCATE (exp_evals_re(k), exp_evals_im(k))
826 17074 : exp_evals_re(:) = COS(-qs_ot_env%rot_mat_evals(:))
827 17074 : exp_evals_im(:) = SIN(-qs_ot_env%rot_mat_evals(:))
828 :
829 : ! Compute: rot_mat_u = \sum_ij exp_evals_ij * |rot_mat_evec_i> <rot_mat_evec_j|
830 : ! Note that we need only two matrix multiplications because rot_mat_u is real.
831 3194 : CALL dbcsr_copy(buf_1, qs_ot_env%rot_mat_evec_re, name="buf_1")
832 3194 : CALL dbcsr_scale_by_vector(buf_1, alpha=exp_evals_re, side='right')
833 3194 : CALL dbcsr_copy(buf_2, qs_ot_env%rot_mat_evec_im, name="buf_2")
834 3194 : CALL dbcsr_scale_by_vector(buf_2, alpha=exp_evals_im, side='right')
835 3194 : CALL dbcsr_add(buf_1, buf_2, alpha_scalar=+1.0_dp, beta_scalar=-1.0_dp)
836 3194 : 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)
837 :
838 3194 : CALL dbcsr_copy(buf_1, qs_ot_env%rot_mat_evec_im)
839 3194 : CALL dbcsr_scale_by_vector(buf_1, alpha=exp_evals_re, side='right')
840 3194 : CALL dbcsr_copy(buf_2, qs_ot_env%rot_mat_evec_re)
841 3194 : CALL dbcsr_scale_by_vector(buf_2, alpha=exp_evals_im, side='right')
842 3194 : CALL dbcsr_add(buf_1, buf_2, alpha_scalar=+1.0_dp, beta_scalar=+1.0_dp)
843 3194 : 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)
844 :
845 : ! Clean up.
846 3194 : CALL dbcsr_release(buf_1)
847 3194 : CALL dbcsr_release(buf_2)
848 3194 : DEALLOCATE (exp_evals_re, exp_evals_im)
849 : END IF
850 :
851 3246 : CALL timestop(handle)
852 :
853 6492 : END SUBROUTINE qs_ot_generate_rotation
854 :
855 : ! **************************************************************************************************
856 : !> \brief computes the derivative fields with respect to rot_mat_x
857 : !> \param qs_ot_env valid qs_ot_env. In particular qs_ot_generate_rotation has to be called before
858 : !> and the rot_mat_dedu matrix has to be up to date
859 : !> \par History
860 : !> 08.2004 created [ Joost VandeVondele ]
861 : !> 12.2024 Rewrite to use only real matrices [Ole Schuett]
862 : ! **************************************************************************************************
863 3284 : SUBROUTINE qs_ot_rot_mat_derivative(qs_ot_env)
864 : TYPE(qs_ot_type) :: qs_ot_env
865 :
866 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_ot_rot_mat_derivative'
867 :
868 : INTEGER :: handle, i, j, k
869 : REAL(KIND=dp) :: e1, e2
870 : TYPE(dbcsr_type) :: outer_deriv_re, outer_deriv_im, mat_buf, &
871 : inner_deriv_re, inner_deriv_im
872 : TYPE(dbcsr_iterator_type) :: iter
873 1642 : INTEGER, DIMENSION(:), POINTER :: row_blk_offset, col_blk_offset
874 1642 : REAL(dp), DIMENSION(:, :), POINTER :: block_in_re, block_in_im, block_out_re, block_out_im
875 : INTEGER :: row, col
876 : LOGICAL :: found
877 : COMPLEX(dp) :: cval_in, cval_out
878 : TYPE(dbcsr_distribution_type) :: dist
879 :
880 1642 : CALL timeset(routineN, handle)
881 :
882 1642 : CALL dbcsr_get_info(qs_ot_env%rot_mat_u, nfullrows_total=k)
883 1642 : IF (k /= 0) THEN
884 1616 : CALL dbcsr_copy(qs_ot_env%matrix_buf1, qs_ot_env%rot_mat_dedu)
885 : ! now we get to the derivative wrt the antisymmetric matrix rot_mat_x
886 1616 : CALL dbcsr_copy(mat_buf, qs_ot_env%rot_mat_dedu, "mat_buf")
887 :
888 : ! inner_deriv_ij = <rot_mat_evec_i| rot_mat_dedu |rot_mat_evec_j>
889 1616 : CALL dbcsr_copy(inner_deriv_re, qs_ot_env%rot_mat_dedu, "inner_deriv_re") ! TODO just create
890 1616 : CALL dbcsr_copy(inner_deriv_im, qs_ot_env%rot_mat_dedu, "inner_deriv_im") ! TODO just create
891 :
892 1616 : 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)
893 1616 : CALL dbcsr_multiply('T', 'N', +1.0_dp, qs_ot_env%rot_mat_evec_im, mat_buf, 0.0_dp, inner_deriv_re)
894 1616 : CALL dbcsr_multiply('T', 'N', +1.0_dp, qs_ot_env%rot_mat_evec_re, mat_buf, 0.0_dp, inner_deriv_im)
895 :
896 1616 : 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)
897 1616 : CALL dbcsr_multiply('T', 'N', +1.0_dp, qs_ot_env%rot_mat_evec_re, mat_buf, 1.0_dp, inner_deriv_re)
898 1616 : CALL dbcsr_multiply('T', 'N', -1.0_dp, qs_ot_env%rot_mat_evec_im, mat_buf, 1.0_dp, inner_deriv_im)
899 :
900 : ! outer_deriv_ij = cint(eval_i, eval_j) * inner_deriv_ij
901 1616 : CALL dbcsr_copy(outer_deriv_re, qs_ot_env%rot_mat_dedu, "outer_deriv_re") ! TODO just create
902 1616 : CALL dbcsr_copy(outer_deriv_im, qs_ot_env%rot_mat_dedu, "outer_deriv_im") ! TODO just create
903 :
904 1616 : CALL dbcsr_get_info(qs_ot_env%rot_mat_dedu, row_blk_offset=row_blk_offset, col_blk_offset=col_blk_offset)
905 1616 : CALL dbcsr_iterator_start(iter, qs_ot_env%rot_mat_dedu)
906 2424 : DO WHILE (dbcsr_iterator_blocks_left(iter))
907 808 : CALL dbcsr_iterator_next_block(iter, row, col)
908 808 : CALL dbcsr_get_block_p(inner_deriv_re, row, col, block_in_re, found)
909 808 : CALL dbcsr_get_block_p(inner_deriv_im, row, col, block_in_im, found)
910 808 : CALL dbcsr_get_block_p(outer_deriv_re, row, col, block_out_re, found)
911 808 : CALL dbcsr_get_block_p(outer_deriv_im, row, col, block_out_im, found)
912 :
913 6081 : DO i = 1, SIZE(block_in_re, 1)
914 25820 : DO j = 1, SIZE(block_in_re, 2)
915 21355 : e1 = qs_ot_env%rot_mat_evals(row_blk_offset(row) + i - 1)
916 21355 : e2 = qs_ot_env%rot_mat_evals(col_blk_offset(col) + j - 1)
917 21355 : cval_in = CMPLX(block_in_re(i, j), block_in_im(i, j), dp)
918 21355 : cval_out = cval_in*cint(e1, e2)
919 21355 : block_out_re(i, j) = REAL(cval_out)
920 25012 : block_out_im(i, j) = AIMAG(cval_out)
921 : END DO
922 : END DO
923 : END DO
924 1616 : CALL dbcsr_iterator_stop(iter)
925 1616 : CALL dbcsr_release(inner_deriv_re)
926 1616 : CALL dbcsr_release(inner_deriv_im)
927 :
928 : ! Compute: matrix_buf1 = \sum_i outer_deriv_ij * |rot_mat_evec_i> <rot_mat_evec_j|
929 1616 : CALL dbcsr_multiply('N', 'N', +1.0_dp, qs_ot_env%rot_mat_evec_re, outer_deriv_re, 0.0_dp, mat_buf)
930 1616 : CALL dbcsr_multiply('N', 'N', -1.0_dp, qs_ot_env%rot_mat_evec_im, outer_deriv_im, 1.0_dp, mat_buf)
931 1616 : 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)
932 :
933 1616 : CALL dbcsr_multiply('N', 'N', +1.0_dp, qs_ot_env%rot_mat_evec_re, outer_deriv_im, 0.0_dp, mat_buf)
934 1616 : CALL dbcsr_multiply('N', 'N', +1.0_dp, qs_ot_env%rot_mat_evec_im, outer_deriv_re, 1.0_dp, mat_buf)
935 1616 : 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)
936 :
937 : ! Account for anti-symmetry of rot_mat_x.
938 1616 : CALL dbcsr_get_info(qs_ot_env%matrix_buf3, distribution=dist)
939 : CALL dbcsr_transposed(qs_ot_env%matrix_buf2, qs_ot_env%matrix_buf1, &
940 : shallow_data_copy=.FALSE., use_distribution=dist, &
941 1616 : transpose_distribution=.FALSE.)
942 :
943 : ! rot_mat_gx = matrix_buf1^T - matrix_buf1
944 1616 : CALL dbcsr_add(qs_ot_env%matrix_buf1, qs_ot_env%matrix_buf2, alpha_scalar=-1.0_dp, beta_scalar=+1.0_dp)
945 1616 : CALL dbcsr_copy(qs_ot_env%rot_mat_gx, qs_ot_env%matrix_buf1)
946 :
947 1616 : CALL dbcsr_release(mat_buf)
948 1616 : CALL dbcsr_release(outer_deriv_re)
949 1616 : CALL dbcsr_release(outer_deriv_im)
950 : END IF
951 3284 : CALL timestop(handle)
952 : CONTAINS
953 :
954 : ! **************************************************************************************************
955 : !> \brief ...
956 : !> \param e1 ...
957 : !> \param e2 ...
958 : !> \return ...
959 : ! **************************************************************************************************
960 21355 : FUNCTION cint(e1, e2)
961 : REAL(KIND=dp) :: e1, e2
962 : COMPLEX(KIND=dp) :: cint
963 :
964 : COMPLEX(KIND=dp) :: l1, l2, x
965 : INTEGER :: I
966 :
967 21355 : l1 = (0.0_dp, -1.0_dp)*e1
968 21355 : l2 = (0.0_dp, -1.0_dp)*e2
969 21355 : IF (ABS(l1 - l2) > 0.5_dp) THEN
970 994 : cint = (EXP(l1) - EXP(l2))/(l1 - l2)
971 : ELSE
972 : x = 1.0_dp
973 : cint = 0.0_dp
974 346137 : DO I = 1, 16
975 325776 : cint = cint + x
976 346137 : x = x*(l1 - l2)/REAL(I + 1, KIND=dp)
977 : END DO
978 20361 : cint = cint*EXP(l2)
979 : END IF
980 21355 : END FUNCTION cint
981 : END SUBROUTINE qs_ot_rot_mat_derivative
982 :
983 : ! **************************************************************************************************
984 : !> \brief decide strategy
985 : !> tries to decide if the taylor expansion of cos(sqrt(xsx)) converges rapidly enough
986 : !> to make a taylor expansion of the functions cos(sqrt(xsx)) and sin(sqrt(xsx))/sqrt(xsx)
987 : !> and their derivatives faster than their computation based on diagonalization since xsx can
988 : !> be very small, especially during dynamics, only a few terms might indeed be needed we find
989 : !> the necessary order N to have largest_eval_upper_bound**(N+1)/(2(N+1))! < eps_taylor
990 : !> \param qs_ot_env ...
991 : ! **************************************************************************************************
992 104893 : SUBROUTINE decide_strategy(qs_ot_env)
993 : TYPE(qs_ot_type) :: qs_ot_env
994 :
995 : INTEGER :: N
996 : REAL(KIND=dp) :: num_error
997 :
998 104893 : qs_ot_env%do_taylor = .FALSE.
999 104893 : N = 0
1000 104893 : num_error = qs_ot_env%largest_eval_upper_bound/(2.0_dp)
1001 447667 : DO WHILE (num_error > qs_ot_env%settings%eps_taylor .AND. N <= qs_ot_env%settings%max_taylor)
1002 342774 : N = N + 1
1003 386259 : num_error = num_error*qs_ot_env%largest_eval_upper_bound/REAL((2*N + 1)*(2*N + 2), KIND=dp)
1004 : END DO
1005 104893 : qs_ot_env%taylor_order = N
1006 104893 : IF (qs_ot_env%taylor_order <= qs_ot_env%settings%max_taylor) THEN
1007 56314 : qs_ot_env%do_taylor = .TRUE.
1008 : END IF
1009 :
1010 104893 : END SUBROUTINE decide_strategy
1011 :
1012 : ! **************************************************************************************************
1013 : !> \brief c=(c0*cos(p^0.5)+x*sin(p^0.5)*p^(-0.5)) x rot_mat_u
1014 : !> this assumes that x is already ortho to S*C0, and that p is x*S*x
1015 : !> rot_mat_u is an optional rotation matrix
1016 : !> \param matrix_c ...
1017 : !> \param matrix_x ...
1018 : !> \param qs_ot_env ...
1019 : ! **************************************************************************************************
1020 194550 : SUBROUTINE qs_ot_get_orbitals(matrix_c, matrix_x, qs_ot_env)
1021 :
1022 : TYPE(dbcsr_type), POINTER :: matrix_c, matrix_x
1023 : TYPE(qs_ot_type) :: qs_ot_env
1024 :
1025 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_ot_get_orbitals'
1026 : REAL(KIND=dp), PARAMETER :: rone = 1.0_dp, rzero = 0.0_dp
1027 :
1028 : INTEGER :: handle, k, n
1029 : TYPE(dbcsr_type), POINTER :: matrix_kk
1030 :
1031 97275 : CALL timeset(routineN, handle)
1032 :
1033 97275 : CALL dbcsr_get_info(matrix_x, nfullrows_total=n, nfullcols_total=k)
1034 :
1035 : ! rotate the multiplying matrices cosp and sinp instead of the result,
1036 : ! this should be cheaper for large basis sets
1037 97275 : IF (qs_ot_env%settings%do_rotation) THEN
1038 3032 : matrix_kk => qs_ot_env%matrix_buf1
1039 : CALL dbcsr_multiply('N', 'N', rone, qs_ot_env%matrix_cosp, &
1040 3032 : qs_ot_env%rot_mat_u, rzero, matrix_kk)
1041 : ELSE
1042 94243 : matrix_kk => qs_ot_env%matrix_cosp
1043 : END IF
1044 :
1045 : CALL dbcsr_multiply('N', 'N', rone, qs_ot_env%matrix_c0, matrix_kk, &
1046 97275 : rzero, matrix_c)
1047 :
1048 97275 : IF (qs_ot_env%settings%do_rotation) THEN
1049 3032 : matrix_kk => qs_ot_env%matrix_buf1
1050 : CALL dbcsr_multiply('N', 'N', rone, qs_ot_env%matrix_sinp, &
1051 3032 : qs_ot_env%rot_mat_u, rzero, matrix_kk)
1052 : ELSE
1053 94243 : matrix_kk => qs_ot_env%matrix_sinp
1054 : END IF
1055 : CALL dbcsr_multiply('N', 'N', rone, matrix_x, matrix_kk, &
1056 97275 : rone, matrix_c)
1057 :
1058 97275 : CALL timestop(handle)
1059 :
1060 97275 : END SUBROUTINE qs_ot_get_orbitals
1061 :
1062 : ! **************************************************************************************************
1063 : !> \brief this routines computes dE/dx=dx, with dx ortho to sc0
1064 : !> needs dE/dC=hc,C0,X,SX,p
1065 : !> if preconditioned it will not be the derivative, but the lagrangian multiplier
1066 : !> is changed so that P*dE/dx is the right derivative (i.e. in the allowed subspace)
1067 : !> \param matrix_hc ...
1068 : !> \param matrix_x ...
1069 : !> \param matrix_sx ...
1070 : !> \param matrix_gx ...
1071 : !> \param qs_ot_env ...
1072 : ! **************************************************************************************************
1073 224283 : SUBROUTINE qs_ot_get_derivative(matrix_hc, matrix_x, matrix_sx, matrix_gx, &
1074 : qs_ot_env)
1075 : TYPE(dbcsr_type), POINTER :: matrix_hc, matrix_x, matrix_sx, matrix_gx
1076 : TYPE(qs_ot_type) :: qs_ot_env
1077 :
1078 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_ot_get_derivative'
1079 : REAL(KIND=dp), PARAMETER :: rone = 1.0_dp, rzero = 0.0_dp
1080 :
1081 : INTEGER :: handle, k, n, ortho_k
1082 : TYPE(dbcsr_type), POINTER :: matrix_hc_local, matrix_target
1083 :
1084 74761 : CALL timeset(routineN, handle)
1085 :
1086 74761 : NULLIFY (matrix_hc_local)
1087 :
1088 74761 : CALL dbcsr_get_info(matrix_x, nfullrows_total=n, nfullcols_total=k)
1089 :
1090 : ! could in principle be taken inside qs_ot_get_derivative_* for increased efficiency
1091 : ! create a local rotated version of matrix_hc leaving matrix_hc untouched (needed
1092 : ! for lagrangian multipliers)
1093 74761 : IF (qs_ot_env%settings%do_rotation) THEN
1094 1642 : CALL dbcsr_copy(matrix_gx, matrix_hc) ! use gx as temporary
1095 1642 : CALL dbcsr_init_p(matrix_hc_local)
1096 1642 : CALL dbcsr_copy(matrix_hc_local, matrix_hc, name='matrix_hc_local')
1097 1642 : CALL dbcsr_set(matrix_hc_local, 0.0_dp)
1098 1642 : CALL dbcsr_multiply('N', 'T', rone, matrix_gx, qs_ot_env%rot_mat_u, rzero, matrix_hc_local)
1099 : ELSE
1100 73119 : matrix_hc_local => matrix_hc
1101 : END IF
1102 :
1103 74761 : IF (qs_ot_env%do_taylor) THEN
1104 41203 : CALL qs_ot_get_derivative_taylor(matrix_hc_local, matrix_x, matrix_sx, matrix_gx, qs_ot_env)
1105 : ELSE
1106 33558 : CALL qs_ot_get_derivative_diag(matrix_hc_local, matrix_x, matrix_sx, matrix_gx, qs_ot_env)
1107 : END IF
1108 :
1109 : ! and make it orthogonal
1110 74761 : CALL dbcsr_get_info(qs_ot_env%matrix_sc0, nfullcols_total=ortho_k)
1111 :
1112 74761 : IF (ASSOCIATED(qs_ot_env%preconditioner)) THEN
1113 64275 : matrix_target => qs_ot_env%matrix_psc0
1114 : ELSE
1115 10486 : matrix_target => qs_ot_env%matrix_sc0
1116 : END IF
1117 : ! first make the matrix os if not yet valid
1118 74761 : IF (.NOT. qs_ot_env%os_valid) THEN
1119 : ! this assumes that the preconditioner is a single matrix
1120 : ! that maps sc0 onto psc0
1121 :
1122 8150 : IF (ASSOCIATED(qs_ot_env%preconditioner)) THEN
1123 : CALL apply_preconditioner(qs_ot_env%preconditioner, qs_ot_env%matrix_sc0, &
1124 7168 : qs_ot_env%matrix_psc0)
1125 : END IF
1126 : CALL dbcsr_multiply('T', 'N', rone, &
1127 : qs_ot_env%matrix_sc0, matrix_target, &
1128 8150 : rzero, qs_ot_env%matrix_os)
1129 : CALL cp_dbcsr_cholesky_decompose(qs_ot_env%matrix_os, &
1130 8150 : para_env=qs_ot_env%para_env, blacs_env=qs_ot_env%blacs_env)
1131 : CALL cp_dbcsr_cholesky_invert(qs_ot_env%matrix_os, &
1132 : para_env=qs_ot_env%para_env, blacs_env=qs_ot_env%blacs_env, &
1133 8150 : uplo_to_full=.TRUE.)
1134 8150 : qs_ot_env%os_valid = .TRUE.
1135 : END IF
1136 : CALL dbcsr_multiply('T', 'N', rone, matrix_target, matrix_gx, &
1137 74761 : rzero, qs_ot_env%matrix_buf1_ortho)
1138 : CALL dbcsr_multiply('N', 'N', rone, qs_ot_env%matrix_os, &
1139 74761 : qs_ot_env%matrix_buf1_ortho, rzero, qs_ot_env%matrix_buf2_ortho)
1140 : CALL dbcsr_multiply('N', 'N', -rone, qs_ot_env%matrix_sc0, &
1141 74761 : qs_ot_env%matrix_buf2_ortho, rone, matrix_gx)
1142 : ! also treat the rot_mat gradient here
1143 74761 : IF (qs_ot_env%settings%do_rotation) THEN
1144 1642 : CALL qs_ot_rot_mat_derivative(qs_ot_env)
1145 : END IF
1146 :
1147 74761 : IF (qs_ot_env%settings%do_rotation) THEN
1148 1642 : CALL dbcsr_release_p(matrix_hc_local)
1149 : END IF
1150 :
1151 74761 : CALL timestop(handle)
1152 :
1153 74761 : END SUBROUTINE qs_ot_get_derivative
1154 :
1155 : ! **************************************************************************************************
1156 : !> \brief ...
1157 : !> \param matrix_hc ...
1158 : !> \param matrix_x ...
1159 : !> \param matrix_sx ...
1160 : !> \param matrix_gx ...
1161 : !> \param qs_ot_env ...
1162 : ! **************************************************************************************************
1163 100674 : SUBROUTINE qs_ot_get_derivative_diag(matrix_hc, matrix_x, matrix_sx, matrix_gx, &
1164 : qs_ot_env)
1165 :
1166 : TYPE(dbcsr_type), POINTER :: matrix_hc, matrix_x, matrix_sx, matrix_gx
1167 : TYPE(qs_ot_type) :: qs_ot_env
1168 :
1169 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_ot_get_derivative_diag'
1170 : REAL(KIND=dp), PARAMETER :: rone = 1.0_dp, rzero = 0.0_dp
1171 :
1172 : INTEGER :: handle, k, n
1173 : TYPE(dbcsr_distribution_type) :: dist
1174 :
1175 33558 : CALL timeset(routineN, handle)
1176 :
1177 33558 : CALL dbcsr_get_info(matrix_x, nfullrows_total=n, nfullcols_total=k)
1178 :
1179 : ! go for the derivative now
1180 : ! this de/dc*(dX/dx)*sinp
1181 33558 : CALL dbcsr_multiply('N', 'N', rone, matrix_hc, qs_ot_env%matrix_sinp, rzero, matrix_gx)
1182 : ! overlap hc*x
1183 33558 : CALL dbcsr_multiply('T', 'N', rone, matrix_hc, matrix_x, rzero, qs_ot_env%matrix_buf2)
1184 : ! get it in the basis of the eigenvectors
1185 : CALL dbcsr_multiply('N', 'N', rone, qs_ot_env%matrix_buf2, qs_ot_env%matrix_r, &
1186 33558 : rzero, qs_ot_env%matrix_buf1)
1187 : CALL dbcsr_multiply('T', 'N', rone, qs_ot_env%matrix_r, qs_ot_env%matrix_buf1, &
1188 33558 : rzero, qs_ot_env%matrix_buf2)
1189 :
1190 : ! get the schur product of O_uv*B_uv
1191 : CALL dbcsr_hadamard_product(qs_ot_env%matrix_buf2, qs_ot_env%matrix_sinp_b, &
1192 33558 : qs_ot_env%matrix_buf3)
1193 :
1194 : ! overlap hc*c0
1195 : CALL dbcsr_multiply('T', 'N', rone, matrix_hc, qs_ot_env%matrix_c0, rzero, &
1196 33558 : qs_ot_env%matrix_buf2)
1197 : ! get it in the basis of the eigenvectors
1198 : CALL dbcsr_multiply('N', 'N', rone, qs_ot_env%matrix_buf2, qs_ot_env%matrix_r, &
1199 33558 : rzero, qs_ot_env%matrix_buf1)
1200 : CALL dbcsr_multiply('T', 'N', rone, qs_ot_env%matrix_r, qs_ot_env%matrix_buf1, &
1201 33558 : rzero, qs_ot_env%matrix_buf2)
1202 :
1203 : CALL dbcsr_hadamard_product(qs_ot_env%matrix_buf2, qs_ot_env%matrix_cosp_b, &
1204 33558 : qs_ot_env%matrix_buf4)
1205 :
1206 : ! add the two bs and compute b+b^T
1207 : CALL dbcsr_add(qs_ot_env%matrix_buf3, qs_ot_env%matrix_buf4, &
1208 33558 : alpha_scalar=rone, beta_scalar=rone)
1209 :
1210 : ! get the b in the eigenvector basis
1211 : CALL dbcsr_multiply('N', 'T', rone, qs_ot_env%matrix_buf3, qs_ot_env%matrix_r, &
1212 33558 : rzero, qs_ot_env%matrix_buf1)
1213 : CALL dbcsr_multiply('N', 'N', rone, qs_ot_env%matrix_r, qs_ot_env%matrix_buf1, &
1214 33558 : rzero, qs_ot_env%matrix_buf3)
1215 33558 : CALL dbcsr_get_info(qs_ot_env%matrix_buf3, distribution=dist)
1216 : CALL dbcsr_transposed(qs_ot_env%matrix_buf1, qs_ot_env%matrix_buf3, &
1217 : shallow_data_copy=.FALSE., use_distribution=dist, &
1218 33558 : transpose_distribution=.FALSE.)
1219 : CALL dbcsr_add(qs_ot_env%matrix_buf3, qs_ot_env%matrix_buf1, &
1220 33558 : alpha_scalar=rone, beta_scalar=rone)
1221 :
1222 : ! and add to the derivative
1223 : CALL dbcsr_multiply('N', 'N', rone, matrix_sx, qs_ot_env%matrix_buf3, &
1224 33558 : rone, matrix_gx)
1225 33558 : CALL timestop(handle)
1226 :
1227 33558 : END SUBROUTINE qs_ot_get_derivative_diag
1228 :
1229 : ! **************************************************************************************************
1230 : !> \brief compute the derivative of the taylor expansion below
1231 : !> \param matrix_hc ...
1232 : !> \param matrix_x ...
1233 : !> \param matrix_sx ...
1234 : !> \param matrix_gx ...
1235 : !> \param qs_ot_env ...
1236 : ! **************************************************************************************************
1237 147928 : 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 41203 : CALL timeset(routineN, handle)
1252 :
1253 41203 : 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 41203 : CALL dbcsr_multiply('N', 'N', rone, matrix_hc, qs_ot_env%matrix_sinp, rzero, matrix_gx)
1258 :
1259 41203 : IF (qs_ot_env%taylor_order <= 0) THEN
1260 8442 : CALL timestop(handle)
1261 8442 : RETURN
1262 : END IF
1263 :
1264 : ! we store the matrix that will multiply sx in matrix_r
1265 32761 : 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 32761 : matrix_left => qs_ot_env%matrix_cosp_b
1269 32761 : matrix_right => qs_ot_env%matrix_sinp_b
1270 :
1271 : ! overlap hc*x and add its transpose to matrix_left
1272 32761 : CALL dbcsr_multiply('T', 'N', rone, matrix_hc, matrix_x, rzero, matrix_left)
1273 32761 : 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 32761 : transpose_distribution=.FALSE.)
1277 : CALL dbcsr_add(matrix_left, qs_ot_env%matrix_buf1, &
1278 32761 : alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
1279 32761 : CALL dbcsr_copy(matrix_right, matrix_left)
1280 :
1281 : ! first order
1282 32761 : sinfactor = -1.0_dp/(2.0_dp*3.0_dp)
1283 : CALL dbcsr_add(qs_ot_env%matrix_r, matrix_left, &
1284 32761 : alpha_scalar=1.0_dp, beta_scalar=sinfactor)
1285 :
1286 : ! M
1287 : ! OM+MO
1288 : ! OOM+OMO+MOO
1289 : ! ...
1290 68989 : DO i = 2, qs_ot_env%taylor_order
1291 36228 : sinfactor = sinfactor*(-1.0_dp)/REAL(2*i*(2*i + 1), KIND=dp)
1292 36228 : CALL dbcsr_multiply('N', 'N', rone, qs_ot_env%matrix_p, matrix_left, rzero, qs_ot_env%matrix_buf1)
1293 36228 : CALL dbcsr_multiply('N', 'N', rone, matrix_right, qs_ot_env%matrix_p, rzero, matrix_left)
1294 36228 : CALL dbcsr_copy(matrix_right, matrix_left)
1295 : CALL dbcsr_add(matrix_left, qs_ot_env%matrix_buf1, &
1296 36228 : 1.0_dp, 1.0_dp)
1297 : CALL dbcsr_add(qs_ot_env%matrix_r, matrix_left, &
1298 68989 : alpha_scalar=1.0_dp, beta_scalar=sinfactor)
1299 : END DO
1300 :
1301 : ! overlap hc*c0 and add its transpose to matrix_left
1302 32761 : CALL dbcsr_multiply('T', 'N', rone, matrix_hc, qs_ot_env%matrix_c0, rzero, matrix_left)
1303 32761 : 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 32761 : transpose_distribution=.FALSE.)
1307 32761 : CALL dbcsr_add(matrix_left, qs_ot_env%matrix_buf1, 1.0_dp, 1.0_dp)
1308 32761 : CALL dbcsr_copy(matrix_right, matrix_left)
1309 :
1310 : ! first order
1311 32761 : cosfactor = -1.0_dp/(1.0_dp*2.0_dp)
1312 : CALL dbcsr_add(qs_ot_env%matrix_r, matrix_left, &
1313 32761 : alpha_scalar=1.0_dp, beta_scalar=cosfactor)
1314 :
1315 : ! M
1316 : ! OM+MO
1317 : ! OOM+OMO+MOO
1318 : ! ...
1319 68989 : DO i = 2, qs_ot_env%taylor_order
1320 36228 : cosfactor = cosfactor*(-1.0_dp)/REAL(2*i*(2*i - 1), KIND=dp)
1321 36228 : CALL dbcsr_multiply('N', 'N', rone, qs_ot_env%matrix_p, matrix_left, rzero, qs_ot_env%matrix_buf1)
1322 36228 : CALL dbcsr_multiply('N', 'N', rone, matrix_right, qs_ot_env%matrix_p, rzero, matrix_left)
1323 36228 : CALL dbcsr_copy(matrix_right, matrix_left)
1324 36228 : 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 68989 : alpha_scalar=1.0_dp, beta_scalar=cosfactor)
1327 : END DO
1328 :
1329 : ! and add to the derivative
1330 32761 : CALL dbcsr_multiply('N', 'N', rone, matrix_sx, qs_ot_env%matrix_r, rone, matrix_gx)
1331 :
1332 32761 : CALL timestop(handle)
1333 :
1334 41203 : END SUBROUTINE qs_ot_get_derivative_taylor
1335 :
1336 : ! *************************************************************************************************
1337 : !> \brief computes a taylor expansion.
1338 : !> \param qs_ot_env ...
1339 : ! **************************************************************************************************
1340 91127 : SUBROUTINE qs_ot_p2m_taylor(qs_ot_env)
1341 : TYPE(qs_ot_type) :: qs_ot_env
1342 :
1343 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_ot_p2m_taylor'
1344 : REAL(KIND=dp), PARAMETER :: rone = 1.0_dp, rzero = 0.0_dp
1345 :
1346 : INTEGER :: handle, i, k
1347 : REAL(KIND=dp) :: cosfactor, sinfactor
1348 :
1349 56314 : CALL timeset(routineN, handle)
1350 :
1351 : ! zeroth order
1352 56314 : CALL dbcsr_set(qs_ot_env%matrix_cosp, rzero)
1353 56314 : CALL dbcsr_set(qs_ot_env%matrix_sinp, rzero)
1354 56314 : CALL dbcsr_add_on_diag(qs_ot_env%matrix_cosp, rone)
1355 56314 : CALL dbcsr_add_on_diag(qs_ot_env%matrix_sinp, rone)
1356 :
1357 56314 : IF (qs_ot_env%taylor_order <= 0) THEN
1358 9082 : CALL timestop(handle)
1359 21501 : RETURN
1360 : END IF
1361 :
1362 : ! first order
1363 47232 : cosfactor = -1.0_dp/(1.0_dp*2.0_dp)
1364 47232 : sinfactor = -1.0_dp/(2.0_dp*3.0_dp)
1365 47232 : CALL dbcsr_add(qs_ot_env%matrix_cosp, qs_ot_env%matrix_p, alpha_scalar=1.0_dp, beta_scalar=cosfactor)
1366 47232 : CALL dbcsr_add(qs_ot_env%matrix_sinp, qs_ot_env%matrix_p, alpha_scalar=1.0_dp, beta_scalar=sinfactor)
1367 47232 : IF (qs_ot_env%taylor_order <= 1) THEN
1368 12419 : CALL timestop(handle)
1369 12419 : RETURN
1370 : END IF
1371 :
1372 : ! other orders
1373 34813 : CALL dbcsr_get_info(qs_ot_env%matrix_p, nfullrows_total=k)
1374 34813 : CALL dbcsr_copy(qs_ot_env%matrix_r, qs_ot_env%matrix_p)
1375 :
1376 87460 : DO i = 2, qs_ot_env%taylor_order
1377 : ! new power of p
1378 : CALL dbcsr_multiply('N', 'N', rone, qs_ot_env%matrix_p, qs_ot_env%matrix_r, &
1379 52647 : rzero, qs_ot_env%matrix_buf1)
1380 52647 : CALL dbcsr_copy(qs_ot_env%matrix_r, qs_ot_env%matrix_buf1)
1381 : ! add to the taylor expansion so far
1382 52647 : cosfactor = cosfactor*(-1.0_dp)/REAL(2*i*(2*i - 1), KIND=dp)
1383 52647 : sinfactor = sinfactor*(-1.0_dp)/REAL(2*i*(2*i + 1), KIND=dp)
1384 : CALL dbcsr_add(qs_ot_env%matrix_cosp, qs_ot_env%matrix_r, &
1385 52647 : alpha_scalar=1.0_dp, beta_scalar=cosfactor)
1386 : CALL dbcsr_add(qs_ot_env%matrix_sinp, qs_ot_env%matrix_r, &
1387 87460 : alpha_scalar=1.0_dp, beta_scalar=sinfactor)
1388 : END DO
1389 :
1390 34813 : CALL timestop(handle)
1391 :
1392 : END SUBROUTINE qs_ot_p2m_taylor
1393 :
1394 : ! **************************************************************************************************
1395 : !> \brief 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 : !> \param qs_ot_env ...
1401 : ! **************************************************************************************************
1402 194316 : SUBROUTINE qs_ot_p2m_diag(qs_ot_env)
1403 :
1404 : TYPE(qs_ot_type) :: qs_ot_env
1405 :
1406 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_ot_p2m_diag'
1407 : REAL(KIND=dp), PARAMETER :: rone = 1.0_dp, rzero = 0.0_dp
1408 :
1409 : INTEGER :: col, col_offset, col_size, handle, i, j, &
1410 : k, row, row_offset, row_size
1411 48579 : REAL(dp), DIMENSION(:, :), POINTER :: block
1412 : REAL(KIND=dp) :: a, b
1413 : TYPE(dbcsr_iterator_type) :: iter
1414 :
1415 48579 : CALL timeset(routineN, handle)
1416 :
1417 48579 : CALL dbcsr_get_info(qs_ot_env%matrix_p, nfullrows_total=k)
1418 48579 : CALL dbcsr_copy(qs_ot_env%matrix_buf1, qs_ot_env%matrix_p)
1419 : CALL cp_dbcsr_syevd(qs_ot_env%matrix_buf1, qs_ot_env%matrix_r, qs_ot_env%evals, &
1420 48579 : qs_ot_env%para_env, qs_ot_env%blacs_env)
1421 518740 : DO i = 1, k
1422 518740 : qs_ot_env%evals(i) = MAX(0.0_dp, qs_ot_env%evals(i))
1423 : END DO
1424 :
1425 48579 : !$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(i) SHARED(k,qs_ot_env)
1426 : DO i = 1, k
1427 : qs_ot_env%dum(i) = COS(SQRT(qs_ot_env%evals(i)))
1428 : END DO
1429 48579 : CALL dbcsr_copy(qs_ot_env%matrix_buf1, qs_ot_env%matrix_r)
1430 48579 : CALL dbcsr_scale_by_vector(qs_ot_env%matrix_buf1, alpha=qs_ot_env%dum, side='right')
1431 : CALL dbcsr_multiply('N', 'T', rone, qs_ot_env%matrix_r, qs_ot_env%matrix_buf1, &
1432 48579 : rzero, qs_ot_env%matrix_cosp)
1433 :
1434 48579 : !$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(i) SHARED(k,qs_ot_env)
1435 : DO i = 1, k
1436 : qs_ot_env%dum(i) = qs_ot_sinc(SQRT(qs_ot_env%evals(i)))
1437 : END DO
1438 48579 : CALL dbcsr_copy(qs_ot_env%matrix_buf1, qs_ot_env%matrix_r)
1439 48579 : CALL dbcsr_scale_by_vector(qs_ot_env%matrix_buf1, alpha=qs_ot_env%dum, side='right')
1440 : CALL dbcsr_multiply('N', 'T', rone, qs_ot_env%matrix_r, qs_ot_env%matrix_buf1, &
1441 48579 : rzero, qs_ot_env%matrix_sinp)
1442 :
1443 48579 : CALL dbcsr_copy(qs_ot_env%matrix_cosp_b, qs_ot_env%matrix_cosp)
1444 48579 : CALL dbcsr_iterator_start(iter, qs_ot_env%matrix_cosp_b)
1445 83542 : DO WHILE (dbcsr_iterator_blocks_left(iter))
1446 : CALL dbcsr_iterator_next_block(iter, row, col, block, &
1447 : row_size=row_size, col_size=col_size, &
1448 34963 : row_offset=row_offset, col_offset=col_offset)
1449 540687 : DO j = 1, col_size
1450 10647570 : DO i = 1, row_size
1451 : a = (SQRT(qs_ot_env%evals(row_offset + i - 1)) &
1452 10155462 : - SQRT(qs_ot_env%evals(col_offset + j - 1)))/2.0_dp
1453 : b = (SQRT(qs_ot_env%evals(row_offset + i - 1)) &
1454 10155462 : + SQRT(qs_ot_env%evals(col_offset + j - 1)))/2.0_dp
1455 10612607 : block(i, j) = -0.5_dp*qs_ot_sinc(a)*qs_ot_sinc(b)
1456 : END DO
1457 : END DO
1458 : END DO
1459 48579 : CALL dbcsr_iterator_stop(iter)
1460 :
1461 48579 : CALL dbcsr_copy(qs_ot_env%matrix_sinp_b, qs_ot_env%matrix_sinp)
1462 48579 : CALL dbcsr_iterator_start(iter, qs_ot_env%matrix_sinp_b)
1463 83542 : DO WHILE (dbcsr_iterator_blocks_left(iter))
1464 : CALL dbcsr_iterator_next_block(iter, row, col, block, &
1465 : row_size=row_size, col_size=col_size, &
1466 34963 : row_offset=row_offset, col_offset=col_offset)
1467 540687 : DO j = 1, col_size
1468 10647570 : DO i = 1, row_size
1469 10155462 : a = SQRT(qs_ot_env%evals(row_offset + i - 1))
1470 10155462 : b = SQRT(qs_ot_env%evals(col_offset + j - 1))
1471 10612607 : block(i, j) = qs_ot_sincf(a, b)
1472 : END DO
1473 : END DO
1474 : END DO
1475 48579 : CALL dbcsr_iterator_stop(iter)
1476 :
1477 48579 : CALL timestop(handle)
1478 :
1479 48579 : END SUBROUTINE qs_ot_p2m_diag
1480 :
1481 : ! **************************************************************************************************
1482 : !> \brief computes sin(x)/x for all values of the argument
1483 : !> \param x ...
1484 : !> \return ...
1485 : ! **************************************************************************************************
1486 28486491 : FUNCTION qs_ot_sinc(x)
1487 :
1488 : REAL(KIND=dp), INTENT(IN) :: x
1489 : REAL(KIND=dp) :: qs_ot_sinc
1490 :
1491 : REAL(KIND=dp), PARAMETER :: q1 = 1.0_dp, q2 = -q1/(2.0_dp*3.0_dp), q3 = -q2/(4.0_dp*5.0_dp), &
1492 : q4 = -q3/(6.0_dp*7.0_dp), q5 = -q4/(8.0_dp*9.0_dp), q6 = -q5/(10.0_dp*11.0_dp), &
1493 : q7 = -q6/(12.0_dp*13.0_dp), q8 = -q7/(14.0_dp*15.0_dp), q9 = -q8/(16.0_dp*17.0_dp), &
1494 : q10 = -q9/(18.0_dp*19.0_dp)
1495 :
1496 : REAL(KIND=dp) :: y
1497 :
1498 28486491 : IF (ABS(x) > 0.5_dp) THEN
1499 8665293 : qs_ot_sinc = SIN(x)/x
1500 : ELSE
1501 19821198 : y = x*x
1502 19821198 : qs_ot_sinc = q1 + y*(q2 + y*(q3 + y*(q4 + y*(q5 + y*(q6 + y*(q7 + y*(q8 + y*(q9 + y*(q10)))))))))
1503 : END IF
1504 28486491 : END FUNCTION qs_ot_sinc
1505 :
1506 : ! **************************************************************************************************
1507 : !> \brief computes (1/(x^2-y^2))*(sinc(x)-sinc(y)) for all positive values of the arguments
1508 : !> \param xa ...
1509 : !> \param ya ...
1510 : !> \return ...
1511 : ! **************************************************************************************************
1512 10155462 : FUNCTION qs_ot_sincf(xa, ya)
1513 :
1514 : REAL(KIND=dp), INTENT(IN) :: xa, ya
1515 : REAL(KIND=dp) :: qs_ot_sincf
1516 :
1517 : INTEGER :: i
1518 : REAL(KIND=dp) :: a, b, rs, sf, x, xs, y, ybx, ybxs
1519 :
1520 : ! this is currently a limit of the routine, could be removed rather easily
1521 10155462 : IF (xa < 0) CPABORT("x is negative")
1522 10155462 : IF (ya < 0) CPABORT("y is negative")
1523 :
1524 10155462 : IF (xa < ya) THEN
1525 4869320 : x = ya
1526 4869320 : y = xa
1527 : ELSE
1528 5286142 : x = xa
1529 5286142 : y = ya
1530 : END IF
1531 :
1532 10155462 : IF (x < 0.5_dp) THEN ! use series, keeping in mind that x,y,x+y,x-y can all be zero
1533 :
1534 6302759 : qs_ot_sincf = 0.0_dp
1535 6302759 : IF (x > 0.0_dp) THEN
1536 6104901 : ybx = y/x
1537 : ELSE ! should be irrelevant !?
1538 : ybx = 0.0_dp
1539 : END IF
1540 :
1541 6302759 : sf = -1.0_dp/((1.0_dp + ybx)*6.0_dp)
1542 6302759 : rs = 1.0_dp
1543 6302759 : ybxs = ybx
1544 6302759 : xs = 1.0_dp
1545 :
1546 69330349 : DO i = 1, 10
1547 63027590 : qs_ot_sincf = qs_ot_sincf + sf*rs*xs*(1.0_dp + ybxs)
1548 63027590 : sf = -sf/(REAL((2*i + 2), dp)*REAL((2*i + 3), dp))
1549 63027590 : rs = rs + ybxs
1550 63027590 : ybxs = ybxs*ybx
1551 69330349 : xs = xs*x*x
1552 : END DO
1553 :
1554 : ELSE ! no series expansion
1555 3852703 : IF (x - y > 0.1_dp) THEN ! safe to use the normal form
1556 3567730 : qs_ot_sincf = (qs_ot_sinc(x) - qs_ot_sinc(y))/((x + y)*(x - y))
1557 : ELSE
1558 284973 : a = (x + y)/2.0_dp
1559 284973 : b = (x - y)/2.0_dp ! might be close to zero
1560 : ! y (=(a-b)) can not be close to zero since it is close to x>0.5
1561 284973 : qs_ot_sincf = (qs_ot_sinc(b)*COS(a) - qs_ot_sinc(a)*COS(b))/(2*x*y)
1562 : END IF
1563 : END IF
1564 :
1565 10155462 : END FUNCTION qs_ot_sincf
1566 :
1567 : END MODULE qs_ot
|
