Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief 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 : ! **************************************************************************************************
15 : MODULE preconditioner
16 : USE cp_blacs_env, ONLY: cp_blacs_env_type
17 : USE cp_cfm_basic_linalg, ONLY: cp_cfm_triangular_invert,&
18 : cp_cfm_triangular_multiply
19 : USE cp_cfm_cholesky, ONLY: cp_cfm_cholesky_decompose
20 : USE cp_cfm_diag, ONLY: cp_cfm_diag_workspace_type,&
21 : cp_cfm_heevd_local,&
22 : cp_cfm_local_workspace_prepare
23 : USE cp_cfm_types, ONLY: cp_cfm_create,&
24 : cp_cfm_release,&
25 : cp_cfm_to_cfm,&
26 : cp_cfm_type,&
27 : cp_fm_to_cfm
28 : USE cp_control_types, ONLY: dft_control_type
29 : USE cp_dbcsr_api, ONLY: dbcsr_get_info,&
30 : dbcsr_get_matrix_type,&
31 : dbcsr_p_type,&
32 : dbcsr_release_p,&
33 : dbcsr_type,&
34 : dbcsr_type_no_symmetry
35 : USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
36 : copy_dbcsr_to_local,&
37 : cp_dbcsr_sm_fm_multiply
38 : USE cp_fm_basic_linalg, ONLY: cp_fm_uplo_to_full
39 : USE cp_fm_cholesky, ONLY: cp_fm_cholesky_decompose,&
40 : cp_fm_cholesky_invert
41 : USE cp_fm_diag, ONLY: choose_eigv_solver,&
42 : diag_check_requested,&
43 : diag_lib_explicit,&
44 : direct_generalized_diagonalization
45 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
46 : cp_fm_struct_release,&
47 : cp_fm_struct_type
48 : USE cp_fm_types, ONLY: cp_fm_create,&
49 : cp_fm_get_info,&
50 : cp_fm_p_type,&
51 : cp_fm_release,&
52 : cp_fm_set_all,&
53 : cp_fm_to_fm,&
54 : cp_fm_type
55 : USE cp_log_handling, ONLY: cp_get_default_logger,&
56 : cp_logger_get_default_io_unit,&
57 : cp_logger_type
58 : USE input_constants, ONLY: &
59 : cholesky_reduce, ot_lattice_fft_auto, ot_lattice_fft_off, ot_lattice_fft_on, &
60 : ot_low_rank_base_lattice_fft, ot_low_rank_base_overlap, ot_precond_fermi_low_rank, &
61 : ot_precond_full_all, ot_precond_full_all_covariant, ot_precond_full_kinetic, &
62 : ot_precond_full_single, ot_precond_full_single_inverse, ot_precond_none, &
63 : ot_precond_s_inverse, ot_precond_solver_default, ot_precond_solver_inv_chol, &
64 : ot_precond_solver_update
65 : USE input_section_types, ONLY: section_vals_get_subs_vals,&
66 : section_vals_type,&
67 : section_vals_val_get
68 : USE kinds, ONLY: default_string_length,&
69 : dp
70 : USE lattice_preconditioner_operator, ONLY: &
71 : build_lattice_inverse, build_lattice_inverse_from_blocks, build_lattice_local_correction, &
72 : lattice_fft_auto_max_cost_ratio, lattice_fft_auto_max_storage_ratio, &
73 : lattice_fft_auto_min_cells, lattice_fft_auto_select, lattice_fft_cost_ratio, &
74 : lattice_fft_storage_ratio, lattice_grid, lattice_local_cost_ratio, &
75 : lattice_local_storage_ratio
76 : USE low_rank_preconditioner_model, ONLY: low_rank_inverse_weight,&
77 : low_rank_select_rank
78 : USE message_passing, ONLY: mp_para_env_type
79 : USE parallel_gemm_api, ONLY: parallel_gemm
80 : USE preconditioner_apply, ONLY: apply_preconditioner_cfm_complex,&
81 : apply_preconditioner_complex_batch,&
82 : apply_preconditioner_dbcsr,&
83 : apply_preconditioner_dbcsr_complex,&
84 : apply_preconditioner_fm
85 : USE preconditioner_makes, ONLY: &
86 : make_complex_fermi_low_rank, make_complex_full_all, make_complex_full_all_covariant, &
87 : make_complex_full_kinetic, make_complex_full_s_inverse, make_complex_full_single, &
88 : make_complex_full_single_inverse, make_preconditioner_matrix, prepare_complex_full_all
89 : USE preconditioner_solvers, ONLY: solve_preconditioner,&
90 : transfer_dbcsr_to_fm,&
91 : transfer_fm_to_dbcsr
92 : USE preconditioner_types, ONLY: destroy_preconditioner,&
93 : init_preconditioner,&
94 : preconditioner_p_type,&
95 : preconditioner_type
96 : USE qs_environment_types, ONLY: get_qs_env,&
97 : qs_environment_type
98 : USE qs_mo_methods, ONLY: calculate_subspace_eigenvalues
99 : USE qs_mo_types, ONLY: deallocate_mo_set,&
100 : duplicate_mo_set,&
101 : get_mo_set,&
102 : mo_set_type,&
103 : set_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 :
117 : PUBLIC :: make_preconditioner, make_preconditioner_complex_fermi_low_rank, &
118 : make_preconditioner_complex_full_all, make_preconditioner_complex_full_all_batch, &
119 : make_preconditioner_complex_full_all_covariant, &
120 : make_preconditioner_complex_full_kinetic, &
121 : make_preconditioner_complex_full_s_inverse, &
122 : make_preconditioner_complex_full_single, &
123 : make_preconditioner_complex_full_single_inverse, restart_preconditioner, &
124 : dbcsr_pair_to_cfm
125 : PUBLIC :: apply_preconditioner, prepare_preconditioner, apply_preconditioner_complex_batch
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
134 : INTERFACE apply_preconditioner
135 : MODULE PROCEDURE apply_preconditioner_cfm_complex
136 : MODULE PROCEDURE apply_preconditioner_dbcsr
137 : MODULE PROCEDURE apply_preconditioner_dbcsr_complex
138 : MODULE PROCEDURE apply_preconditioner_fm
139 : END INTERFACE
140 :
141 : ! **************************************************************************************************
142 :
143 : CONTAINS
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 2 : 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 2 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: index_to_cell
164 : INTEGER, DIMENSION(3) :: delta
165 2 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
166 : LOGICAL :: used_fft
167 : REAL(KIND=dp) :: denominator, difference, numerator, value
168 2 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: dense_operator, xkp
169 2 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: blocks
170 : REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :), &
171 2 : 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 0 : CPASSERT(ASSOCIATED(operator_matrix))
177 2 : CPASSERT(preconditioner_env%lattice_fft_mode /= ot_lattice_fft_off)
178 2 : preconditioner_env%lattice_fft_active = .FALSE.
179 2 : preconditioner_env%lattice_projection_error = 0.0_dp
180 2 : preconditioner_env%lattice_corrected_projection_error = 0.0_dp
181 2 : preconditioner_env%lattice_cost_ratio = 0.0_dp
182 2 : preconditioner_env%lattice_storage_ratio = 0.0_dp
183 :
184 8 : IF (ANY(preconditioner_env%lattice_dims <= 0)) CPABORT("Lattice dimensions must be positive")
185 :
186 8 : IF (preconditioner_env%lattice_fft_mode == ot_lattice_fft_auto .AND. &
187 : PRODUCT(preconditioner_env%lattice_dims) < lattice_fft_auto_min_cells) THEN
188 0 : IF (preconditioner_env%para_env%is_source()) THEN
189 0 : logger => cp_get_default_logger()
190 0 : output_unit = cp_logger_get_default_io_unit(logger)
191 0 : IF (output_unit > 0) THEN
192 : WRITE (output_unit, '(T2,A,A,A,I0,A,I0)') &
193 0 : "OT| Lattice FFT AUTO kept the conventional ", TRIM(operator_name), &
194 0 : " solver: ", PRODUCT(preconditioner_env%lattice_dims), &
195 0 : " cells are fewer than ", lattice_fft_auto_min_cells
196 : END IF
197 : END IF
198 0 : RETURN
199 : END IF
200 :
201 2 : NULLIFY (fm_struct)
202 2 : CALL dbcsr_get_info(operator_matrix, nfullrows_total=n)
203 8 : ncell = PRODUCT(preconditioner_env%lattice_dims)
204 2 : IF (MODULO(n, ncell) /= 0) CPABORT("AO dimension does not match replicated cells")
205 2 : block_size = n/ncell
206 2 : 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) + &
210 : lattice_local_cost_ratio( &
211 : preconditioner_env%lattice_dims, &
212 2 : preconditioner_env%lattice_local_cells)
213 : preconditioner_env%lattice_storage_ratio = lattice_fft_storage_ratio( &
214 : preconditioner_env%lattice_dims, 1) + &
215 : lattice_local_storage_ratio( &
216 : preconditioner_env%lattice_dims, 1, &
217 2 : preconditioner_env%lattice_local_cells)
218 2 : 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 0 : IF (preconditioner_env%para_env%is_source()) THEN
222 0 : logger => cp_get_default_logger()
223 0 : output_unit = cp_logger_get_default_io_unit(logger)
224 0 : IF (output_unit > 0) THEN
225 : WRITE (output_unit, '(T2,A,A,A,F8.3,A,F8.3,A)') &
226 0 : "OT| Lattice FFT AUTO kept the conventional ", TRIM(operator_name), &
227 0 : " solver: estimated work ratio ", preconditioner_env%lattice_cost_ratio, &
228 0 : " storage ratio ", preconditioner_env%lattice_storage_ratio, " selected: F"
229 : END IF
230 : END IF
231 0 : RETURN
232 : END IF
233 : CALL cp_fm_struct_create(fm_struct, nrow_global=n, ncol_global=n, &
234 : context=preconditioner_env%ctxt, &
235 2 : para_env=preconditioner_env%para_env)
236 2 : CALL cp_fm_create(operator_fm, fm_struct, name="lattice preconditioner operator")
237 2 : CALL cp_fm_struct_release(fm_struct)
238 2 : CALL cp_fm_set_all(operator_fm, alpha=0.0_dp, beta=0.0_dp)
239 2 : CALL copy_dbcsr_to_fm(operator_matrix, operator_fm)
240 2 : CALL cp_fm_create(work, operator_fm%matrix_struct, name="lattice operator work")
241 2 : CALL cp_fm_uplo_to_full(operator_fm, work)
242 2 : CALL cp_fm_release(work)
243 : CALL cp_fm_get_info(operator_fm, nrow_local=nrow_local, ncol_local=ncol_local, &
244 2 : row_indices=row_indices, col_indices=col_indices, local_data=local_data)
245 10 : ALLOCATE (blocks(block_size, block_size, ncell), source=0.0_dp)
246 10 : ALLOCATE (index_to_cell(3, ncell), xkp(3, ncell))
247 2 : CALL lattice_grid(preconditioner_env%lattice_dims, index_to_cell, xkp)
248 2 : IF (preconditioner_env%lattice_local_cells > 0) THEN
249 0 : ALLOCATE (dense_operator(n, n), source=0.0_dp)
250 : END IF
251 2 : denominator = 0.0_dp
252 34 : DO j = 1, ncol_local
253 290 : DO i = 1, nrow_local
254 288 : IF (row_indices(i) <= col_indices(j)) THEN
255 136 : value = local_data(i, j)
256 136 : row_cell = (row_indices(i) - 1)/block_size + 1
257 136 : col_cell = (col_indices(j) - 1)/block_size + 1
258 136 : alpha = MODULO(row_indices(i) - 1, block_size) + 1
259 136 : beta = MODULO(col_indices(j) - 1, block_size) + 1
260 : delta = MODULO(index_to_cell(:, row_cell) - index_to_cell(:, col_cell), &
261 544 : preconditioner_env%lattice_dims)
262 : delta_index = delta(1) + preconditioner_env%lattice_dims(1)*(delta(2) + &
263 136 : preconditioner_env%lattice_dims(2)*delta(3)) + 1
264 : blocks(alpha, beta, delta_index) = blocks(alpha, beta, delta_index) + &
265 136 : value/REAL(ncell, KIND=dp)
266 136 : denominator = denominator + value*value
267 136 : IF (ALLOCATED(dense_operator)) THEN
268 0 : dense_operator(row_indices(i), col_indices(j)) = value
269 : END IF
270 136 : IF (row_indices(i) < col_indices(j)) THEN
271 : delta = MODULO(index_to_cell(:, col_cell) - index_to_cell(:, row_cell), &
272 480 : preconditioner_env%lattice_dims)
273 : delta_index = delta(1) + preconditioner_env%lattice_dims(1)*(delta(2) + &
274 120 : preconditioner_env%lattice_dims(2)*delta(3)) + 1
275 : blocks(beta, alpha, delta_index) = blocks(beta, alpha, delta_index) + &
276 120 : value/REAL(ncell, KIND=dp)
277 120 : denominator = denominator + value*value
278 : END IF
279 : END IF
280 : END DO
281 : END DO
282 2 : CALL preconditioner_env%para_env%sum(blocks)
283 2 : numerator = 0.0_dp
284 34 : DO j = 1, ncol_local
285 290 : DO i = 1, nrow_local
286 288 : IF (row_indices(i) <= col_indices(j)) THEN
287 136 : value = local_data(i, j)
288 136 : row_cell = (row_indices(i) - 1)/block_size + 1
289 136 : col_cell = (col_indices(j) - 1)/block_size + 1
290 136 : alpha = MODULO(row_indices(i) - 1, block_size) + 1
291 136 : beta = MODULO(col_indices(j) - 1, block_size) + 1
292 : delta = MODULO(index_to_cell(:, row_cell) - index_to_cell(:, col_cell), &
293 544 : preconditioner_env%lattice_dims)
294 : delta_index = delta(1) + preconditioner_env%lattice_dims(1)*(delta(2) + &
295 136 : preconditioner_env%lattice_dims(2)*delta(3)) + 1
296 136 : difference = value - blocks(alpha, beta, delta_index)
297 136 : numerator = numerator + difference*difference
298 136 : IF (row_indices(i) < col_indices(j)) THEN
299 : delta = MODULO(index_to_cell(:, col_cell) - index_to_cell(:, row_cell), &
300 480 : preconditioner_env%lattice_dims)
301 : delta_index = delta(1) + preconditioner_env%lattice_dims(1)*(delta(2) + &
302 120 : preconditioner_env%lattice_dims(2)*delta(3)) + 1
303 120 : difference = value - blocks(beta, alpha, delta_index)
304 120 : numerator = numerator + difference*difference
305 : END IF
306 : END IF
307 : END DO
308 : END DO
309 2 : CALL preconditioner_env%para_env%sum(numerator)
310 2 : CALL preconditioner_env%para_env%sum(denominator)
311 : preconditioner_env%lattice_projection_error = &
312 2 : SQRT(numerator/MAX(denominator, TINY(denominator)))
313 : CALL build_lattice_inverse_from_blocks( &
314 2 : blocks, preconditioner_env%lattice_dims, preconditioner_env%lattice_inverse_k, used_fft, info)
315 : preconditioner_env%lattice_corrected_projection_error = &
316 2 : preconditioner_env%lattice_projection_error
317 2 : IF (ALLOCATED(dense_operator)) THEN
318 0 : CALL preconditioner_env%para_env%sum(dense_operator)
319 0 : DO j = 1, n
320 0 : DO i = j + 1, n
321 0 : dense_operator(i, j) = dense_operator(j, i)
322 : END DO
323 : END DO
324 : CALL build_lattice_local_correction( &
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 0 : preconditioner_env%lattice_corrected_projection_error, local_info)
331 0 : IF (local_info /= 0) THEN
332 0 : CPABORT("LATTICE_FFT local correction found a non-positive coarse operator")
333 : END IF
334 0 : DEALLOCATE (dense_operator)
335 : END IF
336 2 : CALL cp_fm_release(operator_fm)
337 2 : DEALLOCATE (blocks, index_to_cell, xkp)
338 4 : SELECT CASE (preconditioner_env%lattice_fft_mode)
339 : CASE (ot_lattice_fft_on)
340 2 : IF (info /= 0) THEN
341 0 : CPABORT("LATTICE_FFT found a non-positive Fourier block for "//TRIM(operator_name))
342 : END IF
343 2 : IF (.NOT. used_fft) THEN
344 0 : CPABORT("LATTICE_FFT could not initialize its FFT for "//TRIM(operator_name))
345 : END IF
346 2 : preconditioner_env%lattice_fft_active = .TRUE.
347 : CASE (ot_lattice_fft_auto)
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 0 : preconditioner_env%lattice_projection_error, used_fft, info)
353 : CASE DEFAULT
354 2 : CPABORT("Unknown LATTICE_FFT mode")
355 : END SELECT
356 :
357 2 : IF (.NOT. preconditioner_env%lattice_fft_active .AND. &
358 : ALLOCATED(preconditioner_env%lattice_inverse_k)) THEN
359 0 : DEALLOCATE (preconditioner_env%lattice_inverse_k)
360 : END IF
361 2 : IF (preconditioner_env%para_env%is_source()) THEN
362 1 : logger => cp_get_default_logger()
363 1 : output_unit = cp_logger_get_default_io_unit(logger)
364 1 : IF (output_unit > 0) THEN
365 : WRITE (output_unit, '(T2,A,A,A,3(I0,1X),A,I0)') &
366 1 : "OT| Lattice FFT ", TRIM(operator_name), " cells: ", &
367 2 : preconditioner_env%lattice_dims, "AO block: ", preconditioner_env%lattice_block_size
368 : WRITE (output_unit, '(T2,A,ES12.4,A,L1)') &
369 1 : "OT| Block-circulant projection error: ", &
370 1 : preconditioner_env%lattice_projection_error, &
371 2 : " selected: ", preconditioner_env%lattice_fft_active
372 1 : IF (preconditioner_env%lattice_local_cells > 0) THEN
373 : WRITE (output_unit, '(T2,A,I0,A,ES12.4)') &
374 0 : "OT| Balanced local correction cells: ", &
375 0 : preconditioner_env%lattice_local_cells, " residual error: ", &
376 0 : preconditioner_env%lattice_corrected_projection_error
377 : END IF
378 : WRITE (output_unit, '(T2,A,F8.3,A,F8.3)') &
379 1 : "OT| Estimated lattice/dense work ratio: ", &
380 1 : preconditioner_env%lattice_cost_ratio, " storage ratio: ", &
381 2 : preconditioner_env%lattice_storage_ratio
382 : END IF
383 : END IF
384 :
385 12 : 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 0 : 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 0 : 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 0 : INTEGER, ALLOCATABLE, DIMENSION(:) :: local_indices
404 0 : 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 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: dense_operator, dense_sv, dense_vectors, &
409 0 : local_coarse_inverse, local_operator, &
410 0 : scaled_sv
411 : REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :), &
412 0 : POINTER :: local_data
413 : TYPE(cp_fm_type) :: sv_fm
414 : TYPE(cp_logger_type), POINTER :: logger
415 :
416 0 : CPASSERT(ASSOCIATED(matrix_s))
417 0 : CPASSERT(ASSOCIATED(preconditioner_env%fm))
418 0 : CPASSERT(ASSOCIATED(preconditioner_env%full_evals))
419 0 : CPASSERT(ASSOCIATED(preconditioner_env%occ_evals))
420 0 : CPASSERT(preconditioner_env%lattice_fft_mode /= ot_lattice_fft_off)
421 :
422 0 : preconditioner_env%lattice_fft_active = .FALSE.
423 0 : preconditioner_env%lattice_projection_error = 0.0_dp
424 0 : preconditioner_env%lattice_corrected_projection_error = 0.0_dp
425 0 : preconditioner_env%lattice_cost_ratio = 0.0_dp
426 0 : preconditioner_env%lattice_storage_ratio = 0.0_dp
427 0 : IF (ALLOCATED(preconditioner_env%lattice_state_inverse_k)) THEN
428 0 : DEALLOCATE (preconditioner_env%lattice_state_inverse_k)
429 : END IF
430 0 : IF (ALLOCATED(preconditioner_env%lattice_occ_vectors)) THEN
431 0 : DEALLOCATE (preconditioner_env%lattice_occ_vectors)
432 : END IF
433 0 : IF (ALLOCATED(preconditioner_env%lattice_occ_dual)) THEN
434 0 : DEALLOCATE (preconditioner_env%lattice_occ_dual)
435 : END IF
436 0 : IF (ALLOCATED(preconditioner_env%lattice_local_indices)) THEN
437 0 : DEALLOCATE (preconditioner_env%lattice_local_indices)
438 : END IF
439 0 : IF (ALLOCATED(preconditioner_env%lattice_local_operator)) THEN
440 0 : DEALLOCATE (preconditioner_env%lattice_local_operator)
441 : END IF
442 0 : IF (ALLOCATED(preconditioner_env%lattice_local_coarse_inverse)) THEN
443 0 : DEALLOCATE (preconditioner_env%lattice_local_coarse_inverse)
444 : END IF
445 0 : IF (ALLOCATED(preconditioner_env%lattice_state_local_indices)) THEN
446 0 : DEALLOCATE (preconditioner_env%lattice_state_local_indices)
447 : END IF
448 0 : IF (ALLOCATED(preconditioner_env%lattice_state_local_operator)) THEN
449 0 : DEALLOCATE (preconditioner_env%lattice_state_local_operator)
450 : END IF
451 0 : IF (ALLOCATED(preconditioner_env%lattice_state_local_coarse_inverse)) THEN
452 0 : DEALLOCATE (preconditioner_env%lattice_state_local_coarse_inverse)
453 : END IF
454 0 : IF (ANY(preconditioner_env%lattice_dims <= 0)) CPABORT("Lattice dimensions must be positive")
455 0 : IF (preconditioner_env%lattice_fft_mode == ot_lattice_fft_auto .AND. &
456 : PRODUCT(preconditioner_env%lattice_dims) < lattice_fft_auto_min_cells) THEN
457 0 : IF (preconditioner_env%para_env%is_source()) THEN
458 0 : logger => cp_get_default_logger()
459 0 : output_unit = cp_logger_get_default_io_unit(logger)
460 0 : IF (output_unit > 0) THEN
461 : WRITE (output_unit, '(T2,A,A,A,I0,A,I0)') &
462 0 : "OT| Lattice FFT AUTO kept the conventional ", TRIM(operator_name), &
463 0 : " solver: ", PRODUCT(preconditioner_env%lattice_dims), &
464 0 : " cells are fewer than ", lattice_fft_auto_min_cells
465 : END IF
466 : END IF
467 0 : RETURN
468 : END IF
469 :
470 0 : CALL cp_fm_get_info(preconditioner_env%fm, nrow_global=n, ncol_global=j)
471 0 : CPASSERT(j == n .AND. SIZE(preconditioner_env%full_evals) == n)
472 0 : nstate = SIZE(preconditioner_env%occ_evals)
473 0 : CPASSERT(nstate > 0)
474 0 : ncell = PRODUCT(preconditioner_env%lattice_dims)
475 0 : IF (MODULO(n, ncell) /= 0) CPABORT("AO dimension does not match replicated cells")
476 0 : 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 0 : preconditioner_env%para_env%num_pe)
481 : preconditioner_env%lattice_storage_ratio = lattice_fft_storage_ratio( &
482 0 : preconditioner_env%lattice_dims, nstate)
483 : preconditioner_env%lattice_cost_ratio = preconditioner_env%lattice_cost_ratio + &
484 : lattice_local_cost_ratio( &
485 : preconditioner_env%lattice_dims, &
486 0 : preconditioner_env%lattice_local_cells)
487 : preconditioner_env%lattice_storage_ratio = preconditioner_env%lattice_storage_ratio + &
488 : lattice_local_storage_ratio( &
489 : preconditioner_env%lattice_dims, nstate, &
490 0 : preconditioner_env%lattice_local_cells)
491 0 : 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 0 : IF (preconditioner_env%para_env%is_source()) THEN
495 0 : logger => cp_get_default_logger()
496 0 : output_unit = cp_logger_get_default_io_unit(logger)
497 0 : IF (output_unit > 0) THEN
498 : WRITE (output_unit, '(T2,A,A,A,F8.3,A,F8.3,A)') &
499 0 : "OT| Lattice FFT AUTO kept the conventional ", TRIM(operator_name), &
500 0 : " solver: estimated work ratio ", preconditioner_env%lattice_cost_ratio, &
501 0 : " storage ratio ", preconditioner_env%lattice_storage_ratio, " selected: F"
502 : END IF
503 : END IF
504 0 : RETURN
505 : END IF
506 0 : IF (preconditioner_env%lattice_local_cells > 0) THEN
507 : local_rank = preconditioner_env%lattice_local_cells* &
508 0 : preconditioner_env%lattice_block_size
509 0 : ALLOCATE (preconditioner_env%lattice_state_local_indices(local_rank, nstate), &
510 0 : preconditioner_env%lattice_state_local_operator(n, local_rank, nstate), &
511 0 : preconditioner_env%lattice_state_local_coarse_inverse(local_rank, local_rank, nstate))
512 : END IF
513 0 : CALL cp_fm_create(sv_fm, preconditioner_env%fm%matrix_struct, name="lattice FULL_ALL S V")
514 0 : 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 0 : row_indices=row_indices, col_indices=col_indices, local_data=local_data)
517 0 : ALLOCATE (dense_sv(n, n), dense_vectors(n, nstate), source=0.0_dp)
518 0 : DO j = 1, ncol_local
519 0 : DO i = 1, nrow_local
520 0 : dense_sv(row_indices(i), col_indices(j)) = local_data(i, j)
521 : END DO
522 : END DO
523 0 : CALL preconditioner_env%para_env%sum(dense_sv)
524 0 : 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 0 : row_indices=row_indices, col_indices=col_indices, local_data=local_data)
528 0 : DO j = 1, ncol_local
529 0 : IF (col_indices(j) > nstate) CYCLE
530 0 : DO i = 1, nrow_local
531 0 : dense_vectors(row_indices(i), col_indices(j)) = local_data(i, j)
532 : END DO
533 : END DO
534 0 : CALL preconditioner_env%para_env%sum(dense_vectors)
535 0 : ALLOCATE (preconditioner_env%lattice_occ_vectors(n, nstate), &
536 0 : preconditioner_env%lattice_occ_dual(n, nstate))
537 0 : preconditioner_env%lattice_occ_vectors(:, :) = dense_vectors
538 0 : preconditioner_env%lattice_occ_dual(:, :) = dense_sv(:, 1:nstate)
539 :
540 0 : ALLOCATE (dense_operator(n, n), scaled_sv(n, n))
541 0 : all_selected = .TRUE.
542 0 : DO istate = 1, nstate
543 0 : scaled_sv(:, :) = dense_sv
544 0 : DO j = 1, n
545 : scaled_sv(:, j) = scaled_sv(:, j)*MAX( &
546 : preconditioner_env%energy_gap, &
547 0 : preconditioner_env%full_evals(j) - preconditioner_env%occ_evals(istate))
548 : END DO
549 0 : dense_operator(:, :) = MATMUL(scaled_sv, TRANSPOSE(dense_sv))
550 : CALL build_lattice_inverse(dense_operator, preconditioner_env%lattice_dims, inverse_k, &
551 0 : projection_error, used_fft, info)
552 : preconditioner_env%lattice_projection_error = &
553 0 : MAX(preconditioner_env%lattice_projection_error, projection_error)
554 0 : corrected_projection_error = projection_error
555 0 : IF (preconditioner_env%lattice_local_cells > 0) THEN
556 : CALL build_lattice_local_correction( &
557 : dense_operator, preconditioner_env%lattice_dims, &
558 : preconditioner_env%lattice_local_cells, local_indices, &
559 0 : local_operator, local_coarse_inverse, corrected_projection_error, local_info)
560 0 : IF (local_info /= 0) THEN
561 0 : CPABORT("LATTICE_FFT local correction found a non-positive coarse operator")
562 : END IF
563 0 : preconditioner_env%lattice_state_local_indices(:, istate) = local_indices
564 0 : preconditioner_env%lattice_state_local_operator(:, :, istate) = local_operator
565 0 : preconditioner_env%lattice_state_local_coarse_inverse(:, :, istate) = local_coarse_inverse
566 0 : DEALLOCATE (local_coarse_inverse, local_indices, local_operator)
567 : END IF
568 : preconditioner_env%lattice_corrected_projection_error = &
569 0 : 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 0 : preconditioner_env%para_env%num_pe, projection_error, used_fft, info)
573 0 : IF (preconditioner_env%lattice_fft_mode == ot_lattice_fft_on) THEN
574 0 : IF (info /= 0) THEN
575 0 : CPABORT("LATTICE_FFT found a non-positive Fourier block for "//TRIM(operator_name))
576 : END IF
577 0 : IF (.NOT. used_fft) THEN
578 0 : CPABORT("LATTICE_FFT could not initialize its FFT for "//TRIM(operator_name))
579 : END IF
580 : ELSE
581 0 : all_selected = all_selected .AND. selected
582 : END IF
583 0 : IF (.NOT. ALLOCATED(preconditioner_env%lattice_state_inverse_k)) THEN
584 0 : ALLOCATE (preconditioner_env%lattice_state_inverse_k( &
585 0 : SIZE(inverse_k, 1), SIZE(inverse_k, 2), SIZE(inverse_k, 3), nstate))
586 : END IF
587 0 : IF (info == 0 .AND. used_fft) THEN
588 0 : preconditioner_env%lattice_state_inverse_k(:, :, :, istate) = inverse_k
589 : END IF
590 0 : DEALLOCATE (inverse_k)
591 : END DO
592 :
593 : preconditioner_env%lattice_fft_active = &
594 0 : preconditioner_env%lattice_fft_mode == ot_lattice_fft_on .OR. all_selected
595 0 : IF (.NOT. preconditioner_env%lattice_fft_active .AND. &
596 : ALLOCATED(preconditioner_env%lattice_state_inverse_k)) THEN
597 0 : DEALLOCATE (preconditioner_env%lattice_state_inverse_k)
598 : END IF
599 0 : IF (.NOT. preconditioner_env%lattice_fft_active) THEN
600 0 : DEALLOCATE (preconditioner_env%lattice_occ_vectors, &
601 0 : preconditioner_env%lattice_occ_dual)
602 0 : IF (ALLOCATED(preconditioner_env%lattice_state_local_indices)) THEN
603 : DEALLOCATE (preconditioner_env%lattice_state_local_indices, &
604 0 : preconditioner_env%lattice_state_local_operator, &
605 0 : preconditioner_env%lattice_state_local_coarse_inverse)
606 : END IF
607 : END IF
608 0 : IF (preconditioner_env%para_env%is_source()) THEN
609 0 : logger => cp_get_default_logger()
610 0 : output_unit = cp_logger_get_default_io_unit(logger)
611 0 : IF (output_unit > 0) THEN
612 : WRITE (output_unit, '(T2,A,A,A,3(I0,1X),A,I0,A,I0)') &
613 0 : "OT| Lattice FFT ", TRIM(operator_name), " cells: ", &
614 0 : preconditioner_env%lattice_dims, "AO block: ", &
615 0 : preconditioner_env%lattice_block_size, " states: ", nstate
616 : WRITE (output_unit, '(T2,A,ES12.4,A,L1)') &
617 0 : "OT| Maximum block-circulant projection error: ", &
618 0 : preconditioner_env%lattice_projection_error, &
619 0 : " selected: ", preconditioner_env%lattice_fft_active
620 0 : IF (preconditioner_env%lattice_local_cells > 0) THEN
621 : WRITE (output_unit, '(T2,A,I0,A,ES12.4)') &
622 0 : "OT| Balanced local correction cells/state: ", &
623 0 : preconditioner_env%lattice_local_cells, " maximum residual error: ", &
624 0 : preconditioner_env%lattice_corrected_projection_error
625 : END IF
626 : WRITE (output_unit, '(T2,A,F8.3,A,F8.3)') &
627 0 : "OT| Estimated lattice/dense work ratio: ", &
628 0 : preconditioner_env%lattice_cost_ratio, " storage ratio: ", &
629 0 : preconditioner_env%lattice_storage_ratio
630 : END IF
631 : END IF
632 :
633 0 : DEALLOCATE (dense_operator, dense_sv, dense_vectors, scaled_sv)
634 :
635 0 : 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 0 : 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 0 : CPASSERT(ASSOCIATED(preconditioner_env%fm))
651 0 : CPASSERT(.NOT. ASSOCIATED(preconditioner_env%base_fm))
652 :
653 0 : IF (preconditioner_env%lattice_fft_mode /= ot_lattice_fft_off) THEN
654 0 : IF (.NOT. PRESENT(matrix_s)) THEN
655 0 : CPABORT("LATTICE_FFT requires a non-unit overlap metric for FERMI_LOW_RANK")
656 : END IF
657 0 : CPASSERT(ASSOCIATED(matrix_s))
658 : CALL prepare_lattice_operator(preconditioner_env, matrix_s, "FERMI_LOW_RANK overlap", &
659 0 : SIZE(preconditioner_env%occ_evals))
660 0 : IF (preconditioner_env%lattice_fft_active) THEN
661 0 : preconditioner_env%low_rank_base = ot_low_rank_base_lattice_fft
662 0 : RETURN
663 : END IF
664 0 : preconditioner_env%low_rank_base = ot_low_rank_base_overlap
665 : END IF
666 :
667 0 : IF (preconditioner_env%low_rank_base /= ot_low_rank_base_overlap) THEN
668 0 : CPABORT("Unknown FERMI_LOW_RANK base operator")
669 : END IF
670 0 : ALLOCATE (preconditioner_env%base_fm)
671 : CALL cp_fm_create(preconditioner_env%base_fm, preconditioner_env%fm%matrix_struct, &
672 0 : name="FERMI_LOW_RANK overlap inverse")
673 :
674 0 : IF (PRESENT(matrix_s)) THEN
675 0 : CPASSERT(ASSOCIATED(matrix_s))
676 0 : CALL copy_dbcsr_to_fm(matrix_s, preconditioner_env%base_fm)
677 0 : CALL cp_fm_cholesky_decompose(preconditioner_env%base_fm, info_out=info)
678 0 : IF (info /= 0) CPABORT("FERMI_LOW_RANK requires a positive-definite overlap matrix")
679 0 : CALL cp_fm_cholesky_invert(preconditioner_env%base_fm)
680 0 : CALL cp_fm_create(work, preconditioner_env%base_fm%matrix_struct, name="low-rank base work")
681 0 : CALL cp_fm_uplo_to_full(preconditioner_env%base_fm, work)
682 0 : CALL cp_fm_release(work)
683 : ELSE
684 0 : CALL cp_fm_set_all(preconditioner_env%base_fm, alpha=0.0_dp, beta=1.0_dp)
685 : END IF
686 :
687 0 : 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 9936 : 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 9936 : 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 9936 : CALL timeset(routineN, handle)
742 :
743 9936 : CALL get_mo_set(mo_set=mo_set, mo_coeff=mo_coeff, mo_coeff_b=mo_coeff_b, homo=nhomo)
744 9936 : use_mo_coeff_b = mo_set%use_mo_coeff_b
745 9936 : 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 9936 : my_convert_precond_to_dbcsr = .FALSE.
749 9936 : IF (PRESENT(convert_precond_to_dbcsr)) my_convert_precond_to_dbcsr = convert_precond_to_dbcsr
750 9936 : 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 0 : 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 9936 : IF (preconditioner_env%solver == ot_precond_solver_update) THEN
761 4 : CALL transfer_fm_to_dbcsr(preconditioner_env%fm, preconditioner_env%dbcsr_matrix, matrix_h)
762 : END IF
763 :
764 9936 : needs_full_spectrum = .FALSE.
765 9936 : needs_homo = .FALSE.
766 :
767 3792 : SELECT CASE (precon_type)
768 : CASE (ot_precond_full_all, ot_precond_fermi_low_rank, ot_precond_full_all_covariant)
769 3792 : needs_full_spectrum = .TRUE.
770 : ! both of them need the coefficients as fm's, more matrix mess
771 3792 : IF (use_mo_coeff_b) THEN
772 3488 : CALL copy_dbcsr_to_fm(mo_coeff_b, mo_coeff)
773 : END IF
774 : CASE (ot_precond_full_single)
775 38 : needs_homo = .TRUE.
776 : ! XXXX to be removed if homo estimate only is implemented
777 38 : needs_full_spectrum = .TRUE.
778 : CASE (ot_precond_full_kinetic, ot_precond_s_inverse, ot_precond_full_single_inverse)
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 9936 : CPABORT("The preconditioner is unknown ...")
783 : END SELECT
784 :
785 29520 : ALLOCATE (eigenvalues_ot(k))
786 9936 : energy_homo = 0.0_dp
787 9936 : 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 3830 : 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 3518 : blacs_env=mo_coeff%matrix_struct%context)
799 : ELSE
800 : CALL calculate_subspace_eigenvalues(mo_coeff, matrix_h, &
801 312 : eigenvalues_ot, do_rotation=.FALSE.)
802 : END IF
803 3830 : IF (k > 0) THEN
804 3720 : CPASSERT(nhomo > 0 .AND. nhomo <= k)
805 3720 : energy_homo = eigenvalues_ot(nhomo)
806 : END IF
807 : ELSE
808 6106 : IF (needs_homo) THEN
809 0 : 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 9936 : my_solver_type = solver_type
817 9936 : preconditioner_env%in_use = precon_type
818 9936 : preconditioner_env%cholesky_use = cholesky_reduce
819 9936 : IF (PRESENT(chol_type)) preconditioner_env%cholesky_use = chol_type
820 : preconditioner_env%in_use = precon_type
821 9936 : preconditioner_env%lattice_fft_active = .FALSE.
822 9936 : preconditioner_env%lattice_corrected_projection_error = 0.0_dp
823 9936 : IF (ALLOCATED(preconditioner_env%lattice_inverse_k)) THEN
824 0 : DEALLOCATE (preconditioner_env%lattice_inverse_k)
825 : END IF
826 9936 : IF (ALLOCATED(preconditioner_env%lattice_state_inverse_k)) THEN
827 0 : DEALLOCATE (preconditioner_env%lattice_state_inverse_k)
828 : END IF
829 9936 : IF (ALLOCATED(preconditioner_env%lattice_occ_vectors)) THEN
830 0 : DEALLOCATE (preconditioner_env%lattice_occ_vectors)
831 : END IF
832 9936 : IF (ALLOCATED(preconditioner_env%lattice_occ_dual)) THEN
833 0 : DEALLOCATE (preconditioner_env%lattice_occ_dual)
834 : END IF
835 9936 : IF (ALLOCATED(preconditioner_env%lattice_local_indices)) THEN
836 0 : DEALLOCATE (preconditioner_env%lattice_local_indices)
837 : END IF
838 9936 : IF (ALLOCATED(preconditioner_env%lattice_local_operator)) THEN
839 0 : DEALLOCATE (preconditioner_env%lattice_local_operator)
840 : END IF
841 9936 : IF (ALLOCATED(preconditioner_env%lattice_local_coarse_inverse)) THEN
842 0 : DEALLOCATE (preconditioner_env%lattice_local_coarse_inverse)
843 : END IF
844 9936 : IF (ALLOCATED(preconditioner_env%lattice_state_local_indices)) THEN
845 0 : DEALLOCATE (preconditioner_env%lattice_state_local_indices)
846 : END IF
847 9936 : IF (ALLOCATED(preconditioner_env%lattice_state_local_operator)) THEN
848 0 : DEALLOCATE (preconditioner_env%lattice_state_local_operator)
849 : END IF
850 9936 : IF (ALLOCATED(preconditioner_env%lattice_state_local_coarse_inverse)) THEN
851 0 : DEALLOCATE (preconditioner_env%lattice_state_local_coarse_inverse)
852 : END IF
853 9936 : IF (nhomo == k) THEN
854 : CALL make_preconditioner_matrix(preconditioner_env, matrix_h, matrix_s, matrix_t, mo_coeff, &
855 9848 : 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 88 : para_env=preconditioner_env%para_env)
860 88 : CALL cp_fm_create(mo_occ, fm_struct)
861 88 : CALL cp_fm_to_fm(mo_coeff, mo_occ, nhomo)
862 88 : CALL cp_fm_struct_release(fm_struct)
863 : !
864 : CALL make_preconditioner_matrix(preconditioner_env, matrix_h, matrix_s, matrix_t, mo_occ, &
865 88 : energy_homo, eigenvalues_ot(1:nhomo), energy_gap, my_solver_type)
866 : !
867 88 : CALL cp_fm_release(mo_occ)
868 : END IF
869 :
870 9936 : IF (preconditioner_env%lattice_fft_mode /= ot_lattice_fft_off) THEN
871 0 : SELECT CASE (precon_type)
872 : CASE (ot_precond_full_kinetic)
873 0 : CPASSERT(ASSOCIATED(preconditioner_env%sparse_matrix))
874 : CALL prepare_lattice_operator(preconditioner_env, preconditioner_env%sparse_matrix, &
875 0 : "FULL_KINETIC", nhomo)
876 : CASE (ot_precond_full_single_inverse)
877 0 : CPASSERT(ASSOCIATED(preconditioner_env%sparse_matrix))
878 : CALL prepare_lattice_operator(preconditioner_env, preconditioner_env%sparse_matrix, &
879 0 : "FULL_SINGLE_INVERSE", nhomo)
880 : CASE (ot_precond_s_inverse)
881 2 : CPASSERT(ASSOCIATED(preconditioner_env%sparse_matrix))
882 : CALL prepare_lattice_operator(preconditioner_env, preconditioner_env%sparse_matrix, &
883 2 : "FULL_S_INVERSE", nhomo)
884 : CASE (ot_precond_full_all_covariant)
885 0 : CPASSERT(PRESENT(matrix_s))
886 0 : CPASSERT(ASSOCIATED(matrix_s))
887 : CALL prepare_lattice_state_operator(preconditioner_env, matrix_s, &
888 2 : "FULL_ALL_COVARIANT")
889 : END SELECT
890 : END IF
891 :
892 9936 : IF (.NOT. preconditioner_env%lattice_fft_active) THEN
893 9934 : CALL solve_preconditioner(my_solver_type, preconditioner_env, matrix_s, matrix_h)
894 : ELSE
895 2 : preconditioner_env%solver = ot_precond_solver_inv_chol
896 : END IF
897 :
898 9936 : IF (precon_type == ot_precond_fermi_low_rank) THEN
899 0 : preconditioner_env%spectral_reference = energy_homo
900 : preconditioner_env%spectral_window = &
901 0 : MAX(fermi_low_rank_min_window, fermi_low_rank_gap_factor*preconditioner_env%energy_gap)
902 0 : preconditioner_env%spectral_base_scale = 1.0_dp/preconditioner_env%spectral_window
903 0 : max_rank = preconditioner_env%spectral_max_rank
904 0 : 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 0 : fermi_low_rank_degeneracy_tolerance)
908 0 : CALL prepare_low_rank_base(preconditioner_env, matrix_s)
909 0 : 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 0 : DO WHILE (preconditioner_env%spectral_rank > 0)
913 0 : IF (low_rank_inverse_weight( &
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 0 : preconditioner_env%spectral_rank = preconditioner_env%spectral_rank - 1
918 : END DO
919 : END IF
920 0 : IF (preconditioner_env%para_env%is_source()) THEN
921 0 : logger => cp_get_default_logger()
922 0 : output_unit = cp_logger_get_default_io_unit(logger)
923 0 : IF (output_unit > 0) THEN
924 : WRITE (output_unit, '(T2,A,F14.8,A,I0)') &
925 0 : "OT| FERMI_LOW_RANK common reference [hartree]: ", &
926 0 : preconditioner_env%spectral_reference, " retained rank: ", &
927 0 : 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 9936 : IF (.NOT. preconditioner_env%lattice_fft_active) THEN
936 9934 : IF (my_convert_precond_to_dbcsr) THEN
937 8210 : CALL transfer_fm_to_dbcsr(preconditioner_env%fm, preconditioner_env%dbcsr_matrix, matrix_h)
938 8210 : 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 0 : 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 1724 : preconditioner_env%para_env, preconditioner_env%ctxt)
946 : END IF
947 : END IF
948 :
949 9936 : DEALLOCATE (eigenvalues_ot)
950 :
951 9936 : CALL timestop(handle)
952 :
953 9936 : 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 304 : 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 76 : CALL timeset(routineN, handle)
985 :
986 76 : CPASSERT(ASSOCIATED(matrix_c_re))
987 76 : CPASSERT(ASSOCIATED(matrix_c_im))
988 76 : CPASSERT(ASSOCIATED(matrix_h_re))
989 76 : CPASSERT(ASSOCIATED(matrix_h_im))
990 76 : CPASSERT(ASSOCIATED(matrix_s_re))
991 76 : CPASSERT(ASSOCIATED(matrix_s_im))
992 76 : IF (solver_type /= ot_precond_solver_default) THEN
993 0 : 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 76 : 'complex FULL_ALL C', matrix_c)
997 : CALL dbcsr_pair_to_cfm(matrix_h_re, matrix_h_im, preconditioner_env, &
998 76 : 'complex FULL_ALL H', matrix_h)
999 : CALL dbcsr_pair_to_cfm(matrix_s_re, matrix_s_im, preconditioner_env, &
1000 76 : 'complex FULL_ALL S', matrix_s)
1001 : CALL make_complex_full_all(preconditioner_env, matrix_c, matrix_h, matrix_s, &
1002 76 : energy_gap=energy_gap)
1003 :
1004 76 : CALL cp_cfm_release(matrix_c)
1005 76 : CALL cp_cfm_release(matrix_s)
1006 76 : CALL cp_cfm_release(matrix_h)
1007 :
1008 76 : CALL timestop(handle)
1009 :
1010 76 : END SUBROUTINE make_preconditioner_complex_full_all
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 50 : 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 50 : TYPE(full_all_work_type), ALLOCATABLE :: work(:)
1037 : TYPE(preconditioner_type), POINTER :: p
1038 :
1039 50 : handled = .FALSE.
1040 50 : nchannels = SIZE(preconditioners)
1041 50 : IF (nchannels == 0 .OR. solver_type /= ot_precond_solver_default) RETURN
1042 50 : IF (diag_lib_explicit .OR. diag_check_requested() .OR. direct_generalized_diagonalization) RETURN
1043 109 : DO i = 1, nchannels
1044 109 : IF (preconditioners(i)%preconditioner%para_env%num_pe /= 1) RETURN
1045 : END DO
1046 30 : CALL timeset(routineN, handle)
1047 30 : nworkers = 1
1048 30 : !$ nworkers = MIN(nchannels, omp_get_max_threads())
1049 120 : ALLOCATE (work(nworkers))
1050 89 : DO first = 1, nchannels, nworkers
1051 59 : nbatch = MIN(nworkers, nchannels - first + 1)
1052 118 : DO j = 1, nbatch
1053 59 : i = first + j - 1
1054 59 : p => preconditioners(i)%preconditioner
1055 59 : CALL dbcsr_pair_to_cfm(matrix_c(1, i)%matrix, matrix_c(2, i)%matrix, p, 'FULL_ALL C', work(j)%c)
1056 59 : CALL cp_cfm_create(work(j)%h, matrix_h(1, i)%matrix%matrix_struct, name='FULL_ALL H')
1057 59 : CALL cp_cfm_create(work(j)%s, matrix_s(1, i)%matrix%matrix_struct, name='FULL_ALL S')
1058 59 : CALL cp_fm_to_cfm(matrix_h(1, i)%matrix, matrix_h(2, i)%matrix, work(j)%h)
1059 59 : 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 59 : energy_gap=energy_gap)
1062 59 : n = SIZE(p%full_evals)
1063 : ! Match the existing Cholesky reduction and back transformation exactly.
1064 59 : CALL cp_cfm_cholesky_decompose(work(j)%b)
1065 59 : CALL cp_cfm_triangular_invert(work(j)%b)
1066 59 : CALL cp_cfm_triangular_multiply(work(j)%b, work(j)%a, side='R')
1067 59 : CALL cp_cfm_triangular_multiply(work(j)%b, work(j)%a, transa_tr='C')
1068 59 : CALL cp_cfm_local_workspace_prepare(work(j)%solver, n, work(j)%a%local_data)
1069 59 : CALL cp_cfm_release(work(j)%s)
1070 118 : CALL cp_cfm_release(work(j)%h)
1071 : END DO
1072 59 : !$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 148 : DO j = 1, nbatch
1078 59 : IF (work(j)%info /= 0) CPABORT('Local FULL_ALL generalized diagonalization failed')
1079 59 : i = first + j - 1
1080 59 : p => preconditioners(i)%preconditioner
1081 2536 : p%full_evals(:) = work(j)%solver%evals(:)
1082 59 : CALL cp_cfm_triangular_multiply(work(j)%b, work(j)%a)
1083 59 : CALL cp_cfm_to_cfm(work(j)%a, p%complex_fm)
1084 59 : k = SIZE(p%occ_evals)
1085 643 : p%full_evals(1:k) = p%occ_evals(:)
1086 59 : CALL cp_cfm_to_cfm(work(j)%c, p%complex_fm, k)
1087 59 : CALL cp_cfm_release(work(j)%b)
1088 59 : CALL cp_cfm_release(work(j)%a)
1089 118 : CALL cp_cfm_release(work(j)%c)
1090 : END DO
1091 : END DO
1092 30 : handled = .TRUE.
1093 30 : CALL timestop(handle)
1094 110 : END SUBROUTINE make_preconditioner_complex_full_all_batch
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 32 : 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 8 : REAL(KIND=dp), DIMENSION(:), POINTER :: eigenvalues
1121 : TYPE(cp_cfm_type) :: matrix_h, matrix_s
1122 :
1123 8 : NULLIFY (eigenvalues)
1124 0 : CPASSERT(ASSOCIATED(matrix_h_re))
1125 8 : CPASSERT(ASSOCIATED(matrix_h_im))
1126 8 : CPASSERT(ASSOCIATED(matrix_s_re))
1127 8 : CPASSERT(ASSOCIATED(matrix_s_im))
1128 8 : IF (solver_type /= ot_precond_solver_default) THEN
1129 0 : CPABORT('Complex FULL_SINGLE supports only PRECOND_SOLVER DEFAULT')
1130 : END IF
1131 8 : CALL get_mo_set(mo_set, homo=homo, eigenvalues=eigenvalues)
1132 8 : CPASSERT(ASSOCIATED(eigenvalues))
1133 8 : CPASSERT(homo > 0 .AND. homo <= SIZE(eigenvalues))
1134 :
1135 : CALL dbcsr_pair_to_cfm(matrix_h_re, matrix_h_im, preconditioner_env, &
1136 8 : 'complex FULL_SINGLE H', matrix_h)
1137 : CALL dbcsr_pair_to_cfm(matrix_s_re, matrix_s_im, preconditioner_env, &
1138 8 : 'complex FULL_SINGLE S', matrix_s)
1139 : CALL make_complex_full_single(preconditioner_env, matrix_h, matrix_s, &
1140 8 : eigenvalues(homo), energy_gap)
1141 8 : CALL cp_cfm_release(matrix_s)
1142 8 : CALL cp_cfm_release(matrix_h)
1143 :
1144 8 : END SUBROUTINE make_preconditioner_complex_full_single
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 : ! **************************************************************************************************
1159 184 : SUBROUTINE make_preconditioner_complex_fermi_low_rank( &
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 46 : IF (solver_type /= ot_precond_solver_default) THEN
1172 0 : 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 46 : 'complex FERMI_LOW_RANK C', matrix_c)
1177 : CALL dbcsr_pair_to_cfm(matrix_h_re, matrix_h_im, preconditioner_env, &
1178 46 : 'complex FERMI_LOW_RANK H', matrix_h)
1179 : CALL dbcsr_pair_to_cfm(matrix_s_re, matrix_s_im, preconditioner_env, &
1180 46 : 'complex FERMI_LOW_RANK S', matrix_s)
1181 : CALL make_complex_fermi_low_rank( &
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 46 : fermi_low_rank_degeneracy_tolerance)
1185 :
1186 46 : CALL cp_cfm_release(matrix_s)
1187 46 : CALL cp_cfm_release(matrix_h)
1188 46 : CALL cp_cfm_release(matrix_c)
1189 :
1190 46 : END SUBROUTINE make_preconditioner_complex_fermi_low_rank
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 : ! **************************************************************************************************
1205 200 : SUBROUTINE make_preconditioner_complex_full_all_covariant( &
1206 : preconditioner_env, matrix_c_re, matrix_c_im, matrix_h_re, matrix_h_im, &
1207 50 : 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 50 : IF (solver_type /= ot_precond_solver_default) THEN
1220 0 : 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 50 : 'complex FULL_ALL_COVARIANT C', matrix_c)
1224 : CALL dbcsr_pair_to_cfm(matrix_h_re, matrix_h_im, preconditioner_env, &
1225 50 : 'complex FULL_ALL_COVARIANT H', matrix_h)
1226 : CALL dbcsr_pair_to_cfm(matrix_s_re, matrix_s_im, preconditioner_env, &
1227 50 : 'complex FULL_ALL_COVARIANT S', matrix_s)
1228 : CALL make_complex_full_all_covariant(preconditioner_env, matrix_c, matrix_h, matrix_s, &
1229 92 : energy_gap, occupation_signature)
1230 :
1231 50 : CALL cp_cfm_release(matrix_s)
1232 50 : CALL cp_cfm_release(matrix_h)
1233 50 : CALL cp_cfm_release(matrix_c)
1234 :
1235 50 : END SUBROUTINE make_preconditioner_complex_full_all_covariant
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 : ! **************************************************************************************************
1249 484 : SUBROUTINE make_preconditioner_complex_full_single_inverse(preconditioner_env, &
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 44 : NULLIFY (fm_struct_mo, fm_struct_square)
1269 0 : CPASSERT(ASSOCIATED(matrix_c_re))
1270 44 : CPASSERT(ASSOCIATED(matrix_c_im))
1271 44 : CPASSERT(ASSOCIATED(matrix_h_re))
1272 44 : CPASSERT(ASSOCIATED(matrix_h_im))
1273 44 : CPASSERT(ASSOCIATED(matrix_s_re))
1274 44 : CPASSERT(ASSOCIATED(matrix_s_im))
1275 44 : IF (solver_type /= ot_precond_solver_default .AND. &
1276 : solver_type /= ot_precond_solver_inv_chol) THEN
1277 0 : CPABORT("Complex FULL_SINGLE_INVERSE supports only PRECOND_SOLVER DEFAULT/INVERSE_CHOLESKY")
1278 : END IF
1279 44 : CALL dbcsr_get_info(matrix_c_re, nfullrows_total=n, nfullcols_total=k)
1280 44 : 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 44 : 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 44 : para_env=preconditioner_env%para_env)
1288 44 : CALL cp_fm_create(matrix_c_re_fm, fm_struct_mo, name='complex preconditioner C real')
1289 44 : CALL cp_fm_create(matrix_c_im_fm, fm_struct_mo, name='complex preconditioner C imaginary')
1290 44 : CALL cp_fm_create(matrix_h_re_fm, fm_struct_square, name='complex preconditioner H real')
1291 44 : CALL cp_fm_create(matrix_h_im_fm, fm_struct_square, name='complex preconditioner H imaginary')
1292 44 : CALL cp_fm_create(matrix_s_re_fm, fm_struct_square, name='complex preconditioner S real')
1293 44 : CALL cp_fm_create(matrix_s_im_fm, fm_struct_square, name='complex preconditioner S imaginary')
1294 44 : CALL cp_fm_struct_release(fm_struct_mo)
1295 44 : CALL cp_fm_struct_release(fm_struct_square)
1296 44 : CALL copy_dbcsr_to_fm(matrix_c_re, matrix_c_re_fm)
1297 44 : CALL copy_dbcsr_to_fm(matrix_c_im, matrix_c_im_fm)
1298 44 : CALL copy_dbcsr_to_fm(matrix_h_re, matrix_h_re_fm)
1299 44 : CALL copy_dbcsr_to_fm(matrix_h_im, matrix_h_im_fm)
1300 44 : CALL copy_dbcsr_to_fm(matrix_s_re, matrix_s_re_fm)
1301 44 : CALL copy_dbcsr_to_fm(matrix_s_im, matrix_s_im_fm)
1302 :
1303 44 : CALL cp_cfm_create(matrix_c, matrix_c_re_fm%matrix_struct, name='complex preconditioner C')
1304 44 : CALL cp_cfm_create(matrix_h, matrix_h_re_fm%matrix_struct, name='complex preconditioner H')
1305 44 : CALL cp_cfm_create(matrix_s, matrix_s_re_fm%matrix_struct, name='complex preconditioner S')
1306 44 : CALL cp_fm_to_cfm(matrix_c_re_fm, matrix_c_im_fm, matrix_c)
1307 44 : CALL cp_fm_to_cfm(matrix_h_re_fm, matrix_h_im_fm, matrix_h)
1308 44 : 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 44 : energy_gap)
1311 44 : preconditioner_env%solver = solver_type
1312 :
1313 44 : CALL cp_cfm_release(matrix_s)
1314 44 : CALL cp_cfm_release(matrix_h)
1315 44 : CALL cp_cfm_release(matrix_c)
1316 44 : CALL cp_fm_release(matrix_s_im_fm)
1317 44 : CALL cp_fm_release(matrix_s_re_fm)
1318 44 : CALL cp_fm_release(matrix_h_im_fm)
1319 44 : CALL cp_fm_release(matrix_h_re_fm)
1320 44 : CALL cp_fm_release(matrix_c_im_fm)
1321 44 : CALL cp_fm_release(matrix_c_re_fm)
1322 :
1323 44 : END SUBROUTINE make_preconditioner_complex_full_single_inverse
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 164 : 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 82 : IF (solver_type /= ot_precond_solver_default .AND. &
1342 : solver_type /= ot_precond_solver_inv_chol) THEN
1343 0 : 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 82 : 'complex preconditioner S', matrix_s)
1347 82 : CALL make_complex_full_s_inverse(preconditioner_env, matrix_s)
1348 82 : preconditioner_env%solver = solver_type
1349 82 : CALL cp_cfm_release(matrix_s)
1350 :
1351 82 : END SUBROUTINE make_preconditioner_complex_full_s_inverse
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 1188 : 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 396 : IF (solver_type /= ot_precond_solver_default .AND. &
1377 : solver_type /= ot_precond_solver_inv_chol) THEN
1378 0 : 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 396 : 'complex preconditioner T', matrix_t)
1382 : CALL dbcsr_pair_to_cfm(matrix_s_re, matrix_s_im, preconditioner_env, &
1383 396 : 'complex preconditioner S', matrix_s)
1384 396 : CALL make_complex_full_kinetic(preconditioner_env, matrix_t, matrix_s, energy_gap)
1385 396 : preconditioner_env%solver = solver_type
1386 396 : CALL cp_cfm_release(matrix_s)
1387 396 : CALL cp_cfm_release(matrix_t)
1388 :
1389 396 : END SUBROUTINE make_preconditioner_complex_full_kinetic
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 1569 : 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 1569 : 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 1569 : NULLIFY (fm_struct)
1412 0 : CPASSERT(ASSOCIATED(matrix_re))
1413 1569 : CPASSERT(ASSOCIATED(matrix_im))
1414 1569 : 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 1569 : para_env=preconditioner_env%para_env)
1418 1569 : CALL cp_cfm_create(matrix, fm_struct, name=matrix_name)
1419 : IF (preconditioner_env%para_env%num_pe == 1 .AND. &
1420 1569 : dbcsr_get_matrix_type(matrix_re) == dbcsr_type_no_symmetry .AND. &
1421 : dbcsr_get_matrix_type(matrix_im) == dbcsr_type_no_symmetry) THEN
1422 3556 : ALLOCATE (buffer(nrow, ncol))
1423 889 : CALL copy_dbcsr_to_local(matrix_re, buffer)
1424 434769 : matrix%local_data(1:nrow, 1:ncol) = CMPLX(buffer, 0.0_dp, KIND=dp)
1425 889 : CALL copy_dbcsr_to_local(matrix_im, buffer)
1426 : matrix%local_data(1:nrow, 1:ncol) = &
1427 434769 : matrix%local_data(1:nrow, 1:ncol) + CMPLX(0.0_dp, buffer, KIND=dp)
1428 : ELSE
1429 680 : CALL cp_fm_create(matrix_re_fm, fm_struct, name=TRIM(matrix_name)//' real')
1430 680 : CALL cp_fm_create(matrix_im_fm, fm_struct, name=TRIM(matrix_name)//' imaginary')
1431 680 : CALL copy_dbcsr_to_fm(matrix_re, matrix_re_fm)
1432 680 : CALL copy_dbcsr_to_fm(matrix_im, matrix_im_fm)
1433 680 : CALL cp_fm_to_cfm(matrix_re_fm, matrix_im_fm, matrix)
1434 680 : CALL cp_fm_release(matrix_im_fm)
1435 2249 : CALL cp_fm_release(matrix_re_fm)
1436 : END IF
1437 1569 : CALL cp_fm_struct_release(fm_struct)
1438 :
1439 3138 : 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 7565 : 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 7565 : NULLIFY (para_env, blacs_env)
1460 7565 : CALL get_qs_env(qs_env, para_env=para_env, blacs_env=blacs_env)
1461 :
1462 7565 : IF (ASSOCIATED(preconditioner)) THEN
1463 6700 : SELECT CASE (prec_type)
1464 : CASE (ot_precond_full_all, ot_precond_fermi_low_rank, &
1465 : ot_precond_full_all_covariant, &
1466 : ot_precond_full_single) ! these depend on the ks matrix
1467 3428 : DO ispin = 1, SIZE(preconditioner)
1468 1916 : CALL destroy_preconditioner(preconditioner(ispin)%preconditioner)
1469 3428 : DEALLOCATE (preconditioner(ispin)%preconditioner)
1470 : END DO
1471 1512 : DEALLOCATE (preconditioner)
1472 : CASE (ot_precond_none, ot_precond_full_kinetic, ot_precond_s_inverse, &
1473 : ot_precond_full_single_inverse) ! these are 'independent'
1474 : ! do nothing
1475 : CASE DEFAULT
1476 5188 : CPABORT("Unknown preconditioner type")
1477 : END SELECT
1478 : END IF
1479 :
1480 : ! add an OT preconditioner if none is present
1481 7565 : IF (.NOT. ASSOCIATED(preconditioner)) THEN
1482 7167 : SELECT CASE (prec_type)
1483 : CASE (ot_precond_full_all, ot_precond_fermi_low_rank, ot_precond_full_all_covariant, &
1484 : ot_precond_full_single_inverse)
1485 14077 : ALLOCATE (preconditioner(nspins))
1486 : CASE DEFAULT
1487 4500 : ALLOCATE (preconditioner(1))
1488 : END SELECT
1489 8743 : DO ispin = 1, SIZE(preconditioner)
1490 24270 : ALLOCATE (preconditioner(ispin)%preconditioner)
1491 : CALL init_preconditioner(preconditioner(ispin)%preconditioner, &
1492 : para_env=para_env, &
1493 8743 : blacs_env=blacs_env)
1494 : END DO
1495 : END IF
1496 :
1497 7565 : 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 7565 : 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 7565 : 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 7565 : 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 7565 : 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 7565 : CALL timeset(routineN, handle)
1560 7565 : NULLIFY (matrix_t, mo_coeff_b, mo_coeff, kinetic, dft_control, para_env, blacs_env, &
1561 7565 : input, subsys_section, multiple_unit_cell)
1562 7565 : my_has_unit_metric = .FALSE.
1563 7565 : IF (PRESENT(has_unit_metric)) my_has_unit_metric = has_unit_metric
1564 7565 : my_convert_to_dbcsr = .TRUE.
1565 7565 : IF (PRESENT(convert_to_dbcsr)) my_convert_to_dbcsr = convert_to_dbcsr
1566 7565 : my_full_mo_set = .FALSE.
1567 7565 : IF (PRESENT(full_mo_set)) my_full_mo_set = full_mo_set
1568 7565 : my_chebyshev_degree = 8
1569 7565 : IF (PRESENT(chebyshev_degree)) my_chebyshev_degree = chebyshev_degree
1570 7565 : my_low_rank_base = ot_low_rank_base_overlap
1571 7565 : IF (PRESENT(low_rank_base)) my_low_rank_base = low_rank_base
1572 7565 : my_fermi_low_rank_max_rank = 48
1573 7565 : IF (PRESENT(fermi_low_rank_max_rank)) my_fermi_low_rank_max_rank = fermi_low_rank_max_rank
1574 7565 : my_lattice_fft = ot_lattice_fft_off
1575 7565 : IF (PRESENT(lattice_fft)) my_lattice_fft = lattice_fft
1576 7565 : my_lattice_fft_local_cells = 0
1577 7565 : IF (PRESENT(lattice_fft_local_cells)) my_lattice_fft_local_cells = lattice_fft_local_cells
1578 7547 : IF (my_lattice_fft_local_cells < 0) CPABORT("LATTICE_FFT_LOCAL_CELLS must not be negative")
1579 7565 : do_kpoints = .FALSE.
1580 7565 : 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 7565 : blacs_env=blacs_env)
1588 :
1589 7565 : IF (my_low_rank_base == ot_low_rank_base_lattice_fft) THEN
1590 0 : IF (prec_type /= ot_precond_fermi_low_rank) THEN
1591 0 : CPABORT("FERMI_LOW_RANK_BASE LATTICE_FFT requires PRECONDITIONER FERMI_LOW_RANK")
1592 : END IF
1593 : my_lattice_fft = ot_lattice_fft_on
1594 7565 : ELSE IF (my_low_rank_base /= ot_low_rank_base_overlap) THEN
1595 0 : CPABORT("Unknown FERMI_LOW_RANK_BASE setting")
1596 : END IF
1597 :
1598 : IF (my_lattice_fft /= ot_lattice_fft_off .AND. &
1599 7565 : my_lattice_fft /= ot_lattice_fft_on .AND. &
1600 : my_lattice_fft /= ot_lattice_fft_auto) THEN
1601 0 : CPABORT("Unknown LATTICE_FFT setting")
1602 : END IF
1603 7565 : IF (my_lattice_fft_local_cells > 0 .AND. my_lattice_fft /= ot_lattice_fft_on) THEN
1604 0 : CPABORT("LATTICE_FFT_LOCAL_CELLS requires LATTICE_FFT ON")
1605 : END IF
1606 7565 : IF (my_lattice_fft /= ot_lattice_fft_off) THEN
1607 0 : SELECT CASE (prec_type)
1608 : CASE (ot_precond_fermi_low_rank, ot_precond_full_all_covariant, &
1609 : ot_precond_full_kinetic, ot_precond_full_single_inverse, ot_precond_s_inverse)
1610 : CASE DEFAULT
1611 2 : CPABORT("LATTICE_FFT does not support the selected OT preconditioner")
1612 : END SELECT
1613 2 : IF (solver_type /= ot_precond_solver_default .AND. &
1614 : solver_type /= ot_precond_solver_inv_chol) THEN
1615 0 : CPABORT("LATTICE_FFT requires PRECOND_SOLVER DEFAULT or INVERSE_CHOLESKY")
1616 : END IF
1617 2 : IF (my_has_unit_metric .OR. do_kpoints .OR. .NOT. ASSOCIATED(input)) THEN
1618 0 : IF (my_lattice_fft == ot_lattice_fft_on) THEN
1619 0 : 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 2 : subsys_section => section_vals_get_subs_vals(input, "SUBSYS")
1626 2 : IF (.NOT. ASSOCIATED(subsys_section)) THEN
1627 0 : 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 2 : i_vals=multiple_unit_cell)
1633 2 : IF (.NOT. ASSOCIATED(multiple_unit_cell)) THEN
1634 0 : IF (my_lattice_fft == ot_lattice_fft_on) THEN
1635 0 : 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 8 : IF (SIZE(multiple_unit_cell) /= 3 .OR. ANY(multiple_unit_cell <= 0)) THEN
1640 0 : IF (my_lattice_fft == ot_lattice_fft_on) THEN
1641 0 : CPABORT("LATTICE_FFT ON requires a valid CELL%MULTIPLE_UNIT_CELL")
1642 : END IF
1643 : my_lattice_fft = ot_lattice_fft_off
1644 8 : ELSE IF (PRODUCT(multiple_unit_cell) <= 1 .AND. my_lattice_fft == ot_lattice_fft_on) THEN
1645 0 : CPABORT("LATTICE_FFT ON requires an explicitly replicated multiple unit cell")
1646 : END IF
1647 8 : IF (my_lattice_fft_local_cells >= PRODUCT(multiple_unit_cell)) THEN
1648 0 : 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 7565 : IF (dft_control%qs_control%semi_empirical .OR. dft_control%qs_control%dftb .OR. &
1655 : dft_control%qs_control%xtb) THEN
1656 2146 : IF (prec_type == ot_precond_full_kinetic) THEN
1657 0 : msg = "Full_kinetic not available for semi-empirical methods"
1658 0 : CPABORT(TRIM(msg))
1659 : END IF
1660 2146 : matrix_t => matrix_s(1)%matrix
1661 : ELSE
1662 5419 : CPASSERT(.NOT. my_has_unit_metric)
1663 5419 : CALL get_qs_env(qs_env, kinetic=kinetic)
1664 5419 : matrix_t => kinetic(1)%matrix
1665 : END IF
1666 :
1667 : ! use full set of MOs or just occupied MOs
1668 7565 : nocc = 0
1669 7565 : norb = 0
1670 7565 : IF (my_full_mo_set) THEN
1671 38 : DO ispin = 1, nspins
1672 20 : CALL get_mo_set(mo_set=mos(ispin), homo=nocc(ispin), nmo=norb(ispin))
1673 38 : 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 7565 : n_loops = 1
1678 7565 : IF (prec_type == ot_precond_full_single_inverse) n_loops = nspins
1679 : ! check whether we need the ev and rotate the MOs
1680 2576 : SELECT CASE (prec_type)
1681 : CASE (ot_precond_full_all)
1682 : ! FULL_ALL needs canonical occupied eigenstates and one preconditioner per spin.
1683 2576 : n_loops = nspins
1684 :
1685 2576 : do_co_rotate = ASSOCIATED(qs_env%mo_derivs)
1686 5956 : DO ispin = 1, nspins
1687 3380 : CALL get_mo_set(mo_set=mos(ispin), mo_coeff_b=mo_coeff_b, mo_coeff=mo_coeff)
1688 3380 : use_mo_coeff_b = mos(ispin)%use_mo_coeff_b
1689 5956 : 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 3360 : blacs_env=blacs_env)
1695 20 : 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 20 : blacs_env=blacs_env)
1700 : ELSE
1701 : CALL calculate_subspace_eigenvalues(mo_coeff, matrix_ks(ispin)%matrix, &
1702 0 : do_rotation=.TRUE.)
1703 : END IF
1704 : END DO
1705 : CASE (ot_precond_fermi_low_rank)
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 0 : n_loops = nspins
1709 : CASE (ot_precond_full_all_covariant)
1710 : ! Diagonalize a private occupied copy. The physical MOs and their
1711 : ! derivatives remain in the user-selected gauge.
1712 0 : n_loops = nspins
1713 0 : ALLOCATE (canonical_mos(nspins))
1714 7565 : DO ispin = 1, nspins
1715 0 : CALL duplicate_mo_set(canonical_mos(ispin), mos(ispin))
1716 0 : CALL get_mo_set(mo_set=canonical_mos(ispin), mo_coeff_b=mo_coeff_b, mo_coeff=mo_coeff)
1717 0 : IF (canonical_mos(ispin)%use_mo_coeff_b) THEN
1718 0 : CALL copy_dbcsr_to_fm(mo_coeff_b, mo_coeff)
1719 0 : 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 0 : 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 718 : SELECT CASE (prec_type)
1731 : CASE (ot_precond_none)
1732 1436 : DO ispin = 1, SIZE(ot_preconditioner)
1733 1436 : ot_preconditioner(ispin)%preconditioner%in_use = 0
1734 : END DO
1735 : CASE DEFAULT
1736 22548 : DO icall = 1, n_loops
1737 8136 : ot_preconditioner(icall)%preconditioner%polynomial_degree = my_chebyshev_degree
1738 8136 : ot_preconditioner(icall)%preconditioner%low_rank_base = my_low_rank_base
1739 8136 : ot_preconditioner(icall)%preconditioner%lattice_fft_mode = my_lattice_fft
1740 8136 : ot_preconditioner(icall)%preconditioner%lattice_local_cells = my_lattice_fft_local_cells
1741 32544 : ot_preconditioner(icall)%preconditioner%lattice_dims = [1, 1, 1]
1742 8136 : IF (my_lattice_fft /= ot_lattice_fft_off) THEN
1743 2 : CPASSERT(ASSOCIATED(multiple_unit_cell))
1744 16 : ot_preconditioner(icall)%preconditioner%lattice_dims = multiple_unit_cell
1745 : END IF
1746 8136 : ot_preconditioner(icall)%preconditioner%spectral_max_rank = my_fermi_low_rank_max_rank
1747 14983 : IF (my_has_unit_metric) THEN
1748 474 : 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 0 : 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 474 : convert_precond_to_dbcsr=my_convert_to_dbcsr)
1762 : END IF
1763 : ELSE
1764 7662 : 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 0 : 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 7662 : 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 7565 : IF (ALLOCATED(canonical_mos)) THEN
1785 0 : DO ispin = 1, SIZE(canonical_mos)
1786 0 : CALL deallocate_mo_set(canonical_mos(ispin))
1787 : END DO
1788 0 : DEALLOCATE (canonical_mos)
1789 : END IF
1790 :
1791 : ! reset homo values
1792 7565 : IF (my_full_mo_set) THEN
1793 38 : DO ispin = 1, nspins
1794 38 : CALL set_mo_set(mo_set=mos(ispin), homo=nocc(ispin))
1795 : END DO
1796 : END IF
1797 :
1798 7565 : CALL timestop(handle)
1799 :
1800 15130 : 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 0 : 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 0 : 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 0 : CALL timeset(routineN, handle)
1824 0 : CALL cp_fm_get_info(mo_coeff, nrow_global=nao)
1825 0 : CPASSERT(nocc > 0)
1826 :
1827 : CALL cp_fm_struct_create(fm_struct, nrow_global=nao, ncol_global=nocc, &
1828 : context=preconditioner_env%ctxt, &
1829 0 : para_env=preconditioner_env%para_env)
1830 0 : CALL cp_fm_create(canonical_coeff, fm_struct, name=routineN)
1831 0 : CALL cp_fm_create(weighted_coeff, fm_struct, name=routineN)
1832 0 : CALL cp_fm_struct_release(fm_struct)
1833 0 : 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 0 : para_env=preconditioner_env%para_env)
1838 0 : CALL cp_fm_create(h_block, fm_struct, name=routineN)
1839 0 : IF (ASSOCIATED(preconditioner_env%occ_rotation_dbcsr)) THEN
1840 0 : CALL dbcsr_release_p(preconditioner_env%occ_rotation_dbcsr)
1841 : END IF
1842 0 : IF (ASSOCIATED(preconditioner_env%occ_rotation)) THEN
1843 0 : CALL cp_fm_release(preconditioner_env%occ_rotation)
1844 0 : DEALLOCATE (preconditioner_env%occ_rotation)
1845 : END IF
1846 0 : ALLOCATE (preconditioner_env%occ_rotation)
1847 0 : CALL cp_fm_create(preconditioner_env%occ_rotation, fm_struct, name="occupied rotation")
1848 0 : CALL cp_fm_struct_release(fm_struct)
1849 :
1850 0 : 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 0 : weighted_coeff, 0.0_dp, h_block)
1853 0 : ALLOCATE (occupied_evals(nocc))
1854 0 : 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 0 : preconditioner_env%occ_rotation, 0.0_dp, weighted_coeff)
1858 0 : CALL cp_fm_to_fm(weighted_coeff, mo_coeff, nocc, 1, 1)
1859 :
1860 0 : DEALLOCATE (occupied_evals)
1861 0 : CALL cp_fm_release(canonical_coeff)
1862 0 : CALL cp_fm_release(h_block)
1863 0 : CALL cp_fm_release(weighted_coeff)
1864 0 : CALL timestop(handle)
1865 :
1866 0 : END SUBROUTINE canonicalize_full_all_copy
1867 :
1868 0 : END MODULE preconditioner
|
