(git:3f0d10b)
Loading...
Searching...
No Matches
preconditioner.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \brief computes preconditioners, and implements methods to apply them
10!> currently used in qs_ot
11!> \par History
12!> - [UB] 2009-05-13 Adding stable approximate inverse (full and sparse)
13!> \author Joost VandeVondele (09.2002)
14! **************************************************************************************************
23 USE cp_cfm_types, ONLY: cp_cfm_create,&
29 USE cp_dbcsr_api, ONLY: dbcsr_get_info,&
34 dbcsr_type_no_symmetry
48 USE cp_fm_types, ONLY: cp_fm_create,&
58 USE input_constants, ONLY: &
68 USE kinds, ONLY: default_string_length,&
69 dp
85 USE preconditioner_makes, ONLY: &
101 get_mo_set,&
104
105!$ USE OMP_LIB, ONLY: omp_get_max_threads
106#include "./base/base_uses.f90"
107
108 IMPLICIT NONE
109
110 PRIVATE
111
112 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'preconditioner'
113 REAL(KIND=dp), PARAMETER, PRIVATE :: fermi_low_rank_gap_factor = 4.0_dp, &
114 fermi_low_rank_degeneracy_tolerance = 1.0e-8_dp, &
115 fermi_low_rank_min_window = 1.0_dp
116
126
127 TYPE full_all_work_type
128 INTEGER :: info = 0
129 TYPE(cp_cfm_type) :: c, h, s, a, b
130 TYPE(cp_cfm_diag_workspace_type) :: solver
131 END TYPE full_all_work_type
132
133! The public interface for apply preconditioner, the routines can be found in preconditioner_apply.F
135 MODULE PROCEDURE apply_preconditioner_cfm_complex
136 MODULE PROCEDURE apply_preconditioner_dbcsr
138 MODULE PROCEDURE apply_preconditioner_fm
139 END INTERFACE
140
141! **************************************************************************************************
142
143CONTAINS
144
145! **************************************************************************************************
146
147! **************************************************************************************************
148!> \brief Build the block-circulant inverse of a replicated-cell SPD operator.
149!> \param preconditioner_env ...
150!> \param operator_matrix ...
151!> \param operator_name ...
152!> \param rhs_count number of occupied orbital columns to precondition
153! **************************************************************************************************
154 SUBROUTINE prepare_lattice_operator(preconditioner_env, operator_matrix, operator_name, rhs_count)
155
156 TYPE(preconditioner_type) :: preconditioner_env
157 TYPE(dbcsr_type), POINTER :: operator_matrix
158 CHARACTER(LEN=*), INTENT(IN) :: operator_name
159 INTEGER, INTENT(IN) :: rhs_count
160
161 INTEGER :: alpha, beta, block_size, col_cell, delta_index, i, info, j, local_info, n, ncell, &
162 ncol_local, nrow_local, output_unit, row_cell
163 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: index_to_cell
164 INTEGER, DIMENSION(3) :: delta
165 INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
166 LOGICAL :: used_fft
167 REAL(KIND=dp) :: denominator, difference, numerator, value
168 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: dense_operator, xkp
169 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: blocks
170 REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :), &
171 POINTER :: local_data
172 TYPE(cp_fm_struct_type), POINTER :: fm_struct
173 TYPE(cp_fm_type) :: operator_fm, work
174 TYPE(cp_logger_type), POINTER :: logger
175
176 cpassert(ASSOCIATED(operator_matrix))
177 cpassert(preconditioner_env%lattice_fft_mode /= ot_lattice_fft_off)
178 preconditioner_env%lattice_fft_active = .false.
179 preconditioner_env%lattice_projection_error = 0.0_dp
180 preconditioner_env%lattice_corrected_projection_error = 0.0_dp
181 preconditioner_env%lattice_cost_ratio = 0.0_dp
182 preconditioner_env%lattice_storage_ratio = 0.0_dp
183
184 IF (any(preconditioner_env%lattice_dims <= 0)) cpabort("Lattice dimensions must be positive")
185
186 IF (preconditioner_env%lattice_fft_mode == ot_lattice_fft_auto .AND. &
187 product(preconditioner_env%lattice_dims) < lattice_fft_auto_min_cells) THEN
188 IF (preconditioner_env%para_env%is_source()) THEN
189 logger => cp_get_default_logger()
190 output_unit = cp_logger_get_default_io_unit(logger)
191 IF (output_unit > 0) THEN
192 WRITE (output_unit, '(T2,A,A,A,I0,A,I0)') &
193 "OT| Lattice FFT AUTO kept the conventional ", trim(operator_name), &
194 " solver: ", product(preconditioner_env%lattice_dims), &
195 " cells are fewer than ", lattice_fft_auto_min_cells
196 END IF
197 END IF
198 RETURN
199 END IF
200
201 NULLIFY (fm_struct)
202 CALL dbcsr_get_info(operator_matrix, nfullrows_total=n)
203 ncell = product(preconditioner_env%lattice_dims)
204 IF (modulo(n, ncell) /= 0) cpabort("AO dimension does not match replicated cells")
205 block_size = n/ncell
206 preconditioner_env%lattice_block_size = block_size
207 preconditioner_env%lattice_cost_ratio = lattice_fft_cost_ratio( &
208 preconditioner_env%lattice_dims, block_size, rhs_count, &
209 preconditioner_env%para_env%num_pe) + &
211 preconditioner_env%lattice_dims, &
212 preconditioner_env%lattice_local_cells)
213 preconditioner_env%lattice_storage_ratio = lattice_fft_storage_ratio( &
214 preconditioner_env%lattice_dims, 1) + &
216 preconditioner_env%lattice_dims, 1, &
217 preconditioner_env%lattice_local_cells)
218 IF (preconditioner_env%lattice_fft_mode == ot_lattice_fft_auto .AND. &
219 (preconditioner_env%lattice_cost_ratio > lattice_fft_auto_max_cost_ratio .OR. &
220 preconditioner_env%lattice_storage_ratio > lattice_fft_auto_max_storage_ratio)) THEN
221 IF (preconditioner_env%para_env%is_source()) THEN
222 logger => cp_get_default_logger()
223 output_unit = cp_logger_get_default_io_unit(logger)
224 IF (output_unit > 0) THEN
225 WRITE (output_unit, '(T2,A,A,A,F8.3,A,F8.3,A)') &
226 "OT| Lattice FFT AUTO kept the conventional ", trim(operator_name), &
227 " solver: estimated work ratio ", preconditioner_env%lattice_cost_ratio, &
228 " storage ratio ", preconditioner_env%lattice_storage_ratio, " selected: F"
229 END IF
230 END IF
231 RETURN
232 END IF
233 CALL cp_fm_struct_create(fm_struct, nrow_global=n, ncol_global=n, &
234 context=preconditioner_env%ctxt, &
235 para_env=preconditioner_env%para_env)
236 CALL cp_fm_create(operator_fm, fm_struct, name="lattice preconditioner operator")
237 CALL cp_fm_struct_release(fm_struct)
238 CALL cp_fm_set_all(operator_fm, alpha=0.0_dp, beta=0.0_dp)
239 CALL copy_dbcsr_to_fm(operator_matrix, operator_fm)
240 CALL cp_fm_create(work, operator_fm%matrix_struct, name="lattice operator work")
241 CALL cp_fm_uplo_to_full(operator_fm, work)
242 CALL cp_fm_release(work)
243 CALL cp_fm_get_info(operator_fm, nrow_local=nrow_local, ncol_local=ncol_local, &
244 row_indices=row_indices, col_indices=col_indices, local_data=local_data)
245 ALLOCATE (blocks(block_size, block_size, ncell), source=0.0_dp)
246 ALLOCATE (index_to_cell(3, ncell), xkp(3, ncell))
247 CALL lattice_grid(preconditioner_env%lattice_dims, index_to_cell, xkp)
248 IF (preconditioner_env%lattice_local_cells > 0) THEN
249 ALLOCATE (dense_operator(n, n), source=0.0_dp)
250 END IF
251 denominator = 0.0_dp
252 DO j = 1, ncol_local
253 DO i = 1, nrow_local
254 IF (row_indices(i) <= col_indices(j)) THEN
255 value = local_data(i, j)
256 row_cell = (row_indices(i) - 1)/block_size + 1
257 col_cell = (col_indices(j) - 1)/block_size + 1
258 alpha = modulo(row_indices(i) - 1, block_size) + 1
259 beta = modulo(col_indices(j) - 1, block_size) + 1
260 delta = modulo(index_to_cell(:, row_cell) - index_to_cell(:, col_cell), &
261 preconditioner_env%lattice_dims)
262 delta_index = delta(1) + preconditioner_env%lattice_dims(1)*(delta(2) + &
263 preconditioner_env%lattice_dims(2)*delta(3)) + 1
264 blocks(alpha, beta, delta_index) = blocks(alpha, beta, delta_index) + &
265 value/real(ncell, kind=dp)
266 denominator = denominator + value*value
267 IF (ALLOCATED(dense_operator)) THEN
268 dense_operator(row_indices(i), col_indices(j)) = value
269 END IF
270 IF (row_indices(i) < col_indices(j)) THEN
271 delta = modulo(index_to_cell(:, col_cell) - index_to_cell(:, row_cell), &
272 preconditioner_env%lattice_dims)
273 delta_index = delta(1) + preconditioner_env%lattice_dims(1)*(delta(2) + &
274 preconditioner_env%lattice_dims(2)*delta(3)) + 1
275 blocks(beta, alpha, delta_index) = blocks(beta, alpha, delta_index) + &
276 value/real(ncell, kind=dp)
277 denominator = denominator + value*value
278 END IF
279 END IF
280 END DO
281 END DO
282 CALL preconditioner_env%para_env%sum(blocks)
283 numerator = 0.0_dp
284 DO j = 1, ncol_local
285 DO i = 1, nrow_local
286 IF (row_indices(i) <= col_indices(j)) THEN
287 value = local_data(i, j)
288 row_cell = (row_indices(i) - 1)/block_size + 1
289 col_cell = (col_indices(j) - 1)/block_size + 1
290 alpha = modulo(row_indices(i) - 1, block_size) + 1
291 beta = modulo(col_indices(j) - 1, block_size) + 1
292 delta = modulo(index_to_cell(:, row_cell) - index_to_cell(:, col_cell), &
293 preconditioner_env%lattice_dims)
294 delta_index = delta(1) + preconditioner_env%lattice_dims(1)*(delta(2) + &
295 preconditioner_env%lattice_dims(2)*delta(3)) + 1
296 difference = value - blocks(alpha, beta, delta_index)
297 numerator = numerator + difference*difference
298 IF (row_indices(i) < col_indices(j)) THEN
299 delta = modulo(index_to_cell(:, col_cell) - index_to_cell(:, row_cell), &
300 preconditioner_env%lattice_dims)
301 delta_index = delta(1) + preconditioner_env%lattice_dims(1)*(delta(2) + &
302 preconditioner_env%lattice_dims(2)*delta(3)) + 1
303 difference = value - blocks(beta, alpha, delta_index)
304 numerator = numerator + difference*difference
305 END IF
306 END IF
307 END DO
308 END DO
309 CALL preconditioner_env%para_env%sum(numerator)
310 CALL preconditioner_env%para_env%sum(denominator)
311 preconditioner_env%lattice_projection_error = &
312 sqrt(numerator/max(denominator, tiny(denominator)))
314 blocks, preconditioner_env%lattice_dims, preconditioner_env%lattice_inverse_k, used_fft, info)
315 preconditioner_env%lattice_corrected_projection_error = &
316 preconditioner_env%lattice_projection_error
317 IF (ALLOCATED(dense_operator)) THEN
318 CALL preconditioner_env%para_env%sum(dense_operator)
319 DO j = 1, n
320 DO i = j + 1, n
321 dense_operator(i, j) = dense_operator(j, i)
322 END DO
323 END DO
325 dense_operator, preconditioner_env%lattice_dims, &
326 preconditioner_env%lattice_local_cells, &
327 preconditioner_env%lattice_local_indices, &
328 preconditioner_env%lattice_local_operator, &
329 preconditioner_env%lattice_local_coarse_inverse, &
330 preconditioner_env%lattice_corrected_projection_error, local_info)
331 IF (local_info /= 0) THEN
332 cpabort("LATTICE_FFT local correction found a non-positive coarse operator")
333 END IF
334 DEALLOCATE (dense_operator)
335 END IF
336 CALL cp_fm_release(operator_fm)
337 DEALLOCATE (blocks, index_to_cell, xkp)
338 SELECT CASE (preconditioner_env%lattice_fft_mode)
339 CASE (ot_lattice_fft_on)
340 IF (info /= 0) THEN
341 cpabort("LATTICE_FFT found a non-positive Fourier block for "//trim(operator_name))
342 END IF
343 IF (.NOT. used_fft) THEN
344 cpabort("LATTICE_FFT could not initialize its FFT for "//trim(operator_name))
345 END IF
346 preconditioner_env%lattice_fft_active = .true.
348 preconditioner_env%lattice_fft_active = lattice_fft_auto_select( &
349 preconditioner_env%lattice_dims, &
350 preconditioner_env%lattice_block_size, rhs_count, 1, &
351 preconditioner_env%para_env%num_pe, &
352 preconditioner_env%lattice_projection_error, used_fft, info)
353 CASE DEFAULT
354 cpabort("Unknown LATTICE_FFT mode")
355 END SELECT
356
357 IF (.NOT. preconditioner_env%lattice_fft_active .AND. &
358 ALLOCATED(preconditioner_env%lattice_inverse_k)) THEN
359 DEALLOCATE (preconditioner_env%lattice_inverse_k)
360 END IF
361 IF (preconditioner_env%para_env%is_source()) THEN
362 logger => cp_get_default_logger()
363 output_unit = cp_logger_get_default_io_unit(logger)
364 IF (output_unit > 0) THEN
365 WRITE (output_unit, '(T2,A,A,A,3(I0,1X),A,I0)') &
366 "OT| Lattice FFT ", trim(operator_name), " cells: ", &
367 preconditioner_env%lattice_dims, "AO block: ", preconditioner_env%lattice_block_size
368 WRITE (output_unit, '(T2,A,ES12.4,A,L1)') &
369 "OT| Block-circulant projection error: ", &
370 preconditioner_env%lattice_projection_error, &
371 " selected: ", preconditioner_env%lattice_fft_active
372 IF (preconditioner_env%lattice_local_cells > 0) THEN
373 WRITE (output_unit, '(T2,A,I0,A,ES12.4)') &
374 "OT| Balanced local correction cells: ", &
375 preconditioner_env%lattice_local_cells, " residual error: ", &
376 preconditioner_env%lattice_corrected_projection_error
377 END IF
378 WRITE (output_unit, '(T2,A,F8.3,A,F8.3)') &
379 "OT| Estimated lattice/dense work ratio: ", &
380 preconditioner_env%lattice_cost_ratio, " storage ratio: ", &
381 preconditioner_env%lattice_storage_ratio
382 END IF
383 END IF
384
385 END SUBROUTINE prepare_lattice_operator
386
387! **************************************************************************************************
388!> \brief Build state-dependent block-circulant inverses for rotation-covariant FULL_ALL.
389!> \param preconditioner_env ...
390!> \param matrix_s overlap matrix
391!> \param operator_name ...
392! **************************************************************************************************
393 SUBROUTINE prepare_lattice_state_operator(preconditioner_env, matrix_s, operator_name)
394
395 TYPE(preconditioner_type) :: preconditioner_env
396 TYPE(dbcsr_type), POINTER :: matrix_s
397 CHARACTER(LEN=*), INTENT(IN) :: operator_name
398
399 COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: inverse_k
400 INTEGER :: i, info, istate, j, local_info, &
401 local_rank, n, ncell, ncol_local, &
402 nrow_local, nstate, output_unit
403 INTEGER, ALLOCATABLE, DIMENSION(:) :: local_indices
404 INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
405 LOGICAL :: all_selected, selected, used_fft
406 REAL(KIND=dp) :: corrected_projection_error, &
407 projection_error
408 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: dense_operator, dense_sv, dense_vectors, &
409 local_coarse_inverse, local_operator, &
410 scaled_sv
411 REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :), &
412 POINTER :: local_data
413 TYPE(cp_fm_type) :: sv_fm
414 TYPE(cp_logger_type), POINTER :: logger
415
416 cpassert(ASSOCIATED(matrix_s))
417 cpassert(ASSOCIATED(preconditioner_env%fm))
418 cpassert(ASSOCIATED(preconditioner_env%full_evals))
419 cpassert(ASSOCIATED(preconditioner_env%occ_evals))
420 cpassert(preconditioner_env%lattice_fft_mode /= ot_lattice_fft_off)
421
422 preconditioner_env%lattice_fft_active = .false.
423 preconditioner_env%lattice_projection_error = 0.0_dp
424 preconditioner_env%lattice_corrected_projection_error = 0.0_dp
425 preconditioner_env%lattice_cost_ratio = 0.0_dp
426 preconditioner_env%lattice_storage_ratio = 0.0_dp
427 IF (ALLOCATED(preconditioner_env%lattice_state_inverse_k)) THEN
428 DEALLOCATE (preconditioner_env%lattice_state_inverse_k)
429 END IF
430 IF (ALLOCATED(preconditioner_env%lattice_occ_vectors)) THEN
431 DEALLOCATE (preconditioner_env%lattice_occ_vectors)
432 END IF
433 IF (ALLOCATED(preconditioner_env%lattice_occ_dual)) THEN
434 DEALLOCATE (preconditioner_env%lattice_occ_dual)
435 END IF
436 IF (ALLOCATED(preconditioner_env%lattice_local_indices)) THEN
437 DEALLOCATE (preconditioner_env%lattice_local_indices)
438 END IF
439 IF (ALLOCATED(preconditioner_env%lattice_local_operator)) THEN
440 DEALLOCATE (preconditioner_env%lattice_local_operator)
441 END IF
442 IF (ALLOCATED(preconditioner_env%lattice_local_coarse_inverse)) THEN
443 DEALLOCATE (preconditioner_env%lattice_local_coarse_inverse)
444 END IF
445 IF (ALLOCATED(preconditioner_env%lattice_state_local_indices)) THEN
446 DEALLOCATE (preconditioner_env%lattice_state_local_indices)
447 END IF
448 IF (ALLOCATED(preconditioner_env%lattice_state_local_operator)) THEN
449 DEALLOCATE (preconditioner_env%lattice_state_local_operator)
450 END IF
451 IF (ALLOCATED(preconditioner_env%lattice_state_local_coarse_inverse)) THEN
452 DEALLOCATE (preconditioner_env%lattice_state_local_coarse_inverse)
453 END IF
454 IF (any(preconditioner_env%lattice_dims <= 0)) cpabort("Lattice dimensions must be positive")
455 IF (preconditioner_env%lattice_fft_mode == ot_lattice_fft_auto .AND. &
456 product(preconditioner_env%lattice_dims) < lattice_fft_auto_min_cells) THEN
457 IF (preconditioner_env%para_env%is_source()) THEN
458 logger => cp_get_default_logger()
459 output_unit = cp_logger_get_default_io_unit(logger)
460 IF (output_unit > 0) THEN
461 WRITE (output_unit, '(T2,A,A,A,I0,A,I0)') &
462 "OT| Lattice FFT AUTO kept the conventional ", trim(operator_name), &
463 " solver: ", product(preconditioner_env%lattice_dims), &
464 " cells are fewer than ", lattice_fft_auto_min_cells
465 END IF
466 END IF
467 RETURN
468 END IF
469
470 CALL cp_fm_get_info(preconditioner_env%fm, nrow_global=n, ncol_global=j)
471 cpassert(j == n .AND. SIZE(preconditioner_env%full_evals) == n)
472 nstate = SIZE(preconditioner_env%occ_evals)
473 cpassert(nstate > 0)
474 ncell = product(preconditioner_env%lattice_dims)
475 IF (modulo(n, ncell) /= 0) cpabort("AO dimension does not match replicated cells")
476 preconditioner_env%lattice_block_size = n/ncell
477 preconditioner_env%lattice_cost_ratio = lattice_fft_cost_ratio( &
478 preconditioner_env%lattice_dims, &
479 preconditioner_env%lattice_block_size, nstate, &
480 preconditioner_env%para_env%num_pe)
481 preconditioner_env%lattice_storage_ratio = lattice_fft_storage_ratio( &
482 preconditioner_env%lattice_dims, nstate)
483 preconditioner_env%lattice_cost_ratio = preconditioner_env%lattice_cost_ratio + &
485 preconditioner_env%lattice_dims, &
486 preconditioner_env%lattice_local_cells)
487 preconditioner_env%lattice_storage_ratio = preconditioner_env%lattice_storage_ratio + &
489 preconditioner_env%lattice_dims, nstate, &
490 preconditioner_env%lattice_local_cells)
491 IF (preconditioner_env%lattice_fft_mode == ot_lattice_fft_auto .AND. &
492 (preconditioner_env%lattice_cost_ratio > lattice_fft_auto_max_cost_ratio .OR. &
493 preconditioner_env%lattice_storage_ratio > lattice_fft_auto_max_storage_ratio)) THEN
494 IF (preconditioner_env%para_env%is_source()) THEN
495 logger => cp_get_default_logger()
496 output_unit = cp_logger_get_default_io_unit(logger)
497 IF (output_unit > 0) THEN
498 WRITE (output_unit, '(T2,A,A,A,F8.3,A,F8.3,A)') &
499 "OT| Lattice FFT AUTO kept the conventional ", trim(operator_name), &
500 " solver: estimated work ratio ", preconditioner_env%lattice_cost_ratio, &
501 " storage ratio ", preconditioner_env%lattice_storage_ratio, " selected: F"
502 END IF
503 END IF
504 RETURN
505 END IF
506 IF (preconditioner_env%lattice_local_cells > 0) THEN
507 local_rank = preconditioner_env%lattice_local_cells* &
508 preconditioner_env%lattice_block_size
509 ALLOCATE (preconditioner_env%lattice_state_local_indices(local_rank, nstate), &
510 preconditioner_env%lattice_state_local_operator(n, local_rank, nstate), &
511 preconditioner_env%lattice_state_local_coarse_inverse(local_rank, local_rank, nstate))
512 END IF
513 CALL cp_fm_create(sv_fm, preconditioner_env%fm%matrix_struct, name="lattice FULL_ALL S V")
514 CALL cp_dbcsr_sm_fm_multiply(matrix_s, preconditioner_env%fm, sv_fm, n)
515 CALL cp_fm_get_info(sv_fm, nrow_local=nrow_local, ncol_local=ncol_local, &
516 row_indices=row_indices, col_indices=col_indices, local_data=local_data)
517 ALLOCATE (dense_sv(n, n), dense_vectors(n, nstate), source=0.0_dp)
518 DO j = 1, ncol_local
519 DO i = 1, nrow_local
520 dense_sv(row_indices(i), col_indices(j)) = local_data(i, j)
521 END DO
522 END DO
523 CALL preconditioner_env%para_env%sum(dense_sv)
524 CALL cp_fm_release(sv_fm)
525
526 CALL cp_fm_get_info(preconditioner_env%fm, nrow_local=nrow_local, ncol_local=ncol_local, &
527 row_indices=row_indices, col_indices=col_indices, local_data=local_data)
528 DO j = 1, ncol_local
529 IF (col_indices(j) > nstate) cycle
530 DO i = 1, nrow_local
531 dense_vectors(row_indices(i), col_indices(j)) = local_data(i, j)
532 END DO
533 END DO
534 CALL preconditioner_env%para_env%sum(dense_vectors)
535 ALLOCATE (preconditioner_env%lattice_occ_vectors(n, nstate), &
536 preconditioner_env%lattice_occ_dual(n, nstate))
537 preconditioner_env%lattice_occ_vectors(:, :) = dense_vectors
538 preconditioner_env%lattice_occ_dual(:, :) = dense_sv(:, 1:nstate)
539
540 ALLOCATE (dense_operator(n, n), scaled_sv(n, n))
541 all_selected = .true.
542 DO istate = 1, nstate
543 scaled_sv(:, :) = dense_sv
544 DO j = 1, n
545 scaled_sv(:, j) = scaled_sv(:, j)*max( &
546 preconditioner_env%energy_gap, &
547 preconditioner_env%full_evals(j) - preconditioner_env%occ_evals(istate))
548 END DO
549 dense_operator(:, :) = matmul(scaled_sv, transpose(dense_sv))
550 CALL build_lattice_inverse(dense_operator, preconditioner_env%lattice_dims, inverse_k, &
551 projection_error, used_fft, info)
552 preconditioner_env%lattice_projection_error = &
553 max(preconditioner_env%lattice_projection_error, projection_error)
554 corrected_projection_error = projection_error
555 IF (preconditioner_env%lattice_local_cells > 0) THEN
557 dense_operator, preconditioner_env%lattice_dims, &
558 preconditioner_env%lattice_local_cells, local_indices, &
559 local_operator, local_coarse_inverse, corrected_projection_error, local_info)
560 IF (local_info /= 0) THEN
561 cpabort("LATTICE_FFT local correction found a non-positive coarse operator")
562 END IF
563 preconditioner_env%lattice_state_local_indices(:, istate) = local_indices
564 preconditioner_env%lattice_state_local_operator(:, :, istate) = local_operator
565 preconditioner_env%lattice_state_local_coarse_inverse(:, :, istate) = local_coarse_inverse
566 DEALLOCATE (local_coarse_inverse, local_indices, local_operator)
567 END IF
568 preconditioner_env%lattice_corrected_projection_error = &
569 max(preconditioner_env%lattice_corrected_projection_error, corrected_projection_error)
570 selected = lattice_fft_auto_select(preconditioner_env%lattice_dims, &
571 preconditioner_env%lattice_block_size, nstate, nstate, &
572 preconditioner_env%para_env%num_pe, projection_error, used_fft, info)
573 IF (preconditioner_env%lattice_fft_mode == ot_lattice_fft_on) THEN
574 IF (info /= 0) THEN
575 cpabort("LATTICE_FFT found a non-positive Fourier block for "//trim(operator_name))
576 END IF
577 IF (.NOT. used_fft) THEN
578 cpabort("LATTICE_FFT could not initialize its FFT for "//trim(operator_name))
579 END IF
580 ELSE
581 all_selected = all_selected .AND. selected
582 END IF
583 IF (.NOT. ALLOCATED(preconditioner_env%lattice_state_inverse_k)) THEN
584 ALLOCATE (preconditioner_env%lattice_state_inverse_k( &
585 SIZE(inverse_k, 1), SIZE(inverse_k, 2), SIZE(inverse_k, 3), nstate))
586 END IF
587 IF (info == 0 .AND. used_fft) THEN
588 preconditioner_env%lattice_state_inverse_k(:, :, :, istate) = inverse_k
589 END IF
590 DEALLOCATE (inverse_k)
591 END DO
592
593 preconditioner_env%lattice_fft_active = &
594 preconditioner_env%lattice_fft_mode == ot_lattice_fft_on .OR. all_selected
595 IF (.NOT. preconditioner_env%lattice_fft_active .AND. &
596 ALLOCATED(preconditioner_env%lattice_state_inverse_k)) THEN
597 DEALLOCATE (preconditioner_env%lattice_state_inverse_k)
598 END IF
599 IF (.NOT. preconditioner_env%lattice_fft_active) THEN
600 DEALLOCATE (preconditioner_env%lattice_occ_vectors, &
601 preconditioner_env%lattice_occ_dual)
602 IF (ALLOCATED(preconditioner_env%lattice_state_local_indices)) THEN
603 DEALLOCATE (preconditioner_env%lattice_state_local_indices, &
604 preconditioner_env%lattice_state_local_operator, &
605 preconditioner_env%lattice_state_local_coarse_inverse)
606 END IF
607 END IF
608 IF (preconditioner_env%para_env%is_source()) THEN
609 logger => cp_get_default_logger()
610 output_unit = cp_logger_get_default_io_unit(logger)
611 IF (output_unit > 0) THEN
612 WRITE (output_unit, '(T2,A,A,A,3(I0,1X),A,I0,A,I0)') &
613 "OT| Lattice FFT ", trim(operator_name), " cells: ", &
614 preconditioner_env%lattice_dims, "AO block: ", &
615 preconditioner_env%lattice_block_size, " states: ", nstate
616 WRITE (output_unit, '(T2,A,ES12.4,A,L1)') &
617 "OT| Maximum block-circulant projection error: ", &
618 preconditioner_env%lattice_projection_error, &
619 " selected: ", preconditioner_env%lattice_fft_active
620 IF (preconditioner_env%lattice_local_cells > 0) THEN
621 WRITE (output_unit, '(T2,A,I0,A,ES12.4)') &
622 "OT| Balanced local correction cells/state: ", &
623 preconditioner_env%lattice_local_cells, " maximum residual error: ", &
624 preconditioner_env%lattice_corrected_projection_error
625 END IF
626 WRITE (output_unit, '(T2,A,F8.3,A,F8.3)') &
627 "OT| Estimated lattice/dense work ratio: ", &
628 preconditioner_env%lattice_cost_ratio, " storage ratio: ", &
629 preconditioner_env%lattice_storage_ratio
630 END IF
631 END IF
632
633 DEALLOCATE (dense_operator, dense_sv, dense_vectors, scaled_sv)
634
635 END SUBROUTINE prepare_lattice_state_operator
636
637! **************************************************************************************************
638!> \brief Build the overlap-inverse base of the bounded low-rank spectral preconditioner.
639!> \param preconditioner_env ...
640!> \param matrix_s ...
641! **************************************************************************************************
642 SUBROUTINE prepare_low_rank_base(preconditioner_env, matrix_s)
643
644 TYPE(preconditioner_type) :: preconditioner_env
645 TYPE(dbcsr_type), OPTIONAL, POINTER :: matrix_s
646
647 INTEGER :: info
648 TYPE(cp_fm_type) :: work
649
650 cpassert(ASSOCIATED(preconditioner_env%fm))
651 cpassert(.NOT. ASSOCIATED(preconditioner_env%base_fm))
652
653 IF (preconditioner_env%lattice_fft_mode /= ot_lattice_fft_off) THEN
654 IF (.NOT. PRESENT(matrix_s)) THEN
655 cpabort("LATTICE_FFT requires a non-unit overlap metric for FERMI_LOW_RANK")
656 END IF
657 cpassert(ASSOCIATED(matrix_s))
658 CALL prepare_lattice_operator(preconditioner_env, matrix_s, "FERMI_LOW_RANK overlap", &
659 SIZE(preconditioner_env%occ_evals))
660 IF (preconditioner_env%lattice_fft_active) THEN
661 preconditioner_env%low_rank_base = ot_low_rank_base_lattice_fft
662 RETURN
663 END IF
664 preconditioner_env%low_rank_base = ot_low_rank_base_overlap
665 END IF
666
667 IF (preconditioner_env%low_rank_base /= ot_low_rank_base_overlap) THEN
668 cpabort("Unknown FERMI_LOW_RANK base operator")
669 END IF
670 ALLOCATE (preconditioner_env%base_fm)
671 CALL cp_fm_create(preconditioner_env%base_fm, preconditioner_env%fm%matrix_struct, &
672 name="FERMI_LOW_RANK overlap inverse")
673
674 IF (PRESENT(matrix_s)) THEN
675 cpassert(ASSOCIATED(matrix_s))
676 CALL copy_dbcsr_to_fm(matrix_s, preconditioner_env%base_fm)
677 CALL cp_fm_cholesky_decompose(preconditioner_env%base_fm, info_out=info)
678 IF (info /= 0) cpabort("FERMI_LOW_RANK requires a positive-definite overlap matrix")
679 CALL cp_fm_cholesky_invert(preconditioner_env%base_fm)
680 CALL cp_fm_create(work, preconditioner_env%base_fm%matrix_struct, name="low-rank base work")
681 CALL cp_fm_uplo_to_full(preconditioner_env%base_fm, work)
682 CALL cp_fm_release(work)
683 ELSE
684 CALL cp_fm_set_all(preconditioner_env%base_fm, alpha=0.0_dp, beta=1.0_dp)
685 END IF
686
687 END SUBROUTINE prepare_low_rank_base
688
689! **************************************************************************************************
690
691! creates a preconditioner for the system (H-energy_homo S)
692! this preconditioner is (must be) symmetric positive definite.
693! currently uses a atom-block-diagonal form
694! each block will be ....
695! might overwrite matrix_h, matrix_t
696
697! **************************************************************************************************
698!> \brief ...
699!> \param preconditioner_env ...
700!> \param precon_type ...
701!> \param solver_type ...
702!> \param matrix_h ...
703!> \param matrix_s ...
704!> \param matrix_t ...
705!> \param mo_set ...
706!> \param energy_gap ...
707!> \param convert_precond_to_dbcsr ...
708!> \param chol_type ...
709!> \par History
710!> 09.2014 removed some unused or unfinished methods
711!> removed sparse preconditioners and the
712!> sparse approximate inverse at rev 14341 [Florian Schiffmann]
713! **************************************************************************************************
714 SUBROUTINE make_preconditioner(preconditioner_env, precon_type, solver_type, matrix_h, matrix_s, &
715 matrix_t, mo_set, energy_gap, convert_precond_to_dbcsr, chol_type)
716
717 TYPE(preconditioner_type) :: preconditioner_env
718 INTEGER, INTENT(IN) :: precon_type, solver_type
719 TYPE(dbcsr_type), POINTER :: matrix_h
720 TYPE(dbcsr_type), OPTIONAL, POINTER :: matrix_s, matrix_t
721 TYPE(mo_set_type), INTENT(IN) :: mo_set
722 REAL(kind=dp) :: energy_gap
723 LOGICAL, INTENT(IN), OPTIONAL :: convert_precond_to_dbcsr
724 INTEGER, INTENT(IN), OPTIONAL :: chol_type
725
726 CHARACTER(len=*), PARAMETER :: routinen = 'make_preconditioner'
727
728 INTEGER :: handle, k, max_rank, my_solver_type, &
729 nao, nhomo, output_unit
730 LOGICAL :: my_convert_precond_to_dbcsr, &
731 needs_full_spectrum, needs_homo, &
732 use_mo_coeff_b
733 REAL(kind=dp) :: energy_homo
734 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: eigenvalues_ot
735 TYPE(cp_fm_struct_type), POINTER :: fm_struct
736 TYPE(cp_fm_type) :: mo_occ
737 TYPE(cp_fm_type), POINTER :: mo_coeff
738 TYPE(cp_logger_type), POINTER :: logger
739 TYPE(dbcsr_type), POINTER :: mo_coeff_b
740
741 CALL timeset(routinen, handle)
742
743 CALL get_mo_set(mo_set=mo_set, mo_coeff=mo_coeff, mo_coeff_b=mo_coeff_b, homo=nhomo)
744 use_mo_coeff_b = mo_set%use_mo_coeff_b
745 CALL cp_fm_get_info(mo_coeff, ncol_global=k, nrow_global=nao)
746
747 ! Starting some matrix mess, check where to store the result in preconditioner_env, fm or dbcsr_matrix
748 my_convert_precond_to_dbcsr = .false.
749 IF (PRESENT(convert_precond_to_dbcsr)) my_convert_precond_to_dbcsr = convert_precond_to_dbcsr
750 IF (precon_type == ot_precond_fermi_low_rank .AND. &
751 preconditioner_env%low_rank_base == ot_low_rank_base_lattice_fft) THEN
752 ! The FFT base is an operator rather than a stored DBCSR matrix. Keep the
753 ! complementary vectors in their distributed full-matrix representation.
754 my_convert_precond_to_dbcsr = .false.
755 END IF
756
757 ! Thanks to the mess with the matrices we need to make sure in this case that the
758 ! Previous inverse is properly stored as a sparse matrix, fm gets deallocated here
759 ! if it wasn't anyway
760 IF (preconditioner_env%solver == ot_precond_solver_update) THEN
761 CALL transfer_fm_to_dbcsr(preconditioner_env%fm, preconditioner_env%dbcsr_matrix, matrix_h)
762 END IF
763
764 needs_full_spectrum = .false.
765 needs_homo = .false.
766
767 SELECT CASE (precon_type)
769 needs_full_spectrum = .true.
770 ! both of them need the coefficients as fm's, more matrix mess
771 IF (use_mo_coeff_b) THEN
772 CALL copy_dbcsr_to_fm(mo_coeff_b, mo_coeff)
773 END IF
775 needs_homo = .true.
776 ! XXXX to be removed if homo estimate only is implemented
777 needs_full_spectrum = .true.
779 ! these should be happy without an estimate for the homo energy
780 ! preconditioning can not depend on an absolute eigenvalue, only on eigenvalue differences
781 CASE DEFAULT
782 cpabort("The preconditioner is unknown ...")
783 END SELECT
784
785 ALLOCATE (eigenvalues_ot(k))
786 energy_homo = 0.0_dp
787 IF (needs_full_spectrum) THEN
788 ! XXXXXXXXXXXXXXXX do not touch the initial MOs, could be harmful for either
789 ! the case of non-equivalent MOs but also for the derivate
790 ! we could already have all eigenvalues e.g. full_all and we could skip this
791 ! to be optimised later.
792 ! one flaw is that not all SCF methods (i.e. that go over mo_derivs directly)
793 ! have a 'valid' matrix_h... (we even don't know what evals are in that case)
794 IF (use_mo_coeff_b) THEN
795 CALL calculate_subspace_eigenvalues(mo_coeff_b, matrix_h, &
796 eigenvalues_ot, do_rotation=.false., &
797 para_env=mo_coeff%matrix_struct%para_env, &
798 blacs_env=mo_coeff%matrix_struct%context)
799 ELSE
800 CALL calculate_subspace_eigenvalues(mo_coeff, matrix_h, &
801 eigenvalues_ot, do_rotation=.false.)
802 END IF
803 IF (k > 0) THEN
804 cpassert(nhomo > 0 .AND. nhomo <= k)
805 energy_homo = eigenvalues_ot(nhomo)
806 END IF
807 ELSE
808 IF (needs_homo) THEN
809 cpabort("Not yet implemented")
810 END IF
811 END IF
812
813 ! After all bits and pieces of checking and initialization, here comes the
814 ! part where the preconditioner matrix gets created and solved.
815 ! This will give the matrices for later use
816 my_solver_type = solver_type
817 preconditioner_env%in_use = precon_type
818 preconditioner_env%cholesky_use = cholesky_reduce
819 IF (PRESENT(chol_type)) preconditioner_env%cholesky_use = chol_type
820 preconditioner_env%in_use = precon_type
821 preconditioner_env%lattice_fft_active = .false.
822 preconditioner_env%lattice_corrected_projection_error = 0.0_dp
823 IF (ALLOCATED(preconditioner_env%lattice_inverse_k)) THEN
824 DEALLOCATE (preconditioner_env%lattice_inverse_k)
825 END IF
826 IF (ALLOCATED(preconditioner_env%lattice_state_inverse_k)) THEN
827 DEALLOCATE (preconditioner_env%lattice_state_inverse_k)
828 END IF
829 IF (ALLOCATED(preconditioner_env%lattice_occ_vectors)) THEN
830 DEALLOCATE (preconditioner_env%lattice_occ_vectors)
831 END IF
832 IF (ALLOCATED(preconditioner_env%lattice_occ_dual)) THEN
833 DEALLOCATE (preconditioner_env%lattice_occ_dual)
834 END IF
835 IF (ALLOCATED(preconditioner_env%lattice_local_indices)) THEN
836 DEALLOCATE (preconditioner_env%lattice_local_indices)
837 END IF
838 IF (ALLOCATED(preconditioner_env%lattice_local_operator)) THEN
839 DEALLOCATE (preconditioner_env%lattice_local_operator)
840 END IF
841 IF (ALLOCATED(preconditioner_env%lattice_local_coarse_inverse)) THEN
842 DEALLOCATE (preconditioner_env%lattice_local_coarse_inverse)
843 END IF
844 IF (ALLOCATED(preconditioner_env%lattice_state_local_indices)) THEN
845 DEALLOCATE (preconditioner_env%lattice_state_local_indices)
846 END IF
847 IF (ALLOCATED(preconditioner_env%lattice_state_local_operator)) THEN
848 DEALLOCATE (preconditioner_env%lattice_state_local_operator)
849 END IF
850 IF (ALLOCATED(preconditioner_env%lattice_state_local_coarse_inverse)) THEN
851 DEALLOCATE (preconditioner_env%lattice_state_local_coarse_inverse)
852 END IF
853 IF (nhomo == k) THEN
854 CALL make_preconditioner_matrix(preconditioner_env, matrix_h, matrix_s, matrix_t, mo_coeff, &
855 energy_homo, eigenvalues_ot, energy_gap, my_solver_type)
856 ELSE
857 CALL cp_fm_struct_create(fm_struct, nrow_global=nao, ncol_global=nhomo, &
858 context=preconditioner_env%ctxt, &
859 para_env=preconditioner_env%para_env)
860 CALL cp_fm_create(mo_occ, fm_struct)
861 CALL cp_fm_to_fm(mo_coeff, mo_occ, nhomo)
862 CALL cp_fm_struct_release(fm_struct)
863 !
864 CALL make_preconditioner_matrix(preconditioner_env, matrix_h, matrix_s, matrix_t, mo_occ, &
865 energy_homo, eigenvalues_ot(1:nhomo), energy_gap, my_solver_type)
866 !
867 CALL cp_fm_release(mo_occ)
868 END IF
869
870 IF (preconditioner_env%lattice_fft_mode /= ot_lattice_fft_off) THEN
871 SELECT CASE (precon_type)
873 cpassert(ASSOCIATED(preconditioner_env%sparse_matrix))
874 CALL prepare_lattice_operator(preconditioner_env, preconditioner_env%sparse_matrix, &
875 "FULL_KINETIC", nhomo)
877 cpassert(ASSOCIATED(preconditioner_env%sparse_matrix))
878 CALL prepare_lattice_operator(preconditioner_env, preconditioner_env%sparse_matrix, &
879 "FULL_SINGLE_INVERSE", nhomo)
881 cpassert(ASSOCIATED(preconditioner_env%sparse_matrix))
882 CALL prepare_lattice_operator(preconditioner_env, preconditioner_env%sparse_matrix, &
883 "FULL_S_INVERSE", nhomo)
885 cpassert(PRESENT(matrix_s))
886 cpassert(ASSOCIATED(matrix_s))
887 CALL prepare_lattice_state_operator(preconditioner_env, matrix_s, &
888 "FULL_ALL_COVARIANT")
889 END SELECT
890 END IF
891
892 IF (.NOT. preconditioner_env%lattice_fft_active) THEN
893 CALL solve_preconditioner(my_solver_type, preconditioner_env, matrix_s, matrix_h)
894 ELSE
895 preconditioner_env%solver = ot_precond_solver_inv_chol
896 END IF
897
898 IF (precon_type == ot_precond_fermi_low_rank) THEN
899 preconditioner_env%spectral_reference = energy_homo
900 preconditioner_env%spectral_window = &
901 max(fermi_low_rank_min_window, fermi_low_rank_gap_factor*preconditioner_env%energy_gap)
902 preconditioner_env%spectral_base_scale = 1.0_dp/preconditioner_env%spectral_window
903 max_rank = preconditioner_env%spectral_max_rank
904 IF (max_rank == -1) max_rank = SIZE(preconditioner_env%full_evals) - nhomo
905 preconditioner_env%spectral_rank = low_rank_select_rank( &
906 preconditioner_env%full_evals, nhomo, max_rank, &
907 fermi_low_rank_degeneracy_tolerance)
908 CALL prepare_low_rank_base(preconditioner_env, matrix_s)
909 IF (preconditioner_env%lattice_fft_active) THEN
910 ! With an approximate base, retain only non-negative spectral updates.
911 ! This preserves positive definiteness independently of the projection error.
912 DO WHILE (preconditioner_env%spectral_rank > 0)
914 preconditioner_env%full_evals(nhomo + preconditioner_env%spectral_rank), &
915 preconditioner_env%spectral_reference, preconditioner_env%energy_gap) >= &
916 preconditioner_env%spectral_base_scale) EXIT
917 preconditioner_env%spectral_rank = preconditioner_env%spectral_rank - 1
918 END DO
919 END IF
920 IF (preconditioner_env%para_env%is_source()) THEN
921 logger => cp_get_default_logger()
922 output_unit = cp_logger_get_default_io_unit(logger)
923 IF (output_unit > 0) THEN
924 WRITE (output_unit, '(T2,A,F14.8,A,I0)') &
925 "OT| FERMI_LOW_RANK common reference [hartree]: ", &
926 preconditioner_env%spectral_reference, " retained rank: ", &
927 preconditioner_env%spectral_rank
928 END IF
929 END IF
930 END IF
931
932 ! Here comes more matrix mess, make sure to output the correct matrix format,
933 ! A bit pointless to convert the cholesky factorized version as it doesn't work in
934 ! dbcsr form and will crash later,...
935 IF (.NOT. preconditioner_env%lattice_fft_active) THEN
936 IF (my_convert_precond_to_dbcsr) THEN
937 CALL transfer_fm_to_dbcsr(preconditioner_env%fm, preconditioner_env%dbcsr_matrix, matrix_h)
938 IF (precon_type == ot_precond_fermi_low_rank .AND. &
939 preconditioner_env%low_rank_base == ot_low_rank_base_overlap) THEN
940 CALL transfer_fm_to_dbcsr(preconditioner_env%base_fm, &
941 preconditioner_env%base_dbcsr_matrix, matrix_h)
942 END IF
943 ELSE
944 CALL transfer_dbcsr_to_fm(preconditioner_env%dbcsr_matrix, preconditioner_env%fm, &
945 preconditioner_env%para_env, preconditioner_env%ctxt)
946 END IF
947 END IF
948
949 DEALLOCATE (eigenvalues_ot)
950
951 CALL timestop(handle)
952
953 END SUBROUTINE make_preconditioner
954
955! **************************************************************************************************
956!> \brief Construct FULL_ALL directly from one complex H(k), S(k), and C(k) channel.
957!> \param preconditioner_env preconditioner storage
958!> \param matrix_c_re real part of the active OT reference orbitals
959!> \param matrix_c_im imaginary part of the active OT reference orbitals
960!> \param matrix_h_re real part of H(k)
961!> \param matrix_h_im imaginary part of H(k)
962!> \param matrix_s_re real part of S(k)
963!> \param matrix_s_im imaginary part of S(k)
964!> \param energy_gap denominator floor
965!> \param solver_type inverse application selected in the OT input
966! **************************************************************************************************
967 SUBROUTINE make_preconditioner_complex_full_all(preconditioner_env, &
968 matrix_c_re, matrix_c_im, &
969 matrix_h_re, matrix_h_im, &
970 matrix_s_re, matrix_s_im, &
971 energy_gap, solver_type)
972
973 TYPE(preconditioner_type) :: preconditioner_env
974 TYPE(dbcsr_type), POINTER :: matrix_c_re, matrix_c_im, matrix_h_re, &
975 matrix_h_im, matrix_s_re, matrix_s_im
976 REAL(kind=dp), INTENT(IN) :: energy_gap
977 INTEGER, INTENT(IN) :: solver_type
978
979 CHARACTER(len=*), PARAMETER :: routinen = 'make_preconditioner_complex_full_all'
980
981 INTEGER :: handle
982 TYPE(cp_cfm_type) :: matrix_c, matrix_h, matrix_s
983
984 CALL timeset(routinen, handle)
985
986 cpassert(ASSOCIATED(matrix_c_re))
987 cpassert(ASSOCIATED(matrix_c_im))
988 cpassert(ASSOCIATED(matrix_h_re))
989 cpassert(ASSOCIATED(matrix_h_im))
990 cpassert(ASSOCIATED(matrix_s_re))
991 cpassert(ASSOCIATED(matrix_s_im))
992 IF (solver_type /= ot_precond_solver_default) THEN
993 cpabort('Complex FULL_ALL supports only PRECOND_SOLVER DEFAULT')
994 END IF
995 CALL dbcsr_pair_to_cfm(matrix_c_re, matrix_c_im, preconditioner_env, &
996 'complex FULL_ALL C', matrix_c)
997 CALL dbcsr_pair_to_cfm(matrix_h_re, matrix_h_im, preconditioner_env, &
998 'complex FULL_ALL H', matrix_h)
999 CALL dbcsr_pair_to_cfm(matrix_s_re, matrix_s_im, preconditioner_env, &
1000 'complex FULL_ALL S', matrix_s)
1001 CALL make_complex_full_all(preconditioner_env, matrix_c, matrix_h, matrix_s, &
1002 energy_gap=energy_gap)
1003
1004 CALL cp_cfm_release(matrix_c)
1005 CALL cp_cfm_release(matrix_s)
1006 CALL cp_cfm_release(matrix_h)
1007
1008 CALL timestop(handle)
1009
1011
1012! **************************************************************************************************
1013!> \brief Build local FULL_ALL preconditioners from cached operators in bounded batches.
1014!> Only the independent eigensolves run inside OpenMP.
1015!> \param preconditioners local channel preconditioners
1016!> \param matrix_c real/imaginary DBCSR reference orbitals
1017!> \param matrix_h cached real/imaginary Hamiltonians
1018!> \param matrix_s cached real/imaginary overlaps
1019!> \param energy_gap denominator floor
1020!> \param solver_type inverse application selected in the OT input
1021!> \param handled false when the existing distributed or explicitly selected solver is required
1022! **************************************************************************************************
1023 SUBROUTINE make_preconditioner_complex_full_all_batch(preconditioners, matrix_c, matrix_h, matrix_s, &
1024 energy_gap, solver_type, handled)
1025 TYPE(preconditioner_p_type), INTENT(IN) :: preconditioners(:)
1026 TYPE(dbcsr_p_type), INTENT(IN) :: matrix_c(:, :)
1027 TYPE(cp_fm_p_type), INTENT(IN) :: matrix_h(:, :), matrix_s(:, :)
1028 REAL(dp), INTENT(IN) :: energy_gap
1029 INTEGER, INTENT(IN) :: solver_type
1030 LOGICAL, INTENT(OUT) :: handled
1031
1032 CHARACTER(len=*), PARAMETER :: routinen = 'make_preconditioner_complex_full_all_batch'
1033
1034 INTEGER :: first, handle, i, j, k, n, nbatch, &
1035 nchannels, nworkers
1036 TYPE(full_all_work_type), ALLOCATABLE :: work(:)
1037 TYPE(preconditioner_type), POINTER :: p
1038
1039 handled = .false.
1040 nchannels = SIZE(preconditioners)
1041 IF (nchannels == 0 .OR. solver_type /= ot_precond_solver_default) RETURN
1043 DO i = 1, nchannels
1044 IF (preconditioners(i)%preconditioner%para_env%num_pe /= 1) RETURN
1045 END DO
1046 CALL timeset(routinen, handle)
1047 nworkers = 1
1048!$ nworkers = MIN(nchannels, omp_get_max_threads())
1049 ALLOCATE (work(nworkers))
1050 DO first = 1, nchannels, nworkers
1051 nbatch = min(nworkers, nchannels - first + 1)
1052 DO j = 1, nbatch
1053 i = first + j - 1
1054 p => preconditioners(i)%preconditioner
1055 CALL dbcsr_pair_to_cfm(matrix_c(1, i)%matrix, matrix_c(2, i)%matrix, p, 'FULL_ALL C', work(j)%c)
1056 CALL cp_cfm_create(work(j)%h, matrix_h(1, i)%matrix%matrix_struct, name='FULL_ALL H')
1057 CALL cp_cfm_create(work(j)%s, matrix_s(1, i)%matrix%matrix_struct, name='FULL_ALL S')
1058 CALL cp_fm_to_cfm(matrix_h(1, i)%matrix, matrix_h(2, i)%matrix, work(j)%h)
1059 CALL cp_fm_to_cfm(matrix_s(1, i)%matrix, matrix_s(2, i)%matrix, work(j)%s)
1060 CALL prepare_complex_full_all(p, work(j)%c, work(j)%h, work(j)%s, work(j)%a, work(j)%b, &
1061 energy_gap=energy_gap)
1062 n = SIZE(p%full_evals)
1063 ! Match the existing Cholesky reduction and back transformation exactly.
1064 CALL cp_cfm_cholesky_decompose(work(j)%b)
1065 CALL cp_cfm_triangular_invert(work(j)%b)
1066 CALL cp_cfm_triangular_multiply(work(j)%b, work(j)%a, side='R')
1067 CALL cp_cfm_triangular_multiply(work(j)%b, work(j)%a, transa_tr='C')
1068 CALL cp_cfm_local_workspace_prepare(work(j)%solver, n, work(j)%a%local_data)
1069 CALL cp_cfm_release(work(j)%s)
1070 CALL cp_cfm_release(work(j)%h)
1071 END DO
1072!$OMP PARALLEL DO DEFAULT(NONE) SCHEDULE(STATIC) NUM_THREADS(nbatch) IF(nbatch > 1) SHARED(work,nbatch)
1073 DO j = 1, nbatch
1074 CALL cp_cfm_heevd_local(work(j)%a%local_data, info=work(j)%info, ws=work(j)%solver)
1075 END DO
1076!$OMP END PARALLEL DO
1077 DO j = 1, nbatch
1078 IF (work(j)%info /= 0) cpabort('Local FULL_ALL generalized diagonalization failed')
1079 i = first + j - 1
1080 p => preconditioners(i)%preconditioner
1081 p%full_evals(:) = work(j)%solver%evals(:)
1082 CALL cp_cfm_triangular_multiply(work(j)%b, work(j)%a)
1083 CALL cp_cfm_to_cfm(work(j)%a, p%complex_fm)
1084 k = SIZE(p%occ_evals)
1085 p%full_evals(1:k) = p%occ_evals(:)
1086 CALL cp_cfm_to_cfm(work(j)%c, p%complex_fm, k)
1087 CALL cp_cfm_release(work(j)%b)
1088 CALL cp_cfm_release(work(j)%a)
1089 CALL cp_cfm_release(work(j)%c)
1090 END DO
1091 END DO
1092 handled = .true.
1093 CALL timestop(handle)
1095
1096! **************************************************************************************************
1097!> \brief Construct a complex FULL_SINGLE preconditioner from H(k) and S(k).
1098!> \param preconditioner_env preconditioner storage
1099!> \param matrix_h_re real part of H(k)
1100!> \param matrix_h_im imaginary part of H(k)
1101!> \param matrix_s_re real part of S(k)
1102!> \param matrix_s_im imaginary part of S(k)
1103!> \param mo_set orbital energy labels defining the occupied edge
1104!> \param energy_gap denominator floor
1105!> \param solver_type requested inverse solver
1106! **************************************************************************************************
1107 SUBROUTINE make_preconditioner_complex_full_single(preconditioner_env, &
1108 matrix_h_re, matrix_h_im, &
1109 matrix_s_re, matrix_s_im, &
1110 mo_set, energy_gap, solver_type)
1111
1112 TYPE(preconditioner_type) :: preconditioner_env
1113 TYPE(dbcsr_type), POINTER :: matrix_h_re, matrix_h_im, matrix_s_re, &
1114 matrix_s_im
1115 TYPE(mo_set_type), INTENT(IN) :: mo_set
1116 REAL(kind=dp), INTENT(IN) :: energy_gap
1117 INTEGER, INTENT(IN) :: solver_type
1118
1119 INTEGER :: homo
1120 REAL(kind=dp), DIMENSION(:), POINTER :: eigenvalues
1121 TYPE(cp_cfm_type) :: matrix_h, matrix_s
1122
1123 NULLIFY (eigenvalues)
1124 cpassert(ASSOCIATED(matrix_h_re))
1125 cpassert(ASSOCIATED(matrix_h_im))
1126 cpassert(ASSOCIATED(matrix_s_re))
1127 cpassert(ASSOCIATED(matrix_s_im))
1128 IF (solver_type /= ot_precond_solver_default) THEN
1129 cpabort('Complex FULL_SINGLE supports only PRECOND_SOLVER DEFAULT')
1130 END IF
1131 CALL get_mo_set(mo_set, homo=homo, eigenvalues=eigenvalues)
1132 cpassert(ASSOCIATED(eigenvalues))
1133 cpassert(homo > 0 .AND. homo <= SIZE(eigenvalues))
1134
1135 CALL dbcsr_pair_to_cfm(matrix_h_re, matrix_h_im, preconditioner_env, &
1136 'complex FULL_SINGLE H', matrix_h)
1137 CALL dbcsr_pair_to_cfm(matrix_s_re, matrix_s_im, preconditioner_env, &
1138 'complex FULL_SINGLE S', matrix_s)
1139 CALL make_complex_full_single(preconditioner_env, matrix_h, matrix_s, &
1140 eigenvalues(homo), energy_gap)
1141 CALL cp_cfm_release(matrix_s)
1142 CALL cp_cfm_release(matrix_h)
1143
1145
1146! **************************************************************************************************
1147!> \brief Construct FERMI_LOW_RANK directly for one complex k-point channel.
1148!> \param preconditioner_env preconditioner storage
1149!> \param matrix_c_re real part of the active OT reference orbitals
1150!> \param matrix_c_im imaginary part of the active OT reference orbitals
1151!> \param matrix_h_re real part of H(k)
1152!> \param matrix_h_im imaginary part of H(k)
1153!> \param matrix_s_re real part of S(k)
1154!> \param matrix_s_im imaginary part of S(k)
1155!> \param energy_gap denominator floor
1156!> \param max_rank maximum complementary rank, or -1 for the full complement
1157!> \param solver_type inverse application selected in the OT input
1158! **************************************************************************************************
1160 preconditioner_env, matrix_c_re, matrix_c_im, matrix_h_re, matrix_h_im, &
1161 matrix_s_re, matrix_s_im, energy_gap, max_rank, solver_type)
1162
1163 TYPE(preconditioner_type) :: preconditioner_env
1164 TYPE(dbcsr_type), POINTER :: matrix_c_re, matrix_c_im, matrix_h_re, &
1165 matrix_h_im, matrix_s_re, matrix_s_im
1166 REAL(kind=dp), INTENT(IN) :: energy_gap
1167 INTEGER, INTENT(IN) :: max_rank, solver_type
1168
1169 TYPE(cp_cfm_type) :: matrix_c, matrix_h, matrix_s
1170
1171 IF (solver_type /= ot_precond_solver_default) THEN
1172 cpabort('Complex FERMI_LOW_RANK supports only PRECOND_SOLVER DEFAULT')
1173 END IF
1174
1175 CALL dbcsr_pair_to_cfm(matrix_c_re, matrix_c_im, preconditioner_env, &
1176 'complex FERMI_LOW_RANK C', matrix_c)
1177 CALL dbcsr_pair_to_cfm(matrix_h_re, matrix_h_im, preconditioner_env, &
1178 'complex FERMI_LOW_RANK H', matrix_h)
1179 CALL dbcsr_pair_to_cfm(matrix_s_re, matrix_s_im, preconditioner_env, &
1180 'complex FERMI_LOW_RANK S', matrix_s)
1182 preconditioner_env, matrix_c, matrix_h, matrix_s, energy_gap, &
1183 max_rank, max(fermi_low_rank_min_window, fermi_low_rank_gap_factor*energy_gap), &
1184 fermi_low_rank_degeneracy_tolerance)
1185
1186 CALL cp_cfm_release(matrix_s)
1187 CALL cp_cfm_release(matrix_h)
1188 CALL cp_cfm_release(matrix_c)
1189
1191
1192! **************************************************************************************************
1193!> \brief Construct rotation-covariant FULL_ALL for one complex k-point channel.
1194!> \param preconditioner_env preconditioner storage
1195!> \param matrix_c_re real part of the active OT reference orbitals
1196!> \param matrix_c_im imaginary part of the active OT reference orbitals
1197!> \param matrix_h_re real part of H(k)
1198!> \param matrix_h_im imaginary part of H(k)
1199!> \param matrix_s_re real part of S(k)
1200!> \param matrix_s_im imaginary part of S(k)
1201!> \param energy_gap denominator floor
1202!> \param solver_type inverse application selected in the OT input
1203!> \param occupation_signature occupations that define the admissible private gauge
1204! **************************************************************************************************
1206 preconditioner_env, matrix_c_re, matrix_c_im, matrix_h_re, matrix_h_im, &
1207 matrix_s_re, matrix_s_im, energy_gap, solver_type, occupation_signature)
1208
1209 TYPE(preconditioner_type) :: preconditioner_env
1210 TYPE(dbcsr_type), POINTER :: matrix_c_re, matrix_c_im, matrix_h_re, &
1211 matrix_h_im, matrix_s_re, matrix_s_im
1212 REAL(kind=dp), INTENT(IN) :: energy_gap
1213 INTEGER, INTENT(IN) :: solver_type
1214 REAL(kind=dp), DIMENSION(:, :), INTENT(IN), &
1215 OPTIONAL :: occupation_signature
1216
1217 TYPE(cp_cfm_type) :: matrix_c, matrix_h, matrix_s
1218
1219 IF (solver_type /= ot_precond_solver_default) THEN
1220 cpabort('Complex FULL_ALL_COVARIANT supports only PRECOND_SOLVER DEFAULT')
1221 END IF
1222 CALL dbcsr_pair_to_cfm(matrix_c_re, matrix_c_im, preconditioner_env, &
1223 'complex FULL_ALL_COVARIANT C', matrix_c)
1224 CALL dbcsr_pair_to_cfm(matrix_h_re, matrix_h_im, preconditioner_env, &
1225 'complex FULL_ALL_COVARIANT H', matrix_h)
1226 CALL dbcsr_pair_to_cfm(matrix_s_re, matrix_s_im, preconditioner_env, &
1227 'complex FULL_ALL_COVARIANT S', matrix_s)
1228 CALL make_complex_full_all_covariant(preconditioner_env, matrix_c, matrix_h, matrix_s, &
1229 energy_gap, occupation_signature)
1230
1231 CALL cp_cfm_release(matrix_s)
1232 CALL cp_cfm_release(matrix_h)
1233 CALL cp_cfm_release(matrix_c)
1234
1236
1237! **************************************************************************************************
1238!> \brief Construct a complex FULL_SINGLE_INVERSE preconditioner without discarding Im(H,S,C).
1239!> \param preconditioner_env preconditioner storage
1240!> \param matrix_c_re real part of the occupied reference orbitals
1241!> \param matrix_c_im imaginary part of the occupied reference orbitals
1242!> \param matrix_h_re real part of H(k)
1243!> \param matrix_h_im imaginary part of H(k)
1244!> \param matrix_s_re real part of S(k)
1245!> \param matrix_s_im imaginary part of S(k)
1246!> \param energy_gap lower spectral bound
1247!> \param solver_type requested inverse solver
1248! **************************************************************************************************
1250 matrix_c_re, matrix_c_im, &
1251 matrix_h_re, matrix_h_im, &
1252 matrix_s_re, matrix_s_im, &
1253 energy_gap, solver_type)
1254
1255 TYPE(preconditioner_type) :: preconditioner_env
1256 TYPE(dbcsr_type), POINTER :: matrix_c_re, matrix_c_im, matrix_h_re, &
1257 matrix_h_im, matrix_s_re, matrix_s_im
1258 REAL(kind=dp), INTENT(IN) :: energy_gap
1259 INTEGER, INTENT(IN) :: solver_type
1260
1261 INTEGER :: k, n
1262 TYPE(cp_cfm_type) :: matrix_c, matrix_h, matrix_s
1263 TYPE(cp_fm_struct_type), POINTER :: fm_struct_mo, fm_struct_square
1264 TYPE(cp_fm_type) :: matrix_c_im_fm, matrix_c_re_fm, &
1265 matrix_h_im_fm, matrix_h_re_fm, &
1266 matrix_s_im_fm, matrix_s_re_fm
1267
1268 NULLIFY (fm_struct_mo, fm_struct_square)
1269 cpassert(ASSOCIATED(matrix_c_re))
1270 cpassert(ASSOCIATED(matrix_c_im))
1271 cpassert(ASSOCIATED(matrix_h_re))
1272 cpassert(ASSOCIATED(matrix_h_im))
1273 cpassert(ASSOCIATED(matrix_s_re))
1274 cpassert(ASSOCIATED(matrix_s_im))
1275 IF (solver_type /= ot_precond_solver_default .AND. &
1276 solver_type /= ot_precond_solver_inv_chol) THEN
1277 cpabort("Complex FULL_SINGLE_INVERSE supports only PRECOND_SOLVER DEFAULT/INVERSE_CHOLESKY")
1278 END IF
1279 CALL dbcsr_get_info(matrix_c_re, nfullrows_total=n, nfullcols_total=k)
1280 cpassert(n > 0 .AND. k > 0 .AND. k <= n)
1281
1282 CALL cp_fm_struct_create(fm_struct_mo, nrow_global=n, ncol_global=k, &
1283 context=preconditioner_env%ctxt, &
1284 para_env=preconditioner_env%para_env)
1285 CALL cp_fm_struct_create(fm_struct_square, nrow_global=n, ncol_global=n, &
1286 context=preconditioner_env%ctxt, &
1287 para_env=preconditioner_env%para_env)
1288 CALL cp_fm_create(matrix_c_re_fm, fm_struct_mo, name='complex preconditioner C real')
1289 CALL cp_fm_create(matrix_c_im_fm, fm_struct_mo, name='complex preconditioner C imaginary')
1290 CALL cp_fm_create(matrix_h_re_fm, fm_struct_square, name='complex preconditioner H real')
1291 CALL cp_fm_create(matrix_h_im_fm, fm_struct_square, name='complex preconditioner H imaginary')
1292 CALL cp_fm_create(matrix_s_re_fm, fm_struct_square, name='complex preconditioner S real')
1293 CALL cp_fm_create(matrix_s_im_fm, fm_struct_square, name='complex preconditioner S imaginary')
1294 CALL cp_fm_struct_release(fm_struct_mo)
1295 CALL cp_fm_struct_release(fm_struct_square)
1296 CALL copy_dbcsr_to_fm(matrix_c_re, matrix_c_re_fm)
1297 CALL copy_dbcsr_to_fm(matrix_c_im, matrix_c_im_fm)
1298 CALL copy_dbcsr_to_fm(matrix_h_re, matrix_h_re_fm)
1299 CALL copy_dbcsr_to_fm(matrix_h_im, matrix_h_im_fm)
1300 CALL copy_dbcsr_to_fm(matrix_s_re, matrix_s_re_fm)
1301 CALL copy_dbcsr_to_fm(matrix_s_im, matrix_s_im_fm)
1302
1303 CALL cp_cfm_create(matrix_c, matrix_c_re_fm%matrix_struct, name='complex preconditioner C')
1304 CALL cp_cfm_create(matrix_h, matrix_h_re_fm%matrix_struct, name='complex preconditioner H')
1305 CALL cp_cfm_create(matrix_s, matrix_s_re_fm%matrix_struct, name='complex preconditioner S')
1306 CALL cp_fm_to_cfm(matrix_c_re_fm, matrix_c_im_fm, matrix_c)
1307 CALL cp_fm_to_cfm(matrix_h_re_fm, matrix_h_im_fm, matrix_h)
1308 CALL cp_fm_to_cfm(matrix_s_re_fm, matrix_s_im_fm, matrix_s)
1309 CALL make_complex_full_single_inverse(preconditioner_env, matrix_c, matrix_h, matrix_s, &
1310 energy_gap)
1311 preconditioner_env%solver = solver_type
1312
1313 CALL cp_cfm_release(matrix_s)
1314 CALL cp_cfm_release(matrix_h)
1315 CALL cp_cfm_release(matrix_c)
1316 CALL cp_fm_release(matrix_s_im_fm)
1317 CALL cp_fm_release(matrix_s_re_fm)
1318 CALL cp_fm_release(matrix_h_im_fm)
1319 CALL cp_fm_release(matrix_h_re_fm)
1320 CALL cp_fm_release(matrix_c_im_fm)
1321 CALL cp_fm_release(matrix_c_re_fm)
1322
1324
1325! **************************************************************************************************
1326!> \brief Construct a complex FULL_S_INVERSE preconditioner.
1327!> \param preconditioner_env preconditioner storage
1328!> \param matrix_s_re real part of S(k)
1329!> \param matrix_s_im imaginary part of S(k)
1330!> \param solver_type requested inverse solver
1331! **************************************************************************************************
1332 SUBROUTINE make_preconditioner_complex_full_s_inverse(preconditioner_env, &
1333 matrix_s_re, matrix_s_im, solver_type)
1334
1335 TYPE(preconditioner_type) :: preconditioner_env
1336 TYPE(dbcsr_type), POINTER :: matrix_s_re, matrix_s_im
1337 INTEGER, INTENT(IN) :: solver_type
1338
1339 TYPE(cp_cfm_type) :: matrix_s
1340
1341 IF (solver_type /= ot_precond_solver_default .AND. &
1342 solver_type /= ot_precond_solver_inv_chol) THEN
1343 cpabort("Complex FULL_S_INVERSE supports only PRECOND_SOLVER DEFAULT/INVERSE_CHOLESKY")
1344 END IF
1345 CALL dbcsr_pair_to_cfm(matrix_s_re, matrix_s_im, preconditioner_env, &
1346 'complex preconditioner S', matrix_s)
1347 CALL make_complex_full_s_inverse(preconditioner_env, matrix_s)
1348 preconditioner_env%solver = solver_type
1349 CALL cp_cfm_release(matrix_s)
1350
1352
1353! **************************************************************************************************
1354!> \brief Construct a complex FULL_KINETIC preconditioner.
1355!> \param preconditioner_env preconditioner storage
1356!> \param matrix_t_re real part of T(k)
1357!> \param matrix_t_im imaginary part of T(k)
1358!> \param matrix_s_re real part of S(k)
1359!> \param matrix_s_im imaginary part of S(k)
1360!> \param energy_gap non-negative overlap shift
1361!> \param solver_type requested inverse solver
1362! **************************************************************************************************
1363 SUBROUTINE make_preconditioner_complex_full_kinetic(preconditioner_env, &
1364 matrix_t_re, matrix_t_im, &
1365 matrix_s_re, matrix_s_im, &
1366 energy_gap, solver_type)
1367
1368 TYPE(preconditioner_type) :: preconditioner_env
1369 TYPE(dbcsr_type), POINTER :: matrix_t_re, matrix_t_im, matrix_s_re, &
1370 matrix_s_im
1371 REAL(kind=dp), INTENT(IN) :: energy_gap
1372 INTEGER, INTENT(IN) :: solver_type
1373
1374 TYPE(cp_cfm_type) :: matrix_s, matrix_t
1375
1376 IF (solver_type /= ot_precond_solver_default .AND. &
1377 solver_type /= ot_precond_solver_inv_chol) THEN
1378 cpabort("Complex FULL_KINETIC supports only PRECOND_SOLVER DEFAULT/INVERSE_CHOLESKY")
1379 END IF
1380 CALL dbcsr_pair_to_cfm(matrix_t_re, matrix_t_im, preconditioner_env, &
1381 'complex preconditioner T', matrix_t)
1382 CALL dbcsr_pair_to_cfm(matrix_s_re, matrix_s_im, preconditioner_env, &
1383 'complex preconditioner S', matrix_s)
1384 CALL make_complex_full_kinetic(preconditioner_env, matrix_t, matrix_s, energy_gap)
1385 preconditioner_env%solver = solver_type
1386 CALL cp_cfm_release(matrix_s)
1387 CALL cp_cfm_release(matrix_t)
1388
1390
1391! **************************************************************************************************
1392!> \brief Copy a real/imaginary DBCSR pair to one distributed complex full matrix.
1393!> \param matrix_re real matrix component
1394!> \param matrix_im imaginary matrix component
1395!> \param preconditioner_env source of the distribution context
1396!> \param matrix_name matrix label
1397!> \param matrix complex output matrix
1398! **************************************************************************************************
1399 SUBROUTINE dbcsr_pair_to_cfm(matrix_re, matrix_im, preconditioner_env, matrix_name, matrix)
1400
1401 TYPE(dbcsr_type), POINTER :: matrix_re, matrix_im
1402 TYPE(preconditioner_type) :: preconditioner_env
1403 CHARACTER(LEN=*), INTENT(IN) :: matrix_name
1404 TYPE(cp_cfm_type), INTENT(OUT) :: matrix
1405
1406 INTEGER :: ncol, nrow
1407 REAL(dp), ALLOCATABLE :: buffer(:, :)
1408 TYPE(cp_fm_struct_type), POINTER :: fm_struct
1409 TYPE(cp_fm_type) :: matrix_im_fm, matrix_re_fm
1410
1411 NULLIFY (fm_struct)
1412 cpassert(ASSOCIATED(matrix_re))
1413 cpassert(ASSOCIATED(matrix_im))
1414 CALL dbcsr_get_info(matrix_re, nfullrows_total=nrow, nfullcols_total=ncol)
1415 CALL cp_fm_struct_create(fm_struct, nrow_global=nrow, ncol_global=ncol, &
1416 context=preconditioner_env%ctxt, &
1417 para_env=preconditioner_env%para_env)
1418 CALL cp_cfm_create(matrix, fm_struct, name=matrix_name)
1419 IF (preconditioner_env%para_env%num_pe == 1 .AND. &
1420 dbcsr_get_matrix_type(matrix_re) == dbcsr_type_no_symmetry .AND. &
1421 dbcsr_get_matrix_type(matrix_im) == dbcsr_type_no_symmetry) THEN
1422 ALLOCATE (buffer(nrow, ncol))
1423 CALL copy_dbcsr_to_local(matrix_re, buffer)
1424 matrix%local_data(1:nrow, 1:ncol) = cmplx(buffer, 0.0_dp, kind=dp)
1425 CALL copy_dbcsr_to_local(matrix_im, buffer)
1426 matrix%local_data(1:nrow, 1:ncol) = &
1427 matrix%local_data(1:nrow, 1:ncol) + cmplx(0.0_dp, buffer, kind=dp)
1428 ELSE
1429 CALL cp_fm_create(matrix_re_fm, fm_struct, name=trim(matrix_name)//' real')
1430 CALL cp_fm_create(matrix_im_fm, fm_struct, name=trim(matrix_name)//' imaginary')
1431 CALL copy_dbcsr_to_fm(matrix_re, matrix_re_fm)
1432 CALL copy_dbcsr_to_fm(matrix_im, matrix_im_fm)
1433 CALL cp_fm_to_cfm(matrix_re_fm, matrix_im_fm, matrix)
1434 CALL cp_fm_release(matrix_im_fm)
1435 CALL cp_fm_release(matrix_re_fm)
1436 END IF
1437 CALL cp_fm_struct_release(fm_struct)
1438
1439 END SUBROUTINE dbcsr_pair_to_cfm
1440
1441! **************************************************************************************************
1442!> \brief Allows for a restart of the preconditioner
1443!> depending on the method it purges all arrays or keeps them
1444!> \param qs_env ...
1445!> \param preconditioner ...
1446!> \param prec_type ...
1447!> \param nspins ...
1448! **************************************************************************************************
1449 SUBROUTINE restart_preconditioner(qs_env, preconditioner, prec_type, nspins)
1450
1451 TYPE(qs_environment_type), POINTER :: qs_env
1452 TYPE(preconditioner_p_type), DIMENSION(:), POINTER :: preconditioner
1453 INTEGER, INTENT(IN) :: prec_type, nspins
1454
1455 INTEGER :: ispin
1456 TYPE(cp_blacs_env_type), POINTER :: blacs_env
1457 TYPE(mp_para_env_type), POINTER :: para_env
1458
1459 NULLIFY (para_env, blacs_env)
1460 CALL get_qs_env(qs_env, para_env=para_env, blacs_env=blacs_env)
1461
1462 IF (ASSOCIATED(preconditioner)) THEN
1463 SELECT CASE (prec_type)
1466 ot_precond_full_single) ! these depend on the ks matrix
1467 DO ispin = 1, SIZE(preconditioner)
1468 CALL destroy_preconditioner(preconditioner(ispin)%preconditioner)
1469 DEALLOCATE (preconditioner(ispin)%preconditioner)
1470 END DO
1471 DEALLOCATE (preconditioner)
1473 ot_precond_full_single_inverse) ! these are 'independent'
1474 ! do nothing
1475 CASE DEFAULT
1476 cpabort("Unknown preconditioner type")
1477 END SELECT
1478 END IF
1479
1480 ! add an OT preconditioner if none is present
1481 IF (.NOT. ASSOCIATED(preconditioner)) THEN
1482 SELECT CASE (prec_type)
1485 ALLOCATE (preconditioner(nspins))
1486 CASE DEFAULT
1487 ALLOCATE (preconditioner(1))
1488 END SELECT
1489 DO ispin = 1, SIZE(preconditioner)
1490 ALLOCATE (preconditioner(ispin)%preconditioner)
1491 CALL init_preconditioner(preconditioner(ispin)%preconditioner, &
1492 para_env=para_env, &
1493 blacs_env=blacs_env)
1494 END DO
1495 END IF
1496
1497 END SUBROUTINE restart_preconditioner
1498
1499! **************************************************************************************************
1500!> \brief ...
1501!> \param qs_env ...
1502!> \param mos ...
1503!> \param matrix_ks ...
1504!> \param matrix_s ...
1505!> \param ot_preconditioner ...
1506!> \param prec_type ...
1507!> \param solver_type ...
1508!> \param energy_gap ...
1509!> \param nspins ...
1510!> \param has_unit_metric ...
1511!> \param convert_to_dbcsr ...
1512!> \param chol_type ...
1513!> \param full_mo_set ...
1514!> \param chebyshev_degree ...
1515!> \param low_rank_base ...
1516!> \param fermi_low_rank_max_rank ...
1517!> \param lattice_fft ...
1518!> \param lattice_fft_local_cells ...
1519! **************************************************************************************************
1520 SUBROUTINE prepare_preconditioner(qs_env, mos, matrix_ks, matrix_s, &
1521 ot_preconditioner, prec_type, solver_type, &
1522 energy_gap, nspins, has_unit_metric, &
1523 convert_to_dbcsr, chol_type, full_mo_set, chebyshev_degree, &
1524 low_rank_base, fermi_low_rank_max_rank, lattice_fft, &
1525 lattice_fft_local_cells)
1526
1527 TYPE(qs_environment_type), POINTER :: qs_env
1528 TYPE(mo_set_type), DIMENSION(:), INTENT(INOUT) :: mos
1529 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s
1530 TYPE(preconditioner_p_type), DIMENSION(:), POINTER :: ot_preconditioner
1531 INTEGER, INTENT(IN) :: prec_type, solver_type
1532 REAL(dp), INTENT(IN) :: energy_gap
1533 INTEGER, INTENT(IN) :: nspins
1534 LOGICAL, INTENT(IN), OPTIONAL :: has_unit_metric, convert_to_dbcsr
1535 INTEGER, INTENT(IN), OPTIONAL :: chol_type
1536 LOGICAL, INTENT(IN), OPTIONAL :: full_mo_set
1537 INTEGER, INTENT(IN), OPTIONAL :: chebyshev_degree, low_rank_base, &
1538 fermi_low_rank_max_rank, lattice_fft, &
1539 lattice_fft_local_cells
1540
1541 CHARACTER(LEN=*), PARAMETER :: routinen = 'prepare_preconditioner'
1542
1543 CHARACTER(LEN=default_string_length) :: msg
1544 INTEGER :: handle, icall, ispin, my_chebyshev_degree, my_fermi_low_rank_max_rank, &
1545 my_lattice_fft, my_lattice_fft_local_cells, my_low_rank_base, n_loops
1546 INTEGER, DIMENSION(5) :: nocc, norb
1547 INTEGER, DIMENSION(:), POINTER :: multiple_unit_cell
1548 LOGICAL :: do_co_rotate, do_kpoints, my_convert_to_dbcsr, my_full_mo_set, &
1549 my_has_unit_metric, use_covariant_full_all, use_mo_coeff_b
1550 TYPE(cp_blacs_env_type), POINTER :: blacs_env
1551 TYPE(cp_fm_type), POINTER :: mo_coeff
1552 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: kinetic
1553 TYPE(dbcsr_type), POINTER :: matrix_t, mo_coeff_b
1554 TYPE(dft_control_type), POINTER :: dft_control
1555 TYPE(mo_set_type), ALLOCATABLE, DIMENSION(:) :: canonical_mos
1556 TYPE(mp_para_env_type), POINTER :: para_env
1557 TYPE(section_vals_type), POINTER :: input, subsys_section
1558
1559 CALL timeset(routinen, handle)
1560 NULLIFY (matrix_t, mo_coeff_b, mo_coeff, kinetic, dft_control, para_env, blacs_env, &
1561 input, subsys_section, multiple_unit_cell)
1562 my_has_unit_metric = .false.
1563 IF (PRESENT(has_unit_metric)) my_has_unit_metric = has_unit_metric
1564 my_convert_to_dbcsr = .true.
1565 IF (PRESENT(convert_to_dbcsr)) my_convert_to_dbcsr = convert_to_dbcsr
1566 my_full_mo_set = .false.
1567 IF (PRESENT(full_mo_set)) my_full_mo_set = full_mo_set
1568 my_chebyshev_degree = 8
1569 IF (PRESENT(chebyshev_degree)) my_chebyshev_degree = chebyshev_degree
1570 my_low_rank_base = ot_low_rank_base_overlap
1571 IF (PRESENT(low_rank_base)) my_low_rank_base = low_rank_base
1572 my_fermi_low_rank_max_rank = 48
1573 IF (PRESENT(fermi_low_rank_max_rank)) my_fermi_low_rank_max_rank = fermi_low_rank_max_rank
1574 my_lattice_fft = ot_lattice_fft_off
1575 IF (PRESENT(lattice_fft)) my_lattice_fft = lattice_fft
1576 my_lattice_fft_local_cells = 0
1577 IF (PRESENT(lattice_fft_local_cells)) my_lattice_fft_local_cells = lattice_fft_local_cells
1578 IF (my_lattice_fft_local_cells < 0) cpabort("LATTICE_FFT_LOCAL_CELLS must not be negative")
1579 do_kpoints = .false.
1580 use_covariant_full_all = prec_type == ot_precond_full_all_covariant
1581
1582 CALL get_qs_env(qs_env, &
1583 dft_control=dft_control, &
1584 do_kpoints=do_kpoints, &
1585 input=input, &
1586 para_env=para_env, &
1587 blacs_env=blacs_env)
1588
1589 IF (my_low_rank_base == ot_low_rank_base_lattice_fft) THEN
1590 IF (prec_type /= ot_precond_fermi_low_rank) THEN
1591 cpabort("FERMI_LOW_RANK_BASE LATTICE_FFT requires PRECONDITIONER FERMI_LOW_RANK")
1592 END IF
1593 my_lattice_fft = ot_lattice_fft_on
1594 ELSE IF (my_low_rank_base /= ot_low_rank_base_overlap) THEN
1595 cpabort("Unknown FERMI_LOW_RANK_BASE setting")
1596 END IF
1597
1598 IF (my_lattice_fft /= ot_lattice_fft_off .AND. &
1599 my_lattice_fft /= ot_lattice_fft_on .AND. &
1600 my_lattice_fft /= ot_lattice_fft_auto) THEN
1601 cpabort("Unknown LATTICE_FFT setting")
1602 END IF
1603 IF (my_lattice_fft_local_cells > 0 .AND. my_lattice_fft /= ot_lattice_fft_on) THEN
1604 cpabort("LATTICE_FFT_LOCAL_CELLS requires LATTICE_FFT ON")
1605 END IF
1606 IF (my_lattice_fft /= ot_lattice_fft_off) THEN
1607 SELECT CASE (prec_type)
1610 CASE DEFAULT
1611 cpabort("LATTICE_FFT does not support the selected OT preconditioner")
1612 END SELECT
1613 IF (solver_type /= ot_precond_solver_default .AND. &
1614 solver_type /= ot_precond_solver_inv_chol) THEN
1615 cpabort("LATTICE_FFT requires PRECOND_SOLVER DEFAULT or INVERSE_CHOLESKY")
1616 END IF
1617 IF (my_has_unit_metric .OR. do_kpoints .OR. .NOT. ASSOCIATED(input)) THEN
1618 IF (my_lattice_fft == ot_lattice_fft_on) THEN
1619 cpabort("LATTICE_FFT ON requires a non-unit metric in a real Gamma-point calculation")
1620 END IF
1621 my_lattice_fft = ot_lattice_fft_off
1622 END IF
1623 END IF
1624 IF (my_lattice_fft /= ot_lattice_fft_off) THEN
1625 subsys_section => section_vals_get_subs_vals(input, "SUBSYS")
1626 IF (.NOT. ASSOCIATED(subsys_section)) THEN
1627 IF (my_lattice_fft == ot_lattice_fft_on) cpabort("LATTICE_FFT ON requires a SUBSYS section")
1628 my_lattice_fft = ot_lattice_fft_off
1629 END IF
1630 IF (my_lattice_fft /= ot_lattice_fft_off) THEN
1631 CALL section_vals_val_get(subsys_section, "CELL%MULTIPLE_UNIT_CELL", &
1632 i_vals=multiple_unit_cell)
1633 IF (.NOT. ASSOCIATED(multiple_unit_cell)) THEN
1634 IF (my_lattice_fft == ot_lattice_fft_on) THEN
1635 cpabort("LATTICE_FFT ON requires a valid CELL%MULTIPLE_UNIT_CELL")
1636 END IF
1637 my_lattice_fft = ot_lattice_fft_off
1638 ELSE
1639 IF (SIZE(multiple_unit_cell) /= 3 .OR. any(multiple_unit_cell <= 0)) THEN
1640 IF (my_lattice_fft == ot_lattice_fft_on) THEN
1641 cpabort("LATTICE_FFT ON requires a valid CELL%MULTIPLE_UNIT_CELL")
1642 END IF
1643 my_lattice_fft = ot_lattice_fft_off
1644 ELSE IF (product(multiple_unit_cell) <= 1 .AND. my_lattice_fft == ot_lattice_fft_on) THEN
1645 cpabort("LATTICE_FFT ON requires an explicitly replicated multiple unit cell")
1646 END IF
1647 IF (my_lattice_fft_local_cells >= product(multiple_unit_cell)) THEN
1648 cpabort("LATTICE_FFT_LOCAL_CELLS must be smaller than the replicated cell count")
1649 END IF
1650 END IF
1651 END IF
1652 END IF
1653
1654 IF (dft_control%qs_control%semi_empirical .OR. dft_control%qs_control%dftb .OR. &
1655 dft_control%qs_control%xtb) THEN
1656 IF (prec_type == ot_precond_full_kinetic) THEN
1657 msg = "Full_kinetic not available for semi-empirical methods"
1658 cpabort(trim(msg))
1659 END IF
1660 matrix_t => matrix_s(1)%matrix
1661 ELSE
1662 cpassert(.NOT. my_has_unit_metric)
1663 CALL get_qs_env(qs_env, kinetic=kinetic)
1664 matrix_t => kinetic(1)%matrix
1665 END IF
1666
1667 ! use full set of MOs or just occupied MOs
1668 nocc = 0
1669 norb = 0
1670 IF (my_full_mo_set) THEN
1671 DO ispin = 1, nspins
1672 CALL get_mo_set(mo_set=mos(ispin), homo=nocc(ispin), nmo=norb(ispin))
1673 CALL set_mo_set(mo_set=mos(ispin), homo=norb(ispin))
1674 END DO
1675 END IF
1676 !determines how often make preconditioner is called, spin dependent methods have to be called twice
1677 n_loops = 1
1678 IF (prec_type == ot_precond_full_single_inverse) n_loops = nspins
1679 ! check whether we need the ev and rotate the MOs
1680 SELECT CASE (prec_type)
1681 CASE (ot_precond_full_all)
1682 ! FULL_ALL needs canonical occupied eigenstates and one preconditioner per spin.
1683 n_loops = nspins
1684
1685 do_co_rotate = ASSOCIATED(qs_env%mo_derivs)
1686 DO ispin = 1, nspins
1687 CALL get_mo_set(mo_set=mos(ispin), mo_coeff_b=mo_coeff_b, mo_coeff=mo_coeff)
1688 use_mo_coeff_b = mos(ispin)%use_mo_coeff_b
1689 IF (use_mo_coeff_b .AND. do_co_rotate) THEN
1690 CALL calculate_subspace_eigenvalues(mo_coeff_b, matrix_ks(ispin)%matrix, &
1691 do_rotation=.true., &
1692 co_rotate=qs_env%mo_derivs(ispin)%matrix, &
1693 para_env=para_env, &
1694 blacs_env=blacs_env)
1695 ELSE IF (use_mo_coeff_b) THEN
1696 CALL calculate_subspace_eigenvalues(mo_coeff_b, matrix_ks(ispin)%matrix, &
1697 do_rotation=.true., &
1698 para_env=para_env, &
1699 blacs_env=blacs_env)
1700 ELSE
1701 CALL calculate_subspace_eigenvalues(mo_coeff, matrix_ks(ispin)%matrix, &
1702 do_rotation=.true.)
1703 END IF
1704 END DO
1706 ! The common-reference low-rank operator depends only on the occupied projector.
1707 ! Keep the input gauge unchanged so that OT ROTATION remains well defined.
1708 n_loops = nspins
1710 ! Diagonalize a private occupied copy. The physical MOs and their
1711 ! derivatives remain in the user-selected gauge.
1712 n_loops = nspins
1713 ALLOCATE (canonical_mos(nspins))
1714 DO ispin = 1, nspins
1715 CALL duplicate_mo_set(canonical_mos(ispin), mos(ispin))
1716 CALL get_mo_set(mo_set=canonical_mos(ispin), mo_coeff_b=mo_coeff_b, mo_coeff=mo_coeff)
1717 IF (canonical_mos(ispin)%use_mo_coeff_b) THEN
1718 CALL copy_dbcsr_to_fm(mo_coeff_b, mo_coeff)
1719 canonical_mos(ispin)%use_mo_coeff_b = .false.
1720 END IF
1721 CALL canonicalize_full_all_copy(ot_preconditioner(ispin)%preconditioner, &
1722 mo_coeff, matrix_ks(ispin)%matrix, &
1723 canonical_mos(ispin)%homo)
1724 END DO
1725 CASE DEFAULT
1726 ! No need to rotate the MOs
1727 END SELECT
1728
1729 ! check whether we have a preconditioner
1730 SELECT CASE (prec_type)
1731 CASE (ot_precond_none)
1732 DO ispin = 1, SIZE(ot_preconditioner)
1733 ot_preconditioner(ispin)%preconditioner%in_use = 0
1734 END DO
1735 CASE DEFAULT
1736 DO icall = 1, n_loops
1737 ot_preconditioner(icall)%preconditioner%polynomial_degree = my_chebyshev_degree
1738 ot_preconditioner(icall)%preconditioner%low_rank_base = my_low_rank_base
1739 ot_preconditioner(icall)%preconditioner%lattice_fft_mode = my_lattice_fft
1740 ot_preconditioner(icall)%preconditioner%lattice_local_cells = my_lattice_fft_local_cells
1741 ot_preconditioner(icall)%preconditioner%lattice_dims = [1, 1, 1]
1742 IF (my_lattice_fft /= ot_lattice_fft_off) THEN
1743 cpassert(ASSOCIATED(multiple_unit_cell))
1744 ot_preconditioner(icall)%preconditioner%lattice_dims = multiple_unit_cell
1745 END IF
1746 ot_preconditioner(icall)%preconditioner%spectral_max_rank = my_fermi_low_rank_max_rank
1747 IF (my_has_unit_metric) THEN
1748 IF (use_covariant_full_all) THEN
1749 CALL make_preconditioner(ot_preconditioner(icall)%preconditioner, &
1750 prec_type, solver_type, &
1751 matrix_h=matrix_ks(icall)%matrix, &
1752 mo_set=canonical_mos(icall), &
1753 energy_gap=energy_gap, &
1754 convert_precond_to_dbcsr=my_convert_to_dbcsr)
1755 ELSE
1756 CALL make_preconditioner(ot_preconditioner(icall)%preconditioner, &
1757 prec_type, solver_type, &
1758 matrix_h=matrix_ks(icall)%matrix, &
1759 mo_set=mos(icall), &
1760 energy_gap=energy_gap, &
1761 convert_precond_to_dbcsr=my_convert_to_dbcsr)
1762 END IF
1763 ELSE
1764 IF (use_covariant_full_all) THEN
1765 CALL make_preconditioner(ot_preconditioner(icall)%preconditioner, &
1766 prec_type, solver_type, &
1767 matrix_h=matrix_ks(icall)%matrix, &
1768 matrix_s=matrix_s(1)%matrix, matrix_t=matrix_t, &
1769 mo_set=canonical_mos(icall), &
1770 energy_gap=energy_gap, &
1771 convert_precond_to_dbcsr=my_convert_to_dbcsr, chol_type=chol_type)
1772 ELSE
1773 CALL make_preconditioner(ot_preconditioner(icall)%preconditioner, &
1774 prec_type, solver_type, &
1775 matrix_h=matrix_ks(icall)%matrix, &
1776 matrix_s=matrix_s(1)%matrix, matrix_t=matrix_t, &
1777 mo_set=mos(icall), energy_gap=energy_gap, &
1778 convert_precond_to_dbcsr=my_convert_to_dbcsr, chol_type=chol_type)
1779 END IF
1780 END IF
1781 END DO
1782 END SELECT
1783
1784 IF (ALLOCATED(canonical_mos)) THEN
1785 DO ispin = 1, SIZE(canonical_mos)
1786 CALL deallocate_mo_set(canonical_mos(ispin))
1787 END DO
1788 DEALLOCATE (canonical_mos)
1789 END IF
1790
1791 ! reset homo values
1792 IF (my_full_mo_set) THEN
1793 DO ispin = 1, nspins
1794 CALL set_mo_set(mo_set=mos(ispin), homo=nocc(ispin))
1795 END DO
1796 END IF
1797
1798 CALL timestop(handle)
1799
1800 END SUBROUTINE prepare_preconditioner
1801
1802! **************************************************************************************************
1803!> \brief Canonicalize an occupied copy and retain the occupied-space rotation.
1804!> \param preconditioner_env ...
1805!> \param mo_coeff ...
1806!> \param matrix_h ...
1807!> \param nocc ...
1808! **************************************************************************************************
1809 SUBROUTINE canonicalize_full_all_copy(preconditioner_env, mo_coeff, matrix_h, nocc)
1810
1811 TYPE(preconditioner_type) :: preconditioner_env
1812 TYPE(cp_fm_type), POINTER :: mo_coeff
1813 TYPE(dbcsr_type), POINTER :: matrix_h
1814 INTEGER, INTENT(IN) :: nocc
1815
1816 CHARACTER(LEN=*), PARAMETER :: routinen = 'canonicalize_full_all_copy'
1817
1818 INTEGER :: handle, nao
1819 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: occupied_evals
1820 TYPE(cp_fm_struct_type), POINTER :: fm_struct
1821 TYPE(cp_fm_type) :: canonical_coeff, h_block, weighted_coeff
1822
1823 CALL timeset(routinen, handle)
1824 CALL cp_fm_get_info(mo_coeff, nrow_global=nao)
1825 cpassert(nocc > 0)
1826
1827 CALL cp_fm_struct_create(fm_struct, nrow_global=nao, ncol_global=nocc, &
1828 context=preconditioner_env%ctxt, &
1829 para_env=preconditioner_env%para_env)
1830 CALL cp_fm_create(canonical_coeff, fm_struct, name=routinen)
1831 CALL cp_fm_create(weighted_coeff, fm_struct, name=routinen)
1832 CALL cp_fm_struct_release(fm_struct)
1833 CALL cp_fm_to_fm(mo_coeff, canonical_coeff, nocc)
1834
1835 CALL cp_fm_struct_create(fm_struct, nrow_global=nocc, ncol_global=nocc, &
1836 context=preconditioner_env%ctxt, &
1837 para_env=preconditioner_env%para_env)
1838 CALL cp_fm_create(h_block, fm_struct, name=routinen)
1839 IF (ASSOCIATED(preconditioner_env%occ_rotation_dbcsr)) THEN
1840 CALL dbcsr_release_p(preconditioner_env%occ_rotation_dbcsr)
1841 END IF
1842 IF (ASSOCIATED(preconditioner_env%occ_rotation)) THEN
1843 CALL cp_fm_release(preconditioner_env%occ_rotation)
1844 DEALLOCATE (preconditioner_env%occ_rotation)
1845 END IF
1846 ALLOCATE (preconditioner_env%occ_rotation)
1847 CALL cp_fm_create(preconditioner_env%occ_rotation, fm_struct, name="occupied rotation")
1848 CALL cp_fm_struct_release(fm_struct)
1849
1850 CALL cp_dbcsr_sm_fm_multiply(matrix_h, canonical_coeff, weighted_coeff, nocc)
1851 CALL parallel_gemm('T', 'N', nocc, nocc, nao, 1.0_dp, canonical_coeff, &
1852 weighted_coeff, 0.0_dp, h_block)
1853 ALLOCATE (occupied_evals(nocc))
1854 CALL choose_eigv_solver(h_block, preconditioner_env%occ_rotation, occupied_evals)
1855
1856 CALL parallel_gemm('N', 'N', nao, nocc, nocc, 1.0_dp, canonical_coeff, &
1857 preconditioner_env%occ_rotation, 0.0_dp, weighted_coeff)
1858 CALL cp_fm_to_fm(weighted_coeff, mo_coeff, nocc, 1, 1)
1859
1860 DEALLOCATE (occupied_evals)
1861 CALL cp_fm_release(canonical_coeff)
1862 CALL cp_fm_release(h_block)
1863 CALL cp_fm_release(weighted_coeff)
1864 CALL timestop(handle)
1865
1866 END SUBROUTINE canonicalize_full_all_copy
1867
1868END MODULE preconditioner
static GRID_HOST_DEVICE int modulo(int a, int m)
Equivalent of Fortran's MODULO, which always return a positive number. https://gcc....
methods related to the blacs parallel environment
Basic linear algebra operations for complex full matrices.
subroutine, public cp_cfm_triangular_multiply(triangular_matrix, matrix_b, side, transa_tr, invert_tr, uplo_tr, unit_diag_tr, n_rows, n_cols, alpha)
Multiplies in place by a triangular matrix: matrix_b = alpha op(triangular_matrix) matrix_b or (if si...
subroutine, public cp_cfm_triangular_invert(matrix_a, uplo, info_out)
Inverts a triangular matrix.
various cholesky decomposition related routines
subroutine, public cp_cfm_cholesky_decompose(matrix, n, info_out)
Used to replace a symmetric positive definite matrix M with its Cholesky decomposition U: M = U^T * U...
used for collecting diagonalization schemes available for cp_cfm_type
Definition cp_cfm_diag.F:14
subroutine, public cp_cfm_heevd_local(matrix, vectors, info, ws)
Local LAPACK ZHEEVD with prepared scratch, without allocation or communication.
subroutine, public cp_cfm_local_workspace_prepare(ws, n, a, b)
Query LAPACK once per matrix order/solver and retain private scratch across k points.
Represents a complex full matrix distributed on many processors.
subroutine, public cp_cfm_release(matrix)
Releases a full matrix.
subroutine, public cp_fm_to_cfm(msourcer, msourcei, mtarget)
Construct a complex full matrix by taking its real and imaginary parts from two separate real-value f...
subroutine, public cp_cfm_create(matrix, matrix_struct, name, nrow, ncol, set_zero)
Creates a new full matrix with the given structure.
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
subroutine, public dbcsr_release_p(matrix)
...
character function, public dbcsr_get_matrix_type(matrix)
...
subroutine, public dbcsr_get_info(matrix, nblkrows_total, nblkcols_total, nfullrows_total, nfullcols_total, nblkrows_local, nblkcols_local, nfullrows_local, nfullcols_local, my_prow, my_pcol, local_rows, local_cols, proc_row_dist, proc_col_dist, row_blk_size, col_blk_size, row_blk_offset, col_blk_offset, distribution, name, matrix_type, group)
...
DBCSR operations in CP2K.
subroutine, public cp_dbcsr_sm_fm_multiply(matrix, fm_in, fm_out, ncol, alpha, beta)
multiply a dbcsr with a fm matrix
subroutine, public copy_dbcsr_to_fm(matrix, fm, plan, check_plan)
Copy a DBCSR matrix to a BLACS matrix.
subroutine, public copy_dbcsr_to_local(matrix, dense)
Load a single-rank, nonsymmetric matrix without redistribution.
Basic linear algebra operations for full matrices.
subroutine, public cp_fm_uplo_to_full(matrix, work, uplo)
given a triangular matrix according to uplo, computes the corresponding full matrix
various cholesky decomposition related routines
subroutine, public cp_fm_cholesky_invert(matrix, n, info_out, work)
Replace the upper Cholesky factor by the upper triangle of the inverse. When work is provided,...
subroutine, public cp_fm_cholesky_decompose(matrix, n, info_out)
used to replace a symmetric positive def. matrix M with its cholesky decomposition U: M = U^T * U,...
used for collecting some of the diagonalization schemes available for cp_fm_type. cp_fm_power also mo...
Definition cp_fm_diag.F:17
logical, save, public direct_generalized_diagonalization
Definition cp_fm_diag.F:106
subroutine, public choose_eigv_solver(matrix, eigenvectors, eigenvalues, info)
Choose the Eigensolver depending on which library is available ELPA seems to be unstable for small sy...
Definition cp_fm_diag.F:262
logical function, public diag_check_requested()
Return whether diagonalization checks should be performed.
Definition cp_fm_diag.F:316
logical, save, public diag_lib_explicit
Definition cp_fm_diag.F:96
represent the structure of a full matrix
subroutine, public cp_fm_struct_create(fmstruct, para_env, context, nrow_global, ncol_global, nrow_block, ncol_block, descriptor, first_p_pos, local_leading_dimension, template_fmstruct, square_blocks, force_block)
allocates and initializes a full matrix structure
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
represent a full matrix distributed on many processors
Definition cp_fm_types.F:15
subroutine, public cp_fm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, nrow_locals, ncol_locals, matrix_struct, para_env)
returns all kind of information about the full matrix
subroutine, public cp_fm_set_all(matrix, alpha, beta)
set all elements of a matrix to the same value, and optionally the diagonal to a different one
subroutine, public cp_fm_create(matrix, matrix_struct, name, nrow, ncol, set_zero)
creates a new full matrix with the given structure
various routines to log and control the output. The idea is that decisions about where to log should ...
integer function, public cp_logger_get_default_io_unit(logger)
returns the unit nr for the ionode (-1 on all other processors) skips as well checks if the procs cal...
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
collects all constants needed in input so that they can be used without circular dependencies
integer, parameter, public ot_low_rank_base_overlap
integer, parameter, public ot_precond_full_all_covariant
integer, parameter, public ot_low_rank_base_lattice_fft
integer, parameter, public ot_lattice_fft_on
integer, parameter, public ot_precond_full_kinetic
integer, parameter, public cholesky_reduce
integer, parameter, public ot_lattice_fft_auto
integer, parameter, public ot_precond_solver_default
integer, parameter, public ot_lattice_fft_off
integer, parameter, public ot_precond_full_single
integer, parameter, public ot_precond_solver_inv_chol
integer, parameter, public ot_precond_none
integer, parameter, public ot_precond_full_single_inverse
integer, parameter, public ot_precond_fermi_low_rank
integer, parameter, public ot_precond_s_inverse
integer, parameter, public ot_precond_solver_update
integer, parameter, public ot_precond_full_all
objects that represent the structure of input sections and the data contained in an input section
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Defines the basic variable types.
Definition kinds.F:23
integer, parameter, public dp
Definition kinds.F:34
integer, parameter, public default_string_length
Definition kinds.F:57
Block-circulant inverse operators for replicated-cell OT preconditioners.
real(kind=dp), parameter, public lattice_fft_auto_max_cost_ratio
subroutine, public build_lattice_local_correction(operator_matrix, lattice_dims, local_cell_count, correction_indices, operator_columns, coarse_inverse, corrected_projection_error, info)
Build a balanced local correction for selected non-circulant cell blocks.
subroutine, public build_lattice_inverse(operator_matrix, lattice_dims, inverse_k, projection_error, used_fft, info)
Project a cell-major SPD operator onto lattice translations and invert its Fourier blocks.
pure real(kind=dp) function, public lattice_local_storage_ratio(lattice_dims, operator_count, local_cell_count)
Estimate balanced-local-correction storage relative to one dense real operator.
pure real(kind=dp) function, public lattice_fft_storage_ratio(lattice_dims, operator_count)
Estimate reciprocal-block storage relative to one dense real operator.
pure logical function, public lattice_fft_auto_select(lattice_dims, block_size, rhs_count, operator_count, num_pe, projection_error, used_fft, info)
Conservative error-, work-, and storage-aware AUTO selector.
subroutine, public lattice_grid(lattice_dims, index_to_cell, xkp)
Enumerate a regular direct/reciprocal lattice with the first direction varying fastest.
integer, parameter, public lattice_fft_auto_min_cells
subroutine, public build_lattice_inverse_from_blocks(blocks, lattice_dims, inverse_k, used_fft, info)
Transform and invert already projected real-space lattice blocks.
pure real(kind=dp) function, public lattice_local_cost_ratio(lattice_dims, local_cell_count)
Estimate balanced-local-correction work relative to a dense inverse application.
pure real(kind=dp) function, public lattice_fft_cost_ratio(lattice_dims, block_size, rhs_count, num_pe)
Estimate lattice-FFT application work relative to a dense inverse application.
real(kind=dp), parameter, public lattice_fft_auto_max_storage_ratio
Small algebraic helpers for the rotationally covariant low-rank OT preconditioner.
pure elemental real(kind=dp) function, public low_rank_inverse_weight(eigenvalue, reference, energy_gap)
Spectral inverse weight relative to a common occupied reference level.
pure integer function, public low_rank_select_rank(eigenvalues, nocc, max_rank, degeneracy_tolerance)
Select a bounded spectral rank without cutting a degenerate boundary manifold.
Interface to the message passing library MPI.
basic linear algebra operations for full matrixes
computes preconditioners, and implements methods to apply them currently used in qs_ot
subroutine, public apply_preconditioner_dbcsr_complex(preconditioner_env, matrix_in_re, matrix_in_im, matrix_out_re, matrix_out_im)
Apply a complex orbital preconditioner to paired real/imaginary DBCSR matrices.
subroutine, public apply_preconditioner_fm(preconditioner_env, matrix_in, matrix_out)
applies a previously created preconditioner to a full matrix
subroutine, public apply_preconditioner_cfm_complex(preconditioner_env, matrix_in, matrix_out)
Apply a complex k-point orbital preconditioner.
subroutine, public apply_preconditioner_dbcsr(preconditioner_env, matrix_in, matrix_out)
...
subroutine, public apply_preconditioner_complex_batch(preconditioners, matrix_in, matrix_out, handled, output_scale)
Apply single-rank complex preconditioners in bounded OpenMP batches. Packing and DBCSR writeback stay...
computes preconditioners, and implements methods to apply them currently used in qs_ot
subroutine, public make_complex_full_single_inverse(preconditioner_env, matrix_c0, matrix_h, matrix_s, energy_gap)
Build a gauge-covariant FULL_SINGLE_INVERSE operator for a complex k-point channel.
subroutine, public make_complex_full_kinetic(preconditioner_env, matrix_t, matrix_s, energy_gap)
Build the inverse complex kinetic-plus-overlap preconditioner.
subroutine, public make_complex_full_s_inverse(preconditioner_env, matrix_s)
Build the inverse complex overlap preconditioner.
subroutine, public make_complex_full_all_covariant(preconditioner_env, matrix_c0, matrix_h, matrix_s, energy_gap, occupation_signature)
Build the gauge-covariant FULL_ALL inverse for a complex k-point channel.
subroutine, public make_complex_full_single(preconditioner_env, matrix_h, matrix_s, energy_homo, energy_gap)
Build the complex spectral FULL_SINGLE preconditioner.
subroutine, public prepare_complex_full_all(preconditioner_env, matrix_c0, matrix_h, matrix_s, matrix_tmp, matrix_s_copy, c0_evals, energy_gap)
Prepare the FULL_ALL eigenproblem using the same construction for every solver.
subroutine, public make_complex_fermi_low_rank(preconditioner_env, matrix_c0, matrix_h, matrix_s, energy_gap, max_rank, spectral_window, degeneracy_tolerance)
Build the bounded low-rank spectral inverse for a complex k-point channel.
subroutine, public make_preconditioner_matrix(preconditioner_env, matrix_h, matrix_s, matrix_t, mo_coeff, energy_homo, eigenvalues_ot, energy_gap, my_solver_type)
...
subroutine, public make_complex_full_all(preconditioner_env, matrix_c0, matrix_h, matrix_s, c0_evals, energy_gap)
Build the state-selective FULL_ALL operator for a complex k-point channel. The occupied/reference sub...
solves the preconditioner, contains to utility function for fm<->dbcsr transfers, should be moved soo...
subroutine, public transfer_dbcsr_to_fm(dbcsr_matrix, fm_matrix, para_env, context)
transfers a dbcsr to a full matrix
subroutine, public solve_preconditioner(my_solver_type, preconditioner_env, matrix_s, matrix_h)
...
subroutine, public transfer_fm_to_dbcsr(fm_matrix, dbcsr_matrix, template_mat)
transfers a full matrix to a dbcsr
types of preconditioners
subroutine, public init_preconditioner(preconditioner_env, para_env, blacs_env)
...
subroutine, public destroy_preconditioner(preconditioner_env)
...
computes preconditioners, and implements methods to apply them currently used in qs_ot
subroutine, public make_preconditioner_complex_full_s_inverse(preconditioner_env, matrix_s_re, matrix_s_im, solver_type)
Construct a complex FULL_S_INVERSE preconditioner.
subroutine, public make_preconditioner_complex_full_single(preconditioner_env, matrix_h_re, matrix_h_im, matrix_s_re, matrix_s_im, mo_set, energy_gap, solver_type)
Construct a complex FULL_SINGLE preconditioner from H(k) and S(k).
subroutine, public dbcsr_pair_to_cfm(matrix_re, matrix_im, preconditioner_env, matrix_name, matrix)
Copy a real/imaginary DBCSR pair to one distributed complex full matrix.
subroutine, public make_preconditioner_complex_full_single_inverse(preconditioner_env, matrix_c_re, matrix_c_im, matrix_h_re, matrix_h_im, matrix_s_re, matrix_s_im, energy_gap, solver_type)
Construct a complex FULL_SINGLE_INVERSE preconditioner without discarding Im(H,S,C).
subroutine, public make_preconditioner_complex_full_all_batch(preconditioners, matrix_c, matrix_h, matrix_s, energy_gap, solver_type, handled)
Build local FULL_ALL preconditioners from cached operators in bounded batches. Only the independent e...
subroutine, public make_preconditioner_complex_full_all_covariant(preconditioner_env, matrix_c_re, matrix_c_im, matrix_h_re, matrix_h_im, matrix_s_re, matrix_s_im, energy_gap, solver_type, occupation_signature)
Construct rotation-covariant FULL_ALL for one complex k-point channel.
subroutine, public restart_preconditioner(qs_env, preconditioner, prec_type, nspins)
Allows for a restart of the preconditioner depending on the method it purges all arrays or keeps them...
subroutine, public prepare_preconditioner(qs_env, mos, matrix_ks, matrix_s, ot_preconditioner, prec_type, solver_type, energy_gap, nspins, has_unit_metric, convert_to_dbcsr, chol_type, full_mo_set, chebyshev_degree, low_rank_base, fermi_low_rank_max_rank, lattice_fft, lattice_fft_local_cells)
...
subroutine, public make_preconditioner_complex_full_all(preconditioner_env, matrix_c_re, matrix_c_im, matrix_h_re, matrix_h_im, matrix_s_re, matrix_s_im, energy_gap, solver_type)
Construct FULL_ALL directly from one complex H(k), S(k), and C(k) channel.
subroutine, public make_preconditioner_complex_fermi_low_rank(preconditioner_env, matrix_c_re, matrix_c_im, matrix_h_re, matrix_h_im, matrix_s_re, matrix_s_im, energy_gap, max_rank, solver_type)
Construct FERMI_LOW_RANK directly for one complex k-point channel.
subroutine, public make_preconditioner_complex_full_kinetic(preconditioner_env, matrix_t_re, matrix_t_im, matrix_s_re, matrix_s_im, energy_gap, solver_type)
Construct a complex FULL_KINETIC preconditioner.
subroutine, public make_preconditioner(preconditioner_env, precon_type, solver_type, matrix_h, matrix_s, matrix_t, mo_set, energy_gap, convert_precond_to_dbcsr, chol_type)
...
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, mimic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, sab_cneo, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, matrix_vhxc, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, xcint_weights, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, rhoz_cneo_set, ecoul_1c, rho0_s_rs, rho0_s_gs, rhoz_cneo_s_rs, rhoz_cneo_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, harris_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs, do_rixs, tb_tblite)
Get the QUICKSTEP environment.
collects routines that perform operations directly related to MOs
Definition and initialisation of the mo data type.
Definition qs_mo_types.F:22
subroutine, public duplicate_mo_set(mo_set_new, mo_set_old)
allocate a new mo_set, and copy the old data
subroutine, public set_mo_set(mo_set, maxocc, homo, lfomo, nao, nelectron, n_el_f, nmo, eigenvalues, occupation_numbers, uniform_occupation, kts, mu, flexible_electron_count)
Set the components of a MO set data structure.
subroutine, public deallocate_mo_set(mo_set)
Deallocate a wavefunction data structure.
subroutine, public get_mo_set(mo_set, maxocc, homo, lfomo, nao, nelectron, n_el_f, nmo, eigenvalues, occupation_numbers, mo_coeff, mo_coeff_b, uniform_occupation, kts, mu, flexible_electron_count, cmo_coeff)
Get the components of a MO set data structure.
represent a blacs multidimensional parallel environment (for the mpi corrispective see cp_paratypes/m...
Represent a complex full matrix.
keeps the information about the structure of a full matrix
just to build arrays of pointers to matrices
represent a full matrix
type of a logger, at the moment it contains just a print level starting at which level it should be l...
stores all the informations relevant to an mpi environment