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