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