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