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_type
25 : USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm
26 : USE cp_fm_basic_linalg, ONLY: cp_fm_uplo_to_full
27 : USE cp_fm_cholesky, ONLY: cp_fm_cholesky_decompose,&
28 : cp_fm_cholesky_invert
29 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
30 : cp_fm_struct_release,&
31 : cp_fm_struct_type
32 : USE cp_fm_types, ONLY: cp_fm_create,&
33 : cp_fm_get_info,&
34 : cp_fm_release,&
35 : cp_fm_set_all,&
36 : cp_fm_to_fm,&
37 : cp_fm_type
38 : USE cp_log_handling, ONLY: cp_get_default_logger,&
39 : cp_logger_get_default_io_unit,&
40 : cp_logger_type
41 : USE input_constants, ONLY: &
42 : cholesky_reduce, ot_low_rank_base_lattice_fft, ot_low_rank_base_overlap, &
43 : ot_precond_fermi_low_rank, ot_precond_full_all, ot_precond_full_kinetic, &
44 : ot_precond_full_single, ot_precond_full_single_inverse, ot_precond_none, &
45 : ot_precond_s_inverse, ot_precond_solver_default, ot_precond_solver_inv_chol, &
46 : ot_precond_solver_update
47 : USE input_section_types, ONLY: section_vals_get_subs_vals,&
48 : section_vals_type,&
49 : section_vals_val_get
50 : USE kinds, ONLY: default_string_length,&
51 : dp
52 : USE lattice_low_rank_preconditioner, ONLY: build_lattice_inverse
53 : USE low_rank_preconditioner_model, ONLY: low_rank_inverse_weight,&
54 : low_rank_select_rank
55 : USE message_passing, ONLY: mp_para_env_type
56 : USE preconditioner_apply, ONLY: apply_preconditioner_cfm_complex,&
57 : apply_preconditioner_dbcsr,&
58 : apply_preconditioner_dbcsr_complex,&
59 : apply_preconditioner_fm
60 : USE preconditioner_makes, ONLY: make_complex_full_all,&
61 : make_complex_full_kinetic,&
62 : make_complex_full_s_inverse,&
63 : make_complex_full_single,&
64 : make_complex_full_single_inverse,&
65 : make_preconditioner_matrix
66 : USE preconditioner_solvers, ONLY: solve_preconditioner,&
67 : transfer_dbcsr_to_fm,&
68 : transfer_fm_to_dbcsr
69 : USE preconditioner_types, ONLY: destroy_preconditioner,&
70 : init_preconditioner,&
71 : preconditioner_p_type,&
72 : preconditioner_type
73 : USE qs_environment_types, ONLY: get_qs_env,&
74 : qs_environment_type
75 : USE qs_mo_methods, ONLY: calculate_subspace_eigenvalues
76 : USE qs_mo_types, ONLY: get_mo_set,&
77 : mo_set_type,&
78 : set_mo_set
79 : #include "./base/base_uses.f90"
80 :
81 : IMPLICIT NONE
82 :
83 : PRIVATE
84 :
85 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'preconditioner'
86 : INTEGER, PARAMETER, PRIVATE :: fermi_low_rank_max_rank = 48
87 : REAL(KIND=dp), PARAMETER, PRIVATE :: fermi_low_rank_gap_factor = 4.0_dp, &
88 : fermi_low_rank_degeneracy_tolerance = 1.0E-8_dp, &
89 : fermi_low_rank_min_window = 1.0_dp
90 :
91 : PUBLIC :: make_preconditioner, make_preconditioner_complex_full_all, &
92 : make_preconditioner_complex_full_kinetic, &
93 : make_preconditioner_complex_full_s_inverse, &
94 : make_preconditioner_complex_full_single, &
95 : make_preconditioner_complex_full_single_inverse, restart_preconditioner, &
96 : dbcsr_pair_to_cfm
97 : PUBLIC :: apply_preconditioner, prepare_preconditioner
98 :
99 : ! The public interface for apply preconditioner, the routines can be found in preconditioner_apply.F
100 : INTERFACE apply_preconditioner
101 : MODULE PROCEDURE apply_preconditioner_cfm_complex
102 : MODULE PROCEDURE apply_preconditioner_dbcsr
103 : MODULE PROCEDURE apply_preconditioner_dbcsr_complex
104 : MODULE PROCEDURE apply_preconditioner_fm
105 : END INTERFACE
106 :
107 : ! **************************************************************************************************
108 :
109 : CONTAINS
110 :
111 : ! **************************************************************************************************
112 :
113 : ! **************************************************************************************************
114 : !> \brief Build the overlap-inverse base of the bounded low-rank spectral preconditioner.
115 : !> \param preconditioner_env ...
116 : !> \param matrix_s ...
117 : ! **************************************************************************************************
118 4 : SUBROUTINE prepare_low_rank_base(preconditioner_env, matrix_s)
119 :
120 : TYPE(preconditioner_type) :: preconditioner_env
121 : TYPE(dbcsr_type), OPTIONAL, POINTER :: matrix_s
122 :
123 : INTEGER :: i, info, j, n, ncol_local, nrow_local, &
124 : output_unit
125 4 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
126 : LOGICAL :: used_fft
127 4 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: overlap
128 : REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :), &
129 4 : POINTER :: local_data
130 : TYPE(cp_fm_type) :: overlap_fm, work
131 : TYPE(cp_logger_type), POINTER :: logger
132 :
133 0 : CPASSERT(ASSOCIATED(preconditioner_env%fm))
134 4 : CPASSERT(.NOT. ASSOCIATED(preconditioner_env%base_fm))
135 :
136 6 : SELECT CASE (preconditioner_env%low_rank_base)
137 : CASE (ot_low_rank_base_overlap)
138 2 : ALLOCATE (preconditioner_env%base_fm)
139 : CALL cp_fm_create(preconditioner_env%base_fm, preconditioner_env%fm%matrix_struct, &
140 2 : name="FERMI_LOW_RANK overlap inverse")
141 :
142 2 : IF (PRESENT(matrix_s)) THEN
143 2 : CPASSERT(ASSOCIATED(matrix_s))
144 2 : CALL copy_dbcsr_to_fm(matrix_s, preconditioner_env%base_fm)
145 2 : CALL cp_fm_cholesky_decompose(preconditioner_env%base_fm, info_out=info)
146 2 : IF (info /= 0) CPABORT("FERMI_LOW_RANK requires a positive-definite overlap matrix")
147 2 : CALL cp_fm_cholesky_invert(preconditioner_env%base_fm)
148 2 : CALL cp_fm_create(work, preconditioner_env%base_fm%matrix_struct, name="low-rank base work")
149 2 : CALL cp_fm_uplo_to_full(preconditioner_env%base_fm, work)
150 2 : CALL cp_fm_release(work)
151 : ELSE
152 0 : CALL cp_fm_set_all(preconditioner_env%base_fm, alpha=0.0_dp, beta=1.0_dp)
153 : END IF
154 :
155 : CASE (ot_low_rank_base_lattice_fft)
156 2 : IF (.NOT. PRESENT(matrix_s)) THEN
157 0 : CPABORT("FERMI_LOW_RANK_BASE LATTICE_FFT requires a non-unit overlap metric")
158 : END IF
159 2 : CPASSERT(ASSOCIATED(matrix_s))
160 2 : CALL cp_fm_get_info(preconditioner_env%fm, nrow_global=n)
161 : CALL cp_fm_create(overlap_fm, preconditioner_env%fm%matrix_struct, &
162 2 : name="low-rank lattice overlap")
163 2 : CALL cp_fm_set_all(overlap_fm, alpha=0.0_dp, beta=0.0_dp)
164 2 : CALL copy_dbcsr_to_fm(matrix_s, overlap_fm)
165 2 : CALL cp_fm_create(work, overlap_fm%matrix_struct, name="lattice overlap work")
166 2 : CALL cp_fm_uplo_to_full(overlap_fm, work)
167 2 : CALL cp_fm_release(work)
168 : CALL cp_fm_get_info(overlap_fm, nrow_local=nrow_local, ncol_local=ncol_local, &
169 2 : row_indices=row_indices, col_indices=col_indices, local_data=local_data)
170 8 : ALLOCATE (overlap(n, n), source=0.0_dp)
171 134 : DO j = 1, ncol_local
172 4490 : DO i = 1, nrow_local
173 4488 : IF (row_indices(i) <= col_indices(j)) THEN
174 2211 : overlap(row_indices(i), col_indices(j)) = local_data(i, j)
175 : END IF
176 : END DO
177 : END DO
178 2 : CALL preconditioner_env%para_env%sum(overlap)
179 134 : DO j = 1, n
180 4424 : DO i = j + 1, n
181 4422 : overlap(i, j) = overlap(j, i)
182 : END DO
183 : END DO
184 : CALL build_lattice_inverse(overlap, preconditioner_env%lattice_dims, &
185 : preconditioner_env%lattice_inverse_k, &
186 2 : preconditioner_env%lattice_projection_error, used_fft, info)
187 2 : CALL cp_fm_release(overlap_fm)
188 2 : DEALLOCATE (overlap)
189 2 : IF (info /= 0) THEN
190 0 : CPABORT("FERMI_LOW_RANK_BASE LATTICE_FFT found a non-positive Fourier overlap block")
191 : END IF
192 2 : IF (.NOT. used_fft) THEN
193 0 : CPABORT("FERMI_LOW_RANK_BASE LATTICE_FFT could not initialize its lattice FFT")
194 : END IF
195 : preconditioner_env%lattice_block_size = &
196 8 : n/PRODUCT(preconditioner_env%lattice_dims)
197 9 : IF (preconditioner_env%para_env%is_source()) THEN
198 1 : logger => cp_get_default_logger()
199 1 : output_unit = cp_logger_get_default_io_unit(logger)
200 1 : IF (output_unit > 0) THEN
201 : WRITE (output_unit, '(T2,A,3(I0,1X),A,I0)') &
202 1 : "OT| FERMI_LOW_RANK lattice FFT cells: ", preconditioner_env%lattice_dims, &
203 2 : "AO block: ", preconditioner_env%lattice_block_size
204 : WRITE (output_unit, '(T2,A,ES12.4)') &
205 1 : "OT| Block-circulant overlap projection error: ", &
206 2 : preconditioner_env%lattice_projection_error
207 : END IF
208 : END IF
209 :
210 : CASE DEFAULT
211 4 : CPABORT("Unknown FERMI_LOW_RANK base operator")
212 : END SELECT
213 :
214 8 : END SUBROUTINE prepare_low_rank_base
215 :
216 : ! **************************************************************************************************
217 :
218 : ! creates a preconditioner for the system (H-energy_homo S)
219 : ! this preconditioner is (must be) symmetric positive definite.
220 : ! currently uses a atom-block-diagonal form
221 : ! each block will be ....
222 : ! might overwrite matrix_h, matrix_t
223 :
224 : ! **************************************************************************************************
225 : !> \brief ...
226 : !> \param preconditioner_env ...
227 : !> \param precon_type ...
228 : !> \param solver_type ...
229 : !> \param matrix_h ...
230 : !> \param matrix_s ...
231 : !> \param matrix_t ...
232 : !> \param mo_set ...
233 : !> \param energy_gap ...
234 : !> \param convert_precond_to_dbcsr ...
235 : !> \param chol_type ...
236 : !> \par History
237 : !> 09.2014 removed some unused or unfinished methods
238 : !> removed sparse preconditioners and the
239 : !> sparse approximate inverse at rev 14341 [Florian Schiffmann]
240 : ! **************************************************************************************************
241 9704 : SUBROUTINE make_preconditioner(preconditioner_env, precon_type, solver_type, matrix_h, matrix_s, &
242 : matrix_t, mo_set, energy_gap, convert_precond_to_dbcsr, chol_type)
243 :
244 : TYPE(preconditioner_type) :: preconditioner_env
245 : INTEGER, INTENT(IN) :: precon_type, solver_type
246 : TYPE(dbcsr_type), POINTER :: matrix_h
247 : TYPE(dbcsr_type), OPTIONAL, POINTER :: matrix_s, matrix_t
248 : TYPE(mo_set_type), INTENT(IN) :: mo_set
249 : REAL(KIND=dp) :: energy_gap
250 : LOGICAL, INTENT(IN), OPTIONAL :: convert_precond_to_dbcsr
251 : INTEGER, INTENT(IN), OPTIONAL :: chol_type
252 :
253 : CHARACTER(len=*), PARAMETER :: routineN = 'make_preconditioner'
254 :
255 : INTEGER :: handle, k, my_solver_type, nao, nhomo, &
256 : output_unit
257 : LOGICAL :: my_convert_precond_to_dbcsr, &
258 : needs_full_spectrum, needs_homo, &
259 : use_mo_coeff_b
260 : REAL(KIND=dp) :: energy_homo
261 9704 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eigenvalues_ot
262 : TYPE(cp_fm_struct_type), POINTER :: fm_struct
263 : TYPE(cp_fm_type) :: mo_occ
264 : TYPE(cp_fm_type), POINTER :: mo_coeff
265 : TYPE(cp_logger_type), POINTER :: logger
266 : TYPE(dbcsr_type), POINTER :: mo_coeff_b
267 :
268 9704 : CALL timeset(routineN, handle)
269 :
270 9704 : CALL get_mo_set(mo_set=mo_set, mo_coeff=mo_coeff, mo_coeff_b=mo_coeff_b, homo=nhomo)
271 9704 : use_mo_coeff_b = mo_set%use_mo_coeff_b
272 9704 : CALL cp_fm_get_info(mo_coeff, ncol_global=k, nrow_global=nao)
273 :
274 : ! Starting some matrix mess, check where to store the result in preconditioner_env, fm or dbcsr_matrix
275 9704 : my_convert_precond_to_dbcsr = .FALSE.
276 9704 : IF (PRESENT(convert_precond_to_dbcsr)) my_convert_precond_to_dbcsr = convert_precond_to_dbcsr
277 9704 : IF (precon_type == ot_precond_fermi_low_rank .AND. &
278 : preconditioner_env%low_rank_base == ot_low_rank_base_lattice_fft) THEN
279 : ! The FFT base is an operator rather than a stored DBCSR matrix. Keep the
280 : ! complementary vectors in their distributed full-matrix representation.
281 2 : my_convert_precond_to_dbcsr = .FALSE.
282 : END IF
283 :
284 : ! Thanks to the mess with the matrices we need to make sure in this case that the
285 : ! Previous inverse is properly stored as a sparse matrix, fm gets deallocated here
286 : ! if it wasn't anyway
287 9704 : IF (preconditioner_env%solver == ot_precond_solver_update) THEN
288 4 : CALL transfer_fm_to_dbcsr(preconditioner_env%fm, preconditioner_env%dbcsr_matrix, matrix_h)
289 : END IF
290 :
291 9704 : needs_full_spectrum = .FALSE.
292 9704 : needs_homo = .FALSE.
293 :
294 3766 : SELECT CASE (precon_type)
295 : CASE (ot_precond_full_all, ot_precond_fermi_low_rank)
296 3766 : needs_full_spectrum = .TRUE.
297 : ! both of them need the coefficients as fm's, more matrix mess
298 3766 : IF (use_mo_coeff_b) THEN
299 3462 : CALL copy_dbcsr_to_fm(mo_coeff_b, mo_coeff)
300 : END IF
301 : CASE (ot_precond_full_single)
302 38 : needs_homo = .TRUE.
303 : ! XXXX to be removed if homo estimate only is implemented
304 38 : needs_full_spectrum = .TRUE.
305 : CASE (ot_precond_full_kinetic, ot_precond_s_inverse, ot_precond_full_single_inverse)
306 : ! these should be happy without an estimate for the homo energy
307 : ! preconditioning can not depend on an absolute eigenvalue, only on eigenvalue differences
308 : CASE DEFAULT
309 9704 : CPABORT("The preconditioner is unknown ...")
310 : END SELECT
311 :
312 28928 : ALLOCATE (eigenvalues_ot(k))
313 9704 : energy_homo = 0.0_dp
314 9704 : IF (needs_full_spectrum) THEN
315 : ! XXXXXXXXXXXXXXXX do not touch the initial MOs, could be harmful for either
316 : ! the case of non-equivalent MOs but also for the derivate
317 : ! we could already have all eigenvalues e.g. full_all and we could skip this
318 : ! to be optimised later.
319 : ! one flaw is that not all SCF methods (i.e. that go over mo_derivs directly)
320 : ! have a 'valid' matrix_h... (we even don't know what evals are in that case)
321 3804 : IF (use_mo_coeff_b) THEN
322 : CALL calculate_subspace_eigenvalues(mo_coeff_b, matrix_h, &
323 : eigenvalues_ot, do_rotation=.FALSE., &
324 : para_env=mo_coeff%matrix_struct%para_env, &
325 3492 : blacs_env=mo_coeff%matrix_struct%context)
326 : ELSE
327 : CALL calculate_subspace_eigenvalues(mo_coeff, matrix_h, &
328 312 : eigenvalues_ot, do_rotation=.FALSE.)
329 : END IF
330 3804 : IF (k > 0) THEN
331 3694 : CPASSERT(nhomo > 0 .AND. nhomo <= k)
332 3694 : energy_homo = eigenvalues_ot(nhomo)
333 : END IF
334 : ELSE
335 5900 : IF (needs_homo) THEN
336 0 : CPABORT("Not yet implemented")
337 : END IF
338 : END IF
339 :
340 : ! After all bits and pieces of checking and initialization, here comes the
341 : ! part where the preconditioner matrix gets created and solved.
342 : ! This will give the matrices for later use
343 9704 : my_solver_type = solver_type
344 9704 : preconditioner_env%in_use = precon_type
345 9704 : preconditioner_env%cholesky_use = cholesky_reduce
346 9704 : IF (PRESENT(chol_type)) preconditioner_env%cholesky_use = chol_type
347 : preconditioner_env%in_use = precon_type
348 9704 : IF (nhomo == k) THEN
349 : CALL make_preconditioner_matrix(preconditioner_env, matrix_h, matrix_s, matrix_t, mo_coeff, &
350 9616 : energy_homo, eigenvalues_ot, energy_gap, my_solver_type)
351 : ELSE
352 : CALL cp_fm_struct_create(fm_struct, nrow_global=nao, ncol_global=nhomo, &
353 : context=preconditioner_env%ctxt, &
354 88 : para_env=preconditioner_env%para_env)
355 88 : CALL cp_fm_create(mo_occ, fm_struct)
356 88 : CALL cp_fm_to_fm(mo_coeff, mo_occ, nhomo)
357 88 : CALL cp_fm_struct_release(fm_struct)
358 : !
359 : CALL make_preconditioner_matrix(preconditioner_env, matrix_h, matrix_s, matrix_t, mo_occ, &
360 88 : energy_homo, eigenvalues_ot(1:nhomo), energy_gap, my_solver_type)
361 : !
362 88 : CALL cp_fm_release(mo_occ)
363 : END IF
364 :
365 9704 : CALL solve_preconditioner(my_solver_type, preconditioner_env, matrix_s, matrix_h)
366 :
367 9704 : IF (precon_type == ot_precond_fermi_low_rank) THEN
368 4 : preconditioner_env%spectral_reference = energy_homo
369 : preconditioner_env%spectral_window = &
370 4 : MAX(fermi_low_rank_min_window, fermi_low_rank_gap_factor*preconditioner_env%energy_gap)
371 4 : preconditioner_env%spectral_base_scale = 1.0_dp/preconditioner_env%spectral_window
372 : preconditioner_env%spectral_rank = low_rank_select_rank( &
373 : preconditioner_env%full_evals, nhomo, fermi_low_rank_max_rank, &
374 4 : fermi_low_rank_degeneracy_tolerance)
375 4 : IF (preconditioner_env%low_rank_base == ot_low_rank_base_lattice_fft) THEN
376 : ! With an approximate base, retain only non-negative spectral updates.
377 : ! This preserves positive definiteness independently of the projection error.
378 26 : DO WHILE (preconditioner_env%spectral_rank > 0)
379 26 : IF (low_rank_inverse_weight( &
380 : preconditioner_env%full_evals(nhomo + preconditioner_env%spectral_rank), &
381 : preconditioner_env%spectral_reference, preconditioner_env%energy_gap) >= &
382 : preconditioner_env%spectral_base_scale) EXIT
383 26 : preconditioner_env%spectral_rank = preconditioner_env%spectral_rank - 1
384 : END DO
385 : END IF
386 4 : IF (preconditioner_env%para_env%is_source()) THEN
387 2 : logger => cp_get_default_logger()
388 2 : output_unit = cp_logger_get_default_io_unit(logger)
389 2 : IF (output_unit > 0) THEN
390 : WRITE (output_unit, '(T2,A,F14.8,A,I0)') &
391 2 : "OT| FERMI_LOW_RANK common reference [hartree]: ", &
392 2 : preconditioner_env%spectral_reference, " retained rank: ", &
393 4 : preconditioner_env%spectral_rank
394 : END IF
395 : END IF
396 4 : CALL prepare_low_rank_base(preconditioner_env, matrix_s)
397 : END IF
398 :
399 : ! Here comes more matrix mess, make sure to output the correct matrix format,
400 : ! A bit pointless to convert the cholesky factorized version as it doesn't work in
401 : ! dbcsr form and will crash later,...
402 9704 : IF (my_convert_precond_to_dbcsr) THEN
403 7972 : CALL transfer_fm_to_dbcsr(preconditioner_env%fm, preconditioner_env%dbcsr_matrix, matrix_h)
404 7972 : IF (precon_type == ot_precond_fermi_low_rank .AND. &
405 : preconditioner_env%low_rank_base == ot_low_rank_base_overlap) THEN
406 : CALL transfer_fm_to_dbcsr(preconditioner_env%base_fm, &
407 2 : preconditioner_env%base_dbcsr_matrix, matrix_h)
408 : END IF
409 : ELSE
410 : CALL transfer_dbcsr_to_fm(preconditioner_env%dbcsr_matrix, preconditioner_env%fm, &
411 1732 : preconditioner_env%para_env, preconditioner_env%ctxt)
412 : END IF
413 :
414 9704 : DEALLOCATE (eigenvalues_ot)
415 :
416 9704 : CALL timestop(handle)
417 :
418 9704 : END SUBROUTINE make_preconditioner
419 :
420 : ! **************************************************************************************************
421 : !> \brief Construct FULL_ALL directly from one complex H(k), S(k), and C(k) channel.
422 : !> \param preconditioner_env preconditioner storage
423 : !> \param matrix_c_re real part of the active OT reference orbitals
424 : !> \param matrix_c_im imaginary part of the active OT reference orbitals
425 : !> \param matrix_h_re real part of H(k)
426 : !> \param matrix_h_im imaginary part of H(k)
427 : !> \param matrix_s_re real part of S(k)
428 : !> \param matrix_s_im imaginary part of S(k)
429 : !> \param mo_set reference-orbital energy labels
430 : !> \param energy_gap denominator floor
431 : !> \param solver_type inverse application selected in the OT input
432 : ! **************************************************************************************************
433 520 : SUBROUTINE make_preconditioner_complex_full_all(preconditioner_env, &
434 : matrix_c_re, matrix_c_im, &
435 : matrix_h_re, matrix_h_im, &
436 : matrix_s_re, matrix_s_im, &
437 : mo_set, energy_gap, solver_type)
438 :
439 : TYPE(preconditioner_type) :: preconditioner_env
440 : TYPE(dbcsr_type), POINTER :: matrix_c_re, matrix_c_im, matrix_h_re, &
441 : matrix_h_im, matrix_s_re, matrix_s_im
442 : TYPE(mo_set_type), INTENT(IN) :: mo_set
443 : REAL(KIND=dp), INTENT(IN) :: energy_gap
444 : INTEGER, INTENT(IN) :: solver_type
445 :
446 : CHARACTER(len=*), PARAMETER :: routineN = 'make_preconditioner_complex_full_all'
447 :
448 : INTEGER :: handle, k
449 104 : REAL(KIND=dp), DIMENSION(:), POINTER :: eigenvalues
450 : TYPE(cp_cfm_type) :: matrix_c, matrix_h, matrix_s
451 :
452 104 : CALL timeset(routineN, handle)
453 :
454 104 : NULLIFY (eigenvalues)
455 104 : CPASSERT(ASSOCIATED(matrix_c_re))
456 104 : CPASSERT(ASSOCIATED(matrix_c_im))
457 104 : CPASSERT(ASSOCIATED(matrix_h_re))
458 104 : CPASSERT(ASSOCIATED(matrix_h_im))
459 104 : CPASSERT(ASSOCIATED(matrix_s_re))
460 104 : CPASSERT(ASSOCIATED(matrix_s_im))
461 104 : IF (solver_type /= ot_precond_solver_default) THEN
462 0 : CPABORT('Complex FULL_ALL supports only PRECOND_SOLVER DEFAULT')
463 : END IF
464 104 : CALL dbcsr_get_info(matrix_c_re, nfullcols_total=k)
465 104 : CPASSERT(k > 0)
466 104 : CALL get_mo_set(mo_set, eigenvalues=eigenvalues)
467 104 : CPASSERT(ASSOCIATED(eigenvalues))
468 104 : CPASSERT(SIZE(eigenvalues) >= k)
469 :
470 : CALL dbcsr_pair_to_cfm(matrix_c_re, matrix_c_im, preconditioner_env, &
471 104 : 'complex FULL_ALL C', matrix_c)
472 : CALL dbcsr_pair_to_cfm(matrix_h_re, matrix_h_im, preconditioner_env, &
473 104 : 'complex FULL_ALL H', matrix_h)
474 : CALL dbcsr_pair_to_cfm(matrix_s_re, matrix_s_im, preconditioner_env, &
475 104 : 'complex FULL_ALL S', matrix_s)
476 : CALL make_complex_full_all(preconditioner_env, matrix_c, matrix_h, matrix_s, &
477 104 : eigenvalues(1:k), energy_gap)
478 :
479 104 : CALL cp_cfm_release(matrix_c)
480 104 : CALL cp_cfm_release(matrix_s)
481 104 : CALL cp_cfm_release(matrix_h)
482 :
483 104 : CALL timestop(handle)
484 :
485 104 : END SUBROUTINE make_preconditioner_complex_full_all
486 :
487 : ! **************************************************************************************************
488 : !> \brief Construct a complex FULL_SINGLE preconditioner from H(k) and S(k).
489 : !> \param preconditioner_env preconditioner storage
490 : !> \param matrix_h_re real part of H(k)
491 : !> \param matrix_h_im imaginary part of H(k)
492 : !> \param matrix_s_re real part of S(k)
493 : !> \param matrix_s_im imaginary part of S(k)
494 : !> \param mo_set orbital energy labels defining the occupied edge
495 : !> \param energy_gap denominator floor
496 : !> \param solver_type requested inverse solver
497 : ! **************************************************************************************************
498 104 : SUBROUTINE make_preconditioner_complex_full_single(preconditioner_env, &
499 : matrix_h_re, matrix_h_im, &
500 : matrix_s_re, matrix_s_im, &
501 : mo_set, energy_gap, solver_type)
502 :
503 : TYPE(preconditioner_type) :: preconditioner_env
504 : TYPE(dbcsr_type), POINTER :: matrix_h_re, matrix_h_im, matrix_s_re, &
505 : matrix_s_im
506 : TYPE(mo_set_type), INTENT(IN) :: mo_set
507 : REAL(KIND=dp), INTENT(IN) :: energy_gap
508 : INTEGER, INTENT(IN) :: solver_type
509 :
510 : INTEGER :: homo
511 26 : REAL(KIND=dp), DIMENSION(:), POINTER :: eigenvalues
512 : TYPE(cp_cfm_type) :: matrix_h, matrix_s
513 :
514 26 : NULLIFY (eigenvalues)
515 0 : CPASSERT(ASSOCIATED(matrix_h_re))
516 26 : CPASSERT(ASSOCIATED(matrix_h_im))
517 26 : CPASSERT(ASSOCIATED(matrix_s_re))
518 26 : CPASSERT(ASSOCIATED(matrix_s_im))
519 26 : IF (solver_type /= ot_precond_solver_default) THEN
520 0 : CPABORT('Complex FULL_SINGLE supports only PRECOND_SOLVER DEFAULT')
521 : END IF
522 26 : CALL get_mo_set(mo_set, homo=homo, eigenvalues=eigenvalues)
523 26 : CPASSERT(ASSOCIATED(eigenvalues))
524 26 : CPASSERT(homo > 0 .AND. homo <= SIZE(eigenvalues))
525 :
526 : CALL dbcsr_pair_to_cfm(matrix_h_re, matrix_h_im, preconditioner_env, &
527 26 : 'complex FULL_SINGLE H', matrix_h)
528 : CALL dbcsr_pair_to_cfm(matrix_s_re, matrix_s_im, preconditioner_env, &
529 26 : 'complex FULL_SINGLE S', matrix_s)
530 : CALL make_complex_full_single(preconditioner_env, matrix_h, matrix_s, &
531 26 : eigenvalues(homo), energy_gap)
532 26 : CALL cp_cfm_release(matrix_s)
533 26 : CALL cp_cfm_release(matrix_h)
534 :
535 26 : END SUBROUTINE make_preconditioner_complex_full_single
536 :
537 : ! **************************************************************************************************
538 : !> \brief Construct a complex FULL_SINGLE_INVERSE preconditioner without discarding Im(H,S,C).
539 : !> \param preconditioner_env preconditioner storage
540 : !> \param matrix_c_re real part of the occupied reference orbitals
541 : !> \param matrix_c_im imaginary part of the occupied reference orbitals
542 : !> \param matrix_h_re real part of H(k)
543 : !> \param matrix_h_im imaginary part of H(k)
544 : !> \param matrix_s_re real part of S(k)
545 : !> \param matrix_s_im imaginary part of S(k)
546 : !> \param energy_gap lower spectral bound
547 : !> \param solver_type requested inverse solver
548 : ! **************************************************************************************************
549 484 : SUBROUTINE make_preconditioner_complex_full_single_inverse(preconditioner_env, &
550 : matrix_c_re, matrix_c_im, &
551 : matrix_h_re, matrix_h_im, &
552 : matrix_s_re, matrix_s_im, &
553 : energy_gap, solver_type)
554 :
555 : TYPE(preconditioner_type) :: preconditioner_env
556 : TYPE(dbcsr_type), POINTER :: matrix_c_re, matrix_c_im, matrix_h_re, &
557 : matrix_h_im, matrix_s_re, matrix_s_im
558 : REAL(KIND=dp), INTENT(IN) :: energy_gap
559 : INTEGER, INTENT(IN) :: solver_type
560 :
561 : INTEGER :: k, n
562 : TYPE(cp_cfm_type) :: matrix_c, matrix_h, matrix_s
563 : TYPE(cp_fm_struct_type), POINTER :: fm_struct_mo, fm_struct_square
564 : TYPE(cp_fm_type) :: matrix_c_im_fm, matrix_c_re_fm, &
565 : matrix_h_im_fm, matrix_h_re_fm, &
566 : matrix_s_im_fm, matrix_s_re_fm
567 :
568 44 : NULLIFY (fm_struct_mo, fm_struct_square)
569 0 : CPASSERT(ASSOCIATED(matrix_c_re))
570 44 : CPASSERT(ASSOCIATED(matrix_c_im))
571 44 : CPASSERT(ASSOCIATED(matrix_h_re))
572 44 : CPASSERT(ASSOCIATED(matrix_h_im))
573 44 : CPASSERT(ASSOCIATED(matrix_s_re))
574 44 : CPASSERT(ASSOCIATED(matrix_s_im))
575 44 : IF (solver_type /= ot_precond_solver_default .AND. &
576 : solver_type /= ot_precond_solver_inv_chol) THEN
577 0 : CPABORT("Complex FULL_SINGLE_INVERSE supports only PRECOND_SOLVER DEFAULT/INVERSE_CHOLESKY")
578 : END IF
579 44 : CALL dbcsr_get_info(matrix_c_re, nfullrows_total=n, nfullcols_total=k)
580 44 : CPASSERT(n > 0 .AND. k > 0 .AND. k <= n)
581 :
582 : CALL cp_fm_struct_create(fm_struct_mo, nrow_global=n, ncol_global=k, &
583 : context=preconditioner_env%ctxt, &
584 44 : para_env=preconditioner_env%para_env)
585 : CALL cp_fm_struct_create(fm_struct_square, nrow_global=n, ncol_global=n, &
586 : context=preconditioner_env%ctxt, &
587 44 : para_env=preconditioner_env%para_env)
588 44 : CALL cp_fm_create(matrix_c_re_fm, fm_struct_mo, name='complex preconditioner C real')
589 44 : CALL cp_fm_create(matrix_c_im_fm, fm_struct_mo, name='complex preconditioner C imaginary')
590 44 : CALL cp_fm_create(matrix_h_re_fm, fm_struct_square, name='complex preconditioner H real')
591 44 : CALL cp_fm_create(matrix_h_im_fm, fm_struct_square, name='complex preconditioner H imaginary')
592 44 : CALL cp_fm_create(matrix_s_re_fm, fm_struct_square, name='complex preconditioner S real')
593 44 : CALL cp_fm_create(matrix_s_im_fm, fm_struct_square, name='complex preconditioner S imaginary')
594 44 : CALL cp_fm_struct_release(fm_struct_mo)
595 44 : CALL cp_fm_struct_release(fm_struct_square)
596 44 : CALL copy_dbcsr_to_fm(matrix_c_re, matrix_c_re_fm)
597 44 : CALL copy_dbcsr_to_fm(matrix_c_im, matrix_c_im_fm)
598 44 : CALL copy_dbcsr_to_fm(matrix_h_re, matrix_h_re_fm)
599 44 : CALL copy_dbcsr_to_fm(matrix_h_im, matrix_h_im_fm)
600 44 : CALL copy_dbcsr_to_fm(matrix_s_re, matrix_s_re_fm)
601 44 : CALL copy_dbcsr_to_fm(matrix_s_im, matrix_s_im_fm)
602 :
603 44 : CALL cp_cfm_create(matrix_c, matrix_c_re_fm%matrix_struct, name='complex preconditioner C')
604 44 : CALL cp_cfm_create(matrix_h, matrix_h_re_fm%matrix_struct, name='complex preconditioner H')
605 44 : CALL cp_cfm_create(matrix_s, matrix_s_re_fm%matrix_struct, name='complex preconditioner S')
606 44 : CALL cp_fm_to_cfm(matrix_c_re_fm, matrix_c_im_fm, matrix_c)
607 44 : CALL cp_fm_to_cfm(matrix_h_re_fm, matrix_h_im_fm, matrix_h)
608 44 : CALL cp_fm_to_cfm(matrix_s_re_fm, matrix_s_im_fm, matrix_s)
609 : CALL make_complex_full_single_inverse(preconditioner_env, matrix_c, matrix_h, matrix_s, &
610 44 : energy_gap)
611 44 : preconditioner_env%solver = solver_type
612 :
613 44 : CALL cp_cfm_release(matrix_s)
614 44 : CALL cp_cfm_release(matrix_h)
615 44 : CALL cp_cfm_release(matrix_c)
616 44 : CALL cp_fm_release(matrix_s_im_fm)
617 44 : CALL cp_fm_release(matrix_s_re_fm)
618 44 : CALL cp_fm_release(matrix_h_im_fm)
619 44 : CALL cp_fm_release(matrix_h_re_fm)
620 44 : CALL cp_fm_release(matrix_c_im_fm)
621 44 : CALL cp_fm_release(matrix_c_re_fm)
622 :
623 44 : END SUBROUTINE make_preconditioner_complex_full_single_inverse
624 :
625 : ! **************************************************************************************************
626 : !> \brief Construct a complex FULL_S_INVERSE preconditioner.
627 : !> \param preconditioner_env preconditioner storage
628 : !> \param matrix_s_re real part of S(k)
629 : !> \param matrix_s_im imaginary part of S(k)
630 : !> \param solver_type requested inverse solver
631 : ! **************************************************************************************************
632 180 : SUBROUTINE make_preconditioner_complex_full_s_inverse(preconditioner_env, &
633 : matrix_s_re, matrix_s_im, solver_type)
634 :
635 : TYPE(preconditioner_type) :: preconditioner_env
636 : TYPE(dbcsr_type), POINTER :: matrix_s_re, matrix_s_im
637 : INTEGER, INTENT(IN) :: solver_type
638 :
639 : TYPE(cp_cfm_type) :: matrix_s
640 :
641 90 : IF (solver_type /= ot_precond_solver_default .AND. &
642 : solver_type /= ot_precond_solver_inv_chol) THEN
643 0 : CPABORT("Complex FULL_S_INVERSE supports only PRECOND_SOLVER DEFAULT/INVERSE_CHOLESKY")
644 : END IF
645 : CALL dbcsr_pair_to_cfm(matrix_s_re, matrix_s_im, preconditioner_env, &
646 90 : 'complex preconditioner S', matrix_s)
647 90 : CALL make_complex_full_s_inverse(preconditioner_env, matrix_s)
648 90 : preconditioner_env%solver = solver_type
649 90 : CALL cp_cfm_release(matrix_s)
650 :
651 90 : END SUBROUTINE make_preconditioner_complex_full_s_inverse
652 :
653 : ! **************************************************************************************************
654 : !> \brief Construct a complex FULL_KINETIC preconditioner.
655 : !> \param preconditioner_env preconditioner storage
656 : !> \param matrix_t_re real part of T(k)
657 : !> \param matrix_t_im imaginary part of T(k)
658 : !> \param matrix_s_re real part of S(k)
659 : !> \param matrix_s_im imaginary part of S(k)
660 : !> \param energy_gap non-negative overlap shift
661 : !> \param solver_type requested inverse solver
662 : ! **************************************************************************************************
663 366 : SUBROUTINE make_preconditioner_complex_full_kinetic(preconditioner_env, &
664 : matrix_t_re, matrix_t_im, &
665 : matrix_s_re, matrix_s_im, &
666 : energy_gap, solver_type)
667 :
668 : TYPE(preconditioner_type) :: preconditioner_env
669 : TYPE(dbcsr_type), POINTER :: matrix_t_re, matrix_t_im, matrix_s_re, &
670 : matrix_s_im
671 : REAL(KIND=dp), INTENT(IN) :: energy_gap
672 : INTEGER, INTENT(IN) :: solver_type
673 :
674 : TYPE(cp_cfm_type) :: matrix_s, matrix_t
675 :
676 122 : IF (solver_type /= ot_precond_solver_default .AND. &
677 : solver_type /= ot_precond_solver_inv_chol) THEN
678 0 : CPABORT("Complex FULL_KINETIC supports only PRECOND_SOLVER DEFAULT/INVERSE_CHOLESKY")
679 : END IF
680 : CALL dbcsr_pair_to_cfm(matrix_t_re, matrix_t_im, preconditioner_env, &
681 122 : 'complex preconditioner T', matrix_t)
682 : CALL dbcsr_pair_to_cfm(matrix_s_re, matrix_s_im, preconditioner_env, &
683 122 : 'complex preconditioner S', matrix_s)
684 122 : CALL make_complex_full_kinetic(preconditioner_env, matrix_t, matrix_s, energy_gap)
685 122 : preconditioner_env%solver = solver_type
686 122 : CALL cp_cfm_release(matrix_s)
687 122 : CALL cp_cfm_release(matrix_t)
688 :
689 122 : END SUBROUTINE make_preconditioner_complex_full_kinetic
690 :
691 : ! **************************************************************************************************
692 : !> \brief Copy a real/imaginary DBCSR pair to one distributed complex full matrix.
693 : !> \param matrix_re real matrix component
694 : !> \param matrix_im imaginary matrix component
695 : !> \param preconditioner_env source of the distribution context
696 : !> \param matrix_name matrix label
697 : !> \param matrix complex output matrix
698 : ! **************************************************************************************************
699 2454 : SUBROUTINE dbcsr_pair_to_cfm(matrix_re, matrix_im, preconditioner_env, matrix_name, matrix)
700 :
701 : TYPE(dbcsr_type), POINTER :: matrix_re, matrix_im
702 : TYPE(preconditioner_type) :: preconditioner_env
703 : CHARACTER(LEN=*), INTENT(IN) :: matrix_name
704 : TYPE(cp_cfm_type), INTENT(OUT) :: matrix
705 :
706 : INTEGER :: ncol, nrow
707 : TYPE(cp_fm_struct_type), POINTER :: fm_struct
708 : TYPE(cp_fm_type) :: matrix_im_fm, matrix_re_fm
709 :
710 818 : NULLIFY (fm_struct)
711 0 : CPASSERT(ASSOCIATED(matrix_re))
712 818 : CPASSERT(ASSOCIATED(matrix_im))
713 818 : CALL dbcsr_get_info(matrix_re, nfullrows_total=nrow, nfullcols_total=ncol)
714 : CALL cp_fm_struct_create(fm_struct, nrow_global=nrow, ncol_global=ncol, &
715 : context=preconditioner_env%ctxt, &
716 818 : para_env=preconditioner_env%para_env)
717 818 : CALL cp_fm_create(matrix_re_fm, fm_struct, name=TRIM(matrix_name)//' real')
718 818 : CALL cp_fm_create(matrix_im_fm, fm_struct, name=TRIM(matrix_name)//' imaginary')
719 818 : CALL cp_fm_struct_release(fm_struct)
720 818 : CALL copy_dbcsr_to_fm(matrix_re, matrix_re_fm)
721 818 : CALL copy_dbcsr_to_fm(matrix_im, matrix_im_fm)
722 818 : CALL cp_cfm_create(matrix, matrix_re_fm%matrix_struct, name=matrix_name)
723 818 : CALL cp_fm_to_cfm(matrix_re_fm, matrix_im_fm, matrix)
724 818 : CALL cp_fm_release(matrix_im_fm)
725 818 : CALL cp_fm_release(matrix_re_fm)
726 :
727 818 : END SUBROUTINE dbcsr_pair_to_cfm
728 :
729 : ! **************************************************************************************************
730 : !> \brief Allows for a restart of the preconditioner
731 : !> depending on the method it purges all arrays or keeps them
732 : !> \param qs_env ...
733 : !> \param preconditioner ...
734 : !> \param prec_type ...
735 : !> \param nspins ...
736 : ! **************************************************************************************************
737 7481 : SUBROUTINE restart_preconditioner(qs_env, preconditioner, prec_type, nspins)
738 :
739 : TYPE(qs_environment_type), POINTER :: qs_env
740 : TYPE(preconditioner_p_type), DIMENSION(:), POINTER :: preconditioner
741 : INTEGER, INTENT(IN) :: prec_type, nspins
742 :
743 : INTEGER :: ispin
744 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
745 : TYPE(mp_para_env_type), POINTER :: para_env
746 :
747 7481 : NULLIFY (para_env, blacs_env)
748 7481 : CALL get_qs_env(qs_env, para_env=para_env, blacs_env=blacs_env)
749 :
750 7481 : IF (ASSOCIATED(preconditioner)) THEN
751 6668 : SELECT CASE (prec_type)
752 : CASE (ot_precond_full_all, ot_precond_fermi_low_rank, &
753 : ot_precond_full_single) ! these depend on the ks matrix
754 3460 : DO ispin = 1, SIZE(preconditioner)
755 1930 : CALL destroy_preconditioner(preconditioner(ispin)%preconditioner)
756 3460 : DEALLOCATE (preconditioner(ispin)%preconditioner)
757 : END DO
758 1530 : DEALLOCATE (preconditioner)
759 : CASE (ot_precond_none, ot_precond_full_kinetic, ot_precond_s_inverse, &
760 : ot_precond_full_single_inverse) ! these are 'independent'
761 : ! do nothing
762 : CASE DEFAULT
763 5138 : CPABORT("Unknown preconditioner type")
764 : END SELECT
765 : END IF
766 :
767 : ! add an OT preconditioner if none is present
768 7481 : IF (.NOT. ASSOCIATED(preconditioner)) THEN
769 7105 : SELECT CASE (prec_type)
770 : CASE (ot_precond_full_all, ot_precond_fermi_low_rank, ot_precond_full_single_inverse)
771 13845 : ALLOCATE (preconditioner(nspins))
772 : CASE DEFAULT
773 4514 : ALLOCATE (preconditioner(1))
774 : END SELECT
775 8663 : DO ispin = 1, SIZE(preconditioner)
776 23950 : ALLOCATE (preconditioner(ispin)%preconditioner)
777 : CALL init_preconditioner(preconditioner(ispin)%preconditioner, &
778 : para_env=para_env, &
779 8663 : blacs_env=blacs_env)
780 : END DO
781 : END IF
782 :
783 7481 : END SUBROUTINE restart_preconditioner
784 :
785 : ! **************************************************************************************************
786 : !> \brief ...
787 : !> \param qs_env ...
788 : !> \param mos ...
789 : !> \param matrix_ks ...
790 : !> \param matrix_s ...
791 : !> \param ot_preconditioner ...
792 : !> \param prec_type ...
793 : !> \param solver_type ...
794 : !> \param energy_gap ...
795 : !> \param nspins ...
796 : !> \param has_unit_metric ...
797 : !> \param convert_to_dbcsr ...
798 : !> \param chol_type ...
799 : !> \param full_mo_set ...
800 : !> \param chebyshev_degree ...
801 : !> \param low_rank_base ...
802 : ! **************************************************************************************************
803 7481 : SUBROUTINE prepare_preconditioner(qs_env, mos, matrix_ks, matrix_s, &
804 : ot_preconditioner, prec_type, solver_type, &
805 : energy_gap, nspins, has_unit_metric, &
806 : convert_to_dbcsr, chol_type, full_mo_set, chebyshev_degree, &
807 : low_rank_base)
808 :
809 : TYPE(qs_environment_type), POINTER :: qs_env
810 : TYPE(mo_set_type), DIMENSION(:), INTENT(INOUT) :: mos
811 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s
812 : TYPE(preconditioner_p_type), DIMENSION(:), POINTER :: ot_preconditioner
813 : INTEGER, INTENT(IN) :: prec_type, solver_type
814 : REAL(dp), INTENT(IN) :: energy_gap
815 : INTEGER, INTENT(IN) :: nspins
816 : LOGICAL, INTENT(IN), OPTIONAL :: has_unit_metric, convert_to_dbcsr
817 : INTEGER, INTENT(IN), OPTIONAL :: chol_type
818 : LOGICAL, INTENT(IN), OPTIONAL :: full_mo_set
819 : INTEGER, INTENT(IN), OPTIONAL :: chebyshev_degree, low_rank_base
820 :
821 : CHARACTER(LEN=*), PARAMETER :: routineN = 'prepare_preconditioner'
822 :
823 : CHARACTER(LEN=default_string_length) :: msg
824 : INTEGER :: handle, icall, ispin, &
825 : my_chebyshev_degree, my_low_rank_base, &
826 : n_loops
827 : INTEGER, DIMENSION(5) :: nocc, norb
828 7481 : INTEGER, DIMENSION(:), POINTER :: multiple_unit_cell
829 : LOGICAL :: do_co_rotate, do_kpoints, &
830 : my_convert_to_dbcsr, my_full_mo_set, &
831 : my_has_unit_metric, use_mo_coeff_b
832 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
833 : TYPE(cp_fm_type), POINTER :: mo_coeff
834 7481 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: kinetic
835 : TYPE(dbcsr_type), POINTER :: matrix_t, mo_coeff_b
836 : TYPE(dft_control_type), POINTER :: dft_control
837 : TYPE(mp_para_env_type), POINTER :: para_env
838 : TYPE(section_vals_type), POINTER :: input, subsys_section
839 :
840 7481 : CALL timeset(routineN, handle)
841 7481 : NULLIFY (matrix_t, mo_coeff_b, mo_coeff, kinetic, dft_control, para_env, blacs_env, &
842 7481 : input, subsys_section, multiple_unit_cell)
843 7481 : my_has_unit_metric = .FALSE.
844 7481 : IF (PRESENT(has_unit_metric)) my_has_unit_metric = has_unit_metric
845 7481 : my_convert_to_dbcsr = .TRUE.
846 7481 : IF (PRESENT(convert_to_dbcsr)) my_convert_to_dbcsr = convert_to_dbcsr
847 7481 : my_full_mo_set = .FALSE.
848 7481 : IF (PRESENT(full_mo_set)) my_full_mo_set = full_mo_set
849 7481 : my_chebyshev_degree = 8
850 7481 : IF (PRESENT(chebyshev_degree)) my_chebyshev_degree = chebyshev_degree
851 7481 : my_low_rank_base = ot_low_rank_base_overlap
852 7481 : IF (PRESENT(low_rank_base)) my_low_rank_base = low_rank_base
853 7481 : do_kpoints = .FALSE.
854 :
855 : CALL get_qs_env(qs_env, &
856 : dft_control=dft_control, &
857 : do_kpoints=do_kpoints, &
858 : input=input, &
859 : para_env=para_env, &
860 7481 : blacs_env=blacs_env)
861 :
862 7481 : IF (my_low_rank_base == ot_low_rank_base_lattice_fft) THEN
863 2 : IF (prec_type /= ot_precond_fermi_low_rank) THEN
864 0 : CPABORT("FERMI_LOW_RANK_BASE LATTICE_FFT requires PRECONDITIONER FERMI_LOW_RANK")
865 : END IF
866 2 : IF (my_has_unit_metric) THEN
867 0 : CPABORT("FERMI_LOW_RANK_BASE LATTICE_FFT is unavailable with a unit overlap metric")
868 : END IF
869 2 : IF (do_kpoints) THEN
870 0 : CPABORT("FERMI_LOW_RANK_BASE LATTICE_FFT currently supports only real Gamma-point OT")
871 : END IF
872 2 : IF (.NOT. ASSOCIATED(input)) THEN
873 0 : CPABORT("FERMI_LOW_RANK_BASE LATTICE_FFT could not access the input tree")
874 : END IF
875 2 : subsys_section => section_vals_get_subs_vals(input, "SUBSYS")
876 2 : IF (.NOT. ASSOCIATED(subsys_section)) THEN
877 0 : CPABORT("FERMI_LOW_RANK_BASE LATTICE_FFT requires a SUBSYS section")
878 : END IF
879 : CALL section_vals_val_get(subsys_section, "CELL%MULTIPLE_UNIT_CELL", &
880 2 : i_vals=multiple_unit_cell)
881 2 : IF (.NOT. ASSOCIATED(multiple_unit_cell)) THEN
882 0 : CPABORT("FERMI_LOW_RANK_BASE LATTICE_FFT could not read MULTIPLE_UNIT_CELL")
883 : END IF
884 16 : IF (SIZE(multiple_unit_cell) /= 3 .OR. ANY(multiple_unit_cell <= 0) .OR. &
885 : PRODUCT(multiple_unit_cell) <= 1) THEN
886 0 : CPABORT("FERMI_LOW_RANK_BASE LATTICE_FFT requires an explicit replicated multiple unit cell")
887 : END IF
888 2 : my_convert_to_dbcsr = .FALSE.
889 7479 : ELSE IF (my_low_rank_base /= ot_low_rank_base_overlap) THEN
890 0 : CPABORT("Unknown FERMI_LOW_RANK_BASE setting")
891 : END IF
892 :
893 7481 : IF (dft_control%qs_control%semi_empirical .OR. dft_control%qs_control%dftb .OR. &
894 : dft_control%qs_control%xtb) THEN
895 2148 : IF (prec_type == ot_precond_full_kinetic) THEN
896 0 : msg = "Full_kinetic not available for semi-empirical methods"
897 0 : CPABORT(TRIM(msg))
898 : END IF
899 2148 : matrix_t => matrix_s(1)%matrix
900 : ELSE
901 5333 : CPASSERT(.NOT. my_has_unit_metric)
902 5333 : CALL get_qs_env(qs_env, kinetic=kinetic)
903 5333 : matrix_t => kinetic(1)%matrix
904 : END IF
905 :
906 : ! use full set of MOs or just occupied MOs
907 7481 : nocc = 0
908 7481 : norb = 0
909 7481 : IF (my_full_mo_set) THEN
910 38 : DO ispin = 1, nspins
911 20 : CALL get_mo_set(mo_set=mos(ispin), homo=nocc(ispin), nmo=norb(ispin))
912 38 : CALL set_mo_set(mo_set=mos(ispin), homo=norb(ispin))
913 : END DO
914 : END IF
915 : !determines how often make preconditioner is called, spin dependent methods have to be called twice
916 7481 : n_loops = 1
917 7481 : IF (prec_type == ot_precond_full_single_inverse) n_loops = nspins
918 : ! check whether we need the ev and rotate the MOs
919 2556 : SELECT CASE (prec_type)
920 : CASE (ot_precond_full_all)
921 : ! FULL_ALL needs canonical occupied eigenstates and one preconditioner per spin.
922 2556 : n_loops = nspins
923 :
924 2556 : do_co_rotate = ASSOCIATED(qs_env%mo_derivs)
925 5906 : DO ispin = 1, nspins
926 3350 : CALL get_mo_set(mo_set=mos(ispin), mo_coeff_b=mo_coeff_b, mo_coeff=mo_coeff)
927 3350 : use_mo_coeff_b = mos(ispin)%use_mo_coeff_b
928 5906 : IF (use_mo_coeff_b .AND. do_co_rotate) THEN
929 : CALL calculate_subspace_eigenvalues(mo_coeff_b, matrix_ks(ispin)%matrix, &
930 : do_rotation=.TRUE., &
931 : co_rotate=qs_env%mo_derivs(ispin)%matrix, &
932 : para_env=para_env, &
933 3330 : blacs_env=blacs_env)
934 20 : ELSE IF (use_mo_coeff_b) THEN
935 : CALL calculate_subspace_eigenvalues(mo_coeff_b, matrix_ks(ispin)%matrix, &
936 : do_rotation=.TRUE., &
937 : para_env=para_env, &
938 20 : blacs_env=blacs_env)
939 : ELSE
940 : CALL calculate_subspace_eigenvalues(mo_coeff, matrix_ks(ispin)%matrix, &
941 0 : do_rotation=.TRUE.)
942 : END IF
943 : END DO
944 : CASE (ot_precond_fermi_low_rank)
945 : ! The common-reference low-rank operator depends only on the occupied projector.
946 : ! Keep the input gauge unchanged so that OT ROTATION remains well defined.
947 7481 : n_loops = nspins
948 : CASE DEFAULT
949 : ! No need to rotate the MOs
950 : END SELECT
951 :
952 : ! check whether we have a preconditioner
953 724 : SELECT CASE (prec_type)
954 : CASE (ot_precond_none)
955 1448 : DO ispin = 1, SIZE(ot_preconditioner)
956 1448 : ot_preconditioner(ispin)%preconditioner%in_use = 0
957 : END DO
958 : CASE DEFAULT
959 22140 : DO icall = 1, n_loops
960 7902 : ot_preconditioner(icall)%preconditioner%polynomial_degree = my_chebyshev_degree
961 7902 : ot_preconditioner(icall)%preconditioner%low_rank_base = my_low_rank_base
962 7902 : IF (my_low_rank_base == ot_low_rank_base_lattice_fft) THEN
963 16 : ot_preconditioner(icall)%preconditioner%lattice_dims = multiple_unit_cell
964 : END IF
965 14659 : IF (my_has_unit_metric) THEN
966 : CALL make_preconditioner(ot_preconditioner(icall)%preconditioner, &
967 : prec_type, &
968 : solver_type, &
969 : matrix_h=matrix_ks(icall)%matrix, &
970 : mo_set=mos(icall), &
971 : energy_gap=energy_gap, &
972 474 : convert_precond_to_dbcsr=my_convert_to_dbcsr)
973 : ELSE
974 : CALL make_preconditioner(ot_preconditioner(icall)%preconditioner, &
975 : prec_type, &
976 : solver_type, &
977 : matrix_h=matrix_ks(icall)%matrix, &
978 : matrix_s=matrix_s(1)%matrix, &
979 : matrix_t=matrix_t, &
980 : mo_set=mos(icall), &
981 : energy_gap=energy_gap, &
982 7428 : convert_precond_to_dbcsr=my_convert_to_dbcsr, chol_type=chol_type)
983 : END IF
984 : END DO
985 : END SELECT
986 :
987 : ! reset homo values
988 7481 : IF (my_full_mo_set) THEN
989 38 : DO ispin = 1, nspins
990 38 : CALL set_mo_set(mo_set=mos(ispin), homo=nocc(ispin))
991 : END DO
992 : END IF
993 :
994 7481 : CALL timestop(handle)
995 :
996 7481 : END SUBROUTINE prepare_preconditioner
997 :
998 : END MODULE preconditioner
|
