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