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_types, ONLY: cp_cfm_create,&
18 : cp_cfm_release,&
19 : cp_cfm_type,&
20 : cp_fm_to_cfm
21 : USE cp_control_types, ONLY: dft_control_type
22 : USE cp_dbcsr_api, ONLY: dbcsr_get_info,&
23 : dbcsr_p_type,&
24 : dbcsr_release_p,&
25 : dbcsr_type
26 : USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
27 : cp_dbcsr_sm_fm_multiply
28 : USE cp_fm_basic_linalg, ONLY: cp_fm_uplo_to_full
29 : USE cp_fm_cholesky, ONLY: cp_fm_cholesky_decompose,&
30 : cp_fm_cholesky_invert
31 : USE cp_fm_diag, ONLY: choose_eigv_solver
32 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
33 : cp_fm_struct_release,&
34 : cp_fm_struct_type
35 : USE cp_fm_types, ONLY: cp_fm_create,&
36 : cp_fm_get_info,&
37 : cp_fm_release,&
38 : cp_fm_set_all,&
39 : cp_fm_to_fm,&
40 : cp_fm_type
41 : USE cp_log_handling, ONLY: cp_get_default_logger,&
42 : cp_logger_get_default_io_unit,&
43 : cp_logger_type
44 : USE input_constants, ONLY: &
45 : cholesky_reduce, ot_lattice_fft_auto, ot_lattice_fft_off, ot_lattice_fft_on, &
46 : ot_low_rank_base_lattice_fft, ot_low_rank_base_overlap, ot_precond_fermi_low_rank, &
47 : ot_precond_full_all, ot_precond_full_all_covariant, ot_precond_full_kinetic, &
48 : ot_precond_full_single, ot_precond_full_single_inverse, ot_precond_none, &
49 : ot_precond_s_inverse, ot_precond_solver_default, ot_precond_solver_inv_chol, &
50 : ot_precond_solver_update
51 : USE input_section_types, ONLY: section_vals_get_subs_vals,&
52 : section_vals_type,&
53 : section_vals_val_get
54 : USE kinds, ONLY: default_string_length,&
55 : dp
56 : USE lattice_preconditioner_operator, ONLY: &
57 : build_lattice_inverse, build_lattice_inverse_from_blocks, build_lattice_local_correction, &
58 : lattice_fft_auto_max_cost_ratio, lattice_fft_auto_max_storage_ratio, &
59 : lattice_fft_auto_min_cells, lattice_fft_auto_select, lattice_fft_cost_ratio, &
60 : lattice_fft_storage_ratio, lattice_grid, lattice_local_cost_ratio, &
61 : lattice_local_storage_ratio
62 : USE low_rank_preconditioner_model, ONLY: low_rank_inverse_weight,&
63 : low_rank_select_rank
64 : USE message_passing, ONLY: mp_para_env_type
65 : USE parallel_gemm_api, ONLY: parallel_gemm
66 : USE preconditioner_apply, ONLY: apply_preconditioner_cfm_complex,&
67 : apply_preconditioner_dbcsr,&
68 : apply_preconditioner_dbcsr_complex,&
69 : apply_preconditioner_fm
70 : USE preconditioner_makes, ONLY: make_complex_fermi_low_rank,&
71 : make_complex_full_all,&
72 : make_complex_full_all_covariant,&
73 : make_complex_full_kinetic,&
74 : make_complex_full_s_inverse,&
75 : make_complex_full_single,&
76 : make_complex_full_single_inverse,&
77 : make_preconditioner_matrix
78 : USE preconditioner_solvers, ONLY: solve_preconditioner,&
79 : transfer_dbcsr_to_fm,&
80 : transfer_fm_to_dbcsr
81 : USE preconditioner_types, ONLY: destroy_preconditioner,&
82 : init_preconditioner,&
83 : preconditioner_p_type,&
84 : preconditioner_type
85 : USE qs_environment_types, ONLY: get_qs_env,&
86 : qs_environment_type
87 : USE qs_mo_methods, ONLY: calculate_subspace_eigenvalues
88 : USE qs_mo_types, ONLY: deallocate_mo_set,&
89 : duplicate_mo_set,&
90 : get_mo_set,&
91 : mo_set_type,&
92 : set_mo_set
93 : #include "./base/base_uses.f90"
94 :
95 : IMPLICIT NONE
96 :
97 : PRIVATE
98 :
99 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'preconditioner'
100 : REAL(KIND=dp), PARAMETER, PRIVATE :: fermi_low_rank_gap_factor = 4.0_dp, &
101 : fermi_low_rank_degeneracy_tolerance = 1.0E-8_dp, &
102 : fermi_low_rank_min_window = 1.0_dp
103 :
104 : PUBLIC :: make_preconditioner, make_preconditioner_complex_fermi_low_rank, &
105 : make_preconditioner_complex_full_all, &
106 : make_preconditioner_complex_full_all_covariant, &
107 : make_preconditioner_complex_full_kinetic, &
108 : make_preconditioner_complex_full_s_inverse, &
109 : make_preconditioner_complex_full_single, &
110 : make_preconditioner_complex_full_single_inverse, restart_preconditioner, &
111 : dbcsr_pair_to_cfm
112 : PUBLIC :: apply_preconditioner, prepare_preconditioner
113 :
114 : ! The public interface for apply preconditioner, the routines can be found in preconditioner_apply.F
115 : INTERFACE apply_preconditioner
116 : MODULE PROCEDURE apply_preconditioner_cfm_complex
117 : MODULE PROCEDURE apply_preconditioner_dbcsr
118 : MODULE PROCEDURE apply_preconditioner_dbcsr_complex
119 : MODULE PROCEDURE apply_preconditioner_fm
120 : END INTERFACE
121 :
122 : ! **************************************************************************************************
123 :
124 : CONTAINS
125 :
126 : ! **************************************************************************************************
127 :
128 : ! **************************************************************************************************
129 : !> \brief Build the block-circulant inverse of a replicated-cell SPD operator.
130 : !> \param preconditioner_env ...
131 : !> \param operator_matrix ...
132 : !> \param operator_name ...
133 : !> \param rhs_count number of occupied orbital columns to precondition
134 : ! **************************************************************************************************
135 16 : SUBROUTINE prepare_lattice_operator(preconditioner_env, operator_matrix, operator_name, rhs_count)
136 :
137 : TYPE(preconditioner_type) :: preconditioner_env
138 : TYPE(dbcsr_type), POINTER :: operator_matrix
139 : CHARACTER(LEN=*), INTENT(IN) :: operator_name
140 : INTEGER, INTENT(IN) :: rhs_count
141 :
142 : INTEGER :: alpha, beta, block_size, col_cell, delta_index, i, info, j, local_info, n, ncell, &
143 : ncol_local, nrow_local, output_unit, row_cell
144 16 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: index_to_cell
145 : INTEGER, DIMENSION(3) :: delta
146 16 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
147 : LOGICAL :: used_fft
148 : REAL(KIND=dp) :: denominator, difference, numerator, value
149 16 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: dense_operator, xkp
150 16 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: blocks
151 : REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :), &
152 16 : POINTER :: local_data
153 : TYPE(cp_fm_struct_type), POINTER :: fm_struct
154 : TYPE(cp_fm_type) :: operator_fm, work
155 : TYPE(cp_logger_type), POINTER :: logger
156 :
157 0 : CPASSERT(ASSOCIATED(operator_matrix))
158 16 : CPASSERT(preconditioner_env%lattice_fft_mode /= ot_lattice_fft_off)
159 16 : preconditioner_env%lattice_fft_active = .FALSE.
160 16 : preconditioner_env%lattice_projection_error = 0.0_dp
161 16 : preconditioner_env%lattice_corrected_projection_error = 0.0_dp
162 16 : preconditioner_env%lattice_cost_ratio = 0.0_dp
163 16 : preconditioner_env%lattice_storage_ratio = 0.0_dp
164 :
165 64 : IF (ANY(preconditioner_env%lattice_dims <= 0)) CPABORT("Lattice dimensions must be positive")
166 :
167 64 : IF (preconditioner_env%lattice_fft_mode == ot_lattice_fft_auto .AND. &
168 : PRODUCT(preconditioner_env%lattice_dims) < lattice_fft_auto_min_cells) THEN
169 0 : IF (preconditioner_env%para_env%is_source()) THEN
170 0 : logger => cp_get_default_logger()
171 0 : output_unit = cp_logger_get_default_io_unit(logger)
172 0 : IF (output_unit > 0) THEN
173 : WRITE (output_unit, '(T2,A,A,A,I0,A,I0)') &
174 0 : "OT| Lattice FFT AUTO kept the conventional ", TRIM(operator_name), &
175 0 : " solver: ", PRODUCT(preconditioner_env%lattice_dims), &
176 0 : " cells are fewer than ", lattice_fft_auto_min_cells
177 : END IF
178 : END IF
179 0 : RETURN
180 : END IF
181 :
182 16 : NULLIFY (fm_struct)
183 16 : CALL dbcsr_get_info(operator_matrix, nfullrows_total=n)
184 64 : ncell = PRODUCT(preconditioner_env%lattice_dims)
185 16 : IF (MODULO(n, ncell) /= 0) CPABORT("AO dimension does not match replicated cells")
186 16 : block_size = n/ncell
187 16 : preconditioner_env%lattice_block_size = block_size
188 : preconditioner_env%lattice_cost_ratio = lattice_fft_cost_ratio( &
189 : preconditioner_env%lattice_dims, block_size, rhs_count, &
190 : preconditioner_env%para_env%num_pe) + &
191 : lattice_local_cost_ratio( &
192 : preconditioner_env%lattice_dims, &
193 16 : preconditioner_env%lattice_local_cells)
194 : preconditioner_env%lattice_storage_ratio = lattice_fft_storage_ratio( &
195 : preconditioner_env%lattice_dims, 1) + &
196 : lattice_local_storage_ratio( &
197 : preconditioner_env%lattice_dims, 1, &
198 16 : preconditioner_env%lattice_local_cells)
199 16 : IF (preconditioner_env%lattice_fft_mode == ot_lattice_fft_auto .AND. &
200 : (preconditioner_env%lattice_cost_ratio > lattice_fft_auto_max_cost_ratio .OR. &
201 : preconditioner_env%lattice_storage_ratio > lattice_fft_auto_max_storage_ratio)) THEN
202 2 : IF (preconditioner_env%para_env%is_source()) THEN
203 1 : logger => cp_get_default_logger()
204 1 : output_unit = cp_logger_get_default_io_unit(logger)
205 1 : IF (output_unit > 0) THEN
206 : WRITE (output_unit, '(T2,A,A,A,F8.3,A,F8.3,A)') &
207 1 : "OT| Lattice FFT AUTO kept the conventional ", TRIM(operator_name), &
208 1 : " solver: estimated work ratio ", preconditioner_env%lattice_cost_ratio, &
209 2 : " storage ratio ", preconditioner_env%lattice_storage_ratio, " selected: F"
210 : END IF
211 : END IF
212 2 : RETURN
213 : END IF
214 : CALL cp_fm_struct_create(fm_struct, nrow_global=n, ncol_global=n, &
215 : context=preconditioner_env%ctxt, &
216 14 : para_env=preconditioner_env%para_env)
217 14 : CALL cp_fm_create(operator_fm, fm_struct, name="lattice preconditioner operator")
218 14 : CALL cp_fm_struct_release(fm_struct)
219 14 : CALL cp_fm_set_all(operator_fm, alpha=0.0_dp, beta=0.0_dp)
220 14 : CALL copy_dbcsr_to_fm(operator_matrix, operator_fm)
221 14 : CALL cp_fm_create(work, operator_fm%matrix_struct, name="lattice operator work")
222 14 : CALL cp_fm_uplo_to_full(operator_fm, work)
223 14 : CALL cp_fm_release(work)
224 : CALL cp_fm_get_info(operator_fm, nrow_local=nrow_local, ncol_local=ncol_local, &
225 14 : row_indices=row_indices, col_indices=col_indices, local_data=local_data)
226 70 : ALLOCATE (blocks(block_size, block_size, ncell), source=0.0_dp)
227 70 : ALLOCATE (index_to_cell(3, ncell), xkp(3, ncell))
228 14 : CALL lattice_grid(preconditioner_env%lattice_dims, index_to_cell, xkp)
229 14 : IF (preconditioner_env%lattice_local_cells > 0) THEN
230 8 : ALLOCATE (dense_operator(n, n), source=0.0_dp)
231 : END IF
232 14 : denominator = 0.0_dp
233 338 : DO j = 1, ncol_local
234 7094 : DO i = 1, nrow_local
235 7080 : IF (row_indices(i) <= col_indices(j)) THEN
236 3459 : value = local_data(i, j)
237 3459 : row_cell = (row_indices(i) - 1)/block_size + 1
238 3459 : col_cell = (col_indices(j) - 1)/block_size + 1
239 3459 : alpha = MODULO(row_indices(i) - 1, block_size) + 1
240 3459 : beta = MODULO(col_indices(j) - 1, block_size) + 1
241 : delta = MODULO(index_to_cell(:, row_cell) - index_to_cell(:, col_cell), &
242 13836 : preconditioner_env%lattice_dims)
243 : delta_index = delta(1) + preconditioner_env%lattice_dims(1)*(delta(2) + &
244 3459 : preconditioner_env%lattice_dims(2)*delta(3)) + 1
245 : blocks(alpha, beta, delta_index) = blocks(alpha, beta, delta_index) + &
246 3459 : value/REAL(ncell, KIND=dp)
247 3459 : denominator = denominator + value*value
248 3459 : IF (ALLOCATED(dense_operator)) THEN
249 36 : dense_operator(row_indices(i), col_indices(j)) = value
250 : END IF
251 3459 : IF (row_indices(i) < col_indices(j)) THEN
252 : delta = MODULO(index_to_cell(:, col_cell) - index_to_cell(:, row_cell), &
253 13188 : preconditioner_env%lattice_dims)
254 : delta_index = delta(1) + preconditioner_env%lattice_dims(1)*(delta(2) + &
255 3297 : preconditioner_env%lattice_dims(2)*delta(3)) + 1
256 : blocks(beta, alpha, delta_index) = blocks(beta, alpha, delta_index) + &
257 3297 : value/REAL(ncell, KIND=dp)
258 3297 : denominator = denominator + value*value
259 : END IF
260 : END IF
261 : END DO
262 : END DO
263 14 : CALL preconditioner_env%para_env%sum(blocks)
264 14 : numerator = 0.0_dp
265 338 : DO j = 1, ncol_local
266 7094 : DO i = 1, nrow_local
267 7080 : IF (row_indices(i) <= col_indices(j)) THEN
268 3459 : value = local_data(i, j)
269 3459 : row_cell = (row_indices(i) - 1)/block_size + 1
270 3459 : col_cell = (col_indices(j) - 1)/block_size + 1
271 3459 : alpha = MODULO(row_indices(i) - 1, block_size) + 1
272 3459 : beta = MODULO(col_indices(j) - 1, block_size) + 1
273 : delta = MODULO(index_to_cell(:, row_cell) - index_to_cell(:, col_cell), &
274 13836 : preconditioner_env%lattice_dims)
275 : delta_index = delta(1) + preconditioner_env%lattice_dims(1)*(delta(2) + &
276 3459 : preconditioner_env%lattice_dims(2)*delta(3)) + 1
277 3459 : difference = value - blocks(alpha, beta, delta_index)
278 3459 : numerator = numerator + difference*difference
279 3459 : IF (row_indices(i) < col_indices(j)) THEN
280 : delta = MODULO(index_to_cell(:, col_cell) - index_to_cell(:, row_cell), &
281 13188 : preconditioner_env%lattice_dims)
282 : delta_index = delta(1) + preconditioner_env%lattice_dims(1)*(delta(2) + &
283 3297 : preconditioner_env%lattice_dims(2)*delta(3)) + 1
284 3297 : difference = value - blocks(beta, alpha, delta_index)
285 3297 : numerator = numerator + difference*difference
286 : END IF
287 : END IF
288 : END DO
289 : END DO
290 14 : CALL preconditioner_env%para_env%sum(numerator)
291 14 : CALL preconditioner_env%para_env%sum(denominator)
292 : preconditioner_env%lattice_projection_error = &
293 14 : SQRT(numerator/MAX(denominator, TINY(denominator)))
294 : CALL build_lattice_inverse_from_blocks( &
295 14 : blocks, preconditioner_env%lattice_dims, preconditioner_env%lattice_inverse_k, used_fft, info)
296 : preconditioner_env%lattice_corrected_projection_error = &
297 14 : preconditioner_env%lattice_projection_error
298 14 : IF (ALLOCATED(dense_operator)) THEN
299 2 : CALL preconditioner_env%para_env%sum(dense_operator)
300 18 : DO j = 1, n
301 74 : DO i = j + 1, n
302 72 : dense_operator(i, j) = dense_operator(j, i)
303 : END DO
304 : END DO
305 : CALL build_lattice_local_correction( &
306 : dense_operator, preconditioner_env%lattice_dims, &
307 : preconditioner_env%lattice_local_cells, &
308 : preconditioner_env%lattice_local_indices, &
309 : preconditioner_env%lattice_local_operator, &
310 : preconditioner_env%lattice_local_coarse_inverse, &
311 2 : preconditioner_env%lattice_corrected_projection_error, local_info)
312 2 : IF (local_info /= 0) THEN
313 0 : CPABORT("LATTICE_FFT local correction found a non-positive coarse operator")
314 : END IF
315 2 : DEALLOCATE (dense_operator)
316 : END IF
317 14 : CALL cp_fm_release(operator_fm)
318 14 : DEALLOCATE (blocks, index_to_cell, xkp)
319 24 : SELECT CASE (preconditioner_env%lattice_fft_mode)
320 : CASE (ot_lattice_fft_on)
321 10 : IF (info /= 0) THEN
322 0 : CPABORT("LATTICE_FFT found a non-positive Fourier block for "//TRIM(operator_name))
323 : END IF
324 10 : IF (.NOT. used_fft) THEN
325 0 : CPABORT("LATTICE_FFT could not initialize its FFT for "//TRIM(operator_name))
326 : END IF
327 10 : preconditioner_env%lattice_fft_active = .TRUE.
328 : CASE (ot_lattice_fft_auto)
329 : preconditioner_env%lattice_fft_active = lattice_fft_auto_select( &
330 : preconditioner_env%lattice_dims, &
331 : preconditioner_env%lattice_block_size, rhs_count, 1, &
332 : preconditioner_env%para_env%num_pe, &
333 4 : preconditioner_env%lattice_projection_error, used_fft, info)
334 : CASE DEFAULT
335 14 : CPABORT("Unknown LATTICE_FFT mode")
336 : END SELECT
337 :
338 14 : IF (.NOT. preconditioner_env%lattice_fft_active .AND. &
339 : ALLOCATED(preconditioner_env%lattice_inverse_k)) THEN
340 2 : DEALLOCATE (preconditioner_env%lattice_inverse_k)
341 : END IF
342 14 : IF (preconditioner_env%para_env%is_source()) THEN
343 7 : logger => cp_get_default_logger()
344 7 : output_unit = cp_logger_get_default_io_unit(logger)
345 7 : IF (output_unit > 0) THEN
346 : WRITE (output_unit, '(T2,A,A,A,3(I0,1X),A,I0)') &
347 7 : "OT| Lattice FFT ", TRIM(operator_name), " cells: ", &
348 14 : preconditioner_env%lattice_dims, "AO block: ", preconditioner_env%lattice_block_size
349 : WRITE (output_unit, '(T2,A,ES12.4,A,L1)') &
350 7 : "OT| Block-circulant projection error: ", &
351 7 : preconditioner_env%lattice_projection_error, &
352 14 : " selected: ", preconditioner_env%lattice_fft_active
353 7 : IF (preconditioner_env%lattice_local_cells > 0) THEN
354 : WRITE (output_unit, '(T2,A,I0,A,ES12.4)') &
355 1 : "OT| Balanced local correction cells: ", &
356 1 : preconditioner_env%lattice_local_cells, " residual error: ", &
357 2 : preconditioner_env%lattice_corrected_projection_error
358 : END IF
359 : WRITE (output_unit, '(T2,A,F8.3,A,F8.3)') &
360 7 : "OT| Estimated lattice/dense work ratio: ", &
361 7 : preconditioner_env%lattice_cost_ratio, " storage ratio: ", &
362 14 : preconditioner_env%lattice_storage_ratio
363 : END IF
364 : END IF
365 :
366 90 : END SUBROUTINE prepare_lattice_operator
367 :
368 : ! **************************************************************************************************
369 : !> \brief Build state-dependent block-circulant inverses for rotation-covariant FULL_ALL.
370 : !> \param preconditioner_env ...
371 : !> \param matrix_s overlap matrix
372 : !> \param operator_name ...
373 : ! **************************************************************************************************
374 2 : SUBROUTINE prepare_lattice_state_operator(preconditioner_env, matrix_s, operator_name)
375 :
376 : TYPE(preconditioner_type) :: preconditioner_env
377 : TYPE(dbcsr_type), POINTER :: matrix_s
378 : CHARACTER(LEN=*), INTENT(IN) :: operator_name
379 :
380 2 : COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: inverse_k
381 : INTEGER :: i, info, istate, j, local_info, &
382 : local_rank, n, ncell, ncol_local, &
383 : nrow_local, nstate, output_unit
384 2 : INTEGER, ALLOCATABLE, DIMENSION(:) :: local_indices
385 2 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
386 : LOGICAL :: all_selected, selected, used_fft
387 : REAL(KIND=dp) :: corrected_projection_error, &
388 : projection_error
389 2 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: dense_operator, dense_sv, dense_vectors, &
390 2 : local_coarse_inverse, local_operator, &
391 2 : scaled_sv
392 : REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :), &
393 2 : POINTER :: local_data
394 : TYPE(cp_fm_type) :: sv_fm
395 : TYPE(cp_logger_type), POINTER :: logger
396 :
397 0 : CPASSERT(ASSOCIATED(matrix_s))
398 2 : CPASSERT(ASSOCIATED(preconditioner_env%fm))
399 2 : CPASSERT(ASSOCIATED(preconditioner_env%full_evals))
400 2 : CPASSERT(ASSOCIATED(preconditioner_env%occ_evals))
401 2 : CPASSERT(preconditioner_env%lattice_fft_mode /= ot_lattice_fft_off)
402 :
403 2 : preconditioner_env%lattice_fft_active = .FALSE.
404 2 : preconditioner_env%lattice_projection_error = 0.0_dp
405 2 : preconditioner_env%lattice_corrected_projection_error = 0.0_dp
406 2 : preconditioner_env%lattice_cost_ratio = 0.0_dp
407 2 : preconditioner_env%lattice_storage_ratio = 0.0_dp
408 2 : IF (ALLOCATED(preconditioner_env%lattice_state_inverse_k)) THEN
409 0 : DEALLOCATE (preconditioner_env%lattice_state_inverse_k)
410 : END IF
411 2 : IF (ALLOCATED(preconditioner_env%lattice_occ_vectors)) THEN
412 0 : DEALLOCATE (preconditioner_env%lattice_occ_vectors)
413 : END IF
414 2 : IF (ALLOCATED(preconditioner_env%lattice_occ_dual)) THEN
415 0 : DEALLOCATE (preconditioner_env%lattice_occ_dual)
416 : END IF
417 2 : IF (ALLOCATED(preconditioner_env%lattice_local_indices)) THEN
418 0 : DEALLOCATE (preconditioner_env%lattice_local_indices)
419 : END IF
420 2 : IF (ALLOCATED(preconditioner_env%lattice_local_operator)) THEN
421 0 : DEALLOCATE (preconditioner_env%lattice_local_operator)
422 : END IF
423 2 : IF (ALLOCATED(preconditioner_env%lattice_local_coarse_inverse)) THEN
424 0 : DEALLOCATE (preconditioner_env%lattice_local_coarse_inverse)
425 : END IF
426 2 : IF (ALLOCATED(preconditioner_env%lattice_state_local_indices)) THEN
427 0 : DEALLOCATE (preconditioner_env%lattice_state_local_indices)
428 : END IF
429 2 : IF (ALLOCATED(preconditioner_env%lattice_state_local_operator)) THEN
430 0 : DEALLOCATE (preconditioner_env%lattice_state_local_operator)
431 : END IF
432 2 : IF (ALLOCATED(preconditioner_env%lattice_state_local_coarse_inverse)) THEN
433 0 : DEALLOCATE (preconditioner_env%lattice_state_local_coarse_inverse)
434 : END IF
435 8 : IF (ANY(preconditioner_env%lattice_dims <= 0)) CPABORT("Lattice dimensions must be positive")
436 8 : IF (preconditioner_env%lattice_fft_mode == ot_lattice_fft_auto .AND. &
437 : PRODUCT(preconditioner_env%lattice_dims) < lattice_fft_auto_min_cells) THEN
438 0 : IF (preconditioner_env%para_env%is_source()) THEN
439 0 : logger => cp_get_default_logger()
440 0 : output_unit = cp_logger_get_default_io_unit(logger)
441 0 : IF (output_unit > 0) THEN
442 : WRITE (output_unit, '(T2,A,A,A,I0,A,I0)') &
443 0 : "OT| Lattice FFT AUTO kept the conventional ", TRIM(operator_name), &
444 0 : " solver: ", PRODUCT(preconditioner_env%lattice_dims), &
445 0 : " cells are fewer than ", lattice_fft_auto_min_cells
446 : END IF
447 : END IF
448 0 : RETURN
449 : END IF
450 :
451 2 : CALL cp_fm_get_info(preconditioner_env%fm, nrow_global=n, ncol_global=j)
452 2 : CPASSERT(j == n .AND. SIZE(preconditioner_env%full_evals) == n)
453 2 : nstate = SIZE(preconditioner_env%occ_evals)
454 2 : CPASSERT(nstate > 0)
455 8 : ncell = PRODUCT(preconditioner_env%lattice_dims)
456 2 : IF (MODULO(n, ncell) /= 0) CPABORT("AO dimension does not match replicated cells")
457 2 : preconditioner_env%lattice_block_size = n/ncell
458 : preconditioner_env%lattice_cost_ratio = lattice_fft_cost_ratio( &
459 : preconditioner_env%lattice_dims, &
460 : preconditioner_env%lattice_block_size, nstate, &
461 2 : preconditioner_env%para_env%num_pe)
462 : preconditioner_env%lattice_storage_ratio = lattice_fft_storage_ratio( &
463 2 : preconditioner_env%lattice_dims, nstate)
464 : preconditioner_env%lattice_cost_ratio = preconditioner_env%lattice_cost_ratio + &
465 : lattice_local_cost_ratio( &
466 : preconditioner_env%lattice_dims, &
467 2 : preconditioner_env%lattice_local_cells)
468 : preconditioner_env%lattice_storage_ratio = preconditioner_env%lattice_storage_ratio + &
469 : lattice_local_storage_ratio( &
470 : preconditioner_env%lattice_dims, nstate, &
471 2 : preconditioner_env%lattice_local_cells)
472 2 : IF (preconditioner_env%lattice_fft_mode == ot_lattice_fft_auto .AND. &
473 : (preconditioner_env%lattice_cost_ratio > lattice_fft_auto_max_cost_ratio .OR. &
474 : preconditioner_env%lattice_storage_ratio > lattice_fft_auto_max_storage_ratio)) THEN
475 0 : IF (preconditioner_env%para_env%is_source()) THEN
476 0 : logger => cp_get_default_logger()
477 0 : output_unit = cp_logger_get_default_io_unit(logger)
478 0 : IF (output_unit > 0) THEN
479 : WRITE (output_unit, '(T2,A,A,A,F8.3,A,F8.3,A)') &
480 0 : "OT| Lattice FFT AUTO kept the conventional ", TRIM(operator_name), &
481 0 : " solver: estimated work ratio ", preconditioner_env%lattice_cost_ratio, &
482 0 : " storage ratio ", preconditioner_env%lattice_storage_ratio, " selected: F"
483 : END IF
484 : END IF
485 0 : RETURN
486 : END IF
487 2 : IF (preconditioner_env%lattice_local_cells > 0) THEN
488 : local_rank = preconditioner_env%lattice_local_cells* &
489 2 : preconditioner_env%lattice_block_size
490 0 : ALLOCATE (preconditioner_env%lattice_state_local_indices(local_rank, nstate), &
491 0 : preconditioner_env%lattice_state_local_operator(n, local_rank, nstate), &
492 24 : preconditioner_env%lattice_state_local_coarse_inverse(local_rank, local_rank, nstate))
493 : END IF
494 2 : CALL cp_fm_create(sv_fm, preconditioner_env%fm%matrix_struct, name="lattice FULL_ALL S V")
495 2 : CALL cp_dbcsr_sm_fm_multiply(matrix_s, preconditioner_env%fm, sv_fm, n)
496 : CALL cp_fm_get_info(sv_fm, nrow_local=nrow_local, ncol_local=ncol_local, &
497 2 : row_indices=row_indices, col_indices=col_indices, local_data=local_data)
498 14 : ALLOCATE (dense_sv(n, n), dense_vectors(n, nstate), source=0.0_dp)
499 10 : DO j = 1, ncol_local
500 26 : DO i = 1, nrow_local
501 24 : dense_sv(row_indices(i), col_indices(j)) = local_data(i, j)
502 : END DO
503 : END DO
504 2 : CALL preconditioner_env%para_env%sum(dense_sv)
505 2 : CALL cp_fm_release(sv_fm)
506 :
507 : CALL cp_fm_get_info(preconditioner_env%fm, nrow_local=nrow_local, ncol_local=ncol_local, &
508 2 : row_indices=row_indices, col_indices=col_indices, local_data=local_data)
509 10 : DO j = 1, ncol_local
510 8 : IF (col_indices(j) > nstate) CYCLE
511 14 : DO i = 1, nrow_local
512 16 : dense_vectors(row_indices(i), col_indices(j)) = local_data(i, j)
513 : END DO
514 : END DO
515 2 : CALL preconditioner_env%para_env%sum(dense_vectors)
516 0 : ALLOCATE (preconditioner_env%lattice_occ_vectors(n, nstate), &
517 12 : preconditioner_env%lattice_occ_dual(n, nstate))
518 22 : preconditioner_env%lattice_occ_vectors(:, :) = dense_vectors
519 22 : preconditioner_env%lattice_occ_dual(:, :) = dense_sv(:, 1:nstate)
520 :
521 10 : ALLOCATE (dense_operator(n, n), scaled_sv(n, n))
522 2 : all_selected = .TRUE.
523 6 : DO istate = 1, nstate
524 84 : scaled_sv(:, :) = dense_sv
525 20 : DO j = 1, n
526 : scaled_sv(:, j) = scaled_sv(:, j)*MAX( &
527 : preconditioner_env%energy_gap, &
528 84 : preconditioner_env%full_evals(j) - preconditioner_env%occ_evals(istate))
529 : END DO
530 596 : dense_operator(:, :) = MATMUL(scaled_sv, TRANSPOSE(dense_sv))
531 : CALL build_lattice_inverse(dense_operator, preconditioner_env%lattice_dims, inverse_k, &
532 4 : projection_error, used_fft, info)
533 : preconditioner_env%lattice_projection_error = &
534 4 : MAX(preconditioner_env%lattice_projection_error, projection_error)
535 4 : corrected_projection_error = projection_error
536 4 : IF (preconditioner_env%lattice_local_cells > 0) THEN
537 : CALL build_lattice_local_correction( &
538 : dense_operator, preconditioner_env%lattice_dims, &
539 : preconditioner_env%lattice_local_cells, local_indices, &
540 4 : local_operator, local_coarse_inverse, corrected_projection_error, local_info)
541 4 : IF (local_info /= 0) THEN
542 0 : CPABORT("LATTICE_FFT local correction found a non-positive coarse operator")
543 : END IF
544 12 : preconditioner_env%lattice_state_local_indices(:, istate) = local_indices
545 44 : preconditioner_env%lattice_state_local_operator(:, :, istate) = local_operator
546 28 : preconditioner_env%lattice_state_local_coarse_inverse(:, :, istate) = local_coarse_inverse
547 4 : DEALLOCATE (local_coarse_inverse, local_indices, local_operator)
548 : END IF
549 : preconditioner_env%lattice_corrected_projection_error = &
550 4 : MAX(preconditioner_env%lattice_corrected_projection_error, corrected_projection_error)
551 : selected = lattice_fft_auto_select(preconditioner_env%lattice_dims, &
552 : preconditioner_env%lattice_block_size, nstate, nstate, &
553 4 : preconditioner_env%para_env%num_pe, projection_error, used_fft, info)
554 4 : IF (preconditioner_env%lattice_fft_mode == ot_lattice_fft_on) THEN
555 4 : IF (info /= 0) THEN
556 0 : CPABORT("LATTICE_FFT found a non-positive Fourier block for "//TRIM(operator_name))
557 : END IF
558 4 : IF (.NOT. used_fft) THEN
559 0 : CPABORT("LATTICE_FFT could not initialize its FFT for "//TRIM(operator_name))
560 : END IF
561 : ELSE
562 0 : all_selected = all_selected .AND. selected
563 : END IF
564 4 : IF (.NOT. ALLOCATED(preconditioner_env%lattice_state_inverse_k)) THEN
565 0 : ALLOCATE (preconditioner_env%lattice_state_inverse_k( &
566 12 : SIZE(inverse_k, 1), SIZE(inverse_k, 2), SIZE(inverse_k, 3), nstate))
567 : END IF
568 4 : IF (info == 0 .AND. used_fft) THEN
569 60 : preconditioner_env%lattice_state_inverse_k(:, :, :, istate) = inverse_k
570 : END IF
571 6 : DEALLOCATE (inverse_k)
572 : END DO
573 :
574 : preconditioner_env%lattice_fft_active = &
575 2 : preconditioner_env%lattice_fft_mode == ot_lattice_fft_on .OR. all_selected
576 2 : IF (.NOT. preconditioner_env%lattice_fft_active .AND. &
577 : ALLOCATED(preconditioner_env%lattice_state_inverse_k)) THEN
578 0 : DEALLOCATE (preconditioner_env%lattice_state_inverse_k)
579 : END IF
580 2 : IF (.NOT. preconditioner_env%lattice_fft_active) THEN
581 0 : DEALLOCATE (preconditioner_env%lattice_occ_vectors, &
582 0 : preconditioner_env%lattice_occ_dual)
583 0 : IF (ALLOCATED(preconditioner_env%lattice_state_local_indices)) THEN
584 : DEALLOCATE (preconditioner_env%lattice_state_local_indices, &
585 0 : preconditioner_env%lattice_state_local_operator, &
586 0 : preconditioner_env%lattice_state_local_coarse_inverse)
587 : END IF
588 : END IF
589 2 : IF (preconditioner_env%para_env%is_source()) THEN
590 1 : logger => cp_get_default_logger()
591 1 : output_unit = cp_logger_get_default_io_unit(logger)
592 1 : IF (output_unit > 0) THEN
593 : WRITE (output_unit, '(T2,A,A,A,3(I0,1X),A,I0,A,I0)') &
594 1 : "OT| Lattice FFT ", TRIM(operator_name), " cells: ", &
595 1 : preconditioner_env%lattice_dims, "AO block: ", &
596 2 : preconditioner_env%lattice_block_size, " states: ", nstate
597 : WRITE (output_unit, '(T2,A,ES12.4,A,L1)') &
598 1 : "OT| Maximum block-circulant projection error: ", &
599 1 : preconditioner_env%lattice_projection_error, &
600 2 : " selected: ", preconditioner_env%lattice_fft_active
601 1 : IF (preconditioner_env%lattice_local_cells > 0) THEN
602 : WRITE (output_unit, '(T2,A,I0,A,ES12.4)') &
603 1 : "OT| Balanced local correction cells/state: ", &
604 1 : preconditioner_env%lattice_local_cells, " maximum residual error: ", &
605 2 : preconditioner_env%lattice_corrected_projection_error
606 : END IF
607 : WRITE (output_unit, '(T2,A,F8.3,A,F8.3)') &
608 1 : "OT| Estimated lattice/dense work ratio: ", &
609 1 : preconditioner_env%lattice_cost_ratio, " storage ratio: ", &
610 2 : preconditioner_env%lattice_storage_ratio
611 : END IF
612 : END IF
613 :
614 2 : DEALLOCATE (dense_operator, dense_sv, dense_vectors, scaled_sv)
615 :
616 12 : END SUBROUTINE prepare_lattice_state_operator
617 :
618 : ! **************************************************************************************************
619 : !> \brief Build the overlap-inverse base of the bounded low-rank spectral preconditioner.
620 : !> \param preconditioner_env ...
621 : !> \param matrix_s ...
622 : ! **************************************************************************************************
623 6 : SUBROUTINE prepare_low_rank_base(preconditioner_env, matrix_s)
624 :
625 : TYPE(preconditioner_type) :: preconditioner_env
626 : TYPE(dbcsr_type), OPTIONAL, POINTER :: matrix_s
627 :
628 : INTEGER :: info
629 : TYPE(cp_fm_type) :: work
630 :
631 0 : CPASSERT(ASSOCIATED(preconditioner_env%fm))
632 6 : CPASSERT(.NOT. ASSOCIATED(preconditioner_env%base_fm))
633 :
634 6 : IF (preconditioner_env%lattice_fft_mode /= ot_lattice_fft_off) THEN
635 2 : IF (.NOT. PRESENT(matrix_s)) THEN
636 0 : CPABORT("LATTICE_FFT requires a non-unit overlap metric for FERMI_LOW_RANK")
637 : END IF
638 2 : CPASSERT(ASSOCIATED(matrix_s))
639 : CALL prepare_lattice_operator(preconditioner_env, matrix_s, "FERMI_LOW_RANK overlap", &
640 2 : SIZE(preconditioner_env%occ_evals))
641 2 : IF (preconditioner_env%lattice_fft_active) THEN
642 2 : preconditioner_env%low_rank_base = ot_low_rank_base_lattice_fft
643 2 : RETURN
644 : END IF
645 0 : preconditioner_env%low_rank_base = ot_low_rank_base_overlap
646 : END IF
647 :
648 4 : IF (preconditioner_env%low_rank_base /= ot_low_rank_base_overlap) THEN
649 0 : CPABORT("Unknown FERMI_LOW_RANK base operator")
650 : END IF
651 4 : ALLOCATE (preconditioner_env%base_fm)
652 : CALL cp_fm_create(preconditioner_env%base_fm, preconditioner_env%fm%matrix_struct, &
653 4 : name="FERMI_LOW_RANK overlap inverse")
654 :
655 4 : IF (PRESENT(matrix_s)) THEN
656 4 : CPASSERT(ASSOCIATED(matrix_s))
657 4 : CALL copy_dbcsr_to_fm(matrix_s, preconditioner_env%base_fm)
658 4 : CALL cp_fm_cholesky_decompose(preconditioner_env%base_fm, info_out=info)
659 4 : IF (info /= 0) CPABORT("FERMI_LOW_RANK requires a positive-definite overlap matrix")
660 4 : CALL cp_fm_cholesky_invert(preconditioner_env%base_fm)
661 4 : CALL cp_fm_create(work, preconditioner_env%base_fm%matrix_struct, name="low-rank base work")
662 4 : CALL cp_fm_uplo_to_full(preconditioner_env%base_fm, work)
663 4 : CALL cp_fm_release(work)
664 : ELSE
665 0 : CALL cp_fm_set_all(preconditioner_env%base_fm, alpha=0.0_dp, beta=1.0_dp)
666 : END IF
667 :
668 6 : END SUBROUTINE prepare_low_rank_base
669 :
670 : ! **************************************************************************************************
671 :
672 : ! creates a preconditioner for the system (H-energy_homo S)
673 : ! this preconditioner is (must be) symmetric positive definite.
674 : ! currently uses a atom-block-diagonal form
675 : ! each block will be ....
676 : ! might overwrite matrix_h, matrix_t
677 :
678 : ! **************************************************************************************************
679 : !> \brief ...
680 : !> \param preconditioner_env ...
681 : !> \param precon_type ...
682 : !> \param solver_type ...
683 : !> \param matrix_h ...
684 : !> \param matrix_s ...
685 : !> \param matrix_t ...
686 : !> \param mo_set ...
687 : !> \param energy_gap ...
688 : !> \param convert_precond_to_dbcsr ...
689 : !> \param chol_type ...
690 : !> \par History
691 : !> 09.2014 removed some unused or unfinished methods
692 : !> removed sparse preconditioners and the
693 : !> sparse approximate inverse at rev 14341 [Florian Schiffmann]
694 : ! **************************************************************************************************
695 9940 : SUBROUTINE make_preconditioner(preconditioner_env, precon_type, solver_type, matrix_h, matrix_s, &
696 : matrix_t, mo_set, energy_gap, convert_precond_to_dbcsr, chol_type)
697 :
698 : TYPE(preconditioner_type) :: preconditioner_env
699 : INTEGER, INTENT(IN) :: precon_type, solver_type
700 : TYPE(dbcsr_type), POINTER :: matrix_h
701 : TYPE(dbcsr_type), OPTIONAL, POINTER :: matrix_s, matrix_t
702 : TYPE(mo_set_type), INTENT(IN) :: mo_set
703 : REAL(KIND=dp) :: energy_gap
704 : LOGICAL, INTENT(IN), OPTIONAL :: convert_precond_to_dbcsr
705 : INTEGER, INTENT(IN), OPTIONAL :: chol_type
706 :
707 : CHARACTER(len=*), PARAMETER :: routineN = 'make_preconditioner'
708 :
709 : INTEGER :: handle, k, max_rank, my_solver_type, &
710 : nao, nhomo, output_unit
711 : LOGICAL :: my_convert_precond_to_dbcsr, &
712 : needs_full_spectrum, needs_homo, &
713 : use_mo_coeff_b
714 : REAL(KIND=dp) :: energy_homo
715 9940 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eigenvalues_ot
716 : TYPE(cp_fm_struct_type), POINTER :: fm_struct
717 : TYPE(cp_fm_type) :: mo_occ
718 : TYPE(cp_fm_type), POINTER :: mo_coeff
719 : TYPE(cp_logger_type), POINTER :: logger
720 : TYPE(dbcsr_type), POINTER :: mo_coeff_b
721 :
722 9940 : CALL timeset(routineN, handle)
723 :
724 9940 : CALL get_mo_set(mo_set=mo_set, mo_coeff=mo_coeff, mo_coeff_b=mo_coeff_b, homo=nhomo)
725 9940 : use_mo_coeff_b = mo_set%use_mo_coeff_b
726 9940 : CALL cp_fm_get_info(mo_coeff, ncol_global=k, nrow_global=nao)
727 :
728 : ! Starting some matrix mess, check where to store the result in preconditioner_env, fm or dbcsr_matrix
729 9940 : my_convert_precond_to_dbcsr = .FALSE.
730 9940 : IF (PRESENT(convert_precond_to_dbcsr)) my_convert_precond_to_dbcsr = convert_precond_to_dbcsr
731 9940 : IF (precon_type == ot_precond_fermi_low_rank .AND. &
732 : preconditioner_env%low_rank_base == ot_low_rank_base_lattice_fft) THEN
733 : ! The FFT base is an operator rather than a stored DBCSR matrix. Keep the
734 : ! complementary vectors in their distributed full-matrix representation.
735 0 : my_convert_precond_to_dbcsr = .FALSE.
736 : END IF
737 :
738 : ! Thanks to the mess with the matrices we need to make sure in this case that the
739 : ! Previous inverse is properly stored as a sparse matrix, fm gets deallocated here
740 : ! if it wasn't anyway
741 9940 : IF (preconditioner_env%solver == ot_precond_solver_update) THEN
742 4 : CALL transfer_fm_to_dbcsr(preconditioner_env%fm, preconditioner_env%dbcsr_matrix, matrix_h)
743 : END IF
744 :
745 9940 : needs_full_spectrum = .FALSE.
746 9940 : needs_homo = .FALSE.
747 :
748 3784 : SELECT CASE (precon_type)
749 : CASE (ot_precond_full_all, ot_precond_fermi_low_rank, ot_precond_full_all_covariant)
750 3784 : needs_full_spectrum = .TRUE.
751 : ! both of them need the coefficients as fm's, more matrix mess
752 3784 : IF (use_mo_coeff_b) THEN
753 3472 : CALL copy_dbcsr_to_fm(mo_coeff_b, mo_coeff)
754 : END IF
755 : CASE (ot_precond_full_single)
756 38 : needs_homo = .TRUE.
757 : ! XXXX to be removed if homo estimate only is implemented
758 38 : needs_full_spectrum = .TRUE.
759 : CASE (ot_precond_full_kinetic, ot_precond_s_inverse, ot_precond_full_single_inverse)
760 : ! these should be happy without an estimate for the homo energy
761 : ! preconditioning can not depend on an absolute eigenvalue, only on eigenvalue differences
762 : CASE DEFAULT
763 9940 : CPABORT("The preconditioner is unknown ...")
764 : END SELECT
765 :
766 29532 : ALLOCATE (eigenvalues_ot(k))
767 9940 : energy_homo = 0.0_dp
768 9940 : IF (needs_full_spectrum) THEN
769 : ! XXXXXXXXXXXXXXXX do not touch the initial MOs, could be harmful for either
770 : ! the case of non-equivalent MOs but also for the derivate
771 : ! we could already have all eigenvalues e.g. full_all and we could skip this
772 : ! to be optimised later.
773 : ! one flaw is that not all SCF methods (i.e. that go over mo_derivs directly)
774 : ! have a 'valid' matrix_h... (we even don't know what evals are in that case)
775 3822 : IF (use_mo_coeff_b) THEN
776 : CALL calculate_subspace_eigenvalues(mo_coeff_b, matrix_h, &
777 : eigenvalues_ot, do_rotation=.FALSE., &
778 : para_env=mo_coeff%matrix_struct%para_env, &
779 3502 : blacs_env=mo_coeff%matrix_struct%context)
780 : ELSE
781 : CALL calculate_subspace_eigenvalues(mo_coeff, matrix_h, &
782 320 : eigenvalues_ot, do_rotation=.FALSE.)
783 : END IF
784 3822 : IF (k > 0) THEN
785 3712 : CPASSERT(nhomo > 0 .AND. nhomo <= k)
786 3712 : energy_homo = eigenvalues_ot(nhomo)
787 : END IF
788 : ELSE
789 6118 : IF (needs_homo) THEN
790 0 : CPABORT("Not yet implemented")
791 : END IF
792 : END IF
793 :
794 : ! After all bits and pieces of checking and initialization, here comes the
795 : ! part where the preconditioner matrix gets created and solved.
796 : ! This will give the matrices for later use
797 9940 : my_solver_type = solver_type
798 9940 : preconditioner_env%in_use = precon_type
799 9940 : preconditioner_env%cholesky_use = cholesky_reduce
800 9940 : IF (PRESENT(chol_type)) preconditioner_env%cholesky_use = chol_type
801 : preconditioner_env%in_use = precon_type
802 9940 : preconditioner_env%lattice_fft_active = .FALSE.
803 9940 : preconditioner_env%lattice_corrected_projection_error = 0.0_dp
804 9940 : IF (ALLOCATED(preconditioner_env%lattice_inverse_k)) THEN
805 0 : DEALLOCATE (preconditioner_env%lattice_inverse_k)
806 : END IF
807 9940 : IF (ALLOCATED(preconditioner_env%lattice_state_inverse_k)) THEN
808 0 : DEALLOCATE (preconditioner_env%lattice_state_inverse_k)
809 : END IF
810 9940 : IF (ALLOCATED(preconditioner_env%lattice_occ_vectors)) THEN
811 0 : DEALLOCATE (preconditioner_env%lattice_occ_vectors)
812 : END IF
813 9940 : IF (ALLOCATED(preconditioner_env%lattice_occ_dual)) THEN
814 0 : DEALLOCATE (preconditioner_env%lattice_occ_dual)
815 : END IF
816 9940 : IF (ALLOCATED(preconditioner_env%lattice_local_indices)) THEN
817 0 : DEALLOCATE (preconditioner_env%lattice_local_indices)
818 : END IF
819 9940 : IF (ALLOCATED(preconditioner_env%lattice_local_operator)) THEN
820 0 : DEALLOCATE (preconditioner_env%lattice_local_operator)
821 : END IF
822 9940 : IF (ALLOCATED(preconditioner_env%lattice_local_coarse_inverse)) THEN
823 0 : DEALLOCATE (preconditioner_env%lattice_local_coarse_inverse)
824 : END IF
825 9940 : IF (ALLOCATED(preconditioner_env%lattice_state_local_indices)) THEN
826 0 : DEALLOCATE (preconditioner_env%lattice_state_local_indices)
827 : END IF
828 9940 : IF (ALLOCATED(preconditioner_env%lattice_state_local_operator)) THEN
829 0 : DEALLOCATE (preconditioner_env%lattice_state_local_operator)
830 : END IF
831 9940 : IF (ALLOCATED(preconditioner_env%lattice_state_local_coarse_inverse)) THEN
832 0 : DEALLOCATE (preconditioner_env%lattice_state_local_coarse_inverse)
833 : END IF
834 9940 : IF (nhomo == k) THEN
835 : CALL make_preconditioner_matrix(preconditioner_env, matrix_h, matrix_s, matrix_t, mo_coeff, &
836 9852 : energy_homo, eigenvalues_ot, energy_gap, my_solver_type)
837 : ELSE
838 : CALL cp_fm_struct_create(fm_struct, nrow_global=nao, ncol_global=nhomo, &
839 : context=preconditioner_env%ctxt, &
840 88 : para_env=preconditioner_env%para_env)
841 88 : CALL cp_fm_create(mo_occ, fm_struct)
842 88 : CALL cp_fm_to_fm(mo_coeff, mo_occ, nhomo)
843 88 : CALL cp_fm_struct_release(fm_struct)
844 : !
845 : CALL make_preconditioner_matrix(preconditioner_env, matrix_h, matrix_s, matrix_t, mo_occ, &
846 88 : energy_homo, eigenvalues_ot(1:nhomo), energy_gap, my_solver_type)
847 : !
848 88 : CALL cp_fm_release(mo_occ)
849 : END IF
850 :
851 9940 : IF (preconditioner_env%lattice_fft_mode /= ot_lattice_fft_off) THEN
852 8 : SELECT CASE (precon_type)
853 : CASE (ot_precond_full_kinetic)
854 8 : CPASSERT(ASSOCIATED(preconditioner_env%sparse_matrix))
855 : CALL prepare_lattice_operator(preconditioner_env, preconditioner_env%sparse_matrix, &
856 8 : "FULL_KINETIC", nhomo)
857 : CASE (ot_precond_full_single_inverse)
858 2 : CPASSERT(ASSOCIATED(preconditioner_env%sparse_matrix))
859 : CALL prepare_lattice_operator(preconditioner_env, preconditioner_env%sparse_matrix, &
860 2 : "FULL_SINGLE_INVERSE", nhomo)
861 : CASE (ot_precond_s_inverse)
862 4 : CPASSERT(ASSOCIATED(preconditioner_env%sparse_matrix))
863 : CALL prepare_lattice_operator(preconditioner_env, preconditioner_env%sparse_matrix, &
864 4 : "FULL_S_INVERSE", nhomo)
865 : CASE (ot_precond_full_all_covariant)
866 2 : CPASSERT(PRESENT(matrix_s))
867 2 : CPASSERT(ASSOCIATED(matrix_s))
868 : CALL prepare_lattice_state_operator(preconditioner_env, matrix_s, &
869 20 : "FULL_ALL_COVARIANT")
870 : END SELECT
871 : END IF
872 :
873 9940 : IF (.NOT. preconditioner_env%lattice_fft_active) THEN
874 9928 : CALL solve_preconditioner(my_solver_type, preconditioner_env, matrix_s, matrix_h)
875 : ELSE
876 12 : preconditioner_env%solver = ot_precond_solver_inv_chol
877 : END IF
878 :
879 9940 : IF (precon_type == ot_precond_fermi_low_rank) THEN
880 6 : preconditioner_env%spectral_reference = energy_homo
881 : preconditioner_env%spectral_window = &
882 6 : MAX(fermi_low_rank_min_window, fermi_low_rank_gap_factor*preconditioner_env%energy_gap)
883 6 : preconditioner_env%spectral_base_scale = 1.0_dp/preconditioner_env%spectral_window
884 6 : max_rank = preconditioner_env%spectral_max_rank
885 6 : IF (max_rank == -1) max_rank = SIZE(preconditioner_env%full_evals) - nhomo
886 : preconditioner_env%spectral_rank = low_rank_select_rank( &
887 : preconditioner_env%full_evals, nhomo, max_rank, &
888 6 : fermi_low_rank_degeneracy_tolerance)
889 6 : CALL prepare_low_rank_base(preconditioner_env, matrix_s)
890 6 : IF (preconditioner_env%lattice_fft_active) THEN
891 : ! With an approximate base, retain only non-negative spectral updates.
892 : ! This preserves positive definiteness independently of the projection error.
893 26 : DO WHILE (preconditioner_env%spectral_rank > 0)
894 26 : IF (low_rank_inverse_weight( &
895 : preconditioner_env%full_evals(nhomo + preconditioner_env%spectral_rank), &
896 : preconditioner_env%spectral_reference, preconditioner_env%energy_gap) >= &
897 : preconditioner_env%spectral_base_scale) EXIT
898 26 : preconditioner_env%spectral_rank = preconditioner_env%spectral_rank - 1
899 : END DO
900 : END IF
901 6 : IF (preconditioner_env%para_env%is_source()) THEN
902 3 : logger => cp_get_default_logger()
903 3 : output_unit = cp_logger_get_default_io_unit(logger)
904 3 : IF (output_unit > 0) THEN
905 : WRITE (output_unit, '(T2,A,F14.8,A,I0)') &
906 3 : "OT| FERMI_LOW_RANK common reference [hartree]: ", &
907 3 : preconditioner_env%spectral_reference, " retained rank: ", &
908 6 : preconditioner_env%spectral_rank
909 : END IF
910 : END IF
911 : END IF
912 :
913 : ! Here comes more matrix mess, make sure to output the correct matrix format,
914 : ! A bit pointless to convert the cholesky factorized version as it doesn't work in
915 : ! dbcsr form and will crash later,...
916 9940 : IF (.NOT. preconditioner_env%lattice_fft_active) THEN
917 9926 : IF (my_convert_precond_to_dbcsr) THEN
918 8196 : CALL transfer_fm_to_dbcsr(preconditioner_env%fm, preconditioner_env%dbcsr_matrix, matrix_h)
919 8196 : IF (precon_type == ot_precond_fermi_low_rank .AND. &
920 : preconditioner_env%low_rank_base == ot_low_rank_base_overlap) THEN
921 : CALL transfer_fm_to_dbcsr(preconditioner_env%base_fm, &
922 4 : preconditioner_env%base_dbcsr_matrix, matrix_h)
923 : END IF
924 : ELSE
925 : CALL transfer_dbcsr_to_fm(preconditioner_env%dbcsr_matrix, preconditioner_env%fm, &
926 1730 : preconditioner_env%para_env, preconditioner_env%ctxt)
927 : END IF
928 : END IF
929 :
930 9940 : DEALLOCATE (eigenvalues_ot)
931 :
932 9940 : CALL timestop(handle)
933 :
934 9940 : END SUBROUTINE make_preconditioner
935 :
936 : ! **************************************************************************************************
937 : !> \brief Construct FULL_ALL directly from one complex H(k), S(k), and C(k) channel.
938 : !> \param preconditioner_env preconditioner storage
939 : !> \param matrix_c_re real part of the active OT reference orbitals
940 : !> \param matrix_c_im imaginary part of the active OT reference orbitals
941 : !> \param matrix_h_re real part of H(k)
942 : !> \param matrix_h_im imaginary part of H(k)
943 : !> \param matrix_s_re real part of S(k)
944 : !> \param matrix_s_im imaginary part of S(k)
945 : !> \param mo_set ...
946 : !> \param energy_gap denominator floor
947 : !> \param solver_type inverse application selected in the OT input
948 : ! **************************************************************************************************
949 565 : SUBROUTINE make_preconditioner_complex_full_all(preconditioner_env, &
950 : matrix_c_re, matrix_c_im, &
951 : matrix_h_re, matrix_h_im, &
952 : matrix_s_re, matrix_s_im, &
953 : mo_set, energy_gap, solver_type)
954 :
955 : TYPE(preconditioner_type) :: preconditioner_env
956 : TYPE(dbcsr_type), POINTER :: matrix_c_re, matrix_c_im, matrix_h_re, &
957 : matrix_h_im, matrix_s_re, matrix_s_im
958 : TYPE(mo_set_type), INTENT(IN) :: mo_set
959 : REAL(KIND=dp), INTENT(IN) :: energy_gap
960 : INTEGER, INTENT(IN) :: solver_type
961 :
962 : CHARACTER(len=*), PARAMETER :: routineN = 'make_preconditioner_complex_full_all'
963 :
964 : INTEGER :: handle, k
965 113 : REAL(KIND=dp), DIMENSION(:), POINTER :: eigenvalues
966 : TYPE(cp_cfm_type) :: matrix_c, matrix_h, matrix_s
967 :
968 113 : CALL timeset(routineN, handle)
969 :
970 113 : NULLIFY (eigenvalues)
971 113 : CPASSERT(ASSOCIATED(matrix_c_re))
972 113 : CPASSERT(ASSOCIATED(matrix_c_im))
973 113 : CPASSERT(ASSOCIATED(matrix_h_re))
974 113 : CPASSERT(ASSOCIATED(matrix_h_im))
975 113 : CPASSERT(ASSOCIATED(matrix_s_re))
976 113 : CPASSERT(ASSOCIATED(matrix_s_im))
977 113 : IF (solver_type /= ot_precond_solver_default) THEN
978 0 : CPABORT('Complex FULL_ALL supports only PRECOND_SOLVER DEFAULT')
979 : END IF
980 113 : CALL dbcsr_get_info(matrix_c_re, nfullcols_total=k)
981 113 : CPASSERT(k > 0)
982 113 : CALL get_mo_set(mo_set, eigenvalues=eigenvalues)
983 113 : CPASSERT(ASSOCIATED(eigenvalues))
984 113 : CPASSERT(SIZE(eigenvalues) >= k)
985 :
986 : CALL dbcsr_pair_to_cfm(matrix_c_re, matrix_c_im, preconditioner_env, &
987 113 : 'complex FULL_ALL C', matrix_c)
988 : CALL dbcsr_pair_to_cfm(matrix_h_re, matrix_h_im, preconditioner_env, &
989 113 : 'complex FULL_ALL H', matrix_h)
990 : CALL dbcsr_pair_to_cfm(matrix_s_re, matrix_s_im, preconditioner_env, &
991 113 : 'complex FULL_ALL S', matrix_s)
992 : CALL make_complex_full_all(preconditioner_env, matrix_c, matrix_h, matrix_s, &
993 113 : eigenvalues(1:k), energy_gap)
994 :
995 113 : CALL cp_cfm_release(matrix_c)
996 113 : CALL cp_cfm_release(matrix_s)
997 113 : CALL cp_cfm_release(matrix_h)
998 :
999 113 : CALL timestop(handle)
1000 :
1001 113 : END SUBROUTINE make_preconditioner_complex_full_all
1002 :
1003 : ! **************************************************************************************************
1004 : !> \brief Construct a complex FULL_SINGLE preconditioner from H(k) and S(k).
1005 : !> \param preconditioner_env preconditioner storage
1006 : !> \param matrix_h_re real part of H(k)
1007 : !> \param matrix_h_im imaginary part of H(k)
1008 : !> \param matrix_s_re real part of S(k)
1009 : !> \param matrix_s_im imaginary part of S(k)
1010 : !> \param mo_set orbital energy labels defining the occupied edge
1011 : !> \param energy_gap denominator floor
1012 : !> \param solver_type requested inverse solver
1013 : ! **************************************************************************************************
1014 104 : SUBROUTINE make_preconditioner_complex_full_single(preconditioner_env, &
1015 : matrix_h_re, matrix_h_im, &
1016 : matrix_s_re, matrix_s_im, &
1017 : mo_set, energy_gap, solver_type)
1018 :
1019 : TYPE(preconditioner_type) :: preconditioner_env
1020 : TYPE(dbcsr_type), POINTER :: matrix_h_re, matrix_h_im, matrix_s_re, &
1021 : matrix_s_im
1022 : TYPE(mo_set_type), INTENT(IN) :: mo_set
1023 : REAL(KIND=dp), INTENT(IN) :: energy_gap
1024 : INTEGER, INTENT(IN) :: solver_type
1025 :
1026 : INTEGER :: homo
1027 26 : REAL(KIND=dp), DIMENSION(:), POINTER :: eigenvalues
1028 : TYPE(cp_cfm_type) :: matrix_h, matrix_s
1029 :
1030 26 : NULLIFY (eigenvalues)
1031 0 : CPASSERT(ASSOCIATED(matrix_h_re))
1032 26 : CPASSERT(ASSOCIATED(matrix_h_im))
1033 26 : CPASSERT(ASSOCIATED(matrix_s_re))
1034 26 : CPASSERT(ASSOCIATED(matrix_s_im))
1035 26 : IF (solver_type /= ot_precond_solver_default) THEN
1036 0 : CPABORT('Complex FULL_SINGLE supports only PRECOND_SOLVER DEFAULT')
1037 : END IF
1038 26 : CALL get_mo_set(mo_set, homo=homo, eigenvalues=eigenvalues)
1039 26 : CPASSERT(ASSOCIATED(eigenvalues))
1040 26 : CPASSERT(homo > 0 .AND. homo <= SIZE(eigenvalues))
1041 :
1042 : CALL dbcsr_pair_to_cfm(matrix_h_re, matrix_h_im, preconditioner_env, &
1043 26 : 'complex FULL_SINGLE H', matrix_h)
1044 : CALL dbcsr_pair_to_cfm(matrix_s_re, matrix_s_im, preconditioner_env, &
1045 26 : 'complex FULL_SINGLE S', matrix_s)
1046 : CALL make_complex_full_single(preconditioner_env, matrix_h, matrix_s, &
1047 26 : eigenvalues(homo), energy_gap)
1048 26 : CALL cp_cfm_release(matrix_s)
1049 26 : CALL cp_cfm_release(matrix_h)
1050 :
1051 26 : END SUBROUTINE make_preconditioner_complex_full_single
1052 :
1053 : ! **************************************************************************************************
1054 : !> \brief Construct FERMI_LOW_RANK directly for one complex k-point channel.
1055 : !> \param preconditioner_env preconditioner storage
1056 : !> \param matrix_c_re real part of the active OT reference orbitals
1057 : !> \param matrix_c_im imaginary part of the active OT reference orbitals
1058 : !> \param matrix_h_re real part of H(k)
1059 : !> \param matrix_h_im imaginary part of H(k)
1060 : !> \param matrix_s_re real part of S(k)
1061 : !> \param matrix_s_im imaginary part of S(k)
1062 : !> \param energy_gap denominator floor
1063 : !> \param max_rank maximum complementary rank, or -1 for the full complement
1064 : !> \param solver_type inverse application selected in the OT input
1065 : ! **************************************************************************************************
1066 184 : SUBROUTINE make_preconditioner_complex_fermi_low_rank( &
1067 : preconditioner_env, matrix_c_re, matrix_c_im, matrix_h_re, matrix_h_im, &
1068 : matrix_s_re, matrix_s_im, energy_gap, max_rank, solver_type)
1069 :
1070 : TYPE(preconditioner_type) :: preconditioner_env
1071 : TYPE(dbcsr_type), POINTER :: matrix_c_re, matrix_c_im, matrix_h_re, &
1072 : matrix_h_im, matrix_s_re, matrix_s_im
1073 : REAL(KIND=dp), INTENT(IN) :: energy_gap
1074 : INTEGER, INTENT(IN) :: max_rank, solver_type
1075 :
1076 : TYPE(cp_cfm_type) :: matrix_c, matrix_h, matrix_s
1077 :
1078 46 : IF (solver_type /= ot_precond_solver_default) THEN
1079 0 : CPABORT('Complex FERMI_LOW_RANK supports only PRECOND_SOLVER DEFAULT')
1080 : END IF
1081 :
1082 : CALL dbcsr_pair_to_cfm(matrix_c_re, matrix_c_im, preconditioner_env, &
1083 46 : 'complex FERMI_LOW_RANK C', matrix_c)
1084 : CALL dbcsr_pair_to_cfm(matrix_h_re, matrix_h_im, preconditioner_env, &
1085 46 : 'complex FERMI_LOW_RANK H', matrix_h)
1086 : CALL dbcsr_pair_to_cfm(matrix_s_re, matrix_s_im, preconditioner_env, &
1087 46 : 'complex FERMI_LOW_RANK S', matrix_s)
1088 : CALL make_complex_fermi_low_rank( &
1089 : preconditioner_env, matrix_c, matrix_h, matrix_s, energy_gap, &
1090 : max_rank, MAX(fermi_low_rank_min_window, fermi_low_rank_gap_factor*energy_gap), &
1091 46 : fermi_low_rank_degeneracy_tolerance)
1092 :
1093 46 : CALL cp_cfm_release(matrix_s)
1094 46 : CALL cp_cfm_release(matrix_h)
1095 46 : CALL cp_cfm_release(matrix_c)
1096 :
1097 46 : END SUBROUTINE make_preconditioner_complex_fermi_low_rank
1098 :
1099 : ! **************************************************************************************************
1100 : !> \brief Construct rotation-covariant FULL_ALL for one complex k-point channel.
1101 : !> \param preconditioner_env preconditioner storage
1102 : !> \param matrix_c_re real part of the active OT reference orbitals
1103 : !> \param matrix_c_im imaginary part of the active OT reference orbitals
1104 : !> \param matrix_h_re real part of H(k)
1105 : !> \param matrix_h_im imaginary part of H(k)
1106 : !> \param matrix_s_re real part of S(k)
1107 : !> \param matrix_s_im imaginary part of S(k)
1108 : !> \param energy_gap denominator floor
1109 : !> \param solver_type inverse application selected in the OT input
1110 : !> \param occupation_signature occupations that define the admissible private gauge
1111 : ! **************************************************************************************************
1112 208 : SUBROUTINE make_preconditioner_complex_full_all_covariant( &
1113 : preconditioner_env, matrix_c_re, matrix_c_im, matrix_h_re, matrix_h_im, &
1114 52 : matrix_s_re, matrix_s_im, energy_gap, solver_type, occupation_signature)
1115 :
1116 : TYPE(preconditioner_type) :: preconditioner_env
1117 : TYPE(dbcsr_type), POINTER :: matrix_c_re, matrix_c_im, matrix_h_re, &
1118 : matrix_h_im, matrix_s_re, matrix_s_im
1119 : REAL(KIND=dp), INTENT(IN) :: energy_gap
1120 : INTEGER, INTENT(IN) :: solver_type
1121 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN), &
1122 : OPTIONAL :: occupation_signature
1123 :
1124 : TYPE(cp_cfm_type) :: matrix_c, matrix_h, matrix_s
1125 :
1126 52 : IF (solver_type /= ot_precond_solver_default) THEN
1127 0 : CPABORT('Complex FULL_ALL_COVARIANT supports only PRECOND_SOLVER DEFAULT')
1128 : END IF
1129 : CALL dbcsr_pair_to_cfm(matrix_c_re, matrix_c_im, preconditioner_env, &
1130 52 : 'complex FULL_ALL_COVARIANT C', matrix_c)
1131 : CALL dbcsr_pair_to_cfm(matrix_h_re, matrix_h_im, preconditioner_env, &
1132 52 : 'complex FULL_ALL_COVARIANT H', matrix_h)
1133 : CALL dbcsr_pair_to_cfm(matrix_s_re, matrix_s_im, preconditioner_env, &
1134 52 : 'complex FULL_ALL_COVARIANT S', matrix_s)
1135 : CALL make_complex_full_all_covariant(preconditioner_env, matrix_c, matrix_h, matrix_s, &
1136 96 : energy_gap, occupation_signature)
1137 :
1138 52 : CALL cp_cfm_release(matrix_s)
1139 52 : CALL cp_cfm_release(matrix_h)
1140 52 : CALL cp_cfm_release(matrix_c)
1141 :
1142 52 : END SUBROUTINE make_preconditioner_complex_full_all_covariant
1143 :
1144 : ! **************************************************************************************************
1145 : !> \brief Construct a complex FULL_SINGLE_INVERSE preconditioner without discarding Im(H,S,C).
1146 : !> \param preconditioner_env preconditioner storage
1147 : !> \param matrix_c_re real part of the occupied reference orbitals
1148 : !> \param matrix_c_im imaginary part of the occupied reference orbitals
1149 : !> \param matrix_h_re real part of H(k)
1150 : !> \param matrix_h_im imaginary part of H(k)
1151 : !> \param matrix_s_re real part of S(k)
1152 : !> \param matrix_s_im imaginary part of S(k)
1153 : !> \param energy_gap lower spectral bound
1154 : !> \param solver_type requested inverse solver
1155 : ! **************************************************************************************************
1156 484 : SUBROUTINE make_preconditioner_complex_full_single_inverse(preconditioner_env, &
1157 : matrix_c_re, matrix_c_im, &
1158 : matrix_h_re, matrix_h_im, &
1159 : matrix_s_re, matrix_s_im, &
1160 : energy_gap, solver_type)
1161 :
1162 : TYPE(preconditioner_type) :: preconditioner_env
1163 : TYPE(dbcsr_type), POINTER :: matrix_c_re, matrix_c_im, matrix_h_re, &
1164 : matrix_h_im, matrix_s_re, matrix_s_im
1165 : REAL(KIND=dp), INTENT(IN) :: energy_gap
1166 : INTEGER, INTENT(IN) :: solver_type
1167 :
1168 : INTEGER :: k, n
1169 : TYPE(cp_cfm_type) :: matrix_c, matrix_h, matrix_s
1170 : TYPE(cp_fm_struct_type), POINTER :: fm_struct_mo, fm_struct_square
1171 : TYPE(cp_fm_type) :: matrix_c_im_fm, matrix_c_re_fm, &
1172 : matrix_h_im_fm, matrix_h_re_fm, &
1173 : matrix_s_im_fm, matrix_s_re_fm
1174 :
1175 44 : NULLIFY (fm_struct_mo, fm_struct_square)
1176 0 : CPASSERT(ASSOCIATED(matrix_c_re))
1177 44 : CPASSERT(ASSOCIATED(matrix_c_im))
1178 44 : CPASSERT(ASSOCIATED(matrix_h_re))
1179 44 : CPASSERT(ASSOCIATED(matrix_h_im))
1180 44 : CPASSERT(ASSOCIATED(matrix_s_re))
1181 44 : CPASSERT(ASSOCIATED(matrix_s_im))
1182 44 : IF (solver_type /= ot_precond_solver_default .AND. &
1183 : solver_type /= ot_precond_solver_inv_chol) THEN
1184 0 : CPABORT("Complex FULL_SINGLE_INVERSE supports only PRECOND_SOLVER DEFAULT/INVERSE_CHOLESKY")
1185 : END IF
1186 44 : CALL dbcsr_get_info(matrix_c_re, nfullrows_total=n, nfullcols_total=k)
1187 44 : CPASSERT(n > 0 .AND. k > 0 .AND. k <= n)
1188 :
1189 : CALL cp_fm_struct_create(fm_struct_mo, nrow_global=n, ncol_global=k, &
1190 : context=preconditioner_env%ctxt, &
1191 44 : para_env=preconditioner_env%para_env)
1192 : CALL cp_fm_struct_create(fm_struct_square, nrow_global=n, ncol_global=n, &
1193 : context=preconditioner_env%ctxt, &
1194 44 : para_env=preconditioner_env%para_env)
1195 44 : CALL cp_fm_create(matrix_c_re_fm, fm_struct_mo, name='complex preconditioner C real')
1196 44 : CALL cp_fm_create(matrix_c_im_fm, fm_struct_mo, name='complex preconditioner C imaginary')
1197 44 : CALL cp_fm_create(matrix_h_re_fm, fm_struct_square, name='complex preconditioner H real')
1198 44 : CALL cp_fm_create(matrix_h_im_fm, fm_struct_square, name='complex preconditioner H imaginary')
1199 44 : CALL cp_fm_create(matrix_s_re_fm, fm_struct_square, name='complex preconditioner S real')
1200 44 : CALL cp_fm_create(matrix_s_im_fm, fm_struct_square, name='complex preconditioner S imaginary')
1201 44 : CALL cp_fm_struct_release(fm_struct_mo)
1202 44 : CALL cp_fm_struct_release(fm_struct_square)
1203 44 : CALL copy_dbcsr_to_fm(matrix_c_re, matrix_c_re_fm)
1204 44 : CALL copy_dbcsr_to_fm(matrix_c_im, matrix_c_im_fm)
1205 44 : CALL copy_dbcsr_to_fm(matrix_h_re, matrix_h_re_fm)
1206 44 : CALL copy_dbcsr_to_fm(matrix_h_im, matrix_h_im_fm)
1207 44 : CALL copy_dbcsr_to_fm(matrix_s_re, matrix_s_re_fm)
1208 44 : CALL copy_dbcsr_to_fm(matrix_s_im, matrix_s_im_fm)
1209 :
1210 44 : CALL cp_cfm_create(matrix_c, matrix_c_re_fm%matrix_struct, name='complex preconditioner C')
1211 44 : CALL cp_cfm_create(matrix_h, matrix_h_re_fm%matrix_struct, name='complex preconditioner H')
1212 44 : CALL cp_cfm_create(matrix_s, matrix_s_re_fm%matrix_struct, name='complex preconditioner S')
1213 44 : CALL cp_fm_to_cfm(matrix_c_re_fm, matrix_c_im_fm, matrix_c)
1214 44 : CALL cp_fm_to_cfm(matrix_h_re_fm, matrix_h_im_fm, matrix_h)
1215 44 : CALL cp_fm_to_cfm(matrix_s_re_fm, matrix_s_im_fm, matrix_s)
1216 : CALL make_complex_full_single_inverse(preconditioner_env, matrix_c, matrix_h, matrix_s, &
1217 44 : energy_gap)
1218 44 : preconditioner_env%solver = solver_type
1219 :
1220 44 : CALL cp_cfm_release(matrix_s)
1221 44 : CALL cp_cfm_release(matrix_h)
1222 44 : CALL cp_cfm_release(matrix_c)
1223 44 : CALL cp_fm_release(matrix_s_im_fm)
1224 44 : CALL cp_fm_release(matrix_s_re_fm)
1225 44 : CALL cp_fm_release(matrix_h_im_fm)
1226 44 : CALL cp_fm_release(matrix_h_re_fm)
1227 44 : CALL cp_fm_release(matrix_c_im_fm)
1228 44 : CALL cp_fm_release(matrix_c_re_fm)
1229 :
1230 44 : END SUBROUTINE make_preconditioner_complex_full_single_inverse
1231 :
1232 : ! **************************************************************************************************
1233 : !> \brief Construct a complex FULL_S_INVERSE preconditioner.
1234 : !> \param preconditioner_env preconditioner storage
1235 : !> \param matrix_s_re real part of S(k)
1236 : !> \param matrix_s_im imaginary part of S(k)
1237 : !> \param solver_type requested inverse solver
1238 : ! **************************************************************************************************
1239 164 : SUBROUTINE make_preconditioner_complex_full_s_inverse(preconditioner_env, &
1240 : matrix_s_re, matrix_s_im, solver_type)
1241 :
1242 : TYPE(preconditioner_type) :: preconditioner_env
1243 : TYPE(dbcsr_type), POINTER :: matrix_s_re, matrix_s_im
1244 : INTEGER, INTENT(IN) :: solver_type
1245 :
1246 : TYPE(cp_cfm_type) :: matrix_s
1247 :
1248 82 : IF (solver_type /= ot_precond_solver_default .AND. &
1249 : solver_type /= ot_precond_solver_inv_chol) THEN
1250 0 : CPABORT("Complex FULL_S_INVERSE supports only PRECOND_SOLVER DEFAULT/INVERSE_CHOLESKY")
1251 : END IF
1252 : CALL dbcsr_pair_to_cfm(matrix_s_re, matrix_s_im, preconditioner_env, &
1253 82 : 'complex preconditioner S', matrix_s)
1254 82 : CALL make_complex_full_s_inverse(preconditioner_env, matrix_s)
1255 82 : preconditioner_env%solver = solver_type
1256 82 : CALL cp_cfm_release(matrix_s)
1257 :
1258 82 : END SUBROUTINE make_preconditioner_complex_full_s_inverse
1259 :
1260 : ! **************************************************************************************************
1261 : !> \brief Construct a complex FULL_KINETIC preconditioner.
1262 : !> \param preconditioner_env preconditioner storage
1263 : !> \param matrix_t_re real part of T(k)
1264 : !> \param matrix_t_im imaginary part of T(k)
1265 : !> \param matrix_s_re real part of S(k)
1266 : !> \param matrix_s_im imaginary part of S(k)
1267 : !> \param energy_gap non-negative overlap shift
1268 : !> \param solver_type requested inverse solver
1269 : ! **************************************************************************************************
1270 360 : SUBROUTINE make_preconditioner_complex_full_kinetic(preconditioner_env, &
1271 : matrix_t_re, matrix_t_im, &
1272 : matrix_s_re, matrix_s_im, &
1273 : energy_gap, solver_type)
1274 :
1275 : TYPE(preconditioner_type) :: preconditioner_env
1276 : TYPE(dbcsr_type), POINTER :: matrix_t_re, matrix_t_im, matrix_s_re, &
1277 : matrix_s_im
1278 : REAL(KIND=dp), INTENT(IN) :: energy_gap
1279 : INTEGER, INTENT(IN) :: solver_type
1280 :
1281 : TYPE(cp_cfm_type) :: matrix_s, matrix_t
1282 :
1283 120 : IF (solver_type /= ot_precond_solver_default .AND. &
1284 : solver_type /= ot_precond_solver_inv_chol) THEN
1285 0 : CPABORT("Complex FULL_KINETIC supports only PRECOND_SOLVER DEFAULT/INVERSE_CHOLESKY")
1286 : END IF
1287 : CALL dbcsr_pair_to_cfm(matrix_t_re, matrix_t_im, preconditioner_env, &
1288 120 : 'complex preconditioner T', matrix_t)
1289 : CALL dbcsr_pair_to_cfm(matrix_s_re, matrix_s_im, preconditioner_env, &
1290 120 : 'complex preconditioner S', matrix_s)
1291 120 : CALL make_complex_full_kinetic(preconditioner_env, matrix_t, matrix_s, energy_gap)
1292 120 : preconditioner_env%solver = solver_type
1293 120 : CALL cp_cfm_release(matrix_s)
1294 120 : CALL cp_cfm_release(matrix_t)
1295 :
1296 120 : END SUBROUTINE make_preconditioner_complex_full_kinetic
1297 :
1298 : ! **************************************************************************************************
1299 : !> \brief Copy a real/imaginary DBCSR pair to one distributed complex full matrix.
1300 : !> \param matrix_re real matrix component
1301 : !> \param matrix_im imaginary matrix component
1302 : !> \param preconditioner_env source of the distribution context
1303 : !> \param matrix_name matrix label
1304 : !> \param matrix complex output matrix
1305 : ! **************************************************************************************************
1306 3357 : SUBROUTINE dbcsr_pair_to_cfm(matrix_re, matrix_im, preconditioner_env, matrix_name, matrix)
1307 :
1308 : TYPE(dbcsr_type), POINTER :: matrix_re, matrix_im
1309 : TYPE(preconditioner_type) :: preconditioner_env
1310 : CHARACTER(LEN=*), INTENT(IN) :: matrix_name
1311 : TYPE(cp_cfm_type), INTENT(OUT) :: matrix
1312 :
1313 : INTEGER :: ncol, nrow
1314 : TYPE(cp_fm_struct_type), POINTER :: fm_struct
1315 : TYPE(cp_fm_type) :: matrix_im_fm, matrix_re_fm
1316 :
1317 1119 : NULLIFY (fm_struct)
1318 0 : CPASSERT(ASSOCIATED(matrix_re))
1319 1119 : CPASSERT(ASSOCIATED(matrix_im))
1320 1119 : CALL dbcsr_get_info(matrix_re, nfullrows_total=nrow, nfullcols_total=ncol)
1321 : CALL cp_fm_struct_create(fm_struct, nrow_global=nrow, ncol_global=ncol, &
1322 : context=preconditioner_env%ctxt, &
1323 1119 : para_env=preconditioner_env%para_env)
1324 1119 : CALL cp_fm_create(matrix_re_fm, fm_struct, name=TRIM(matrix_name)//' real')
1325 1119 : CALL cp_fm_create(matrix_im_fm, fm_struct, name=TRIM(matrix_name)//' imaginary')
1326 1119 : CALL cp_fm_struct_release(fm_struct)
1327 1119 : CALL copy_dbcsr_to_fm(matrix_re, matrix_re_fm)
1328 1119 : CALL copy_dbcsr_to_fm(matrix_im, matrix_im_fm)
1329 1119 : CALL cp_cfm_create(matrix, matrix_re_fm%matrix_struct, name=matrix_name)
1330 1119 : CALL cp_fm_to_cfm(matrix_re_fm, matrix_im_fm, matrix)
1331 1119 : CALL cp_fm_release(matrix_im_fm)
1332 1119 : CALL cp_fm_release(matrix_re_fm)
1333 :
1334 1119 : END SUBROUTINE dbcsr_pair_to_cfm
1335 :
1336 : ! **************************************************************************************************
1337 : !> \brief Allows for a restart of the preconditioner
1338 : !> depending on the method it purges all arrays or keeps them
1339 : !> \param qs_env ...
1340 : !> \param preconditioner ...
1341 : !> \param prec_type ...
1342 : !> \param nspins ...
1343 : ! **************************************************************************************************
1344 7611 : SUBROUTINE restart_preconditioner(qs_env, preconditioner, prec_type, nspins)
1345 :
1346 : TYPE(qs_environment_type), POINTER :: qs_env
1347 : TYPE(preconditioner_p_type), DIMENSION(:), POINTER :: preconditioner
1348 : INTEGER, INTENT(IN) :: prec_type, nspins
1349 :
1350 : INTEGER :: ispin
1351 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
1352 : TYPE(mp_para_env_type), POINTER :: para_env
1353 :
1354 7611 : NULLIFY (para_env, blacs_env)
1355 7611 : CALL get_qs_env(qs_env, para_env=para_env, blacs_env=blacs_env)
1356 :
1357 7611 : IF (ASSOCIATED(preconditioner)) THEN
1358 6776 : SELECT CASE (prec_type)
1359 : CASE (ot_precond_full_all, ot_precond_fermi_low_rank, &
1360 : ot_precond_full_all_covariant, &
1361 : ot_precond_full_single) ! these depend on the ks matrix
1362 3476 : DO ispin = 1, SIZE(preconditioner)
1363 1938 : CALL destroy_preconditioner(preconditioner(ispin)%preconditioner)
1364 3476 : DEALLOCATE (preconditioner(ispin)%preconditioner)
1365 : END DO
1366 1538 : DEALLOCATE (preconditioner)
1367 : CASE (ot_precond_none, ot_precond_full_kinetic, ot_precond_s_inverse, &
1368 : ot_precond_full_single_inverse) ! these are 'independent'
1369 : ! do nothing
1370 : CASE DEFAULT
1371 5238 : CPABORT("Unknown preconditioner type")
1372 : END SELECT
1373 : END IF
1374 :
1375 : ! add an OT preconditioner if none is present
1376 7611 : IF (.NOT. ASSOCIATED(preconditioner)) THEN
1377 7169 : SELECT CASE (prec_type)
1378 : CASE (ot_precond_full_all, ot_precond_fermi_low_rank, ot_precond_full_all_covariant, &
1379 : ot_precond_full_single_inverse)
1380 13959 : ALLOCATE (preconditioner(nspins))
1381 : CASE DEFAULT
1382 4564 : ALLOCATE (preconditioner(1))
1383 : END SELECT
1384 8749 : DO ispin = 1, SIZE(preconditioner)
1385 24190 : ALLOCATE (preconditioner(ispin)%preconditioner)
1386 : CALL init_preconditioner(preconditioner(ispin)%preconditioner, &
1387 : para_env=para_env, &
1388 8749 : blacs_env=blacs_env)
1389 : END DO
1390 : END IF
1391 :
1392 7611 : END SUBROUTINE restart_preconditioner
1393 :
1394 : ! **************************************************************************************************
1395 : !> \brief ...
1396 : !> \param qs_env ...
1397 : !> \param mos ...
1398 : !> \param matrix_ks ...
1399 : !> \param matrix_s ...
1400 : !> \param ot_preconditioner ...
1401 : !> \param prec_type ...
1402 : !> \param solver_type ...
1403 : !> \param energy_gap ...
1404 : !> \param nspins ...
1405 : !> \param has_unit_metric ...
1406 : !> \param convert_to_dbcsr ...
1407 : !> \param chol_type ...
1408 : !> \param full_mo_set ...
1409 : !> \param chebyshev_degree ...
1410 : !> \param low_rank_base ...
1411 : !> \param fermi_low_rank_max_rank ...
1412 : !> \param lattice_fft ...
1413 : !> \param lattice_fft_local_cells ...
1414 : ! **************************************************************************************************
1415 7611 : SUBROUTINE prepare_preconditioner(qs_env, mos, matrix_ks, matrix_s, &
1416 : ot_preconditioner, prec_type, solver_type, &
1417 : energy_gap, nspins, has_unit_metric, &
1418 : convert_to_dbcsr, chol_type, full_mo_set, chebyshev_degree, &
1419 : low_rank_base, fermi_low_rank_max_rank, lattice_fft, &
1420 : lattice_fft_local_cells)
1421 :
1422 : TYPE(qs_environment_type), POINTER :: qs_env
1423 : TYPE(mo_set_type), DIMENSION(:), INTENT(INOUT) :: mos
1424 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s
1425 : TYPE(preconditioner_p_type), DIMENSION(:), POINTER :: ot_preconditioner
1426 : INTEGER, INTENT(IN) :: prec_type, solver_type
1427 : REAL(dp), INTENT(IN) :: energy_gap
1428 : INTEGER, INTENT(IN) :: nspins
1429 : LOGICAL, INTENT(IN), OPTIONAL :: has_unit_metric, convert_to_dbcsr
1430 : INTEGER, INTENT(IN), OPTIONAL :: chol_type
1431 : LOGICAL, INTENT(IN), OPTIONAL :: full_mo_set
1432 : INTEGER, INTENT(IN), OPTIONAL :: chebyshev_degree, low_rank_base, &
1433 : fermi_low_rank_max_rank, lattice_fft, &
1434 : lattice_fft_local_cells
1435 :
1436 : CHARACTER(LEN=*), PARAMETER :: routineN = 'prepare_preconditioner'
1437 :
1438 : CHARACTER(LEN=default_string_length) :: msg
1439 : INTEGER :: handle, icall, ispin, my_chebyshev_degree, my_fermi_low_rank_max_rank, &
1440 : my_lattice_fft, my_lattice_fft_local_cells, my_low_rank_base, n_loops
1441 : INTEGER, DIMENSION(5) :: nocc, norb
1442 7611 : INTEGER, DIMENSION(:), POINTER :: multiple_unit_cell
1443 : LOGICAL :: do_co_rotate, do_kpoints, my_convert_to_dbcsr, my_full_mo_set, &
1444 : my_has_unit_metric, use_covariant_full_all, use_mo_coeff_b
1445 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
1446 : TYPE(cp_fm_type), POINTER :: mo_coeff
1447 7611 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: kinetic
1448 : TYPE(dbcsr_type), POINTER :: matrix_t, mo_coeff_b
1449 : TYPE(dft_control_type), POINTER :: dft_control
1450 7611 : TYPE(mo_set_type), ALLOCATABLE, DIMENSION(:) :: canonical_mos
1451 : TYPE(mp_para_env_type), POINTER :: para_env
1452 : TYPE(section_vals_type), POINTER :: input, subsys_section
1453 :
1454 7611 : CALL timeset(routineN, handle)
1455 7611 : NULLIFY (matrix_t, mo_coeff_b, mo_coeff, kinetic, dft_control, para_env, blacs_env, &
1456 7611 : input, subsys_section, multiple_unit_cell)
1457 7611 : my_has_unit_metric = .FALSE.
1458 7611 : IF (PRESENT(has_unit_metric)) my_has_unit_metric = has_unit_metric
1459 7611 : my_convert_to_dbcsr = .TRUE.
1460 7611 : IF (PRESENT(convert_to_dbcsr)) my_convert_to_dbcsr = convert_to_dbcsr
1461 7611 : my_full_mo_set = .FALSE.
1462 7611 : IF (PRESENT(full_mo_set)) my_full_mo_set = full_mo_set
1463 7611 : my_chebyshev_degree = 8
1464 7611 : IF (PRESENT(chebyshev_degree)) my_chebyshev_degree = chebyshev_degree
1465 7611 : my_low_rank_base = ot_low_rank_base_overlap
1466 7611 : IF (PRESENT(low_rank_base)) my_low_rank_base = low_rank_base
1467 7611 : my_fermi_low_rank_max_rank = 48
1468 7611 : IF (PRESENT(fermi_low_rank_max_rank)) my_fermi_low_rank_max_rank = fermi_low_rank_max_rank
1469 7611 : my_lattice_fft = ot_lattice_fft_off
1470 7611 : IF (PRESENT(lattice_fft)) my_lattice_fft = lattice_fft
1471 7611 : my_lattice_fft_local_cells = 0
1472 7611 : IF (PRESENT(lattice_fft_local_cells)) my_lattice_fft_local_cells = lattice_fft_local_cells
1473 7593 : IF (my_lattice_fft_local_cells < 0) CPABORT("LATTICE_FFT_LOCAL_CELLS must not be negative")
1474 7611 : do_kpoints = .FALSE.
1475 7611 : use_covariant_full_all = prec_type == ot_precond_full_all_covariant
1476 :
1477 : CALL get_qs_env(qs_env, &
1478 : dft_control=dft_control, &
1479 : do_kpoints=do_kpoints, &
1480 : input=input, &
1481 : para_env=para_env, &
1482 7611 : blacs_env=blacs_env)
1483 :
1484 7611 : IF (my_low_rank_base == ot_low_rank_base_lattice_fft) THEN
1485 0 : IF (prec_type /= ot_precond_fermi_low_rank) THEN
1486 0 : CPABORT("FERMI_LOW_RANK_BASE LATTICE_FFT requires PRECONDITIONER FERMI_LOW_RANK")
1487 : END IF
1488 : my_lattice_fft = ot_lattice_fft_on
1489 7611 : ELSE IF (my_low_rank_base /= ot_low_rank_base_overlap) THEN
1490 0 : CPABORT("Unknown FERMI_LOW_RANK_BASE setting")
1491 : END IF
1492 :
1493 : IF (my_lattice_fft /= ot_lattice_fft_off .AND. &
1494 7611 : my_lattice_fft /= ot_lattice_fft_on .AND. &
1495 : my_lattice_fft /= ot_lattice_fft_auto) THEN
1496 0 : CPABORT("Unknown LATTICE_FFT setting")
1497 : END IF
1498 7611 : IF (my_lattice_fft_local_cells > 0 .AND. my_lattice_fft /= ot_lattice_fft_on) THEN
1499 0 : CPABORT("LATTICE_FFT_LOCAL_CELLS requires LATTICE_FFT ON")
1500 : END IF
1501 7611 : IF (my_lattice_fft /= ot_lattice_fft_off) THEN
1502 0 : SELECT CASE (prec_type)
1503 : CASE (ot_precond_fermi_low_rank, ot_precond_full_all_covariant, &
1504 : ot_precond_full_kinetic, ot_precond_full_single_inverse, ot_precond_s_inverse)
1505 : CASE DEFAULT
1506 18 : CPABORT("LATTICE_FFT does not support the selected OT preconditioner")
1507 : END SELECT
1508 18 : IF (solver_type /= ot_precond_solver_default .AND. &
1509 : solver_type /= ot_precond_solver_inv_chol) THEN
1510 0 : CPABORT("LATTICE_FFT requires PRECOND_SOLVER DEFAULT or INVERSE_CHOLESKY")
1511 : END IF
1512 18 : IF (my_has_unit_metric .OR. do_kpoints .OR. .NOT. ASSOCIATED(input)) THEN
1513 0 : IF (my_lattice_fft == ot_lattice_fft_on) THEN
1514 0 : CPABORT("LATTICE_FFT ON requires a non-unit metric in a real Gamma-point calculation")
1515 : END IF
1516 : my_lattice_fft = ot_lattice_fft_off
1517 : END IF
1518 : END IF
1519 : IF (my_lattice_fft /= ot_lattice_fft_off) THEN
1520 18 : subsys_section => section_vals_get_subs_vals(input, "SUBSYS")
1521 18 : IF (.NOT. ASSOCIATED(subsys_section)) THEN
1522 0 : IF (my_lattice_fft == ot_lattice_fft_on) CPABORT("LATTICE_FFT ON requires a SUBSYS section")
1523 : my_lattice_fft = ot_lattice_fft_off
1524 : END IF
1525 : IF (my_lattice_fft /= ot_lattice_fft_off) THEN
1526 : CALL section_vals_val_get(subsys_section, "CELL%MULTIPLE_UNIT_CELL", &
1527 18 : i_vals=multiple_unit_cell)
1528 18 : IF (.NOT. ASSOCIATED(multiple_unit_cell)) THEN
1529 0 : IF (my_lattice_fft == ot_lattice_fft_on) THEN
1530 0 : CPABORT("LATTICE_FFT ON requires a valid CELL%MULTIPLE_UNIT_CELL")
1531 : END IF
1532 : my_lattice_fft = ot_lattice_fft_off
1533 : ELSE
1534 72 : IF (SIZE(multiple_unit_cell) /= 3 .OR. ANY(multiple_unit_cell <= 0)) THEN
1535 0 : IF (my_lattice_fft == ot_lattice_fft_on) THEN
1536 0 : CPABORT("LATTICE_FFT ON requires a valid CELL%MULTIPLE_UNIT_CELL")
1537 : END IF
1538 : my_lattice_fft = ot_lattice_fft_off
1539 72 : ELSE IF (PRODUCT(multiple_unit_cell) <= 1 .AND. my_lattice_fft == ot_lattice_fft_on) THEN
1540 0 : CPABORT("LATTICE_FFT ON requires an explicitly replicated multiple unit cell")
1541 : END IF
1542 72 : IF (my_lattice_fft_local_cells >= PRODUCT(multiple_unit_cell)) THEN
1543 0 : CPABORT("LATTICE_FFT_LOCAL_CELLS must be smaller than the replicated cell count")
1544 : END IF
1545 : END IF
1546 : END IF
1547 : END IF
1548 :
1549 7611 : IF (dft_control%qs_control%semi_empirical .OR. dft_control%qs_control%dftb .OR. &
1550 : dft_control%qs_control%xtb) THEN
1551 2152 : IF (prec_type == ot_precond_full_kinetic) THEN
1552 0 : msg = "Full_kinetic not available for semi-empirical methods"
1553 0 : CPABORT(TRIM(msg))
1554 : END IF
1555 2152 : matrix_t => matrix_s(1)%matrix
1556 : ELSE
1557 5459 : CPASSERT(.NOT. my_has_unit_metric)
1558 5459 : CALL get_qs_env(qs_env, kinetic=kinetic)
1559 5459 : matrix_t => kinetic(1)%matrix
1560 : END IF
1561 :
1562 : ! use full set of MOs or just occupied MOs
1563 7611 : nocc = 0
1564 7611 : norb = 0
1565 7611 : IF (my_full_mo_set) THEN
1566 38 : DO ispin = 1, nspins
1567 20 : CALL get_mo_set(mo_set=mos(ispin), homo=nocc(ispin), nmo=norb(ispin))
1568 38 : CALL set_mo_set(mo_set=mos(ispin), homo=norb(ispin))
1569 : END DO
1570 : END IF
1571 : !determines how often make preconditioner is called, spin dependent methods have to be called twice
1572 7611 : n_loops = 1
1573 7611 : IF (prec_type == ot_precond_full_single_inverse) n_loops = nspins
1574 : ! check whether we need the ev and rotate the MOs
1575 2562 : SELECT CASE (prec_type)
1576 : CASE (ot_precond_full_all)
1577 : ! FULL_ALL needs canonical occupied eigenstates and one preconditioner per spin.
1578 2562 : n_loops = nspins
1579 :
1580 2562 : do_co_rotate = ASSOCIATED(qs_env%mo_derivs)
1581 5920 : DO ispin = 1, nspins
1582 3358 : CALL get_mo_set(mo_set=mos(ispin), mo_coeff_b=mo_coeff_b, mo_coeff=mo_coeff)
1583 3358 : use_mo_coeff_b = mos(ispin)%use_mo_coeff_b
1584 5920 : IF (use_mo_coeff_b .AND. do_co_rotate) THEN
1585 : CALL calculate_subspace_eigenvalues(mo_coeff_b, matrix_ks(ispin)%matrix, &
1586 : do_rotation=.TRUE., &
1587 : co_rotate=qs_env%mo_derivs(ispin)%matrix, &
1588 : para_env=para_env, &
1589 3338 : blacs_env=blacs_env)
1590 20 : ELSE IF (use_mo_coeff_b) THEN
1591 : CALL calculate_subspace_eigenvalues(mo_coeff_b, matrix_ks(ispin)%matrix, &
1592 : do_rotation=.TRUE., &
1593 : para_env=para_env, &
1594 20 : blacs_env=blacs_env)
1595 : ELSE
1596 : CALL calculate_subspace_eigenvalues(mo_coeff, matrix_ks(ispin)%matrix, &
1597 0 : do_rotation=.TRUE.)
1598 : END IF
1599 : END DO
1600 : CASE (ot_precond_fermi_low_rank)
1601 : ! The common-reference low-rank operator depends only on the occupied projector.
1602 : ! Keep the input gauge unchanged so that OT ROTATION remains well defined.
1603 6 : n_loops = nspins
1604 : CASE (ot_precond_full_all_covariant)
1605 : ! Diagonalize a private occupied copy. The physical MOs and their
1606 : ! derivatives remain in the user-selected gauge.
1607 8 : n_loops = nspins
1608 32 : ALLOCATE (canonical_mos(nspins))
1609 7627 : DO ispin = 1, nspins
1610 8 : CALL duplicate_mo_set(canonical_mos(ispin), mos(ispin))
1611 8 : CALL get_mo_set(mo_set=canonical_mos(ispin), mo_coeff_b=mo_coeff_b, mo_coeff=mo_coeff)
1612 8 : IF (canonical_mos(ispin)%use_mo_coeff_b) THEN
1613 8 : CALL copy_dbcsr_to_fm(mo_coeff_b, mo_coeff)
1614 8 : canonical_mos(ispin)%use_mo_coeff_b = .FALSE.
1615 : END IF
1616 : CALL canonicalize_full_all_copy(ot_preconditioner(ispin)%preconditioner, &
1617 : mo_coeff, matrix_ks(ispin)%matrix, &
1618 16 : canonical_mos(ispin)%homo)
1619 : END DO
1620 : CASE DEFAULT
1621 : ! No need to rotate the MOs
1622 : END SELECT
1623 :
1624 : ! check whether we have a preconditioner
1625 724 : SELECT CASE (prec_type)
1626 : CASE (ot_precond_none)
1627 1448 : DO ispin = 1, SIZE(ot_preconditioner)
1628 1448 : ot_preconditioner(ispin)%preconditioner%in_use = 0
1629 : END DO
1630 : CASE DEFAULT
1631 22636 : DO icall = 1, n_loops
1632 8138 : ot_preconditioner(icall)%preconditioner%polynomial_degree = my_chebyshev_degree
1633 8138 : ot_preconditioner(icall)%preconditioner%low_rank_base = my_low_rank_base
1634 8138 : ot_preconditioner(icall)%preconditioner%lattice_fft_mode = my_lattice_fft
1635 8138 : ot_preconditioner(icall)%preconditioner%lattice_local_cells = my_lattice_fft_local_cells
1636 32552 : ot_preconditioner(icall)%preconditioner%lattice_dims = [1, 1, 1]
1637 8138 : IF (my_lattice_fft /= ot_lattice_fft_off) THEN
1638 18 : CPASSERT(ASSOCIATED(multiple_unit_cell))
1639 144 : ot_preconditioner(icall)%preconditioner%lattice_dims = multiple_unit_cell
1640 : END IF
1641 8138 : ot_preconditioner(icall)%preconditioner%spectral_max_rank = my_fermi_low_rank_max_rank
1642 15025 : IF (my_has_unit_metric) THEN
1643 474 : IF (use_covariant_full_all) THEN
1644 : CALL make_preconditioner(ot_preconditioner(icall)%preconditioner, &
1645 : prec_type, solver_type, &
1646 : matrix_h=matrix_ks(icall)%matrix, &
1647 : mo_set=canonical_mos(icall), &
1648 : energy_gap=energy_gap, &
1649 0 : convert_precond_to_dbcsr=my_convert_to_dbcsr)
1650 : ELSE
1651 : CALL make_preconditioner(ot_preconditioner(icall)%preconditioner, &
1652 : prec_type, solver_type, &
1653 : matrix_h=matrix_ks(icall)%matrix, &
1654 : mo_set=mos(icall), &
1655 : energy_gap=energy_gap, &
1656 474 : convert_precond_to_dbcsr=my_convert_to_dbcsr)
1657 : END IF
1658 : ELSE
1659 7664 : IF (use_covariant_full_all) THEN
1660 : CALL make_preconditioner(ot_preconditioner(icall)%preconditioner, &
1661 : prec_type, solver_type, &
1662 : matrix_h=matrix_ks(icall)%matrix, &
1663 : matrix_s=matrix_s(1)%matrix, matrix_t=matrix_t, &
1664 : mo_set=canonical_mos(icall), &
1665 : energy_gap=energy_gap, &
1666 8 : convert_precond_to_dbcsr=my_convert_to_dbcsr, chol_type=chol_type)
1667 : ELSE
1668 : CALL make_preconditioner(ot_preconditioner(icall)%preconditioner, &
1669 : prec_type, solver_type, &
1670 : matrix_h=matrix_ks(icall)%matrix, &
1671 : matrix_s=matrix_s(1)%matrix, matrix_t=matrix_t, &
1672 : mo_set=mos(icall), energy_gap=energy_gap, &
1673 7656 : convert_precond_to_dbcsr=my_convert_to_dbcsr, chol_type=chol_type)
1674 : END IF
1675 : END IF
1676 : END DO
1677 : END SELECT
1678 :
1679 7611 : IF (ALLOCATED(canonical_mos)) THEN
1680 16 : DO ispin = 1, SIZE(canonical_mos)
1681 16 : CALL deallocate_mo_set(canonical_mos(ispin))
1682 : END DO
1683 8 : DEALLOCATE (canonical_mos)
1684 : END IF
1685 :
1686 : ! reset homo values
1687 7611 : IF (my_full_mo_set) THEN
1688 38 : DO ispin = 1, nspins
1689 38 : CALL set_mo_set(mo_set=mos(ispin), homo=nocc(ispin))
1690 : END DO
1691 : END IF
1692 :
1693 7611 : CALL timestop(handle)
1694 :
1695 15222 : END SUBROUTINE prepare_preconditioner
1696 :
1697 : ! **************************************************************************************************
1698 : !> \brief Canonicalize an occupied copy and retain the occupied-space rotation.
1699 : !> \param preconditioner_env ...
1700 : !> \param mo_coeff ...
1701 : !> \param matrix_h ...
1702 : !> \param nocc ...
1703 : ! **************************************************************************************************
1704 8 : SUBROUTINE canonicalize_full_all_copy(preconditioner_env, mo_coeff, matrix_h, nocc)
1705 :
1706 : TYPE(preconditioner_type) :: preconditioner_env
1707 : TYPE(cp_fm_type), POINTER :: mo_coeff
1708 : TYPE(dbcsr_type), POINTER :: matrix_h
1709 : INTEGER, INTENT(IN) :: nocc
1710 :
1711 : CHARACTER(LEN=*), PARAMETER :: routineN = 'canonicalize_full_all_copy'
1712 :
1713 : INTEGER :: handle, nao
1714 8 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: occupied_evals
1715 : TYPE(cp_fm_struct_type), POINTER :: fm_struct
1716 : TYPE(cp_fm_type) :: canonical_coeff, h_block, weighted_coeff
1717 :
1718 8 : CALL timeset(routineN, handle)
1719 8 : CALL cp_fm_get_info(mo_coeff, nrow_global=nao)
1720 8 : CPASSERT(nocc > 0)
1721 :
1722 : CALL cp_fm_struct_create(fm_struct, nrow_global=nao, ncol_global=nocc, &
1723 : context=preconditioner_env%ctxt, &
1724 8 : para_env=preconditioner_env%para_env)
1725 8 : CALL cp_fm_create(canonical_coeff, fm_struct, name=routineN)
1726 8 : CALL cp_fm_create(weighted_coeff, fm_struct, name=routineN)
1727 8 : CALL cp_fm_struct_release(fm_struct)
1728 8 : CALL cp_fm_to_fm(mo_coeff, canonical_coeff, nocc)
1729 :
1730 : CALL cp_fm_struct_create(fm_struct, nrow_global=nocc, ncol_global=nocc, &
1731 : context=preconditioner_env%ctxt, &
1732 8 : para_env=preconditioner_env%para_env)
1733 8 : CALL cp_fm_create(h_block, fm_struct, name=routineN)
1734 8 : IF (ASSOCIATED(preconditioner_env%occ_rotation_dbcsr)) THEN
1735 0 : CALL dbcsr_release_p(preconditioner_env%occ_rotation_dbcsr)
1736 : END IF
1737 8 : IF (ASSOCIATED(preconditioner_env%occ_rotation)) THEN
1738 0 : CALL cp_fm_release(preconditioner_env%occ_rotation)
1739 0 : DEALLOCATE (preconditioner_env%occ_rotation)
1740 : END IF
1741 8 : ALLOCATE (preconditioner_env%occ_rotation)
1742 8 : CALL cp_fm_create(preconditioner_env%occ_rotation, fm_struct, name="occupied rotation")
1743 8 : CALL cp_fm_struct_release(fm_struct)
1744 :
1745 8 : CALL cp_dbcsr_sm_fm_multiply(matrix_h, canonical_coeff, weighted_coeff, nocc)
1746 : CALL parallel_gemm('T', 'N', nocc, nocc, nao, 1.0_dp, canonical_coeff, &
1747 8 : weighted_coeff, 0.0_dp, h_block)
1748 24 : ALLOCATE (occupied_evals(nocc))
1749 8 : CALL choose_eigv_solver(h_block, preconditioner_env%occ_rotation, occupied_evals)
1750 :
1751 : CALL parallel_gemm('N', 'N', nao, nocc, nocc, 1.0_dp, canonical_coeff, &
1752 8 : preconditioner_env%occ_rotation, 0.0_dp, weighted_coeff)
1753 8 : CALL cp_fm_to_fm(weighted_coeff, mo_coeff, nocc, 1, 1)
1754 :
1755 8 : DEALLOCATE (occupied_evals)
1756 8 : CALL cp_fm_release(canonical_coeff)
1757 8 : CALL cp_fm_release(h_block)
1758 8 : CALL cp_fm_release(weighted_coeff)
1759 8 : CALL timestop(handle)
1760 :
1761 24 : END SUBROUTINE canonicalize_full_all_copy
1762 :
1763 : END MODULE preconditioner
|
