(git:e7e05ae)
almo_scf_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 Subroutines for ALMO SCF
10 !> \par History
11 !> 2011.06 created [Rustam Z Khaliullin]
12 !> 2018.09 smearing support [Ruben Staub]
13 !> \author Rustam Z Khaliullin
14 ! **************************************************************************************************
15 MODULE almo_scf_methods
16  USE almo_scf_types, ONLY: almo_scf_env_type,&
17  almo_scf_history_type
18  USE bibliography, ONLY: kolafa2004,&
19  kuhne2007,&
20  cite_reference
21  USE cp_blacs_env, ONLY: cp_blacs_env_type
22  USE cp_dbcsr_cholesky, ONLY: cp_dbcsr_cholesky_decompose,&
23  cp_dbcsr_cholesky_invert
24  USE cp_log_handling, ONLY: cp_get_default_logger,&
25  cp_logger_get_default_unit_nr,&
26  cp_logger_type
27  USE dbcsr_api, ONLY: &
28  dbcsr_add, dbcsr_add_on_diag, dbcsr_copy, dbcsr_create, dbcsr_distribution_get, &
29  dbcsr_distribution_type, dbcsr_filter, dbcsr_finalize, dbcsr_frobenius_norm, &
30  dbcsr_get_block_p, dbcsr_get_diag, dbcsr_get_info, dbcsr_get_stored_coordinates, &
31  dbcsr_init_random, dbcsr_iterator_blocks_left, dbcsr_iterator_next_block, &
32  dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, dbcsr_multiply, &
33  dbcsr_nblkcols_total, dbcsr_nblkrows_total, dbcsr_print, dbcsr_release, &
34  dbcsr_reserve_block2d, dbcsr_scale_by_vector, dbcsr_set, dbcsr_set_diag, dbcsr_transposed, &
35  dbcsr_type, dbcsr_type_no_symmetry, dbcsr_type_symmetric, dbcsr_work_create
36  USE domain_submatrix_methods, ONLY: &
37  add_submatrices, construct_dbcsr_from_submatrices, construct_submatrices, &
38  copy_submatrices, copy_submatrix_data, init_submatrices, multiply_submatrices, &
39  print_submatrices, release_submatrices
40  USE domain_submatrix_types, ONLY: domain_map_type,&
41  domain_submatrix_type,&
42  select_row,&
43  select_row_col
44  USE fermi_utils, ONLY: fermifixed
45 !! for smearing
46  USE input_constants, ONLY: almo_domain_layout_molecular,&
47  almo_mat_distr_atomic,&
48  almo_scf_diag,&
49  spd_inversion_dense_cholesky,&
50  spd_inversion_ls_hotelling,&
51  spd_inversion_ls_taylor
52  USE iterate_matrix, ONLY: invert_hotelling,&
53  invert_taylor,&
54  matrix_sqrt_newton_schulz
55  USE kinds, ONLY: dp
56  USE mathlib, ONLY: binomial
57  USE message_passing, ONLY: mp_comm_type,&
58  mp_para_env_type
59  USE util, ONLY: sort
60 #include "./base/base_uses.f90"
61 
62  IMPLICIT NONE
63 
64  PRIVATE
65 
66  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'almo_scf_methods'
67 
68  PUBLIC almo_scf_ks_to_ks_blk, almo_scf_p_blk_to_t_blk, &
69  almo_scf_t_to_proj, almo_scf_ks_blk_to_tv_blk, &
70  almo_scf_ks_xx_to_tv_xx, almo_scf_t_rescaling, &
71  apply_projector, get_overlap, &
72  generator_to_unitary, &
73  orthogonalize_mos, &
74  pseudo_invert_diagonal_blk, construct_test, &
75  construct_domain_preconditioner, &
76  apply_domain_operators, &
77  construct_domain_s_inv, &
78  construct_domain_s_sqrt, &
79  distribute_domains, &
80  almo_scf_ks_to_ks_xx, &
81  construct_domain_r_down, &
82  xalmo_initial_guess, &
83  fill_matrix_with_ones
84 
85 CONTAINS
86 
87 ! **************************************************************************************************
88 !> \brief Fill all matrix blocks with 1.0_dp
89 !> \param matrix ...
90 !> \par History
91 !> 2019.09 created [Rustam Z Khaliullin]
92 !> \author Rustam Z Khaliullin
93 ! **************************************************************************************************
94  SUBROUTINE fill_matrix_with_ones(matrix)
95 
96  TYPE(dbcsr_type), INTENT(INOUT) :: matrix
97 
98  INTEGER :: col, hold, iblock_col, iblock_row, &
99  mynode, nblkcols_tot, nblkrows_tot, row
100  LOGICAL :: tr
101  REAL(kind=dp), DIMENSION(:, :), POINTER :: p_new_block
102  TYPE(dbcsr_distribution_type) :: dist
103 
104  CALL dbcsr_get_info(matrix, distribution=dist)
105  CALL dbcsr_distribution_get(dist, mynode=mynode)
106  CALL dbcsr_work_create(matrix, work_mutable=.true.)
107  nblkrows_tot = dbcsr_nblkrows_total(matrix)
108  nblkcols_tot = dbcsr_nblkcols_total(matrix)
109  DO row = 1, nblkrows_tot
110  DO col = 1, nblkcols_tot
111  tr = .false.
112  iblock_row = row
113  iblock_col = col
114  CALL dbcsr_get_stored_coordinates(matrix, &
115  iblock_row, iblock_col, hold)
116  IF (hold .EQ. mynode) THEN
117  NULLIFY (p_new_block)
118  CALL dbcsr_reserve_block2d(matrix, &
119  iblock_row, iblock_col, p_new_block)
120  cpassert(ASSOCIATED(p_new_block))
121  p_new_block(:, :) = 1.0_dp
122  END IF
123  END DO
124  END DO
125  CALL dbcsr_finalize(matrix)
126 
127  END SUBROUTINE fill_matrix_with_ones
128 
129 ! **************************************************************************************************
130 !> \brief builds projected KS matrices for the overlapping domains
131 !> also computes the DIIS error vector as a by-product
132 !> \param almo_scf_env ...
133 !> \par History
134 !> 2013.03 created [Rustam Z Khaliullin]
135 !> \author Rustam Z Khaliullin
136 ! **************************************************************************************************
137  SUBROUTINE almo_scf_ks_to_ks_xx(almo_scf_env)
138 
139  TYPE(almo_scf_env_type), INTENT(INOUT) :: almo_scf_env
140 
141  CHARACTER(LEN=*), PARAMETER :: routinen = 'almo_scf_ks_to_ks_xx'
142 
143  INTEGER :: handle, ispin, ndomains
144  REAL(kind=dp) :: eps_multiply
145  TYPE(dbcsr_type) :: matrix_tmp1, matrix_tmp2, matrix_tmp3, matrix_tmp4, matrix_tmp5, &
146  matrix_tmp6, matrix_tmp7, matrix_tmp8, matrix_tmp9
147  TYPE(domain_submatrix_type), ALLOCATABLE, &
148  DIMENSION(:) :: subm_tmp1, subm_tmp2, subm_tmp3
149 
150  CALL timeset(routinen, handle)
151 
152  eps_multiply = almo_scf_env%eps_filter
153 
154  DO ispin = 1, almo_scf_env%nspins
155 
156  ndomains = dbcsr_nblkcols_total(almo_scf_env%quench_t(ispin))
157 
158  ! 0. Create KS_xx
159  CALL construct_submatrices( &
160  almo_scf_env%matrix_ks(ispin), &
161  almo_scf_env%domain_ks_xx(:, ispin), &
162  almo_scf_env%quench_t(ispin), &
163  almo_scf_env%domain_map(ispin), &
164  almo_scf_env%cpu_of_domain, &
165  select_row_col)
166 
167  !!!!! RZK-warning MAKE SURE THAT YOU NEED BLOCKS OUTSIDE QUENCH_T
168  !!!!! FOR ALL NO-MATRICES NOT COMPUTING THEM CAN SAVE LOTS OF TIME
169 
170  ! 1. TMP1=KS.T
171  ! Cost: NOn
172  !matrix_tmp1 = create NxO, full
173  CALL dbcsr_create(matrix_tmp1, &
174  template=almo_scf_env%matrix_t(ispin))
175  CALL dbcsr_multiply("N", "N", 1.0_dp, almo_scf_env%matrix_ks(ispin), &
176  almo_scf_env%matrix_t(ispin), &
177  0.0_dp, matrix_tmp1, &
178  filter_eps=eps_multiply)
179 
180  ! 2. TMP2=TMP1.SigInv=KS.T.SigInv
181  ! Cost: NOO
182  !matrix_tmp2 = create NxO, full
183  CALL dbcsr_create(matrix_tmp2, &
184  template=almo_scf_env%matrix_t(ispin))
185  CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_tmp1, &
186  almo_scf_env%matrix_sigma_inv(ispin), &
187  0.0_dp, matrix_tmp2, &
188  filter_eps=eps_multiply)
189 
190  ! 3. TMP1=S.T
191  ! Cost: NOn
192  CALL dbcsr_multiply("N", "N", 1.0_dp, almo_scf_env%matrix_s(1), &
193  almo_scf_env%matrix_t(ispin), &
194  0.0_dp, matrix_tmp1, &
195  filter_eps=eps_multiply)
196 
197  ! 4. TMP4=TMP2.tr(TMP1)=KS.T.SigInv.tr(T).S
198  ! Cost: NNO
199  !matrix_tmp4 = create NxN
200  CALL dbcsr_create(matrix_tmp4, &
201  template=almo_scf_env%matrix_s(1), &
202  matrix_type=dbcsr_type_no_symmetry)
203  CALL dbcsr_multiply("N", "T", 1.0_dp, matrix_tmp2, &
204  matrix_tmp1, &
205  0.0_dp, matrix_tmp4, &
206  filter_eps=eps_multiply)
207 
208  ! 5. KS_xx=KS_xx-TMP4_xx-tr(TMP4_xx)
209  ALLOCATE (subm_tmp1(ndomains))
210  CALL init_submatrices(subm_tmp1)
211  CALL construct_submatrices( &
212  matrix_tmp4, &
213  subm_tmp1, &
214  almo_scf_env%quench_t(ispin), &
215  almo_scf_env%domain_map(ispin), &
216  almo_scf_env%cpu_of_domain, &
217  select_row_col)
218  CALL add_submatrices(1.0_dp, almo_scf_env%domain_ks_xx(:, ispin), &
219  -1.0_dp, subm_tmp1, 'N')
220  CALL add_submatrices(1.0_dp, almo_scf_env%domain_ks_xx(:, ispin), &
221  -1.0_dp, subm_tmp1, 'T')
222 
223  ! 6. TMP3=tr(TMP4).T=S.T.SigInv.tr(T).KS.T
224  ! Cost: NOn
225  !matrix_tmp3 = create NxO, full
226  CALL dbcsr_create(matrix_tmp3, &
227  template=almo_scf_env%matrix_t(ispin), &
228  matrix_type=dbcsr_type_no_symmetry)
229  CALL dbcsr_multiply("T", "N", 1.0_dp, &
230  matrix_tmp4, &
231  almo_scf_env%matrix_t(ispin), &
232  0.0_dp, matrix_tmp3, &
233  filter_eps=eps_multiply)
234  CALL dbcsr_release(matrix_tmp4)
235 
236  ! 8. TMP6=TMP3.SigInv=S.T.SigInv.tr(T).KS.T.SigInv
237  ! Cost: NOO
238  !matrix_tmp6 = create NxO, full
239  CALL dbcsr_create(matrix_tmp6, &
240  template=almo_scf_env%matrix_t(ispin), &
241  matrix_type=dbcsr_type_no_symmetry)
242  CALL dbcsr_multiply("N", "N", 1.0_dp, &
243  matrix_tmp3, &
244  almo_scf_env%matrix_sigma_inv(ispin), &
245  0.0_dp, matrix_tmp6, &
246  filter_eps=eps_multiply)
247 
248  ! 8A. Use intermediate matrices to evaluate the gradient/error
249  ! Err=(TMP2-TMP6)_q=(KS.T.SigInv-S.T.SigInv.tr(T).KS.T.SigInv)_q
250  ! error vector in AO-MO basis
251  CALL dbcsr_copy(almo_scf_env%matrix_err_xx(ispin), &
252  almo_scf_env%quench_t(ispin))
253  CALL dbcsr_copy(almo_scf_env%matrix_err_xx(ispin), &
254  matrix_tmp2, keep_sparsity=.true.)
255  CALL dbcsr_create(matrix_tmp4, &
256  template=almo_scf_env%matrix_t(ispin), &
257  matrix_type=dbcsr_type_no_symmetry)
258  CALL dbcsr_copy(matrix_tmp4, &
259  almo_scf_env%quench_t(ispin))
260  CALL dbcsr_copy(matrix_tmp4, &
261  matrix_tmp6, keep_sparsity=.true.)
262  CALL dbcsr_add(almo_scf_env%matrix_err_xx(ispin), &
263  matrix_tmp4, 1.0_dp, -1.0_dp)
264  CALL dbcsr_release(matrix_tmp4)
265  !
266  ! error vector in AO-AO basis
267  ! RZK-warning tmp4 can be created using the sparsity pattern,
268  ! then retain_sparsity can be used to perform the multiply
269  ! this will save some time
270  CALL dbcsr_copy(matrix_tmp3, &
271  matrix_tmp2)
272  CALL dbcsr_add(matrix_tmp3, &
273  matrix_tmp6, 1.0_dp, -1.0_dp)
274  CALL dbcsr_create(matrix_tmp4, &
275  template=almo_scf_env%matrix_s(1), &
276  matrix_type=dbcsr_type_no_symmetry)
277  CALL dbcsr_multiply("N", "T", 1.0_dp, &
278  matrix_tmp3, &
279  almo_scf_env%matrix_t(ispin), &
280  0.0_dp, matrix_tmp4, &
281  filter_eps=eps_multiply)
282  CALL construct_submatrices( &
283  matrix_tmp4, &
284  almo_scf_env%domain_err(:, ispin), &
285  almo_scf_env%quench_t(ispin), &
286  almo_scf_env%domain_map(ispin), &
287  almo_scf_env%cpu_of_domain, &
288  select_row_col)
289  CALL dbcsr_release(matrix_tmp4)
290  ! domain_err submatrices are in down-up representation
291  ! bring them into the orthogonalized basis
292  ALLOCATE (subm_tmp2(ndomains))
293  CALL init_submatrices(subm_tmp2)
294  CALL multiply_submatrices('N', 'N', 1.0_dp, &
295  almo_scf_env%domain_err(:, ispin), &
296  almo_scf_env%domain_s_sqrt(:, ispin), 0.0_dp, subm_tmp2)
297  CALL multiply_submatrices('N', 'N', 1.0_dp, &
298  almo_scf_env%domain_s_sqrt_inv(:, ispin), &
299  subm_tmp2, 0.0_dp, almo_scf_env%domain_err(:, ispin))
300 
301  ! 9. TMP5=TMP6.tr(TMP1)=S.T.SigInv.tr(T).KS.T.SigInv.tr(T).S
302  ! Cost: NNO
303  ! matrix_tmp5 = create NxN, full
304  ! RZK-warning tmp5 can be created using the sparsity pattern,
305  ! then retain_sparsity can be used to perform the multiply
306  ! this will save some time
307  CALL dbcsr_create(matrix_tmp5, &
308  template=almo_scf_env%matrix_s(1), &
309  matrix_type=dbcsr_type_no_symmetry)
310  CALL dbcsr_multiply("N", "T", 1.0_dp, &
311  matrix_tmp6, &
312  matrix_tmp1, &
313  0.0_dp, matrix_tmp5, &
314  filter_eps=eps_multiply)
315 
316  ! 10. KS_xx=KS_xx+TMP5_xx
317  CALL construct_submatrices( &
318  matrix_tmp5, &
319  subm_tmp1, &
320  almo_scf_env%quench_t(ispin), &
321  almo_scf_env%domain_map(ispin), &
322  almo_scf_env%cpu_of_domain, &
323  select_row_col)
324  CALL dbcsr_release(matrix_tmp5)
325  CALL add_submatrices(1.0_dp, almo_scf_env%domain_ks_xx(:, ispin), &
326  1.0_dp, subm_tmp1, 'N')
327 
328  ! 11. KS_xx=KS_xx + [S.T]_xx.[SigInv.tr(T).KS.(1-T.SigInv.tr(T).S)]_xx + transposed
329  ALLOCATE (subm_tmp3(ndomains))
330  CALL init_submatrices(subm_tmp3)
331  CALL construct_submatrices( &
332  matrix_tmp2, &
333  subm_tmp2, &
334  almo_scf_env%quench_t(ispin), &
335  almo_scf_env%domain_map(ispin), &
336  almo_scf_env%cpu_of_domain, &
337  select_row)
338  CALL construct_submatrices( &
339  matrix_tmp6, &
340  subm_tmp3, &
341  almo_scf_env%quench_t(ispin), &
342  almo_scf_env%domain_map(ispin), &
343  almo_scf_env%cpu_of_domain, &
344  select_row)
345  CALL dbcsr_release(matrix_tmp6)
346  CALL add_submatrices(1.0_dp, subm_tmp2, &
347  -1.0_dp, subm_tmp3, 'N')
348  CALL construct_submatrices( &
349  matrix_tmp1, &
350  subm_tmp3, &
351  almo_scf_env%quench_t(ispin), &
352  almo_scf_env%domain_map(ispin), &
353  almo_scf_env%cpu_of_domain, &
354  select_row)
355  CALL multiply_submatrices('N', 'T', 1.0_dp, subm_tmp2, &
356  subm_tmp3, 0.0_dp, subm_tmp1)
357  CALL add_submatrices(1.0_dp, almo_scf_env%domain_ks_xx(:, ispin), &
358  1.0_dp, subm_tmp1, 'N')
359  CALL add_submatrices(1.0_dp, almo_scf_env%domain_ks_xx(:, ispin), &
360  1.0_dp, subm_tmp1, 'T')
361 
362  ! 12. TMP7=tr(T).KS.T.SigInv
363  CALL dbcsr_create(matrix_tmp7, &
364  template=almo_scf_env%matrix_sigma_blk(ispin), &
365  matrix_type=dbcsr_type_no_symmetry)
366  CALL dbcsr_multiply("T", "N", 1.0_dp, &
367  almo_scf_env%matrix_t(ispin), &
368  matrix_tmp2, &
369  0.0_dp, matrix_tmp7, &
370  filter_eps=eps_multiply)
371 
372  ! 13. TMP8=[SigInv.tr(T).KS.T.SigInv]_xx
373  CALL dbcsr_create(matrix_tmp8, &
374  template=almo_scf_env%matrix_sigma_blk(ispin), &
375  matrix_type=dbcsr_type_symmetric)
376  CALL dbcsr_copy(matrix_tmp8, almo_scf_env%matrix_sigma_blk(ispin))
377  CALL dbcsr_multiply("N", "N", 1.0_dp, &
378  almo_scf_env%matrix_sigma_inv(ispin), &
379  matrix_tmp7, &
380  0.0_dp, matrix_tmp8, &
381  retain_sparsity=.true., &
382  filter_eps=eps_multiply)
383  CALL dbcsr_release(matrix_tmp7)
384 
385  ! 13. TMP9=[S.T]_xx
386  CALL dbcsr_create(matrix_tmp9, &
387  template=almo_scf_env%matrix_t(ispin), &
388  matrix_type=dbcsr_type_no_symmetry)
389  CALL dbcsr_copy(matrix_tmp9, almo_scf_env%quench_t(ispin))
390  CALL dbcsr_copy(matrix_tmp9, matrix_tmp1, keep_sparsity=.true.)
391 
392  ! 14. TMP3=TMP9.TMP8=[S.T]_xx.[SigInv.tr(T).KS.T.SigInv]_xx
393  CALL dbcsr_multiply("N", "N", 1.0_dp, &
394  matrix_tmp9, &
395  matrix_tmp8, &
396  0.0_dp, matrix_tmp3, &
397  filter_eps=eps_multiply)
398  CALL dbcsr_release(matrix_tmp8)
399  CALL dbcsr_release(matrix_tmp9)
400 
401  ! 15. KS_xx=KS_xx+[S.T]_xx.[SigInv.tr(T).KS.T.SigInv]_xx.[tr(T).S]_xx
402  CALL construct_submatrices( &
403  matrix_tmp3, &
404  subm_tmp2, &
405  almo_scf_env%quench_t(ispin), &
406  almo_scf_env%domain_map(ispin), &
407  almo_scf_env%cpu_of_domain, &
408  select_row)
409  CALL multiply_submatrices('N', 'T', 1.0_dp, subm_tmp2, &
410  subm_tmp3, 0.0_dp, subm_tmp1)
411  CALL add_submatrices(1.0_dp, almo_scf_env%domain_ks_xx(:, ispin), &
412  1.0_dp, subm_tmp1, 'N')
413 
414  !!!!!!! use intermediate matrices to get the error vector !!!!!!!
415  !!!!!!! make sure s_blk_sqrt and its inverse exist (i.e. we use diag algorithm)
416  !CPPrecondition(almo_scf_env%almo_update_algorithm.eq.almo_scf_diag,cp_failure_level,routineP,failure)
417  !! tmp_err = (1-S.T_blk.SigInv.tr(T_blk)).F.T_blk.SigInv
418  !CALL dbcsr_init(matrix_tmp_err)
419  !CALL dbcsr_create(matrix_tmp_err,&
420  ! template=almo_scf_env%matrix_t(ispin))
421  !CALL dbcsr_copy(matrix_tmp_err,&
422  ! matrix_tmp2)
423  !CALL dbcsr_add(matrix_tmp_err,matrix_tmp3,&
424  ! 1.0_dp,-1.0_dp)
425  !! err_blk = tmp_err.tr(T_blk)
426  !CALL dbcsr_copy(almo_scf_env%matrix_err_blk(ispin),&
427  ! almo_scf_env%matrix_s_blk_sqrt(1))
428  !CALL dbcsr_multiply("N", "T", 1.0_dp, matrix_tmp_err,&
429  ! almo_scf_env%matrix_t(ispin),&
430  ! 0.0_dp, almo_scf_env%matrix_err_blk(ispin),&
431  ! retain_sparsity=.TRUE.,&
432  ! filter_eps=eps_multiply)
433  !CALL dbcsr_release(matrix_tmp_err)
434  !! bring to the orthogonal basis
435  !! err_blk = (S_blk^-1/2).err_blk.(S_blk^1/2)
436  !CALL dbcsr_init(matrix_tmp_err)
437  !CALL dbcsr_create(matrix_tmp_err,&
438  ! template=almo_scf_env%matrix_err_blk(ispin))
439  !CALL dbcsr_multiply("N", "N", 1.0_dp,&
440  ! almo_scf_env%matrix_err_blk(ispin),&
441  ! almo_scf_env%matrix_s_blk_sqrt(1),&
442  ! 0.0_dp, matrix_tmp_err,&
443  ! filter_eps=eps_multiply)
444  !CALL dbcsr_multiply("N", "N", 1.0_dp,&
445  ! almo_scf_env%matrix_s_blk_sqrt_inv(1),&
446  ! matrix_tmp_err,&
447  ! 0.0_dp, almo_scf_env%matrix_err_blk(ispin),&
448  ! filter_eps=eps_multiply)
449  !! subtract transpose
450  !CALL dbcsr_transposed(matrix_tmp_err,&
451  ! almo_scf_env%matrix_err_blk(ispin))
452  !CALL dbcsr_add(almo_scf_env%matrix_err_blk(ispin),&
453  ! matrix_tmp_err,&
454  ! 1.0_dp,-1.0_dp)
455  !CALL dbcsr_release(matrix_tmp_err)
456  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
457 
458  CALL release_submatrices(subm_tmp3)
459  CALL release_submatrices(subm_tmp2)
460  CALL release_submatrices(subm_tmp1)
461  DEALLOCATE (subm_tmp3)
462  DEALLOCATE (subm_tmp2)
463  DEALLOCATE (subm_tmp1)
464  CALL dbcsr_release(matrix_tmp3)
465  CALL dbcsr_release(matrix_tmp2)
466  CALL dbcsr_release(matrix_tmp1)
467 
468  END DO ! spins
469 
470  CALL timestop(handle)
471 
472  END SUBROUTINE almo_scf_ks_to_ks_xx
473 
474 ! **************************************************************************************************
475 !> \brief computes the projected KS from the total KS matrix
476 !> also computes the DIIS error vector as a by-product
477 !> \param almo_scf_env ...
478 !> \par History
479 !> 2011.06 created [Rustam Z Khaliullin]
480 !> \author Rustam Z Khaliullin
481 ! **************************************************************************************************
482  SUBROUTINE almo_scf_ks_to_ks_blk(almo_scf_env)
483 
484  TYPE(almo_scf_env_type), INTENT(INOUT) :: almo_scf_env
485 
486  CHARACTER(LEN=*), PARAMETER :: routinen = 'almo_scf_ks_to_ks_blk'
487 
488  INTEGER :: handle, ispin
489  REAL(kind=dp) :: eps_multiply
490  TYPE(dbcsr_type) :: matrix_tmp1, matrix_tmp2, matrix_tmp3, matrix_tmp4, matrix_tmp5, &
491  matrix_tmp6, matrix_tmp7, matrix_tmp8, matrix_tmp9, matrix_tmp_err
492 
493  CALL timeset(routinen, handle)
494 
495  eps_multiply = almo_scf_env%eps_filter
496 
497  DO ispin = 1, almo_scf_env%nspins
498 
499  ! 1. TMP1=KS.T_blk
500  ! Cost: NOn
501  !matrix_tmp1 = create NxO, full
502  CALL dbcsr_create(matrix_tmp1, &
503  template=almo_scf_env%matrix_t(ispin))
504  CALL dbcsr_multiply("N", "N", 1.0_dp, almo_scf_env%matrix_ks(ispin), &
505  almo_scf_env%matrix_t_blk(ispin), &
506  0.0_dp, matrix_tmp1, &
507  filter_eps=eps_multiply)
508  ! 2. TMP2=TMP1.SigInv=KS.T_blk.SigInv
509  ! Cost: NOO
510  !matrix_tmp2 = create NxO, full
511  CALL dbcsr_create(matrix_tmp2, &
512  template=almo_scf_env%matrix_t(ispin))
513  CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_tmp1, &
514  almo_scf_env%matrix_sigma_inv(ispin), &
515  0.0_dp, matrix_tmp2, &
516  filter_eps=eps_multiply)
517 
518  !!!!!! use intermediate matrices to get the error vector !!!!!!!
519  !CALL dbcsr_copy(almo_scf_env%matrix_err_blk(ispin),&
520  ! almo_scf_env%matrix_t_blk(ispin))
521  !CALL dbcsr_copy(almo_scf_env%matrix_err_blk(ispin),&
522  ! matrix_tmp2,&
523  ! keep_sparsity=.TRUE.)
524  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
525 
526  ! 3. TMP1=S.T_blk
527  ! Cost: NOn
528  CALL dbcsr_multiply("N", "N", 1.0_dp, almo_scf_env%matrix_s(1), &
529  almo_scf_env%matrix_t_blk(ispin), &
530  0.0_dp, matrix_tmp1, &
531  filter_eps=eps_multiply)
532 
533  ! 4. TMP4_blk=TMP2.tr(TMP1)=KS.T_blk.SigInv.tr(T_blk).S
534  ! Cost: NnO
535  !matrix_tmp4 = create NxN, blk
536  CALL dbcsr_create(matrix_tmp4, &
537  template=almo_scf_env%matrix_s_blk(1), &
538  matrix_type=dbcsr_type_no_symmetry)
539  CALL dbcsr_copy(matrix_tmp4, almo_scf_env%matrix_s_blk(1))
540  CALL dbcsr_multiply("N", "T", 1.0_dp, matrix_tmp2, &
541  matrix_tmp1, &
542  0.0_dp, matrix_tmp4, &
543  retain_sparsity=.true., &
544  filter_eps=eps_multiply)
545 
546  ! 5. KS_blk=KS_blk-TMP4_blk
547  CALL dbcsr_copy(almo_scf_env%matrix_ks_blk(ispin), &
548  almo_scf_env%matrix_ks(ispin), keep_sparsity=.true.)
549  CALL dbcsr_add(almo_scf_env%matrix_ks_blk(ispin), &
550  matrix_tmp4, &
551  1.0_dp, -1.0_dp)
552 
553  ! 6. TMP5_blk=tr(TMP4_blk)
554  ! KS_blk=KS_blk-tr(TMP4_blk)
555  !matrix_tmp5 = create NxN, blk
556  CALL dbcsr_create(matrix_tmp5, &
557  template=almo_scf_env%matrix_s_blk(1), &
558  matrix_type=dbcsr_type_no_symmetry)
559  CALL dbcsr_transposed(matrix_tmp5, matrix_tmp4)
560  CALL dbcsr_add(almo_scf_env%matrix_ks_blk(ispin), matrix_tmp5, &
561  1.0_dp, -1.0_dp)
562 
563  ! 7. TMP3=tr(T_blk).TMP2=tr(T_blk).KS.T_blk.SigInv
564  ! Cost: OOn
565  !matrix_tmp3 = create OxO, full
566  CALL dbcsr_create(matrix_tmp3, &
567  template=almo_scf_env%matrix_sigma_inv(ispin), &
568  matrix_type=dbcsr_type_no_symmetry)
569  CALL dbcsr_multiply("T", "N", 1.0_dp, &
570  almo_scf_env%matrix_t_blk(ispin), &
571  matrix_tmp2, &
572  0.0_dp, matrix_tmp3, &
573  filter_eps=eps_multiply)
574 
575  ! 8. TMP6=SigInv.TMP3=SigInv.tr(T_blk).KS.T_blk.SigInv
576  ! Cost: OOO
577  !matrix_tmp6 = create OxO, full
578  CALL dbcsr_create(matrix_tmp6, &
579  template=almo_scf_env%matrix_sigma_inv(ispin), &
580  matrix_type=dbcsr_type_no_symmetry)
581  CALL dbcsr_multiply("N", "N", 1.0_dp, &
582  almo_scf_env%matrix_sigma_inv(ispin), &
583  matrix_tmp3, &
584  0.0_dp, matrix_tmp6, &
585  filter_eps=eps_multiply)
586 
587  ! 9. TMP3=TMP1.TMP6=S.T_blk.SigInv.tr(T_blk).KS.T_blk.SigInv
588  ! Cost: NOO
589  !matrix_tmp3 = re-create NxO, full
590  CALL dbcsr_release(matrix_tmp3)
591  CALL dbcsr_create(matrix_tmp3, &
592  template=almo_scf_env%matrix_t(ispin))
593  CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_tmp1, &
594  matrix_tmp6, &
595  0.0_dp, matrix_tmp3, &
596  filter_eps=eps_multiply)
597 
598  !!!!!! use intermediate matrices to get the error vector !!!!!!!
599  !CALL dbcsr_init(matrix_tmp_err)
600  !CALL dbcsr_create(matrix_tmp_err,&
601  ! template=almo_scf_env%matrix_t_blk(ispin))
602  !CALL dbcsr_copy(matrix_tmp_err,&
603  ! almo_scf_env%matrix_t_blk(ispin))
604  !CALL dbcsr_copy(matrix_tmp_err,matrix_tmp3,&
605  ! keep_sparsity=.TRUE.)
606  !CALL dbcsr_add(almo_scf_env%matrix_err_blk(ispin),matrix_tmp_err,&
607  ! 1.0_dp,-1.0_dp)
608  !CALL dbcsr_release(matrix_tmp_err)
609  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
610 
611  !!!!!! use intermediate matrices to get the error vector !!!!!!!
612  !!!!!! make sure s_blk_sqrt and its inverse exist (i.e. we use diag algorithm)
613  cpassert(almo_scf_env%almo_update_algorithm .EQ. almo_scf_diag)
614  ! tmp_err = (1-S.T_blk.SigInv.tr(T_blk)).F.T_blk.SigInv
615  CALL dbcsr_create(matrix_tmp_err, &
616  template=almo_scf_env%matrix_t_blk(ispin))
617  CALL dbcsr_copy(matrix_tmp_err, &
618  matrix_tmp2)
619  CALL dbcsr_add(matrix_tmp_err, matrix_tmp3, &
620  1.0_dp, -1.0_dp)
621  ! err_blk = tmp_err.tr(T_blk)
622  CALL dbcsr_copy(almo_scf_env%matrix_err_blk(ispin), &
623  almo_scf_env%matrix_s_blk_sqrt(1))
624  CALL dbcsr_multiply("N", "T", 1.0_dp, matrix_tmp_err, &
625  almo_scf_env%matrix_t_blk(ispin), &
626  0.0_dp, almo_scf_env%matrix_err_blk(ispin), &
627  retain_sparsity=.true., &
628  filter_eps=eps_multiply)
629  CALL dbcsr_release(matrix_tmp_err)
630  ! bring to the orthogonal basis
631  ! err_blk = (S_blk^-1/2).err_blk.(S_blk^1/2)
632  CALL dbcsr_create(matrix_tmp_err, &
633  template=almo_scf_env%matrix_err_blk(ispin))
634  CALL dbcsr_multiply("N", "N", 1.0_dp, &
635  almo_scf_env%matrix_err_blk(ispin), &
636  almo_scf_env%matrix_s_blk_sqrt(1), &
637  0.0_dp, matrix_tmp_err, &
638  filter_eps=eps_multiply)
639  CALL dbcsr_multiply("N", "N", 1.0_dp, &
640  almo_scf_env%matrix_s_blk_sqrt_inv(1), &
641  matrix_tmp_err, &
642  0.0_dp, almo_scf_env%matrix_err_blk(ispin), &
643  filter_eps=eps_multiply)
644 
645  ! subtract transpose
646  CALL dbcsr_transposed(matrix_tmp_err, &
647  almo_scf_env%matrix_err_blk(ispin))
648  CALL dbcsr_add(almo_scf_env%matrix_err_blk(ispin), &
649  matrix_tmp_err, &
650  1.0_dp, -1.0_dp)
651  CALL dbcsr_release(matrix_tmp_err)
652  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
653 
654  ! later we will need only the blk version of TMP6
655  ! create it here and release TMP6
656  !matrix_tmp9 = create OxO, blk
657  !matrix_tmp9 = copy data from matrix_tmp6, retain sparsity
658  !matrix_tmp6 = release
659  CALL dbcsr_create(matrix_tmp9, &
660  template=almo_scf_env%matrix_sigma_blk(ispin), &
661  matrix_type=dbcsr_type_no_symmetry)
662  CALL dbcsr_copy(matrix_tmp9, almo_scf_env%matrix_sigma_blk(ispin))
663  CALL dbcsr_copy(matrix_tmp9, matrix_tmp6, keep_sparsity=.true.)
664  CALL dbcsr_release(matrix_tmp6)
665 
666  !10. KS_blk=KS_blk+TMP3.tr(TMP1)=
667  ! =KS_blk+S.T_blk.SigInv.tr.(T_blk).KS.T_blk.SigInv.tr(T_blk).S
668  ! Cost: NnO
669  CALL dbcsr_multiply("N", "T", 1.0_dp, matrix_tmp3, &
670  matrix_tmp1, &
671  1.0_dp, almo_scf_env%matrix_ks_blk(ispin), &
672  retain_sparsity=.true., &
673  filter_eps=eps_multiply)
674 
675  ! 11. TMP4_blk=TMP7_blk.tr(TMP8_blk)
676  ! Cost: Nnn
677  !matrix_tmp7 = create NxO, blk
678  !matrix_tmp7 = copy data from matrix_tmp3, retain sparsity
679  !matrix_tmp3 = release
680  !matrix_tmp8 = create NxO, blk
681  !matrix_tmp8 = copy data from matrix_tmp1, retain sparsity
682  !matrix_tmp1 = release
683  CALL dbcsr_create(matrix_tmp7, &
684  template=almo_scf_env%matrix_t_blk(ispin))
685  ! transfer only the ALMO blocks from tmp3 into tmp7:
686  ! first, copy t_blk into tmp7 to transfer the blk structure,
687  ! then copy tmp3 into tmp7 with retain_sparsity
688  CALL dbcsr_copy(matrix_tmp7, almo_scf_env%matrix_t_blk(ispin))
689  CALL dbcsr_copy(matrix_tmp7, matrix_tmp3, keep_sparsity=.true.)
690  CALL dbcsr_release(matrix_tmp3)
691  ! do the same for tmp1->tmp8
692  CALL dbcsr_create(matrix_tmp8, &
693  template=almo_scf_env%matrix_t_blk(ispin))
694  CALL dbcsr_copy(matrix_tmp8, almo_scf_env%matrix_t_blk(ispin))
695  CALL dbcsr_copy(matrix_tmp8, matrix_tmp1, keep_sparsity=.true.)
696  CALL dbcsr_release(matrix_tmp1)
697  CALL dbcsr_multiply("N", "T", 1.0_dp, matrix_tmp7, &
698  matrix_tmp8, &
699  0.0_dp, matrix_tmp4, &
700  filter_eps=eps_multiply, &
701  retain_sparsity=.true.)
702 
703  ! 12. KS_blk=KS_blk-TMP4_blk
704  CALL dbcsr_add(almo_scf_env%matrix_ks_blk(ispin), matrix_tmp4, &
705  1.0_dp, -1.0_dp)
706 
707  ! 13. TMP5_blk=tr(TMP5_blk)
708  ! KS_blk=KS_blk-tr(TMP4_blk)
709  CALL dbcsr_transposed(matrix_tmp5, matrix_tmp4)
710  CALL dbcsr_add(almo_scf_env%matrix_ks_blk(ispin), matrix_tmp5, &
711  1.0_dp, -1.0_dp)
712 
713  ! 14. TMP4_blk=TMP7_blk.tr(TMP8_blk)
714  ! Cost: Nnn
715  CALL dbcsr_copy(matrix_tmp7, matrix_tmp2, keep_sparsity=.true.)
716  CALL dbcsr_release(matrix_tmp2)
717  CALL dbcsr_multiply("N", "T", 1.0_dp, matrix_tmp7, &
718  matrix_tmp8, &
719  0.0_dp, matrix_tmp4, &
720  retain_sparsity=.true., &
721  filter_eps=eps_multiply)
722  ! 15. KS_blk=KS_blk+TMP4_blk
723  CALL dbcsr_add(almo_scf_env%matrix_ks_blk(ispin), matrix_tmp4, &
724  1.0_dp, 1.0_dp)
725 
726  ! 16. KS_blk=KS_blk+tr(TMP4_blk)
727  CALL dbcsr_transposed(matrix_tmp5, matrix_tmp4)
728  CALL dbcsr_release(matrix_tmp4)
729  CALL dbcsr_add(almo_scf_env%matrix_ks_blk(ispin), matrix_tmp5, &
730  1.0_dp, 1.0_dp)
731  CALL dbcsr_release(matrix_tmp5)
732 
733  ! 17. TMP10_blk=TMP8_blk.TMP9_blk
734  ! Cost: Noo
735  CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_tmp8, &
736  matrix_tmp9, &
737  0.0_dp, matrix_tmp7, &
738  retain_sparsity=.true., &
739  filter_eps=eps_multiply)
740  CALL dbcsr_release(matrix_tmp9)
741 
742  ! 18. KS_blk=TMP7_blk.tr(TMP8_blk)
743  ! Cost: Nno
744  CALL dbcsr_multiply("N", "T", 1.0_dp, matrix_tmp7, &
745  matrix_tmp8, &
746  1.0_dp, almo_scf_env%matrix_ks_blk(ispin), &
747  retain_sparsity=.true., &
748  filter_eps=eps_multiply)
749  CALL dbcsr_release(matrix_tmp7)
750  CALL dbcsr_release(matrix_tmp8)
751 
752  END DO ! spins
753 
754  CALL timestop(handle)
755 
756  END SUBROUTINE almo_scf_ks_to_ks_blk
757 
758 ! **************************************************************************************************
759 !> \brief ALMOs by diagonalizing the KS domain submatrices
760 !> computes both the occupied and virtual orbitals
761 !> \param almo_scf_env ...
762 !> \par History
763 !> 2013.03 created [Rustam Z Khaliullin]
764 !> 2018.09 smearing support [Ruben Staub]
765 !> \author Rustam Z Khaliullin
766 ! **************************************************************************************************
767  SUBROUTINE almo_scf_ks_xx_to_tv_xx(almo_scf_env)
768 
769  TYPE(almo_scf_env_type), INTENT(INOUT) :: almo_scf_env
770 
771  CHARACTER(LEN=*), PARAMETER :: routinen = 'almo_scf_ks_xx_to_tv_xx'
772 
773  INTEGER :: handle, iblock_size, idomain, info, &
774  ispin, lwork, ndomains
775  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: eigenvalues, work
776  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: data_copy
777  TYPE(domain_submatrix_type), ALLOCATABLE, &
778  DIMENSION(:) :: subm_ks_xx_orthog, subm_t, subm_tmp
779 
780  CALL timeset(routinen, handle)
781 
782  IF (almo_scf_env%domain_layout_aos == almo_domain_layout_molecular .AND. &
783  almo_scf_env%mat_distr_aos == almo_mat_distr_atomic) THEN
784  cpabort("a domain must be located entirely on a CPU")
785  END IF
786 
787  ndomains = almo_scf_env%ndomains
788  ALLOCATE (subm_tmp(ndomains))
789  ALLOCATE (subm_ks_xx_orthog(ndomains))
790  ALLOCATE (subm_t(ndomains))
791 
792  DO ispin = 1, almo_scf_env%nspins
793 
794  CALL init_submatrices(subm_tmp)
795  CALL init_submatrices(subm_ks_xx_orthog)
796 
797  ! TRY: project out T0-occupied space for each domain
798  ! F=(1-R_du).F.(1-tr(R_du))
799  !CALL copy_submatrices(almo_scf_env%domain_ks_xx(:,ispin),&
800  ! subm_ks_xx_orthog,copy_data=.TRUE.)
801  !CALL multiply_submatrices('N','N',1.0_dp,&
802  ! almo_scf_env%domain_r_down_up(:,ispin),&
803  ! almo_scf_env%domain_ks_xx(:,ispin),0.0_dp,subm_tmp)
804  !CALL add_submatrices(1.0_dp,subm_ks_xx_orthog,-1.0_dp,subm_tmp,'N')
805  !CALL add_submatrices(1.0_dp,subm_ks_xx_orthog,-1.0_dp,subm_tmp,'T')
806  !!CALL multiply_submatrices('N','T',1.0_dp,subm_tmp,&
807  !! almo_scf_env%domain_r_down_up(:,ispin),&
808  !! 1.0_dp,subm_ks_xx_orthog)
809 
810  ! convert blocks to the orthogonal basis set
811  ! TRY: replace one multiply
812  !CALL multiply_submatrices('N','N',1.0_dp,subm_ks_xx_orthog,&
813  ! almo_scf_env%domain_s_sqrt_inv(:,ispin),0.0_dp,subm_tmp)
814  CALL multiply_submatrices('N', 'N', 1.0_dp, almo_scf_env%domain_ks_xx(:, ispin), &
815  almo_scf_env%domain_s_sqrt_inv(:, ispin), 0.0_dp, subm_tmp)
816  CALL multiply_submatrices('N', 'N', 1.0_dp, almo_scf_env%domain_s_sqrt_inv(:, ispin), &
817  subm_tmp, 0.0_dp, subm_ks_xx_orthog)
818  CALL release_submatrices(subm_tmp)
819 
820  ! create temporary matrices for occupied and virtual orbitals
821  ! represented in the orthogonalized basis set
822  CALL init_submatrices(subm_t)
823 
824  ! loop over domains - perform diagonalization
825  DO idomain = 1, ndomains
826 
827  ! check if the submatrix exists
828  IF (subm_ks_xx_orthog(idomain)%domain .GT. 0) THEN
829 
830  iblock_size = subm_ks_xx_orthog(idomain)%nrows
831 
832  ! Prepare data
833  ALLOCATE (eigenvalues(iblock_size))
834  ALLOCATE (data_copy(iblock_size, iblock_size))
835  data_copy(:, :) = subm_ks_xx_orthog(idomain)%mdata(:, :)
836 
837  ! Query the optimal workspace for dsyev
838  lwork = -1
839  ALLOCATE (work(max(1, lwork)))
840  CALL dsyev('V', 'L', iblock_size, data_copy, iblock_size, eigenvalues, work, lwork, info)
841  lwork = int(work(1))
842  DEALLOCATE (work)
843 
844  ! Allocate the workspace and solve the eigenproblem
845  ALLOCATE (work(max(1, lwork)))
846  CALL dsyev('V', 'L', iblock_size, data_copy, iblock_size, eigenvalues, work, lwork, info)
847  IF (info .NE. 0) THEN
848  cpabort("DSYEV failed")
849  END IF
850 
851  ! Copy occupied eigenvectors
852  IF (almo_scf_env%domain_t(idomain, ispin)%ncols .NE. &
853  almo_scf_env%nocc_of_domain(idomain, ispin)) THEN
854  cpabort("wrong domain structure")
855  END IF
856  CALL copy_submatrices(almo_scf_env%domain_t(idomain, ispin), &
857  subm_t(idomain), .false.)
858  CALL copy_submatrix_data(data_copy(:, 1:almo_scf_env%nocc_of_domain(idomain, ispin)), &
859  subm_t(idomain))
860  !! Copy occupied eigenvalues if smearing requested
861  IF (almo_scf_env%smear) THEN
862  almo_scf_env%mo_energies(1 + sum(almo_scf_env%nocc_of_domain(:idomain - 1, ispin)) &
863  :sum(almo_scf_env%nocc_of_domain(:idomain, ispin)), ispin) &
864  = eigenvalues(1:almo_scf_env%nocc_of_domain(idomain, ispin))
865  END IF
866 
867  DEALLOCATE (work)
868  DEALLOCATE (data_copy)
869  DEALLOCATE (eigenvalues)
870 
871  END IF ! submatrix for the domain exists
872 
873  END DO ! loop over domains
874 
875  CALL release_submatrices(subm_ks_xx_orthog)
876 
877  ! convert orbitals to the AO basis set (from orthogonalized AOs)
878  CALL multiply_submatrices('N', 'N', 1.0_dp, almo_scf_env%domain_s_sqrt_inv(:, ispin), &
879  subm_t, 0.0_dp, almo_scf_env%domain_t(:, ispin))
880  CALL release_submatrices(subm_t)
881 
882  ! convert domain orbitals to a dbcsr matrix
883  CALL construct_dbcsr_from_submatrices( &
884  almo_scf_env%matrix_t(ispin), &
885  almo_scf_env%domain_t(:, ispin), &
886  almo_scf_env%quench_t(ispin))
887  CALL dbcsr_filter(almo_scf_env%matrix_t(ispin), &
888  almo_scf_env%eps_filter)
889 
890  ! TRY: add T0 component
891  !!CALL dbcsr_add(almo_scf_env%matrix_t(ispin),&
892  !! almo_scf_env%matrix_t_blk(ispin),1.0_dp,1.0_dp)
893 
894  END DO ! spins
895 
896  DEALLOCATE (subm_tmp)
897  DEALLOCATE (subm_ks_xx_orthog)
898  DEALLOCATE (subm_t)
899 
900  CALL timestop(handle)
901 
902  END SUBROUTINE almo_scf_ks_xx_to_tv_xx
903 
904 ! **************************************************************************************************
905 !> \brief computes ALMOs by diagonalizing the projected blocked KS matrix
906 !> uses the diagonalization code for blocks
907 !> computes both the occupied and virtual orbitals
908 !> \param almo_scf_env ...
909 !> \par History
910 !> 2011.07 created [Rustam Z Khaliullin]
911 !> 2018.09 smearing support [Ruben Staub]
912 !> \author Rustam Z Khaliullin
913 ! **************************************************************************************************
914  SUBROUTINE almo_scf_ks_blk_to_tv_blk(almo_scf_env)
915 
916  TYPE(almo_scf_env_type), INTENT(INOUT) :: almo_scf_env
917 
918  CHARACTER(LEN=*), PARAMETER :: routinen = 'almo_scf_ks_blk_to_tv_blk'
919 
920  INTEGER :: handle, iblock_col, iblock_row, &
921  iblock_size, info, ispin, lwork, &
922  nocc_of_block, nvirt_of_block, orbital
923  LOGICAL :: block_needed
924  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: eigenvalues, work
925  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: data_copy
926  REAL(kind=dp), DIMENSION(:, :), POINTER :: data_p, p_new_block
927  TYPE(dbcsr_iterator_type) :: iter
928  TYPE(dbcsr_type) :: matrix_ks_blk_orthog, &
929  matrix_t_blk_orthog, matrix_tmp, &
930  matrix_v_blk_orthog
931 
932  CALL timeset(routinen, handle)
933 
934  IF (almo_scf_env%domain_layout_aos == almo_domain_layout_molecular .AND. &
935  almo_scf_env%mat_distr_aos == almo_mat_distr_atomic) THEN
936  cpabort("a domain must be located entirely on a CPU")
937  END IF
938 
939  DO ispin = 1, almo_scf_env%nspins
940 
941  CALL dbcsr_create(matrix_tmp, template=almo_scf_env%matrix_ks_blk(ispin), &
942  matrix_type=dbcsr_type_no_symmetry)
943  CALL dbcsr_create(matrix_ks_blk_orthog, template=almo_scf_env%matrix_ks_blk(ispin), &
944  matrix_type=dbcsr_type_no_symmetry)
945 
946  ! convert blocks to the orthogonal basis set
947  CALL dbcsr_multiply("N", "N", 1.0_dp, almo_scf_env%matrix_ks_blk(ispin), &
948  almo_scf_env%matrix_s_blk_sqrt_inv(1), 0.0_dp, matrix_tmp, &
949  filter_eps=almo_scf_env%eps_filter)
950  CALL dbcsr_multiply("N", "N", 1.0_dp, almo_scf_env%matrix_s_blk_sqrt_inv(1), &
951  matrix_tmp, 0.0_dp, matrix_ks_blk_orthog, &
952  filter_eps=almo_scf_env%eps_filter)
953 
954  CALL dbcsr_release(matrix_tmp)
955 
956  ! create temporary matrices for occupied and virtual orbitals
957  ! represented in the orthogonalized AOs basis set
958  CALL dbcsr_create(matrix_t_blk_orthog, template=almo_scf_env%matrix_t_blk(ispin))
959  CALL dbcsr_create(matrix_v_blk_orthog, template=almo_scf_env%matrix_v_full_blk(ispin))
960  CALL dbcsr_work_create(matrix_t_blk_orthog, work_mutable=.true.)
961  CALL dbcsr_work_create(matrix_v_blk_orthog, work_mutable=.true.)
962 
963  CALL dbcsr_work_create(almo_scf_env%matrix_eoo(ispin), work_mutable=.true.)
964  CALL dbcsr_work_create(almo_scf_env%matrix_evv_full(ispin), work_mutable=.true.)
965 
966  CALL dbcsr_iterator_start(iter, matrix_ks_blk_orthog)
967 
968  DO WHILE (dbcsr_iterator_blocks_left(iter))
969  CALL dbcsr_iterator_next_block(iter, iblock_row, iblock_col, data_p, row_size=iblock_size)
970 
971  IF (iblock_row .NE. iblock_col) THEN
972  cpabort("off-diagonal block found")
973  END IF
974 
975  block_needed = .true.
976  IF (almo_scf_env%nocc_of_domain(iblock_col, ispin) .EQ. 0 .AND. &
977  almo_scf_env%nvirt_of_domain(iblock_col, ispin) .EQ. 0) THEN
978  block_needed = .false.
979  END IF
980 
981  IF (block_needed) THEN
982 
983  ! Prepare data
984  ALLOCATE (eigenvalues(iblock_size))
985  ALLOCATE (data_copy(iblock_size, iblock_size))
986  data_copy(:, :) = data_p(:, :)
987 
988  ! Query the optimal workspace for dsyev
989  lwork = -1
990  ALLOCATE (work(max(1, lwork)))
991  CALL dsyev('V', 'L', iblock_size, data_copy, iblock_size, eigenvalues, work, lwork, info)
992  lwork = int(work(1))
993  DEALLOCATE (work)
994 
995  ! Allocate the workspace and solve the eigenproblem
996  ALLOCATE (work(max(1, lwork)))
997  CALL dsyev('V', 'L', iblock_size, data_copy, iblock_size, eigenvalues, work, lwork, info)
998  IF (info .NE. 0) THEN
999  cpabort("DSYEV failed")
1000  END IF
1001 
1002  !!! RZK-warning !!!
1003  !!! IT IS EXTREMELY IMPORTANT THAT THE DIAGONAL BLOCKS OF THE !!!
1004  !!! FOLLOWING MATRICES ARE LOCATED ON THE SAME NODES WITH !!!
1005  !!! THE CORRESPONDING DIAGONAL BLOCKS OF THE FOCK MATRIX: !!!
1006  !!! T, V, E_o, E_v
1007 
1008  ! copy eigenvectors into two dbcsr matrices - occupied and virtuals
1009  nocc_of_block = almo_scf_env%nocc_of_domain(iblock_col, ispin)
1010  IF (nocc_of_block .GT. 0) THEN
1011  NULLIFY (p_new_block)
1012  CALL dbcsr_reserve_block2d(matrix_t_blk_orthog, iblock_row, iblock_col, p_new_block)
1013  cpassert(ASSOCIATED(p_new_block))
1014  p_new_block(:, :) = data_copy(:, 1:nocc_of_block)
1015  ! copy eigenvalues into diagonal dbcsr matrix - Eoo
1016  NULLIFY (p_new_block)
1017  CALL dbcsr_reserve_block2d(almo_scf_env%matrix_eoo(ispin), iblock_row, iblock_col, p_new_block)
1018  cpassert(ASSOCIATED(p_new_block))
1019  p_new_block(:, :) = 0.0_dp
1020  DO orbital = 1, nocc_of_block
1021  p_new_block(orbital, orbital) = eigenvalues(orbital)
1022  END DO
1023  !! Retrieve occupied MOs energies for smearing purpose, if requested
1024  !! RS-WARNING: Hack to retrieve the occupied energies, since matrix_eoo seems to be ill-defined
1025  !! for multiprocessing (any idea for fix?)
1026  !! RS-WARNING: This section is not suitable for parallel run !!!
1027  !! (but usually fails less than retrieving the diagonal of matrix_eoo)
1028  !! RS-WARNING: This method will likely keep the energies of the initial guess if run in parallel
1029  !! (which is still a reasonable smearing in most cases...)
1030  IF (almo_scf_env%smear) THEN
1031  DO orbital = 1, nocc_of_block
1032  almo_scf_env%mo_energies(sum(almo_scf_env%nocc_of_domain(:iblock_row - 1, ispin)) + orbital, &
1033  ispin) = eigenvalues(orbital)
1034  END DO
1035  END IF
1036  END IF
1037 
1038  ! now virtuals
1039  nvirt_of_block = almo_scf_env%nvirt_of_domain(iblock_col, ispin)
1040  IF (nvirt_of_block .GT. 0) THEN
1041  NULLIFY (p_new_block)
1042  CALL dbcsr_reserve_block2d(matrix_v_blk_orthog, iblock_row, iblock_col, p_new_block)
1043  cpassert(ASSOCIATED(p_new_block))
1044  p_new_block(:, :) = data_copy(:, (nocc_of_block + 1):(nocc_of_block + nvirt_of_block))
1045  ! virtual energies
1046  NULLIFY (p_new_block)
1047  CALL dbcsr_reserve_block2d(almo_scf_env%matrix_evv_full(ispin), iblock_row, iblock_col, p_new_block)
1048  cpassert(ASSOCIATED(p_new_block))
1049  p_new_block(:, :) = 0.0_dp
1050  DO orbital = 1, nvirt_of_block
1051  p_new_block(orbital, orbital) = eigenvalues(nocc_of_block + orbital)
1052  END DO
1053  END IF
1054 
1055  DEALLOCATE (work)
1056  DEALLOCATE (data_copy)
1057  DEALLOCATE (eigenvalues)
1058 
1059  END IF
1060 
1061  END DO
1062  CALL dbcsr_iterator_stop(iter)
1063 
1064  CALL dbcsr_finalize(matrix_t_blk_orthog)
1065  CALL dbcsr_finalize(matrix_v_blk_orthog)
1066  CALL dbcsr_finalize(almo_scf_env%matrix_eoo(ispin))
1067  CALL dbcsr_finalize(almo_scf_env%matrix_evv_full(ispin))
1068 
1069  !! RS-WARNING: When matrix_eoo will be well-defined with multiprocessing,
1070  !! the following should be the preferred way to retrieve the occupied energies:
1071  !! Retrieve occupied MOs energies for smearing purpose, if requested
1072  !! IF (almo_scf_env%smear) THEN
1073  !! CALL dbcsr_get_diag(almo_scf_env%matrix_eoo(ispin), almo_scf_env%mo_energies(:, ispin))
1074  !! END IF
1075 
1076  CALL dbcsr_filter(matrix_t_blk_orthog, almo_scf_env%eps_filter)
1077  CALL dbcsr_filter(matrix_v_blk_orthog, almo_scf_env%eps_filter)
1078 
1079  CALL dbcsr_release(matrix_ks_blk_orthog)
1080 
1081  ! convert orbitals to the AO basis set (from orthogonalized AOs)
1082  CALL dbcsr_multiply("N", "N", 1.0_dp, almo_scf_env%matrix_s_blk_sqrt_inv(1), &
1083  matrix_t_blk_orthog, 0.0_dp, almo_scf_env%matrix_t_blk(ispin), &
1084  filter_eps=almo_scf_env%eps_filter)
1085  CALL dbcsr_multiply("N", "N", 1.0_dp, almo_scf_env%matrix_s_blk_sqrt_inv(1), &
1086  matrix_v_blk_orthog, 0.0_dp, almo_scf_env%matrix_v_full_blk(ispin), &
1087  filter_eps=almo_scf_env%eps_filter)
1088 
1089  CALL dbcsr_release(matrix_t_blk_orthog)
1090  CALL dbcsr_release(matrix_v_blk_orthog)
1091 
1092  END DO ! spins
1093 
1094  CALL timestop(handle)
1095 
1096  END SUBROUTINE almo_scf_ks_blk_to_tv_blk
1097 
1098 ! **************************************************************************************************
1099 !> \brief inverts block-diagonal blocks of a dbcsr_matrix
1100 !> \param matrix_in ...
1101 !> \param matrix_out ...
1102 !> \param nocc ...
1103 !> \par History
1104 !> 2012.05 created [Rustam Z Khaliullin]
1105 !> \author Rustam Z Khaliullin
1106 ! **************************************************************************************************
1107  SUBROUTINE pseudo_invert_diagonal_blk(matrix_in, matrix_out, nocc)
1108 
1109  TYPE(dbcsr_type), INTENT(IN) :: matrix_in
1110  TYPE(dbcsr_type), INTENT(INOUT) :: matrix_out
1111  INTEGER, DIMENSION(:) :: nocc
1112 
1113  CHARACTER(LEN=*), PARAMETER :: routinen = 'pseudo_invert_diagonal_blk'
1114 
1115  INTEGER :: handle, iblock_col, iblock_row, &
1116  iblock_size, methodid
1117  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: data_copy
1118  REAL(kind=dp), DIMENSION(:, :), POINTER :: data_p, p_new_block
1119  TYPE(dbcsr_iterator_type) :: iter
1120 
1121  CALL timeset(routinen, handle)
1122 
1123  CALL dbcsr_create(matrix_out, template=matrix_in)
1124  CALL dbcsr_work_create(matrix_out, work_mutable=.true.)
1125 
1126  CALL dbcsr_iterator_start(iter, matrix_in)
1127 
1128  DO WHILE (dbcsr_iterator_blocks_left(iter))
1129 
1130  CALL dbcsr_iterator_next_block(iter, iblock_row, iblock_col, data_p, row_size=iblock_size)
1131 
1132  IF (iblock_row == iblock_col) THEN
1133 
1134  ! Prepare data
1135  ALLOCATE (data_copy(iblock_size, iblock_size))
1136  !data_copy(:,:)=data_p(:,:)
1137 
1138  ! 0. Cholesky factorization
1139  ! 1. Diagonalization
1140  methodid = 1
1141  CALL pseudo_invert_matrix(data_p, data_copy, iblock_size, &
1142  methodid, &
1143  range1=nocc(iblock_row), range2=nocc(iblock_row), &
1144  !range1_thr,range2_thr,&
1145  shift=1.0e-5_dp)
1146  !!! IT IS EXTREMELY IMPORTANT THAT THE BLOCKS OF THE "OUT" !!!
1147  !!! MATRIX ARE DISTRIBUTED AS THE BLOCKS OF THE "IN" MATRIX !!!
1148 
1149  NULLIFY (p_new_block)
1150  CALL dbcsr_reserve_block2d(matrix_out, iblock_row, iblock_col, p_new_block)
1151  cpassert(ASSOCIATED(p_new_block))
1152  p_new_block(:, :) = data_copy(:, :)
1153 
1154  DEALLOCATE (data_copy)
1155 
1156  END IF
1157 
1158  END DO
1159  CALL dbcsr_iterator_stop(iter)
1160 
1161  CALL dbcsr_finalize(matrix_out)
1162 
1163  CALL timestop(handle)
1164 
1165  END SUBROUTINE pseudo_invert_diagonal_blk
1166 
1167 ! **************************************************************************************************
1168 !> \brief computes occupied ALMOs from the superimposed atomic density blocks
1169 !> \param almo_scf_env ...
1170 !> \param ionic ...
1171 !> \par History
1172 !> 2011.06 created [Rustam Z Khaliullin]
1173 !> \author Rustam Z Khaliullin
1174 ! **************************************************************************************************
1175  SUBROUTINE almo_scf_p_blk_to_t_blk(almo_scf_env, ionic)
1176 
1177  TYPE(almo_scf_env_type), INTENT(INOUT) :: almo_scf_env
1178  LOGICAL, INTENT(IN) :: ionic
1179 
1180  CHARACTER(LEN=*), PARAMETER :: routinen = 'almo_scf_p_blk_to_t_blk'
1181 
1182  INTEGER :: handle, iblock_col, iblock_row, &
1183  iblock_size, info, ispin, lwork, &
1184  nocc_of_block, unit_nr
1185  LOGICAL :: block_needed
1186  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: eigenvalues, work
1187  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: data_copy
1188  REAL(kind=dp), DIMENSION(:, :), POINTER :: data_p, p_new_block
1189  TYPE(cp_logger_type), POINTER :: logger
1190  TYPE(dbcsr_iterator_type) :: iter
1191  TYPE(dbcsr_type) :: matrix_t_blk_tmp
1192 
1193  CALL timeset(routinen, handle)
1194 
1195  ! get a useful unit_nr
1196  logger => cp_get_default_logger()
1197  IF (logger%para_env%is_source()) THEN
1198  unit_nr = cp_logger_get_default_unit_nr(logger, local=.true.)
1199  ELSE
1200  unit_nr = -1
1201  END IF
1202 
1203  DO ispin = 1, almo_scf_env%nspins
1204 
1205  IF (ionic) THEN
1206 
1207  ! create a temporary matrix to keep the eigenvectors
1208  CALL dbcsr_create(matrix_t_blk_tmp, &
1209  template=almo_scf_env%matrix_t_blk(ispin))
1210  CALL dbcsr_work_create(matrix_t_blk_tmp, &
1211  work_mutable=.true.)
1212 
1213  CALL dbcsr_iterator_start(iter, almo_scf_env%matrix_p_blk(ispin))
1214  DO WHILE (dbcsr_iterator_blocks_left(iter))
1215  CALL dbcsr_iterator_next_block(iter, iblock_row, iblock_col, data_p, row_size=iblock_size)
1216 
1217  block_needed = .false.
1218 
1219  IF (iblock_row == iblock_col) THEN
1220  block_needed = .true.
1221  END IF
1222 
1223  IF (.NOT. block_needed) THEN
1224  cpabort("off-diag block found")
1225  END IF
1226 
1227  ! Prepare data
1228  ALLOCATE (eigenvalues(iblock_size))
1229  ALLOCATE (data_copy(iblock_size, iblock_size))
1230  data_copy(:, :) = data_p(:, :)
1231 
1232  ! Query the optimal workspace for dsyev
1233  lwork = -1
1234  ALLOCATE (work(max(1, lwork)))
1235  CALL dsyev('V', 'L', iblock_size, data_copy, iblock_size, eigenvalues, &
1236  work, lwork, info)
1237  lwork = int(work(1))
1238  DEALLOCATE (work)
1239 
1240  ! Allocate the workspace and solve the eigenproblem
1241  ALLOCATE (work(max(1, lwork)))
1242  CALL dsyev('V', 'L', iblock_size, data_copy, iblock_size, eigenvalues, work, lwork, info)
1243  IF (info .NE. 0) THEN
1244  IF (unit_nr > 0) THEN
1245  WRITE (unit_nr, *) 'BLOCK = ', iblock_row
1246  WRITE (unit_nr, *) 'INFO =', info
1247  WRITE (unit_nr, *) data_p(:, :)
1248  END IF
1249  cpabort("DSYEV failed")
1250  END IF
1251 
1252  !!! IT IS EXTREMELY IMPORTANT THAT THE DIAGONAL BLOCKS OF THE !!!
1253  !!! P AND T MATRICES ARE LOCATED ON THE SAME NODES !!!
1254 
1255  ! copy eigenvectors into two dbcsr matrices - occupied and virtuals
1256  NULLIFY (p_new_block)
1257  CALL dbcsr_reserve_block2d(matrix_t_blk_tmp, &
1258  iblock_row, iblock_col, p_new_block)
1259  nocc_of_block = SIZE(p_new_block, 2)
1260  cpassert(ASSOCIATED(p_new_block))
1261  cpassert(nocc_of_block .GT. 0)
1262  p_new_block(:, :) = data_copy(:, iblock_size - nocc_of_block + 1:)
1263 
1264  DEALLOCATE (work)
1265  DEALLOCATE (data_copy)
1266  DEALLOCATE (eigenvalues)
1267 
1268  END DO
1269  CALL dbcsr_iterator_stop(iter)
1270 
1271  CALL dbcsr_finalize(matrix_t_blk_tmp)
1272  CALL dbcsr_filter(matrix_t_blk_tmp, &
1273  almo_scf_env%eps_filter)
1274  CALL dbcsr_copy(almo_scf_env%matrix_t_blk(ispin), &
1275  matrix_t_blk_tmp)
1276  CALL dbcsr_release(matrix_t_blk_tmp)
1277 
1278  ELSE
1279 
1280  !! generate a random set of ALMOs
1281  !! matrix_t_blk should already be initiated to the proper domain structure
1282  CALL dbcsr_init_random(almo_scf_env%matrix_t_blk(ispin), &
1283  keep_sparsity=.true.)
1284 
1285  CALL dbcsr_create(matrix_t_blk_tmp, &
1286  template=almo_scf_env%matrix_t_blk(ispin), &
1287  matrix_type=dbcsr_type_no_symmetry)
1288 
1289  ! use current ALMOs in matrix_t_blk and project them onto the blocked dm
1290  ! compute T_new = R_blk S_blk T_random
1291  CALL dbcsr_multiply("N", "N", 1.0_dp, almo_scf_env%matrix_s_blk(1), &
1292  almo_scf_env%matrix_t_blk(ispin), &
1293  0.0_dp, matrix_t_blk_tmp, &
1294  filter_eps=almo_scf_env%eps_filter)
1295 
1296  CALL dbcsr_multiply("N", "N", 1.0_dp, &
1297  almo_scf_env%matrix_p_blk(ispin), matrix_t_blk_tmp, &
1298  0.0_dp, almo_scf_env%matrix_t_blk(ispin), &
1299  filter_eps=almo_scf_env%eps_filter)
1300 
1301  CALL dbcsr_release(matrix_t_blk_tmp)
1302 
1303  END IF
1304 
1305  END DO
1306 
1307  CALL timestop(handle)
1308 
1309  END SUBROUTINE almo_scf_p_blk_to_t_blk
1310 
1311 ! **************************************************************************************************
1312 !> \brief Apply an occupation-rescaling trick to ALMOs for smearing.
1313 !> Partially occupied orbitals are considered full and rescaled by SQRT(occupation_number)
1314 !> (this was designed to be used with smearing only)
1315 !> \param matrix_t ...
1316 !> \param mo_energies ...
1317 !> \param mu_of_domain ...
1318 !> \param real_ne_of_domain ...
1319 !> \param spin_kTS ...
1320 !> \param smear_e_temp ...
1321 !> \param ndomains ...
1322 !> \param nocc_of_domain ...
1323 !> \par History
1324 !> 2018.09 created [Ruben Staub]
1325 !> \author Ruben Staub
1326 ! **************************************************************************************************
1327  SUBROUTINE almo_scf_t_rescaling(matrix_t, mo_energies, mu_of_domain, real_ne_of_domain, &
1328  spin_kTS, smear_e_temp, ndomains, nocc_of_domain)
1329 
1330  TYPE(dbcsr_type), INTENT(INOUT) :: matrix_t
1331  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: mo_energies
1332  REAL(kind=dp), DIMENSION(:), INTENT(INOUT) :: mu_of_domain, real_ne_of_domain
1333  REAL(kind=dp), INTENT(INOUT) :: spin_kts
1334  REAL(kind=dp), INTENT(IN) :: smear_e_temp
1335  INTEGER, INTENT(IN) :: ndomains
1336  INTEGER, DIMENSION(:), INTENT(IN) :: nocc_of_domain
1337 
1338  CHARACTER(LEN=*), PARAMETER :: routinen = 'almo_scf_t_rescaling'
1339 
1340  INTEGER :: handle, idomain, neigenval_used, nmo
1341  REAL(kind=dp) :: kts
1342  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: occupation_numbers, rescaling_factors
1343 
1344  CALL timeset(routinen, handle)
1345 
1346  !!
1347  !! Initialization
1348  !!
1349  nmo = SIZE(mo_energies)
1350  ALLOCATE (occupation_numbers(nmo))
1351  ALLOCATE (rescaling_factors(nmo))
1352 
1353  !!
1354  !! Set occupation numbers for smearing
1355  !!
1356  !! occupation numbers are obtained using Fermi-Dirac smearing with orbital energies stored in mo_energies
1357  !! nocc_of_domain is the number of partially occupied orbitals, while real_ne_of_domain is the real number of electrons
1358  neigenval_used = 0 !! this is used as an offset to copy sub-arrays
1359 
1360  !! Reset electronic entropy term
1361  spin_kts = 0.0_dp
1362 
1363  !! Apply Fermi-Dirac smearing for each domain and store associated occupations for the whole system
1364  DO idomain = 1, ndomains
1365  CALL fermifixed(occupation_numbers(1 + neigenval_used:nocc_of_domain(idomain) + neigenval_used), &
1366  mu_of_domain(idomain), &
1367  kts, &
1368  mo_energies(1 + neigenval_used:nocc_of_domain(idomain) + neigenval_used), &
1369  real_ne_of_domain(idomain), &
1370  smear_e_temp, &
1371  1.0_dp) !! Warning, maxocc is set to 1 since we don't want to interfere with the spin_factor rescaling
1372  spin_kts = spin_kts + kts !! Add up electronic entropy contributions
1373  neigenval_used = neigenval_used + nocc_of_domain(idomain) !! Update eigenvalues index offset
1374  END DO
1375  rescaling_factors(:) = sqrt(occupation_numbers) !! scale = sqrt(occupation_number)
1376 
1377  !!
1378  !! Rescaling electronic entropy contribution by spin_factor (deprecated)
1379  !! (currently, entropy is rescaled by spin_factor with the density matrix)
1380  !!
1381  !!IF (almo_scf_env%nspins == 1) THEN
1382  !! spin_kTS = spin_kTS*2.0_dp
1383  !!END IF
1384 
1385  !!
1386  !! Rescaling of T (i.e. ALMOs)
1387  !!
1388  CALL dbcsr_scale_by_vector(matrix_t, rescaling_factors, side='right') !! Apply occupation-rescaling trick
1389 
1390  !!
1391  !! Debug tools (for debug purpose only)
1392  !!
1393  !! WRITE (*,*) "occupations", occupation_numbers(:) !! debug
1394  !! WRITE (*,*) "eigenvalues", mo_energies(:) !! debug
1395  !! WRITE (*,*) "kTS (spin_factor excluded) = ", spin_kTS !! debug
1396 
1397  !!
1398  !! Cleaning up before exit
1399  !!
1400  DEALLOCATE (occupation_numbers)
1401  DEALLOCATE (rescaling_factors)
1402 
1403  CALL timestop(handle)
1404 
1405  END SUBROUTINE almo_scf_t_rescaling
1406 
1407 ! **************************************************************************************************
1408 !> \brief Computes the overlap matrix of MO orbitals
1409 !> \param bra ...
1410 !> \param ket ...
1411 !> \param overlap ...
1412 !> \param metric ...
1413 !> \param retain_overlap_sparsity ...
1414 !> \param eps_filter ...
1415 !> \param smear ...
1416 !> \par History
1417 !> 2011.08 created [Rustam Z Khaliullin]
1418 !> 2018.09 smearing support [Ruben Staub]
1419 !> \author Rustam Z Khaliullin
1420 ! **************************************************************************************************
1421  SUBROUTINE get_overlap(bra, ket, overlap, metric, retain_overlap_sparsity, &
1422  eps_filter, smear)
1423 
1424  TYPE(dbcsr_type), INTENT(IN) :: bra, ket
1425  TYPE(dbcsr_type), INTENT(INOUT) :: overlap
1426  TYPE(dbcsr_type), INTENT(IN) :: metric
1427  LOGICAL, INTENT(IN), OPTIONAL :: retain_overlap_sparsity
1428  REAL(kind=dp) :: eps_filter
1429  LOGICAL, INTENT(IN), OPTIONAL :: smear
1430 
1431  CHARACTER(LEN=*), PARAMETER :: routinen = 'get_overlap'
1432 
1433  INTEGER :: dim0, handle
1434  LOGICAL :: local_retain_sparsity, smearing
1435  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: diag_correction
1436  TYPE(dbcsr_type) :: tmp
1437 
1438  CALL timeset(routinen, handle)
1439 
1440  IF (.NOT. PRESENT(retain_overlap_sparsity)) THEN
1441  local_retain_sparsity = .false.
1442  ELSE
1443  local_retain_sparsity = retain_overlap_sparsity
1444  END IF
1445 
1446  IF (.NOT. PRESENT(smear)) THEN
1447  smearing = .false.
1448  ELSE
1449  smearing = smear
1450  END IF
1451 
1452  CALL dbcsr_create(tmp, template=ket, &
1453  matrix_type=dbcsr_type_no_symmetry)
1454 
1455  ! TMP=metric*ket
1456  CALL dbcsr_multiply("N", "N", 1.0_dp, &
1457  metric, ket, 0.0_dp, tmp, &
1458  filter_eps=eps_filter)
1459 
1460  ! OVERLAP=tr(bra)*TMP
1461  CALL dbcsr_multiply("T", "N", 1.0_dp, &
1462  bra, tmp, 0.0_dp, overlap, &
1463  retain_sparsity=local_retain_sparsity, &
1464  filter_eps=eps_filter)
1465 
1466  CALL dbcsr_release(tmp)
1467 
1468  !! If smearing ALMO is requested, apply correction of the occupation-rescaling trick
1469  !! (i.e. converting rescaled orbitals into selfish orbitals)
1470  !! (i.e. the diagonal blocks of the rescaled overlap must remain unscaled)
1471  !! Since we have here orthonormal MOs within a fragment, diagonal blocks are identity matrices
1472  !! Therefore, one only need to restore the diagonal to 1
1473  !! RS-WARNING: Assume orthonormal MOs within a fragment
1474  IF (smearing) THEN
1475  CALL dbcsr_get_info(overlap, nfullrows_total=dim0)
1476  ALLOCATE (diag_correction(dim0))
1477  diag_correction = 1.0_dp
1478  CALL dbcsr_set_diag(overlap, diag_correction)
1479  DEALLOCATE (diag_correction)
1480  END IF
1481 
1482  CALL timestop(handle)
1483 
1484  END SUBROUTINE get_overlap
1485 
1486 ! **************************************************************************************************
1487 !> \brief orthogonalize MOs
1488 !> \param ket ...
1489 !> \param overlap ...
1490 !> \param metric ...
1491 !> \param retain_locality ...
1492 !> \param only_normalize ...
1493 !> \param nocc_of_domain ...
1494 !> \param eps_filter ...
1495 !> \param order_lanczos ...
1496 !> \param eps_lanczos ...
1497 !> \param max_iter_lanczos ...
1498 !> \param overlap_sqrti ...
1499 !> \param smear ...
1500 !> \par History
1501 !> 2012.03 created [Rustam Z Khaliullin]
1502 !> 2018.09 smearing support [Ruben Staub]
1503 !> \author Rustam Z Khaliullin
1504 ! **************************************************************************************************
1505  SUBROUTINE orthogonalize_mos(ket, overlap, metric, retain_locality, only_normalize, &
1506  nocc_of_domain, eps_filter, order_lanczos, eps_lanczos, &
1507  max_iter_lanczos, overlap_sqrti, smear)
1508 
1509  TYPE(dbcsr_type), INTENT(INOUT) :: ket, overlap
1510  TYPE(dbcsr_type), INTENT(IN) :: metric
1511  LOGICAL, INTENT(IN) :: retain_locality, only_normalize
1512  INTEGER, DIMENSION(:), INTENT(IN) :: nocc_of_domain
1513  REAL(kind=dp) :: eps_filter
1514  INTEGER, INTENT(IN) :: order_lanczos
1515  REAL(kind=dp), INTENT(IN) :: eps_lanczos
1516  INTEGER, INTENT(IN) :: max_iter_lanczos
1517  TYPE(dbcsr_type), INTENT(INOUT), OPTIONAL :: overlap_sqrti
1518  LOGICAL, INTENT(IN), OPTIONAL :: smear
1519 
1520  CHARACTER(LEN=*), PARAMETER :: routinen = 'orthogonalize_mos'
1521 
1522  INTEGER :: dim0, handle
1523  LOGICAL :: smearing
1524  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: diagonal
1525  TYPE(dbcsr_type) :: matrix_sigma_blk_sqrt, &
1526  matrix_sigma_blk_sqrt_inv, &
1527  matrix_t_blk_tmp
1528 
1529  CALL timeset(routinen, handle)
1530 
1531  IF (.NOT. PRESENT(smear)) THEN
1532  smearing = .false.
1533  ELSE
1534  smearing = smear
1535  END IF
1536 
1537  ! create block-diagonal sparsity pattern for the overlap
1538  ! in case retain_locality is set to true
1539  ! RZK-warning this will fail if distribution blocks are smaller than domains!!!
1540  CALL dbcsr_set(overlap, 0.0_dp)
1541  CALL dbcsr_add_on_diag(overlap, 1.0_dp)
1542  CALL dbcsr_filter(overlap, eps_filter)
1543 
1544  CALL get_overlap(ket, ket, overlap, metric, retain_locality, &
1545  eps_filter, smear=smearing)
1546 
1547  IF (only_normalize) THEN
1548 
1549  CALL dbcsr_get_info(overlap, nfullrows_total=dim0)
1550  ALLOCATE (diagonal(dim0))
1551  CALL dbcsr_get_diag(overlap, diagonal)
1552  CALL dbcsr_set(overlap, 0.0_dp)
1553  CALL dbcsr_set_diag(overlap, diagonal)
1554  DEALLOCATE (diagonal)
1555  CALL dbcsr_filter(overlap, eps_filter)
1556 
1557  END IF
1558 
1559  CALL dbcsr_create(matrix_sigma_blk_sqrt, template=overlap, &
1560  matrix_type=dbcsr_type_no_symmetry)
1561  CALL dbcsr_create(matrix_sigma_blk_sqrt_inv, template=overlap, &
1562  matrix_type=dbcsr_type_no_symmetry)
1563 
1564  ! compute sqrt and sqrt_inv of the blocked MO overlap
1565  CALL set_zero_electron_blocks_in_mo_mo_matrix(overlap, nocc_of_domain, 1.0_dp)
1566  CALL matrix_sqrt_newton_schulz(matrix_sigma_blk_sqrt, matrix_sigma_blk_sqrt_inv, &
1567  overlap, threshold=eps_filter, &
1568  order=order_lanczos, &
1569  eps_lanczos=eps_lanczos, &
1570  max_iter_lanczos=max_iter_lanczos)
1571  CALL set_zero_electron_blocks_in_mo_mo_matrix(overlap, nocc_of_domain, 0.0_dp)
1572  !CALL set_zero_electron_blocks_in_mo_mo_matrix(matrix_sigma_blk_sqrt,nocc_of_domain,0.0_dp)
1573  CALL set_zero_electron_blocks_in_mo_mo_matrix(matrix_sigma_blk_sqrt_inv, nocc_of_domain, 0.0_dp)
1574 
1575  CALL dbcsr_create(matrix_t_blk_tmp, &
1576  template=ket, &
1577  matrix_type=dbcsr_type_no_symmetry)
1578 
1579  CALL dbcsr_multiply("N", "N", 1.0_dp, &
1580  ket, &
1581  matrix_sigma_blk_sqrt_inv, &
1582  0.0_dp, matrix_t_blk_tmp, &
1583  filter_eps=eps_filter)
1584 
1585  ! update the orbitals with the orthonormalized MOs
1586  CALL dbcsr_copy(ket, matrix_t_blk_tmp)
1587 
1588  ! return overlap SQRT_INV if necessary
1589  IF (PRESENT(overlap_sqrti)) THEN
1590  CALL dbcsr_copy(overlap_sqrti, &
1591  matrix_sigma_blk_sqrt_inv)
1592  END IF
1593 
1594  CALL dbcsr_release(matrix_t_blk_tmp)
1595  CALL dbcsr_release(matrix_sigma_blk_sqrt)
1596  CALL dbcsr_release(matrix_sigma_blk_sqrt_inv)
1597 
1598  CALL timestop(handle)
1599 
1600  END SUBROUTINE orthogonalize_mos
1601 
1602 ! **************************************************************************************************
1603 !> \brief computes the idempotent density matrix from MOs
1604 !> MOs can be either orthogonal or non-orthogonal
1605 !> \param t ...
1606 !> \param p ...
1607 !> \param eps_filter ...
1608 !> \param orthog_orbs ...
1609 !> \param nocc_of_domain ...
1610 !> \param s ...
1611 !> \param sigma ...
1612 !> \param sigma_inv ...
1613 !> \param use_guess ...
1614 !> \param smear ...
1615 !> \param algorithm to inver sigma: 0 - Hotelling (linear), 1 - Cholesky (cubic, low prefactor)
1616 !> \param para_env ...
1617 !> \param blacs_env ...
1618 !> \param eps_lanczos ...
1619 !> \param max_iter_lanczos ...
1620 !> \param inverse_accelerator ...
1621 !> \param inv_eps_factor ...
1622 !> \par History
1623 !> 2011.07 created [Rustam Z Khaliullin]
1624 !> 2018.09 smearing support [Ruben Staub]
1625 !> \author Rustam Z Khaliullin
1626 ! **************************************************************************************************
1627  SUBROUTINE almo_scf_t_to_proj(t, p, eps_filter, orthog_orbs, nocc_of_domain, s, sigma, sigma_inv, &
1628  use_guess, smear, algorithm, para_env, blacs_env, eps_lanczos, &
1629  max_iter_lanczos, inverse_accelerator, inv_eps_factor)
1630 
1631  TYPE(dbcsr_type), INTENT(IN) :: t
1632  TYPE(dbcsr_type), INTENT(INOUT) :: p
1633  REAL(kind=dp), INTENT(IN) :: eps_filter
1634  LOGICAL, INTENT(IN) :: orthog_orbs
1635  INTEGER, DIMENSION(:), INTENT(IN), OPTIONAL :: nocc_of_domain
1636  TYPE(dbcsr_type), INTENT(IN), OPTIONAL :: s
1637  TYPE(dbcsr_type), INTENT(INOUT), OPTIONAL :: sigma, sigma_inv
1638  LOGICAL, INTENT(IN), OPTIONAL :: use_guess, smear
1639  INTEGER, INTENT(IN), OPTIONAL :: algorithm
1640  TYPE(mp_para_env_type), OPTIONAL, POINTER :: para_env
1641  TYPE(cp_blacs_env_type), OPTIONAL, POINTER :: blacs_env
1642  REAL(kind=dp), INTENT(IN), OPTIONAL :: eps_lanczos
1643  INTEGER, INTENT(IN), OPTIONAL :: max_iter_lanczos, inverse_accelerator
1644  REAL(kind=dp), INTENT(IN), OPTIONAL :: inv_eps_factor
1645 
1646  CHARACTER(LEN=*), PARAMETER :: routinen = 'almo_scf_t_to_proj'
1647 
1648  INTEGER :: dim0, handle, my_accelerator, &
1649  my_algorithm
1650  LOGICAL :: smearing, use_sigma_inv_guess
1651  REAL(kind=dp) :: my_inv_eps_factor
1652  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: diag_correction
1653  TYPE(dbcsr_type) :: t_tmp
1654 
1655  CALL timeset(routinen, handle)
1656 
1657  ! make sure that S, sigma and sigma_inv are present for non-orthogonal orbitals
1658  IF (.NOT. orthog_orbs) THEN
1659  IF ((.NOT. PRESENT(s)) .OR. (.NOT. PRESENT(sigma)) .OR. &
1660  (.NOT. PRESENT(sigma_inv)) .OR. (.NOT. PRESENT(nocc_of_domain))) THEN
1661  cpabort("Nonorthogonal orbitals need more input")
1662  END IF
1663  END IF
1664 
1665  my_algorithm = 0
1666  IF (PRESENT(algorithm)) my_algorithm = algorithm
1667 
1668  IF (my_algorithm == 1 .AND. (.NOT. PRESENT(para_env) .OR. .NOT. PRESENT(blacs_env))) &
1669  cpabort("PARA and BLACS env are necessary for cholesky algorithm")
1670 
1671  use_sigma_inv_guess = .false.
1672  IF (PRESENT(use_guess)) THEN
1673  use_sigma_inv_guess = use_guess
1674  END IF
1675 
1676  IF (.NOT. PRESENT(smear)) THEN
1677  smearing = .false.
1678  ELSE
1679  smearing = smear
1680  END IF
1681 
1682  my_accelerator = 1
1683  IF (PRESENT(inverse_accelerator)) my_accelerator = inverse_accelerator
1684 
1685  my_inv_eps_factor = 10.0_dp
1686  IF (PRESENT(inv_eps_factor)) my_inv_eps_factor = inv_eps_factor
1687 
1688  IF (orthog_orbs) THEN
1689 
1690  CALL dbcsr_multiply("N", "T", 1.0_dp, t, t, &
1691  0.0_dp, p, filter_eps=eps_filter)
1692 
1693  ELSE
1694 
1695  CALL dbcsr_create(t_tmp, template=t)
1696 
1697  ! TMP=S.T
1698  CALL dbcsr_multiply("N", "N", 1.0_dp, s, t, 0.0_dp, t_tmp, &
1699  filter_eps=eps_filter)
1700 
1701  ! Sig=tr(T).TMP - get MO overlap
1702  CALL dbcsr_multiply("T", "N", 1.0_dp, t, t_tmp, 0.0_dp, sigma, &
1703  filter_eps=eps_filter)
1704 
1705  !! If smearing ALMO is requested, apply correction of the occupation-rescaling trick
1706  !! (i.e. converting rescaled orbitals into selfish orbitals)
1707  !! (i.e. the diagonal blocks of the rescaled overlap must remain unscaled)
1708  !! Since we have here orthonormal MOs within a fragment, diagonal blocks are identity matrices
1709  !! Therefore, one only need to restore the diagonal to 1
1710  !! RS-WARNING: Assume orthonormal MOs within a fragment
1711  IF (smearing) THEN
1712  CALL dbcsr_get_info(sigma, nfullrows_total=dim0)
1713  ALLOCATE (diag_correction(dim0))
1714  diag_correction = 1.0_dp
1715  CALL dbcsr_set_diag(sigma, diag_correction)
1716  DEALLOCATE (diag_correction)
1717  END IF
1718 
1719  ! invert MO overlap
1720  CALL set_zero_electron_blocks_in_mo_mo_matrix(sigma, nocc_of_domain, 1.0_dp)
1721  SELECT CASE (my_algorithm)
1722  CASE (spd_inversion_ls_taylor)
1723 
1724  CALL invert_taylor( &
1725  matrix_inverse=sigma_inv, &
1726  matrix=sigma, &
1727  use_inv_as_guess=use_sigma_inv_guess, &
1728  threshold=eps_filter*my_inv_eps_factor, &
1729  filter_eps=eps_filter, &
1730  !accelerator_order=my_accelerator, &
1731  !eps_lanczos=eps_lanczos, &
1732  !max_iter_lanczos=max_iter_lanczos, &
1733  silent=.false.)
1734 
1735  CASE (spd_inversion_ls_hotelling)
1736 
1737  CALL invert_hotelling( &
1738  matrix_inverse=sigma_inv, &
1739  matrix=sigma, &
1740  use_inv_as_guess=use_sigma_inv_guess, &
1741  threshold=eps_filter*my_inv_eps_factor, &
1742  filter_eps=eps_filter, &
1743  accelerator_order=my_accelerator, &
1744  eps_lanczos=eps_lanczos, &
1745  max_iter_lanczos=max_iter_lanczos, &
1746  silent=.false.)
1747 
1748  CASE (spd_inversion_dense_cholesky)
1749 
1750  ! invert using cholesky
1751  CALL dbcsr_copy(sigma_inv, sigma)
1752  CALL cp_dbcsr_cholesky_decompose(sigma_inv, &
1753  para_env=para_env, &
1754  blacs_env=blacs_env)
1755  CALL cp_dbcsr_cholesky_invert(sigma_inv, &
1756  para_env=para_env, &
1757  blacs_env=blacs_env, &
1758  upper_to_full=.true.)
1759  CALL dbcsr_filter(sigma_inv, eps_filter)
1760  CASE DEFAULT
1761  cpabort("Illegal MO overalp inversion algorithm")
1762  END SELECT
1763  CALL set_zero_electron_blocks_in_mo_mo_matrix(sigma, nocc_of_domain, 0.0_dp)
1764  CALL set_zero_electron_blocks_in_mo_mo_matrix(sigma_inv, nocc_of_domain, 0.0_dp)
1765 
1766  ! TMP=T.SigInv
1767  CALL dbcsr_multiply("N", "N", 1.0_dp, t, sigma_inv, 0.0_dp, t_tmp, &
1768  filter_eps=eps_filter)
1769 
1770  ! P=TMP.tr(T_blk)
1771  CALL dbcsr_multiply("N", "T", 1.0_dp, t_tmp, t, 0.0_dp, p, &
1772  filter_eps=eps_filter)
1773 
1774  CALL dbcsr_release(t_tmp)
1775 
1776  END IF
1777 
1778  CALL timestop(handle)
1779 
1780  END SUBROUTINE almo_scf_t_to_proj
1781 
1782 ! **************************************************************************************************
1783 !> \brief self-explanatory
1784 !> \param matrix ...
1785 !> \param nocc_of_domain ...
1786 !> \param value ...
1787 !> \param
1788 !> \par History
1789 !> 2016.12 created [Rustam Z Khaliullin]
1790 !> \author Rustam Z Khaliullin
1791 ! **************************************************************************************************
1792  SUBROUTINE set_zero_electron_blocks_in_mo_mo_matrix(matrix, nocc_of_domain, value)
1793 
1794  TYPE(dbcsr_type), INTENT(INOUT) :: matrix
1795  INTEGER, DIMENSION(:), INTENT(IN) :: nocc_of_domain
1796  REAL(kind=dp), INTENT(IN) :: value
1797 
1798  INTEGER :: hold, iblock_col, iblock_row, mynode, &
1799  nblkrows_tot, row
1800  LOGICAL :: found, tr
1801  REAL(kind=dp), DIMENSION(:, :), POINTER :: p_new_block
1802  TYPE(dbcsr_distribution_type) :: dist
1803 
1804  CALL dbcsr_get_info(matrix, distribution=dist)
1805  CALL dbcsr_distribution_get(dist, mynode=mynode)
1806  !mynode = dbcsr_mp_mynode(dbcsr_distribution_mp(dbcsr_distribution(matrix)))
1807  CALL dbcsr_work_create(matrix, work_mutable=.true.)
1808 
1809  nblkrows_tot = dbcsr_nblkrows_total(matrix)
1810 
1811  DO row = 1, nblkrows_tot
1812  IF (nocc_of_domain(row) == 0) THEN
1813  tr = .false.
1814  iblock_row = row
1815  iblock_col = row
1816  CALL dbcsr_get_stored_coordinates(matrix, iblock_row, iblock_col, hold)
1817  IF (hold .EQ. mynode) THEN
1818  NULLIFY (p_new_block)
1819  CALL dbcsr_get_block_p(matrix, iblock_row, iblock_col, p_new_block, found)
1820  IF (found) THEN
1821  p_new_block(1, 1) = value
1822  ELSE
1823  CALL dbcsr_reserve_block2d(matrix, iblock_row, iblock_col, p_new_block)
1824  cpassert(ASSOCIATED(p_new_block))
1825  p_new_block(1, 1) = value
1826  END IF
1827  END IF ! mynode
1828  END IF !zero-electron block
1829  END DO
1830 
1831  CALL dbcsr_finalize(matrix)
1832 
1833  END SUBROUTINE set_zero_electron_blocks_in_mo_mo_matrix
1834 
1835 ! **************************************************************************************************
1836 !> \brief applies projector to the orbitals
1837 !> |psi_out> = P |psi_in> OR |psi_out> = (1-P) |psi_in>,
1838 !> where P = |psi_proj> (<psi_proj|psi_roj>)^{-1} <psi_proj|
1839 !> \param psi_in ...
1840 !> \param psi_out ...
1841 !> \param psi_projector ...
1842 !> \param metric ...
1843 !> \param project_out ...
1844 !> \param psi_projector_orthogonal ...
1845 !> \param proj_in_template ...
1846 !> \param eps_filter ...
1847 !> \param sig_inv_projector ...
1848 !> \param sig_inv_template ...
1849 !> \par History
1850 !> 2011.10 created [Rustam Z Khaliullin]
1851 !> \author Rustam Z Khaliullin
1852 ! **************************************************************************************************
1853  SUBROUTINE apply_projector(psi_in, psi_out, psi_projector, metric, project_out, &
1854  psi_projector_orthogonal, proj_in_template, eps_filter, sig_inv_projector, &
1855  sig_inv_template)
1856 
1857  TYPE(dbcsr_type), INTENT(IN) :: psi_in
1858  TYPE(dbcsr_type), INTENT(INOUT) :: psi_out
1859  TYPE(dbcsr_type), INTENT(IN) :: psi_projector, metric
1860  LOGICAL, INTENT(IN) :: project_out, psi_projector_orthogonal
1861  TYPE(dbcsr_type), INTENT(IN) :: proj_in_template
1862  REAL(kind=dp), INTENT(IN) :: eps_filter
1863  TYPE(dbcsr_type), INTENT(IN), OPTIONAL :: sig_inv_projector, sig_inv_template
1864 
1865  CHARACTER(LEN=*), PARAMETER :: routinen = 'apply_projector'
1866 
1867  INTEGER :: handle
1868  TYPE(dbcsr_type) :: tmp_no, tmp_ov, tmp_ov2, tmp_sig, &
1869  tmp_sig_inv
1870 
1871  CALL timeset(routinen, handle)
1872 
1873  ! =S*PSI_proj
1874  CALL dbcsr_create(tmp_no, template=psi_projector)
1875  CALL dbcsr_multiply("N", "N", 1.0_dp, &
1876  metric, psi_projector, &
1877  0.0_dp, tmp_no, &
1878  filter_eps=eps_filter)
1879 
1880  ! =tr(S.PSI_proj)*PSI_in
1881  CALL dbcsr_create(tmp_ov, template=proj_in_template)
1882  CALL dbcsr_multiply("T", "N", 1.0_dp, &
1883  tmp_no, psi_in, &
1884  0.0_dp, tmp_ov, &
1885  filter_eps=eps_filter)
1886 
1887  IF (.NOT. psi_projector_orthogonal) THEN
1888  ! =SigInv_proj*Sigma_OV
1889  CALL dbcsr_create(tmp_ov2, &
1890  template=proj_in_template)
1891  IF (PRESENT(sig_inv_projector)) THEN
1892  CALL dbcsr_create(tmp_sig_inv, &
1893  template=sig_inv_projector)
1894  CALL dbcsr_copy(tmp_sig_inv, sig_inv_projector)
1895  ELSE
1896  IF (.NOT. PRESENT(sig_inv_template)) THEN
1897  cpabort("PROGRAMMING ERROR: provide either template or sig_inv")
1898  END IF
1899  ! compute inverse overlap of the projector orbitals
1900  CALL dbcsr_create(tmp_sig, &
1901  template=sig_inv_template, &
1902  matrix_type=dbcsr_type_no_symmetry)
1903  CALL dbcsr_multiply("T", "N", 1.0_dp, &
1904  psi_projector, tmp_no, 0.0_dp, tmp_sig, &
1905  filter_eps=eps_filter)
1906  CALL dbcsr_create(tmp_sig_inv, &
1907  template=sig_inv_template, &
1908  matrix_type=dbcsr_type_no_symmetry)
1909  CALL invert_hotelling(tmp_sig_inv, tmp_sig, &
1910  threshold=eps_filter)
1911  CALL dbcsr_release(tmp_sig)
1912  END IF
1913  CALL dbcsr_multiply("N", "N", 1.0_dp, &
1914  tmp_sig_inv, tmp_ov, 0.0_dp, tmp_ov2, &
1915  filter_eps=eps_filter)
1916  CALL dbcsr_release(tmp_sig_inv)
1917  CALL dbcsr_copy(tmp_ov, tmp_ov2)
1918  CALL dbcsr_release(tmp_ov2)
1919  END IF
1920  CALL dbcsr_release(tmp_no)
1921 
1922  ! =PSI_proj*TMP_OV
1923  CALL dbcsr_multiply("N", "N", 1.0_dp, &
1924  psi_projector, tmp_ov, 0.0_dp, psi_out, &
1925  filter_eps=eps_filter)
1926  CALL dbcsr_release(tmp_ov)
1927 
1928  ! V_out=V_in-V_out
1929  IF (project_out) THEN
1930  CALL dbcsr_add(psi_out, psi_in, -1.0_dp, +1.0_dp)
1931  END IF
1932 
1933  CALL timestop(handle)
1934 
1935  END SUBROUTINE apply_projector
1936 
1937 !! **************************************************************************************************
1938 !!> \brief projects the occupied space out from the provided orbitals
1939 !!> \par History
1940 !!> 2011.07 created [Rustam Z Khaliullin]
1941 !!> \author Rustam Z Khaliullin
1942 !! **************************************************************************************************
1943 ! SUBROUTINE almo_scf_p_out_from_v(v_in,v_out,ov_template,ispin,almo_scf_env)
1944 !
1945 ! TYPE(dbcsr_type), INTENT(IN) :: v_in, ov_template
1946 ! TYPE(dbcsr_type), INTENT(INOUT) :: v_out
1947 ! INTEGER, INTENT(IN) :: ispin
1948 ! TYPE(almo_scf_env_type), INTENT(INOUT) :: almo_scf_env
1949 !
1950 ! CHARACTER(LEN=*), PARAMETER :: &
1951 ! routineN = 'almo_scf_p_out_from_v', &
1952 ! routineP = moduleN//':'//routineN
1953 !
1954 ! TYPE(dbcsr_type) :: tmp_on, tmp_ov, tmp_ov2
1955 ! INTEGER :: handle
1956 ! LOGICAL :: failure
1957 !
1958 ! CALL timeset(routineN,handle)
1959 !
1960 ! ! =tr(T_blk)*S
1961 ! CALL dbcsr_init(tmp_on)
1962 ! CALL dbcsr_create(tmp_on,&
1963 ! template=almo_scf_env%matrix_t_tr(ispin))
1964 ! CALL dbcsr_multiply("T","N",1.0_dp,&
1965 ! almo_scf_env%matrix_t_blk(ispin),&
1966 ! almo_scf_env%matrix_s(1),&
1967 ! 0.0_dp,tmp_on,&
1968 ! filter_eps=almo_scf_env%eps_filter)
1969 !
1970 ! ! =tr(T_blk).S*V_in
1971 ! CALL dbcsr_init(tmp_ov)
1972 ! CALL dbcsr_create(tmp_ov,template=ov_template)
1973 ! CALL dbcsr_multiply("N","N",1.0_dp,&
1974 ! tmp_on,v_in,0.0_dp,tmp_ov,&
1975 ! filter_eps=almo_scf_env%eps_filter)
1976 ! CALL dbcsr_release(tmp_on)
1977 !
1978 ! ! =SigmaInv*Sigma_OV
1979 ! CALL dbcsr_init(tmp_ov2)
1980 ! CALL dbcsr_create(tmp_ov2,template=ov_template)
1981 ! CALL dbcsr_multiply("N","N",1.0_dp,&
1982 ! almo_scf_env%matrix_sigma_inv(ispin),&
1983 ! tmp_ov,0.0_dp,tmp_ov2,&
1984 ! filter_eps=almo_scf_env%eps_filter)
1985 ! CALL dbcsr_release(tmp_ov)
1986 !
1987 ! ! =T_blk*SigmaInv.Sigma_OV
1988 ! CALL dbcsr_multiply("N","N",1.0_dp,&
1989 ! almo_scf_env%matrix_t_blk(ispin),&
1990 ! tmp_ov2,0.0_dp,v_out,&
1991 ! filter_eps=almo_scf_env%eps_filter)
1992 ! CALL dbcsr_release(tmp_ov2)
1993 !
1994 ! ! V_out=V_in-V_out=
1995 ! CALL dbcsr_add(v_out,v_in,-1.0_dp,+1.0_dp)
1996 !
1997 ! CALL timestop(handle)
1998 !
1999 ! END SUBROUTINE almo_scf_p_out_from_v
2000 
2001 ! **************************************************************************************************
2002 !> \brief computes a unitary matrix from an arbitrary "generator" matrix
2003 !> U = ( 1 - X + tr(X) ) ( 1 + X - tr(X) )^(-1)
2004 !> \param X ...
2005 !> \param U ...
2006 !> \param eps_filter ...
2007 !> \par History
2008 !> 2011.08 created [Rustam Z Khaliullin]
2009 !> \author Rustam Z Khaliullin
2010 ! **************************************************************************************************
2011  SUBROUTINE generator_to_unitary(X, U, eps_filter)
2012 
2013  TYPE(dbcsr_type), INTENT(IN) :: x
2014  TYPE(dbcsr_type), INTENT(INOUT) :: u
2015  REAL(kind=dp), INTENT(IN) :: eps_filter
2016 
2017  CHARACTER(LEN=*), PARAMETER :: routinen = 'generator_to_unitary'
2018 
2019  INTEGER :: handle, unit_nr
2020  LOGICAL :: safe_mode
2021  REAL(kind=dp) :: frob_matrix, frob_matrix_base
2022  TYPE(cp_logger_type), POINTER :: logger
2023  TYPE(dbcsr_type) :: delta, t1, t2, tmp1
2024 
2025  CALL timeset(routinen, handle)
2026 
2027  safe_mode = .true.
2028 
2029  ! get a useful output_unit
2030  logger => cp_get_default_logger()
2031  IF (logger%para_env%is_source()) THEN
2032  unit_nr = cp_logger_get_default_unit_nr(logger, local=.true.)
2033  ELSE
2034  unit_nr = -1
2035  END IF
2036 
2037  CALL dbcsr_create(t1, template=x, &
2038  matrix_type=dbcsr_type_no_symmetry)
2039  CALL dbcsr_create(t2, template=x, &
2040  matrix_type=dbcsr_type_no_symmetry)
2041 
2042  ! create antisymmetric Delta = -X + tr(X)
2043  CALL dbcsr_create(delta, template=x, &
2044  matrix_type=dbcsr_type_no_symmetry)
2045  CALL dbcsr_transposed(delta, x)
2046 ! check that transposed is added correctly
2047  CALL dbcsr_add(delta, x, 1.0_dp, -1.0_dp)
2048 
2049  ! compute (1 - Delta)^(-1)
2050  CALL dbcsr_add_on_diag(t1, 1.0_dp)
2051  CALL dbcsr_add(t1, delta, 1.0_dp, -1.0_dp)
2052  CALL invert_hotelling(t2, t1, threshold=eps_filter)
2053 
2054  IF (safe_mode) THEN
2055 
2056  CALL dbcsr_create(tmp1, template=x, &
2057  matrix_type=dbcsr_type_no_symmetry)
2058  CALL dbcsr_multiply("N", "N", 1.0_dp, t2, t1, 0.0_dp, tmp1, &
2059  filter_eps=eps_filter)
2060  frob_matrix_base = dbcsr_frobenius_norm(tmp1)
2061  CALL dbcsr_add_on_diag(tmp1, -1.0_dp)
2062  frob_matrix = dbcsr_frobenius_norm(tmp1)
2063  IF (unit_nr > 0) WRITE (unit_nr, *) "Error for (inv(A)*A-I)", frob_matrix/frob_matrix_base
2064  CALL dbcsr_release(tmp1)
2065  END IF
2066 
2067  CALL dbcsr_multiply("N", "N", 1.0_dp, delta, t2, 0.0_dp, u, &
2068  filter_eps=eps_filter)
2069  CALL dbcsr_add(u, t2, 1.0_dp, 1.0_dp)
2070 
2071  IF (safe_mode) THEN
2072 
2073  CALL dbcsr_create(tmp1, template=x, &
2074  matrix_type=dbcsr_type_no_symmetry)
2075  CALL dbcsr_multiply("T", "N", 1.0_dp, u, u, 0.0_dp, tmp1, &
2076  filter_eps=eps_filter)
2077  frob_matrix_base = dbcsr_frobenius_norm(tmp1)
2078  CALL dbcsr_add_on_diag(tmp1, -1.0_dp)
2079  frob_matrix = dbcsr_frobenius_norm(tmp1)
2080  IF (unit_nr > 0) WRITE (unit_nr, *) "Error for (trn(U)*U-I)", frob_matrix/frob_matrix_base
2081  CALL dbcsr_release(tmp1)
2082  END IF
2083 
2084  CALL timestop(handle)
2085 
2086  END SUBROUTINE generator_to_unitary
2087 
2088 ! **************************************************************************************************
2089 !> \brief Parallel code for domain specific operations (my_action)
2090 !> 0. out = op1 * in
2091 !> 1. out = in - op2 * op1 * in
2092 !> \param matrix_in ...
2093 !> \param matrix_out ...
2094 !> \param operator1 ...
2095 !> \param operator2 ...
2096 !> \param dpattern ...
2097 !> \param map ...
2098 !> \param node_of_domain ...
2099 !> \param my_action ...
2100 !> \param filter_eps ...
2101 !> \param matrix_trimmer ...
2102 !> \param use_trimmer ...
2103 !> \par History
2104 !> 2013.01 created [Rustam Z. Khaliullin]
2105 !> \author Rustam Z. Khaliullin
2106 ! **************************************************************************************************
2107  SUBROUTINE apply_domain_operators(matrix_in, matrix_out, operator1, operator2, &
2108  dpattern, map, node_of_domain, my_action, filter_eps, matrix_trimmer, use_trimmer)
2109 
2110  TYPE(dbcsr_type), INTENT(IN) :: matrix_in
2111  TYPE(dbcsr_type), INTENT(INOUT) :: matrix_out
2112  TYPE(domain_submatrix_type), DIMENSION(:), &
2113  INTENT(IN) :: operator1
2114  TYPE(domain_submatrix_type), DIMENSION(:), &
2115  INTENT(IN), OPTIONAL :: operator2
2116  TYPE(dbcsr_type), INTENT(IN) :: dpattern
2117  TYPE(domain_map_type), INTENT(IN) :: map
2118  INTEGER, DIMENSION(:), INTENT(IN) :: node_of_domain
2119  INTEGER, INTENT(IN) :: my_action
2120  REAL(kind=dp) :: filter_eps
2121  TYPE(dbcsr_type), INTENT(IN), OPTIONAL :: matrix_trimmer
2122  LOGICAL, INTENT(IN), OPTIONAL :: use_trimmer
2123 
2124  CHARACTER(len=*), PARAMETER :: routinen = 'apply_domain_operators'
2125 
2126  INTEGER :: handle, ndomains
2127  LOGICAL :: matrix_trimmer_required, my_use_trimmer, &
2128  operator2_required
2129  TYPE(domain_submatrix_type), ALLOCATABLE, &
2130  DIMENSION(:) :: subm_in, subm_out, subm_temp
2131 
2132  CALL timeset(routinen, handle)
2133 
2134  my_use_trimmer = .false.
2135  IF (PRESENT(use_trimmer)) THEN
2136  my_use_trimmer = use_trimmer
2137  END IF
2138 
2139  operator2_required = .false.
2140  matrix_trimmer_required = .false.
2141 
2142  IF (my_action .EQ. 1) operator2_required = .true.
2143 
2144  IF (my_use_trimmer) THEN
2145  matrix_trimmer_required = .true.
2146  cpabort("TRIMMED PROJECTOR DISABLED!")
2147  END IF
2148 
2149  IF (.NOT. PRESENT(operator2) .AND. operator2_required) THEN
2150  cpabort("SECOND OPERATOR IS REQUIRED")
2151  END IF
2152  IF (.NOT. PRESENT(matrix_trimmer) .AND. matrix_trimmer_required) THEN
2153  cpabort("TRIMMER MATRIX IS REQUIRED")
2154  END IF
2155 
2156  ndomains = dbcsr_nblkcols_total(dpattern)
2157 
2158  ALLOCATE (subm_in(ndomains))
2159  ALLOCATE (subm_temp(ndomains))
2160  ALLOCATE (subm_out(ndomains))
2161  !!!TRIM ALLOCATE(subm_trimmer(ndomains))
2162  CALL init_submatrices(subm_in)
2163  CALL init_submatrices(subm_temp)
2164  CALL init_submatrices(subm_out)
2165 
2166  CALL construct_submatrices(matrix_in, subm_in, &
2167  dpattern, map, node_of_domain, select_row)
2168 
2169  !!!TRIM IF (matrix_trimmer_required) THEN
2170  !!!TRIM CALL construct_submatrices(matrix_trimmer,subm_trimmer,&
2171  !!!TRIM dpattern,map,node_of_domain,select_row)
2172  !!!TRIM ENDIF
2173 
2174  IF (my_action .EQ. 0) THEN
2175  ! for example, apply preconditioner
2176  CALL multiply_submatrices('N', 'N', 1.0_dp, operator1, &
2177  subm_in, 0.0_dp, subm_out)
2178  ELSE IF (my_action .EQ. 1) THEN
2179  ! use for projectors
2180  CALL copy_submatrices(subm_in, subm_out, .true.)
2181  CALL multiply_submatrices('N', 'N', 1.0_dp, operator1, &
2182  subm_in, 0.0_dp, subm_temp)
2183  CALL multiply_submatrices('N', 'N', -1.0_dp, operator2, &
2184  subm_temp, 1.0_dp, subm_out)
2185  ELSE
2186  cpabort("ILLEGAL ACTION")
2187  END IF
2188 
2189  CALL construct_dbcsr_from_submatrices(matrix_out, subm_out, dpattern)
2190  CALL dbcsr_filter(matrix_out, filter_eps)
2191 
2192  CALL release_submatrices(subm_out)
2193  CALL release_submatrices(subm_temp)
2194  CALL release_submatrices(subm_in)
2195 
2196  DEALLOCATE (subm_out)
2197  DEALLOCATE (subm_temp)
2198  DEALLOCATE (subm_in)
2199 
2200  CALL timestop(handle)
2201 
2202  END SUBROUTINE apply_domain_operators
2203 
2204 ! **************************************************************************************************
2205 !> \brief Constructs preconditioners for each domain
2206 !> -1. projected preconditioner
2207 !> 0. simple preconditioner
2208 !> \param matrix_main ...
2209 !> \param subm_s_inv ...
2210 !> \param subm_s_inv_half ...
2211 !> \param subm_s_half ...
2212 !> \param subm_r_down ...
2213 !> \param matrix_trimmer ...
2214 !> \param dpattern ...
2215 !> \param map ...
2216 !> \param node_of_domain ...
2217 !> \param preconditioner ...
2218 !> \param bad_modes_projector_down ...
2219 !> \param use_trimmer ...
2220 !> \param eps_zero_eigenvalues ...
2221 !> \param my_action ...
2222 !> \param skip_inversion ...
2223 !> \par History
2224 !> 2013.01 created [Rustam Z. Khaliullin]
2225 !> \author Rustam Z. Khaliullin
2226 ! **************************************************************************************************
2227  SUBROUTINE construct_domain_preconditioner(matrix_main, subm_s_inv, &
2228  subm_s_inv_half, subm_s_half, &
2229  subm_r_down, matrix_trimmer, &
2230  dpattern, map, node_of_domain, &
2231  preconditioner, &
2232  bad_modes_projector_down, &
2233  use_trimmer, &
2234  eps_zero_eigenvalues, &
2235  my_action, skip_inversion)
2236 
2237  TYPE(dbcsr_type), INTENT(INOUT) :: matrix_main
2238  TYPE(domain_submatrix_type), DIMENSION(:), &
2239  INTENT(IN), OPTIONAL :: subm_s_inv, subm_s_inv_half, &
2240  subm_s_half, subm_r_down
2241  TYPE(dbcsr_type), INTENT(IN), OPTIONAL :: matrix_trimmer
2242  TYPE(dbcsr_type), INTENT(IN) :: dpattern
2243  TYPE(domain_map_type), INTENT(IN) :: map
2244  INTEGER, DIMENSION(:), INTENT(IN) :: node_of_domain
2245  TYPE(domain_submatrix_type), DIMENSION(:), &
2246  INTENT(INOUT) :: preconditioner
2247  TYPE(domain_submatrix_type), DIMENSION(:), &
2248  INTENT(INOUT), OPTIONAL :: bad_modes_projector_down
2249  LOGICAL, INTENT(IN), OPTIONAL :: use_trimmer
2250  REAL(kind=dp), INTENT(IN), OPTIONAL :: eps_zero_eigenvalues
2251  INTEGER, INTENT(IN) :: my_action
2252  LOGICAL, INTENT(IN), OPTIONAL :: skip_inversion
2253 
2254  CHARACTER(len=*), PARAMETER :: routinen = 'construct_domain_preconditioner'
2255 
2256  INTEGER :: handle, idomain, index1_end, &
2257  index1_start, n_domain_mos, naos, &
2258  nblkrows_tot, ndomains, neighbor, row
2259  INTEGER, DIMENSION(:), POINTER :: nmos
2260  LOGICAL :: eps_zero_eigenvalues_required, matrix_r_required, matrix_s_half_required, &
2261  matrix_s_inv_half_required, matrix_s_inv_required, matrix_trimmer_required, &
2262  my_skip_inversion, my_use_trimmer
2263  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: minv, proj_array
2264  TYPE(domain_submatrix_type), ALLOCATABLE, &
2265  DIMENSION(:) :: subm_main, subm_tmp, subm_tmp2
2266 
2267  CALL timeset(routinen, handle)
2268 
2269  my_use_trimmer = .false.
2270  IF (PRESENT(use_trimmer)) THEN
2271  my_use_trimmer = use_trimmer
2272  END IF
2273 
2274  my_skip_inversion = .false.
2275  IF (PRESENT(skip_inversion)) THEN
2276  my_skip_inversion = skip_inversion
2277  END IF
2278 
2279  matrix_s_inv_half_required = .false.
2280  matrix_s_half_required = .false.
2281  eps_zero_eigenvalues_required = .false.
2282  matrix_s_inv_required = .false.
2283  matrix_trimmer_required = .false.
2284  matrix_r_required = .false.
2285 
2286  IF (my_action .EQ. -1) matrix_s_inv_required = .true.
2287  IF (my_action .EQ. -1) matrix_r_required = .true.
2288  IF (my_use_trimmer) THEN
2289  matrix_trimmer_required = .true.
2290  cpabort("TRIMMED PRECONDITIONER DISABLED!")
2291  END IF
2292  ! tie the following optional arguments together to prevent bad calls
2293  IF (PRESENT(bad_modes_projector_down)) THEN
2294  matrix_s_inv_half_required = .true.
2295  matrix_s_half_required = .true.
2296  eps_zero_eigenvalues_required = .true.
2297  END IF
2298 
2299  ! check if all required optional arguments are provided
2300  IF (.NOT. PRESENT(subm_s_inv_half) .AND. matrix_s_inv_half_required) THEN
2301  cpabort("S_inv_half SUBMATRICES ARE REQUIRED")
2302  END IF
2303  IF (.NOT. PRESENT(subm_s_half) .AND. matrix_s_half_required) THEN
2304  cpabort("S_half SUBMATRICES ARE REQUIRED")
2305  END IF
2306  IF (.NOT. PRESENT(eps_zero_eigenvalues) .AND. eps_zero_eigenvalues_required) THEN
2307  cpabort("EPS_ZERO_EIGENVALUES IS REQUIRED")
2308  END IF
2309  IF (.NOT. PRESENT(subm_s_inv) .AND. matrix_s_inv_required) THEN
2310  cpabort("S_inv SUBMATRICES ARE REQUIRED")
2311  END IF
2312  IF (.NOT. PRESENT(subm_r_down) .AND. matrix_r_required) THEN
2313  cpabort("R SUBMATRICES ARE REQUIRED")
2314  END IF
2315  IF (.NOT. PRESENT(matrix_trimmer) .AND. matrix_trimmer_required) THEN
2316  cpabort("TRIMMER MATRIX IS REQUIRED")
2317  END IF
2318 
2319  CALL dbcsr_get_info(dpattern, &
2320  nblkcols_total=ndomains, &
2321  nblkrows_total=nblkrows_tot, &
2322  col_blk_size=nmos)
2323 
2324  ALLOCATE (subm_main(ndomains))
2325  CALL init_submatrices(subm_main)
2326  !!!TRIM ALLOCATE(subm_trimmer(ndomains))
2327 
2328  CALL construct_submatrices(matrix_main, subm_main, &
2329  dpattern, map, node_of_domain, select_row_col)
2330 
2331  !!!TRIM IF (matrix_trimmer_required) THEN
2332  !!!TRIM CALL construct_submatrices(matrix_trimmer,subm_trimmer,&
2333  !!!TRIM dpattern,map,node_of_domain,select_row)
2334  !!!TRIM ENDIF
2335 
2336  IF (my_action .EQ. -1) THEN
2337  ! project out the local occupied space
2338  !tmp=MATMUL(subm_r(idomain)%mdata,Minv)
2339  !Minv=MATMUL(tmp,subm_main(idomain)%mdata)
2340  !subm_main(idomain)%mdata=subm_main(idomain)%mdata-&
2341  ! Minv-TRANSPOSE(Minv)+MATMUL(Minv,TRANSPOSE(tmp))
2342  ALLOCATE (subm_tmp(ndomains))
2343  ALLOCATE (subm_tmp2(ndomains))
2344  CALL init_submatrices(subm_tmp)
2345  CALL init_submatrices(subm_tmp2)
2346  CALL multiply_submatrices('N', 'N', 1.0_dp, subm_r_down, &
2347  subm_s_inv, 0.0_dp, subm_tmp)
2348  CALL multiply_submatrices('N', 'N', 1.0_dp, subm_tmp, &
2349  subm_main, 0.0_dp, subm_tmp2)
2350  CALL add_submatrices(1.0_dp, subm_main, -1.0_dp, subm_tmp2, 'N')
2351  CALL add_submatrices(1.0_dp, subm_main, -1.0_dp, subm_tmp2, 'T')
2352  CALL multiply_submatrices('N', 'T', 1.0_dp, subm_tmp2, &
2353  subm_tmp, 1.0_dp, subm_main)
2354  CALL release_submatrices(subm_tmp)
2355  CALL release_submatrices(subm_tmp2)
2356  DEALLOCATE (subm_tmp2)
2357  DEALLOCATE (subm_tmp)
2358  END IF
2359 
2360  IF (my_skip_inversion) THEN
2361  CALL copy_submatrices(subm_main, preconditioner, .true.)
2362  ELSE
2363  ! loop over domains - perform inversion
2364  DO idomain = 1, ndomains
2365 
2366  ! check if the submatrix exists
2367  IF (subm_main(idomain)%domain .GT. 0) THEN
2368 
2369  ! find sizes of MO submatrices
2370  IF (idomain .EQ. 1) THEN
2371  index1_start = 1
2372  ELSE
2373  index1_start = map%index1(idomain - 1)
2374  END IF
2375  index1_end = map%index1(idomain) - 1
2376 
2377  n_domain_mos = 0
2378  DO row = index1_start, index1_end
2379  neighbor = map%pairs(row, 1)
2380  n_domain_mos = n_domain_mos + nmos(neighbor)
2381  END DO
2382 
2383  naos = subm_main(idomain)%nrows
2384  !WRITE(*,*) "Domain, mo_self_and_neig, ao_domain: ", idomain, n_domain_mos, naos
2385 
2386  ALLOCATE (minv(naos, naos))
2387 
2388  !!!TRIM IF (my_use_trimmer) THEN
2389  !!!TRIM ! THIS IS SUPER EXPENSIVE (ELIMINATE)
2390  !!!TRIM ! trim the main matrix before inverting
2391  !!!TRIM ! assume that the trimmer columns are different (i.e. the main matrix is different for each MO)
2392  !!!TRIM allocate(tmp(naos,nmos(idomain)))
2393  !!!TRIM DO ii=1, nmos(idomain)
2394  !!!TRIM ! transform the main matrix using the trimmer for the current MO
2395  !!!TRIM DO jj=1, naos
2396  !!!TRIM DO kk=1, naos
2397  !!!TRIM Mstore(jj,kk)=sumb_main(idomain)%mdata(jj,kk)*&
2398  !!!TRIM subm_trimmer(idomain)%mdata(jj,ii)*&
2399  !!!TRIM subm_trimmer(idomain)%mdata(kk,ii)
2400  !!!TRIM ENDDO
2401  !!!TRIM ENDDO
2402  !!!TRIM ! invert the main matrix (exclude some eigenvalues, shift some)
2403  !!!TRIM CALL pseudo_invert_matrix(A=Mstore,Ainv=Minv,N=naos,method=1,&
2404  !!!TRIM !range1_thr=1.0E-9_dp,range2_thr=1.0E-9_dp,&
2405  !!!TRIM shift=1.0E-5_dp,&
2406  !!!TRIM range1=nmos(idomain),range2=nmos(idomain),&
2407  !!!TRIM
2408  !!!TRIM ! apply the inverted matrix
2409  !!!TRIM ! RZK-warning this is only possible when the preconditioner is applied
2410  !!!TRIM tmp(:,ii)=MATMUL(Minv,subm_in(idomain)%mdata(:,ii))
2411  !!!TRIM ENDDO
2412  !!!TRIM subm_out=MATMUL(tmp,sigma)
2413  !!!TRIM deallocate(tmp)
2414  !!!TRIM ELSE
2415 
2416  IF (PRESENT(bad_modes_projector_down)) THEN
2417  ALLOCATE (proj_array(naos, naos))
2418  CALL pseudo_invert_matrix(a=subm_main(idomain)%mdata, ainv=minv, n=naos, method=1, &
2419  range1=nmos(idomain), range2=n_domain_mos, &
2420  range1_thr=eps_zero_eigenvalues, &
2421  bad_modes_projector_down=proj_array, &
2422  s_inv_half=subm_s_inv_half(idomain)%mdata, &
2423  s_half=subm_s_half(idomain)%mdata &
2424  )
2425  ELSE
2426  CALL pseudo_invert_matrix(a=subm_main(idomain)%mdata, ainv=minv, n=naos, method=1, &
2427  range1=nmos(idomain), range2=n_domain_mos)
2428  END IF
2429  !!!TRIM ENDIF
2430 
2431  CALL copy_submatrices(subm_main(idomain), preconditioner(idomain), .false.)
2432  CALL copy_submatrix_data(minv, preconditioner(idomain))
2433  DEALLOCATE (minv)
2434 
2435  IF (PRESENT(bad_modes_projector_down)) THEN
2436  CALL copy_submatrices(subm_main(idomain), bad_modes_projector_down(idomain), .false.)
2437  CALL copy_submatrix_data(proj_array, bad_modes_projector_down(idomain))
2438  DEALLOCATE (proj_array)
2439  END IF
2440 
2441  END IF ! submatrix for the domain exists
2442 
2443  END DO ! loop over domains
2444 
2445  END IF
2446 
2447  CALL release_submatrices(subm_main)
2448  DEALLOCATE (subm_main)
2449  !DEALLOCATE(subm_s)
2450  !DEALLOCATE(subm_r)
2451 
2452  !IF (matrix_r_required) THEN
2453  ! CALL dbcsr_release(m_tmp_no_1)
2454  ! CALL dbcsr_release(m_tmp_no_2)
2455  ! CALL dbcsr_release(matrix_r)
2456  !ENDIF
2457 
2458  CALL timestop(handle)
2459 
2460  END SUBROUTINE construct_domain_preconditioner
2461 
2462 ! **************************************************************************************************
2463 !> \brief Constructs S^(+1/2) and S^(-1/2) submatrices for each domain
2464 !> \param matrix_s ...
2465 !> \param subm_s_sqrt ...
2466 !> \param subm_s_sqrt_inv ...
2467 !> \param dpattern ...
2468 !> \param map ...
2469 !> \param node_of_domain ...
2470 !> \par History
2471 !> 2013.03 created [Rustam Z. Khaliullin]
2472 !> \author Rustam Z. Khaliullin
2473 ! **************************************************************************************************
2474  SUBROUTINE construct_domain_s_sqrt(matrix_s, subm_s_sqrt, subm_s_sqrt_inv, &
2475  dpattern, map, node_of_domain)
2476 
2477  TYPE(dbcsr_type), INTENT(IN) :: matrix_s
2478  TYPE(domain_submatrix_type), DIMENSION(:), &
2479  INTENT(INOUT) :: subm_s_sqrt, subm_s_sqrt_inv
2480  TYPE(dbcsr_type), INTENT(IN) :: dpattern
2481  TYPE(domain_map_type), INTENT(IN) :: map
2482  INTEGER, DIMENSION(:), INTENT(IN) :: node_of_domain
2483 
2484  CHARACTER(len=*), PARAMETER :: routinen = 'construct_domain_s_sqrt'
2485 
2486  INTEGER :: handle, idomain, naos, ndomains
2487  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: ssqrt, ssqrtinv
2488  TYPE(domain_submatrix_type), ALLOCATABLE, &
2489  DIMENSION(:) :: subm_s
2490 
2491  CALL timeset(routinen, handle)
2492 
2493  ndomains = dbcsr_nblkcols_total(dpattern)
2494  cpassert(SIZE(subm_s_sqrt) .EQ. ndomains)
2495  cpassert(SIZE(subm_s_sqrt_inv) .EQ. ndomains)
2496  ALLOCATE (subm_s(ndomains))
2497  CALL init_submatrices(subm_s)
2498 
2499  CALL construct_submatrices(matrix_s, subm_s, &
2500  dpattern, map, node_of_domain, select_row_col)
2501 
2502  ! loop over domains - perform inversion
2503  DO idomain = 1, ndomains
2504 
2505  ! check if the submatrix exists
2506  IF (subm_s(idomain)%domain .GT. 0) THEN
2507 
2508  naos = subm_s(idomain)%nrows
2509 
2510  ALLOCATE (ssqrt(naos, naos))
2511  ALLOCATE (ssqrtinv(naos, naos))
2512 
2513  CALL matrix_sqrt(a=subm_s(idomain)%mdata, asqrt=ssqrt, asqrtinv=ssqrtinv, &
2514  n=naos)
2515 
2516  CALL copy_submatrices(subm_s(idomain), subm_s_sqrt(idomain), .false.)
2517  CALL copy_submatrix_data(ssqrt, subm_s_sqrt(idomain))
2518 
2519  CALL copy_submatrices(subm_s(idomain), subm_s_sqrt_inv(idomain), .false.)
2520  CALL copy_submatrix_data(ssqrtinv, subm_s_sqrt_inv(idomain))
2521 
2522  DEALLOCATE (ssqrtinv)
2523  DEALLOCATE (ssqrt)
2524 
2525  END IF ! submatrix for the domain exists
2526 
2527  END DO ! loop over domains
2528 
2529  CALL release_submatrices(subm_s)
2530  DEALLOCATE (subm_s)
2531 
2532  CALL timestop(handle)
2533 
2534  END SUBROUTINE construct_domain_s_sqrt
2535 
2536 ! **************************************************************************************************
2537 !> \brief Constructs S_inv block for each domain
2538 !> \param matrix_s ...
2539 !> \param subm_s_inv ...
2540 !> \param dpattern ...
2541 !> \param map ...
2542 !> \param node_of_domain ...
2543 !> \par History
2544 !> 2013.02 created [Rustam Z. Khaliullin]
2545 !> \author Rustam Z. Khaliullin
2546 ! **************************************************************************************************
2547  SUBROUTINE construct_domain_s_inv(matrix_s, subm_s_inv, dpattern, map, &
2548  node_of_domain)
2549  TYPE(dbcsr_type), INTENT(IN) :: matrix_s
2550  TYPE(domain_submatrix_type), DIMENSION(:), &
2551  INTENT(INOUT) :: subm_s_inv
2552  TYPE(dbcsr_type), INTENT(IN) :: dpattern
2553  TYPE(domain_map_type), INTENT(IN) :: map
2554  INTEGER, DIMENSION(:), INTENT(IN) :: node_of_domain
2555 
2556  CHARACTER(len=*), PARAMETER :: routinen = 'construct_domain_s_inv'
2557 
2558  INTEGER :: handle, idomain, naos, ndomains
2559  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: sinv
2560  TYPE(domain_submatrix_type), ALLOCATABLE, &
2561  DIMENSION(:) :: subm_s
2562 
2563  CALL timeset(routinen, handle)
2564 
2565  ndomains = dbcsr_nblkcols_total(dpattern)
2566 
2567  cpassert(SIZE(subm_s_inv) .EQ. ndomains)
2568  ALLOCATE (subm_s(ndomains))
2569  CALL init_submatrices(subm_s)
2570 
2571  CALL construct_submatrices(matrix_s, subm_s, &
2572  dpattern, map, node_of_domain, select_row_col)
2573 
2574  ! loop over domains - perform inversion
2575  DO idomain = 1, ndomains
2576 
2577  ! check if the submatrix exists
2578  IF (subm_s(idomain)%domain .GT. 0) THEN
2579 
2580  naos = subm_s(idomain)%nrows
2581 
2582  ALLOCATE (sinv(naos, naos))
2583 
2584  CALL pseudo_invert_matrix(a=subm_s(idomain)%mdata, ainv=sinv, n=naos, &
2585  method=0)
2586 
2587  CALL copy_submatrices(subm_s(idomain), subm_s_inv(idomain), .false.)
2588  CALL copy_submatrix_data(sinv, subm_s_inv(idomain))
2589 
2590  DEALLOCATE (sinv)
2591 
2592  END IF ! submatrix for the domain exists
2593 
2594  END DO ! loop over domains
2595 
2596  CALL release_submatrices(subm_s)
2597  DEALLOCATE (subm_s)
2598 
2599  CALL timestop(handle)
2600 
2601  END SUBROUTINE construct_domain_s_inv
2602 
2603 ! **************************************************************************************************
2604 !> \brief Constructs subblocks of the covariant-covariant projectors (i.e. DM without spin factor)
2605 !> \param matrix_t ...
2606 !> \param matrix_sigma_inv ...
2607 !> \param matrix_s ...
2608 !> \param subm_r_down ...
2609 !> \param dpattern ...
2610 !> \param map ...
2611 !> \param node_of_domain ...
2612 !> \param filter_eps ...
2613 !> \par History
2614 !> 2013.02 created [Rustam Z. Khaliullin]
2615 !> \author Rustam Z. Khaliullin
2616 ! **************************************************************************************************
2617  SUBROUTINE construct_domain_r_down(matrix_t, matrix_sigma_inv, matrix_s, &
2618  subm_r_down, dpattern, map, node_of_domain, filter_eps)
2619 
2620  TYPE(dbcsr_type), INTENT(IN) :: matrix_t, matrix_sigma_inv, matrix_s
2621  TYPE(domain_submatrix_type), DIMENSION(:), &
2622  INTENT(INOUT) :: subm_r_down
2623  TYPE(dbcsr_type), INTENT(IN) :: dpattern
2624  TYPE(domain_map_type), INTENT(IN) :: map
2625  INTEGER, DIMENSION(:), INTENT(IN) :: node_of_domain
2626  REAL(kind=dp) :: filter_eps
2627 
2628  CHARACTER(len=*), PARAMETER :: routinen = 'construct_domain_r_down'
2629 
2630  INTEGER :: handle, ndomains
2631  TYPE(dbcsr_type) :: m_tmp_no_1, m_tmp_no_2, matrix_r
2632 
2633  CALL timeset(routinen, handle)
2634 
2635  ! compute the density matrix in the COVARIANT representation
2636  CALL dbcsr_create(matrix_r, &
2637  template=matrix_s, &
2638  matrix_type=dbcsr_type_symmetric)
2639  CALL dbcsr_create(m_tmp_no_1, &
2640  template=matrix_t, &
2641  matrix_type=dbcsr_type_no_symmetry)
2642  CALL dbcsr_create(m_tmp_no_2, &
2643  template=matrix_t, &
2644  matrix_type=dbcsr_type_no_symmetry)
2645 
2646  CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_s, matrix_t, &
2647  0.0_dp, m_tmp_no_1, filter_eps=filter_eps)
2648  CALL dbcsr_multiply("N", "N", 1.0_dp, m_tmp_no_1, matrix_sigma_inv, &
2649  0.0_dp, m_tmp_no_2, filter_eps=filter_eps)
2650  CALL dbcsr_multiply("N", "T", 1.0_dp, m_tmp_no_2, m_tmp_no_1, &
2651  0.0_dp, matrix_r, filter_eps=filter_eps)
2652 
2653  CALL dbcsr_release(m_tmp_no_1)
2654  CALL dbcsr_release(m_tmp_no_2)
2655 
2656  ndomains = dbcsr_nblkcols_total(dpattern)
2657  cpassert(SIZE(subm_r_down) .EQ. ndomains)
2658 
2659  CALL construct_submatrices(matrix_r, subm_r_down, &
2660  dpattern, map, node_of_domain, select_row_col)
2661 
2662  CALL dbcsr_release(matrix_r)
2663 
2664  CALL timestop(handle)
2665 
2666  END SUBROUTINE construct_domain_r_down
2667 
2668 ! **************************************************************************************************
2669 !> \brief Finds the square root of a matrix and its inverse
2670 !> \param A ...
2671 !> \param Asqrt ...
2672 !> \param Asqrtinv ...
2673 !> \param N ...
2674 !> \par History
2675 !> 2013.03 created [Rustam Z. Khaliullin]
2676 !> \author Rustam Z. Khaliullin
2677 ! **************************************************************************************************
2678  SUBROUTINE matrix_sqrt(A, Asqrt, Asqrtinv, N)
2679 
2680  REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: a
2681  REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: asqrt, asqrtinv
2682  INTEGER, INTENT(IN) :: n
2683 
2684  CHARACTER(len=*), PARAMETER :: routinen = 'matrix_sqrt'
2685 
2686  INTEGER :: handle, info, jj, lwork, unit_nr
2687  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: eigenvalues, work
2688  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: test, testn
2689  TYPE(cp_logger_type), POINTER :: logger
2690 
2691  CALL timeset(routinen, handle)
2692 
2693  asqrtinv = a
2694  info = 0
2695 
2696  ! get a useful unit_nr
2697  logger => cp_get_default_logger()
2698  IF (logger%para_env%is_source()) THEN
2699  unit_nr = cp_logger_get_default_unit_nr(logger, local=.true.)
2700  ELSE
2701  unit_nr = -1
2702  END IF
2703 
2704  !CALL DPOTRF('L', N, Ainv, N, INFO )
2705  !IF( INFO.NE.0 ) THEN
2706  ! CPErrorMessage(cp_failure_level,routineP,"DPOTRF failed")
2707  ! CPPrecondition(.FALSE.,cp_failure_level,routineP,failure)
2708  !END IF
2709  !CALL DPOTRI('L', N, Ainv, N, INFO )
2710  !IF( INFO.NE.0 ) THEN
2711  ! CPErrorMessage(cp_failure_level,routineP,"DPOTRI failed")
2712  ! CPPrecondition(.FALSE.,cp_failure_level,routineP,failure)
2713  !END IF
2714  !! complete the matrix
2715  !DO ii=1,N
2716  ! DO jj=ii+1,N
2717  ! Ainv(ii,jj)=Ainv(jj,ii)
2718  ! ENDDO
2719  ! !WRITE(*,'(100F13.9)') Ainv(ii,:)
2720  !ENDDO
2721 
2722  ! diagonalize first
2723  ALLOCATE (eigenvalues(n))
2724  ! Query the optimal workspace for dsyev
2725  lwork = -1
2726  ALLOCATE (work(max(1, lwork)))
2727  CALL dsyev('V', 'L', n, asqrtinv, n, eigenvalues, work, lwork, info)
2728  lwork = int(work(1))
2729  DEALLOCATE (work)
2730  ! Allocate the workspace and solve the eigenproblem
2731  ALLOCATE (work(max(1, lwork)))
2732  CALL dsyev('V', 'L', n, asqrtinv, n, eigenvalues, work, lwork, info)
2733  IF (info .NE. 0) THEN
2734  IF (unit_nr > 0) WRITE (unit_nr, *) 'DSYEV ERROR MESSAGE: ', info
2735  cpabort("DSYEV failed")
2736  END IF
2737  DEALLOCATE (work)
2738 
2739  ! take functions of eigenvalues and use eigenvectors to compute the matrix function
2740  ! first sqrt
2741  ALLOCATE (test(n, n))
2742  DO jj = 1, n
2743  test(jj, :) = asqrtinv(:, jj)*sqrt(eigenvalues(jj))
2744  END DO
2745  ALLOCATE (testn(n, n))
2746  testn(:, :) = matmul(asqrtinv, test)
2747  asqrt = testn
2748  ! now, sqrt_inv
2749  DO jj = 1, n
2750  test(jj, :) = asqrtinv(:, jj)/sqrt(eigenvalues(jj))
2751  END DO
2752  testn(:, :) = matmul(asqrtinv, test)
2753  asqrtinv = testn
2754  DEALLOCATE (test, testn)
2755 
2756  DEALLOCATE (eigenvalues)
2757 
2758 !!! ! compute the error
2759 !!! allocate(test(N,N))
2760 !!! test=MATMUL(Ainv,A)
2761 !!! DO ii=1,N
2762 !!! test(ii,ii)=test(ii,ii)-1.0_dp
2763 !!! ENDDO
2764 !!! test_error=0.0_dp
2765 !!! DO ii=1,N
2766 !!! DO jj=1,N
2767 !!! test_error=test_error+test(jj,ii)*test(jj,ii)
2768 !!! ENDDO
2769 !!! ENDDO
2770 !!! WRITE(*,*) "Inversion error: ", SQRT(test_error)
2771 !!! deallocate(test)
2772 
2773  CALL timestop(handle)
2774 
2775  END SUBROUTINE matrix_sqrt
2776 
2777 ! **************************************************************************************************
2778 !> \brief Inverts a matrix using a requested method
2779 !> 0. Cholesky factorization
2780 !> 1. Diagonalization
2781 !> \param A ...
2782 !> \param Ainv ...
2783 !> \param N ...
2784 !> \param method ...
2785 !> \param range1 ...
2786 !> \param range2 ...
2787 !> \param range1_thr ...
2788 !> \param shift ...
2789 !> \param bad_modes_projector_down ...
2790 !> \param s_inv_half ...
2791 !> \param s_half ...
2792 !> \par History
2793 !> 2012.04 created [Rustam Z. Khaliullin]
2794 !> \author Rustam Z. Khaliullin
2795 ! **************************************************************************************************
2796  SUBROUTINE pseudo_invert_matrix(A, Ainv, N, method, range1, range2, range1_thr, &
2797  shift, bad_modes_projector_down, s_inv_half, s_half)
2798 
2799  REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: a
2800  REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: ainv
2801  INTEGER, INTENT(IN) :: n, method
2802  INTEGER, INTENT(IN), OPTIONAL :: range1, range2
2803  REAL(kind=dp), INTENT(IN), OPTIONAL :: range1_thr, shift
2804  REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT), &
2805  OPTIONAL :: bad_modes_projector_down
2806  REAL(kind=dp), DIMENSION(:, :), INTENT(IN), &
2807  OPTIONAL :: s_inv_half, s_half
2808 
2809  CHARACTER(len=*), PARAMETER :: routinen = 'pseudo_invert_matrix'
2810 
2811  INTEGER :: handle, ii, info, jj, lwork, range1_eiv, &
2812  range2_eiv, range3_eiv, unit_nr
2813  LOGICAL :: use_both, use_ranges_only, use_thr_only
2814  REAL(kind=dp) :: my_shift
2815  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: eigenvalues, work
2816  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: temp1, temp2, temp3, temp4
2817  TYPE(cp_logger_type), POINTER :: logger
2818 
2819  CALL timeset(routinen, handle)
2820 
2821  ! get a useful unit_nr
2822  logger => cp_get_default_logger()
2823  IF (logger%para_env%is_source()) THEN
2824  unit_nr = cp_logger_get_default_unit_nr(logger, local=.true.)
2825  ELSE
2826  unit_nr = -1
2827  END IF
2828 
2829  IF (method .EQ. 1) THEN
2830 
2831  IF ((PRESENT(range1) .AND. (.NOT. PRESENT(range2))) .OR. (PRESENT(range2) .AND. (.NOT. PRESENT(range1)))) THEN
2832  cpabort("range1 and range2 must be provided together")
2833  END IF
2834 
2835  IF (PRESENT(range1) .AND. PRESENT(range1_thr)) THEN
2836  use_both = .true.
2837  use_thr_only = .false.
2838  use_ranges_only = .false.
2839  ELSE
2840  use_both = .false.
2841 
2842  IF (PRESENT(range1)) THEN
2843  use_ranges_only = .true.
2844  ELSE
2845  use_ranges_only = .false.
2846  END IF
2847 
2848  IF (PRESENT(range1_thr)) THEN
2849  use_thr_only = .true.
2850  ELSE
2851  use_thr_only = .false.
2852  END IF
2853 
2854  END IF
2855 
2856  IF ((PRESENT(s_half) .AND. (.NOT. PRESENT(s_inv_half))) .OR. (PRESENT(s_inv_half) .AND. (.NOT. PRESENT(s_half)))) THEN
2857  cpabort("Domain overlap matrix missing")
2858  END IF
2859  END IF
2860 
2861  my_shift = 0.0_dp
2862  IF (PRESENT(shift)) THEN
2863  my_shift = shift
2864  END IF
2865 
2866  ainv = a
2867  info = 0
2868 
2869  SELECT CASE (method)
2870  CASE (0) ! Inversion via cholesky factorization
2871 
2872  CALL dpotrf('L', n, ainv, n, info)
2873  IF (info .NE. 0) THEN
2874  cpabort("DPOTRF failed")
2875  END IF
2876  CALL dpotri('L', n, ainv, n, info)
2877  IF (info .NE. 0) THEN
2878  cpabort("DPOTRI failed")
2879  END IF
2880  ! complete the matrix
2881  DO ii = 1, n
2882  DO jj = ii + 1, n
2883  ainv(ii, jj) = ainv(jj, ii)
2884  END DO
2885  !WRITE(*,'(100F13.9)') Ainv(ii,:)
2886  END DO
2887 
2888  CASE (1)
2889 
2890  ! diagonalize first
2891  ALLOCATE (eigenvalues(n))
2892  ALLOCATE (temp1(n, n))
2893  ALLOCATE (temp4(n, n))
2894  IF (PRESENT(s_inv_half)) THEN
2895  CALL dsymm('L', 'U', n, n, 1.0_dp, s_inv_half, n, a, n, 0.0_dp, temp1, n)
2896  CALL dsymm('R', 'U', n, n, 1.0_dp, s_inv_half, n, temp1, n, 0.0_dp, ainv, n)
2897  END IF
2898  ! Query the optimal workspace for dsyev
2899  lwork = -1
2900  ALLOCATE (work(max(1, lwork)))
2901  CALL dsyev('V', 'L', n, ainv, n, eigenvalues, work, lwork, info)
2902 
2903  lwork = int(work(1))
2904  DEALLOCATE (work)
2905  ! Allocate the workspace and solve the eigenproblem
2906  ALLOCATE (work(max(1, lwork)))
2907  CALL dsyev('V', 'L', n, ainv, n, eigenvalues, work, lwork, info)
2908 
2909  IF (info .NE. 0) THEN
2910  IF (unit_nr > 0) WRITE (unit_nr, *) 'EIGENSYSTEM ERROR MESSAGE: ', info
2911  cpabort("Eigenproblem routine failed")
2912  END IF
2913  DEALLOCATE (work)
2914 
2915  !WRITE(*,*) "EIGENVALS: "
2916  !WRITE(*,'(4F13.9)') eigenvalues(:)
2917 
2918  ! invert eigenvalues and use eigenvectors to compute pseudo Ainv
2919  ! project out near-zero eigenvalue modes
2920  ALLOCATE (temp2(n, n))
2921  IF (PRESENT(bad_modes_projector_down)) ALLOCATE (temp3(n, n))
2922  temp2(1:n, 1:n) = ainv(1:n, 1:n)
2923 
2924  range1_eiv = 0
2925  range2_eiv = 0
2926  range3_eiv = 0
2927 
2928  IF (use_both) THEN
2929  DO jj = 1, n
2930  IF ((jj .LE. range2) .AND. (eigenvalues(jj) .LT. range1_thr)) THEN
2931  temp1(jj, :) = temp2(:, jj)*0.0_dp
2932  IF (PRESENT(bad_modes_projector_down)) temp3(jj, :) = ainv(:, jj)*1.0_dp
2933  range1_eiv = range1_eiv + 1
2934  ELSE
2935  temp1(jj, :) = temp2(:, jj)/(eigenvalues(jj) + my_shift)
2936  IF (PRESENT(bad_modes_projector_down)) temp3(jj, :) = ainv(:, jj)*0.0_dp
2937  range2_eiv = range2_eiv + 1
2938  END IF
2939  END DO
2940  ELSE
2941  IF (use_ranges_only) THEN
2942  DO jj = 1, n
2943  IF (jj .LE. range1) THEN
2944  temp1(jj, :) = temp2(:, jj)*0.0_dp
2945  IF (PRESENT(bad_modes_projector_down)) temp3(jj, :) = ainv(:, jj)*1.0_dp
2946  range1_eiv = range1_eiv + 1
2947  ELSE IF (jj .LE. range2) THEN
2948  temp1(jj, :) = temp2(:, jj)*1.0_dp
2949  IF (PRESENT(bad_modes_projector_down)) temp3(jj, :) = ainv(:, jj)*1.0_dp
2950  range2_eiv = range2_eiv + 1
2951  ELSE
2952  temp1(jj, :) = temp2(:, jj)/(eigenvalues(jj) + my_shift)
2953  IF (PRESENT(bad_modes_projector_down)) temp3(jj, :) = ainv(:, jj)*0.0_dp
2954  range3_eiv = range3_eiv + 1
2955  END IF
2956  END DO
2957  ELSE IF (use_thr_only) THEN
2958  DO jj = 1, n
2959  IF (eigenvalues(jj) .LT. range1_thr) THEN
2960  temp1(jj, :) = temp2(:, jj)*0.0_dp
2961  IF (PRESENT(bad_modes_projector_down)) temp3(jj, :) = ainv(:, jj)*1.0_dp
2962  range1_eiv = range1_eiv + 1
2963  ELSE
2964  temp1(jj, :) = temp2(:, jj)/(eigenvalues(jj) + my_shift)
2965  IF (PRESENT(bad_modes_projector_down)) temp3(jj, :) = ainv(:, jj)*0.0_dp
2966  range2_eiv = range2_eiv + 1
2967  END IF
2968  END DO
2969  ELSE ! no ranges, no thresholds
2970  cpabort("Invert using Cholesky. It would be faster.")
2971  END IF
2972  END IF
2973  !WRITE(*,*) ' EIV RANGES: ', range1_eiv, range2_eiv, range3_eiv
2974  IF (PRESENT(bad_modes_projector_down)) THEN
2975  IF (PRESENT(s_half)) THEN
2976  CALL dsymm('L', 'U', n, n, 1.0_dp, s_half, n, temp2, n, 0.0_dp, ainv, n)
2977  CALL dsymm('R', 'U', n, n, 1.0_dp, s_half, n, temp3, n, 0.0_dp, temp4, n)
2978  CALL dgemm('N', 'N', n, n, n, 1.0_dp, ainv, n, temp4, n, 0.0_dp, bad_modes_projector_down, n)
2979  ELSE
2980  CALL dgemm('N', 'N', n, n, n, 1.0_dp, temp2, n, temp3, n, 0.0_dp, bad_modes_projector_down, n)
2981  END IF
2982  END IF
2983 
2984  IF (PRESENT(s_inv_half)) THEN
2985  CALL dsymm('L', 'U', n, n, 1.0_dp, s_inv_half, n, temp2, n, 0.0_dp, temp4, n)
2986  CALL dsymm('R', 'U', n, n, 1.0_dp, s_inv_half, n, temp1, n, 0.0_dp, temp2, n)
2987  CALL dgemm('N', 'N', n, n, n, 1.0_dp, temp4, n, temp2, n, 0.0_dp, ainv, n)
2988  ELSE
2989  CALL dgemm('N', 'N', n, n, n, 1.0_dp, temp2, n, temp1, n, 0.0_dp, ainv, n)
2990  END IF
2991  DEALLOCATE (temp1, temp2, temp4)
2992  IF (PRESENT(bad_modes_projector_down)) DEALLOCATE (temp3)
2993  DEALLOCATE (eigenvalues)
2994 
2995  CASE DEFAULT
2996 
2997  cpabort("Illegal method selected for matrix inversion")
2998 
2999  END SELECT
3000 
3001  !! compute the inversion error
3002  !allocate(temp1(N,N))
3003  !temp1=MATMUL(Ainv,A)
3004  !DO ii=1,N
3005  ! temp1(ii,ii)=temp1(ii,ii)-1.0_dp
3006  !ENDDO
3007  !temp1_error=0.0_dp
3008  !DO ii=1,N
3009  ! DO jj=1,N
3010  ! temp1_error=temp1_error+temp1(jj,ii)*temp1(jj,ii)
3011  ! ENDDO
3012  !ENDDO
3013  !WRITE(*,*) "Inversion error: ", SQRT(temp1_error)
3014  !deallocate(temp1)
3015 
3016  CALL timestop(handle)
3017 
3018  END SUBROUTINE pseudo_invert_matrix
3019 
3020 ! **************************************************************************************************
3021 !> \brief Find matrix power using diagonalization
3022 !> \param A ...
3023 !> \param Apow ...
3024 !> \param power ...
3025 !> \param N ...
3026 !> \param range1 ...
3027 !> \param range1_thr ...
3028 !> \param shift ...
3029 !> \par History
3030 !> 2012.04 created [Rustam Z. Khaliullin]
3031 !> \author Rustam Z. Khaliullin
3032 ! **************************************************************************************************
3033 
3034  SUBROUTINE pseudo_matrix_power(A, Apow, power, N, range1, range1_thr, shift)
3035 
3036  REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: a
3037  REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: apow
3038  REAL(kind=dp), INTENT(IN) :: power
3039  INTEGER, INTENT(IN) :: n
3040  INTEGER, INTENT(IN), OPTIONAL :: range1
3041  REAL(kind=dp), INTENT(IN), OPTIONAL :: range1_thr, shift
3042 
3043  CHARACTER(len=*), PARAMETER :: routinen = 'pseudo_matrix_power'
3044 
3045  INTEGER :: handle, info, jj, lwork, range1_eiv, &
3046  range2_eiv, unit_nr
3047  LOGICAL :: use_both, use_ranges, use_thr
3048  REAL(kind=dp) :: my_shift
3049  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: eigenvalues, work
3050  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: temp1, temp2
3051  TYPE(cp_logger_type), POINTER :: logger
3052 
3053  CALL timeset(routinen, handle)
3054 
3055  ! get a useful unit_nr
3056  logger => cp_get_default_logger()
3057  IF (logger%para_env%is_source()) THEN
3058  unit_nr = cp_logger_get_default_unit_nr(logger, local=.true.)
3059  ELSE
3060  unit_nr = -1
3061  END IF
3062 
3063  IF (PRESENT(range1) .AND. PRESENT(range1_thr)) THEN
3064  use_both = .true.
3065  ELSE
3066  use_both = .false.
3067  IF (PRESENT(range1)) THEN
3068  use_ranges = .true.
3069  ELSE
3070  use_ranges = .false.
3071  IF (PRESENT(range1_thr)) THEN
3072  use_thr = .true.
3073  ELSE
3074  use_thr = .false.
3075  END IF
3076  END IF
3077  END IF
3078 
3079  my_shift = 0.0_dp
3080  IF (PRESENT(shift)) THEN
3081  my_shift = shift
3082  END IF
3083 
3084  apow = a
3085  info = 0
3086 
3087  ! diagonalize first
3088  ALLOCATE (eigenvalues(n))
3089  ALLOCATE (temp1(n, n))
3090 
3091  ! Query the optimal workspace for dsyev
3092  lwork = -1
3093  ALLOCATE (work(max(1, lwork)))
3094  CALL dsyev('V', 'L', n, apow, n, eigenvalues, work, lwork, info)
3095 
3096  lwork = int(work(1))
3097  DEALLOCATE (work)
3098  ! Allocate the workspace and solve the eigenproblem
3099  ALLOCATE (work(max(1, lwork)))
3100 
3101  CALL dsyev('V', 'L', n, apow, n, eigenvalues, work, lwork, info)
3102 
3103  IF (info .NE. 0) THEN
3104  IF (unit_nr > 0) WRITE (unit_nr, *) 'EIGENSYSTEM ERROR MESSAGE: ', info
3105  cpabort("Eigenproblem routine failed")
3106  END IF
3107  DEALLOCATE (work)
3108 
3109  !WRITE(*,*) "EIGENVALS: "
3110  !WRITE(*,'(4F13.9)') eigenvalues(:)
3111 
3112  ! invert eigenvalues and use eigenvectors to compute pseudo Ainv
3113  ! project out near-zero eigenvalue modes
3114  ALLOCATE (temp2(n, n))
3115 
3116  temp2(1:n, 1:n) = apow(1:n, 1:n)
3117 
3118  range1_eiv = 0
3119  range2_eiv = 0
3120 
3121  IF (use_both) THEN
3122  DO jj = 1, n
3123  IF ((jj .LE. range1) .AND. (eigenvalues(jj) .LT. range1_thr)) THEN
3124  temp1(jj, :) = temp2(:, jj)*0.0_dp
3125  range1_eiv = range1_eiv + 1
3126  ELSE
3127  temp1(jj, :) = temp2(:, jj)*((eigenvalues(jj) + my_shift)**power)
3128  END IF
3129  END DO
3130  ELSE
3131  IF (use_ranges) THEN
3132  DO jj = 1, n
3133  IF (jj .LE. range1) THEN
3134  temp1(jj, :) = temp2(:, jj)*0.0_dp
3135  range1_eiv = range1_eiv + 1
3136  ELSE
3137  temp1(jj, :) = temp2(:, jj)*((eigenvalues(jj) + my_shift)**power)
3138  END IF
3139  END DO
3140  ELSE
3141  IF (use_thr) THEN
3142  DO jj = 1, n
3143  IF (eigenvalues(jj) .LT. range1_thr) THEN
3144  temp1(jj, :) = temp2(:, jj)*0.0_dp
3145 
3146  range1_eiv = range1_eiv + 1
3147  ELSE
3148  temp1(jj, :) = temp2(:, jj)*((eigenvalues(jj) + my_shift)**power)
3149  END IF
3150  END DO
3151  ELSE
3152  DO jj = 1, n
3153  temp1(jj, :) = temp2(:, jj)*((eigenvalues(jj) + my_shift)**power)
3154  END DO
3155  END IF
3156  END IF
3157  END IF
3158  !WRITE(*,*) ' EIV RANGES: ', range1_eiv, range2_eiv, range3_eiv
3159  apow = matmul(temp2, temp1)
3160  DEALLOCATE (temp1, temp2)
3161  DEALLOCATE (eigenvalues)
3162 
3163  CALL timestop(handle)
3164 
3165  END SUBROUTINE pseudo_matrix_power
3166 
3167 ! **************************************************************************************************
3168 !> \brief Load balancing of the submatrix computations
3169 !> \param almo_scf_env ...
3170 !> \par History
3171 !> 2013.02 created [Rustam Z. Khaliullin]
3172 !> \author Rustam Z. Khaliullin
3173 ! **************************************************************************************************
3174  SUBROUTINE distribute_domains(almo_scf_env)
3175 
3176  TYPE(almo_scf_env_type), INTENT(INOUT) :: almo_scf_env
3177 
3178  CHARACTER(len=*), PARAMETER :: routinen = 'distribute_domains'
3179 
3180  INTEGER :: handle, idomain, least_loaded, nao, &
3181  ncpus, ndomains
3182  INTEGER, ALLOCATABLE, DIMENSION(:) :: index0
3183  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: cpu_load, domain_load
3184  TYPE(dbcsr_distribution_type) :: dist
3185 
3186  CALL timeset(routinen, handle)
3187 
3188  ndomains = almo_scf_env%ndomains
3189  CALL dbcsr_get_info(almo_scf_env%matrix_s(1), distribution=dist)
3190  CALL dbcsr_distribution_get(dist, numnodes=ncpus)
3191 
3192  ALLOCATE (domain_load(ndomains))
3193  DO idomain = 1, ndomains
3194  nao = almo_scf_env%nbasis_of_domain(idomain)
3195  domain_load(idomain) = (nao*nao*nao)*1.0_dp
3196  END DO
3197 
3198  ALLOCATE (index0(ndomains))
3199 
3200  CALL sort(domain_load, ndomains, index0)
3201 
3202  ALLOCATE (cpu_load(ncpus))
3203  cpu_load(:) = 0.0_dp
3204 
3205  DO idomain = 1, ndomains
3206  least_loaded = minloc(cpu_load, 1)
3207  cpu_load(least_loaded) = cpu_load(least_loaded) + domain_load(idomain)
3208  almo_scf_env%cpu_of_domain(index0(idomain)) = least_loaded - 1
3209  END DO
3210 
3211  DEALLOCATE (cpu_load)
3212  DEALLOCATE (index0)
3213  DEALLOCATE (domain_load)
3214 
3215  CALL timestop(handle)
3216 
3217  END SUBROUTINE distribute_domains
3218 
3219 ! **************************************************************************************************
3220 !> \brief Tests construction and release of domain submatrices
3221 !> \param matrix_no ...
3222 !> \param dpattern ...
3223 !> \param map ...
3224 !> \param node_of_domain ...
3225 !> \par History
3226 !> 2013.01 created [Rustam Z. Khaliullin]
3227 !> \author Rustam Z. Khaliullin
3228 ! **************************************************************************************************
3229  SUBROUTINE construct_test(matrix_no, dpattern, map, node_of_domain)
3230 
3231  TYPE(dbcsr_type), INTENT(IN) :: matrix_no, dpattern
3232  TYPE(domain_map_type), INTENT(IN) :: map
3233  INTEGER, DIMENSION(:), INTENT(IN) :: node_of_domain
3234 
3235  CHARACTER(len=*), PARAMETER :: routinen = 'construct_test'
3236 
3237  INTEGER :: groupid, handle, ndomains
3238  TYPE(dbcsr_type) :: copy1
3239  TYPE(domain_submatrix_type), ALLOCATABLE, &
3240  DIMENSION(:) :: subm_nn, subm_no
3241  TYPE(mp_comm_type) :: group
3242 
3243  CALL timeset(routinen, handle)
3244 
3245  ndomains = dbcsr_nblkcols_total(dpattern)
3246  CALL dbcsr_get_info(dpattern, group=groupid)
3247  CALL group%set_handle(groupid)
3248 
3249  ALLOCATE (subm_no(ndomains), subm_nn(ndomains))
3250  CALL init_submatrices(subm_no)
3251  CALL init_submatrices(subm_nn)
3252 
3253  !CALL dbcsr_print(matrix_nn)
3254  !CALL construct_submatrices(matrix_nn,subm_nn,dpattern,map,select_row_col)
3255  !CALL print_submatrices(subm_nn,Group)
3256 
3257  !CALL dbcsr_print(matrix_no)
3258  CALL construct_submatrices(matrix_no, subm_no, dpattern, map, node_of_domain, select_row)
3259  CALL print_submatrices(subm_no, group)
3260 
3261  CALL dbcsr_create(copy1, template=matrix_no)
3262  CALL dbcsr_copy(copy1, matrix_no)
3263  CALL dbcsr_print(copy1)
3264  CALL construct_dbcsr_from_submatrices(copy1, subm_no, dpattern)
3265  CALL dbcsr_print(copy1)
3266  CALL dbcsr_release(copy1)
3267 
3268  CALL release_submatrices(subm_no)
3269  CALL release_submatrices(subm_nn)
3270  DEALLOCATE (subm_no, subm_nn)
3271 
3272  CALL timestop(handle)
3273 
3274  END SUBROUTINE construct_test
3275 
3276 ! **************************************************************************************************
3277 !> \brief create the initial guess for XALMOs
3278 !> \param m_guess ...
3279 !> \param m_t_in ...
3280 !> \param m_t0 ...
3281 !> \param m_quench_t ...
3282 !> \param m_overlap ...
3283 !> \param m_sigma_tmpl ...
3284 !> \param nspins ...
3285 !> \param xalmo_history ...
3286 !> \param assume_t0_q0x ...
3287 !> \param optimize_theta ...
3288 !> \param envelope_amplitude ...
3289 !> \param eps_filter ...
3290 !> \param order_lanczos ...
3291 !> \param eps_lanczos ...
3292 !> \param max_iter_lanczos ...
3293 !> \param nocc_of_domain ...
3294 !> \par History
3295 !> 2016.11 created [Rustam Z Khaliullin]
3296 !> \author Rustam Z Khaliullin
3297 ! **************************************************************************************************
3298  SUBROUTINE xalmo_initial_guess(m_guess, m_t_in, m_t0, m_quench_t, &
3299  m_overlap, m_sigma_tmpl, nspins, xalmo_history, assume_t0_q0x, &
3300  optimize_theta, envelope_amplitude, eps_filter, order_lanczos, eps_lanczos, &
3301  max_iter_lanczos, nocc_of_domain)
3302 
3303  TYPE(dbcsr_type), DIMENSION(:), INTENT(INOUT) :: m_guess
3304  TYPE(dbcsr_type), DIMENSION(:), INTENT(IN) :: m_t_in, m_t0, m_quench_t
3305  TYPE(dbcsr_type), INTENT(IN) :: m_overlap
3306  TYPE(dbcsr_type), DIMENSION(:), INTENT(IN) :: m_sigma_tmpl
3307  INTEGER, INTENT(IN) :: nspins
3308  TYPE(almo_scf_history_type), INTENT(IN) :: xalmo_history
3309  LOGICAL, INTENT(IN) :: assume_t0_q0x, optimize_theta
3310  REAL(kind=dp), INTENT(IN) :: envelope_amplitude, eps_filter
3311  INTEGER, INTENT(IN) :: order_lanczos
3312  REAL(kind=dp), INTENT(IN) :: eps_lanczos
3313  INTEGER, INTENT(IN) :: max_iter_lanczos
3314  INTEGER, DIMENSION(:, :), INTENT(IN) :: nocc_of_domain
3315 
3316  CHARACTER(len=*), PARAMETER :: routinen = 'xalmo_initial_guess'
3317 
3318  INTEGER :: handle, iaspc, ispin, istore, naspc, &
3319  unit_nr
3320  LOGICAL :: aspc_guess
3321  REAL(kind=dp) :: alpha
3322  TYPE(cp_logger_type), POINTER :: logger
3323  TYPE(dbcsr_type) :: m_extrapolated, m_sigma_tmp
3324 
3325  CALL timeset(routinen, handle)
3326 
3327  ! get a useful output_unit
3328  logger => cp_get_default_logger()
3329  IF (logger%para_env%is_source()) THEN
3330  unit_nr = cp_logger_get_default_unit_nr(logger, local=.true.)
3331  ELSE
3332  unit_nr = -1
3333  END IF
3334 
3335  IF (optimize_theta) THEN
3336  cpwarn("unused option")
3337  ! just not to trigger unused variable
3338  alpha = envelope_amplitude
3339  END IF
3340 
3341  ! if extrapolation order is zero then the standard guess is used
3342  ! ... the number of stored history points will remain zero if extrapolation order is zero
3343  IF (xalmo_history%istore == 0) THEN
3344  aspc_guess = .false.
3345  ELSE
3346  aspc_guess = .true.
3347  END IF
3348 
3349  ! create initial guess
3350  IF (.NOT. aspc_guess) THEN
3351 
3352  DO ispin = 1, nspins
3353 
3354  ! zero initial guess for the delocalization amplitudes
3355  ! or the supplied guess for orbitals
3356  IF (assume_t0_q0x) THEN
3357  CALL dbcsr_set(m_guess(ispin), 0.0_dp)
3358  ELSE
3359  ! copy coefficients to m_guess
3360  CALL dbcsr_copy(m_guess(ispin), m_t_in(ispin))
3361  END IF
3362 
3363  END DO !ispins
3364 
3365  ELSE !aspc_guess
3366 
3367  CALL cite_reference(kolafa2004)
3368  CALL cite_reference(kuhne2007)
3369 
3370  naspc = min(xalmo_history%istore, xalmo_history%nstore)
3371  IF (unit_nr > 0) THEN
3372  WRITE (unit_nr, fmt="(/,T2,A,/,/,T3,A,I0)") &
3373  "Parameters for the always stable predictor-corrector (ASPC) method:", &
3374  "ASPC order: ", naspc
3375  END IF
3376 
3377  DO ispin = 1, nspins
3378 
3379  CALL dbcsr_create(m_extrapolated, &
3380  template=m_quench_t(ispin), matrix_type=dbcsr_type_no_symmetry)
3381  CALL dbcsr_create(m_sigma_tmp, &
3382  template=m_sigma_tmpl(ispin), matrix_type=dbcsr_type_no_symmetry)
3383 
3384  ! set to zero before accumulation
3385  CALL dbcsr_set(m_guess(ispin), 0.0_dp)
3386 
3387  ! extrapolation
3388  DO iaspc = 1, naspc
3389 
3390  istore = mod(xalmo_history%istore - iaspc, xalmo_history%nstore) + 1
3391  alpha = (-1.0_dp)**(iaspc + 1)*real(iaspc, kind=dp)* &
3392  binomial(2*naspc, naspc - iaspc)/binomial(2*naspc - 2, naspc - 1)
3393  IF (unit_nr > 0) THEN
3394  WRITE (unit_nr, fmt="(T3,A2,I0,A4,F10.6)") &
3395  "B(", iaspc, ") = ", alpha
3396  END IF
3397 
3398  ! m_extrapolated - initialize the correct sparsity pattern
3399  ! it must be kept throughout extrapolation
3400  CALL dbcsr_copy(m_extrapolated, m_quench_t(ispin))
3401 
3402  ! project t0 onto the previous DMs
3403  ! note that t0 is projected instead of any other matrix (e.g.
3404  ! t_SCF from the prev step or random t)
3405  ! this is done to keep orbitals phase (i.e. sign) the same as in
3406  ! t0. if this is not done then subtracting t0 on the next step
3407  ! will produce a terrible guess and extrapolation will fail
3408  CALL dbcsr_multiply("N", "N", 1.0_dp, &
3409  xalmo_history%matrix_p_up_down(ispin, istore), &
3410  m_t0(ispin), &
3411  0.0_dp, m_extrapolated, &
3412  retain_sparsity=.true.)
3413  ! normalize MOs
3414  CALL orthogonalize_mos(ket=m_extrapolated, &
3415  overlap=m_sigma_tmp, &
3416  metric=m_overlap, &
3417  retain_locality=.true., &
3418  only_normalize=.false., &
3419  nocc_of_domain=nocc_of_domain(:, ispin), &
3420  eps_filter=eps_filter, &
3421  order_lanczos=order_lanczos, &
3422  eps_lanczos=eps_lanczos, &
3423  max_iter_lanczos=max_iter_lanczos)
3424 
3425  ! now accumulate. correct sparsity is ensured
3426  CALL dbcsr_add(m_guess(ispin), m_extrapolated, &
3427  1.0_dp, (1.0_dp*alpha)/naspc)
3428 
3429  END DO !iaspc
3430 
3431  CALL dbcsr_release(m_extrapolated)
3432 
3433  ! normalize MOs
3434  CALL orthogonalize_mos(ket=m_guess(ispin), &
3435  overlap=m_sigma_tmp, &
3436  metric=m_overlap, &
3437  retain_locality=.true., &
3438  only_normalize=.false., &
3439  nocc_of_domain=nocc_of_domain(:, ispin), &
3440  eps_filter=eps_filter, &
3441  order_lanczos=order_lanczos, &
3442  eps_lanczos=eps_lanczos, &
3443  max_iter_lanczos=max_iter_lanczos)
3444 
3445  CALL dbcsr_release(m_sigma_tmp)
3446 
3447  ! project the t0 space out from the extrapolated state
3448  ! this can be done outside this subroutine
3449  IF (assume_t0_q0x) THEN
3450  CALL dbcsr_add(m_guess(ispin), m_t0(ispin), &
3451  1.0_dp, -1.0_dp)
3452  END IF !assume_t0_q0x
3453 
3454  END DO !ispin
3455 
3456  END IF !aspc_guess?
3457 
3458  CALL timestop(handle)
3459 
3460  END SUBROUTINE xalmo_initial_guess
3461 
3462 END MODULE almo_scf_methods
3463 
dgemm
static void dgemm(const char transa, const char transb, const int m, const int n, const int k, const double alpha, const double *a, const int lda, const double *b, const int ldb, const double beta, double *c, const int ldc)
Convenient wrapper to hide Fortran nature of dgemm_, swapping a and b.
Definition: grid_cpu_task_list.c:195
almo_scf_methods
Subroutines for ALMO SCF.
Definition: almo_scf_methods.F:15
almo_scf_methods::construct_domain_preconditioner
subroutine, public construct_domain_preconditioner(matrix_main, subm_s_inv, subm_s_inv_half, subm_s_half, subm_r_down, matrix_trimmer, dpattern, map, node_of_domain, preconditioner, bad_modes_projector_down, use_trimmer, eps_zero_eigenvalues, my_action, skip_inversion)
Constructs preconditioners for each domain -1. projected preconditioner 0. simple preconditioner.
Definition: almo_scf_methods.F:2236
almo_scf_methods::almo_scf_ks_xx_to_tv_xx
subroutine, public almo_scf_ks_xx_to_tv_xx(almo_scf_env)
ALMOs by diagonalizing the KS domain submatrices computes both the occupied and virtual orbitals.
Definition: almo_scf_methods.F:768
almo_scf_methods::xalmo_initial_guess
subroutine, public xalmo_initial_guess(m_guess, m_t_in, m_t0, m_quench_t, m_overlap, m_sigma_tmpl, nspins, xalmo_history, assume_t0_q0x, optimize_theta, envelope_amplitude, eps_filter, order_lanczos, eps_lanczos, max_iter_lanczos, nocc_of_domain)
create the initial guess for XALMOs
Definition: almo_scf_methods.F:3302
almo_scf_methods::generator_to_unitary
subroutine, public generator_to_unitary(X, U, eps_filter)
computes a unitary matrix from an arbitrary "generator" matrix U = ( 1 - X + tr(X) ) ( 1 + X - tr(X) ...
Definition: almo_scf_methods.F:2012
almo_scf_methods::construct_test
subroutine, public construct_test(matrix_no, dpattern, map, node_of_domain)
Tests construction and release of domain submatrices.
Definition: almo_scf_methods.F:3230
almo_scf_methods::distribute_domains
subroutine, public distribute_domains(almo_scf_env)
Load balancing of the submatrix computations.
Definition: almo_scf_methods.F:3175
almo_scf_methods::almo_scf_p_blk_to_t_blk
subroutine, public almo_scf_p_blk_to_t_blk(almo_scf_env, ionic)
computes occupied ALMOs from the superimposed atomic density blocks
Definition: almo_scf_methods.F:1176
almo_scf_methods::almo_scf_t_rescaling
subroutine, public almo_scf_t_rescaling(matrix_t, mo_energies, mu_of_domain, real_ne_of_domain, spin_kTS, smear_e_temp, ndomains, nocc_of_domain)
Apply an occupation-rescaling trick to ALMOs for smearing. Partially occupied orbitals are considered...
Definition: almo_scf_methods.F:1329
almo_scf_methods::pseudo_invert_diagonal_blk
subroutine, public pseudo_invert_diagonal_blk(matrix_in, matrix_out, nocc)
inverts block-diagonal blocks of a dbcsr_matrix
Definition: almo_scf_methods.F:1108
almo_scf_methods::almo_scf_ks_blk_to_tv_blk
subroutine, public almo_scf_ks_blk_to_tv_blk(almo_scf_env)
computes ALMOs by diagonalizing the projected blocked KS matrix uses the diagonalization code for blo...
Definition: almo_scf_methods.F:915
almo_scf_methods::apply_domain_operators
subroutine, public apply_domain_operators(matrix_in, matrix_out, operator1, operator2, dpattern, map, node_of_domain, my_action, filter_eps, matrix_trimmer, use_trimmer)
Parallel code for domain specific operations (my_action) 0. out = op1 * in.
Definition: almo_scf_methods.F:2109
almo_scf_methods::construct_domain_r_down
subroutine, public construct_domain_r_down(matrix_t, matrix_sigma_inv, matrix_s, subm_r_down, dpattern, map, node_of_domain, filter_eps)
Constructs subblocks of the covariant-covariant projectors (i.e. DM without spin factor)
Definition: almo_scf_methods.F:2619
almo_scf_methods::almo_scf_t_to_proj
subroutine, public almo_scf_t_to_proj(t, p, eps_filter, orthog_orbs, nocc_of_domain, s, sigma, sigma_inv, use_guess, smear, algorithm, para_env, blacs_env, eps_lanczos, max_iter_lanczos, inverse_accelerator, inv_eps_factor)
computes the idempotent density matrix from MOs MOs can be either orthogonal or non-orthogonal
Definition: almo_scf_methods.F:1630
almo_scf_methods::construct_domain_s_inv
subroutine, public construct_domain_s_inv(matrix_s, subm_s_inv, dpattern, map, node_of_domain)
Constructs S_inv block for each domain.
Definition: almo_scf_methods.F:2549
almo_scf_methods::almo_scf_ks_to_ks_blk
subroutine, public almo_scf_ks_to_ks_blk(almo_scf_env)
computes the projected KS from the total KS matrix also computes the DIIS error vector as a by-produc...
Definition: almo_scf_methods.F:483
almo_scf_methods::get_overlap
subroutine, public get_overlap(bra, ket, overlap, metric, retain_overlap_sparsity, eps_filter, smear)
Computes the overlap matrix of MO orbitals.
Definition: almo_scf_methods.F:1423
almo_scf_methods::fill_matrix_with_ones
subroutine, public fill_matrix_with_ones(matrix)
Fill all matrix blocks with 1.0_dp.
Definition: almo_scf_methods.F:95
almo_scf_methods::apply_projector
subroutine, public apply_projector(psi_in, psi_out, psi_projector, metric, project_out, psi_projector_orthogonal, proj_in_template, eps_filter, sig_inv_projector, sig_inv_template)
applies projector to the orbitals |psi_out> = P |psi_in> OR |psi_out> = (1-P) |psi_in>,...
Definition: almo_scf_methods.F:1856
almo_scf_methods::construct_domain_s_sqrt
subroutine, public construct_domain_s_sqrt(matrix_s, subm_s_sqrt, subm_s_sqrt_inv, dpattern, map, node_of_domain)
Constructs S^(+1/2) and S^(-1/2) submatrices for each domain.
Definition: almo_scf_methods.F:2476
almo_scf_methods::orthogonalize_mos
subroutine, public orthogonalize_mos(ket, overlap, metric, retain_locality, only_normalize, nocc_of_domain, eps_filter, order_lanczos, eps_lanczos, max_iter_lanczos, overlap_sqrti, smear)
orthogonalize MOs
Definition: almo_scf_methods.F:1508
almo_scf_methods::almo_scf_ks_to_ks_xx
subroutine, public almo_scf_ks_to_ks_xx(almo_scf_env)
builds projected KS matrices for the overlapping domains also computes the DIIS error vector as a by-...
Definition: almo_scf_methods.F:138
almo_scf_types
Types for all ALMO-based methods.
Definition: almo_scf_types.F:15
bibliography
collects all references to literature in CP2K as new algorithms / method are included from literature...
Definition: bibliography.F:28
bibliography::kuhne2007
integer, save, public kuhne2007
Definition: bibliography.F:43
bibliography::kolafa2004
integer, save, public kolafa2004
Definition: bibliography.F:43
cp_blacs_env
methods related to the blacs parallel environment
Definition: cp_blacs_env.F:15
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_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
domain_submatrix_methods
Subroutines to handle submatrices.
Definition: domain_submatrix_methods.F:14
domain_submatrix_methods::copy_submatrix_data
subroutine, public copy_submatrix_data(array, copy)
...
Definition: domain_submatrix_methods.F:259
domain_submatrix_methods::print_submatrices
subroutine, public print_submatrices(submatrices, mpgroup)
...
Definition: domain_submatrix_methods.F:1597
domain_submatrix_methods::construct_dbcsr_from_submatrices
subroutine, public construct_dbcsr_from_submatrices(matrix, submatrix, distr_pattern)
Constructs a DBCSR matrix from submatrices.
Definition: domain_submatrix_methods.F:1336
domain_submatrix_methods::construct_submatrices
subroutine, public construct_submatrices(matrix, submatrix, distr_pattern, domain_map, node_of_domain, job_type)
Constructs submatrices for each ALMO domain by collecting distributed DBCSR blocks to local arrays.
Definition: domain_submatrix_methods.F:839
domain_submatrix_types
Types to handle submatrices.
Definition: domain_submatrix_types.F:14
domain_submatrix_types::select_row
integer, parameter, public select_row
Definition: domain_submatrix_types.F:26
domain_submatrix_types::select_row_col
integer, parameter, public select_row_col
Definition: domain_submatrix_types.F:25
fermi_utils
deal with the Fermi distribution, compute it, fix mu, get derivs
Definition: fermi_utils.F:13
fermi_utils::fermifixed
subroutine, public fermifixed(f, mu, kTS, e, N, T, maxocc, estate, festate)
returns occupations according to Fermi-Dirac statistics for a given set of energies and number of ele...
Definition: fermi_utils.F:202
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition: input_constants.F:17
input_constants::spd_inversion_dense_cholesky
integer, parameter, public spd_inversion_dense_cholesky
Definition: input_constants.F:1041
input_constants::almo_domain_layout_molecular
integer, parameter, public almo_domain_layout_molecular
Definition: input_constants.F:976
input_constants::almo_scf_diag
integer, parameter, public almo_scf_diag
Definition: input_constants.F:997
input_constants::almo_mat_distr_atomic
integer, parameter, public almo_mat_distr_atomic
Definition: input_constants.F:980
input_constants::spd_inversion_ls_taylor
integer, parameter, public spd_inversion_ls_taylor
Definition: input_constants.F:1041
input_constants::spd_inversion_ls_hotelling
integer, parameter, public spd_inversion_ls_hotelling
Definition: input_constants.F:1041
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
iterate_matrix::invert_taylor
subroutine, public invert_taylor(matrix_inverse, matrix, threshold, use_inv_as_guess, norm_convergence, filter_eps, accelerator_order, max_iter_lanczos, eps_lanczos, silent)
invert a symmetric positive definite diagonally dominant matrix
Definition: iterate_matrix.F:285
iterate_matrix::invert_hotelling
subroutine, public invert_hotelling(matrix_inverse, matrix, threshold, use_inv_as_guess, norm_convergence, filter_eps, accelerator_order, max_iter_lanczos, eps_lanczos, silent)
invert a symmetric positive definite matrix by Hotelling's method explicit symmetrization makes this ...
Definition: iterate_matrix.F:472
kinds
Defines the basic variable types.
Definition: kinds.F:23
kinds::dp
integer, parameter, public dp
Definition: kinds.F:34
mathlib
Collection of simple mathematical functions and subroutines.
Definition: mathlib.F:15
mathlib::binomial
elemental real(kind=dp) function, public binomial(n, k)
The binomial coefficient n over k for 0 <= k <= n is calculated, otherwise zero is returned.
Definition: mathlib.F:206
message_passing
Interface to the message passing library MPI.
Definition: message_passing.F:23
preconditioner
computes preconditioners, and implements methods to apply them currently used in qs_ot
Definition: preconditioner.F:15
util
All kind of helpful little routines.
Definition: util.F:14
orbital
Orbital angular momentum.
Definition: grid_common.h:125