(git:495eafe)
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))) &
1624 cpabort("PARA and BLACS env are necessary for cholesky algorithm")
1625
1626 use_sigma_inv_guess = .false.
1627 IF (PRESENT(use_guess)) THEN
1628 use_sigma_inv_guess = use_guess
1629 END IF
1630
1631 IF (.NOT. PRESENT(smear)) THEN
1632 smearing = .false.
1633 ELSE
1634 smearing = smear
1635 END IF
1636
1637 my_accelerator = 1
1638 IF (PRESENT(inverse_accelerator)) my_accelerator = inverse_accelerator
1639
1640 my_inv_eps_factor = 10.0_dp
1641 IF (PRESENT(inv_eps_factor)) my_inv_eps_factor = inv_eps_factor
1642
1643 IF (orthog_orbs) THEN
1644
1645 CALL dbcsr_multiply("N", "T", 1.0_dp, t, t, &
1646 0.0_dp, p, filter_eps=eps_filter)
1647
1648 ELSE
1649
1650 CALL dbcsr_create(t_tmp, template=t)
1651
1652 ! TMP=S.T
1653 CALL dbcsr_multiply("N", "N", 1.0_dp, s, t, 0.0_dp, t_tmp, &
1654 filter_eps=eps_filter)
1655
1656 ! Sig=tr(T).TMP - get MO overlap
1657 CALL dbcsr_multiply("T", "N", 1.0_dp, t, t_tmp, 0.0_dp, sigma, &
1658 filter_eps=eps_filter)
1659
1660 !! If smearing ALMO is requested, apply correction of the occupation-rescaling trick
1661 !! (i.e. converting rescaled orbitals into selfish orbitals)
1662 !! (i.e. the diagonal blocks of the rescaled overlap must remain unscaled)
1663 !! Since we have here orthonormal MOs within a fragment, diagonal blocks are identity matrices
1664 !! Therefore, one only need to restore the diagonal to 1
1665 !! RS-WARNING: Assume orthonormal MOs within a fragment
1666 IF (smearing) THEN
1667 CALL dbcsr_get_info(sigma, nfullrows_total=dim0)
1668 ALLOCATE (diag_correction(dim0))
1669 diag_correction = 1.0_dp
1670 CALL dbcsr_set_diag(sigma, diag_correction)
1671 DEALLOCATE (diag_correction)
1672 END IF
1673
1674 ! invert MO overlap
1675 CALL set_zero_electron_blocks_in_mo_mo_matrix(sigma, nocc_of_domain, 1.0_dp)
1676 SELECT CASE (my_algorithm)
1677 CASE (spd_inversion_ls_taylor)
1678
1679 CALL invert_taylor( &
1680 matrix_inverse=sigma_inv, &
1681 matrix=sigma, &
1682 use_inv_as_guess=use_sigma_inv_guess, &
1683 threshold=eps_filter*my_inv_eps_factor, &
1684 filter_eps=eps_filter, &
1685 !accelerator_order=my_accelerator, &
1686 !eps_lanczos=eps_lanczos, &
1687 !max_iter_lanczos=max_iter_lanczos, &
1688 silent=.false.)
1689
1690 CASE (spd_inversion_ls_hotelling)
1691
1692 CALL invert_hotelling( &
1693 matrix_inverse=sigma_inv, &
1694 matrix=sigma, &
1695 use_inv_as_guess=use_sigma_inv_guess, &
1696 threshold=eps_filter*my_inv_eps_factor, &
1697 filter_eps=eps_filter, &
1698 accelerator_order=my_accelerator, &
1699 eps_lanczos=eps_lanczos, &
1700 max_iter_lanczos=max_iter_lanczos, &
1701 silent=.false.)
1702
1703 CASE (spd_inversion_dense_cholesky)
1704
1705 ! invert using cholesky
1706 CALL dbcsr_copy(sigma_inv, sigma)
1707 CALL cp_dbcsr_cholesky_decompose(sigma_inv, &
1708 para_env=para_env, &
1709 blacs_env=blacs_env)
1710 CALL cp_dbcsr_cholesky_invert(sigma_inv, &
1711 para_env=para_env, &
1712 blacs_env=blacs_env, &
1713 uplo_to_full=.true.)
1714 CALL dbcsr_filter(sigma_inv, eps_filter)
1715 CASE DEFAULT
1716 cpabort("Illegal MO overalp inversion algorithm")
1717 END SELECT
1718 CALL set_zero_electron_blocks_in_mo_mo_matrix(sigma, nocc_of_domain, 0.0_dp)
1719 CALL set_zero_electron_blocks_in_mo_mo_matrix(sigma_inv, nocc_of_domain, 0.0_dp)
1720
1721 ! TMP=T.SigInv
1722 CALL dbcsr_multiply("N", "N", 1.0_dp, t, sigma_inv, 0.0_dp, t_tmp, &
1723 filter_eps=eps_filter)
1724
1725 ! P=TMP.tr(T_blk)
1726 CALL dbcsr_multiply("N", "T", 1.0_dp, t_tmp, t, 0.0_dp, p, &
1727 filter_eps=eps_filter)
1728
1729 CALL dbcsr_release(t_tmp)
1730
1731 END IF
1732
1733 CALL timestop(handle)
1734
1735 END SUBROUTINE almo_scf_t_to_proj
1736
1737! **************************************************************************************************
1738!> \brief self-explanatory
1739!> \param matrix ...
1740!> \param nocc_of_domain ...
1741!> \param value ...
1742!> \param
1743!> \par History
1744!> 2016.12 created [Rustam Z Khaliullin]
1745!> \author Rustam Z Khaliullin
1746! **************************************************************************************************
1747 SUBROUTINE set_zero_electron_blocks_in_mo_mo_matrix(matrix, nocc_of_domain, value)
1748
1749 TYPE(dbcsr_type), INTENT(INOUT) :: matrix
1750 INTEGER, DIMENSION(:), INTENT(IN) :: nocc_of_domain
1751 REAL(kind=dp), INTENT(IN) :: value
1752
1753 INTEGER :: col, row
1754 REAL(kind=dp), DIMENSION(:, :), POINTER :: block
1755 TYPE(dbcsr_iterator_type) :: iter
1756
1757 CALL dbcsr_reserve_diag_blocks(matrix)
1758
1759 CALL dbcsr_iterator_start(iter, matrix)
1760 DO WHILE (dbcsr_iterator_blocks_left(iter))
1761 CALL dbcsr_iterator_next_block(iter, row, col, block)
1762 IF (row == col .AND. nocc_of_domain(row) == 0) THEN
1763 block(1, 1) = value
1764 END IF
1765 END DO
1766 CALL dbcsr_iterator_stop(iter)
1767
1768 END SUBROUTINE set_zero_electron_blocks_in_mo_mo_matrix
1769
1770! **************************************************************************************************
1771!> \brief applies projector to the orbitals
1772!> |psi_out> = P |psi_in> OR |psi_out> = (1-P) |psi_in>,
1773!> where P = |psi_proj> (<psi_proj|psi_roj>)^{-1} <psi_proj|
1774!> \param psi_in ...
1775!> \param psi_out ...
1776!> \param psi_projector ...
1777!> \param metric ...
1778!> \param project_out ...
1779!> \param psi_projector_orthogonal ...
1780!> \param proj_in_template ...
1781!> \param eps_filter ...
1782!> \param sig_inv_projector ...
1783!> \param sig_inv_template ...
1784!> \par History
1785!> 2011.10 created [Rustam Z Khaliullin]
1786!> \author Rustam Z Khaliullin
1787! **************************************************************************************************
1788 SUBROUTINE apply_projector(psi_in, psi_out, psi_projector, metric, project_out, &
1789 psi_projector_orthogonal, proj_in_template, eps_filter, sig_inv_projector, &
1790 sig_inv_template)
1791
1792 TYPE(dbcsr_type), INTENT(IN) :: psi_in
1793 TYPE(dbcsr_type), INTENT(INOUT) :: psi_out
1794 TYPE(dbcsr_type), INTENT(IN) :: psi_projector, metric
1795 LOGICAL, INTENT(IN) :: project_out, psi_projector_orthogonal
1796 TYPE(dbcsr_type), INTENT(IN) :: proj_in_template
1797 REAL(kind=dp), INTENT(IN) :: eps_filter
1798 TYPE(dbcsr_type), INTENT(IN), OPTIONAL :: sig_inv_projector, sig_inv_template
1799
1800 CHARACTER(LEN=*), PARAMETER :: routinen = 'apply_projector'
1801
1802 INTEGER :: handle
1803 TYPE(dbcsr_type) :: tmp_no, tmp_ov, tmp_ov2, tmp_sig, &
1804 tmp_sig_inv
1805
1806 CALL timeset(routinen, handle)
1807
1808 ! =S*PSI_proj
1809 CALL dbcsr_create(tmp_no, template=psi_projector)
1810 CALL dbcsr_multiply("N", "N", 1.0_dp, &
1811 metric, psi_projector, &
1812 0.0_dp, tmp_no, &
1813 filter_eps=eps_filter)
1814
1815 ! =tr(S.PSI_proj)*PSI_in
1816 CALL dbcsr_create(tmp_ov, template=proj_in_template)
1817 CALL dbcsr_multiply("T", "N", 1.0_dp, &
1818 tmp_no, psi_in, &
1819 0.0_dp, tmp_ov, &
1820 filter_eps=eps_filter)
1821
1822 IF (.NOT. psi_projector_orthogonal) THEN
1823 ! =SigInv_proj*Sigma_OV
1824 CALL dbcsr_create(tmp_ov2, &
1825 template=proj_in_template)
1826 IF (PRESENT(sig_inv_projector)) THEN
1827 CALL dbcsr_create(tmp_sig_inv, &
1828 template=sig_inv_projector)
1829 CALL dbcsr_copy(tmp_sig_inv, sig_inv_projector)
1830 ELSE
1831 IF (.NOT. PRESENT(sig_inv_template)) THEN
1832 cpabort("PROGRAMMING ERROR: provide either template or sig_inv")
1833 END IF
1834 ! compute inverse overlap of the projector orbitals
1835 CALL dbcsr_create(tmp_sig, &
1836 template=sig_inv_template, &
1837 matrix_type=dbcsr_type_no_symmetry)
1838 CALL dbcsr_multiply("T", "N", 1.0_dp, &
1839 psi_projector, tmp_no, 0.0_dp, tmp_sig, &
1840 filter_eps=eps_filter)
1841 CALL dbcsr_create(tmp_sig_inv, &
1842 template=sig_inv_template, &
1843 matrix_type=dbcsr_type_no_symmetry)
1844 CALL invert_hotelling(tmp_sig_inv, tmp_sig, &
1845 threshold=eps_filter)
1846 CALL dbcsr_release(tmp_sig)
1847 END IF
1848 CALL dbcsr_multiply("N", "N", 1.0_dp, &
1849 tmp_sig_inv, tmp_ov, 0.0_dp, tmp_ov2, &
1850 filter_eps=eps_filter)
1851 CALL dbcsr_release(tmp_sig_inv)
1852 CALL dbcsr_copy(tmp_ov, tmp_ov2)
1853 CALL dbcsr_release(tmp_ov2)
1854 END IF
1855 CALL dbcsr_release(tmp_no)
1856
1857 ! =PSI_proj*TMP_OV
1858 CALL dbcsr_multiply("N", "N", 1.0_dp, &
1859 psi_projector, tmp_ov, 0.0_dp, psi_out, &
1860 filter_eps=eps_filter)
1861 CALL dbcsr_release(tmp_ov)
1862
1863 ! V_out=V_in-V_out
1864 IF (project_out) THEN
1865 CALL dbcsr_add(psi_out, psi_in, -1.0_dp, +1.0_dp)
1866 END IF
1867
1868 CALL timestop(handle)
1869
1870 END SUBROUTINE apply_projector
1871
1872!! **************************************************************************************************
1873!!> \brief projects the occupied space out from the provided orbitals
1874!!> \par History
1875!!> 2011.07 created [Rustam Z Khaliullin]
1876!!> \author Rustam Z Khaliullin
1877!! **************************************************************************************************
1878! SUBROUTINE almo_scf_p_out_from_v(v_in,v_out,ov_template,ispin,almo_scf_env)
1879!
1880! TYPE(dbcsr_type), INTENT(IN) :: v_in, ov_template
1881! TYPE(dbcsr_type), INTENT(INOUT) :: v_out
1882! INTEGER, INTENT(IN) :: ispin
1883! TYPE(almo_scf_env_type), INTENT(INOUT) :: almo_scf_env
1884!
1885! CHARACTER(LEN=*), PARAMETER :: &
1886! routineN = 'almo_scf_p_out_from_v', &
1887! routineP = moduleN//':'//routineN
1888!
1889! TYPE(dbcsr_type) :: tmp_on, tmp_ov, tmp_ov2
1890! INTEGER :: handle
1891! LOGICAL :: failure
1892!
1893! CALL timeset(routineN,handle)
1894!
1895! ! =tr(T_blk)*S
1896! CALL dbcsr_init(tmp_on)
1897! CALL dbcsr_create(tmp_on,&
1898! template=almo_scf_env%matrix_t_tr(ispin))
1899! CALL dbcsr_multiply("T","N",1.0_dp,&
1900! almo_scf_env%matrix_t_blk(ispin),&
1901! almo_scf_env%matrix_s(1),&
1902! 0.0_dp,tmp_on,&
1903! filter_eps=almo_scf_env%eps_filter)
1904!
1905! ! =tr(T_blk).S*V_in
1906! CALL dbcsr_init(tmp_ov)
1907! CALL dbcsr_create(tmp_ov,template=ov_template)
1908! CALL dbcsr_multiply("N","N",1.0_dp,&
1909! tmp_on,v_in,0.0_dp,tmp_ov,&
1910! filter_eps=almo_scf_env%eps_filter)
1911! CALL dbcsr_release(tmp_on)
1912!
1913! ! =SigmaInv*Sigma_OV
1914! CALL dbcsr_init(tmp_ov2)
1915! CALL dbcsr_create(tmp_ov2,template=ov_template)
1916! CALL dbcsr_multiply("N","N",1.0_dp,&
1917! almo_scf_env%matrix_sigma_inv(ispin),&
1918! tmp_ov,0.0_dp,tmp_ov2,&
1919! filter_eps=almo_scf_env%eps_filter)
1920! CALL dbcsr_release(tmp_ov)
1921!
1922! ! =T_blk*SigmaInv.Sigma_OV
1923! CALL dbcsr_multiply("N","N",1.0_dp,&
1924! almo_scf_env%matrix_t_blk(ispin),&
1925! tmp_ov2,0.0_dp,v_out,&
1926! filter_eps=almo_scf_env%eps_filter)
1927! CALL dbcsr_release(tmp_ov2)
1928!
1929! ! V_out=V_in-V_out=
1930! CALL dbcsr_add(v_out,v_in,-1.0_dp,+1.0_dp)
1931!
1932! CALL timestop(handle)
1933!
1934! END SUBROUTINE almo_scf_p_out_from_v
1935
1936! **************************************************************************************************
1937!> \brief computes a unitary matrix from an arbitrary "generator" matrix
1938!> U = ( 1 - X + tr(X) ) ( 1 + X - tr(X) )^(-1)
1939!> \param X ...
1940!> \param U ...
1941!> \param eps_filter ...
1942!> \par History
1943!> 2011.08 created [Rustam Z Khaliullin]
1944!> \author Rustam Z Khaliullin
1945! **************************************************************************************************
1946 SUBROUTINE generator_to_unitary(X, U, eps_filter)
1947
1948 TYPE(dbcsr_type), INTENT(IN) :: x
1949 TYPE(dbcsr_type), INTENT(INOUT) :: u
1950 REAL(kind=dp), INTENT(IN) :: eps_filter
1951
1952 CHARACTER(LEN=*), PARAMETER :: routinen = 'generator_to_unitary'
1953
1954 INTEGER :: handle, unit_nr
1955 LOGICAL :: safe_mode
1956 REAL(kind=dp) :: frob_matrix, frob_matrix_base
1957 TYPE(cp_logger_type), POINTER :: logger
1958 TYPE(dbcsr_type) :: delta, t1, t2, tmp1
1959
1960 CALL timeset(routinen, handle)
1961
1962 safe_mode = .true.
1963
1964 ! get a useful output_unit
1965 logger => cp_get_default_logger()
1966 IF (logger%para_env%is_source()) THEN
1967 unit_nr = cp_logger_get_default_unit_nr(logger, local=.true.)
1968 ELSE
1969 unit_nr = -1
1970 END IF
1971
1972 CALL dbcsr_create(t1, template=x, &
1973 matrix_type=dbcsr_type_no_symmetry)
1974 CALL dbcsr_create(t2, template=x, &
1975 matrix_type=dbcsr_type_no_symmetry)
1976
1977 ! create antisymmetric Delta = -X + tr(X)
1978 CALL dbcsr_create(delta, template=x, &
1979 matrix_type=dbcsr_type_no_symmetry)
1980 CALL dbcsr_transposed(delta, x)
1981! check that transposed is added correctly
1982 CALL dbcsr_add(delta, x, 1.0_dp, -1.0_dp)
1983
1984 ! compute (1 - Delta)^(-1)
1985 CALL dbcsr_add_on_diag(t1, 1.0_dp)
1986 CALL dbcsr_add(t1, delta, 1.0_dp, -1.0_dp)
1987 CALL invert_hotelling(t2, t1, threshold=eps_filter)
1988
1989 IF (safe_mode) THEN
1990
1991 CALL dbcsr_create(tmp1, template=x, &
1992 matrix_type=dbcsr_type_no_symmetry)
1993 CALL dbcsr_multiply("N", "N", 1.0_dp, t2, t1, 0.0_dp, tmp1, &
1994 filter_eps=eps_filter)
1995 frob_matrix_base = dbcsr_frobenius_norm(tmp1)
1996 CALL dbcsr_add_on_diag(tmp1, -1.0_dp)
1997 frob_matrix = dbcsr_frobenius_norm(tmp1)
1998 IF (unit_nr > 0) WRITE (unit_nr, *) "Error for (inv(A)*A-I)", frob_matrix/frob_matrix_base
1999 CALL dbcsr_release(tmp1)
2000 END IF
2001
2002 CALL dbcsr_multiply("N", "N", 1.0_dp, delta, t2, 0.0_dp, u, &
2003 filter_eps=eps_filter)
2004 CALL dbcsr_add(u, t2, 1.0_dp, 1.0_dp)
2005
2006 IF (safe_mode) THEN
2007
2008 CALL dbcsr_create(tmp1, template=x, &
2009 matrix_type=dbcsr_type_no_symmetry)
2010 CALL dbcsr_multiply("T", "N", 1.0_dp, u, u, 0.0_dp, tmp1, &
2011 filter_eps=eps_filter)
2012 frob_matrix_base = dbcsr_frobenius_norm(tmp1)
2013 CALL dbcsr_add_on_diag(tmp1, -1.0_dp)
2014 frob_matrix = dbcsr_frobenius_norm(tmp1)
2015 IF (unit_nr > 0) WRITE (unit_nr, *) "Error for (trn(U)*U-I)", frob_matrix/frob_matrix_base
2016 CALL dbcsr_release(tmp1)
2017 END IF
2018
2019 CALL timestop(handle)
2020
2021 END SUBROUTINE generator_to_unitary
2022
2023! **************************************************************************************************
2024!> \brief Parallel code for domain specific operations (my_action)
2025!> 0. out = op1 * in
2026!> 1. out = in - op2 * op1 * in
2027!> \param matrix_in ...
2028!> \param matrix_out ...
2029!> \param operator1 ...
2030!> \param operator2 ...
2031!> \param dpattern ...
2032!> \param map ...
2033!> \param node_of_domain ...
2034!> \param my_action ...
2035!> \param filter_eps ...
2036!> \param matrix_trimmer ...
2037!> \param use_trimmer ...
2038!> \par History
2039!> 2013.01 created [Rustam Z. Khaliullin]
2040!> \author Rustam Z. Khaliullin
2041! **************************************************************************************************
2042 SUBROUTINE apply_domain_operators(matrix_in, matrix_out, operator1, operator2, &
2043 dpattern, map, node_of_domain, my_action, filter_eps, matrix_trimmer, use_trimmer)
2044
2045 TYPE(dbcsr_type), INTENT(INOUT) :: matrix_in, matrix_out
2046 TYPE(domain_submatrix_type), DIMENSION(:), &
2047 INTENT(IN) :: operator1
2048 TYPE(domain_submatrix_type), DIMENSION(:), &
2049 INTENT(IN), OPTIONAL :: operator2
2050 TYPE(dbcsr_type), INTENT(IN) :: dpattern
2051 TYPE(domain_map_type), INTENT(IN) :: map
2052 INTEGER, DIMENSION(:), INTENT(IN) :: node_of_domain
2053 INTEGER, INTENT(IN) :: my_action
2054 REAL(kind=dp) :: filter_eps
2055 TYPE(dbcsr_type), INTENT(IN), OPTIONAL :: matrix_trimmer
2056 LOGICAL, INTENT(IN), OPTIONAL :: use_trimmer
2057
2058 CHARACTER(len=*), PARAMETER :: routinen = 'apply_domain_operators'
2059
2060 INTEGER :: handle, ndomains
2061 LOGICAL :: matrix_trimmer_required, my_use_trimmer, &
2062 operator2_required
2063 TYPE(domain_submatrix_type), ALLOCATABLE, &
2064 DIMENSION(:) :: subm_in, subm_out, subm_temp
2065
2066 CALL timeset(routinen, handle)
2067
2068 my_use_trimmer = .false.
2069 IF (PRESENT(use_trimmer)) THEN
2070 my_use_trimmer = use_trimmer
2071 END IF
2072
2073 operator2_required = .false.
2074 matrix_trimmer_required = .false.
2075
2076 IF (my_action == 1) operator2_required = .true.
2077
2078 IF (my_use_trimmer) THEN
2079 matrix_trimmer_required = .true.
2080 cpabort("TRIMMED PROJECTOR DISABLED!")
2081 END IF
2082
2083 IF (.NOT. PRESENT(operator2) .AND. operator2_required) THEN
2084 cpabort("SECOND OPERATOR IS REQUIRED")
2085 END IF
2086 IF (.NOT. PRESENT(matrix_trimmer) .AND. matrix_trimmer_required) THEN
2087 cpabort("TRIMMER MATRIX IS REQUIRED")
2088 END IF
2089
2090 CALL dbcsr_get_info(dpattern, nblkcols_total=ndomains)
2091
2092 ALLOCATE (subm_in(ndomains))
2093 ALLOCATE (subm_temp(ndomains))
2094 ALLOCATE (subm_out(ndomains))
2095 !!!TRIM ALLOCATE(subm_trimmer(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 !!!TRIM IF (matrix_trimmer_required) THEN
2104 !!!TRIM CALL construct_submatrices(matrix_trimmer,subm_trimmer,&
2105 !!!TRIM dpattern,map,node_of_domain,select_row)
2106 !!!TRIM ENDIF
2107
2108 IF (my_action == 0) THEN
2109 ! for example, apply preconditioner
2110 CALL multiply_submatrices('N', 'N', 1.0_dp, operator1, &
2111 subm_in, 0.0_dp, subm_out)
2112 ELSE IF (my_action == 1) THEN
2113 ! use for projectors
2114 CALL copy_submatrices(subm_in, subm_out, .true.)
2115 CALL multiply_submatrices('N', 'N', 1.0_dp, operator1, &
2116 subm_in, 0.0_dp, subm_temp)
2117 CALL multiply_submatrices('N', 'N', -1.0_dp, operator2, &
2118 subm_temp, 1.0_dp, subm_out)
2119 ELSE
2120 cpabort("ILLEGAL ACTION")
2121 END IF
2122
2123 CALL construct_dbcsr_from_submatrices(matrix_out, subm_out, dpattern)
2124 CALL dbcsr_filter(matrix_out, filter_eps)
2125
2126 CALL release_submatrices(subm_out)
2127 CALL release_submatrices(subm_temp)
2128 CALL release_submatrices(subm_in)
2129
2130 DEALLOCATE (subm_out)
2131 DEALLOCATE (subm_temp)
2132 DEALLOCATE (subm_in)
2133
2134 CALL timestop(handle)
2135
2136 END SUBROUTINE apply_domain_operators
2137
2138! **************************************************************************************************
2139!> \brief Constructs preconditioners for each domain
2140!> -1. projected preconditioner
2141!> 0. simple preconditioner
2142!> \param matrix_main ...
2143!> \param subm_s_inv ...
2144!> \param subm_s_inv_half ...
2145!> \param subm_s_half ...
2146!> \param subm_r_down ...
2147!> \param matrix_trimmer ...
2148!> \param dpattern ...
2149!> \param map ...
2150!> \param node_of_domain ...
2151!> \param preconditioner ...
2152!> \param bad_modes_projector_down ...
2153!> \param use_trimmer ...
2154!> \param eps_zero_eigenvalues ...
2155!> \param my_action ...
2156!> \param skip_inversion ...
2157!> \par History
2158!> 2013.01 created [Rustam Z. Khaliullin]
2159!> \author Rustam Z. Khaliullin
2160! **************************************************************************************************
2161 SUBROUTINE construct_domain_preconditioner(matrix_main, subm_s_inv, &
2162 subm_s_inv_half, subm_s_half, &
2163 subm_r_down, matrix_trimmer, &
2164 dpattern, map, node_of_domain, &
2165 preconditioner, &
2166 bad_modes_projector_down, &
2167 use_trimmer, &
2168 eps_zero_eigenvalues, &
2169 my_action, skip_inversion)
2170
2171 TYPE(dbcsr_type), INTENT(INOUT) :: matrix_main
2172 TYPE(domain_submatrix_type), DIMENSION(:), &
2173 INTENT(IN), OPTIONAL :: subm_s_inv, subm_s_inv_half, &
2174 subm_s_half, subm_r_down
2175 TYPE(dbcsr_type), INTENT(IN), OPTIONAL :: matrix_trimmer
2176 TYPE(dbcsr_type), INTENT(IN) :: dpattern
2177 TYPE(domain_map_type), INTENT(IN) :: map
2178 INTEGER, DIMENSION(:), INTENT(IN) :: node_of_domain
2179 TYPE(domain_submatrix_type), DIMENSION(:), &
2180 INTENT(INOUT) :: preconditioner
2181 TYPE(domain_submatrix_type), DIMENSION(:), &
2182 INTENT(INOUT), OPTIONAL :: bad_modes_projector_down
2183 LOGICAL, INTENT(IN), OPTIONAL :: use_trimmer
2184 REAL(kind=dp), INTENT(IN), OPTIONAL :: eps_zero_eigenvalues
2185 INTEGER, INTENT(IN) :: my_action
2186 LOGICAL, INTENT(IN), OPTIONAL :: skip_inversion
2187
2188 CHARACTER(len=*), PARAMETER :: routinen = 'construct_domain_preconditioner'
2189
2190 INTEGER :: handle, idomain, index1_end, &
2191 index1_start, n_domain_mos, naos, &
2192 nblkrows_tot, ndomains, neighbor, row
2193 INTEGER, DIMENSION(:), POINTER :: nmos
2194 LOGICAL :: eps_zero_eigenvalues_required, matrix_r_required, matrix_s_half_required, &
2195 matrix_s_inv_half_required, matrix_s_inv_required, matrix_trimmer_required, &
2196 my_skip_inversion, my_use_trimmer
2197 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: minv, proj_array
2198 TYPE(domain_submatrix_type), ALLOCATABLE, &
2199 DIMENSION(:) :: subm_main, subm_tmp, subm_tmp2
2200
2201 CALL timeset(routinen, handle)
2202
2203 my_use_trimmer = .false.
2204 IF (PRESENT(use_trimmer)) THEN
2205 my_use_trimmer = use_trimmer
2206 END IF
2207
2208 my_skip_inversion = .false.
2209 IF (PRESENT(skip_inversion)) THEN
2210 my_skip_inversion = skip_inversion
2211 END IF
2212
2213 matrix_s_inv_half_required = .false.
2214 matrix_s_half_required = .false.
2215 eps_zero_eigenvalues_required = .false.
2216 matrix_s_inv_required = .false.
2217 matrix_trimmer_required = .false.
2218 matrix_r_required = .false.
2219
2220 IF (my_action == -1) matrix_s_inv_required = .true.
2221 IF (my_action == -1) matrix_r_required = .true.
2222 IF (my_use_trimmer) THEN
2223 matrix_trimmer_required = .true.
2224 cpabort("TRIMMED PRECONDITIONER DISABLED!")
2225 END IF
2226 ! tie the following optional arguments together to prevent bad calls
2227 IF (PRESENT(bad_modes_projector_down)) THEN
2228 matrix_s_inv_half_required = .true.
2229 matrix_s_half_required = .true.
2230 eps_zero_eigenvalues_required = .true.
2231 END IF
2232
2233 ! check if all required optional arguments are provided
2234 IF (.NOT. PRESENT(subm_s_inv_half) .AND. matrix_s_inv_half_required) THEN
2235 cpabort("S_inv_half SUBMATRICES ARE REQUIRED")
2236 END IF
2237 IF (.NOT. PRESENT(subm_s_half) .AND. matrix_s_half_required) THEN
2238 cpabort("S_half SUBMATRICES ARE REQUIRED")
2239 END IF
2240 IF (.NOT. PRESENT(eps_zero_eigenvalues) .AND. eps_zero_eigenvalues_required) THEN
2241 cpabort("EPS_ZERO_EIGENVALUES IS REQUIRED")
2242 END IF
2243 IF (.NOT. PRESENT(subm_s_inv) .AND. matrix_s_inv_required) THEN
2244 cpabort("S_inv SUBMATRICES ARE REQUIRED")
2245 END IF
2246 IF (.NOT. PRESENT(subm_r_down) .AND. matrix_r_required) THEN
2247 cpabort("R SUBMATRICES ARE REQUIRED")
2248 END IF
2249 IF (.NOT. PRESENT(matrix_trimmer) .AND. matrix_trimmer_required) THEN
2250 cpabort("TRIMMER MATRIX IS REQUIRED")
2251 END IF
2252
2253 CALL dbcsr_get_info(dpattern, &
2254 nblkcols_total=ndomains, &
2255 nblkrows_total=nblkrows_tot, &
2256 col_blk_size=nmos)
2257
2258 ALLOCATE (subm_main(ndomains))
2259 CALL init_submatrices(subm_main)
2260 !!!TRIM ALLOCATE(subm_trimmer(ndomains))
2261
2262 CALL construct_submatrices(matrix_main, subm_main, &
2263 dpattern, map, node_of_domain, select_row_col)
2264
2265 !!!TRIM IF (matrix_trimmer_required) THEN
2266 !!!TRIM CALL construct_submatrices(matrix_trimmer,subm_trimmer,&
2267 !!!TRIM dpattern,map,node_of_domain,select_row)
2268 !!!TRIM ENDIF
2269
2270 IF (my_action == -1) THEN
2271 ! project out the local occupied space
2272 !tmp=MATMUL(subm_r(idomain)%mdata,Minv)
2273 !Minv=MATMUL(tmp,subm_main(idomain)%mdata)
2274 !subm_main(idomain)%mdata=subm_main(idomain)%mdata-&
2275 ! Minv-TRANSPOSE(Minv)+MATMUL(Minv,TRANSPOSE(tmp))
2276 ALLOCATE (subm_tmp(ndomains))
2277 ALLOCATE (subm_tmp2(ndomains))
2278 CALL init_submatrices(subm_tmp)
2279 CALL init_submatrices(subm_tmp2)
2280 CALL multiply_submatrices('N', 'N', 1.0_dp, subm_r_down, &
2281 subm_s_inv, 0.0_dp, subm_tmp)
2282 CALL multiply_submatrices('N', 'N', 1.0_dp, subm_tmp, &
2283 subm_main, 0.0_dp, subm_tmp2)
2284 CALL add_submatrices(1.0_dp, subm_main, -1.0_dp, subm_tmp2, 'N')
2285 CALL add_submatrices(1.0_dp, subm_main, -1.0_dp, subm_tmp2, 'T')
2286 CALL multiply_submatrices('N', 'T', 1.0_dp, subm_tmp2, &
2287 subm_tmp, 1.0_dp, subm_main)
2288 CALL release_submatrices(subm_tmp)
2289 CALL release_submatrices(subm_tmp2)
2290 DEALLOCATE (subm_tmp2)
2291 DEALLOCATE (subm_tmp)
2292 END IF
2293
2294 IF (my_skip_inversion) THEN
2295 CALL copy_submatrices(subm_main, preconditioner, .true.)
2296 ELSE
2297 ! loop over domains - perform inversion
2298 DO idomain = 1, ndomains
2299
2300 ! check if the submatrix exists
2301 IF (subm_main(idomain)%domain > 0) THEN
2302
2303 ! find sizes of MO submatrices
2304 IF (idomain == 1) THEN
2305 index1_start = 1
2306 ELSE
2307 index1_start = map%index1(idomain - 1)
2308 END IF
2309 index1_end = map%index1(idomain) - 1
2310
2311 n_domain_mos = 0
2312 DO row = index1_start, index1_end
2313 neighbor = map%pairs(row, 1)
2314 n_domain_mos = n_domain_mos + nmos(neighbor)
2315 END DO
2316
2317 naos = subm_main(idomain)%nrows
2318 !WRITE(*,*) "Domain, mo_self_and_neig, ao_domain: ", idomain, n_domain_mos, naos
2319
2320 ALLOCATE (minv(naos, naos))
2321
2322 !!!TRIM IF (my_use_trimmer) THEN
2323 !!!TRIM ! THIS IS SUPER EXPENSIVE (ELIMINATE)
2324 !!!TRIM ! trim the main matrix before inverting
2325 !!!TRIM ! assume that the trimmer columns are different (i.e. the main matrix is different for each MO)
2326 !!!TRIM allocate(tmp(naos,nmos(idomain)))
2327 !!!TRIM DO ii=1, nmos(idomain)
2328 !!!TRIM ! transform the main matrix using the trimmer for the current MO
2329 !!!TRIM DO jj=1, naos
2330 !!!TRIM DO kk=1, naos
2331 !!!TRIM Mstore(jj,kk)=sumb_main(idomain)%mdata(jj,kk)*&
2332 !!!TRIM subm_trimmer(idomain)%mdata(jj,ii)*&
2333 !!!TRIM subm_trimmer(idomain)%mdata(kk,ii)
2334 !!!TRIM ENDDO
2335 !!!TRIM ENDDO
2336 !!!TRIM ! invert the main matrix (exclude some eigenvalues, shift some)
2337 !!!TRIM CALL pseudo_invert_matrix(A=Mstore,Ainv=Minv,N=naos,method=1,&
2338 !!!TRIM !range1_thr=1.0E-9_dp,range2_thr=1.0E-9_dp,&
2339 !!!TRIM shift=1.0E-5_dp,&
2340 !!!TRIM range1=nmos(idomain),range2=nmos(idomain),&
2341 !!!TRIM
2342 !!!TRIM ! apply the inverted matrix
2343 !!!TRIM ! RZK-warning this is only possible when the preconditioner is applied
2344 !!!TRIM tmp(:,ii)=MATMUL(Minv,subm_in(idomain)%mdata(:,ii))
2345 !!!TRIM ENDDO
2346 !!!TRIM subm_out=MATMUL(tmp,sigma)
2347 !!!TRIM deallocate(tmp)
2348 !!!TRIM ELSE
2349
2350 IF (PRESENT(bad_modes_projector_down)) THEN
2351 ALLOCATE (proj_array(naos, naos))
2352 CALL pseudo_invert_matrix(a=subm_main(idomain)%mdata, ainv=minv, n=naos, method=1, &
2353 range1=nmos(idomain), range2=n_domain_mos, &
2354 range1_thr=eps_zero_eigenvalues, &
2355 bad_modes_projector_down=proj_array, &
2356 s_inv_half=subm_s_inv_half(idomain)%mdata, &
2357 s_half=subm_s_half(idomain)%mdata &
2358 )
2359 ELSE
2360 CALL pseudo_invert_matrix(a=subm_main(idomain)%mdata, ainv=minv, n=naos, method=1, &
2361 range1=nmos(idomain), range2=n_domain_mos)
2362 END IF
2363 !!!TRIM ENDIF
2364
2365 CALL copy_submatrices(subm_main(idomain), preconditioner(idomain), .false.)
2366 CALL copy_submatrix_data(minv, preconditioner(idomain))
2367 DEALLOCATE (minv)
2368
2369 IF (PRESENT(bad_modes_projector_down)) THEN
2370 CALL copy_submatrices(subm_main(idomain), bad_modes_projector_down(idomain), .false.)
2371 CALL copy_submatrix_data(proj_array, bad_modes_projector_down(idomain))
2372 DEALLOCATE (proj_array)
2373 END IF
2374
2375 END IF ! submatrix for the domain exists
2376
2377 END DO ! loop over domains
2378
2379 END IF
2380
2381 CALL release_submatrices(subm_main)
2382 DEALLOCATE (subm_main)
2383 !DEALLOCATE(subm_s)
2384 !DEALLOCATE(subm_r)
2385
2386 !IF (matrix_r_required) THEN
2387 ! CALL dbcsr_release(m_tmp_no_1)
2388 ! CALL dbcsr_release(m_tmp_no_2)
2389 ! CALL dbcsr_release(matrix_r)
2390 !ENDIF
2391
2392 CALL timestop(handle)
2393
2394 END SUBROUTINE construct_domain_preconditioner
2395
2396! **************************************************************************************************
2397!> \brief Constructs S^(+1/2) and S^(-1/2) submatrices for each domain
2398!> \param matrix_s ...
2399!> \param subm_s_sqrt ...
2400!> \param subm_s_sqrt_inv ...
2401!> \param dpattern ...
2402!> \param map ...
2403!> \param node_of_domain ...
2404!> \par History
2405!> 2013.03 created [Rustam Z. Khaliullin]
2406!> \author Rustam Z. Khaliullin
2407! **************************************************************************************************
2408 SUBROUTINE construct_domain_s_sqrt(matrix_s, subm_s_sqrt, subm_s_sqrt_inv, &
2409 dpattern, map, node_of_domain)
2410
2411 TYPE(dbcsr_type), INTENT(INOUT) :: matrix_s
2412 TYPE(domain_submatrix_type), DIMENSION(:), &
2413 INTENT(INOUT) :: subm_s_sqrt, subm_s_sqrt_inv
2414 TYPE(dbcsr_type), INTENT(IN) :: dpattern
2415 TYPE(domain_map_type), INTENT(IN) :: map
2416 INTEGER, DIMENSION(:), INTENT(IN) :: node_of_domain
2417
2418 CHARACTER(len=*), PARAMETER :: routinen = 'construct_domain_s_sqrt'
2419
2420 INTEGER :: handle, idomain, naos, ndomains
2421 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: ssqrt, ssqrtinv
2422 TYPE(domain_submatrix_type), ALLOCATABLE, &
2423 DIMENSION(:) :: subm_s
2424
2425 CALL timeset(routinen, handle)
2426
2427 CALL dbcsr_get_info(dpattern, nblkcols_total=ndomains)
2428 cpassert(SIZE(subm_s_sqrt) == ndomains)
2429 cpassert(SIZE(subm_s_sqrt_inv) == ndomains)
2430 ALLOCATE (subm_s(ndomains))
2431 CALL init_submatrices(subm_s)
2432
2433 CALL construct_submatrices(matrix_s, subm_s, &
2434 dpattern, map, node_of_domain, select_row_col)
2435
2436 ! loop over domains - perform inversion
2437 DO idomain = 1, ndomains
2438
2439 ! check if the submatrix exists
2440 IF (subm_s(idomain)%domain > 0) THEN
2441
2442 naos = subm_s(idomain)%nrows
2443
2444 ALLOCATE (ssqrt(naos, naos))
2445 ALLOCATE (ssqrtinv(naos, naos))
2446
2447 CALL matrix_sqrt(a=subm_s(idomain)%mdata, asqrt=ssqrt, asqrtinv=ssqrtinv, &
2448 n=naos)
2449
2450 CALL copy_submatrices(subm_s(idomain), subm_s_sqrt(idomain), .false.)
2451 CALL copy_submatrix_data(ssqrt, subm_s_sqrt(idomain))
2452
2453 CALL copy_submatrices(subm_s(idomain), subm_s_sqrt_inv(idomain), .false.)
2454 CALL copy_submatrix_data(ssqrtinv, subm_s_sqrt_inv(idomain))
2455
2456 DEALLOCATE (ssqrtinv)
2457 DEALLOCATE (ssqrt)
2458
2459 END IF ! submatrix for the domain exists
2460
2461 END DO ! loop over domains
2462
2463 CALL release_submatrices(subm_s)
2464 DEALLOCATE (subm_s)
2465
2466 CALL timestop(handle)
2467
2468 END SUBROUTINE construct_domain_s_sqrt
2469
2470! **************************************************************************************************
2471!> \brief Constructs S_inv block for each domain
2472!> \param matrix_s ...
2473!> \param subm_s_inv ...
2474!> \param dpattern ...
2475!> \param map ...
2476!> \param node_of_domain ...
2477!> \par History
2478!> 2013.02 created [Rustam Z. Khaliullin]
2479!> \author Rustam Z. Khaliullin
2480! **************************************************************************************************
2481 SUBROUTINE construct_domain_s_inv(matrix_s, subm_s_inv, dpattern, map, &
2482 node_of_domain)
2483 TYPE(dbcsr_type), INTENT(INOUT) :: matrix_s
2484 TYPE(domain_submatrix_type), DIMENSION(:), &
2485 INTENT(INOUT) :: subm_s_inv
2486 TYPE(dbcsr_type), INTENT(IN) :: dpattern
2487 TYPE(domain_map_type), INTENT(IN) :: map
2488 INTEGER, DIMENSION(:), INTENT(IN) :: node_of_domain
2489
2490 CHARACTER(len=*), PARAMETER :: routinen = 'construct_domain_s_inv'
2491
2492 INTEGER :: handle, idomain, naos, ndomains
2493 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: sinv
2494 TYPE(domain_submatrix_type), ALLOCATABLE, &
2495 DIMENSION(:) :: subm_s
2496
2497 CALL timeset(routinen, handle)
2498
2499 CALL dbcsr_get_info(dpattern, nblkcols_total=ndomains)
2500
2501 cpassert(SIZE(subm_s_inv) == ndomains)
2502 ALLOCATE (subm_s(ndomains))
2503 CALL init_submatrices(subm_s)
2504
2505 CALL construct_submatrices(matrix_s, subm_s, &
2506 dpattern, map, node_of_domain, select_row_col)
2507
2508 ! loop over domains - perform inversion
2509 DO idomain = 1, ndomains
2510
2511 ! check if the submatrix exists
2512 IF (subm_s(idomain)%domain > 0) THEN
2513
2514 naos = subm_s(idomain)%nrows
2515
2516 ALLOCATE (sinv(naos, naos))
2517
2518 CALL pseudo_invert_matrix(a=subm_s(idomain)%mdata, ainv=sinv, n=naos, &
2519 method=0)
2520
2521 CALL copy_submatrices(subm_s(idomain), subm_s_inv(idomain), .false.)
2522 CALL copy_submatrix_data(sinv, subm_s_inv(idomain))
2523
2524 DEALLOCATE (sinv)
2525
2526 END IF ! submatrix for the domain exists
2527
2528 END DO ! loop over domains
2529
2530 CALL release_submatrices(subm_s)
2531 DEALLOCATE (subm_s)
2532
2533 CALL timestop(handle)
2534
2535 END SUBROUTINE construct_domain_s_inv
2536
2537! **************************************************************************************************
2538!> \brief Constructs subblocks of the covariant-covariant projectors (i.e. DM without spin factor)
2539!> \param matrix_t ...
2540!> \param matrix_sigma_inv ...
2541!> \param matrix_s ...
2542!> \param subm_r_down ...
2543!> \param dpattern ...
2544!> \param map ...
2545!> \param node_of_domain ...
2546!> \param filter_eps ...
2547!> \par History
2548!> 2013.02 created [Rustam Z. Khaliullin]
2549!> \author Rustam Z. Khaliullin
2550! **************************************************************************************************
2551 SUBROUTINE construct_domain_r_down(matrix_t, matrix_sigma_inv, matrix_s, &
2552 subm_r_down, dpattern, map, node_of_domain, filter_eps)
2553
2554 TYPE(dbcsr_type), INTENT(IN) :: matrix_t, matrix_sigma_inv, matrix_s
2555 TYPE(domain_submatrix_type), DIMENSION(:), &
2556 INTENT(INOUT) :: subm_r_down
2557 TYPE(dbcsr_type), INTENT(IN) :: dpattern
2558 TYPE(domain_map_type), INTENT(IN) :: map
2559 INTEGER, DIMENSION(:), INTENT(IN) :: node_of_domain
2560 REAL(kind=dp) :: filter_eps
2561
2562 CHARACTER(len=*), PARAMETER :: routinen = 'construct_domain_r_down'
2563
2564 INTEGER :: handle, ndomains
2565 TYPE(dbcsr_type) :: m_tmp_no_1, m_tmp_no_2, matrix_r
2566
2567 CALL timeset(routinen, handle)
2568
2569 ! compute the density matrix in the COVARIANT representation
2570 CALL dbcsr_create(matrix_r, &
2571 template=matrix_s, &
2572 matrix_type=dbcsr_type_symmetric)
2573 CALL dbcsr_create(m_tmp_no_1, &
2574 template=matrix_t, &
2575 matrix_type=dbcsr_type_no_symmetry)
2576 CALL dbcsr_create(m_tmp_no_2, &
2577 template=matrix_t, &
2578 matrix_type=dbcsr_type_no_symmetry)
2579
2580 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_s, matrix_t, &
2581 0.0_dp, m_tmp_no_1, filter_eps=filter_eps)
2582 CALL dbcsr_multiply("N", "N", 1.0_dp, m_tmp_no_1, matrix_sigma_inv, &
2583 0.0_dp, m_tmp_no_2, filter_eps=filter_eps)
2584 CALL dbcsr_multiply("N", "T", 1.0_dp, m_tmp_no_2, m_tmp_no_1, &
2585 0.0_dp, matrix_r, filter_eps=filter_eps)
2586
2587 CALL dbcsr_release(m_tmp_no_1)
2588 CALL dbcsr_release(m_tmp_no_2)
2589
2590 CALL dbcsr_get_info(dpattern, nblkcols_total=ndomains)
2591 cpassert(SIZE(subm_r_down) == ndomains)
2592
2593 CALL construct_submatrices(matrix_r, subm_r_down, &
2594 dpattern, map, node_of_domain, select_row_col)
2595
2596 CALL dbcsr_release(matrix_r)
2597
2598 CALL timestop(handle)
2599
2600 END SUBROUTINE construct_domain_r_down
2601
2602! **************************************************************************************************
2603!> \brief Finds the square root of a matrix and its inverse
2604!> \param A ...
2605!> \param Asqrt ...
2606!> \param Asqrtinv ...
2607!> \param N ...
2608!> \par History
2609!> 2013.03 created [Rustam Z. Khaliullin]
2610!> \author Rustam Z. Khaliullin
2611! **************************************************************************************************
2612 SUBROUTINE matrix_sqrt(A, Asqrt, Asqrtinv, N)
2613
2614 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: a
2615 REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: asqrt, asqrtinv
2616 INTEGER, INTENT(IN) :: n
2617
2618 CHARACTER(len=*), PARAMETER :: routinen = 'matrix_sqrt'
2619
2620 INTEGER :: handle, info, jj, lwork, unit_nr
2621 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: eigenvalues, work
2622 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: test, testn
2623 TYPE(cp_logger_type), POINTER :: logger
2624
2625 CALL timeset(routinen, handle)
2626
2627 asqrtinv = a
2628 info = 0
2629
2630 ! get a useful unit_nr
2631 logger => cp_get_default_logger()
2632 IF (logger%para_env%is_source()) THEN
2633 unit_nr = cp_logger_get_default_unit_nr(logger, local=.true.)
2634 ELSE
2635 unit_nr = -1
2636 END IF
2637
2638 !CALL dpotrf('L', N, Ainv, N, INFO )
2639 !IF( INFO/=0 ) THEN
2640 ! CPErrorMessage(cp_failure_level,routineP,"DPOTRF failed")
2641 ! CPPrecondition(.FALSE.,cp_failure_level,routineP,failure)
2642 !END IF
2643 !CALL dpotri('L', N, Ainv, N, INFO )
2644 !IF( INFO/=0 ) THEN
2645 ! CPErrorMessage(cp_failure_level,routineP,"DPOTRI failed")
2646 ! CPPrecondition(.FALSE.,cp_failure_level,routineP,failure)
2647 !END IF
2648 !! complete the matrix
2649 !DO ii=1,N
2650 ! DO jj=ii+1,N
2651 ! Ainv(ii,jj)=Ainv(jj,ii)
2652 ! ENDDO
2653 ! !WRITE(*,'(100F13.9)') Ainv(ii,:)
2654 !ENDDO
2655
2656 ! diagonalize first
2657 ALLOCATE (eigenvalues(n))
2658 ! Query the optimal workspace for dsyev
2659 lwork = -1
2660 ALLOCATE (work(max(1, lwork)))
2661 CALL dsyev('V', 'L', n, asqrtinv, n, eigenvalues, work, lwork, info)
2662 lwork = int(work(1))
2663 DEALLOCATE (work)
2664 ! Allocate the workspace and solve the eigenproblem
2665 ALLOCATE (work(max(1, lwork)))
2666 CALL dsyev('V', 'L', n, asqrtinv, n, eigenvalues, work, lwork, info)
2667 IF (info /= 0) THEN
2668 IF (unit_nr > 0) WRITE (unit_nr, *) 'DSYEV ERROR MESSAGE: ', info
2669 cpabort("DSYEV failed")
2670 END IF
2671 DEALLOCATE (work)
2672
2673 ! take functions of eigenvalues and use eigenvectors to compute the matrix function
2674 ! first sqrt
2675 ALLOCATE (test(n, n))
2676 DO jj = 1, n
2677 test(jj, :) = asqrtinv(:, jj)*sqrt(eigenvalues(jj))
2678 END DO
2679 ALLOCATE (testn(n, n))
2680 testn(:, :) = matmul(asqrtinv, test)
2681 asqrt = testn
2682 ! now, sqrt_inv
2683 DO jj = 1, n
2684 test(jj, :) = asqrtinv(:, jj)/sqrt(eigenvalues(jj))
2685 END DO
2686 testn(:, :) = matmul(asqrtinv, test)
2687 asqrtinv = testn
2688 DEALLOCATE (test, testn)
2689
2690 DEALLOCATE (eigenvalues)
2691
2692!!! ! compute the error
2693!!! allocate(test(N,N))
2694!!! test=MATMUL(Ainv,A)
2695!!! DO ii=1,N
2696!!! test(ii,ii)=test(ii,ii)-1.0_dp
2697!!! ENDDO
2698!!! test_error=0.0_dp
2699!!! DO ii=1,N
2700!!! DO jj=1,N
2701!!! test_error=test_error+test(jj,ii)*test(jj,ii)
2702!!! ENDDO
2703!!! ENDDO
2704!!! WRITE(*,*) "Inversion error: ", SQRT(test_error)
2705!!! deallocate(test)
2706
2707 CALL timestop(handle)
2708
2709 END SUBROUTINE matrix_sqrt
2710
2711! **************************************************************************************************
2712!> \brief Inverts a matrix using a requested method
2713!> 0. Cholesky factorization
2714!> 1. Diagonalization
2715!> \param A ...
2716!> \param Ainv ...
2717!> \param N ...
2718!> \param method ...
2719!> \param range1 ...
2720!> \param range2 ...
2721!> \param range1_thr ...
2722!> \param shift ...
2723!> \param bad_modes_projector_down ...
2724!> \param s_inv_half ...
2725!> \param s_half ...
2726!> \par History
2727!> 2012.04 created [Rustam Z. Khaliullin]
2728!> \author Rustam Z. Khaliullin
2729! **************************************************************************************************
2730 SUBROUTINE pseudo_invert_matrix(A, Ainv, N, method, range1, range2, range1_thr, &
2731 shift, bad_modes_projector_down, s_inv_half, s_half)
2732
2733 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: a
2734 REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: ainv
2735 INTEGER, INTENT(IN) :: n, method
2736 INTEGER, INTENT(IN), OPTIONAL :: range1, range2
2737 REAL(kind=dp), INTENT(IN), OPTIONAL :: range1_thr, shift
2738 REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT), &
2739 OPTIONAL :: bad_modes_projector_down
2740 REAL(kind=dp), DIMENSION(:, :), INTENT(IN), &
2741 OPTIONAL :: s_inv_half, s_half
2742
2743 CHARACTER(len=*), PARAMETER :: routinen = 'pseudo_invert_matrix'
2744
2745 INTEGER :: handle, info, jj, lwork, range1_eiv, &
2746 range2_eiv, range3_eiv, unit_nr
2747 LOGICAL :: use_both, use_ranges_only, use_thr_only
2748 REAL(kind=dp) :: my_shift
2749 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: eigenvalues, work
2750 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: temp1, temp2, temp3, temp4
2751 TYPE(cp_logger_type), POINTER :: logger
2752
2753 CALL timeset(routinen, handle)
2754
2755 ! get a useful unit_nr
2756 logger => cp_get_default_logger()
2757 IF (logger%para_env%is_source()) THEN
2758 unit_nr = cp_logger_get_default_unit_nr(logger, local=.true.)
2759 ELSE
2760 unit_nr = -1
2761 END IF
2762
2763 IF (method == 1) THEN
2764
2765 IF ((PRESENT(range1) .AND. (.NOT. PRESENT(range2))) .OR. (PRESENT(range2) .AND. (.NOT. PRESENT(range1)))) THEN
2766 cpabort("range1 and range2 must be provided together")
2767 END IF
2768
2769 IF (PRESENT(range1) .AND. PRESENT(range1_thr)) THEN
2770 use_both = .true.
2771 use_thr_only = .false.
2772 use_ranges_only = .false.
2773 ELSE
2774 use_both = .false.
2775
2776 IF (PRESENT(range1)) THEN
2777 use_ranges_only = .true.
2778 ELSE
2779 use_ranges_only = .false.
2780 END IF
2781
2782 IF (PRESENT(range1_thr)) THEN
2783 use_thr_only = .true.
2784 ELSE
2785 use_thr_only = .false.
2786 END IF
2787
2788 END IF
2789
2790 IF ((PRESENT(s_half) .AND. (.NOT. PRESENT(s_inv_half))) .OR. (PRESENT(s_inv_half) .AND. (.NOT. PRESENT(s_half)))) THEN
2791 cpabort("Domain overlap matrix missing")
2792 END IF
2793 END IF
2794
2795 my_shift = 0.0_dp
2796 IF (PRESENT(shift)) THEN
2797 my_shift = shift
2798 END IF
2799
2800 ainv = a
2801 info = 0
2802
2803 SELECT CASE (method)
2804 CASE (0) ! Inversion via cholesky factorization
2805 CALL invmat_symm(ainv)
2806 CASE (1)
2807
2808 ! diagonalize first
2809 ALLOCATE (eigenvalues(n))
2810 ALLOCATE (temp1(n, n))
2811 ALLOCATE (temp4(n, n))
2812 IF (PRESENT(s_inv_half)) THEN
2813 CALL dsymm('L', 'U', n, n, 1.0_dp, s_inv_half, n, a, n, 0.0_dp, temp1, n)
2814 CALL dsymm('R', 'U', n, n, 1.0_dp, s_inv_half, n, temp1, n, 0.0_dp, ainv, n)
2815 END IF
2816 ! Query the optimal workspace for dsyev
2817 lwork = -1
2818 ALLOCATE (work(max(1, lwork)))
2819 CALL dsyev('V', 'L', n, ainv, n, eigenvalues, work, lwork, info)
2820
2821 lwork = int(work(1))
2822 DEALLOCATE (work)
2823 ! Allocate the workspace and solve the eigenproblem
2824 ALLOCATE (work(max(1, lwork)))
2825 CALL dsyev('V', 'L', n, ainv, n, eigenvalues, work, lwork, info)
2826
2827 IF (info /= 0) THEN
2828 IF (unit_nr > 0) WRITE (unit_nr, *) 'EIGENSYSTEM ERROR MESSAGE: ', info
2829 cpabort("Eigenproblem routine failed")
2830 END IF
2831 DEALLOCATE (work)
2832
2833 !WRITE(*,*) "EIGENVALS: "
2834 !WRITE(*,'(4F13.9)') eigenvalues(:)
2835
2836 ! invert eigenvalues and use eigenvectors to compute pseudo Ainv
2837 ! project out near-zero eigenvalue modes
2838 ALLOCATE (temp2(n, n))
2839 IF (PRESENT(bad_modes_projector_down)) ALLOCATE (temp3(n, n))
2840 temp2(1:n, 1:n) = ainv(1:n, 1:n)
2841
2842 range1_eiv = 0
2843 range2_eiv = 0
2844 range3_eiv = 0
2845
2846 IF (use_both) THEN
2847 DO jj = 1, n
2848 IF ((jj <= range2) .AND. (eigenvalues(jj) < range1_thr)) THEN
2849 temp1(jj, :) = temp2(:, jj)*0.0_dp
2850 IF (PRESENT(bad_modes_projector_down)) temp3(jj, :) = ainv(:, jj)*1.0_dp
2851 range1_eiv = range1_eiv + 1
2852 ELSE
2853 temp1(jj, :) = temp2(:, jj)/(eigenvalues(jj) + my_shift)
2854 IF (PRESENT(bad_modes_projector_down)) temp3(jj, :) = ainv(:, jj)*0.0_dp
2855 range2_eiv = range2_eiv + 1
2856 END IF
2857 END DO
2858 ELSE
2859 IF (use_ranges_only) THEN
2860 DO jj = 1, n
2861 IF (jj <= range1) THEN
2862 temp1(jj, :) = temp2(:, jj)*0.0_dp
2863 IF (PRESENT(bad_modes_projector_down)) temp3(jj, :) = ainv(:, jj)*1.0_dp
2864 range1_eiv = range1_eiv + 1
2865 ELSE IF (jj <= range2) THEN
2866 temp1(jj, :) = temp2(:, jj)*1.0_dp
2867 IF (PRESENT(bad_modes_projector_down)) temp3(jj, :) = ainv(:, jj)*1.0_dp
2868 range2_eiv = range2_eiv + 1
2869 ELSE
2870 temp1(jj, :) = temp2(:, jj)/(eigenvalues(jj) + my_shift)
2871 IF (PRESENT(bad_modes_projector_down)) temp3(jj, :) = ainv(:, jj)*0.0_dp
2872 range3_eiv = range3_eiv + 1
2873 END IF
2874 END DO
2875 ELSE IF (use_thr_only) THEN
2876 DO jj = 1, n
2877 IF (eigenvalues(jj) < range1_thr) THEN
2878 temp1(jj, :) = temp2(:, jj)*0.0_dp
2879 IF (PRESENT(bad_modes_projector_down)) temp3(jj, :) = ainv(:, jj)*1.0_dp
2880 range1_eiv = range1_eiv + 1
2881 ELSE
2882 temp1(jj, :) = temp2(:, jj)/(eigenvalues(jj) + my_shift)
2883 IF (PRESENT(bad_modes_projector_down)) temp3(jj, :) = ainv(:, jj)*0.0_dp
2884 range2_eiv = range2_eiv + 1
2885 END IF
2886 END DO
2887 ELSE ! no ranges, no thresholds
2888 cpabort("Invert using Cholesky. It would be faster.")
2889 END IF
2890 END IF
2891 !WRITE(*,*) ' EIV RANGES: ', range1_eiv, range2_eiv, range3_eiv
2892 IF (PRESENT(bad_modes_projector_down)) THEN
2893 IF (PRESENT(s_half)) THEN
2894 CALL dsymm('L', 'U', n, n, 1.0_dp, s_half, n, temp2, n, 0.0_dp, ainv, n)
2895 CALL dsymm('R', 'U', n, n, 1.0_dp, s_half, n, temp3, n, 0.0_dp, temp4, n)
2896 CALL dgemm('N', 'N', n, n, n, 1.0_dp, ainv, n, temp4, n, 0.0_dp, bad_modes_projector_down, n)
2897 ELSE
2898 CALL dgemm('N', 'N', n, n, n, 1.0_dp, temp2, n, temp3, n, 0.0_dp, bad_modes_projector_down, n)
2899 END IF
2900 END IF
2901
2902 IF (PRESENT(s_inv_half)) THEN
2903 CALL dsymm('L', 'U', n, n, 1.0_dp, s_inv_half, n, temp2, n, 0.0_dp, temp4, n)
2904 CALL dsymm('R', 'U', n, n, 1.0_dp, s_inv_half, n, temp1, n, 0.0_dp, temp2, n)
2905 CALL dgemm('N', 'N', n, n, n, 1.0_dp, temp4, n, temp2, n, 0.0_dp, ainv, n)
2906 ELSE
2907 CALL dgemm('N', 'N', n, n, n, 1.0_dp, temp2, n, temp1, n, 0.0_dp, ainv, n)
2908 END IF
2909 DEALLOCATE (temp1, temp2, temp4)
2910 IF (PRESENT(bad_modes_projector_down)) DEALLOCATE (temp3)
2911 DEALLOCATE (eigenvalues)
2912
2913 CASE DEFAULT
2914
2915 cpabort("Illegal method selected for matrix inversion")
2916
2917 END SELECT
2918
2919 !! compute the inversion error
2920 !allocate(temp1(N,N))
2921 !temp1=MATMUL(Ainv,A)
2922 !DO ii=1,N
2923 ! temp1(ii,ii)=temp1(ii,ii)-1.0_dp
2924 !ENDDO
2925 !temp1_error=0.0_dp
2926 !DO ii=1,N
2927 ! DO jj=1,N
2928 ! temp1_error=temp1_error+temp1(jj,ii)*temp1(jj,ii)
2929 ! ENDDO
2930 !ENDDO
2931 !WRITE(*,*) "Inversion error: ", SQRT(temp1_error)
2932 !deallocate(temp1)
2933
2934 CALL timestop(handle)
2935
2936 END SUBROUTINE pseudo_invert_matrix
2937
2938! **************************************************************************************************
2939!> \brief Find matrix power using diagonalization
2940!> \param A ...
2941!> \param Apow ...
2942!> \param power ...
2943!> \param N ...
2944!> \param range1 ...
2945!> \param range1_thr ...
2946!> \param shift ...
2947!> \par History
2948!> 2012.04 created [Rustam Z. Khaliullin]
2949!> \author Rustam Z. Khaliullin
2950! **************************************************************************************************
2951 SUBROUTINE pseudo_matrix_power(A, Apow, power, N, range1, range1_thr, shift)
2952
2953 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: a
2954 REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: apow
2955 REAL(kind=dp), INTENT(IN) :: power
2956 INTEGER, INTENT(IN) :: n
2957 INTEGER, INTENT(IN), OPTIONAL :: range1
2958 REAL(kind=dp), INTENT(IN), OPTIONAL :: range1_thr, shift
2959
2960 CHARACTER(len=*), PARAMETER :: routinen = 'pseudo_matrix_power'
2961
2962 INTEGER :: handle, info, jj, lwork, range1_eiv, &
2963 range2_eiv, unit_nr
2964 LOGICAL :: use_both, use_ranges, use_thr
2965 REAL(kind=dp) :: my_shift
2966 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: eigenvalues, work
2967 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: temp1, temp2
2968 TYPE(cp_logger_type), POINTER :: logger
2969
2970 CALL timeset(routinen, handle)
2971
2972 ! get a useful unit_nr
2973 logger => cp_get_default_logger()
2974 IF (logger%para_env%is_source()) THEN
2975 unit_nr = cp_logger_get_default_unit_nr(logger, local=.true.)
2976 ELSE
2977 unit_nr = -1
2978 END IF
2979
2980 IF (PRESENT(range1) .AND. PRESENT(range1_thr)) THEN
2981 use_both = .true.
2982 ELSE
2983 use_both = .false.
2984 IF (PRESENT(range1)) THEN
2985 use_ranges = .true.
2986 ELSE
2987 use_ranges = .false.
2988 IF (PRESENT(range1_thr)) THEN
2989 use_thr = .true.
2990 ELSE
2991 use_thr = .false.
2992 END IF
2993 END IF
2994 END IF
2995
2996 my_shift = 0.0_dp
2997 IF (PRESENT(shift)) THEN
2998 my_shift = shift
2999 END IF
3000
3001 apow = a
3002 info = 0
3003
3004 ! diagonalize first
3005 ALLOCATE (eigenvalues(n))
3006 ALLOCATE (temp1(n, n))
3007
3008 ! Query the optimal workspace for dsyev
3009 lwork = -1
3010 ALLOCATE (work(max(1, lwork)))
3011 CALL dsyev('V', 'L', n, apow, n, eigenvalues, work, lwork, info)
3012
3013 lwork = int(work(1))
3014 DEALLOCATE (work)
3015 ! Allocate the workspace and solve the eigenproblem
3016 ALLOCATE (work(max(1, lwork)))
3017 CALL dsyev('V', 'L', n, apow, n, eigenvalues, work, lwork, info)
3018
3019 IF (info /= 0) THEN
3020 IF (unit_nr > 0) WRITE (unit_nr, *) 'EIGENSYSTEM ERROR MESSAGE: ', info
3021 cpabort("Eigenproblem routine failed")
3022 END IF
3023 DEALLOCATE (work)
3024
3025 !WRITE(*,*) "EIGENVALS: "
3026 !WRITE(*,'(4F13.9)') eigenvalues(:)
3027
3028 ! invert eigenvalues and use eigenvectors to compute pseudo Ainv
3029 ! project out near-zero eigenvalue modes
3030 ALLOCATE (temp2(n, n))
3031
3032 temp2(1:n, 1:n) = apow(1:n, 1:n)
3033
3034 range1_eiv = 0
3035 range2_eiv = 0
3036
3037 IF (use_both) THEN
3038 DO jj = 1, n
3039 IF ((jj <= range1) .AND. (eigenvalues(jj) < range1_thr)) THEN
3040 temp1(jj, :) = temp2(:, jj)*0.0_dp
3041 range1_eiv = range1_eiv + 1
3042 ELSE
3043 temp1(jj, :) = temp2(:, jj)*((eigenvalues(jj) + my_shift)**power)
3044 END IF
3045 END DO
3046 ELSE
3047 IF (use_ranges) THEN
3048 DO jj = 1, n
3049 IF (jj <= range1) THEN
3050 temp1(jj, :) = temp2(:, jj)*0.0_dp
3051 range1_eiv = range1_eiv + 1
3052 ELSE
3053 temp1(jj, :) = temp2(:, jj)*((eigenvalues(jj) + my_shift)**power)
3054 END IF
3055 END DO
3056 ELSE
3057 IF (use_thr) THEN
3058 DO jj = 1, n
3059 IF (eigenvalues(jj) < range1_thr) THEN
3060 temp1(jj, :) = temp2(:, jj)*0.0_dp
3061
3062 range1_eiv = range1_eiv + 1
3063 ELSE
3064 temp1(jj, :) = temp2(:, jj)*((eigenvalues(jj) + my_shift)**power)
3065 END IF
3066 END DO
3067 ELSE
3068 DO jj = 1, n
3069 temp1(jj, :) = temp2(:, jj)*((eigenvalues(jj) + my_shift)**power)
3070 END DO
3071 END IF
3072 END IF
3073 END IF
3074 !WRITE(*,*) ' EIV RANGES: ', range1_eiv, range2_eiv, range3_eiv
3075 apow = matmul(temp2, temp1)
3076 DEALLOCATE (temp1, temp2)
3077 DEALLOCATE (eigenvalues)
3078
3079 CALL timestop(handle)
3080
3081 END SUBROUTINE pseudo_matrix_power
3082
3083! **************************************************************************************************
3084!> \brief Load balancing of the submatrix computations
3085!> \param almo_scf_env ...
3086!> \par History
3087!> 2013.02 created [Rustam Z. Khaliullin]
3088!> \author Rustam Z. Khaliullin
3089! **************************************************************************************************
3090 SUBROUTINE distribute_domains(almo_scf_env)
3091
3092 TYPE(almo_scf_env_type), INTENT(INOUT) :: almo_scf_env
3093
3094 CHARACTER(len=*), PARAMETER :: routinen = 'distribute_domains'
3095
3096 INTEGER :: handle, idomain, least_loaded, nao, &
3097 ncpus, ndomains
3098 INTEGER, ALLOCATABLE, DIMENSION(:) :: index0
3099 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: cpu_load, domain_load
3100 TYPE(dbcsr_distribution_type) :: dist
3101
3102 CALL timeset(routinen, handle)
3103
3104 ndomains = almo_scf_env%ndomains
3105 CALL dbcsr_get_info(almo_scf_env%matrix_s(1), distribution=dist)
3106 CALL dbcsr_distribution_get(dist, numnodes=ncpus)
3107
3108 ALLOCATE (domain_load(ndomains))
3109 DO idomain = 1, ndomains
3110 nao = almo_scf_env%nbasis_of_domain(idomain)
3111 domain_load(idomain) = (nao*nao*nao)*1.0_dp
3112 END DO
3113
3114 ALLOCATE (index0(ndomains))
3115
3116 CALL sort(domain_load, ndomains, index0)
3117
3118 ALLOCATE (cpu_load(ncpus))
3119 cpu_load(:) = 0.0_dp
3120
3121 DO idomain = 1, ndomains
3122 least_loaded = minloc(cpu_load, 1)
3123 cpu_load(least_loaded) = cpu_load(least_loaded) + domain_load(idomain)
3124 almo_scf_env%cpu_of_domain(index0(idomain)) = least_loaded - 1
3125 END DO
3126
3127 DEALLOCATE (cpu_load)
3128 DEALLOCATE (index0)
3129 DEALLOCATE (domain_load)
3130
3131 CALL timestop(handle)
3132
3133 END SUBROUTINE distribute_domains
3134
3135! **************************************************************************************************
3136!> \brief Tests construction and release of domain submatrices
3137!> \param matrix_no ...
3138!> \param dpattern ...
3139!> \param map ...
3140!> \param node_of_domain ...
3141!> \par History
3142!> 2013.01 created [Rustam Z. Khaliullin]
3143!> \author Rustam Z. Khaliullin
3144! **************************************************************************************************
3145 SUBROUTINE construct_test(matrix_no, dpattern, map, node_of_domain)
3146
3147 TYPE(dbcsr_type), INTENT(INOUT) :: matrix_no, dpattern
3148 TYPE(domain_map_type), INTENT(IN) :: map
3149 INTEGER, DIMENSION(:), INTENT(IN) :: node_of_domain
3150
3151 CHARACTER(len=*), PARAMETER :: routinen = 'construct_test'
3152
3153 INTEGER :: handle, ndomains
3154 TYPE(dbcsr_type) :: copy1
3155 TYPE(domain_submatrix_type), ALLOCATABLE, &
3156 DIMENSION(:) :: subm_nn, subm_no
3157 TYPE(mp_comm_type) :: group
3158
3159 CALL timeset(routinen, handle)
3160
3161 CALL dbcsr_get_info(dpattern, group=group, nblkcols_total=ndomains)
3162
3163 ALLOCATE (subm_no(ndomains), subm_nn(ndomains))
3164 CALL init_submatrices(subm_no)
3165 CALL init_submatrices(subm_nn)
3166
3167 !CALL dbcsr_print(matrix_nn)
3168 !CALL construct_submatrices(matrix_nn,subm_nn,dpattern,map,select_row_col)
3169 !CALL print_submatrices(subm_nn,Group)
3170
3171 !CALL dbcsr_print(matrix_no)
3172 CALL construct_submatrices(matrix_no, subm_no, dpattern, map, node_of_domain, select_row)
3173 CALL print_submatrices(subm_no, group)
3174
3175 CALL dbcsr_create(copy1, template=matrix_no)
3176 CALL dbcsr_copy(copy1, matrix_no)
3177 CALL dbcsr_print(copy1)
3178 CALL construct_dbcsr_from_submatrices(copy1, subm_no, dpattern)
3179 CALL dbcsr_print(copy1)
3180 CALL dbcsr_release(copy1)
3181
3182 CALL release_submatrices(subm_no)
3183 CALL release_submatrices(subm_nn)
3184 DEALLOCATE (subm_no, subm_nn)
3185
3186 CALL timestop(handle)
3187
3188 END SUBROUTINE construct_test
3189
3190! **************************************************************************************************
3191!> \brief create the initial guess for XALMOs
3192!> \param m_guess ...
3193!> \param m_t_in ...
3194!> \param m_t0 ...
3195!> \param m_quench_t ...
3196!> \param m_overlap ...
3197!> \param m_sigma_tmpl ...
3198!> \param nspins ...
3199!> \param xalmo_history ...
3200!> \param assume_t0_q0x ...
3201!> \param optimize_theta ...
3202!> \param envelope_amplitude ...
3203!> \param eps_filter ...
3204!> \param order_lanczos ...
3205!> \param eps_lanczos ...
3206!> \param max_iter_lanczos ...
3207!> \param nocc_of_domain ...
3208!> \par History
3209!> 2016.11 created [Rustam Z Khaliullin]
3210!> \author Rustam Z Khaliullin
3211! **************************************************************************************************
3212 SUBROUTINE xalmo_initial_guess(m_guess, m_t_in, m_t0, m_quench_t, &
3213 m_overlap, m_sigma_tmpl, nspins, xalmo_history, assume_t0_q0x, &
3214 optimize_theta, envelope_amplitude, eps_filter, order_lanczos, eps_lanczos, &
3215 max_iter_lanczos, nocc_of_domain)
3216
3217 TYPE(dbcsr_type), DIMENSION(:), INTENT(INOUT) :: m_guess
3218 TYPE(dbcsr_type), DIMENSION(:), INTENT(IN) :: m_t_in, m_t0, m_quench_t
3219 TYPE(dbcsr_type), INTENT(IN) :: m_overlap
3220 TYPE(dbcsr_type), DIMENSION(:), INTENT(IN) :: m_sigma_tmpl
3221 INTEGER, INTENT(IN) :: nspins
3222 TYPE(almo_scf_history_type), INTENT(IN) :: xalmo_history
3223 LOGICAL, INTENT(IN) :: assume_t0_q0x, optimize_theta
3224 REAL(kind=dp), INTENT(IN) :: envelope_amplitude, eps_filter
3225 INTEGER, INTENT(IN) :: order_lanczos
3226 REAL(kind=dp), INTENT(IN) :: eps_lanczos
3227 INTEGER, INTENT(IN) :: max_iter_lanczos
3228 INTEGER, DIMENSION(:, :), INTENT(IN) :: nocc_of_domain
3229
3230 CHARACTER(len=*), PARAMETER :: routinen = 'xalmo_initial_guess'
3231
3232 INTEGER :: handle, iaspc, ispin, istore, naspc, &
3233 unit_nr
3234 LOGICAL :: aspc_guess
3235 REAL(kind=dp) :: alpha
3236 TYPE(cp_logger_type), POINTER :: logger
3237 TYPE(dbcsr_type) :: m_extrapolated, m_sigma_tmp
3238
3239 CALL timeset(routinen, handle)
3240
3241 ! get a useful output_unit
3242 logger => cp_get_default_logger()
3243 IF (logger%para_env%is_source()) THEN
3244 unit_nr = cp_logger_get_default_unit_nr(logger, local=.true.)
3245 ELSE
3246 unit_nr = -1
3247 END IF
3248
3249 IF (optimize_theta) THEN
3250 cpwarn("unused option")
3251 ! just not to trigger unused variable
3252 alpha = envelope_amplitude
3253 END IF
3254
3255 ! if extrapolation order is zero then the standard guess is used
3256 ! ... the number of stored history points will remain zero if extrapolation order is zero
3257 IF (xalmo_history%istore == 0) THEN
3258 aspc_guess = .false.
3259 ELSE
3260 aspc_guess = .true.
3261 END IF
3262
3263 ! create initial guess
3264 IF (.NOT. aspc_guess) THEN
3265
3266 DO ispin = 1, nspins
3267
3268 ! zero initial guess for the delocalization amplitudes
3269 ! or the supplied guess for orbitals
3270 IF (assume_t0_q0x) THEN
3271 CALL dbcsr_set(m_guess(ispin), 0.0_dp)
3272 ELSE
3273 ! copy coefficients to m_guess
3274 CALL dbcsr_copy(m_guess(ispin), m_t_in(ispin))
3275 END IF
3276
3277 END DO !ispins
3278
3279 ELSE !aspc_guess
3280
3281 CALL cite_reference(kolafa2004)
3282 CALL cite_reference(kuhne2007)
3283
3284 naspc = min(xalmo_history%istore, xalmo_history%nstore)
3285 IF (unit_nr > 0) THEN
3286 WRITE (unit_nr, fmt="(/,T2,A,/,/,T3,A,I0)") &
3287 "Parameters for the always stable predictor-corrector (ASPC) method:", &
3288 "ASPC order: ", naspc
3289 END IF
3290
3291 DO ispin = 1, nspins
3292
3293 CALL dbcsr_create(m_extrapolated, &
3294 template=m_quench_t(ispin), matrix_type=dbcsr_type_no_symmetry)
3295 CALL dbcsr_create(m_sigma_tmp, &
3296 template=m_sigma_tmpl(ispin), matrix_type=dbcsr_type_no_symmetry)
3297
3298 ! set to zero before accumulation
3299 CALL dbcsr_set(m_guess(ispin), 0.0_dp)
3300
3301 ! extrapolation
3302 DO iaspc = 1, naspc
3303
3304 istore = mod(xalmo_history%istore - iaspc, xalmo_history%nstore) + 1
3305 alpha = (-1.0_dp)**(iaspc + 1)*real(iaspc, kind=dp)* &
3306 binomial(2*naspc, naspc - iaspc)/binomial(2*naspc - 2, naspc - 1)
3307 IF (unit_nr > 0) THEN
3308 WRITE (unit_nr, fmt="(T3,A2,I0,A4,F10.6)") &
3309 "B(", iaspc, ") = ", alpha
3310 END IF
3311
3312 ! m_extrapolated - initialize the correct sparsity pattern
3313 ! it must be kept throughout extrapolation
3314 CALL dbcsr_copy(m_extrapolated, m_quench_t(ispin))
3315
3316 ! project t0 onto the previous DMs
3317 ! note that t0 is projected instead of any other matrix (e.g.
3318 ! t_SCF from the prev step or random t)
3319 ! this is done to keep orbitals phase (i.e. sign) the same as in
3320 ! t0. if this is not done then subtracting t0 on the next step
3321 ! will produce a terrible guess and extrapolation will fail
3322 CALL dbcsr_multiply("N", "N", 1.0_dp, &
3323 xalmo_history%matrix_p_up_down(ispin, istore), &
3324 m_t0(ispin), &
3325 0.0_dp, m_extrapolated, &
3326 retain_sparsity=.true.)
3327 ! normalize MOs
3328 CALL orthogonalize_mos(ket=m_extrapolated, &
3329 overlap=m_sigma_tmp, &
3330 metric=m_overlap, &
3331 retain_locality=.true., &
3332 only_normalize=.false., &
3333 nocc_of_domain=nocc_of_domain(:, ispin), &
3334 eps_filter=eps_filter, &
3335 order_lanczos=order_lanczos, &
3336 eps_lanczos=eps_lanczos, &
3337 max_iter_lanczos=max_iter_lanczos)
3338
3339 ! now accumulate. correct sparsity is ensured
3340 CALL dbcsr_add(m_guess(ispin), m_extrapolated, &
3341 1.0_dp, (1.0_dp*alpha)/naspc)
3342
3343 END DO !iaspc
3344
3345 CALL dbcsr_release(m_extrapolated)
3346
3347 ! normalize MOs
3348 CALL orthogonalize_mos(ket=m_guess(ispin), &
3349 overlap=m_sigma_tmp, &
3350 metric=m_overlap, &
3351 retain_locality=.true., &
3352 only_normalize=.false., &
3353 nocc_of_domain=nocc_of_domain(:, ispin), &
3354 eps_filter=eps_filter, &
3355 order_lanczos=order_lanczos, &
3356 eps_lanczos=eps_lanczos, &
3357 max_iter_lanczos=max_iter_lanczos)
3358
3359 CALL dbcsr_release(m_sigma_tmp)
3360
3361 ! project the t0 space out from the extrapolated state
3362 ! this can be done outside this subroutine
3363 IF (assume_t0_q0x) THEN
3364 CALL dbcsr_add(m_guess(ispin), m_t0(ispin), &
3365 1.0_dp, -1.0_dp)
3366 END IF !assume_t0_q0x
3367
3368 END DO !ispin
3369
3370 END IF !aspc_guess?
3371
3372 CALL timestop(handle)
3373
3374 END SUBROUTINE xalmo_initial_guess
3375
3376END MODULE almo_scf_methods
3377
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:218
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:2170
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:3216
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:3146
almo_scf_methods::distribute_domains
subroutine, public distribute_domains(almo_scf_env)
Load balancing of the submatrix computations.
Definition almo_scf_methods.F:3091
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:1947
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:2044
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:2553
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:2483
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:1791
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:2410
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:1232
cp_dbcsr_api::dbcsr_iterator_next_block
subroutine, public dbcsr_iterator_next_block(iterator, row, column, block, block_number_argument_has_been_removed, row_size, col_size, row_offset, col_offset)
...
Definition cp_dbcsr_api.F:969
cp_dbcsr_api::dbcsr_iterator_blocks_left
logical function, public dbcsr_iterator_blocks_left(iterator)
...
Definition cp_dbcsr_api.F:943
cp_dbcsr_api::dbcsr_iterator_stop
subroutine, public dbcsr_iterator_stop(iterator)
...
Definition cp_dbcsr_api.F:1040
cp_dbcsr_api::dbcsr_copy
subroutine, public dbcsr_copy(matrix_b, matrix_a, name, keep_sparsity, keep_imaginary)
...
Definition cp_dbcsr_api.F:370
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:1086
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:807
cp_dbcsr_api::dbcsr_work_create
subroutine, public dbcsr_work_create(matrix, nblks_guess, sizedata_guess, n, work_mutable)
...
Definition cp_dbcsr_api.F:1297
cp_dbcsr_api::dbcsr_filter
subroutine, public dbcsr_filter(matrix, eps)
...
Definition cp_dbcsr_api.F:657
cp_dbcsr_api::dbcsr_finalize
subroutine, public dbcsr_finalize(matrix)
...
Definition cp_dbcsr_api.F:672
cp_dbcsr_api::dbcsr_iterator_start
subroutine, public dbcsr_iterator_start(iterator, matrix, shared, dynamic, dynamic_byrows)
...
Definition cp_dbcsr_api.F:1002
cp_dbcsr_api::dbcsr_set
subroutine, public dbcsr_set(matrix, alpha)
...
Definition cp_dbcsr_api.F:1194
cp_dbcsr_api::dbcsr_release
subroutine, public dbcsr_release(matrix)
...
Definition cp_dbcsr_api.F:1132
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:1023
cp_dbcsr_api::dbcsr_put_block
subroutine, public dbcsr_put_block(matrix, row, col, block, summation)
...
Definition cp_dbcsr_api.F:1115
cp_dbcsr_api::dbcsr_add
subroutine, public dbcsr_add(matrix_a, matrix_b, alpha_scalar, beta_scalar)
...
Definition cp_dbcsr_api.F:253
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:571
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:567
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition cp_log_handling.F:234
domain_submatrix_methods
Subroutines to handle submatrices.
Definition domain_submatrix_methods.F:14
domain_submatrix_methods::copy_submatrix_data
subroutine, public copy_submatrix_data(array, copy)
...
Definition domain_submatrix_methods.F:258
domain_submatrix_methods::print_submatrices
subroutine, public print_submatrices(submatrices, mpgroup)
...
Definition domain_submatrix_methods.F:1589
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:1332
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:17
input_constants::smear_fermi_dirac
integer, parameter, public smear_fermi_dirac
Definition input_constants.F:910
input_constants::spd_inversion_dense_cholesky
integer, parameter, public spd_inversion_dense_cholesky
Definition input_constants.F:1066
input_constants::almo_domain_layout_molecular
integer, parameter, public almo_domain_layout_molecular
Definition input_constants.F:1001
input_constants::almo_scf_diag
integer, parameter, public almo_scf_diag
Definition input_constants.F:1022
input_constants::almo_mat_distr_atomic
integer, parameter, public almo_mat_distr_atomic
Definition input_constants.F:1005
input_constants::spd_inversion_ls_taylor
integer, parameter, public spd_inversion_ls_taylor
Definition input_constants.F:1066
input_constants::spd_inversion_ls_hotelling
integer, parameter, public spd_inversion_ls_hotelling
Definition input_constants.F:1066
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:1626
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:312
util
All kind of helpful little routines.
Definition util.F:14
almo_scf_types::almo_scf_env_type
Definition almo_scf_types.F:104
almo_scf_types::almo_scf_history_type
Definition almo_scf_types.F:96
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:155
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:743
orbital
Orbital angular momentum.
Definition grid_common.h:133