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