Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 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 8500 : 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 : 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 8500 : CALL timeset(routineN, handle)
128 :
129 8500 : CALL get_mo_set(mo_set=mo_set, mo_coeff=mo_coeff, mo_coeff_b=mo_coeff_b, homo=nhomo)
130 8500 : use_mo_coeff_b = mo_set%use_mo_coeff_b
131 8500 : 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 8500 : my_convert_precond_to_dbcsr = .FALSE.
135 8500 : 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 8500 : IF (preconditioner_env%solver == ot_precond_solver_update) &
141 4 : CALL transfer_fm_to_dbcsr(preconditioner_env%fm, preconditioner_env%dbcsr_matrix, matrix_h)
142 :
143 8500 : needs_full_spectrum = .FALSE.
144 8500 : needs_homo = .FALSE.
145 :
146 11540 : SELECT CASE (precon_type)
147 : CASE (ot_precond_full_all)
148 3040 : needs_full_spectrum = .TRUE.
149 : ! both of them need the coefficients as fm's, more matrix mess
150 3040 : IF (use_mo_coeff_b) THEN
151 2764 : CALL copy_dbcsr_to_fm(mo_coeff_b, mo_coeff)
152 : END IF
153 : CASE (ot_precond_full_single)
154 38 : needs_homo = .TRUE.
155 : ! XXXX to be removed if homo estimate only is implemented
156 38 : needs_full_spectrum = .TRUE.
157 : CASE (ot_precond_full_kinetic, ot_precond_s_inverse, ot_precond_full_single_inverse)
158 : ! these should be happy without an estimate for the homo energy
159 : ! preconditioning can not depend on an absolute eigenvalue, only on eigenvalue differences
160 : CASE DEFAULT
161 8500 : CPABORT("The preconditioner is unknown ...")
162 : END SELECT
163 :
164 25386 : ALLOCATE (eigenvalues_ot(k))
165 8500 : energy_homo = 0.0_dp
166 8500 : IF (needs_full_spectrum) THEN
167 : ! XXXXXXXXXXXXXXXX do not touch the initial MOs, could be harmful for either
168 : ! the case of non-equivalent MOs but also for the derivate
169 : ! we could already have all eigenvalues e.g. full_all and we could skip this
170 : ! to be optimised later.
171 : ! one flaw is that not all SCF methods (i.e. that go over mo_derivs directly)
172 : ! have a 'valid' matrix_h... (we even don't know what evals are in that case)
173 3078 : IF (use_mo_coeff_b) THEN
174 : CALL calculate_subspace_eigenvalues(mo_coeff_b, matrix_h, &
175 : eigenvalues_ot, do_rotation=.FALSE., &
176 : para_env=mo_coeff%matrix_struct%para_env, &
177 2794 : blacs_env=mo_coeff%matrix_struct%context)
178 : ELSE
179 : CALL calculate_subspace_eigenvalues(mo_coeff, matrix_h, &
180 284 : eigenvalues_ot, do_rotation=.FALSE.)
181 : END IF
182 3078 : IF (k > 0) THEN
183 2972 : CPASSERT(nhomo > 0 .AND. nhomo <= k)
184 2972 : energy_homo = eigenvalues_ot(nhomo)
185 : END IF
186 : ELSE
187 5422 : IF (needs_homo) THEN
188 0 : CPABORT("Not yet implemented")
189 : END IF
190 : END IF
191 :
192 : ! After all bits and pieces of checking and initialization, here comes the
193 : ! part where the preconditioner matrix gets created and solved.
194 : ! This will give the matrices for later use
195 8500 : my_solver_type = solver_type
196 8500 : preconditioner_env%in_use = precon_type
197 8500 : preconditioner_env%cholesky_use = cholesky_reduce
198 8500 : IF (PRESENT(chol_type)) preconditioner_env%cholesky_use = chol_type
199 : preconditioner_env%in_use = precon_type
200 8500 : IF (nhomo == k) THEN
201 : CALL make_preconditioner_matrix(preconditioner_env, matrix_h, matrix_s, matrix_t, mo_coeff, &
202 8418 : energy_homo, eigenvalues_ot, energy_gap, my_solver_type)
203 : ELSE
204 : CALL cp_fm_struct_create(fm_struct, nrow_global=nao, ncol_global=nhomo, &
205 : context=preconditioner_env%ctxt, &
206 82 : para_env=preconditioner_env%para_env)
207 82 : CALL cp_fm_create(mo_occ, fm_struct)
208 82 : CALL cp_fm_to_fm(mo_coeff, mo_occ, nhomo)
209 82 : CALL cp_fm_struct_release(fm_struct)
210 : !
211 : CALL make_preconditioner_matrix(preconditioner_env, matrix_h, matrix_s, matrix_t, mo_occ, &
212 82 : energy_homo, eigenvalues_ot(1:nhomo), energy_gap, my_solver_type)
213 : !
214 82 : CALL cp_fm_release(mo_occ)
215 : END IF
216 :
217 8500 : CALL solve_preconditioner(my_solver_type, preconditioner_env, matrix_s, matrix_h)
218 :
219 : ! Here comes more matrix mess, make sure to output the correct matrix format,
220 : ! A bit pointless to convert the cholesky factorized version as it doesn't work in
221 : ! dbcsr form and will crash later,...
222 8500 : IF (my_convert_precond_to_dbcsr) THEN
223 7060 : CALL transfer_fm_to_dbcsr(preconditioner_env%fm, preconditioner_env%dbcsr_matrix, matrix_h)
224 : ELSE
225 : CALL transfer_dbcsr_to_fm(preconditioner_env%dbcsr_matrix, preconditioner_env%fm, &
226 1440 : preconditioner_env%para_env, preconditioner_env%ctxt)
227 : END IF
228 :
229 8500 : DEALLOCATE (eigenvalues_ot)
230 :
231 8500 : CALL timestop(handle)
232 :
233 8500 : END SUBROUTINE make_preconditioner
234 :
235 : ! **************************************************************************************************
236 : !> \brief Allows for a restart of the preconditioner
237 : !> depending on the method it purges all arrays or keeps them
238 : !> \param qs_env ...
239 : !> \param preconditioner ...
240 : !> \param prec_type ...
241 : !> \param nspins ...
242 : ! **************************************************************************************************
243 6713 : SUBROUTINE restart_preconditioner(qs_env, preconditioner, prec_type, nspins)
244 :
245 : TYPE(qs_environment_type), POINTER :: qs_env
246 : TYPE(preconditioner_p_type), DIMENSION(:), POINTER :: preconditioner
247 : INTEGER, INTENT(IN) :: prec_type, nspins
248 :
249 : INTEGER :: ispin
250 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
251 : TYPE(mp_para_env_type), POINTER :: para_env
252 :
253 6713 : NULLIFY (para_env, blacs_env)
254 6713 : CALL get_qs_env(qs_env, para_env=para_env, blacs_env=blacs_env)
255 :
256 6713 : IF (ASSOCIATED(preconditioner)) THEN
257 5910 : SELECT CASE (prec_type)
258 : CASE (ot_precond_full_all, ot_precond_full_single) ! these depend on the ks matrix
259 2742 : DO ispin = 1, SIZE(preconditioner)
260 1554 : CALL destroy_preconditioner(preconditioner(ispin)%preconditioner)
261 2742 : DEALLOCATE (preconditioner(ispin)%preconditioner)
262 : END DO
263 1188 : DEALLOCATE (preconditioner)
264 : CASE (ot_precond_none, ot_precond_full_kinetic, ot_precond_s_inverse, ot_precond_full_single_inverse) ! these are 'independent'
265 : ! do nothing
266 : CASE DEFAULT
267 4722 : CPABORT("")
268 : END SELECT
269 : END IF
270 :
271 : ! add an OT preconditioner if none is present
272 6713 : IF (.NOT. ASSOCIATED(preconditioner)) THEN
273 5743 : SELECT CASE (prec_type)
274 : CASE (ot_precond_full_all, ot_precond_full_single_inverse)
275 11081 : ALLOCATE (preconditioner(nspins))
276 : CASE DEFAULT
277 3794 : ALLOCATE (preconditioner(1))
278 : END SELECT
279 7183 : DO ispin = 1, SIZE(preconditioner)
280 4004 : ALLOCATE (preconditioner(ispin)%preconditioner)
281 : CALL init_preconditioner(preconditioner(ispin)%preconditioner, &
282 : para_env=para_env, &
283 7183 : blacs_env=blacs_env)
284 : END DO
285 : END IF
286 :
287 6713 : END SUBROUTINE restart_preconditioner
288 :
289 : ! **************************************************************************************************
290 : !> \brief ...
291 : !> \param qs_env ...
292 : !> \param mos ...
293 : !> \param matrix_ks ...
294 : !> \param matrix_s ...
295 : !> \param ot_preconditioner ...
296 : !> \param prec_type ...
297 : !> \param solver_type ...
298 : !> \param energy_gap ...
299 : !> \param nspins ...
300 : !> \param has_unit_metric ...
301 : !> \param convert_to_dbcsr ...
302 : !> \param chol_type ...
303 : !> \param full_mo_set ...
304 : ! **************************************************************************************************
305 6713 : SUBROUTINE prepare_preconditioner(qs_env, mos, matrix_ks, matrix_s, &
306 : ot_preconditioner, prec_type, solver_type, &
307 : energy_gap, nspins, has_unit_metric, &
308 : convert_to_dbcsr, chol_type, full_mo_set)
309 :
310 : TYPE(qs_environment_type), POINTER :: qs_env
311 : TYPE(mo_set_type), DIMENSION(:), INTENT(INOUT) :: mos
312 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s
313 : TYPE(preconditioner_p_type), DIMENSION(:), POINTER :: ot_preconditioner
314 : INTEGER, INTENT(IN) :: prec_type, solver_type
315 : REAL(dp), INTENT(IN) :: energy_gap
316 : INTEGER, INTENT(IN) :: nspins
317 : LOGICAL, INTENT(IN), OPTIONAL :: has_unit_metric, convert_to_dbcsr
318 : INTEGER, INTENT(IN), OPTIONAL :: chol_type
319 : LOGICAL, INTENT(IN), OPTIONAL :: full_mo_set
320 :
321 : CHARACTER(LEN=*), PARAMETER :: routineN = 'prepare_preconditioner'
322 :
323 : CHARACTER(LEN=default_string_length) :: msg
324 : INTEGER :: handle, icall, ispin, n_loops
325 : INTEGER, DIMENSION(5) :: nocc, norb
326 : LOGICAL :: do_co_rotate, my_convert_to_dbcsr, &
327 : my_full_mo_set, my_has_unit_metric, &
328 : use_mo_coeff_b
329 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
330 : TYPE(cp_fm_type), POINTER :: mo_coeff
331 6713 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: kinetic
332 : TYPE(dbcsr_type), POINTER :: matrix_t, mo_coeff_b
333 : TYPE(dft_control_type), POINTER :: dft_control
334 : TYPE(mp_para_env_type), POINTER :: para_env
335 :
336 6713 : CALL timeset(routineN, handle)
337 6713 : NULLIFY (matrix_t, mo_coeff_b, mo_coeff, kinetic, dft_control, para_env, blacs_env)
338 6713 : my_has_unit_metric = .FALSE.
339 6713 : IF (PRESENT(has_unit_metric)) my_has_unit_metric = has_unit_metric
340 6713 : my_convert_to_dbcsr = .TRUE.
341 6713 : IF (PRESENT(convert_to_dbcsr)) my_convert_to_dbcsr = convert_to_dbcsr
342 6713 : my_full_mo_set = .FALSE.
343 6713 : IF (PRESENT(full_mo_set)) my_full_mo_set = full_mo_set
344 :
345 : CALL get_qs_env(qs_env, &
346 : dft_control=dft_control, &
347 : para_env=para_env, &
348 6713 : blacs_env=blacs_env)
349 :
350 6713 : IF (dft_control%qs_control%semi_empirical .OR. dft_control%qs_control%dftb .OR. &
351 : dft_control%qs_control%xtb) THEN
352 2088 : IF (prec_type == ot_precond_full_kinetic) THEN
353 0 : msg = "Full_kinetic not available for semi-empirical methods"
354 0 : CPABORT(TRIM(msg))
355 : END IF
356 2088 : matrix_t => matrix_s(1)%matrix
357 : ELSE
358 4625 : CPASSERT(.NOT. my_has_unit_metric)
359 4625 : CALL get_qs_env(qs_env, kinetic=kinetic)
360 4625 : matrix_t => kinetic(1)%matrix
361 : END IF
362 :
363 : ! use full set of MOs or just occupied MOs
364 6713 : nocc = 0
365 6713 : norb = 0
366 6713 : IF (my_full_mo_set) THEN
367 34 : DO ispin = 1, nspins
368 18 : CALL get_mo_set(mo_set=mos(ispin), homo=nocc(ispin), nmo=norb(ispin))
369 34 : CALL set_mo_set(mo_set=mos(ispin), homo=norb(ispin))
370 : END DO
371 : END IF
372 : !determines how often make preconditioner is called, spin dependent methods have to be called twice
373 6713 : n_loops = 1
374 6713 : IF (prec_type == ot_precond_full_single_inverse) n_loops = nspins
375 : ! check whether we need the ev and rotate the MOs
376 1956 : SELECT CASE (prec_type)
377 : CASE (ot_precond_full_all)
378 : ! if one of these preconditioners is used every spin needs to call make_preconditioner
379 1956 : n_loops = nspins
380 :
381 1956 : do_co_rotate = ASSOCIATED(qs_env%mo_derivs)
382 9399 : DO ispin = 1, nspins
383 2686 : CALL get_mo_set(mo_set=mos(ispin), mo_coeff_b=mo_coeff_b, mo_coeff=mo_coeff)
384 2686 : use_mo_coeff_b = mos(ispin)%use_mo_coeff_b
385 4642 : IF (use_mo_coeff_b .AND. do_co_rotate) THEN
386 : CALL calculate_subspace_eigenvalues(mo_coeff_b, matrix_ks(ispin)%matrix, &
387 : do_rotation=.TRUE., &
388 : co_rotate=qs_env%mo_derivs(ispin)%matrix, &
389 : para_env=para_env, &
390 2668 : blacs_env=blacs_env)
391 18 : ELSEIF (use_mo_coeff_b) THEN
392 : CALL calculate_subspace_eigenvalues(mo_coeff_b, matrix_ks(ispin)%matrix, &
393 : do_rotation=.TRUE., &
394 : para_env=para_env, &
395 18 : blacs_env=blacs_env)
396 : ELSE
397 : CALL calculate_subspace_eigenvalues(mo_coeff, matrix_ks(ispin)%matrix, &
398 0 : do_rotation=.TRUE.)
399 : END IF
400 : END DO
401 : CASE DEFAULT
402 : ! No need to rotate the MOs
403 : END SELECT
404 :
405 : ! check whether we have a preconditioner
406 702 : SELECT CASE (prec_type)
407 : CASE (ot_precond_none)
408 1404 : DO ispin = 1, SIZE(ot_preconditioner)
409 1404 : ot_preconditioner(ispin)%preconditioner%in_use = 0
410 : END DO
411 : CASE DEFAULT
412 19710 : DO icall = 1, n_loops
413 12997 : IF (my_has_unit_metric) THEN
414 : CALL make_preconditioner(ot_preconditioner(icall)%preconditioner, &
415 : prec_type, &
416 : solver_type, &
417 : matrix_h=matrix_ks(icall)%matrix, &
418 : mo_set=mos(icall), &
419 : energy_gap=energy_gap, &
420 476 : convert_precond_to_dbcsr=my_convert_to_dbcsr)
421 : ELSE
422 : CALL make_preconditioner(ot_preconditioner(icall)%preconditioner, &
423 : prec_type, &
424 : solver_type, &
425 : matrix_h=matrix_ks(icall)%matrix, &
426 : matrix_s=matrix_s(1)%matrix, &
427 : matrix_t=matrix_t, &
428 : mo_set=mos(icall), &
429 : energy_gap=energy_gap, &
430 6510 : convert_precond_to_dbcsr=my_convert_to_dbcsr, chol_type=chol_type)
431 : END IF
432 : END DO
433 : END SELECT
434 :
435 : ! reset homo values
436 6713 : IF (my_full_mo_set) THEN
437 34 : DO ispin = 1, nspins
438 34 : CALL set_mo_set(mo_set=mos(ispin), homo=nocc(ispin))
439 : END DO
440 : END IF
441 :
442 6713 : CALL timestop(handle)
443 :
444 6713 : END SUBROUTINE prepare_preconditioner
445 :
446 : END MODULE preconditioner
447 :
|
