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