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