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