Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief orbital transformations
10 : !> \par History
11 : !> Added Taylor expansion based computation of the matrix functions (01.2004)
12 : !> added additional rotation variables for non-equivalent occupied orbs (08.2004)
13 : !> \author Joost VandeVondele (06.2002)
14 : ! **************************************************************************************************
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 7486 : 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 7486 : qs_ot_env%preconditioner => preconditioner
79 7486 : qs_ot_env%os_valid = .FALSE.
80 7486 : IF (.NOT. ASSOCIATED(qs_ot_env%matrix_psc0)) THEN
81 7486 : CALL dbcsr_init_p(qs_ot_env%matrix_psc0)
82 7486 : CALL dbcsr_copy(qs_ot_env%matrix_psc0, qs_ot_env%matrix_sc0, 'matrix_psc0')
83 : END IF
84 :
85 7486 : IF (.NOT. qs_ot_env%use_dx) THEN
86 4153 : qs_ot_env%use_dx = .TRUE.
87 4153 : CALL dbcsr_init_p(qs_ot_env%matrix_dx)
88 4153 : CALL dbcsr_copy(qs_ot_env%matrix_dx, qs_ot_env%matrix_gx, 'matrix_dx')
89 4153 : 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 4153 : 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 7486 : 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 .LT. 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) .LE. 0.0_dp) &
299 0 : CPABORT("P not positive definite")
300 1796 : IF (eig(i) .LT. 1.0E-8_dp) THEN
301 0 : fun(i) = 0.0_dp
302 : ELSE
303 1488 : fun(i) = 1.0_dp/SQRT(eig(i))
304 : END IF
305 : END DO
306 308 : CALL dbcsr_copy(W, V)
307 308 : CALL dbcsr_scale_by_vector(V, alpha=fun, side='right')
308 308 : CALL dbcsr_multiply('N', 'T', 1.0_dp, W, V, 0.0_dp, P)
309 : !
310 : ! Update C
311 308 : CALL dbcsr_multiply('N', 'N', 1.0_dp, C_OLD, P, 0.0_dp, C_NEW)
312 : !
313 : ! Update SC if needed
314 308 : IF (update) THEN
315 216 : CALL dbcsr_multiply('N', 'N', 1.0_dp, SC, P, 0.0_dp, C_TMP)
316 216 : CALL dbcsr_copy(SC, C_TMP)
317 : END IF
318 : !
319 308 : DEALLOCATE (eig, fun)
320 : !
321 308 : CALL timestop(handle)
322 308 : END SUBROUTINE qs_ot_ref_lwdn
323 :
324 : ! **************************************************************************************************
325 : !> \brief ...
326 : !> \param qs_ot_env ...
327 : !> \param C_OLD ...
328 : !> \param C_TMP ...
329 : !> \param C_NEW ...
330 : !> \param P ...
331 : !> \param SC ...
332 : !> \param norm_in ...
333 : !> \param update ...
334 : ! **************************************************************************************************
335 7104 : SUBROUTINE qs_ot_ref_poly(qs_ot_env, C_OLD, C_TMP, C_NEW, P, SC, norm_in, update)
336 : !
337 : TYPE(qs_ot_type) :: qs_ot_env
338 : TYPE(dbcsr_type), POINTER :: C_OLD, C_TMP, C_NEW, P
339 : TYPE(dbcsr_type) :: SC
340 : REAL(dp), INTENT(IN) :: norm_in
341 : LOGICAL, INTENT(IN) :: update
342 :
343 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_ot_ref_poly'
344 :
345 : INTEGER :: handle, irefine, k, n
346 : LOGICAL :: quick_exit
347 : REAL(dp) :: norm, norm_fro, norm_gct, occ_in, &
348 : occ_out, rescale
349 : TYPE(dbcsr_type), POINTER :: BUF1, BUF2, BUF_NOSYM, FT, FY
350 :
351 3552 : CALL timeset(routineN, handle)
352 : !
353 3552 : CALL dbcsr_get_info(C_NEW, nfullrows_total=n, nfullcols_total=k)
354 : !
355 3552 : BUF_NOSYM => qs_ot_env%buf1_k_k_nosym ! a buffer
356 3552 : BUF1 => qs_ot_env%buf1_k_k_sym ! a buffer
357 3552 : BUF2 => qs_ot_env%buf2_k_k_sym ! a buffer
358 3552 : FY => qs_ot_env%buf3_k_k_sym ! a buffer
359 3552 : FT => qs_ot_env%buf4_k_k_sym ! a buffer
360 : !
361 : ! initialize the norm (already computed in qs_ot_get_orbitals_ref)
362 3552 : norm = norm_in
363 : !
364 : ! can we do a quick exit?
365 3552 : quick_exit = .FALSE.
366 3552 : IF (norm .LT. qs_ot_env%settings%eps_irac_quick_exit) quick_exit = .TRUE.
367 : !
368 : ! lets refine
369 3552 : rescale = 1.0_dp
370 3986 : DO irefine = 1, qs_ot_env%settings%max_irac
371 : !
372 : ! rescaling
373 3986 : IF (norm .GT. 1.0_dp) THEN
374 12 : CALL dbcsr_scale(P, 1.0_dp/norm)
375 12 : rescale = rescale/SQRT(norm)
376 : END IF
377 : !
378 : ! get the refinement polynomial
379 : CALL qs_ot_refine(P, FY, BUF1, BUF2, qs_ot_env%settings%irac_degree, &
380 3986 : qs_ot_env%settings%eps_irac_filter_matrix)
381 : !
382 : ! collect the transformation
383 3986 : IF (irefine .EQ. 1) THEN
384 3552 : CALL dbcsr_copy(FT, FY, name='FT')
385 : ELSE
386 434 : CALL dbcsr_multiply('N', 'N', 1.0_dp, FT, FY, 0.0_dp, BUF1)
387 434 : IF (qs_ot_env%settings%eps_irac_filter_matrix .GT. 0.0_dp) THEN
388 4 : occ_in = dbcsr_get_occupation(buf1)
389 4 : CALL dbcsr_filter(buf1, qs_ot_env%settings%eps_irac_filter_matrix)
390 4 : occ_out = dbcsr_get_occupation(buf1)
391 : END IF
392 434 : CALL dbcsr_copy(FT, BUF1, name='FT')
393 : END IF
394 : !
395 : ! quick exit if possible
396 3986 : IF (quick_exit) THEN
397 : EXIT
398 : END IF
399 : !
400 : ! P = FY^T * P * FY
401 1712 : CALL dbcsr_multiply('N', 'N', 1.0_dp, P, FY, 0.0_dp, BUF_NOSYM)
402 1712 : IF (qs_ot_env%settings%eps_irac_filter_matrix .GT. 0.0_dp) THEN
403 8 : occ_in = dbcsr_get_occupation(buf_nosym)
404 8 : CALL dbcsr_filter(buf_nosym, qs_ot_env%settings%eps_irac_filter_matrix)
405 8 : occ_out = dbcsr_get_occupation(buf_nosym)
406 : END IF
407 1712 : CALL dbcsr_multiply('N', 'N', 1.0_dp, FY, BUF_NOSYM, 0.0_dp, P)
408 1712 : IF (qs_ot_env%settings%eps_irac_filter_matrix .GT. 0.0_dp) THEN
409 8 : occ_in = dbcsr_get_occupation(p)
410 8 : CALL dbcsr_filter(p, qs_ot_env%settings%eps_irac_filter_matrix)
411 8 : occ_out = dbcsr_get_occupation(p)
412 : END IF
413 : !
414 : ! check ||P-1||_gct
415 1712 : CALL dbcsr_add_on_diag(P, -1.0_dp)
416 1712 : norm_fro = dbcsr_frobenius_norm(P)
417 1712 : norm_gct = dbcsr_gershgorin_norm(P)
418 1712 : CALL dbcsr_add_on_diag(P, 1.0_dp)
419 1712 : norm = MIN(norm_gct, norm_fro)
420 : !
421 : ! printing
422 : !
423 : ! blows up
424 1712 : IF (norm > 1.0E10_dp) THEN
425 : CALL cp_abort(__LOCATION__, &
426 : "Refinement blows up! "// &
427 : "We need you to improve the code, please post your input on "// &
428 0 : "the forum https://www.cp2k.org/")
429 : END IF
430 : !
431 : ! can we do a quick exit next step?
432 1712 : IF (norm .LT. qs_ot_env%settings%eps_irac_quick_exit) quick_exit = .TRUE.
433 : !
434 : ! are we done?
435 3986 : IF (norm .LT. qs_ot_env%settings%eps_irac) EXIT
436 : !
437 : END DO
438 : !
439 : ! C_NEW = C_NEW * FT * rescale
440 3552 : CALL dbcsr_multiply('N', 'N', rescale, C_OLD, FT, 0.0_dp, C_NEW)
441 3552 : IF (qs_ot_env%settings%eps_irac_filter_matrix .GT. 0.0_dp) THEN
442 4 : occ_in = dbcsr_get_occupation(c_new)
443 4 : CALL dbcsr_filter(c_new, qs_ot_env%settings%eps_irac_filter_matrix)
444 4 : occ_out = dbcsr_get_occupation(c_new)
445 : END IF
446 : !
447 : ! update SC = SC * FY * rescale
448 3552 : IF (update) THEN
449 1412 : CALL dbcsr_multiply('N', 'N', rescale, SC, FT, 0.0_dp, C_TMP)
450 1412 : IF (qs_ot_env%settings%eps_irac_filter_matrix .GT. 0.0_dp) THEN
451 4 : occ_in = dbcsr_get_occupation(c_tmp)
452 4 : CALL dbcsr_filter(c_tmp, qs_ot_env%settings%eps_irac_filter_matrix)
453 4 : occ_out = dbcsr_get_occupation(c_tmp)
454 : END IF
455 1412 : CALL dbcsr_copy(SC, C_TMP)
456 : END IF
457 : !
458 3552 : CALL timestop(handle)
459 3552 : END SUBROUTINE qs_ot_ref_poly
460 :
461 : ! **************************************************************************************************
462 : !> \brief ...
463 : !> \param qs_ot_env1 ...
464 : !> \return ...
465 : ! **************************************************************************************************
466 4606 : FUNCTION qs_ot_ref_update(qs_ot_env1) RESULT(update)
467 : !
468 : TYPE(qs_ot_type) :: qs_ot_env1
469 : LOGICAL :: update
470 :
471 4606 : update = .FALSE.
472 4174 : SELECT CASE (qs_ot_env1%settings%ot_method)
473 : CASE ("CG")
474 4174 : SELECT CASE (qs_ot_env1%settings%line_search_method)
475 : CASE ("2PNT")
476 4174 : IF (qs_ot_env1%line_search_count .EQ. 2) update = .TRUE.
477 : CASE DEFAULT
478 4174 : CPABORT("NYI")
479 : END SELECT
480 : CASE ("DIIS")
481 0 : update = .TRUE.
482 : CASE DEFAULT
483 4606 : CPABORT("NYI")
484 : END SELECT
485 4606 : END FUNCTION qs_ot_ref_update
486 :
487 : ! **************************************************************************************************
488 : !> \brief ...
489 : !> \param qs_ot_env1 ...
490 : !> \param norm_in ...
491 : !> \param ortho_irac ...
492 : ! **************************************************************************************************
493 4606 : SUBROUTINE qs_ot_ref_decide(qs_ot_env1, norm_in, ortho_irac)
494 : !
495 : TYPE(qs_ot_type) :: qs_ot_env1
496 : REAL(dp), INTENT(IN) :: norm_in
497 : CHARACTER(LEN=*), INTENT(INOUT) :: ortho_irac
498 :
499 4606 : ortho_irac = qs_ot_env1%settings%ortho_irac
500 4606 : IF (norm_in .LT. qs_ot_env1%settings%eps_irac_switch) ortho_irac = "POLY"
501 4606 : END SUBROUTINE qs_ot_ref_decide
502 :
503 : ! **************************************************************************************************
504 : !> \brief ...
505 : !> \param matrix_c ...
506 : !> \param matrix_s ...
507 : !> \param matrix_x ...
508 : !> \param matrix_sx ...
509 : !> \param matrix_gx_old ...
510 : !> \param matrix_dx ...
511 : !> \param qs_ot_env ...
512 : !> \param qs_ot_env1 ...
513 : ! **************************************************************************************************
514 9212 : SUBROUTINE qs_ot_get_orbitals_ref(matrix_c, matrix_s, matrix_x, matrix_sx, &
515 : matrix_gx_old, matrix_dx, qs_ot_env, qs_ot_env1)
516 : !
517 : TYPE(dbcsr_type), POINTER :: matrix_c, matrix_s, matrix_x, matrix_sx, &
518 : matrix_gx_old, matrix_dx
519 : TYPE(qs_ot_type) :: qs_ot_env, qs_ot_env1
520 :
521 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_ot_get_orbitals_ref'
522 :
523 : CHARACTER(LEN=4) :: ortho_irac
524 : INTEGER :: handle, k, n
525 : LOGICAL :: on_the_fly_loc, update
526 : REAL(dp) :: norm, norm_fro, norm_gct, occ_in, occ_out
527 : TYPE(dbcsr_type), POINTER :: C_NEW, C_OLD, C_TMP, D, G_OLD, P, S, SC
528 :
529 4606 : CALL timeset(routineN, handle)
530 :
531 4606 : CALL dbcsr_get_info(matrix_c, nfullrows_total=n, nfullcols_total=k)
532 : !
533 4606 : C_NEW => matrix_c
534 4606 : C_OLD => matrix_x ! need to be carefully updated for the gradient !
535 4606 : SC => matrix_sx ! need to be carefully updated for the gradient !
536 4606 : G_OLD => matrix_gx_old ! need to be carefully updated for localization !
537 4606 : D => matrix_dx ! need to be carefully updated for localization !
538 4606 : S => matrix_s
539 :
540 4606 : P => qs_ot_env%p_k_k_sym ! a buffer
541 4606 : C_TMP => qs_ot_env%buf1_n_k ! a buffer
542 : !
543 : ! do we need to update C_OLD and SC?
544 4606 : update = qs_ot_ref_update(qs_ot_env1)
545 : !
546 : ! do we want to on the fly localize?
547 : ! for the moment this is set from the input,
548 : ! later we might want to localize every n-step or
549 : ! when the sparsity increases...
550 4606 : on_the_fly_loc = qs_ot_env1%settings%on_the_fly_loc
551 : !
552 : ! compute SC = S*C
553 4606 : IF (ASSOCIATED(S)) THEN
554 4606 : CALL dbcsr_multiply('N', 'N', 1.0_dp, S, C_OLD, 0.0_dp, SC)
555 4606 : IF (qs_ot_env1%settings%eps_irac_filter_matrix .GT. 0.0_dp) THEN
556 4 : occ_in = dbcsr_get_occupation(sc)
557 4 : CALL dbcsr_filter(sc, qs_ot_env1%settings%eps_irac_filter_matrix)
558 4 : occ_out = dbcsr_get_occupation(sc)
559 : END IF
560 : ELSE
561 0 : CALL dbcsr_copy(SC, C_OLD)
562 : END IF
563 : !
564 : ! compute P = C'*SC
565 4606 : CALL dbcsr_multiply('T', 'N', 1.0_dp, C_OLD, SC, 0.0_dp, P)
566 4606 : IF (qs_ot_env1%settings%eps_irac_filter_matrix .GT. 0.0_dp) THEN
567 4 : occ_in = dbcsr_get_occupation(p)
568 4 : CALL dbcsr_filter(p, qs_ot_env1%settings%eps_irac_filter_matrix)
569 4 : occ_out = dbcsr_get_occupation(p)
570 : END IF
571 : !
572 : ! check ||P-1||_f and ||P-1||_gct
573 4606 : CALL dbcsr_add_on_diag(P, -1.0_dp)
574 4606 : norm_fro = dbcsr_frobenius_norm(P)
575 4606 : norm_gct = dbcsr_gershgorin_norm(P)
576 4606 : CALL dbcsr_add_on_diag(P, 1.0_dp)
577 4606 : norm = MIN(norm_gct, norm_fro)
578 4606 : CALL qs_ot_ref_decide(qs_ot_env1, norm, ortho_irac)
579 : !
580 : ! select the orthogonality method
581 746 : SELECT CASE (ortho_irac)
582 : CASE ("CHOL")
583 746 : CALL qs_ot_ref_chol(qs_ot_env, C_OLD, C_TMP, C_NEW, P, SC, update)
584 : CASE ("LWDN")
585 308 : CALL qs_ot_ref_lwdn(qs_ot_env, C_OLD, C_TMP, C_NEW, P, SC, update)
586 : CASE ("POLY")
587 3552 : CALL qs_ot_ref_poly(qs_ot_env, C_OLD, C_TMP, C_NEW, P, SC, norm, update)
588 : CASE DEFAULT
589 4606 : CPABORT("Wrong argument")
590 : END SELECT
591 : !
592 : ! We update the C_i+1 and localization
593 4606 : IF (update) THEN
594 2042 : IF (on_the_fly_loc) THEN
595 84 : CALL qs_ot_on_the_fly_localize(qs_ot_env, C_NEW, SC, G_OLD, D)
596 : END IF
597 2042 : CALL dbcsr_copy(C_OLD, C_NEW)
598 : END IF
599 : !
600 4606 : CALL timestop(handle)
601 4606 : END SUBROUTINE qs_ot_get_orbitals_ref
602 :
603 : ! **************************************************************************************************
604 : !> \brief refinement polynomial of degree 2,3 and 4 (PRB 70, 193102 (2004))
605 : !> \param P ...
606 : !> \param FY ...
607 : !> \param P2 ...
608 : !> \param T ...
609 : !> \param irac_degree ...
610 : !> \param eps_irac_filter_matrix ...
611 : ! **************************************************************************************************
612 7972 : SUBROUTINE qs_ot_refine(P, FY, P2, T, irac_degree, eps_irac_filter_matrix)
613 : TYPE(dbcsr_type), INTENT(inout) :: P, FY, P2, T
614 : INTEGER, INTENT(in) :: irac_degree
615 : REAL(dp), INTENT(in) :: eps_irac_filter_matrix
616 :
617 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_ot_refine'
618 :
619 : INTEGER :: handle, k
620 : REAL(dp) :: occ_in, occ_out, r
621 :
622 3986 : CALL timeset(routineN, handle)
623 :
624 3986 : CALL dbcsr_get_info(P, nfullcols_total=k)
625 3986 : SELECT CASE (irac_degree)
626 : CASE (2)
627 : ! C_out = C_in * ( 15/8 * I - 10/8 * P + 3/8 * P^2)
628 0 : r = 3.0_dp/8.0_dp
629 0 : CALL dbcsr_multiply('N', 'N', r, P, P, 0.0_dp, FY)
630 0 : IF (eps_irac_filter_matrix .GT. 0.0_dp) THEN
631 0 : occ_in = dbcsr_get_occupation(fy)
632 0 : CALL dbcsr_filter(fy, eps_irac_filter_matrix)
633 0 : occ_out = dbcsr_get_occupation(fy)
634 : END IF
635 0 : r = -10.0_dp/8.0_dp
636 0 : CALL dbcsr_add(FY, P, alpha_scalar=1.0_dp, beta_scalar=r)
637 0 : r = 15.0_dp/8.0_dp
638 0 : CALL dbcsr_add_on_diag(FY, alpha=r)
639 : CASE (3)
640 : ! C_out = C_in * ( 35/16 * I - 35/16 * P + 21/16 * P^2 - 5/16 P^3)
641 0 : CALL dbcsr_multiply('N', 'N', 1.0_dp, P, P, 0.0_dp, P2)
642 0 : IF (eps_irac_filter_matrix .GT. 0.0_dp) THEN
643 0 : occ_in = dbcsr_get_occupation(p2)
644 0 : CALL dbcsr_filter(p2, eps_irac_filter_matrix)
645 0 : occ_out = dbcsr_get_occupation(p2)
646 : END IF
647 0 : r = -5.0_dp/16.0_dp
648 0 : CALL dbcsr_multiply('N', 'N', r, P2, P, 0.0_dp, FY)
649 0 : IF (eps_irac_filter_matrix .GT. 0.0_dp) THEN
650 0 : occ_in = dbcsr_get_occupation(fy)
651 0 : CALL dbcsr_filter(fy, eps_irac_filter_matrix)
652 0 : occ_out = dbcsr_get_occupation(fy)
653 : END IF
654 0 : r = 21.0_dp/16.0_dp
655 0 : CALL dbcsr_add(FY, P2, alpha_scalar=1.0_dp, beta_scalar=r)
656 0 : r = -35.0_dp/16.0_dp
657 0 : CALL dbcsr_add(FY, P, alpha_scalar=1.0_dp, beta_scalar=r)
658 0 : r = 35.0_dp/16.0_dp
659 0 : CALL dbcsr_add_on_diag(FY, alpha=r)
660 : CASE (4)
661 : ! C_out = C_in * ( 315/128 * I - 420/128 * P + 378/128 * P^2 - 180/128 P^3 + 35/128 P^4 )
662 : ! = C_in * ( 315/128 * I - 420/128 * P + 378/128 * P^2 + ( - 180/128 * P + 35/128 * P^2 ) * P^2 )
663 3986 : CALL dbcsr_multiply('N', 'N', 1.0_dp, P, P, 0.0_dp, P2) ! P^2
664 3986 : IF (eps_irac_filter_matrix .GT. 0.0_dp) THEN
665 8 : occ_in = dbcsr_get_occupation(p2)
666 8 : CALL dbcsr_filter(p2, eps_irac_filter_matrix)
667 8 : occ_out = dbcsr_get_occupation(p2)
668 : END IF
669 3986 : r = -180.0_dp/128.0_dp
670 3986 : CALL dbcsr_add(T, P, alpha_scalar=0.0_dp, beta_scalar=r) ! T=-180/128*P
671 3986 : r = 35.0_dp/128.0_dp
672 3986 : CALL dbcsr_add(T, P2, alpha_scalar=1.0_dp, beta_scalar=r) ! T=T+35/128*P^2
673 3986 : CALL dbcsr_multiply('N', 'N', 1.0_dp, T, P2, 0.0_dp, FY) ! Y=T*P^2
674 3986 : IF (eps_irac_filter_matrix .GT. 0.0_dp) THEN
675 8 : occ_in = dbcsr_get_occupation(fy)
676 8 : CALL dbcsr_filter(fy, eps_irac_filter_matrix)
677 8 : occ_out = dbcsr_get_occupation(fy)
678 : END IF
679 3986 : r = 378.0_dp/128.0_dp
680 3986 : CALL dbcsr_add(FY, P2, alpha_scalar=1.0_dp, beta_scalar=r) ! Y=Y+378/128*P^2
681 3986 : r = -420.0_dp/128.0_dp
682 3986 : CALL dbcsr_add(FY, P, alpha_scalar=1.0_dp, beta_scalar=r) ! Y=Y-420/128*P
683 3986 : r = 315.0_dp/128.0_dp
684 3986 : CALL dbcsr_add_on_diag(FY, alpha=r) ! Y=Y+315/128*I
685 : CASE DEFAULT
686 3986 : CPABORT("This irac_order NYI")
687 : END SELECT
688 3986 : CALL timestop(handle)
689 3986 : END SUBROUTINE qs_ot_refine
690 :
691 : ! **************************************************************************************************
692 : !> \brief ...
693 : !> \param matrix_hc ...
694 : !> \param matrix_x ...
695 : !> \param matrix_sx ...
696 : !> \param matrix_gx ...
697 : !> \param qs_ot_env ...
698 : ! **************************************************************************************************
699 5992 : SUBROUTINE qs_ot_get_derivative_ref(matrix_hc, matrix_x, matrix_sx, matrix_gx, &
700 : qs_ot_env)
701 : TYPE(dbcsr_type), POINTER :: matrix_hc, matrix_x, matrix_sx, matrix_gx
702 : TYPE(qs_ot_type) :: qs_ot_env
703 :
704 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_ot_get_derivative_ref'
705 :
706 : INTEGER :: handle, k, n
707 : REAL(dp) :: occ_in, occ_out
708 : TYPE(dbcsr_type), POINTER :: C, CHC, G, G_dp, HC, SC
709 :
710 2996 : CALL timeset(routineN, handle)
711 :
712 2996 : CALL dbcsr_get_info(matrix_x, nfullrows_total=n, nfullcols_total=k)
713 : !
714 2996 : C => matrix_x ! NBsf*NOcc
715 2996 : SC => matrix_sx ! NBsf*NOcc need to be up2date
716 2996 : HC => matrix_hc ! NBsf*NOcc
717 2996 : G => matrix_gx ! NBsf*NOcc
718 2996 : CHC => qs_ot_env%buf1_k_k_sym ! buffer
719 2996 : G_dp => qs_ot_env%buf1_n_k_dp ! buffer dp
720 :
721 : ! C'*(H*C)
722 2996 : CALL dbcsr_multiply('T', 'N', 1.0_dp, C, HC, 0.0_dp, CHC)
723 2996 : IF (qs_ot_env%settings%eps_irac_filter_matrix .GT. 0.0_dp) THEN
724 4 : occ_in = dbcsr_get_occupation(chc)
725 4 : CALL dbcsr_filter(chc, qs_ot_env%settings%eps_irac_filter_matrix)
726 4 : occ_out = dbcsr_get_occupation(chc)
727 : END IF
728 : ! (S*C)*(C'*H*C)
729 2996 : CALL dbcsr_multiply('N', 'N', 1.0_dp, SC, CHC, 0.0_dp, G)
730 2996 : IF (qs_ot_env%settings%eps_irac_filter_matrix .GT. 0.0_dp) THEN
731 4 : occ_in = dbcsr_get_occupation(g)
732 4 : CALL dbcsr_filter(g, qs_ot_env%settings%eps_irac_filter_matrix)
733 4 : occ_out = dbcsr_get_occupation(g)
734 : END IF
735 : ! G = 2*(1-S*C*C')*H*C
736 2996 : CALL dbcsr_add(G, HC, alpha_scalar=-1.0_dp, beta_scalar=1.0_dp)
737 : !
738 2996 : CALL timestop(handle)
739 2996 : END SUBROUTINE qs_ot_get_derivative_ref
740 :
741 : ! **************************************************************************************************
742 : !> \brief computes p=x*S*x and the matrix functionals related matrices
743 : !> \param matrix_x ...
744 : !> \param matrix_sx ...
745 : !> \param qs_ot_env ...
746 : ! **************************************************************************************************
747 288915 : SUBROUTINE qs_ot_get_p(matrix_x, matrix_sx, qs_ot_env)
748 :
749 : TYPE(dbcsr_type), POINTER :: matrix_x, matrix_sx
750 : TYPE(qs_ot_type) :: qs_ot_env
751 :
752 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_ot_get_p'
753 : REAL(KIND=dp), PARAMETER :: rone = 1.0_dp, rzero = 0.0_dp
754 :
755 : INTEGER :: handle, k, max_iter, n
756 : LOGICAL :: converged
757 : REAL(KIND=dp) :: max_ev, min_ev, threshold
758 :
759 96305 : CALL timeset(routineN, handle)
760 :
761 96305 : CALL dbcsr_get_info(matrix_x, nfullrows_total=n, nfullcols_total=k)
762 :
763 : ! get the overlap
764 : CALL dbcsr_multiply('T', 'N', rone, matrix_x, matrix_sx, rzero, &
765 96305 : qs_ot_env%matrix_p)
766 :
767 : ! get an upper bound for the largest eigenvalue
768 : ! try using lancos first and fall back to gershgorin norm if it fails
769 96305 : max_iter = 30; threshold = 1.0E-03_dp
770 96305 : CALL arnoldi_extremal(qs_ot_env%matrix_p, max_ev, min_ev, converged, threshold, max_iter)
771 96305 : qs_ot_env%largest_eval_upper_bound = MAX(max_ev, ABS(min_ev))
772 :
773 96305 : IF (.NOT. converged) qs_ot_env%largest_eval_upper_bound = dbcsr_gershgorin_norm(qs_ot_env%matrix_p)
774 96305 : CALL decide_strategy(qs_ot_env)
775 96305 : IF (qs_ot_env%do_taylor) THEN
776 52174 : CALL qs_ot_p2m_taylor(qs_ot_env)
777 : ELSE
778 44131 : CALL qs_ot_p2m_diag(qs_ot_env)
779 : END IF
780 :
781 96305 : IF (qs_ot_env%settings%do_rotation) THEN
782 3206 : CALL qs_ot_generate_rotation(qs_ot_env)
783 : END IF
784 :
785 96305 : CALL timestop(handle)
786 :
787 96305 : END SUBROUTINE qs_ot_get_p
788 :
789 : ! **************************************************************************************************
790 : !> \brief computes the rotation matrix rot_mat_u that is associated to a given
791 : !> rot_mat_x using rot_mat_u=exp(rot_mat_x)
792 : !> \param qs_ot_env a valid qs_ot_env
793 : !> \par History
794 : !> 08.2004 created [Joost VandeVondele]
795 : !> 12.2024 Rewrite to use only real matrices [Ole Schuett]
796 : ! **************************************************************************************************
797 3206 : SUBROUTINE qs_ot_generate_rotation(qs_ot_env)
798 :
799 : TYPE(qs_ot_type) :: qs_ot_env
800 :
801 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_ot_generate_rotation'
802 :
803 : INTEGER :: handle, k
804 3206 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: exp_evals_im, exp_evals_re
805 : TYPE(dbcsr_type) :: buf_1, buf_2
806 :
807 3206 : CALL timeset(routineN, handle)
808 :
809 3206 : CALL dbcsr_get_info(qs_ot_env%rot_mat_x, nfullrows_total=k)
810 :
811 3206 : IF (k /= 0) THEN
812 : ! We want to compute: rot_mat_u = exp(i*rot_mat_x)
813 :
814 : ! Diagonalize: matrix = i*rot_mat_x.
815 : ! Note that matrix is imaginary and hermitian because rot_mat_x is real and anti-symmetric.
816 : CALL cp_dbcsr_heevd(matrix_im=qs_ot_env%rot_mat_x, & ! matrix_re omitted because it's zero
817 : eigenvectors_re=qs_ot_env%rot_mat_evec_re, &
818 : eigenvectors_im=qs_ot_env%rot_mat_evec_im, &
819 : eigenvalues=qs_ot_env%rot_mat_evals, &
820 : para_env=qs_ot_env%para_env, &
821 3154 : blacs_env=qs_ot_env%blacs_env)
822 :
823 : ! Compute: exp_evals = EXP(-i*rot_mat_evals)
824 12616 : ALLOCATE (exp_evals_re(k), exp_evals_im(k))
825 16874 : exp_evals_re(:) = COS(-qs_ot_env%rot_mat_evals(:))
826 16874 : exp_evals_im(:) = SIN(-qs_ot_env%rot_mat_evals(:))
827 :
828 : ! Compute: rot_mat_u = \sum_ij exp_evals_ij * |rot_mat_evec_i> <rot_mat_evec_j|
829 : ! Note that we need only two matrix multiplications because rot_mat_u is real.
830 3154 : CALL dbcsr_copy(buf_1, qs_ot_env%rot_mat_evec_re, name="buf_1")
831 3154 : CALL dbcsr_scale_by_vector(buf_1, alpha=exp_evals_re, side='right')
832 3154 : CALL dbcsr_copy(buf_2, qs_ot_env%rot_mat_evec_im, name="buf_2")
833 3154 : CALL dbcsr_scale_by_vector(buf_2, alpha=exp_evals_im, side='right')
834 3154 : CALL dbcsr_add(buf_1, buf_2, alpha_scalar=+1.0_dp, beta_scalar=-1.0_dp)
835 3154 : CALL dbcsr_multiply('N', 'T', 1.0_dp, buf_1, qs_ot_env%rot_mat_evec_re, 0.0_dp, qs_ot_env%rot_mat_u)
836 :
837 3154 : CALL dbcsr_copy(buf_1, qs_ot_env%rot_mat_evec_im)
838 3154 : CALL dbcsr_scale_by_vector(buf_1, alpha=exp_evals_re, side='right')
839 3154 : CALL dbcsr_copy(buf_2, qs_ot_env%rot_mat_evec_re)
840 3154 : CALL dbcsr_scale_by_vector(buf_2, alpha=exp_evals_im, side='right')
841 3154 : CALL dbcsr_add(buf_1, buf_2, alpha_scalar=+1.0_dp, beta_scalar=+1.0_dp)
842 3154 : CALL dbcsr_multiply('N', 'T', 1.0_dp, buf_1, qs_ot_env%rot_mat_evec_im, 1.0_dp, qs_ot_env%rot_mat_u)
843 :
844 : ! Clean up.
845 3154 : CALL dbcsr_release(buf_1)
846 3154 : CALL dbcsr_release(buf_2)
847 3154 : DEALLOCATE (exp_evals_re, exp_evals_im)
848 : END IF
849 :
850 3206 : CALL timestop(handle)
851 :
852 6412 : END SUBROUTINE qs_ot_generate_rotation
853 :
854 : ! **************************************************************************************************
855 : !> \brief computes the derivative fields with respect to rot_mat_x
856 : !> \param qs_ot_env valid qs_ot_env. In particular qs_ot_generate_rotation has to be called before
857 : !> and the rot_mat_dedu matrix has to be up to date
858 : !> \par History
859 : !> 08.2004 created [ Joost VandeVondele ]
860 : !> 12.2024 Rewrite to use only real matrices [Ole Schuett]
861 : ! **************************************************************************************************
862 3244 : SUBROUTINE qs_ot_rot_mat_derivative(qs_ot_env)
863 : TYPE(qs_ot_type) :: qs_ot_env
864 :
865 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_ot_rot_mat_derivative'
866 :
867 : INTEGER :: handle, i, j, k
868 : REAL(KIND=dp) :: e1, e2
869 : TYPE(dbcsr_type) :: outer_deriv_re, outer_deriv_im, mat_buf, &
870 : inner_deriv_re, inner_deriv_im
871 : TYPE(dbcsr_iterator_type) :: iter
872 1622 : INTEGER, DIMENSION(:), POINTER :: row_blk_offset, col_blk_offset
873 1622 : REAL(dp), DIMENSION(:, :), POINTER :: block_in_re, block_in_im, block_out_re, block_out_im
874 : INTEGER :: row, col
875 : LOGICAL :: found
876 : COMPLEX(dp) :: cval_in, cval_out
877 : TYPE(dbcsr_distribution_type) :: dist
878 :
879 1622 : CALL timeset(routineN, handle)
880 :
881 1622 : CALL dbcsr_get_info(qs_ot_env%rot_mat_u, nfullrows_total=k)
882 1622 : IF (k /= 0) THEN
883 1596 : CALL dbcsr_copy(qs_ot_env%matrix_buf1, qs_ot_env%rot_mat_dedu)
884 : ! now we get to the derivative wrt the antisymmetric matrix rot_mat_x
885 1596 : CALL dbcsr_copy(mat_buf, qs_ot_env%rot_mat_dedu, "mat_buf")
886 :
887 : ! inner_deriv_ij = <rot_mat_evec_i| rot_mat_dedu |rot_mat_evec_j>
888 1596 : CALL dbcsr_copy(inner_deriv_re, qs_ot_env%rot_mat_dedu, "inner_deriv_re") ! TODO just create
889 1596 : CALL dbcsr_copy(inner_deriv_im, qs_ot_env%rot_mat_dedu, "inner_deriv_im") ! TODO just create
890 :
891 1596 : CALL dbcsr_multiply('T', 'N', +1.0_dp, qs_ot_env%rot_mat_dedu, qs_ot_env%rot_mat_evec_im, 0.0_dp, mat_buf)
892 1596 : CALL dbcsr_multiply('T', 'N', +1.0_dp, qs_ot_env%rot_mat_evec_im, mat_buf, 0.0_dp, inner_deriv_re)
893 1596 : CALL dbcsr_multiply('T', 'N', +1.0_dp, qs_ot_env%rot_mat_evec_re, mat_buf, 0.0_dp, inner_deriv_im)
894 :
895 1596 : CALL dbcsr_multiply('T', 'N', +1.0_dp, qs_ot_env%rot_mat_dedu, qs_ot_env%rot_mat_evec_re, 0.0_dp, mat_buf)
896 1596 : CALL dbcsr_multiply('T', 'N', +1.0_dp, qs_ot_env%rot_mat_evec_re, mat_buf, 1.0_dp, inner_deriv_re)
897 1596 : CALL dbcsr_multiply('T', 'N', -1.0_dp, qs_ot_env%rot_mat_evec_im, mat_buf, 1.0_dp, inner_deriv_im)
898 :
899 : ! outer_deriv_ij = cint(eval_i, eval_j) * inner_deriv_ij
900 1596 : CALL dbcsr_copy(outer_deriv_re, qs_ot_env%rot_mat_dedu, "outer_deriv_re") ! TODO just create
901 1596 : CALL dbcsr_copy(outer_deriv_im, qs_ot_env%rot_mat_dedu, "outer_deriv_im") ! TODO just create
902 :
903 1596 : CALL dbcsr_get_info(qs_ot_env%rot_mat_dedu, row_blk_offset=row_blk_offset, col_blk_offset=col_blk_offset)
904 1596 : CALL dbcsr_iterator_start(iter, qs_ot_env%rot_mat_dedu)
905 2394 : DO WHILE (dbcsr_iterator_blocks_left(iter))
906 798 : CALL dbcsr_iterator_next_block(iter, row, col)
907 798 : CALL dbcsr_get_block_p(inner_deriv_re, row, col, block_in_re, found)
908 798 : CALL dbcsr_get_block_p(inner_deriv_im, row, col, block_in_im, found)
909 798 : CALL dbcsr_get_block_p(outer_deriv_re, row, col, block_out_re, found)
910 798 : CALL dbcsr_get_block_p(outer_deriv_im, row, col, block_out_im, found)
911 :
912 6011 : DO i = 1, SIZE(block_in_re, 1)
913 25610 : DO j = 1, SIZE(block_in_re, 2)
914 21195 : e1 = qs_ot_env%rot_mat_evals(row_blk_offset(row) + i - 1)
915 21195 : e2 = qs_ot_env%rot_mat_evals(col_blk_offset(col) + j - 1)
916 21195 : cval_in = CMPLX(block_in_re(i, j), block_in_im(i, j), dp)
917 21195 : cval_out = cval_in*cint(e1, e2)
918 21195 : block_out_re(i, j) = REAL(cval_out)
919 24812 : block_out_im(i, j) = AIMAG(cval_out)
920 : END DO
921 : END DO
922 : END DO
923 1596 : CALL dbcsr_iterator_stop(iter)
924 1596 : CALL dbcsr_release(inner_deriv_re)
925 1596 : CALL dbcsr_release(inner_deriv_im)
926 :
927 : ! Compute: matrix_buf1 = \sum_i outer_deriv_ij * |rot_mat_evec_i> <rot_mat_evec_j|
928 1596 : CALL dbcsr_multiply('N', 'N', +1.0_dp, qs_ot_env%rot_mat_evec_re, outer_deriv_re, 0.0_dp, mat_buf)
929 1596 : CALL dbcsr_multiply('N', 'N', -1.0_dp, qs_ot_env%rot_mat_evec_im, outer_deriv_im, 1.0_dp, mat_buf)
930 1596 : CALL dbcsr_multiply('N', 'T', +1.0_dp, mat_buf, qs_ot_env%rot_mat_evec_re, 0.0_dp, qs_ot_env%matrix_buf1)
931 :
932 1596 : CALL dbcsr_multiply('N', 'N', +1.0_dp, qs_ot_env%rot_mat_evec_re, outer_deriv_im, 0.0_dp, mat_buf)
933 1596 : CALL dbcsr_multiply('N', 'N', +1.0_dp, qs_ot_env%rot_mat_evec_im, outer_deriv_re, 1.0_dp, mat_buf)
934 1596 : CALL dbcsr_multiply('N', 'T', +1.0_dp, mat_buf, qs_ot_env%rot_mat_evec_im, 1.0_dp, qs_ot_env%matrix_buf1)
935 :
936 : ! Account for anti-symmetry of rot_mat_x.
937 1596 : CALL dbcsr_get_info(qs_ot_env%matrix_buf3, distribution=dist)
938 : CALL dbcsr_transposed(qs_ot_env%matrix_buf2, qs_ot_env%matrix_buf1, &
939 : shallow_data_copy=.FALSE., use_distribution=dist, &
940 1596 : transpose_distribution=.FALSE.)
941 :
942 : ! rot_mat_gx = matrix_buf1^T - matrix_buf1
943 1596 : CALL dbcsr_add(qs_ot_env%matrix_buf1, qs_ot_env%matrix_buf2, alpha_scalar=-1.0_dp, beta_scalar=+1.0_dp)
944 1596 : CALL dbcsr_copy(qs_ot_env%rot_mat_gx, qs_ot_env%matrix_buf1)
945 :
946 1596 : CALL dbcsr_release(mat_buf)
947 1596 : CALL dbcsr_release(outer_deriv_re)
948 1596 : CALL dbcsr_release(outer_deriv_im)
949 : END IF
950 3244 : CALL timestop(handle)
951 : CONTAINS
952 :
953 : ! **************************************************************************************************
954 : !> \brief ...
955 : !> \param e1 ...
956 : !> \param e2 ...
957 : !> \return ...
958 : ! **************************************************************************************************
959 21195 : FUNCTION cint(e1, e2)
960 : REAL(KIND=dp) :: e1, e2
961 : COMPLEX(KIND=dp) :: cint
962 :
963 : COMPLEX(KIND=dp) :: l1, l2, x
964 : INTEGER :: I
965 :
966 21195 : l1 = (0.0_dp, -1.0_dp)*e1
967 21195 : l2 = (0.0_dp, -1.0_dp)*e2
968 21195 : IF (ABS(l1 - l2) .GT. 0.5_dp) THEN
969 994 : cint = (EXP(l1) - EXP(l2))/(l1 - l2)
970 : ELSE
971 : x = 1.0_dp
972 : cint = 0.0_dp
973 343417 : DO I = 1, 16
974 323216 : cint = cint + x
975 343417 : x = x*(l1 - l2)/REAL(I + 1, KIND=dp)
976 : END DO
977 20201 : cint = cint*EXP(l2)
978 : END IF
979 21195 : END FUNCTION cint
980 : END SUBROUTINE qs_ot_rot_mat_derivative
981 :
982 : ! **************************************************************************************************
983 : !> \brief decide strategy
984 : !> tries to decide if the taylor expansion of cos(sqrt(xsx)) converges rapidly enough
985 : !> to make a taylor expansion of the functions cos(sqrt(xsx)) and sin(sqrt(xsx))/sqrt(xsx)
986 : !> and their derivatives faster than their computation based on diagonalization since xsx can
987 : !> be very small, especially during dynamics, only a few terms might indeed be needed we find
988 : !> the necessary order N to have largest_eval_upper_bound**(N+1)/(2(N+1))! < eps_taylor
989 : !> \param qs_ot_env ...
990 : ! **************************************************************************************************
991 96305 : SUBROUTINE decide_strategy(qs_ot_env)
992 : TYPE(qs_ot_type) :: qs_ot_env
993 :
994 : INTEGER :: N
995 : REAL(KIND=dp) :: num_error
996 :
997 96305 : qs_ot_env%do_taylor = .FALSE.
998 96305 : N = 0
999 96305 : num_error = qs_ot_env%largest_eval_upper_bound/(2.0_dp)
1000 410937 : DO WHILE (num_error > qs_ot_env%settings%eps_taylor .AND. N <= qs_ot_env%settings%max_taylor)
1001 314632 : N = N + 1
1002 353899 : num_error = num_error*qs_ot_env%largest_eval_upper_bound/REAL((2*N + 1)*(2*N + 2), KIND=dp)
1003 : END DO
1004 96305 : qs_ot_env%taylor_order = N
1005 96305 : IF (qs_ot_env%taylor_order <= qs_ot_env%settings%max_taylor) THEN
1006 52174 : qs_ot_env%do_taylor = .TRUE.
1007 : END IF
1008 :
1009 96305 : END SUBROUTINE decide_strategy
1010 :
1011 : ! **************************************************************************************************
1012 : !> \brief c=(c0*cos(p^0.5)+x*sin(p^0.5)*p^(-0.5)) x rot_mat_u
1013 : !> this assumes that x is already ortho to S*C0, and that p is x*S*x
1014 : !> rot_mat_u is an optional rotation matrix
1015 : !> \param matrix_c ...
1016 : !> \param matrix_x ...
1017 : !> \param qs_ot_env ...
1018 : ! **************************************************************************************************
1019 178918 : SUBROUTINE qs_ot_get_orbitals(matrix_c, matrix_x, qs_ot_env)
1020 :
1021 : TYPE(dbcsr_type), POINTER :: matrix_c, matrix_x
1022 : TYPE(qs_ot_type) :: qs_ot_env
1023 :
1024 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_ot_get_orbitals'
1025 : REAL(KIND=dp), PARAMETER :: rone = 1.0_dp, rzero = 0.0_dp
1026 :
1027 : INTEGER :: handle, k, n
1028 : TYPE(dbcsr_type), POINTER :: matrix_kk
1029 :
1030 89459 : CALL timeset(routineN, handle)
1031 :
1032 89459 : CALL dbcsr_get_info(matrix_x, nfullrows_total=n, nfullcols_total=k)
1033 :
1034 : ! rotate the multiplying matrices cosp and sinp instead of the result,
1035 : ! this should be cheaper for large basis sets
1036 89459 : IF (qs_ot_env%settings%do_rotation) THEN
1037 2996 : matrix_kk => qs_ot_env%matrix_buf1
1038 : CALL dbcsr_multiply('N', 'N', rone, qs_ot_env%matrix_cosp, &
1039 2996 : qs_ot_env%rot_mat_u, rzero, matrix_kk)
1040 : ELSE
1041 86463 : matrix_kk => qs_ot_env%matrix_cosp
1042 : END IF
1043 :
1044 : CALL dbcsr_multiply('N', 'N', rone, qs_ot_env%matrix_c0, matrix_kk, &
1045 89459 : rzero, matrix_c)
1046 :
1047 89459 : IF (qs_ot_env%settings%do_rotation) THEN
1048 2996 : matrix_kk => qs_ot_env%matrix_buf1
1049 : CALL dbcsr_multiply('N', 'N', rone, qs_ot_env%matrix_sinp, &
1050 2996 : qs_ot_env%rot_mat_u, rzero, matrix_kk)
1051 : ELSE
1052 86463 : matrix_kk => qs_ot_env%matrix_sinp
1053 : END IF
1054 : CALL dbcsr_multiply('N', 'N', rone, matrix_x, matrix_kk, &
1055 89459 : rone, matrix_c)
1056 :
1057 89459 : CALL timestop(handle)
1058 :
1059 89459 : END SUBROUTINE qs_ot_get_orbitals
1060 :
1061 : ! **************************************************************************************************
1062 : !> \brief this routines computes dE/dx=dx, with dx ortho to sc0
1063 : !> needs dE/dC=hc,C0,X,SX,p
1064 : !> if preconditioned it will not be the derivative, but the lagrangian multiplier
1065 : !> is changed so that P*dE/dx is the right derivative (i.e. in the allowed subspace)
1066 : !> \param matrix_hc ...
1067 : !> \param matrix_x ...
1068 : !> \param matrix_sx ...
1069 : !> \param matrix_gx ...
1070 : !> \param qs_ot_env ...
1071 : ! **************************************************************************************************
1072 200547 : SUBROUTINE qs_ot_get_derivative(matrix_hc, matrix_x, matrix_sx, matrix_gx, &
1073 : qs_ot_env)
1074 : TYPE(dbcsr_type), POINTER :: matrix_hc, matrix_x, matrix_sx, matrix_gx
1075 : TYPE(qs_ot_type) :: qs_ot_env
1076 :
1077 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_ot_get_derivative'
1078 : REAL(KIND=dp), PARAMETER :: rone = 1.0_dp, rzero = 0.0_dp
1079 :
1080 : INTEGER :: handle, k, n, ortho_k
1081 : TYPE(dbcsr_type), POINTER :: matrix_hc_local, matrix_target
1082 :
1083 66849 : CALL timeset(routineN, handle)
1084 :
1085 66849 : NULLIFY (matrix_hc_local)
1086 :
1087 66849 : CALL dbcsr_get_info(matrix_x, nfullrows_total=n, nfullcols_total=k)
1088 :
1089 : ! could in principle be taken inside qs_ot_get_derivative_* for increased efficiency
1090 : ! create a local rotated version of matrix_hc leaving matrix_hc untouched (needed
1091 : ! for lagrangian multipliers)
1092 66849 : IF (qs_ot_env%settings%do_rotation) THEN
1093 1622 : CALL dbcsr_copy(matrix_gx, matrix_hc) ! use gx as temporary
1094 1622 : CALL dbcsr_init_p(matrix_hc_local)
1095 1622 : CALL dbcsr_copy(matrix_hc_local, matrix_hc, name='matrix_hc_local')
1096 1622 : CALL dbcsr_set(matrix_hc_local, 0.0_dp)
1097 1622 : CALL dbcsr_multiply('N', 'T', rone, matrix_gx, qs_ot_env%rot_mat_u, rzero, matrix_hc_local)
1098 : ELSE
1099 65227 : matrix_hc_local => matrix_hc
1100 : END IF
1101 :
1102 66849 : IF (qs_ot_env%do_taylor) THEN
1103 37395 : CALL qs_ot_get_derivative_taylor(matrix_hc_local, matrix_x, matrix_sx, matrix_gx, qs_ot_env)
1104 : ELSE
1105 29454 : CALL qs_ot_get_derivative_diag(matrix_hc_local, matrix_x, matrix_sx, matrix_gx, qs_ot_env)
1106 : END IF
1107 :
1108 : ! and make it orthogonal
1109 66849 : CALL dbcsr_get_info(qs_ot_env%matrix_sc0, nfullcols_total=ortho_k)
1110 :
1111 66849 : IF (ASSOCIATED(qs_ot_env%preconditioner)) THEN
1112 56027 : matrix_target => qs_ot_env%matrix_psc0
1113 : ELSE
1114 10822 : matrix_target => qs_ot_env%matrix_sc0
1115 : END IF
1116 : ! first make the matrix os if not yet valid
1117 66849 : IF (.NOT. qs_ot_env%os_valid) THEN
1118 : ! this assumes that the preconditioner is a single matrix
1119 : ! that maps sc0 onto psc0
1120 :
1121 7380 : IF (ASSOCIATED(qs_ot_env%preconditioner)) THEN
1122 : CALL apply_preconditioner(qs_ot_env%preconditioner, qs_ot_env%matrix_sc0, &
1123 6410 : qs_ot_env%matrix_psc0)
1124 : END IF
1125 : CALL dbcsr_multiply('T', 'N', rone, &
1126 : qs_ot_env%matrix_sc0, matrix_target, &
1127 7380 : rzero, qs_ot_env%matrix_os)
1128 : CALL cp_dbcsr_cholesky_decompose(qs_ot_env%matrix_os, &
1129 7380 : para_env=qs_ot_env%para_env, blacs_env=qs_ot_env%blacs_env)
1130 : CALL cp_dbcsr_cholesky_invert(qs_ot_env%matrix_os, &
1131 : para_env=qs_ot_env%para_env, blacs_env=qs_ot_env%blacs_env, &
1132 7380 : uplo_to_full=.TRUE.)
1133 7380 : qs_ot_env%os_valid = .TRUE.
1134 : END IF
1135 : CALL dbcsr_multiply('T', 'N', rone, matrix_target, matrix_gx, &
1136 66849 : rzero, qs_ot_env%matrix_buf1_ortho)
1137 : CALL dbcsr_multiply('N', 'N', rone, qs_ot_env%matrix_os, &
1138 66849 : qs_ot_env%matrix_buf1_ortho, rzero, qs_ot_env%matrix_buf2_ortho)
1139 : CALL dbcsr_multiply('N', 'N', -rone, qs_ot_env%matrix_sc0, &
1140 66849 : qs_ot_env%matrix_buf2_ortho, rone, matrix_gx)
1141 : ! also treat the rot_mat gradient here
1142 66849 : IF (qs_ot_env%settings%do_rotation) THEN
1143 1622 : CALL qs_ot_rot_mat_derivative(qs_ot_env)
1144 : END IF
1145 :
1146 66849 : IF (qs_ot_env%settings%do_rotation) THEN
1147 1622 : CALL dbcsr_release_p(matrix_hc_local)
1148 : END IF
1149 :
1150 66849 : CALL timestop(handle)
1151 :
1152 66849 : END SUBROUTINE qs_ot_get_derivative
1153 :
1154 : ! **************************************************************************************************
1155 : !> \brief ...
1156 : !> \param matrix_hc ...
1157 : !> \param matrix_x ...
1158 : !> \param matrix_sx ...
1159 : !> \param matrix_gx ...
1160 : !> \param qs_ot_env ...
1161 : ! **************************************************************************************************
1162 88362 : SUBROUTINE qs_ot_get_derivative_diag(matrix_hc, matrix_x, matrix_sx, matrix_gx, &
1163 : qs_ot_env)
1164 :
1165 : TYPE(dbcsr_type), POINTER :: matrix_hc, matrix_x, matrix_sx, matrix_gx
1166 : TYPE(qs_ot_type) :: qs_ot_env
1167 :
1168 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_ot_get_derivative_diag'
1169 : REAL(KIND=dp), PARAMETER :: rone = 1.0_dp, rzero = 0.0_dp
1170 :
1171 : INTEGER :: handle, k, n
1172 : TYPE(dbcsr_distribution_type) :: dist
1173 :
1174 29454 : CALL timeset(routineN, handle)
1175 :
1176 29454 : CALL dbcsr_get_info(matrix_x, nfullrows_total=n, nfullcols_total=k)
1177 :
1178 : ! go for the derivative now
1179 : ! this de/dc*(dX/dx)*sinp
1180 29454 : CALL dbcsr_multiply('N', 'N', rone, matrix_hc, qs_ot_env%matrix_sinp, rzero, matrix_gx)
1181 : ! overlap hc*x
1182 29454 : CALL dbcsr_multiply('T', 'N', rone, matrix_hc, matrix_x, rzero, qs_ot_env%matrix_buf2)
1183 : ! get it in the basis of the eigenvectors
1184 : CALL dbcsr_multiply('N', 'N', rone, qs_ot_env%matrix_buf2, qs_ot_env%matrix_r, &
1185 29454 : rzero, qs_ot_env%matrix_buf1)
1186 : CALL dbcsr_multiply('T', 'N', rone, qs_ot_env%matrix_r, qs_ot_env%matrix_buf1, &
1187 29454 : rzero, qs_ot_env%matrix_buf2)
1188 :
1189 : ! get the schur product of O_uv*B_uv
1190 : CALL dbcsr_hadamard_product(qs_ot_env%matrix_buf2, qs_ot_env%matrix_sinp_b, &
1191 29454 : qs_ot_env%matrix_buf3)
1192 :
1193 : ! overlap hc*c0
1194 : CALL dbcsr_multiply('T', 'N', rone, matrix_hc, qs_ot_env%matrix_c0, rzero, &
1195 29454 : qs_ot_env%matrix_buf2)
1196 : ! get it in the basis of the eigenvectors
1197 : CALL dbcsr_multiply('N', 'N', rone, qs_ot_env%matrix_buf2, qs_ot_env%matrix_r, &
1198 29454 : rzero, qs_ot_env%matrix_buf1)
1199 : CALL dbcsr_multiply('T', 'N', rone, qs_ot_env%matrix_r, qs_ot_env%matrix_buf1, &
1200 29454 : rzero, qs_ot_env%matrix_buf2)
1201 :
1202 : CALL dbcsr_hadamard_product(qs_ot_env%matrix_buf2, qs_ot_env%matrix_cosp_b, &
1203 29454 : qs_ot_env%matrix_buf4)
1204 :
1205 : ! add the two bs and compute b+b^T
1206 : CALL dbcsr_add(qs_ot_env%matrix_buf3, qs_ot_env%matrix_buf4, &
1207 29454 : alpha_scalar=rone, beta_scalar=rone)
1208 :
1209 : ! get the b in the eigenvector basis
1210 : CALL dbcsr_multiply('N', 'T', rone, qs_ot_env%matrix_buf3, qs_ot_env%matrix_r, &
1211 29454 : rzero, qs_ot_env%matrix_buf1)
1212 : CALL dbcsr_multiply('N', 'N', rone, qs_ot_env%matrix_r, qs_ot_env%matrix_buf1, &
1213 29454 : rzero, qs_ot_env%matrix_buf3)
1214 29454 : CALL dbcsr_get_info(qs_ot_env%matrix_buf3, distribution=dist)
1215 : CALL dbcsr_transposed(qs_ot_env%matrix_buf1, qs_ot_env%matrix_buf3, &
1216 : shallow_data_copy=.FALSE., use_distribution=dist, &
1217 29454 : transpose_distribution=.FALSE.)
1218 : CALL dbcsr_add(qs_ot_env%matrix_buf3, qs_ot_env%matrix_buf1, &
1219 29454 : alpha_scalar=rone, beta_scalar=rone)
1220 :
1221 : ! and add to the derivative
1222 : CALL dbcsr_multiply('N', 'N', rone, matrix_sx, qs_ot_env%matrix_buf3, &
1223 29454 : rone, matrix_gx)
1224 29454 : CALL timestop(handle)
1225 :
1226 29454 : END SUBROUTINE qs_ot_get_derivative_diag
1227 :
1228 : ! **************************************************************************************************
1229 : !> \brief compute the derivative of the taylor expansion below
1230 : !> \param matrix_hc ...
1231 : !> \param matrix_x ...
1232 : !> \param matrix_sx ...
1233 : !> \param matrix_gx ...
1234 : !> \param qs_ot_env ...
1235 : ! **************************************************************************************************
1236 134312 : SUBROUTINE qs_ot_get_derivative_taylor(matrix_hc, matrix_x, matrix_sx, matrix_gx, &
1237 : qs_ot_env)
1238 :
1239 : TYPE(dbcsr_type), POINTER :: matrix_hc, matrix_x, matrix_sx, matrix_gx
1240 : TYPE(qs_ot_type) :: qs_ot_env
1241 :
1242 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_ot_get_derivative_taylor'
1243 : REAL(KIND=dp), PARAMETER :: rone = 1.0_dp, rzero = 0.0_dp
1244 :
1245 : INTEGER :: handle, i, k, n
1246 : REAL(KIND=dp) :: cosfactor, sinfactor
1247 : TYPE(dbcsr_distribution_type) :: dist
1248 : TYPE(dbcsr_type), POINTER :: matrix_left, matrix_right
1249 :
1250 37395 : CALL timeset(routineN, handle)
1251 :
1252 37395 : CALL dbcsr_get_info(matrix_x, nfullrows_total=n, nfullcols_total=k)
1253 :
1254 : ! go for the derivative now
1255 : ! this de/dc*(dX/dx)*sinp i.e. zeroth order
1256 37395 : CALL dbcsr_multiply('N', 'N', rone, matrix_hc, qs_ot_env%matrix_sinp, rzero, matrix_gx)
1257 :
1258 37395 : IF (qs_ot_env%taylor_order .LE. 0) THEN
1259 7634 : CALL timestop(handle)
1260 7634 : RETURN
1261 : END IF
1262 :
1263 : ! we store the matrix that will multiply sx in matrix_r
1264 29761 : CALL dbcsr_set(qs_ot_env%matrix_r, rzero)
1265 :
1266 : ! just better names for matrix_cosp_b and matrix_sinp_b (they are buffer space here)
1267 29761 : matrix_left => qs_ot_env%matrix_cosp_b
1268 29761 : matrix_right => qs_ot_env%matrix_sinp_b
1269 :
1270 : ! overlap hc*x and add its transpose to matrix_left
1271 29761 : CALL dbcsr_multiply('T', 'N', rone, matrix_hc, matrix_x, rzero, matrix_left)
1272 29761 : CALL dbcsr_get_info(matrix_left, distribution=dist)
1273 : CALL dbcsr_transposed(qs_ot_env%matrix_buf1, matrix_left, &
1274 : shallow_data_copy=.FALSE., use_distribution=dist, &
1275 29761 : transpose_distribution=.FALSE.)
1276 : CALL dbcsr_add(matrix_left, qs_ot_env%matrix_buf1, &
1277 29761 : alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
1278 29761 : CALL dbcsr_copy(matrix_right, matrix_left)
1279 :
1280 : ! first order
1281 29761 : sinfactor = -1.0_dp/(2.0_dp*3.0_dp)
1282 : CALL dbcsr_add(qs_ot_env%matrix_r, matrix_left, &
1283 29761 : alpha_scalar=1.0_dp, beta_scalar=sinfactor)
1284 :
1285 : ! M
1286 : ! OM+MO
1287 : ! OOM+OMO+MOO
1288 : ! ...
1289 63341 : DO i = 2, qs_ot_env%taylor_order
1290 33580 : sinfactor = sinfactor*(-1.0_dp)/REAL(2*i*(2*i + 1), KIND=dp)
1291 33580 : CALL dbcsr_multiply('N', 'N', rone, qs_ot_env%matrix_p, matrix_left, rzero, qs_ot_env%matrix_buf1)
1292 33580 : CALL dbcsr_multiply('N', 'N', rone, matrix_right, qs_ot_env%matrix_p, rzero, matrix_left)
1293 33580 : CALL dbcsr_copy(matrix_right, matrix_left)
1294 : CALL dbcsr_add(matrix_left, qs_ot_env%matrix_buf1, &
1295 33580 : 1.0_dp, 1.0_dp)
1296 : CALL dbcsr_add(qs_ot_env%matrix_r, matrix_left, &
1297 63341 : alpha_scalar=1.0_dp, beta_scalar=sinfactor)
1298 : END DO
1299 :
1300 : ! overlap hc*c0 and add its transpose to matrix_left
1301 29761 : CALL dbcsr_multiply('T', 'N', rone, matrix_hc, qs_ot_env%matrix_c0, rzero, matrix_left)
1302 29761 : CALL dbcsr_get_info(matrix_left, distribution=dist)
1303 : CALL dbcsr_transposed(qs_ot_env%matrix_buf1, matrix_left, &
1304 : shallow_data_copy=.FALSE., use_distribution=dist, &
1305 29761 : transpose_distribution=.FALSE.)
1306 29761 : CALL dbcsr_add(matrix_left, qs_ot_env%matrix_buf1, 1.0_dp, 1.0_dp)
1307 29761 : CALL dbcsr_copy(matrix_right, matrix_left)
1308 :
1309 : ! first order
1310 29761 : cosfactor = -1.0_dp/(1.0_dp*2.0_dp)
1311 : CALL dbcsr_add(qs_ot_env%matrix_r, matrix_left, &
1312 29761 : alpha_scalar=1.0_dp, beta_scalar=cosfactor)
1313 :
1314 : ! M
1315 : ! OM+MO
1316 : ! OOM+OMO+MOO
1317 : ! ...
1318 63341 : DO i = 2, qs_ot_env%taylor_order
1319 33580 : cosfactor = cosfactor*(-1.0_dp)/REAL(2*i*(2*i - 1), KIND=dp)
1320 33580 : CALL dbcsr_multiply('N', 'N', rone, qs_ot_env%matrix_p, matrix_left, rzero, qs_ot_env%matrix_buf1)
1321 33580 : CALL dbcsr_multiply('N', 'N', rone, matrix_right, qs_ot_env%matrix_p, rzero, matrix_left)
1322 33580 : CALL dbcsr_copy(matrix_right, matrix_left)
1323 33580 : CALL dbcsr_add(matrix_left, qs_ot_env%matrix_buf1, 1.0_dp, 1.0_dp)
1324 : CALL dbcsr_add(qs_ot_env%matrix_r, matrix_left, &
1325 63341 : alpha_scalar=1.0_dp, beta_scalar=cosfactor)
1326 : END DO
1327 :
1328 : ! and add to the derivative
1329 29761 : CALL dbcsr_multiply('N', 'N', rone, matrix_sx, qs_ot_env%matrix_r, rone, matrix_gx)
1330 :
1331 29761 : CALL timestop(handle)
1332 :
1333 37395 : END SUBROUTINE qs_ot_get_derivative_taylor
1334 :
1335 : ! *************************************************************************************************
1336 : !> \brief computes a taylor expansion.
1337 : !> \param qs_ot_env ...
1338 : ! **************************************************************************************************
1339 85123 : SUBROUTINE qs_ot_p2m_taylor(qs_ot_env)
1340 : TYPE(qs_ot_type) :: qs_ot_env
1341 :
1342 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_ot_p2m_taylor'
1343 : REAL(KIND=dp), PARAMETER :: rone = 1.0_dp, rzero = 0.0_dp
1344 :
1345 : INTEGER :: handle, i, k
1346 : REAL(KIND=dp) :: cosfactor, sinfactor
1347 :
1348 52174 : CALL timeset(routineN, handle)
1349 :
1350 : ! zeroth order
1351 52174 : CALL dbcsr_set(qs_ot_env%matrix_cosp, rzero)
1352 52174 : CALL dbcsr_set(qs_ot_env%matrix_sinp, rzero)
1353 52174 : CALL dbcsr_add_on_diag(qs_ot_env%matrix_cosp, rone)
1354 52174 : CALL dbcsr_add_on_diag(qs_ot_env%matrix_sinp, rone)
1355 :
1356 52174 : IF (qs_ot_env%taylor_order .LE. 0) THEN
1357 8220 : CALL timestop(handle)
1358 19225 : RETURN
1359 : END IF
1360 :
1361 : ! first order
1362 43954 : cosfactor = -1.0_dp/(1.0_dp*2.0_dp)
1363 43954 : sinfactor = -1.0_dp/(2.0_dp*3.0_dp)
1364 43954 : CALL dbcsr_add(qs_ot_env%matrix_cosp, qs_ot_env%matrix_p, alpha_scalar=1.0_dp, beta_scalar=cosfactor)
1365 43954 : CALL dbcsr_add(qs_ot_env%matrix_sinp, qs_ot_env%matrix_p, alpha_scalar=1.0_dp, beta_scalar=sinfactor)
1366 43954 : IF (qs_ot_env%taylor_order .LE. 1) THEN
1367 11005 : CALL timestop(handle)
1368 11005 : RETURN
1369 : END IF
1370 :
1371 : ! other orders
1372 32949 : CALL dbcsr_get_info(qs_ot_env%matrix_p, nfullrows_total=k)
1373 32949 : CALL dbcsr_copy(qs_ot_env%matrix_r, qs_ot_env%matrix_p)
1374 :
1375 82972 : DO i = 2, qs_ot_env%taylor_order
1376 : ! new power of p
1377 : CALL dbcsr_multiply('N', 'N', rone, qs_ot_env%matrix_p, qs_ot_env%matrix_r, &
1378 50023 : rzero, qs_ot_env%matrix_buf1)
1379 50023 : CALL dbcsr_copy(qs_ot_env%matrix_r, qs_ot_env%matrix_buf1)
1380 : ! add to the taylor expansion so far
1381 50023 : cosfactor = cosfactor*(-1.0_dp)/REAL(2*i*(2*i - 1), KIND=dp)
1382 50023 : sinfactor = sinfactor*(-1.0_dp)/REAL(2*i*(2*i + 1), KIND=dp)
1383 : CALL dbcsr_add(qs_ot_env%matrix_cosp, qs_ot_env%matrix_r, &
1384 50023 : alpha_scalar=1.0_dp, beta_scalar=cosfactor)
1385 : CALL dbcsr_add(qs_ot_env%matrix_sinp, qs_ot_env%matrix_r, &
1386 82972 : alpha_scalar=1.0_dp, beta_scalar=sinfactor)
1387 : END DO
1388 :
1389 32949 : CALL timestop(handle)
1390 :
1391 : END SUBROUTINE qs_ot_p2m_taylor
1392 :
1393 : ! **************************************************************************************************
1394 : !> \brief given p, computes - eigenstuff (matrix_r,evals)
1395 : !> - cos(p^0.5),p^(-0.5)*sin(p^0.5)
1396 : !> - the real b matrices, needed for the derivatives of these guys
1397 : !> cosp_b_ij=(1/(2pii) * int(cos(z^1/2)/((z-eval(i))*(z-eval(j))))
1398 : !> sinp_b_ij=(1/(2pii) * int(z^(-1/2)*sin(z^1/2)/((z-eval(i))*(z-eval(j))))
1399 : !> \param qs_ot_env ...
1400 : ! **************************************************************************************************
1401 176524 : SUBROUTINE qs_ot_p2m_diag(qs_ot_env)
1402 :
1403 : TYPE(qs_ot_type) :: qs_ot_env
1404 :
1405 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_ot_p2m_diag'
1406 : REAL(KIND=dp), PARAMETER :: rone = 1.0_dp, rzero = 0.0_dp
1407 :
1408 : INTEGER :: col, col_offset, col_size, handle, i, j, &
1409 : k, row, row_offset, row_size
1410 44131 : REAL(dp), DIMENSION(:, :), POINTER :: block
1411 : REAL(KIND=dp) :: a, b
1412 : TYPE(dbcsr_iterator_type) :: iter
1413 :
1414 44131 : CALL timeset(routineN, handle)
1415 :
1416 44131 : CALL dbcsr_get_info(qs_ot_env%matrix_p, nfullrows_total=k)
1417 44131 : CALL dbcsr_copy(qs_ot_env%matrix_buf1, qs_ot_env%matrix_p)
1418 : CALL cp_dbcsr_syevd(qs_ot_env%matrix_buf1, qs_ot_env%matrix_r, qs_ot_env%evals, &
1419 44131 : qs_ot_env%para_env, qs_ot_env%blacs_env)
1420 504496 : DO i = 1, k
1421 504496 : qs_ot_env%evals(i) = MAX(0.0_dp, qs_ot_env%evals(i))
1422 : END DO
1423 :
1424 44131 : !$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(i) SHARED(k,qs_ot_env)
1425 : DO i = 1, k
1426 : qs_ot_env%dum(i) = COS(SQRT(qs_ot_env%evals(i)))
1427 : END DO
1428 44131 : CALL dbcsr_copy(qs_ot_env%matrix_buf1, qs_ot_env%matrix_r)
1429 44131 : CALL dbcsr_scale_by_vector(qs_ot_env%matrix_buf1, alpha=qs_ot_env%dum, side='right')
1430 : CALL dbcsr_multiply('N', 'T', rone, qs_ot_env%matrix_r, qs_ot_env%matrix_buf1, &
1431 44131 : rzero, qs_ot_env%matrix_cosp)
1432 :
1433 44131 : !$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(i) SHARED(k,qs_ot_env)
1434 : DO i = 1, k
1435 : qs_ot_env%dum(i) = qs_ot_sinc(SQRT(qs_ot_env%evals(i)))
1436 : END DO
1437 44131 : CALL dbcsr_copy(qs_ot_env%matrix_buf1, qs_ot_env%matrix_r)
1438 44131 : CALL dbcsr_scale_by_vector(qs_ot_env%matrix_buf1, alpha=qs_ot_env%dum, side='right')
1439 : CALL dbcsr_multiply('N', 'T', rone, qs_ot_env%matrix_r, qs_ot_env%matrix_buf1, &
1440 44131 : rzero, qs_ot_env%matrix_sinp)
1441 :
1442 44131 : CALL dbcsr_copy(qs_ot_env%matrix_cosp_b, qs_ot_env%matrix_cosp)
1443 44131 : CALL dbcsr_iterator_start(iter, qs_ot_env%matrix_cosp_b)
1444 77688 : DO WHILE (dbcsr_iterator_blocks_left(iter))
1445 : CALL dbcsr_iterator_next_block(iter, row, col, block, &
1446 : row_size=row_size, col_size=col_size, &
1447 33557 : row_offset=row_offset, col_offset=col_offset)
1448 543815 : DO j = 1, col_size
1449 11064784 : DO i = 1, row_size
1450 : a = (SQRT(qs_ot_env%evals(row_offset + i - 1)) &
1451 10565100 : - SQRT(qs_ot_env%evals(col_offset + j - 1)))/2.0_dp
1452 : b = (SQRT(qs_ot_env%evals(row_offset + i - 1)) &
1453 10565100 : + SQRT(qs_ot_env%evals(col_offset + j - 1)))/2.0_dp
1454 11031227 : block(i, j) = -0.5_dp*qs_ot_sinc(a)*qs_ot_sinc(b)
1455 : END DO
1456 : END DO
1457 : END DO
1458 44131 : CALL dbcsr_iterator_stop(iter)
1459 :
1460 44131 : CALL dbcsr_copy(qs_ot_env%matrix_sinp_b, qs_ot_env%matrix_sinp)
1461 44131 : CALL dbcsr_iterator_start(iter, qs_ot_env%matrix_sinp_b)
1462 77688 : DO WHILE (dbcsr_iterator_blocks_left(iter))
1463 : CALL dbcsr_iterator_next_block(iter, row, col, block, &
1464 : row_size=row_size, col_size=col_size, &
1465 33557 : row_offset=row_offset, col_offset=col_offset)
1466 543815 : DO j = 1, col_size
1467 11064784 : DO i = 1, row_size
1468 10565100 : a = SQRT(qs_ot_env%evals(row_offset + i - 1))
1469 10565100 : b = SQRT(qs_ot_env%evals(col_offset + j - 1))
1470 11031227 : block(i, j) = qs_ot_sincf(a, b)
1471 : END DO
1472 : END DO
1473 : END DO
1474 44131 : CALL dbcsr_iterator_stop(iter)
1475 :
1476 44131 : CALL timestop(handle)
1477 :
1478 44131 : END SUBROUTINE qs_ot_p2m_diag
1479 :
1480 : ! **************************************************************************************************
1481 : !> \brief computes sin(x)/x for all values of the argument
1482 : !> \param x ...
1483 : !> \return ...
1484 : ! **************************************************************************************************
1485 29282511 : FUNCTION qs_ot_sinc(x)
1486 :
1487 : REAL(KIND=dp), INTENT(IN) :: x
1488 : REAL(KIND=dp) :: qs_ot_sinc
1489 :
1490 : REAL(KIND=dp), PARAMETER :: q1 = 1.0_dp, q2 = -q1/(2.0_dp*3.0_dp), q3 = -q2/(4.0_dp*5.0_dp), &
1491 : q4 = -q3/(6.0_dp*7.0_dp), q5 = -q4/(8.0_dp*9.0_dp), q6 = -q5/(10.0_dp*11.0_dp), &
1492 : q7 = -q6/(12.0_dp*13.0_dp), q8 = -q7/(14.0_dp*15.0_dp), q9 = -q8/(16.0_dp*17.0_dp), &
1493 : q10 = -q9/(18.0_dp*19.0_dp)
1494 :
1495 : REAL(KIND=dp) :: y
1496 :
1497 29282511 : IF (ABS(x) > 0.5_dp) THEN
1498 8635159 : qs_ot_sinc = SIN(x)/x
1499 : ELSE
1500 20647352 : y = x*x
1501 20647352 : qs_ot_sinc = q1 + y*(q2 + y*(q3 + y*(q4 + y*(q5 + y*(q6 + y*(q7 + y*(q8 + y*(q9 + y*(q10)))))))))
1502 : END IF
1503 29282511 : END FUNCTION qs_ot_sinc
1504 :
1505 : ! **************************************************************************************************
1506 : !> \brief computes (1/(x^2-y^2))*(sinc(x)-sinc(y)) for all positive values of the arguments
1507 : !> \param xa ...
1508 : !> \param ya ...
1509 : !> \return ...
1510 : ! **************************************************************************************************
1511 10565100 : FUNCTION qs_ot_sincf(xa, ya)
1512 :
1513 : REAL(KIND=dp), INTENT(IN) :: xa, ya
1514 : REAL(KIND=dp) :: qs_ot_sincf
1515 :
1516 : INTEGER :: i
1517 : REAL(KIND=dp) :: a, b, rs, sf, x, xs, y, ybx, ybxs
1518 :
1519 : ! this is currently a limit of the routine, could be removed rather easily
1520 10565100 : IF (xa .LT. 0) CPABORT("x is negative")
1521 10565100 : IF (ya .LT. 0) CPABORT("y is negative")
1522 :
1523 10565100 : IF (xa .LT. ya) THEN
1524 5074289 : x = ya
1525 5074289 : y = xa
1526 : ELSE
1527 5490811 : x = xa
1528 5490811 : y = ya
1529 : END IF
1530 :
1531 10565100 : IF (x .LT. 0.5_dp) THEN ! use series, keeping in mind that x,y,x+y,x-y can all be zero
1532 :
1533 6719127 : qs_ot_sincf = 0.0_dp
1534 6719127 : IF (x .GT. 0.0_dp) THEN
1535 6516829 : ybx = y/x
1536 : ELSE ! should be irrelevant !?
1537 : ybx = 0.0_dp
1538 : END IF
1539 :
1540 6719127 : sf = -1.0_dp/((1.0_dp + ybx)*6.0_dp)
1541 6719127 : rs = 1.0_dp
1542 6719127 : ybxs = ybx
1543 6719127 : xs = 1.0_dp
1544 :
1545 73910397 : DO i = 1, 10
1546 67191270 : qs_ot_sincf = qs_ot_sincf + sf*rs*xs*(1.0_dp + ybxs)
1547 67191270 : sf = -sf/(REAL((2*i + 2), dp)*REAL((2*i + 3), dp))
1548 67191270 : rs = rs + ybxs
1549 67191270 : ybxs = ybxs*ybx
1550 73910397 : xs = xs*x*x
1551 : END DO
1552 :
1553 : ELSE ! no series expansion
1554 3845973 : IF (x - y .GT. 0.1_dp) THEN ! safe to use the normal form
1555 3563604 : qs_ot_sincf = (qs_ot_sinc(x) - qs_ot_sinc(y))/((x + y)*(x - y))
1556 : ELSE
1557 282369 : a = (x + y)/2.0_dp
1558 282369 : b = (x - y)/2.0_dp ! might be close to zero
1559 : ! y (=(a-b)) can not be close to zero since it is close to x>0.5
1560 282369 : qs_ot_sincf = (qs_ot_sinc(b)*COS(a) - qs_ot_sinc(a)*COS(b))/(2*x*y)
1561 : END IF
1562 : END IF
1563 :
1564 10565100 : END FUNCTION qs_ot_sincf
1565 :
1566 : END MODULE qs_ot
|
