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_control_types, ONLY: dft_control_type
18 : USE cp_dbcsr_api, ONLY: dbcsr_p_type,&
19 : dbcsr_type
20 : USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm
21 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
22 : cp_fm_struct_release,&
23 : cp_fm_struct_type
24 : USE cp_fm_types, ONLY: cp_fm_create,&
25 : cp_fm_get_info,&
26 : cp_fm_release,&
27 : cp_fm_to_fm,&
28 : cp_fm_type
29 : USE input_constants, ONLY: cholesky_reduce,&
30 : ot_precond_full_all,&
31 : ot_precond_full_kinetic,&
32 : ot_precond_full_single,&
33 : ot_precond_full_single_inverse,&
34 : ot_precond_none,&
35 : ot_precond_s_inverse,&
36 : ot_precond_solver_update
37 : USE kinds, ONLY: default_string_length,&
38 : dp
39 : USE message_passing, ONLY: mp_para_env_type
40 : USE preconditioner_apply, ONLY: apply_preconditioner_dbcsr,&
41 : apply_preconditioner_fm
42 : USE preconditioner_makes, ONLY: make_preconditioner_matrix
43 : USE preconditioner_solvers, ONLY: solve_preconditioner,&
44 : transfer_dbcsr_to_fm,&
45 : transfer_fm_to_dbcsr
46 : USE preconditioner_types, ONLY: destroy_preconditioner,&
47 : init_preconditioner,&
48 : preconditioner_p_type,&
49 : preconditioner_type
50 : USE qs_environment_types, ONLY: get_qs_env,&
51 : qs_environment_type
52 : USE qs_mo_methods, ONLY: calculate_subspace_eigenvalues
53 : USE qs_mo_types, ONLY: get_mo_set,&
54 : mo_set_type,&
55 : set_mo_set
56 : #include "./base/base_uses.f90"
57 :
58 : IMPLICIT NONE
59 :
60 : PRIVATE
61 :
62 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'preconditioner'
63 :
64 : PUBLIC :: make_preconditioner, restart_preconditioner
65 : PUBLIC :: apply_preconditioner, prepare_preconditioner
66 :
67 : ! The public interface for apply preconditioner, the routines can be found in preconditioner_apply.F
68 : INTERFACE apply_preconditioner
69 : MODULE PROCEDURE apply_preconditioner_dbcsr
70 : MODULE PROCEDURE apply_preconditioner_fm
71 : END INTERFACE
72 :
73 : ! **************************************************************************************************
74 :
75 : CONTAINS
76 :
77 : ! **************************************************************************************************
78 :
79 : ! creates a preconditioner for the system (H-energy_homo S)
80 : ! this preconditioner is (must be) symmetric positive definite.
81 : ! currently uses a atom-block-diagonal form
82 : ! each block will be ....
83 : ! might overwrite matrix_h, matrix_t
84 :
85 : ! **************************************************************************************************
86 : !> \brief ...
87 : !> \param preconditioner_env ...
88 : !> \param precon_type ...
89 : !> \param solver_type ...
90 : !> \param matrix_h ...
91 : !> \param matrix_s ...
92 : !> \param matrix_t ...
93 : !> \param mo_set ...
94 : !> \param energy_gap ...
95 : !> \param convert_precond_to_dbcsr ...
96 : !> \param chol_type ...
97 : !> \par History
98 : !> 09.2014 removed some unused or unfinished methods
99 : !> removed sparse preconditioners and the
100 : !> sparse approximate inverse at rev 14341 [Florian Schiffmann]
101 : ! **************************************************************************************************
102 9466 : SUBROUTINE make_preconditioner(preconditioner_env, precon_type, solver_type, matrix_h, matrix_s, &
103 : matrix_t, mo_set, energy_gap, convert_precond_to_dbcsr, chol_type)
104 :
105 : TYPE(preconditioner_type) :: preconditioner_env
106 : INTEGER, INTENT(IN) :: precon_type, solver_type
107 : TYPE(dbcsr_type), POINTER :: matrix_h
108 : TYPE(dbcsr_type), OPTIONAL, POINTER :: matrix_s, matrix_t
109 : TYPE(mo_set_type), INTENT(IN) :: mo_set
110 : REAL(KIND=dp) :: energy_gap
111 : LOGICAL, INTENT(IN), OPTIONAL :: convert_precond_to_dbcsr
112 : INTEGER, INTENT(IN), OPTIONAL :: chol_type
113 :
114 : CHARACTER(len=*), PARAMETER :: routineN = 'make_preconditioner'
115 :
116 : INTEGER :: handle, k, my_solver_type, nao, nhomo
117 : LOGICAL :: my_convert_precond_to_dbcsr, &
118 : needs_full_spectrum, needs_homo, &
119 : use_mo_coeff_b
120 : REAL(KIND=dp) :: energy_homo
121 9466 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eigenvalues_ot
122 : TYPE(cp_fm_struct_type), POINTER :: fm_struct
123 : TYPE(cp_fm_type) :: mo_occ
124 : TYPE(cp_fm_type), POINTER :: mo_coeff
125 : TYPE(dbcsr_type), POINTER :: mo_coeff_b
126 :
127 9466 : CALL timeset(routineN, handle)
128 :
129 9466 : CALL get_mo_set(mo_set=mo_set, mo_coeff=mo_coeff, mo_coeff_b=mo_coeff_b, homo=nhomo)
130 9466 : use_mo_coeff_b = mo_set%use_mo_coeff_b
131 9466 : CALL cp_fm_get_info(mo_coeff, ncol_global=k, nrow_global=nao)
132 :
133 : ! Starting some matrix mess, check where to store the result in preconditioner_env, fm or dbcsr_matrix
134 9466 : my_convert_precond_to_dbcsr = .FALSE.
135 9466 : IF (PRESENT(convert_precond_to_dbcsr)) my_convert_precond_to_dbcsr = convert_precond_to_dbcsr
136 :
137 : ! Thanks to the mess with the matrices we need to make sure in this case that the
138 : ! Previous inverse is properly stored as a sparse matrix, fm gets deallocated here
139 : ! if it wasn't anyway
140 9466 : IF (preconditioner_env%solver == ot_precond_solver_update) THEN
141 4 : CALL transfer_fm_to_dbcsr(preconditioner_env%fm, preconditioner_env%dbcsr_matrix, matrix_h)
142 : END IF
143 :
144 9466 : needs_full_spectrum = .FALSE.
145 9466 : needs_homo = .FALSE.
146 :
147 13114 : SELECT CASE (precon_type)
148 : CASE (ot_precond_full_all)
149 3648 : needs_full_spectrum = .TRUE.
150 : ! both of them need the coefficients as fm's, more matrix mess
151 3648 : IF (use_mo_coeff_b) THEN
152 3362 : CALL copy_dbcsr_to_fm(mo_coeff_b, mo_coeff)
153 : END IF
154 : CASE (ot_precond_full_single)
155 38 : needs_homo = .TRUE.
156 : ! XXXX to be removed if homo estimate only is implemented
157 38 : needs_full_spectrum = .TRUE.
158 : CASE (ot_precond_full_kinetic, ot_precond_s_inverse, ot_precond_full_single_inverse)
159 : ! these should be happy without an estimate for the homo energy
160 : ! preconditioning can not depend on an absolute eigenvalue, only on eigenvalue differences
161 : CASE DEFAULT
162 9466 : CPABORT("The preconditioner is unknown ...")
163 : END SELECT
164 :
165 28282 : ALLOCATE (eigenvalues_ot(k))
166 9466 : energy_homo = 0.0_dp
167 9466 : IF (needs_full_spectrum) THEN
168 : ! XXXXXXXXXXXXXXXX do not touch the initial MOs, could be harmful for either
169 : ! the case of non-equivalent MOs but also for the derivate
170 : ! we could already have all eigenvalues e.g. full_all and we could skip this
171 : ! to be optimised later.
172 : ! one flaw is that not all SCF methods (i.e. that go over mo_derivs directly)
173 : ! have a 'valid' matrix_h... (we even don't know what evals are in that case)
174 3686 : IF (use_mo_coeff_b) THEN
175 : CALL calculate_subspace_eigenvalues(mo_coeff_b, matrix_h, &
176 : eigenvalues_ot, do_rotation=.FALSE., &
177 : para_env=mo_coeff%matrix_struct%para_env, &
178 3392 : blacs_env=mo_coeff%matrix_struct%context)
179 : ELSE
180 : CALL calculate_subspace_eigenvalues(mo_coeff, matrix_h, &
181 294 : eigenvalues_ot, do_rotation=.FALSE.)
182 : END IF
183 3686 : IF (k > 0) THEN
184 3580 : CPASSERT(nhomo > 0 .AND. nhomo <= k)
185 3580 : energy_homo = eigenvalues_ot(nhomo)
186 : END IF
187 : ELSE
188 5780 : IF (needs_homo) THEN
189 0 : CPABORT("Not yet implemented")
190 : END IF
191 : END IF
192 :
193 : ! After all bits and pieces of checking and initialization, here comes the
194 : ! part where the preconditioner matrix gets created and solved.
195 : ! This will give the matrices for later use
196 9466 : my_solver_type = solver_type
197 9466 : preconditioner_env%in_use = precon_type
198 9466 : preconditioner_env%cholesky_use = cholesky_reduce
199 9466 : IF (PRESENT(chol_type)) preconditioner_env%cholesky_use = chol_type
200 : preconditioner_env%in_use = precon_type
201 9466 : IF (nhomo == k) THEN
202 : CALL make_preconditioner_matrix(preconditioner_env, matrix_h, matrix_s, matrix_t, mo_coeff, &
203 9382 : energy_homo, eigenvalues_ot, energy_gap, my_solver_type)
204 : ELSE
205 : CALL cp_fm_struct_create(fm_struct, nrow_global=nao, ncol_global=nhomo, &
206 : context=preconditioner_env%ctxt, &
207 84 : para_env=preconditioner_env%para_env)
208 84 : CALL cp_fm_create(mo_occ, fm_struct)
209 84 : CALL cp_fm_to_fm(mo_coeff, mo_occ, nhomo)
210 84 : CALL cp_fm_struct_release(fm_struct)
211 : !
212 : CALL make_preconditioner_matrix(preconditioner_env, matrix_h, matrix_s, matrix_t, mo_occ, &
213 84 : energy_homo, eigenvalues_ot(1:nhomo), energy_gap, my_solver_type)
214 : !
215 84 : CALL cp_fm_release(mo_occ)
216 : END IF
217 :
218 9466 : CALL solve_preconditioner(my_solver_type, preconditioner_env, matrix_s, matrix_h)
219 :
220 : ! Here comes more matrix mess, make sure to output the correct matrix format,
221 : ! A bit pointless to convert the cholesky factorized version as it doesn't work in
222 : ! dbcsr form and will crash later,...
223 9466 : IF (my_convert_precond_to_dbcsr) THEN
224 7808 : CALL transfer_fm_to_dbcsr(preconditioner_env%fm, preconditioner_env%dbcsr_matrix, matrix_h)
225 : ELSE
226 : CALL transfer_dbcsr_to_fm(preconditioner_env%dbcsr_matrix, preconditioner_env%fm, &
227 1658 : preconditioner_env%para_env, preconditioner_env%ctxt)
228 : END IF
229 :
230 9466 : DEALLOCATE (eigenvalues_ot)
231 :
232 9466 : CALL timestop(handle)
233 :
234 9466 : END SUBROUTINE make_preconditioner
235 :
236 : ! **************************************************************************************************
237 : !> \brief Allows for a restart of the preconditioner
238 : !> depending on the method it purges all arrays or keeps them
239 : !> \param qs_env ...
240 : !> \param preconditioner ...
241 : !> \param prec_type ...
242 : !> \param nspins ...
243 : ! **************************************************************************************************
244 7411 : SUBROUTINE restart_preconditioner(qs_env, preconditioner, prec_type, nspins)
245 :
246 : TYPE(qs_environment_type), POINTER :: qs_env
247 : TYPE(preconditioner_p_type), DIMENSION(:), POINTER :: preconditioner
248 : INTEGER, INTENT(IN) :: prec_type, nspins
249 :
250 : INTEGER :: ispin
251 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
252 : TYPE(mp_para_env_type), POINTER :: para_env
253 :
254 7411 : NULLIFY (para_env, blacs_env)
255 7411 : CALL get_qs_env(qs_env, para_env=para_env, blacs_env=blacs_env)
256 :
257 7411 : IF (ASSOCIATED(preconditioner)) THEN
258 6636 : SELECT CASE (prec_type)
259 : CASE (ot_precond_full_all, ot_precond_full_single) ! these depend on the ks matrix
260 3384 : DO ispin = 1, SIZE(preconditioner)
261 1888 : CALL destroy_preconditioner(preconditioner(ispin)%preconditioner)
262 3384 : DEALLOCATE (preconditioner(ispin)%preconditioner)
263 : END DO
264 1496 : DEALLOCATE (preconditioner)
265 : CASE (ot_precond_none, ot_precond_full_kinetic, ot_precond_s_inverse, &
266 : ot_precond_full_single_inverse) ! these are 'independent'
267 : ! do nothing
268 : CASE DEFAULT
269 5140 : CPABORT("Unknown preconditioner type")
270 : END SELECT
271 : END IF
272 :
273 : ! add an OT preconditioner if none is present
274 7411 : IF (.NOT. ASSOCIATED(preconditioner)) THEN
275 6913 : SELECT CASE (prec_type)
276 : CASE (ot_precond_full_all, ot_precond_full_single_inverse)
277 13467 : ALLOCATE (preconditioner(nspins))
278 : CASE DEFAULT
279 4388 : ALLOCATE (preconditioner(1))
280 : END SELECT
281 8417 : DO ispin = 1, SIZE(preconditioner)
282 4650 : ALLOCATE (preconditioner(ispin)%preconditioner)
283 : CALL init_preconditioner(preconditioner(ispin)%preconditioner, &
284 : para_env=para_env, &
285 8417 : blacs_env=blacs_env)
286 : END DO
287 : END IF
288 :
289 7411 : END SUBROUTINE restart_preconditioner
290 :
291 : ! **************************************************************************************************
292 : !> \brief ...
293 : !> \param qs_env ...
294 : !> \param mos ...
295 : !> \param matrix_ks ...
296 : !> \param matrix_s ...
297 : !> \param ot_preconditioner ...
298 : !> \param prec_type ...
299 : !> \param solver_type ...
300 : !> \param energy_gap ...
301 : !> \param nspins ...
302 : !> \param has_unit_metric ...
303 : !> \param convert_to_dbcsr ...
304 : !> \param chol_type ...
305 : !> \param full_mo_set ...
306 : ! **************************************************************************************************
307 7411 : SUBROUTINE prepare_preconditioner(qs_env, mos, matrix_ks, matrix_s, &
308 : ot_preconditioner, prec_type, solver_type, &
309 : energy_gap, nspins, has_unit_metric, &
310 : convert_to_dbcsr, chol_type, full_mo_set)
311 :
312 : TYPE(qs_environment_type), POINTER :: qs_env
313 : TYPE(mo_set_type), DIMENSION(:), INTENT(INOUT) :: mos
314 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s
315 : TYPE(preconditioner_p_type), DIMENSION(:), POINTER :: ot_preconditioner
316 : INTEGER, INTENT(IN) :: prec_type, solver_type
317 : REAL(dp), INTENT(IN) :: energy_gap
318 : INTEGER, INTENT(IN) :: nspins
319 : LOGICAL, INTENT(IN), OPTIONAL :: has_unit_metric, convert_to_dbcsr
320 : INTEGER, INTENT(IN), OPTIONAL :: chol_type
321 : LOGICAL, INTENT(IN), OPTIONAL :: full_mo_set
322 :
323 : CHARACTER(LEN=*), PARAMETER :: routineN = 'prepare_preconditioner'
324 :
325 : CHARACTER(LEN=default_string_length) :: msg
326 : INTEGER :: handle, icall, ispin, n_loops
327 : INTEGER, DIMENSION(5) :: nocc, norb
328 : LOGICAL :: do_co_rotate, my_convert_to_dbcsr, &
329 : my_full_mo_set, my_has_unit_metric, &
330 : use_mo_coeff_b
331 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
332 : TYPE(cp_fm_type), POINTER :: mo_coeff
333 7411 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: kinetic
334 : TYPE(dbcsr_type), POINTER :: matrix_t, mo_coeff_b
335 : TYPE(dft_control_type), POINTER :: dft_control
336 : TYPE(mp_para_env_type), POINTER :: para_env
337 :
338 7411 : CALL timeset(routineN, handle)
339 7411 : NULLIFY (matrix_t, mo_coeff_b, mo_coeff, kinetic, dft_control, para_env, blacs_env)
340 7411 : my_has_unit_metric = .FALSE.
341 7411 : IF (PRESENT(has_unit_metric)) my_has_unit_metric = has_unit_metric
342 7411 : my_convert_to_dbcsr = .TRUE.
343 7411 : IF (PRESENT(convert_to_dbcsr)) my_convert_to_dbcsr = convert_to_dbcsr
344 7411 : my_full_mo_set = .FALSE.
345 7411 : IF (PRESENT(full_mo_set)) my_full_mo_set = full_mo_set
346 :
347 : CALL get_qs_env(qs_env, &
348 : dft_control=dft_control, &
349 : para_env=para_env, &
350 7411 : blacs_env=blacs_env)
351 :
352 7411 : IF (dft_control%qs_control%semi_empirical .OR. dft_control%qs_control%dftb .OR. &
353 : dft_control%qs_control%xtb) THEN
354 2136 : IF (prec_type == ot_precond_full_kinetic) THEN
355 0 : msg = "Full_kinetic not available for semi-empirical methods"
356 0 : CPABORT(TRIM(msg))
357 : END IF
358 2136 : matrix_t => matrix_s(1)%matrix
359 : ELSE
360 5275 : CPASSERT(.NOT. my_has_unit_metric)
361 5275 : CALL get_qs_env(qs_env, kinetic=kinetic)
362 5275 : matrix_t => kinetic(1)%matrix
363 : END IF
364 :
365 : ! use full set of MOs or just occupied MOs
366 7411 : nocc = 0
367 7411 : norb = 0
368 7411 : IF (my_full_mo_set) THEN
369 38 : DO ispin = 1, nspins
370 20 : CALL get_mo_set(mo_set=mos(ispin), homo=nocc(ispin), nmo=norb(ispin))
371 38 : CALL set_mo_set(mo_set=mos(ispin), homo=norb(ispin))
372 : END DO
373 : END IF
374 : !determines how often make preconditioner is called, spin dependent methods have to be called twice
375 7411 : n_loops = 1
376 7411 : IF (prec_type == ot_precond_full_single_inverse) n_loops = nspins
377 : ! check whether we need the ev and rotate the MOs
378 2484 : SELECT CASE (prec_type)
379 : CASE (ot_precond_full_all)
380 : ! if one of these preconditioners is used every spin needs to call make_preconditioner
381 2484 : n_loops = nspins
382 :
383 2484 : do_co_rotate = ASSOCIATED(qs_env%mo_derivs)
384 10665 : DO ispin = 1, nspins
385 3254 : CALL get_mo_set(mo_set=mos(ispin), mo_coeff_b=mo_coeff_b, mo_coeff=mo_coeff)
386 3254 : use_mo_coeff_b = mos(ispin)%use_mo_coeff_b
387 5738 : IF (use_mo_coeff_b .AND. do_co_rotate) THEN
388 : CALL calculate_subspace_eigenvalues(mo_coeff_b, matrix_ks(ispin)%matrix, &
389 : do_rotation=.TRUE., &
390 : co_rotate=qs_env%mo_derivs(ispin)%matrix, &
391 : para_env=para_env, &
392 3234 : blacs_env=blacs_env)
393 20 : ELSE IF (use_mo_coeff_b) THEN
394 : CALL calculate_subspace_eigenvalues(mo_coeff_b, matrix_ks(ispin)%matrix, &
395 : do_rotation=.TRUE., &
396 : para_env=para_env, &
397 20 : blacs_env=blacs_env)
398 : ELSE
399 : CALL calculate_subspace_eigenvalues(mo_coeff, matrix_ks(ispin)%matrix, &
400 0 : do_rotation=.TRUE.)
401 : END IF
402 : END DO
403 : CASE DEFAULT
404 : ! No need to rotate the MOs
405 : END SELECT
406 :
407 : ! check whether we have a preconditioner
408 724 : SELECT CASE (prec_type)
409 : CASE (ot_precond_none)
410 1448 : DO ispin = 1, SIZE(ot_preconditioner)
411 1448 : ot_preconditioner(ispin)%preconditioner%in_use = 0
412 : END DO
413 : CASE DEFAULT
414 21834 : DO icall = 1, n_loops
415 14423 : IF (my_has_unit_metric) THEN
416 : CALL make_preconditioner(ot_preconditioner(icall)%preconditioner, &
417 : prec_type, &
418 : solver_type, &
419 : matrix_h=matrix_ks(icall)%matrix, &
420 : mo_set=mos(icall), &
421 : energy_gap=energy_gap, &
422 474 : convert_precond_to_dbcsr=my_convert_to_dbcsr)
423 : ELSE
424 : CALL make_preconditioner(ot_preconditioner(icall)%preconditioner, &
425 : prec_type, &
426 : solver_type, &
427 : matrix_h=matrix_ks(icall)%matrix, &
428 : matrix_s=matrix_s(1)%matrix, &
429 : matrix_t=matrix_t, &
430 : mo_set=mos(icall), &
431 : energy_gap=energy_gap, &
432 7262 : convert_precond_to_dbcsr=my_convert_to_dbcsr, chol_type=chol_type)
433 : END IF
434 : END DO
435 : END SELECT
436 :
437 : ! reset homo values
438 7411 : IF (my_full_mo_set) THEN
439 38 : DO ispin = 1, nspins
440 38 : CALL set_mo_set(mo_set=mos(ispin), homo=nocc(ispin))
441 : END DO
442 : END IF
443 :
444 7411 : CALL timestop(handle)
445 :
446 7411 : END SUBROUTINE prepare_preconditioner
447 :
448 : END MODULE preconditioner
449 :
|
