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