(git:ed6f26b)
Loading...
Searching...
No Matches
ct_methods.F
Go to the documentation of this file.
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 Cayley transformation methods
10!> \par History
11!> 2011.06 created [Rustam Z Khaliullin]
12!> \author Rustam Z Khaliullin
13! **************************************************************************************************
14MODULE ct_methods
15 USE cp_dbcsr_api, ONLY: &
16 dbcsr_add, dbcsr_copy, dbcsr_create, dbcsr_desymmetrize, dbcsr_filter, dbcsr_finalize, &
17 dbcsr_get_info, dbcsr_iterator_blocks_left, dbcsr_iterator_next_block, &
18 dbcsr_iterator_readonly_start, dbcsr_iterator_start, dbcsr_iterator_stop, &
19 dbcsr_iterator_type, dbcsr_multiply, dbcsr_put_block, dbcsr_release, dbcsr_scale, &
20 dbcsr_set, dbcsr_transposed, dbcsr_type, dbcsr_type_no_symmetry, dbcsr_work_create
21 USE cp_dbcsr_cholesky, ONLY: cp_dbcsr_cholesky_decompose,&
22 cp_dbcsr_cholesky_invert
23 USE cp_dbcsr_contrib, ONLY: dbcsr_add_on_diag,&
24 dbcsr_dot,&
25 dbcsr_frobenius_norm,&
26 dbcsr_get_diag,&
27 dbcsr_hadamard_product,&
28 dbcsr_maxabs,&
29 dbcsr_reserve_diag_blocks,&
30 dbcsr_set_diag
31 USE cp_dbcsr_diag, ONLY: cp_dbcsr_syevd
32 USE cp_log_handling, ONLY: cp_get_default_logger,&
33 cp_logger_get_default_unit_nr,&
34 cp_logger_type
35 USE ct_types, ONLY: ct_step_env_type
36 USE input_constants, ONLY: &
37 cg_dai_yuan, cg_fletcher, cg_fletcher_reeves, cg_hager_zhang, cg_hestenes_stiefel, &
38 cg_liu_storey, cg_polak_ribiere, cg_zero, tensor_orthogonal, tensor_up_down
39 USE iterate_matrix, ONLY: matrix_sqrt_newton_schulz
40 USE kinds, ONLY: dp
41 USE machine, ONLY: m_walltime
42#include "./base/base_uses.f90"
43
44 IMPLICIT NONE
45
46 PRIVATE
47
48 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'ct_methods'
49
50 ! Public subroutines
51 PUBLIC :: ct_step_execute, analytic_line_search, diagonalize_diagonal_blocks
52
53CONTAINS
54
55! **************************************************************************************************
56!> \brief Performs Cayley transformation
57!> \param cts_env ...
58!> \par History
59!> 2011.06 created [Rustam Z Khaliullin]
60!> \author Rustam Z Khaliullin
61! **************************************************************************************************
62 SUBROUTINE ct_step_execute(cts_env)
63
64 TYPE(ct_step_env_type) :: cts_env
65
66 CHARACTER(len=*), PARAMETER :: routinen = 'ct_step_execute'
67
68 INTEGER :: handle, n, preconditioner_type, unit_nr
69 REAL(kind=dp) :: gap_estimate, safety_margin
70 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: evals
71 TYPE(cp_logger_type), POINTER :: logger
72 TYPE(dbcsr_type) :: matrix_pp, matrix_pq, matrix_qp, &
73 matrix_qp_save, matrix_qq, oo1, &
74 oo1_sqrt, oo1_sqrt_inv, t_corr, tmp1, &
75 u_pp, u_qq
76
77!TYPE(dbcsr_type) :: rst_x1, rst_x2
78!REAL(KIND=dp) :: ener_tmp
79!TYPE(dbcsr_iterator_type) :: iter
80!INTEGER :: iblock_row,iblock_col,&
81! iblock_row_size,iblock_col_size
82!REAL(KIND=dp), DIMENSION(:,:), POINTER :: data_p
83
84 CALL timeset(routinen, handle)
85
86 logger => cp_get_default_logger()
87 IF (logger%para_env%is_source()) THEN
88 unit_nr = cp_logger_get_default_unit_nr(logger, local=.true.)
89 ELSE
90 unit_nr = -1
91 END IF
92
93 ! check if all input is in place and flags are consistent
94 IF (cts_env%update_q .AND. (.NOT. cts_env%update_p)) THEN
95 cpabort("q-update is possible only with p-update")
96 END IF
97
98 IF (cts_env%tensor_type .EQ. tensor_up_down) THEN
99 cpabort("riccati is not implemented for biorthogonal basis")
100 END IF
101
102 IF (.NOT. ASSOCIATED(cts_env%matrix_ks)) THEN
103 cpabort("KS matrix is not associated")
104 END IF
105
106 IF (cts_env%use_virt_orbs .AND. (.NOT. cts_env%use_occ_orbs)) THEN
107 cpabort("virtual orbs can be used only with occupied orbs")
108 END IF
109
110 IF (cts_env%use_occ_orbs) THEN
111 IF (.NOT. ASSOCIATED(cts_env%matrix_t)) THEN
112 cpabort("T matrix is not associated")
113 END IF
114 IF (.NOT. ASSOCIATED(cts_env%matrix_qp_template)) THEN
115 cpabort("QP template is not associated")
116 END IF
117 IF (.NOT. ASSOCIATED(cts_env%matrix_pq_template)) THEN
118 cpabort("PQ template is not associated")
119 END IF
120 END IF
121
122 IF (cts_env%use_virt_orbs) THEN
123 IF (.NOT. ASSOCIATED(cts_env%matrix_v)) THEN
124 cpabort("V matrix is not associated")
125 END IF
126 ELSE
127 IF (.NOT. ASSOCIATED(cts_env%matrix_p)) THEN
128 cpabort("P matrix is not associated")
129 END IF
130 END IF
131
132 IF (cts_env%tensor_type .NE. tensor_up_down .AND. &
133 cts_env%tensor_type .NE. tensor_orthogonal) THEN
134 cpabort("illegal tensor flag")
135 END IF
136
137 ! start real calculations
138 IF (cts_env%use_occ_orbs) THEN
139
140 ! create matrices for various ks blocks
141 CALL dbcsr_create(matrix_pp, &
142 template=cts_env%p_index_up, &
143 matrix_type=dbcsr_type_no_symmetry)
144 CALL dbcsr_create(matrix_qp, &
145 template=cts_env%matrix_qp_template, &
146 matrix_type=dbcsr_type_no_symmetry)
147 CALL dbcsr_create(matrix_qq, &
148 template=cts_env%q_index_up, &
149 matrix_type=dbcsr_type_no_symmetry)
150 CALL dbcsr_create(matrix_pq, &
151 template=cts_env%matrix_pq_template, &
152 matrix_type=dbcsr_type_no_symmetry)
153
154 ! create the residue matrix
155 CALL dbcsr_create(cts_env%matrix_res, &
156 template=cts_env%matrix_qp_template)
157
158 CALL assemble_ks_qp_blocks(cts_env%matrix_ks, &
159 cts_env%matrix_p, &
160 cts_env%matrix_t, &
161 cts_env%matrix_v, &
162 cts_env%q_index_down, &
163 cts_env%p_index_up, &
164 cts_env%q_index_up, &
165 matrix_pp, &
166 matrix_qq, &
167 matrix_qp, &
168 matrix_pq, &
169 cts_env%tensor_type, &
170 cts_env%use_virt_orbs, &
171 cts_env%eps_filter)
172
173 ! create a matrix of single-excitation amplitudes
174 CALL dbcsr_create(cts_env%matrix_x, &
175 template=cts_env%matrix_qp_template)
176 IF (ASSOCIATED(cts_env%matrix_x_guess)) THEN
177 CALL dbcsr_copy(cts_env%matrix_x, &
178 cts_env%matrix_x_guess)
179 IF (cts_env%tensor_type .EQ. tensor_orthogonal) THEN
180 ! bring x from contravariant-covariant representation
181 ! to the orthogonal/cholesky representation
182 ! use res as temporary storage
183 CALL dbcsr_multiply("N", "N", 1.0_dp, cts_env%q_index_down, &
184 cts_env%matrix_x, 0.0_dp, cts_env%matrix_res, &
185 filter_eps=cts_env%eps_filter)
186 CALL dbcsr_multiply("N", "N", 1.0_dp, cts_env%matrix_res, &
187 cts_env%p_index_up, 0.0_dp, &
188 cts_env%matrix_x, &
189 filter_eps=cts_env%eps_filter)
190 END IF
191 ELSE
192 ! set amplitudes to zero
193 CALL dbcsr_set(cts_env%matrix_x, 0.0_dp)
194 END IF
195
196 !SELECT CASE (cts_env%preconditioner_type)
197 !CASE (prec_eigenvector_blocks,prec_eigenvector_full)
198 preconditioner_type = 1
199 safety_margin = 2.0_dp
200 gap_estimate = 0.0001_dp
201 SELECT CASE (preconditioner_type)
202 CASE (1, 2)
203!RZK-warning diagonalization works only with orthogonal tensor!!!
204 ! find a better basis by diagonalizing diagonal blocks
205 ! first pp
206 CALL dbcsr_create(u_pp, template=matrix_pp, &
207 matrix_type=dbcsr_type_no_symmetry)
208 !IF (cts_env%preconditioner_type.eq.prec_eigenvector_full) THEN
209 IF (.true.) THEN
210 CALL dbcsr_get_info(matrix_pp, nfullrows_total=n)
211 ALLOCATE (evals(n))
212 CALL cp_dbcsr_syevd(matrix_pp, u_pp, evals, &
213 cts_env%para_env, cts_env%blacs_env)
214 DEALLOCATE (evals)
215 ELSE
216 CALL diagonalize_diagonal_blocks(matrix_pp, u_pp)
217 END IF
218 ! and now qq
219 CALL dbcsr_create(u_qq, template=matrix_qq, &
220 matrix_type=dbcsr_type_no_symmetry)
221 !IF (cts_env%preconditioner_type.eq.prec_eigenvector_full) THEN
222 IF (.true.) THEN
223 CALL dbcsr_get_info(matrix_qq, nfullrows_total=n)
224 ALLOCATE (evals(n))
225 CALL cp_dbcsr_syevd(matrix_qq, u_qq, evals, &
226 cts_env%para_env, cts_env%blacs_env)
227 DEALLOCATE (evals)
228 ELSE
229 CALL diagonalize_diagonal_blocks(matrix_qq, u_qq)
230 END IF
231
232 ! apply the transformation to all matrices
233 CALL matrix_forward_transform(matrix_pp, u_pp, u_pp, &
234 cts_env%eps_filter)
235 CALL matrix_forward_transform(matrix_qq, u_qq, u_qq, &
236 cts_env%eps_filter)
237 CALL matrix_forward_transform(matrix_qp, u_qq, u_pp, &
238 cts_env%eps_filter)
239 CALL matrix_forward_transform(matrix_pq, u_pp, u_qq, &
240 cts_env%eps_filter)
241 CALL matrix_forward_transform(cts_env%matrix_x, u_qq, u_pp, &
242 cts_env%eps_filter)
243
244 IF (cts_env%max_iter .GE. 0) THEN
245
246 CALL solve_riccati_equation( &
247 pp=matrix_pp, &
248 qq=matrix_qq, &
249 qp=matrix_qp, &
250 pq=matrix_pq, &
251 x=cts_env%matrix_x, &
252 res=cts_env%matrix_res, &
253 neglect_quadratic_term=cts_env%neglect_quadratic_term, &
254 conjugator=cts_env%conjugator, &
255 max_iter=cts_env%max_iter, &
256 eps_convergence=cts_env%eps_convergence, &
257 eps_filter=cts_env%eps_filter, &
258 converged=cts_env%converged)
259
260 IF (cts_env%converged) THEN
261 !IF (unit_nr>0) THEN
262 ! WRITE(unit_nr,*)
263 ! WRITE(unit_nr,'(T6,A)') &
264 ! "RICCATI equations solved"
265 ! CALL m_flush(unit_nr)
266 !ENDIF
267 ELSE
268 cpabort("RICCATI: CG algorithm has NOT converged")
269 END IF
270
271 END IF
272
273 IF (cts_env%calculate_energy_corr) THEN
274
275 CALL dbcsr_dot(matrix_qp, cts_env%matrix_x, cts_env%energy_correction)
276
277 END IF
278
279 CALL dbcsr_release(matrix_pp)
280 CALL dbcsr_release(matrix_qp)
281 CALL dbcsr_release(matrix_qq)
282 CALL dbcsr_release(matrix_pq)
283
284 ! back-transform to the original basis
285 CALL matrix_backward_transform(cts_env%matrix_x, u_qq, &
286 u_pp, cts_env%eps_filter)
287
288 CALL dbcsr_release(u_qq)
289 CALL dbcsr_release(u_pp)
290
291 !CASE (prec_cholesky_inverse)
292 CASE (3)
293
294! RZK-warning implemented only for orthogonal tensors!!!
295! generalization to up_down should be easy
296 CALL dbcsr_create(u_pp, template=matrix_pp, &
297 matrix_type=dbcsr_type_no_symmetry)
298 CALL dbcsr_copy(u_pp, matrix_pp)
299 CALL dbcsr_scale(u_pp, -1.0_dp)
300 CALL dbcsr_add_on_diag(u_pp, &
301 abs(safety_margin*gap_estimate))
302 CALL cp_dbcsr_cholesky_decompose(u_pp, &
303 para_env=cts_env%para_env, &
304 blacs_env=cts_env%blacs_env)
305 CALL cp_dbcsr_cholesky_invert(u_pp, &
306 para_env=cts_env%para_env, &
307 blacs_env=cts_env%blacs_env, &
308 uplo_to_full=.true.)
309 !CALL dbcsr_scale(u_pp,-1.0_dp)
310
311 CALL dbcsr_create(u_qq, template=matrix_qq, &
312 matrix_type=dbcsr_type_no_symmetry)
313 CALL dbcsr_copy(u_qq, matrix_qq)
314 CALL dbcsr_add_on_diag(u_qq, &
315 abs(safety_margin*gap_estimate))
316 CALL cp_dbcsr_cholesky_decompose(u_qq, &
317 para_env=cts_env%para_env, &
318 blacs_env=cts_env%blacs_env)
319 CALL cp_dbcsr_cholesky_invert(u_qq, &
320 para_env=cts_env%para_env, &
321 blacs_env=cts_env%blacs_env, &
322 uplo_to_full=.true.)
323
324 ! transform all riccati matrices (left-right preconditioner)
325 CALL dbcsr_create(tmp1, template=matrix_qq, &
326 matrix_type=dbcsr_type_no_symmetry)
327 CALL dbcsr_multiply("N", "N", 1.0_dp, u_qq, &
328 matrix_qq, 0.0_dp, tmp1, &
329 filter_eps=cts_env%eps_filter)
330 CALL dbcsr_copy(matrix_qq, tmp1)
331 CALL dbcsr_release(tmp1)
332
333 CALL dbcsr_create(tmp1, template=matrix_pp, &
334 matrix_type=dbcsr_type_no_symmetry)
335 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_pp, &
336 u_pp, 0.0_dp, tmp1, &
337 filter_eps=cts_env%eps_filter)
338 CALL dbcsr_copy(matrix_pp, tmp1)
339 CALL dbcsr_release(tmp1)
340
341 CALL dbcsr_create(matrix_qp_save, template=matrix_qp, &
342 matrix_type=dbcsr_type_no_symmetry)
343 CALL dbcsr_copy(matrix_qp_save, matrix_qp)
344
345 CALL dbcsr_create(tmp1, template=matrix_qp, &
346 matrix_type=dbcsr_type_no_symmetry)
347 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_qp, &
348 u_pp, 0.0_dp, tmp1, &
349 filter_eps=cts_env%eps_filter)
350 CALL dbcsr_multiply("N", "N", 1.0_dp, u_qq, tmp1, &
351 0.0_dp, matrix_qp, &
352 filter_eps=cts_env%eps_filter)
353 CALL dbcsr_release(tmp1)
354!CALL dbcsr_print(matrix_qq)
355!CALL dbcsr_print(matrix_qp)
356!CALL dbcsr_print(matrix_pp)
357
358 IF (cts_env%max_iter .GE. 0) THEN
359
360 CALL solve_riccati_equation( &
361 pp=matrix_pp, &
362 qq=matrix_qq, &
363 qp=matrix_qp, &
364 pq=matrix_pq, &
365 oo=u_pp, &
366 vv=u_qq, &
367 x=cts_env%matrix_x, &
368 res=cts_env%matrix_res, &
369 neglect_quadratic_term=cts_env%neglect_quadratic_term, &
370 conjugator=cts_env%conjugator, &
371 max_iter=cts_env%max_iter, &
372 eps_convergence=cts_env%eps_convergence, &
373 eps_filter=cts_env%eps_filter, &
374 converged=cts_env%converged)
375
376 IF (cts_env%converged) THEN
377 !IF (unit_nr>0) THEN
378 ! WRITE(unit_nr,*)
379 ! WRITE(unit_nr,'(T6,A)') &
380 ! "RICCATI equations solved"
381 ! CALL m_flush(unit_nr)
382 !ENDIF
383 ELSE
384 cpabort("RICCATI: CG algorithm has NOT converged")
385 END IF
386
387 END IF
388
389 IF (cts_env%calculate_energy_corr) THEN
390
391 CALL dbcsr_dot(matrix_qp_save, cts_env%matrix_x, cts_env%energy_correction)
392
393 END IF
394 CALL dbcsr_release(matrix_qp_save)
395
396 CALL dbcsr_release(matrix_pp)
397 CALL dbcsr_release(matrix_qp)
398 CALL dbcsr_release(matrix_qq)
399 CALL dbcsr_release(matrix_pq)
400
401 CALL dbcsr_release(u_qq)
402 CALL dbcsr_release(u_pp)
403
404 CASE DEFAULT
405 cpabort("illegal preconditioner type")
406 END SELECT ! preconditioner type
407
408 IF (cts_env%update_p) THEN
409
410 IF (cts_env%tensor_type .EQ. tensor_up_down) THEN
411 cpabort("orbital update is NYI for this tensor type")
412 END IF
413
414 ! transform occupied orbitals
415 ! in a way that preserves the overlap metric
416 CALL dbcsr_create(oo1, &
417 template=cts_env%p_index_up, &
418 matrix_type=dbcsr_type_no_symmetry)
419 CALL dbcsr_create(oo1_sqrt_inv, &
420 template=oo1)
421 CALL dbcsr_create(oo1_sqrt, &
422 template=oo1)
423
424 ! Compute (1+tr(X).X)^(-1/2)_up_down
425 CALL dbcsr_multiply("T", "N", 1.0_dp, cts_env%matrix_x, &
426 cts_env%matrix_x, 0.0_dp, oo1, &
427 filter_eps=cts_env%eps_filter)
428 CALL dbcsr_add_on_diag(oo1, 1.0_dp)
429 CALL matrix_sqrt_newton_schulz(oo1_sqrt, &
430 oo1_sqrt_inv, &
431 oo1, &
432 !if cholesky is used then sqrt
433 !guess cannot be provided
434 !matrix_sqrt_inv_guess=cts_env%p_index_up,&
435 !matrix_sqrt_guess=cts_env%p_index_down,&
436 threshold=cts_env%eps_filter, &
437 order=cts_env%order_lanczos, &
438 eps_lanczos=cts_env%eps_lancsoz, &
439 max_iter_lanczos=cts_env%max_iter_lanczos)
440 CALL dbcsr_multiply("N", "N", 1.0_dp, cts_env%p_index_up, &
441 oo1_sqrt_inv, 0.0_dp, oo1, &
442 filter_eps=cts_env%eps_filter)
443 CALL dbcsr_multiply("N", "N", 1.0_dp, oo1, &
444 cts_env%p_index_down, 0.0_dp, oo1_sqrt, &
445 filter_eps=cts_env%eps_filter)
446 CALL dbcsr_release(oo1)
447 CALL dbcsr_release(oo1_sqrt_inv)
448
449 ! bring x to contravariant-covariant representation now
450 CALL dbcsr_create(matrix_qp, &
451 template=cts_env%matrix_qp_template, &
452 matrix_type=dbcsr_type_no_symmetry)
453 CALL dbcsr_multiply("N", "N", 1.0_dp, cts_env%q_index_up, &
454 cts_env%matrix_x, 0.0_dp, matrix_qp, &
455 filter_eps=cts_env%eps_filter)
456 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_qp, &
457 cts_env%p_index_down, 0.0_dp, &
458 cts_env%matrix_x, &
459 filter_eps=cts_env%eps_filter)
460 CALL dbcsr_release(matrix_qp)
461
462 ! update T=T+X or T=T+V.X (whichever is appropriate)
463 CALL dbcsr_create(t_corr, template=cts_env%matrix_t)
464 IF (cts_env%use_virt_orbs) THEN
465 CALL dbcsr_multiply("N", "N", 1.0_dp, cts_env%matrix_v, &
466 cts_env%matrix_x, 0.0_dp, t_corr, &
467 filter_eps=cts_env%eps_filter)
468 CALL dbcsr_add(cts_env%matrix_t, t_corr, &
469 1.0_dp, 1.0_dp)
470 ELSE
471 CALL dbcsr_add(cts_env%matrix_t, cts_env%matrix_x, &
472 1.0_dp, 1.0_dp)
473 END IF
474 ! adjust T so the metric is preserved: T=(T+X).(1+tr(X).X)^(-1/2)
475 CALL dbcsr_multiply("N", "N", 1.0_dp, cts_env%matrix_t, oo1_sqrt, &
476 0.0_dp, t_corr, filter_eps=cts_env%eps_filter)
477 CALL dbcsr_copy(cts_env%matrix_t, t_corr)
478
479 CALL dbcsr_release(t_corr)
480 CALL dbcsr_release(oo1_sqrt)
481
482 ELSE ! do not update p
483
484 IF (cts_env%tensor_type .EQ. tensor_orthogonal) THEN
485 ! bring x to contravariant-covariant representation
486 CALL dbcsr_create(matrix_qp, &
487 template=cts_env%matrix_qp_template, &
488 matrix_type=dbcsr_type_no_symmetry)
489 CALL dbcsr_multiply("N", "N", 1.0_dp, cts_env%q_index_up, &
490 cts_env%matrix_x, 0.0_dp, matrix_qp, &
491 filter_eps=cts_env%eps_filter)
492 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_qp, &
493 cts_env%p_index_down, 0.0_dp, &
494 cts_env%matrix_x, &
495 filter_eps=cts_env%eps_filter)
496 CALL dbcsr_release(matrix_qp)
497 END IF
498
499 END IF
500
501 ELSE
502 cpabort("illegal occ option")
503 END IF
504
505 CALL timestop(handle)
506
507 END SUBROUTINE ct_step_execute
508
509! **************************************************************************************************
510!> \brief computes oo, ov, vo, and vv blocks of the ks matrix
511!> \param ks ...
512!> \param p ...
513!> \param t ...
514!> \param v ...
515!> \param q_index_down ...
516!> \param p_index_up ...
517!> \param q_index_up ...
518!> \param pp ...
519!> \param qq ...
520!> \param qp ...
521!> \param pq ...
522!> \param tensor_type ...
523!> \param use_virt_orbs ...
524!> \param eps_filter ...
525!> \par History
526!> 2011.06 created [Rustam Z Khaliullin]
527!> \author Rustam Z Khaliullin
528! **************************************************************************************************
529 SUBROUTINE assemble_ks_qp_blocks(ks, p, t, v, q_index_down, &
530 p_index_up, q_index_up, pp, qq, qp, pq, tensor_type, use_virt_orbs, eps_filter)
531
532 TYPE(dbcsr_type), INTENT(IN) :: ks, p, t, v, q_index_down, p_index_up, &
533 q_index_up
534 TYPE(dbcsr_type), INTENT(OUT) :: pp, qq, qp, pq
535 INTEGER, INTENT(IN) :: tensor_type
536 LOGICAL, INTENT(IN) :: use_virt_orbs
537 REAL(kind=dp), INTENT(IN) :: eps_filter
538
539 CHARACTER(len=*), PARAMETER :: routinen = 'assemble_ks_qp_blocks'
540
541 INTEGER :: handle
542 LOGICAL :: library_fixed
543 TYPE(dbcsr_type) :: kst, ksv, no, on, oo, q_index_up_nosym, &
544 sp, spf, t_or, v_or
545
546 CALL timeset(routinen, handle)
547
548 IF (use_virt_orbs) THEN
549
550 ! orthogonalize the orbitals
551 CALL dbcsr_create(t_or, template=t)
552 CALL dbcsr_create(v_or, template=v)
553 CALL dbcsr_multiply("N", "N", 1.0_dp, t, p_index_up, &
554 0.0_dp, t_or, filter_eps=eps_filter)
555 CALL dbcsr_multiply("N", "N", 1.0_dp, v, q_index_up, &
556 0.0_dp, v_or, filter_eps=eps_filter)
557
558 ! KS.T
559 CALL dbcsr_create(kst, template=t)
560 CALL dbcsr_multiply("N", "N", 1.0_dp, ks, t_or, &
561 0.0_dp, kst, filter_eps=eps_filter)
562 ! pp=tr(T)*KS.T
563 CALL dbcsr_multiply("T", "N", 1.0_dp, t_or, kst, &
564 0.0_dp, pp, filter_eps=eps_filter)
565 ! qp=tr(V)*KS.T
566 CALL dbcsr_multiply("T", "N", 1.0_dp, v_or, kst, &
567 0.0_dp, qp, filter_eps=eps_filter)
568 CALL dbcsr_release(kst)
569
570 ! KS.V
571 CALL dbcsr_create(ksv, template=v)
572 CALL dbcsr_multiply("N", "N", 1.0_dp, ks, v_or, &
573 0.0_dp, ksv, filter_eps=eps_filter)
574 ! tr(T)*KS.V
575 CALL dbcsr_multiply("T", "N", 1.0_dp, t_or, ksv, &
576 0.0_dp, pq, filter_eps=eps_filter)
577 ! tr(V)*KS.V
578 CALL dbcsr_multiply("T", "N", 1.0_dp, v_or, ksv, &
579 0.0_dp, qq, filter_eps=eps_filter)
580 CALL dbcsr_release(ksv)
581
582 CALL dbcsr_release(t_or)
583 CALL dbcsr_release(v_or)
584
585 ELSE ! no virtuals, use projected AOs
586
587! THIS PROCEDURE HAS NOT BEEN UPDATED FOR CHOLESKY p/q_index_up/down
588 CALL dbcsr_create(sp, template=q_index_down, &
589 matrix_type=dbcsr_type_no_symmetry)
590 CALL dbcsr_create(spf, template=q_index_down, &
591 matrix_type=dbcsr_type_no_symmetry)
592
593 ! qp=KS*T
594 CALL dbcsr_multiply("N", "N", 1.0_dp, ks, t, 0.0_dp, qp, &
595 filter_eps=eps_filter)
596 ! pp=tr(T)*KS.T
597 CALL dbcsr_multiply("T", "N", 1.0_dp, t, qp, 0.0_dp, pp, &
598 filter_eps=eps_filter)
599 ! sp=-S_*P
600 CALL dbcsr_multiply("N", "N", -1.0_dp, q_index_down, p, 0.0_dp, sp, &
601 filter_eps=eps_filter)
602
603 ! sp=1/S^-S_.P
604 SELECT CASE (tensor_type)
605 CASE (tensor_up_down)
606 CALL dbcsr_add_on_diag(sp, 1.0_dp)
607 CASE (tensor_orthogonal)
608 CALL dbcsr_create(q_index_up_nosym, template=q_index_up, &
609 matrix_type=dbcsr_type_no_symmetry)
610 CALL dbcsr_desymmetrize(q_index_up, q_index_up_nosym)
611 CALL dbcsr_add(sp, q_index_up_nosym, 1.0_dp, 1.0_dp)
612 CALL dbcsr_release(q_index_up_nosym)
613 END SELECT
614
615 ! spf=(1/S^-S_.P)*KS
616 CALL dbcsr_multiply("N", "N", 1.0_dp, sp, ks, 0.0_dp, spf, &
617 filter_eps=eps_filter)
618
619 ! qp=spf*T
620 CALL dbcsr_multiply("N", "N", 1.0_dp, spf, t, 0.0_dp, qp, &
621 filter_eps=eps_filter)
622
623 SELECT CASE (tensor_type)
624 CASE (tensor_up_down)
625 ! pq=tr(qp)
626 CALL dbcsr_transposed(pq, qp, transpose_distribution=.false.)
627 CASE (tensor_orthogonal)
628 ! pq=sig^.tr(qp)
629 CALL dbcsr_multiply("N", "T", 1.0_dp, p_index_up, qp, 0.0_dp, pq, &
630 filter_eps=eps_filter)
631 library_fixed = .false.
632 IF (library_fixed) THEN
633 CALL dbcsr_transposed(qp, pq, transpose_distribution=.false.)
634 ELSE
635 CALL dbcsr_create(no, template=qp, &
636 matrix_type=dbcsr_type_no_symmetry)
637 CALL dbcsr_multiply("N", "N", 1.0_dp, qp, p_index_up, 0.0_dp, no, &
638 filter_eps=eps_filter)
639 CALL dbcsr_copy(qp, no)
640 CALL dbcsr_release(no)
641 END IF
642 END SELECT
643
644 ! qq=spf*tr(sp)
645 CALL dbcsr_multiply("N", "T", 1.0_dp, spf, sp, 0.0_dp, qq, &
646 filter_eps=eps_filter)
647
648 SELECT CASE (tensor_type)
649 CASE (tensor_up_down)
650
651 CALL dbcsr_create(oo, template=pp, &
652 matrix_type=dbcsr_type_no_symmetry)
653 CALL dbcsr_create(no, template=qp, &
654 matrix_type=dbcsr_type_no_symmetry)
655
656 ! first index up
657 CALL dbcsr_multiply("N", "N", 1.0_dp, q_index_up, qq, 0.0_dp, spf, &
658 filter_eps=eps_filter)
659 CALL dbcsr_copy(qq, spf)
660 CALL dbcsr_multiply("N", "N", 1.0_dp, q_index_up, qp, 0.0_dp, no, &
661 filter_eps=eps_filter)
662 CALL dbcsr_copy(qp, no)
663 CALL dbcsr_multiply("N", "N", 1.0_dp, p_index_up, pp, 0.0_dp, oo, &
664 filter_eps=eps_filter)
665 CALL dbcsr_copy(pp, oo)
666 CALL dbcsr_multiply("N", "N", 1.0_dp, p_index_up, pq, 0.0_dp, on, &
667 filter_eps=eps_filter)
668 CALL dbcsr_copy(pq, on)
669
670 CALL dbcsr_release(no)
671 CALL dbcsr_release(oo)
672
673 CASE (tensor_orthogonal)
674
675 CALL dbcsr_create(oo, template=pp, &
676 matrix_type=dbcsr_type_no_symmetry)
677
678 ! both indeces up in the pp block
679 CALL dbcsr_multiply("N", "N", 1.0_dp, p_index_up, pp, 0.0_dp, oo, &
680 filter_eps=eps_filter)
681 CALL dbcsr_multiply("N", "N", 1.0_dp, oo, p_index_up, 0.0_dp, pp, &
682 filter_eps=eps_filter)
683
684 CALL dbcsr_release(oo)
685
686 END SELECT
687
688 CALL dbcsr_release(sp)
689 CALL dbcsr_release(spf)
690
691 END IF
692
693 CALL timestop(handle)
694
695 END SUBROUTINE assemble_ks_qp_blocks
696
697! **************************************************************************************************
698!> \brief Solves the generalized Riccati or Sylvester eqation
699!> using the preconditioned conjugate gradient algorithm
700!> qp + qq.x.oo - vv.x.pp - vv.x.pq.x.oo = 0 [oo and vv are optional]
701!> qp + qq.x - x.pp - x.pq.x = 0
702!> \param pp ...
703!> \param qq ...
704!> \param qp ...
705!> \param pq ...
706!> \param oo ...
707!> \param vv ...
708!> \param x ...
709!> \param res ...
710!> \param neglect_quadratic_term ...
711!> \param conjugator ...
712!> \param max_iter ...
713!> \param eps_convergence ...
714!> \param eps_filter ...
715!> \param converged ...
716!> \par History
717!> 2011.06 created [Rustam Z Khaliullin]
718!> 2011.11 generalized [Rustam Z Khaliullin]
719!> \author Rustam Z Khaliullin
720! **************************************************************************************************
721 RECURSIVE SUBROUTINE solve_riccati_equation(pp, qq, qp, pq, oo, vv, x, res, &
722 neglect_quadratic_term, &
723 conjugator, max_iter, eps_convergence, eps_filter, &
724 converged)
725
726 TYPE(dbcsr_type), INTENT(IN) :: pp, qq
727 TYPE(dbcsr_type), INTENT(INOUT) :: qp
728 TYPE(dbcsr_type), INTENT(IN) :: pq
729 TYPE(dbcsr_type), INTENT(IN), OPTIONAL :: oo, vv
730 TYPE(dbcsr_type), INTENT(INOUT) :: x
731 TYPE(dbcsr_type), INTENT(OUT) :: res
732 LOGICAL, INTENT(IN) :: neglect_quadratic_term
733 INTEGER, INTENT(IN) :: conjugator, max_iter
734 REAL(kind=dp), INTENT(IN) :: eps_convergence, eps_filter
735 LOGICAL, INTENT(OUT) :: converged
736
737 CHARACTER(len=*), PARAMETER :: routinen = 'solve_riccati_equation'
738
739 INTEGER :: handle, istep, iteration, nsteps, &
740 unit_nr, update_prec_freq
741 LOGICAL :: prepare_to_exit, present_oo, present_vv, &
742 quadratic_term, restart_conjugator
743 REAL(kind=dp) :: best_norm, best_step_size, beta, c0, c1, &
744 c2, c3, denom, kappa, numer, &
745 obj_function, t1, t2, tau
746 REAL(kind=dp), DIMENSION(3) :: step_size
747 TYPE(cp_logger_type), POINTER :: logger
748 TYPE(dbcsr_type) :: aux1, aux2, grad, m, n, oo1, oo2, prec, &
749 res_trial, step, step_oo, vv_step
750
751!TYPE(dbcsr_type) :: qqqq, pppp, zero_pq, zero_qp
752
753 CALL timeset(routinen, handle)
754
755 logger => cp_get_default_logger()
756 IF (logger%para_env%is_source()) THEN
757 unit_nr = cp_logger_get_default_unit_nr(logger, local=.true.)
758 ELSE
759 unit_nr = -1
760 END IF
761
762 t1 = m_walltime()
763
764!IF (level.gt.5) THEN
765! CPErrorMessage(cp_failure_level,routineP,"recursion level is too high")
766! CPPrecondition(.FALSE.,cp_failure_level,routineP,failure)
767!ENDIF
768!IF (unit_nr>0) THEN
769! WRITE(unit_nr,*) &
770! "========== LEVEL ",level,"=========="
771!ENDIF
772!CALL dbcsr_print(qq)
773!CALL dbcsr_print(pp)
774!CALL dbcsr_print(qp)
775!!CALL dbcsr_print(pq)
776!IF (unit_nr>0) THEN
777! WRITE(unit_nr,*) &
778! "====== END LEVEL ",level,"=========="
779!ENDIF
780
781 quadratic_term = .NOT. neglect_quadratic_term
782 present_oo = PRESENT(oo)
783 present_vv = PRESENT(vv)
784
785 ! create aux1 matrix and init
786 CALL dbcsr_create(aux1, template=pp)
787 CALL dbcsr_copy(aux1, pp)
788 CALL dbcsr_scale(aux1, -1.0_dp)
789
790 ! create aux2 matrix and init
791 CALL dbcsr_create(aux2, template=qq)
792 CALL dbcsr_copy(aux2, qq)
793
794 ! create the gradient matrix and init
795 CALL dbcsr_create(grad, template=x)
796 CALL dbcsr_set(grad, 0.0_dp)
797
798 ! create a preconditioner
799 ! RZK-warning how to apply it to up_down tensor?
800 CALL dbcsr_create(prec, template=x)
801 !CALL create_preconditioner(prec,aux1,aux2,qp,res,tensor_type,eps_filter)
802 !CALL dbcsr_set(prec,1.0_dp)
803
804 ! create the step matrix and init
805 CALL dbcsr_create(step, template=x)
806 !CALL dbcsr_hadamard_product(prec,grad,step)
807 !CALL dbcsr_scale(step,-1.0_dp)
808
809 CALL dbcsr_create(n, template=x)
810 CALL dbcsr_create(m, template=x)
811 CALL dbcsr_create(oo1, template=pp)
812 CALL dbcsr_create(oo2, template=pp)
813 CALL dbcsr_create(res_trial, template=res)
814 CALL dbcsr_create(vv_step, template=res)
815 CALL dbcsr_create(step_oo, template=res)
816
817 ! start conjugate gradient iterations
818 iteration = 0
819 converged = .false.
820 prepare_to_exit = .false.
821 beta = 0.0_dp
822 best_step_size = 0.0_dp
823 best_norm = 1.0e+100_dp
824 !ecorr=0.0_dp
825 !change_ecorr=0.0_dp
826 restart_conjugator = .false.
827 update_prec_freq = 20
828 DO
829
830 ! (re)-compute the residuals
831 IF (iteration .EQ. 0) THEN
832 CALL dbcsr_copy(res, qp)
833 IF (present_oo) THEN
834 CALL dbcsr_multiply("N", "N", +1.0_dp, qq, x, 0.0_dp, res_trial, &
835 filter_eps=eps_filter)
836 CALL dbcsr_multiply("N", "N", +1.0_dp, res_trial, oo, 1.0_dp, res, &
837 filter_eps=eps_filter)
838 ELSE
839 CALL dbcsr_multiply("N", "N", +1.0_dp, qq, x, 1.0_dp, res, &
840 filter_eps=eps_filter)
841 END IF
842 IF (present_vv) THEN
843 CALL dbcsr_multiply("N", "N", -1.0_dp, x, pp, 0.0_dp, res_trial, &
844 filter_eps=eps_filter)
845 CALL dbcsr_multiply("N", "N", +1.0_dp, vv, res_trial, 1.0_dp, res, &
846 filter_eps=eps_filter)
847 ELSE
848 CALL dbcsr_multiply("N", "N", -1.0_dp, x, pp, 1.0_dp, res, &
849 filter_eps=eps_filter)
850 END IF
851 IF (quadratic_term) THEN
852 IF (present_oo) THEN
853 CALL dbcsr_multiply("N", "N", +1.0_dp, pq, x, 0.0_dp, oo1, &
854 filter_eps=eps_filter)
855 CALL dbcsr_multiply("N", "N", +1.0_dp, oo1, oo, 0.0_dp, oo2, &
856 filter_eps=eps_filter)
857 ELSE
858 CALL dbcsr_multiply("N", "N", +1.0_dp, pq, x, 0.0_dp, oo2, &
859 filter_eps=eps_filter)
860 END IF
861 IF (present_vv) THEN
862 CALL dbcsr_multiply("N", "N", -1.0_dp, x, oo2, 0.0_dp, res_trial, &
863 filter_eps=eps_filter)
864 CALL dbcsr_multiply("N", "N", +1.0_dp, vv, res_trial, 1.0_dp, res, &
865 filter_eps=eps_filter)
866 ELSE
867 CALL dbcsr_multiply("N", "N", -1.0_dp, x, oo2, 1.0_dp, res, &
868 filter_eps=eps_filter)
869 END IF
870 END IF
871 best_norm = dbcsr_maxabs(res)
872 ELSE
873 CALL dbcsr_add(res, m, 1.0_dp, best_step_size)
874 CALL dbcsr_add(res, n, 1.0_dp, -best_step_size*best_step_size)
875 CALL dbcsr_filter(res, eps_filter)
876 END IF
877
878 ! check convergence and other exit criteria
879 converged = (best_norm .LT. eps_convergence)
880 IF (converged .OR. (iteration .GE. max_iter)) THEN
881 prepare_to_exit = .true.
882 END IF
883
884 IF (.NOT. prepare_to_exit) THEN
885
886 ! update aux1=-pp-pq.x.oo and aux2=qq-vv.x.pq
887 IF (quadratic_term) THEN
888 IF (iteration .EQ. 0) THEN
889 IF (present_oo) THEN
890 CALL dbcsr_multiply("N", "N", -1.0_dp, pq, x, 0.0_dp, oo1, &
891 filter_eps=eps_filter)
892 CALL dbcsr_multiply("N", "N", +1.0_dp, oo1, oo, 1.0_dp, aux1, &
893 filter_eps=eps_filter)
894 ELSE
895 CALL dbcsr_multiply("N", "N", -1.0_dp, pq, x, 1.0_dp, aux1, &
896 filter_eps=eps_filter)
897 END IF
898 IF (present_vv) THEN
899 CALL dbcsr_multiply("N", "N", -1.0_dp, vv, x, 0.0_dp, res_trial, &
900 filter_eps=eps_filter)
901 CALL dbcsr_multiply("N", "N", +1.0_dp, res_trial, pq, 1.0_dp, aux2, &
902 filter_eps=eps_filter)
903 ELSE
904 CALL dbcsr_multiply("N", "N", -1.0_dp, x, pq, 1.0_dp, aux2, &
905 filter_eps=eps_filter)
906 END IF
907 ELSE
908 IF (present_oo) THEN
909 CALL dbcsr_multiply("N", "N", -best_step_size, pq, step_oo, 1.0_dp, aux1, &
910 filter_eps=eps_filter)
911 ELSE
912 CALL dbcsr_multiply("N", "N", -best_step_size, pq, step, 1.0_dp, aux1, &
913 filter_eps=eps_filter)
914 END IF
915 IF (present_vv) THEN
916 CALL dbcsr_multiply("N", "N", -best_step_size, vv_step, pq, 1.0_dp, aux2, &
917 filter_eps=eps_filter)
918 ELSE
919 CALL dbcsr_multiply("N", "N", -best_step_size, step, pq, 1.0_dp, aux2, &
920 filter_eps=eps_filter)
921 END IF
922 END IF
923 END IF
924
925 ! recompute the gradient, do not update it yet
926 ! use m matrix as a temporary storage
927 ! grad=t(vv).res.t(aux1)+t(aux2).res.t(oo)
928 IF (present_vv) THEN
929 CALL dbcsr_multiply("N", "T", 1.0_dp, res, aux1, 0.0_dp, res_trial, &
930 filter_eps=eps_filter)
931 CALL dbcsr_multiply("T", "N", 1.0_dp, vv, res_trial, 0.0_dp, m, &
932 filter_eps=eps_filter)
933 ELSE
934 CALL dbcsr_multiply("N", "T", 1.0_dp, res, aux1, 0.0_dp, m, &
935 filter_eps=eps_filter)
936 END IF
937 IF (present_oo) THEN
938 CALL dbcsr_multiply("T", "N", 1.0_dp, aux1, res, 0.0_dp, res_trial, &
939 filter_eps=eps_filter)
940 CALL dbcsr_multiply("N", "T", 1.0_dp, res_trial, oo, 1.0_dp, m, &
941 filter_eps=eps_filter)
942 ELSE
943 CALL dbcsr_multiply("T", "N", 1.0_dp, aux2, res, 1.0_dp, m, &
944 filter_eps=eps_filter)
945 END IF
946
947 ! compute preconditioner
948 !IF (iteration.eq.0.OR.(mod(iteration,update_prec_freq).eq.0)) THEN
949 IF (iteration .EQ. 0) THEN
950 CALL create_preconditioner(prec, aux1, aux2, eps_filter)
951 !restart_conjugator=.TRUE.
952!CALL dbcsr_set(prec,1.0_dp)
953!CALL dbcsr_print(prec)
954 END IF
955
956 ! compute the conjugation coefficient - beta
957 IF ((iteration .EQ. 0) .OR. restart_conjugator) THEN
958 beta = 0.0_dp
959 ELSE
960 restart_conjugator = .false.
961 SELECT CASE (conjugator)
962 CASE (cg_hestenes_stiefel)
963 CALL dbcsr_add(grad, m, -1.0_dp, 1.0_dp)
964 CALL dbcsr_hadamard_product(prec, grad, n)
965 CALL dbcsr_dot(n, m, numer)
966 CALL dbcsr_dot(grad, step, denom)
967 beta = numer/denom
968 CASE (cg_fletcher_reeves)
969 CALL dbcsr_hadamard_product(prec, grad, n)
970 CALL dbcsr_dot(grad, n, denom)
971 CALL dbcsr_hadamard_product(prec, m, n)
972 CALL dbcsr_dot(m, n, numer)
973 beta = numer/denom
974 CASE (cg_polak_ribiere)
975 CALL dbcsr_hadamard_product(prec, grad, n)
976 CALL dbcsr_dot(grad, n, denom)
977 CALL dbcsr_add(grad, m, -1.0_dp, 1.0_dp)
978 CALL dbcsr_hadamard_product(prec, grad, n)
979 CALL dbcsr_dot(n, m, numer)
980 beta = numer/denom
981 CASE (cg_fletcher)
982 CALL dbcsr_hadamard_product(prec, m, n)
983 CALL dbcsr_dot(m, n, numer)
984 CALL dbcsr_dot(grad, step, denom)
985 beta = -1.0_dp*numer/denom
986 CASE (cg_liu_storey)
987 CALL dbcsr_dot(grad, step, denom)
988 CALL dbcsr_add(grad, m, -1.0_dp, 1.0_dp)
989 CALL dbcsr_hadamard_product(prec, grad, n)
990 CALL dbcsr_dot(n, m, numer)
991 beta = -1.0_dp*numer/denom
992 CASE (cg_dai_yuan)
993 CALL dbcsr_hadamard_product(prec, m, n)
994 CALL dbcsr_dot(m, n, numer)
995 CALL dbcsr_add(grad, m, -1.0_dp, 1.0_dp)
996 CALL dbcsr_dot(grad, step, denom)
997 beta = numer/denom
998 CASE (cg_hager_zhang)
999 CALL dbcsr_add(grad, m, -1.0_dp, 1.0_dp)
1000 CALL dbcsr_dot(grad, step, denom)
1001 CALL dbcsr_hadamard_product(prec, grad, n)
1002 CALL dbcsr_dot(n, grad, numer)
1003 kappa = 2.0_dp*numer/denom
1004 CALL dbcsr_dot(n, m, numer)
1005 tau = numer/denom
1006 CALL dbcsr_dot(step, m, numer)
1007 beta = tau - kappa*numer/denom
1008 CASE (cg_zero)
1009 beta = 0.0_dp
1010 CASE DEFAULT
1011 cpabort("illegal conjugator")
1012 END SELECT
1013 END IF ! iteration.eq.0
1014
1015 ! move the current gradient to its storage
1016 CALL dbcsr_copy(grad, m)
1017
1018 ! precondition new gradient (use m as tmp storage)
1019 CALL dbcsr_hadamard_product(prec, grad, m)
1020 CALL dbcsr_filter(m, eps_filter)
1021
1022 ! recompute the step direction
1023 CALL dbcsr_add(step, m, beta, -1.0_dp)
1024 CALL dbcsr_filter(step, eps_filter)
1025
1026!! ALTERNATIVE METHOD TO OBTAIN THE STEP FROM THE GRADIENT
1027!CALL dbcsr_init(qqqq)
1028!CALL dbcsr_create(qqqq,template=qq)
1029!CALL dbcsr_init(pppp)
1030!CALL dbcsr_create(pppp,template=pp)
1031!CALL dbcsr_init(zero_pq)
1032!CALL dbcsr_create(zero_pq,template=pq)
1033!CALL dbcsr_init(zero_qp)
1034!CALL dbcsr_create(zero_qp,template=qp)
1035!CALL dbcsr_multiply("T","N",1.0_dp,aux2,aux2,0.0_dp,qqqq,&
1036! filter_eps=eps_filter)
1037!CALL dbcsr_multiply("N","T",-1.0_dp,aux1,aux1,0.0_dp,pppp,&
1038! filter_eps=eps_filter)
1039!CALL dbcsr_set(zero_qp,0.0_dp)
1040!CALL dbcsr_set(zero_pq,0.0_dp)
1041!CALL solve_riccati_equation(pppp,qqqq,grad,zero_pq,zero_qp,zero_qp,&
1042! .TRUE.,tensor_type,&
1043! conjugator,max_iter,eps_convergence,eps_filter,&
1044! converged,level+1)
1045!CALL dbcsr_release(qqqq)
1046!CALL dbcsr_release(pppp)
1047!CALL dbcsr_release(zero_qp)
1048!CALL dbcsr_release(zero_pq)
1049
1050 ! calculate the optimal step size
1051 ! m=step.aux1+aux2.step
1052 IF (present_vv) THEN
1053 CALL dbcsr_multiply("N", "N", 1.0_dp, vv, step, 0.0_dp, vv_step, &
1054 filter_eps=eps_filter)
1055 CALL dbcsr_multiply("N", "N", 1.0_dp, vv_step, aux1, 0.0_dp, m, &
1056 filter_eps=eps_filter)
1057 ELSE
1058 CALL dbcsr_multiply("N", "N", 1.0_dp, step, aux1, 0.0_dp, m, &
1059 filter_eps=eps_filter)
1060 END IF
1061 IF (present_oo) THEN
1062 CALL dbcsr_multiply("N", "N", 1.0_dp, step, oo, 0.0_dp, step_oo, &
1063 filter_eps=eps_filter)
1064 CALL dbcsr_multiply("N", "N", 1.0_dp, aux2, step_oo, 1.0_dp, m, &
1065 filter_eps=eps_filter)
1066 ELSE
1067 CALL dbcsr_multiply("N", "N", 1.0_dp, aux2, step, 1.0_dp, m, &
1068 filter_eps=eps_filter)
1069 END IF
1070
1071 IF (quadratic_term) THEN
1072 ! n=step.pq.step
1073 IF (present_oo) THEN
1074 CALL dbcsr_multiply("N", "N", 1.0_dp, pq, step, 0.0_dp, oo1, &
1075 filter_eps=eps_filter)
1076 CALL dbcsr_multiply("N", "N", 1.0_dp, oo1, oo, 0.0_dp, oo2, &
1077 filter_eps=eps_filter)
1078 ELSE
1079 CALL dbcsr_multiply("N", "N", 1.0_dp, pq, step, 0.0_dp, oo2, &
1080 filter_eps=eps_filter)
1081 END IF
1082 IF (present_vv) THEN
1083 CALL dbcsr_multiply("N", "N", 1.0_dp, step, oo2, 0.0_dp, res_trial, &
1084 filter_eps=eps_filter)
1085 CALL dbcsr_multiply("N", "N", 1.0_dp, vv, res_trial, 0.0_dp, n, &
1086 filter_eps=eps_filter)
1087 ELSE
1088 CALL dbcsr_multiply("N", "N", 1.0_dp, step, oo2, 0.0_dp, n, &
1089 filter_eps=eps_filter)
1090 END IF
1091
1092 ELSE
1093 CALL dbcsr_set(n, 0.0_dp)
1094 END IF
1095
1096 ! calculate coefficients of the cubic eq for alpha - step size
1097 c0 = 2.0_dp*(dbcsr_frobenius_norm(n))**2
1098
1099 CALL dbcsr_dot(m, n, c1)
1100 c1 = -3.0_dp*c1
1101
1102 CALL dbcsr_dot(res, n, c2)
1103 c2 = -2.0_dp*c2 + (dbcsr_frobenius_norm(m))**2
1104
1105 CALL dbcsr_dot(res, m, c3)
1106
1107 ! find step size
1108 CALL analytic_line_search(c0, c1, c2, c3, step_size, nsteps)
1109
1110 IF (nsteps .EQ. 0) THEN
1111 cpabort("no step sizes!")
1112 END IF
1113 ! if we have several possible step sizes
1114 ! choose one with the lowest objective function
1115 best_norm = 1.0e+100_dp
1116 best_step_size = 0.0_dp
1117 DO istep = 1, nsteps
1118 ! recompute the residues
1119 CALL dbcsr_copy(res_trial, res)
1120 CALL dbcsr_add(res_trial, m, 1.0_dp, step_size(istep))
1121 CALL dbcsr_add(res_trial, n, 1.0_dp, -step_size(istep)*step_size(istep))
1122 CALL dbcsr_filter(res_trial, eps_filter)
1123 ! RZK-warning objective function might be different in the case of
1124 ! tensor_up_down
1125 !obj_function=0.5_dp*(dbcsr_frobenius_norm(res_trial))**2
1126 obj_function = dbcsr_maxabs(res_trial)
1127 IF (obj_function .LT. best_norm) THEN
1128 best_norm = obj_function
1129 best_step_size = step_size(istep)
1130 END IF
1131 END DO
1132
1133 END IF
1134
1135 ! update X along the line
1136 CALL dbcsr_add(x, step, 1.0_dp, best_step_size)
1137 CALL dbcsr_filter(x, eps_filter)
1138
1139 ! evaluate current energy correction
1140 !change_ecorr=ecorr
1141 !CALL dbcsr_dot(qp,x,ecorr,"T","N")
1142 !change_ecorr=ecorr-change_ecorr
1143
1144 ! check convergence and other exit criteria
1145 converged = (best_norm .LT. eps_convergence)
1146 IF (converged .OR. (iteration .GE. max_iter)) THEN
1147 prepare_to_exit = .true.
1148 END IF
1149
1150 t2 = m_walltime()
1151
1152 IF (unit_nr > 0) THEN
1153 WRITE (unit_nr, '(T6,A,1X,I4,1X,E12.3,F8.3)') &
1154 "RICCATI iter ", iteration, best_norm, t2 - t1
1155 !WRITE(unit_nr,'(T6,A,1X,I4,1X,F15.9,F15.9,E12.3,F8.3)') &
1156 ! "RICCATI iter ",iteration,ecorr,change_ecorr,best_norm,t2-t1
1157 END IF
1158
1159 t1 = m_walltime()
1160
1161 iteration = iteration + 1
1162
1163 IF (prepare_to_exit) EXIT
1164
1165 END DO
1166
1167 CALL dbcsr_release(aux1)
1168 CALL dbcsr_release(aux2)
1169 CALL dbcsr_release(grad)
1170 CALL dbcsr_release(step)
1171 CALL dbcsr_release(n)
1172 CALL dbcsr_release(m)
1173 CALL dbcsr_release(oo1)
1174 CALL dbcsr_release(oo2)
1175 CALL dbcsr_release(res_trial)
1176 CALL dbcsr_release(vv_step)
1177 CALL dbcsr_release(step_oo)
1178
1179 CALL timestop(handle)
1180
1181 END SUBROUTINE solve_riccati_equation
1182
1183! **************************************************************************************************
1184!> \brief Computes a preconditioner from diagonal elements of ~f_oo, ~f_vv
1185!> The preconditioner is approximately equal to
1186!> prec_ai ~ (e_a - e_i)^(-2)
1187!> However, the real expression is more complex
1188!> \param prec ...
1189!> \param pp ...
1190!> \param qq ...
1191!> \param eps_filter ...
1192!> \par History
1193!> 2011.07 created [Rustam Z Khaliullin]
1194!> \author Rustam Z Khaliullin
1195! **************************************************************************************************
1196 SUBROUTINE create_preconditioner(prec, pp, qq, eps_filter)
1197
1198 TYPE(dbcsr_type), INTENT(OUT) :: prec
1199 TYPE(dbcsr_type), INTENT(IN) :: pp, qq
1200 REAL(kind=dp), INTENT(IN) :: eps_filter
1201
1202 CHARACTER(len=*), PARAMETER :: routinen = 'create_preconditioner'
1203
1204 INTEGER :: handle, p_nrows, q_nrows
1205 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: p_diagonal, q_diagonal
1206 REAL(kind=dp), DIMENSION(:, :), POINTER :: block
1207 TYPE(dbcsr_iterator_type) :: iter
1208 TYPE(dbcsr_type) :: pp_diag, qq_diag, t1, t2, tmp
1209
1210!LOGICAL, INTENT(IN) :: use_virt_orbs
1211
1212 CALL timeset(routinen, handle)
1213
1214! ! copy diagonal elements
1215! CALL dbcsr_get_info(pp,nfullrows_total=nrows)
1216! CALL dbcsr_init(pp_diag)
1217! CALL dbcsr_create(pp_diag,template=pp)
1218! ALLOCATE(diagonal(nrows))
1219! CALL dbcsr_get_diag(pp,diagonal)
1220! CALL dbcsr_add_on_diag(pp_diag,1.0_dp)
1221! CALL dbcsr_set_diag(pp_diag,diagonal)
1222! DEALLOCATE(diagonal)
1223!
1224 ! initialize a matrix to 1.0
1225 CALL dbcsr_create(tmp, template=prec)
1226 CALL dbcsr_reserve_diag_blocks(tmp)
1227 CALL dbcsr_iterator_start(iter, tmp)
1228 DO WHILE (dbcsr_iterator_blocks_left(iter))
1229 CALL dbcsr_iterator_next_block(iter, block=block)
1230 block(:, :) = 1.0_dp
1231 END DO
1232 CALL dbcsr_iterator_stop(iter)
1233
1234 ! copy diagonal elements of pp into cols of a matrix
1235 CALL dbcsr_get_info(pp, nfullrows_total=p_nrows)
1236 CALL dbcsr_create(pp_diag, template=pp)
1237 ALLOCATE (p_diagonal(p_nrows))
1238 CALL dbcsr_get_diag(pp, p_diagonal)
1239 CALL dbcsr_add_on_diag(pp_diag, 1.0_dp)
1240 CALL dbcsr_set_diag(pp_diag, p_diagonal)
1241 ! RZK-warning is it possible to use dbcsr_scale_by_vector?
1242 ! or even insert elements directly in the prev cycles
1243 CALL dbcsr_create(t2, template=prec)
1244 CALL dbcsr_multiply("N", "N", 1.0_dp, tmp, pp_diag, &
1245 0.0_dp, t2, filter_eps=eps_filter)
1246
1247 ! copy diagonal elements qq into rows of a matrix
1248 CALL dbcsr_get_info(qq, nfullrows_total=q_nrows)
1249 CALL dbcsr_create(qq_diag, template=qq)
1250 ALLOCATE (q_diagonal(q_nrows))
1251 CALL dbcsr_get_diag(qq, q_diagonal)
1252 CALL dbcsr_add_on_diag(qq_diag, 1.0_dp)
1253 CALL dbcsr_set_diag(qq_diag, q_diagonal)
1254 CALL dbcsr_set(tmp, 1.0_dp)
1255 CALL dbcsr_create(t1, template=prec)
1256 CALL dbcsr_multiply("N", "N", 1.0_dp, qq_diag, tmp, &
1257 0.0_dp, t1, filter_eps=eps_filter)
1258
1259 CALL dbcsr_hadamard_product(t1, t2, prec)
1260 CALL dbcsr_release(t1)
1261 CALL dbcsr_scale(prec, 2.0_dp)
1262
1263 ! Get the diagonal of tr(qq).qq
1264 CALL dbcsr_multiply("T", "N", 1.0_dp, qq, qq, &
1265 0.0_dp, qq_diag, retain_sparsity=.true., &
1266 filter_eps=eps_filter)
1267 CALL dbcsr_get_diag(qq_diag, q_diagonal)
1268 CALL dbcsr_set(qq_diag, 0.0_dp)
1269 CALL dbcsr_add_on_diag(qq_diag, 1.0_dp)
1270 CALL dbcsr_set_diag(qq_diag, q_diagonal)
1271 DEALLOCATE (q_diagonal)
1272 CALL dbcsr_set(tmp, 1.0_dp)
1273 CALL dbcsr_multiply("N", "N", 1.0_dp, qq_diag, tmp, &
1274 0.0_dp, t2, filter_eps=eps_filter)
1275 CALL dbcsr_release(qq_diag)
1276 CALL dbcsr_add(prec, t2, 1.0_dp, 1.0_dp)
1277
1278 ! Get the diagonal of pp.tr(pp)
1279 CALL dbcsr_multiply("N", "T", 1.0_dp, pp, pp, &
1280 0.0_dp, pp_diag, retain_sparsity=.true., &
1281 filter_eps=eps_filter)
1282 CALL dbcsr_get_diag(pp_diag, p_diagonal)
1283 CALL dbcsr_set(pp_diag, 0.0_dp)
1284 CALL dbcsr_add_on_diag(pp_diag, 1.0_dp)
1285 CALL dbcsr_set_diag(pp_diag, p_diagonal)
1286 DEALLOCATE (p_diagonal)
1287 CALL dbcsr_set(tmp, 1.0_dp)
1288 CALL dbcsr_multiply("N", "N", 1.0_dp, tmp, pp_diag, &
1289 0.0_dp, t2, filter_eps=eps_filter)
1290 CALL dbcsr_release(tmp)
1291 CALL dbcsr_release(pp_diag)
1292 CALL dbcsr_add(prec, t2, 1.0_dp, 1.0_dp)
1293
1294 ! now add the residual component
1295 !CALL dbcsr_hadamard_product(res,qp,t2)
1296 !CALL dbcsr_add(prec,t2,1.0_dp,-2.0_dp)
1297 CALL dbcsr_release(t2)
1298 CALL inverse_of_elements(prec)
1299 CALL dbcsr_filter(prec, eps_filter)
1300
1301 CALL timestop(handle)
1302
1303 END SUBROUTINE create_preconditioner
1304
1305! **************************************************************************************************
1306!> \brief Computes 1/x of the matrix elements.
1307!> \param matrix ...
1308!> \author Ole Schuett
1309! **************************************************************************************************
1310 SUBROUTINE inverse_of_elements(matrix)
1311 TYPE(dbcsr_type), INTENT(INOUT) :: matrix
1312
1313 CHARACTER(len=*), PARAMETER :: routinen = 'inverse_of_elements'
1314
1315 INTEGER :: handle
1316 REAL(kind=dp), DIMENSION(:, :), POINTER :: block
1317 TYPE(dbcsr_iterator_type) :: iter
1318
1319 CALL timeset(routinen, handle)
1320 CALL dbcsr_iterator_start(iter, matrix)
1321 DO WHILE (dbcsr_iterator_blocks_left(iter))
1322 CALL dbcsr_iterator_next_block(iter, block=block)
1323 block = 1.0_dp/block
1324 END DO
1325 CALL dbcsr_iterator_stop(iter)
1326 CALL timestop(handle)
1327
1328 END SUBROUTINE inverse_of_elements
1329
1330! **************************************************************************************************
1331!> \brief Finds real roots of a cubic equation
1332!> > a*x**3 + b*x**2 + c*x + d = 0
1333!> and returns only those roots for which the derivative is positive
1334!>
1335!> Step 0: Check the true order of the equation. Cubic, quadratic, linear?
1336!> Step 1: Calculate p and q
1337!> p = ( 3*c/a - (b/a)**2 ) / 3
1338!> q = ( 2*(b/a)**3 - 9*b*c/a/a + 27*d/a ) / 27
1339!> Step 2: Calculate discriminant D
1340!> D = (p/3)**3 + (q/2)**2
1341!> Step 3: Depending on the sign of D, we follow different strategy.
1342!> If D<0, three distinct real roots.
1343!> If D=0, three real roots of which at least two are equal.
1344!> If D>0, one real and two complex roots.
1345!> Step 3a: For D>0 and D=0,
1346!> Calculate u and v
1347!> u = cubic_root(-q/2 + sqrt(D))
1348!> v = cubic_root(-q/2 - sqrt(D))
1349!> Find the three transformed roots
1350!> y1 = u + v
1351!> y2 = -(u+v)/2 + i (u-v)*sqrt(3)/2
1352!> y3 = -(u+v)/2 - i (u-v)*sqrt(3)/2
1353!> Step 3b Alternately, for D<0, a trigonometric formulation is more convenient
1354!> y1 = 2 * sqrt(|p|/3) * cos(phi/3)
1355!> y2 = -2 * sqrt(|p|/3) * cos((phi+pi)/3)
1356!> y3 = -2 * sqrt(|p|/3) * cos((phi-pi)/3)
1357!> where phi = acos(-q/2/sqrt(|p|**3/27))
1358!> pi = 3.141592654...
1359!> Step 4 Find the real roots
1360!> x = y - b/a/3
1361!> Step 5 Check the derivative and return only those real roots
1362!> for which the derivative is positive
1363!>
1364!> \param a ...
1365!> \param b ...
1366!> \param c ...
1367!> \param d ...
1368!> \param minima ...
1369!> \param nmins ...
1370!> \par History
1371!> 2011.06 created [Rustam Z Khaliullin]
1372!> \author Rustam Z Khaliullin
1373! **************************************************************************************************
1374 SUBROUTINE analytic_line_search(a, b, c, d, minima, nmins)
1375
1376 REAL(kind=dp), INTENT(IN) :: a, b, c, d
1377 REAL(kind=dp), DIMENSION(3), INTENT(OUT) :: minima
1378 INTEGER, INTENT(OUT) :: nmins
1379
1380 INTEGER :: i, nroots
1381 REAL(kind=dp) :: dd, der, p, phi, pi, q, temp1, temp2, u, &
1382 v, y1, y2, y2i, y2r, y3
1383 REAL(kind=dp), DIMENSION(3) :: x
1384
1385! CALL timeset(routineN,handle)
1386
1387 pi = acos(-1.0_dp)
1388
1389 ! Step 0: Check coefficients and find the true order of the eq
1390 IF (a .EQ. 0.0_dp) THEN
1391 IF (b .EQ. 0.0_dp) THEN
1392 IF (c .EQ. 0.0_dp) THEN
1393 ! Non-equation, no valid solutions
1394 nroots = 0
1395 ELSE
1396 ! Linear equation with one root.
1397 nroots = 1
1398 x(1) = -d/c
1399 END IF
1400 ELSE
1401 ! Quadratic equation with max two roots.
1402 dd = c*c - 4.0_dp*b*d
1403 IF (dd .GT. 0.0_dp) THEN
1404 nroots = 2
1405 x(1) = (-c + sqrt(dd))/2.0_dp/b
1406 x(2) = (-c - sqrt(dd))/2.0_dp/b
1407 ELSE IF (dd .LT. 0.0_dp) THEN
1408 nroots = 0
1409 ELSE
1410 nroots = 1
1411 x(1) = -c/2.0_dp/b
1412 END IF
1413 END IF
1414 ELSE
1415 ! Cubic equation with max three roots
1416 ! Calculate p and q
1417 p = c/a - b*b/a/a/3.0_dp
1418 q = (2.0_dp*b*b*b/a/a/a - 9.0_dp*b*c/a/a + 27.0_dp*d/a)/27.0_dp
1419
1420 ! Calculate DD
1421 dd = p*p*p/27.0_dp + q*q/4.0_dp
1422
1423 IF (dd .LT. 0.0_dp) THEN
1424 ! three real unequal roots -- use the trigonometric formulation
1425 phi = acos(-q/2.0_dp/sqrt(abs(p*p*p)/27.0_dp))
1426 temp1 = 2.0_dp*sqrt(abs(p)/3.0_dp)
1427 y1 = temp1*cos(phi/3.0_dp)
1428 y2 = -temp1*cos((phi + pi)/3.0_dp)
1429 y3 = -temp1*cos((phi - pi)/3.0_dp)
1430 ELSE
1431 ! 1 real & 2 conjugate complex roots OR 3 real roots (some are equal)
1432 temp1 = -q/2.0_dp + sqrt(dd)
1433 temp2 = -q/2.0_dp - sqrt(dd)
1434 u = abs(temp1)**(1.0_dp/3.0_dp)
1435 v = abs(temp2)**(1.0_dp/3.0_dp)
1436 IF (temp1 .LT. 0.0_dp) u = -u
1437 IF (temp2 .LT. 0.0_dp) v = -v
1438 y1 = u + v
1439 y2r = -(u + v)/2.0_dp
1440 y2i = (u - v)*sqrt(3.0_dp)/2.0_dp
1441 END IF
1442
1443 ! Final transformation
1444 temp1 = b/a/3.0_dp
1445 y1 = y1 - temp1
1446 y2 = y2 - temp1
1447 y3 = y3 - temp1
1448 y2r = y2r - temp1
1449
1450 ! Assign answers
1451 IF (dd .LT. 0.0_dp) THEN
1452 nroots = 3
1453 x(1) = y1
1454 x(2) = y2
1455 x(3) = y3
1456 ELSE IF (dd .EQ. 0.0_dp) THEN
1457 nroots = 2
1458 x(1) = y1
1459 x(2) = y2r
1460 !x(3) = cmplx(y2r, 0.)
1461 ELSE
1462 nroots = 1
1463 x(1) = y1
1464 !x(2) = cmplx(y2r, y2i)
1465 !x(3) = cmplx(y2r,-y2i)
1466 END IF
1467
1468 END IF
1469
1470!write(*,'(i2,a)') nroots, ' real root(s)'
1471 nmins = 0
1472 DO i = 1, nroots
1473 ! maximum or minimum? use the derivative
1474 ! 3*a*x**2+2*b*x+c
1475 der = 3.0_dp*a*x(i)*x(i) + 2.0_dp*b*x(i) + c
1476 IF (der .GT. 0.0_dp) THEN
1477 nmins = nmins + 1
1478 minima(nmins) = x(i)
1479!write(*,'(a,i2,a,f10.5)') 'Minimum ', i, ', value: ', x(i)
1480 END IF
1481 END DO
1482
1483! CALL timestop(handle)
1484
1485 END SUBROUTINE analytic_line_search
1486
1487! **************************************************************************************************
1488!> \brief Diagonalizes diagonal blocks of a symmetric dbcsr matrix
1489!> and returs its eigenvectors
1490!> \param matrix ...
1491!> \param c ...
1492!> \param e ...
1493!> \par History
1494!> 2011.07 created [Rustam Z Khaliullin]
1495!> \author Rustam Z Khaliullin
1496! **************************************************************************************************
1497 SUBROUTINE diagonalize_diagonal_blocks(matrix, c, e)
1498
1499 TYPE(dbcsr_type), INTENT(IN) :: matrix
1500 TYPE(dbcsr_type), INTENT(OUT) :: c
1501 TYPE(dbcsr_type), INTENT(OUT), OPTIONAL :: e
1502
1503 CHARACTER(len=*), PARAMETER :: routinen = 'diagonalize_diagonal_blocks'
1504
1505 INTEGER :: handle, iblock_col, iblock_row, &
1506 iblock_size, info, lwork, orbital
1507 LOGICAL :: block_needed, do_eigenvalues
1508 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: eigenvalues, work
1509 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: data_copy, new_block
1510 REAL(kind=dp), DIMENSION(:, :), POINTER :: data_p
1511 TYPE(dbcsr_iterator_type) :: iter
1512
1513 CALL timeset(routinen, handle)
1514
1515 IF (PRESENT(e)) THEN
1516 do_eigenvalues = .true.
1517 ELSE
1518 do_eigenvalues = .false.
1519 END IF
1520
1521 ! create a matrix for eigenvectors
1522 CALL dbcsr_work_create(c, work_mutable=.true.)
1523 IF (do_eigenvalues) &
1524 CALL dbcsr_work_create(e, work_mutable=.true.)
1525
1526 CALL dbcsr_iterator_readonly_start(iter, matrix)
1527
1528 DO WHILE (dbcsr_iterator_blocks_left(iter))
1529
1530 CALL dbcsr_iterator_next_block(iter, iblock_row, iblock_col, data_p, row_size=iblock_size)
1531
1532 block_needed = .false.
1533 IF (iblock_row == iblock_col) block_needed = .true.
1534
1535 IF (block_needed) THEN
1536
1537 ! Prepare data
1538 ALLOCATE (eigenvalues(iblock_size))
1539 ALLOCATE (data_copy(iblock_size, iblock_size))
1540 data_copy(:, :) = data_p(:, :)
1541
1542 ! Query the optimal workspace for dsyev
1543 lwork = -1
1544 ALLOCATE (work(max(1, lwork)))
1545 CALL dsyev('V', 'L', iblock_size, data_copy, iblock_size, eigenvalues, work, lwork, info)
1546 lwork = int(work(1))
1547 DEALLOCATE (work)
1548
1549 ! Allocate the workspace and solve the eigenproblem
1550 ALLOCATE (work(max(1, lwork)))
1551 CALL dsyev('V', 'L', iblock_size, data_copy, iblock_size, eigenvalues, work, lwork, info)
1552 IF (info .NE. 0) cpabort("DSYEV failed")
1553
1554 ! copy eigenvectors into a cp_dbcsr matrix
1555 CALL dbcsr_put_block(c, iblock_row, iblock_col, block=data_copy)
1556
1557 ! if requested copy eigenvalues into a cp_dbcsr matrix
1558 IF (do_eigenvalues) THEN
1559 ALLOCATE (new_block(iblock_size, iblock_size))
1560 new_block(:, :) = 0.0_dp
1561 DO orbital = 1, iblock_size
1562 new_block(orbital, orbital) = eigenvalues(orbital)
1563 END DO
1564 CALL dbcsr_put_block(e, iblock_row, iblock_col, new_block)
1565 DEALLOCATE (new_block)
1566 END IF
1567
1568 DEALLOCATE (work)
1569 DEALLOCATE (data_copy)
1570 DEALLOCATE (eigenvalues)
1571
1572 END IF
1573
1574 END DO
1575
1576 CALL dbcsr_iterator_stop(iter)
1577
1578 CALL dbcsr_finalize(c)
1579 IF (do_eigenvalues) CALL dbcsr_finalize(e)
1580
1581 CALL timestop(handle)
1582
1583 END SUBROUTINE diagonalize_diagonal_blocks
1584
1585! **************************************************************************************************
1586!> \brief Transforms a matrix M_out = tr(U1) * M_in * U2
1587!> \param matrix ...
1588!> \param u1 ...
1589!> \param u2 ...
1590!> \param eps_filter ...
1591!> \par History
1592!> 2011.10 created [Rustam Z Khaliullin]
1593!> \author Rustam Z Khaliullin
1594! **************************************************************************************************
1595 SUBROUTINE matrix_forward_transform(matrix, u1, u2, eps_filter)
1596
1597 TYPE(dbcsr_type), INTENT(INOUT) :: matrix
1598 TYPE(dbcsr_type), INTENT(IN) :: u1, u2
1599 REAL(kind=dp), INTENT(IN) :: eps_filter
1600
1601 CHARACTER(len=*), PARAMETER :: routinen = 'matrix_forward_transform'
1602
1603 INTEGER :: handle
1604 TYPE(dbcsr_type) :: tmp
1605
1606 CALL timeset(routinen, handle)
1607
1608 CALL dbcsr_create(tmp, template=matrix, &
1609 matrix_type=dbcsr_type_no_symmetry)
1610 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix, u2, 0.0_dp, tmp, &
1611 filter_eps=eps_filter)
1612 CALL dbcsr_multiply("T", "N", 1.0_dp, u1, tmp, 0.0_dp, matrix, &
1613 filter_eps=eps_filter)
1614 CALL dbcsr_release(tmp)
1615
1616 CALL timestop(handle)
1617
1618 END SUBROUTINE matrix_forward_transform
1619
1620! **************************************************************************************************
1621!> \brief Transforms a matrix M_out = U1 * M_in * tr(U2)
1622!> \param matrix ...
1623!> \param u1 ...
1624!> \param u2 ...
1625!> \param eps_filter ...
1626!> \par History
1627!> 2011.10 created [Rustam Z Khaliullin]
1628!> \author Rustam Z Khaliullin
1629! **************************************************************************************************
1630 SUBROUTINE matrix_backward_transform(matrix, u1, u2, eps_filter)
1631
1632 TYPE(dbcsr_type), INTENT(INOUT) :: matrix
1633 TYPE(dbcsr_type), INTENT(IN) :: u1, u2
1634 REAL(kind=dp), INTENT(IN) :: eps_filter
1635
1636 CHARACTER(len=*), PARAMETER :: routinen = 'matrix_backward_transform'
1637
1638 INTEGER :: handle
1639 TYPE(dbcsr_type) :: tmp
1640
1641 CALL timeset(routinen, handle)
1642
1643 CALL dbcsr_create(tmp, template=matrix, &
1644 matrix_type=dbcsr_type_no_symmetry)
1645 CALL dbcsr_multiply("N", "T", 1.0_dp, matrix, u2, 0.0_dp, tmp, &
1646 filter_eps=eps_filter)
1647 CALL dbcsr_multiply("N", "N", 1.0_dp, u1, tmp, 0.0_dp, matrix, &
1648 filter_eps=eps_filter)
1649 CALL dbcsr_release(tmp)
1650
1651 CALL timestop(handle)
1652
1653 END SUBROUTINE matrix_backward_transform
1654
1655!! **************************************************************************************************
1656!!> \brief Transforms to a representation in which diagonal blocks
1657!!> of qq and pp matrices are diagonal. This can improve convergence
1658!!> of PCG
1659!!> \par History
1660!!> 2011.07 created [Rustam Z Khaliullin]
1661!!> \author Rustam Z Khaliullin
1662!! **************************************************************************************************
1663! SUBROUTINE transform_matrices_to_blk_diag(matrix_pp,matrix_qq,matrix_qp,&
1664! matrix_pq,eps_filter)
1665!
1666! TYPE(dbcsr_type), INTENT(INOUT) :: matrix_pp, matrix_qq,&
1667! matrix_qp, matrix_pq
1668! REAL(KIND=dp), INTENT(IN) :: eps_filter
1669!
1670! CHARACTER(len=*), PARAMETER :: routineN = 'transform_matrices_to_blk_diag',&
1671! routineP = moduleN//':'//routineN
1672!
1673! TYPE(dbcsr_type) :: tmp_pp, tmp_qq,&
1674! tmp_qp, tmp_pq,&
1675! blk, blk2
1676! INTEGER :: handle
1677!
1678! CALL timeset(routineN,handle)
1679!
1680! ! find a better basis by diagonalizing diagonal blocks
1681! ! first pp
1682! CALL dbcsr_init(blk)
1683! CALL dbcsr_create(blk,template=matrix_pp)
1684! CALL diagonalize_diagonal_blocks(matrix_pp,blk)
1685!
1686! ! convert matrices to the new basis
1687! CALL dbcsr_init(tmp_pp)
1688! CALL dbcsr_create(tmp_pp,template=matrix_pp)
1689! CALL dbcsr_multiply("N","N",1.0_dp,matrix_pp,blk,0.0_dp,tmp_pp,&
1690! filter_eps=eps_filter)
1691! CALL dbcsr_multiply("T","N",1.0_dp,blk,tmp_pp,0.0_dp,matrix_pp,&
1692! filter_eps=eps_filter)
1693! CALL dbcsr_release(tmp_pp)
1694!
1695! ! now qq
1696! CALL dbcsr_init(blk2)
1697! CALL dbcsr_create(blk2,template=matrix_qq)
1698! CALL diagonalize_diagonal_blocks(matrix_qq,blk2)
1699!
1700! CALL dbcsr_init(tmp_qq)
1701! CALL dbcsr_create(tmp_qq,template=matrix_qq)
1702! CALL dbcsr_multiply("N","N",1.0_dp,matrix_qq,blk2,0.0_dp,tmp_qq,&
1703! filter_eps=eps_filter)
1704! CALL dbcsr_multiply("T","N",1.0_dp,blk2,tmp_qq,0.0_dp,matrix_qq,&
1705! filter_eps=eps_filter)
1706! CALL dbcsr_release(tmp_qq)
1707!
1708! ! transform pq
1709! CALL dbcsr_init(tmp_pq)
1710! CALL dbcsr_create(tmp_pq,template=matrix_pq)
1711! CALL dbcsr_multiply("T","N",1.0_dp,blk,matrix_pq,0.0_dp,tmp_pq,&
1712! filter_eps=eps_filter)
1713! CALL dbcsr_multiply("N","N",1.0_dp,tmp_pq,blk2,0.0_dp,matrix_pq,&
1714! filter_eps=eps_filter)
1715! CALL dbcsr_release(tmp_pq)
1716!
1717! ! transform qp
1718! CALL dbcsr_init(tmp_qp)
1719! CALL dbcsr_create(tmp_qp,template=matrix_qp)
1720! CALL dbcsr_multiply("N","N",1.0_dp,matrix_qp,blk,0.0_dp,tmp_qp,&
1721! filter_eps=eps_filter)
1722! CALL dbcsr_multiply("T","N",1.0_dp,blk2,tmp_qp,0.0_dp,matrix_qp,&
1723! filter_eps=eps_filter)
1724! CALL dbcsr_release(tmp_qp)
1725!
1726! CALL dbcsr_release(blk2)
1727! CALL dbcsr_release(blk)
1728!
1729! CALL timestop(handle)
1730!
1731! END SUBROUTINE transform_matrices_to_blk_diag
1732
1733! **************************************************************************************************
1734!> \brief computes oo, ov, vo, and vv blocks of the ks matrix
1735!> \par History
1736!> 2011.06 created [Rustam Z Khaliullin]
1737!> \author Rustam Z Khaliullin
1738! **************************************************************************************************
1739! SUBROUTINE ct_step_env_execute(env)
1740!
1741! TYPE(ct_step_env_type) :: env
1742!
1743! CHARACTER(len=*), PARAMETER :: routineN = 'ct_step_env_execute', &
1744! routineP = moduleN//':'//routineN
1745!
1746! INTEGER :: handle
1747!
1748! CALL timeset(routineN,handle)
1749!
1750!
1751! CALL timestop(handle)
1752!
1753! END SUBROUTINE ct_step_env_execute
1754
1755END MODULE ct_methods
1756
cp_dbcsr_api::dbcsr_create
Definition cp_dbcsr_api.F:192
cp_dbcsr_api
Definition cp_dbcsr_api.F:8
cp_dbcsr_api::dbcsr_transposed
subroutine, public dbcsr_transposed(transposed, normal, shallow_data_copy, transpose_distribution, use_distribution)
...
Definition cp_dbcsr_api.F:1219
cp_dbcsr_api::dbcsr_scale
subroutine, public dbcsr_scale(matrix, alpha_scalar)
...
Definition cp_dbcsr_api.F:1165
cp_dbcsr_api::dbcsr_iterator_next_block
subroutine, public dbcsr_iterator_next_block(iterator, row, column, block, block_number_argument_has_been_removed, row_size, col_size, row_offset, col_offset)
...
Definition cp_dbcsr_api.F:956
cp_dbcsr_api::dbcsr_iterator_blocks_left
logical function, public dbcsr_iterator_blocks_left(iterator)
...
Definition cp_dbcsr_api.F:930
cp_dbcsr_api::dbcsr_iterator_stop
subroutine, public dbcsr_iterator_stop(iterator)
...
Definition cp_dbcsr_api.F:1027
cp_dbcsr_api::dbcsr_desymmetrize
subroutine, public dbcsr_desymmetrize(matrix_a, matrix_b)
...
Definition cp_dbcsr_api.F:517
cp_dbcsr_api::dbcsr_copy
subroutine, public dbcsr_copy(matrix_b, matrix_a, name, keep_sparsity, keep_imaginary)
...
Definition cp_dbcsr_api.F:368
cp_dbcsr_api::dbcsr_multiply
subroutine, public dbcsr_multiply(transa, transb, alpha, matrix_a, matrix_b, beta, matrix_c, first_row, last_row, first_column, last_column, first_k, last_k, retain_sparsity, filter_eps, flop)
...
Definition cp_dbcsr_api.F:1073
cp_dbcsr_api::dbcsr_get_info
subroutine, public dbcsr_get_info(matrix, nblkrows_total, nblkcols_total, nfullrows_total, nfullcols_total, nblkrows_local, nblkcols_local, nfullrows_local, nfullcols_local, my_prow, my_pcol, local_rows, local_cols, proc_row_dist, proc_col_dist, row_blk_size, col_blk_size, row_blk_offset, col_blk_offset, distribution, name, matrix_type, group)
...
Definition cp_dbcsr_api.F:794
cp_dbcsr_api::dbcsr_work_create
subroutine, public dbcsr_work_create(matrix, nblks_guess, sizedata_guess, n, work_mutable)
...
Definition cp_dbcsr_api.F:1284
cp_dbcsr_api::dbcsr_filter
subroutine, public dbcsr_filter(matrix, eps)
...
Definition cp_dbcsr_api.F:644
cp_dbcsr_api::dbcsr_finalize
subroutine, public dbcsr_finalize(matrix)
...
Definition cp_dbcsr_api.F:659
cp_dbcsr_api::dbcsr_iterator_start
subroutine, public dbcsr_iterator_start(iterator, matrix, shared, dynamic, dynamic_byrows)
...
Definition cp_dbcsr_api.F:989
cp_dbcsr_api::dbcsr_set
subroutine, public dbcsr_set(matrix, alpha)
...
Definition cp_dbcsr_api.F:1181
cp_dbcsr_api::dbcsr_release
subroutine, public dbcsr_release(matrix)
...
Definition cp_dbcsr_api.F:1119
cp_dbcsr_api::dbcsr_iterator_readonly_start
subroutine, public dbcsr_iterator_readonly_start(iterator, matrix, shared, dynamic, dynamic_byrows)
Like dbcsr_iterator_start() but with matrix being INTENT(IN). When invoking this routine,...
Definition cp_dbcsr_api.F:1010
cp_dbcsr_api::dbcsr_put_block
subroutine, public dbcsr_put_block(matrix, row, col, block, summation)
...
Definition cp_dbcsr_api.F:1102
cp_dbcsr_api::dbcsr_add
subroutine, public dbcsr_add(matrix_a, matrix_b, alpha_scalar, beta_scalar)
...
Definition cp_dbcsr_api.F:251
cp_dbcsr_cholesky
Interface to (sca)lapack for the Cholesky based procedures.
Definition cp_dbcsr_cholesky.F:17
cp_dbcsr_cholesky::cp_dbcsr_cholesky_decompose
subroutine, public cp_dbcsr_cholesky_decompose(matrix, n, para_env, blacs_env)
used to replace a symmetric positive def. matrix M with its cholesky decomposition U: M = U^T * U,...
Definition cp_dbcsr_cholesky.F:61
cp_dbcsr_cholesky::cp_dbcsr_cholesky_invert
subroutine, public cp_dbcsr_cholesky_invert(matrix, n, para_env, blacs_env, uplo_to_full)
used to replace the cholesky decomposition by the inverse
Definition cp_dbcsr_cholesky.F:114
cp_dbcsr_contrib
Definition cp_dbcsr_contrib.F:8
cp_dbcsr_contrib::dbcsr_set_diag
subroutine, public dbcsr_set_diag(matrix, diag)
Copies the diagonal elements from the given array into the given matrix.
Definition cp_dbcsr_contrib.F:585
cp_dbcsr_contrib::dbcsr_get_diag
subroutine, public dbcsr_get_diag(matrix, diag)
Copies the diagonal elements from the given matrix into the given array.
Definition cp_dbcsr_contrib.F:544
cp_dbcsr_contrib::dbcsr_add_on_diag
subroutine, public dbcsr_add_on_diag(matrix, alpha)
Adds the given scalar to the diagonal of the matrix. Reserves any missing diagonal blocks.
Definition cp_dbcsr_contrib.F:331
cp_dbcsr_contrib::dbcsr_maxabs
real(dp) function, public dbcsr_maxabs(matrix)
Compute the maxabs norm of a dbcsr matrix.
Definition cp_dbcsr_contrib.F:94
cp_dbcsr_contrib::dbcsr_frobenius_norm
real(dp) function, public dbcsr_frobenius_norm(matrix)
Compute the frobenius norm of a dbcsr matrix.
Definition cp_dbcsr_contrib.F:127
cp_dbcsr_contrib::dbcsr_dot
subroutine, public dbcsr_dot(matrix_a, matrix_b, trace)
Computes the dot product of two matrices, also known as the trace of their matrix product.
Definition cp_dbcsr_contrib.F:367
cp_dbcsr_contrib::dbcsr_hadamard_product
subroutine, public dbcsr_hadamard_product(matrix_a, matrix_b, matrix_c)
Hadamard product: C = A . B (C needs to be different from A and B)
Definition cp_dbcsr_contrib.F:53
cp_dbcsr_contrib::dbcsr_reserve_diag_blocks
subroutine, public dbcsr_reserve_diag_blocks(matrix)
Reserves all diagonal blocks.
Definition cp_dbcsr_contrib.F:255
cp_dbcsr_diag
Interface to (sca)lapack for the Cholesky based procedures.
Definition cp_dbcsr_diag.F:17
cp_dbcsr_diag::cp_dbcsr_syevd
subroutine, public cp_dbcsr_syevd(matrix, eigenvectors, eigenvalues, para_env, blacs_env)
...
Definition cp_dbcsr_diag.F:67
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition cp_log_handling.F:41
cp_log_handling::cp_logger_get_default_unit_nr
recursive integer function, public cp_logger_get_default_unit_nr(logger, local, skip_not_ionode)
asks the default unit number of the given logger. try to use cp_logger_get_unit_nr
Definition cp_log_handling.F:567
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition cp_log_handling.F:234
ct_methods
Cayley transformation methods.
Definition ct_methods.F:14
ct_methods::analytic_line_search
subroutine, public analytic_line_search(a, b, c, d, minima, nmins)
Finds real roots of a cubic equation ‍ a*x**3 + b*x**2 + c*x + d = 0 and returns only those roots for...
Definition ct_methods.F:1375
ct_methods::diagonalize_diagonal_blocks
subroutine, public diagonalize_diagonal_blocks(matrix, c, e)
Diagonalizes diagonal blocks of a symmetric dbcsr matrix and returs its eigenvectors.
Definition ct_methods.F:1498
ct_methods::ct_step_execute
subroutine, public ct_step_execute(cts_env)
Performs Cayley transformation.
Definition ct_methods.F:63
ct_types
Types for all cayley transformation methods.
Definition ct_types.F:14
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:17
input_constants::cg_hestenes_stiefel
integer, parameter, public cg_hestenes_stiefel
Definition input_constants.F:1021
input_constants::cg_fletcher
integer, parameter, public cg_fletcher
Definition input_constants.F:1021
input_constants::cg_fletcher_reeves
integer, parameter, public cg_fletcher_reeves
Definition input_constants.F:1021
input_constants::tensor_up_down
integer, parameter, public tensor_up_down
Definition input_constants.F:1046
input_constants::tensor_orthogonal
integer, parameter, public tensor_orthogonal
Definition input_constants.F:1046
input_constants::cg_dai_yuan
integer, parameter, public cg_dai_yuan
Definition input_constants.F:1021
input_constants::cg_liu_storey
integer, parameter, public cg_liu_storey
Definition input_constants.F:1021
input_constants::cg_hager_zhang
integer, parameter, public cg_hager_zhang
Definition input_constants.F:1021
input_constants::cg_zero
integer, parameter, public cg_zero
Definition input_constants.F:1021
input_constants::cg_polak_ribiere
integer, parameter, public cg_polak_ribiere
Definition input_constants.F:1021
iterate_matrix
Routines useful for iterative matrix calculations.
Definition iterate_matrix.F:13
iterate_matrix::matrix_sqrt_newton_schulz
subroutine, public matrix_sqrt_newton_schulz(matrix_sqrt, matrix_sqrt_inv, matrix, threshold, order, eps_lanczos, max_iter_lanczos, symmetrize, converged, iounit)
compute the sqrt of a matrix via the sign function and the corresponding Newton-Schulz iterations the...
Definition iterate_matrix.F:1626
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
kinds::sp
integer, parameter, public sp
Definition kinds.F:33
machine
Machine interface based on Fortran 2003 and POSIX.
Definition machine.F:17
machine::m_walltime
real(kind=dp) function, public m_walltime()
returns time from a real-time clock, protected against rolling early/easily
Definition machine.F:147
cp_dbcsr_api::dbcsr_iterator_type
Definition cp_dbcsr_api.F:186
cp_dbcsr_api::dbcsr_type
Definition cp_dbcsr_api.F:174
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
ct_types::ct_step_env_type
Definition ct_types.F:37
orbital
Orbital angular momentum.
Definition grid_common.h:128