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