(git:b195825)
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 ! **************************************************************************************************
14 MODULE 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 
49 CONTAINS
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 
1750 END 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 w...
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
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
orbital
Orbital angular momentum.
Definition: grid_common.h:125