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