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