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_struct, ONLY: cp_fm_struct_create,&
27 : cp_fm_struct_release,&
28 : cp_fm_struct_type
29 : USE cp_fm_types, ONLY: cp_fm_create,&
30 : cp_fm_get_info,&
31 : cp_fm_release,&
32 : cp_fm_to_fm,&
33 : cp_fm_type
34 : USE input_constants, ONLY: &
35 : cholesky_reduce, ot_precond_full_all, ot_precond_full_kinetic, ot_precond_full_single, &
36 : ot_precond_full_single_inverse, ot_precond_none, ot_precond_s_inverse, &
37 : ot_precond_solver_default, ot_precond_solver_inv_chol, ot_precond_solver_update
38 : USE kinds, ONLY: default_string_length,&
39 : dp
40 : USE message_passing, ONLY: mp_para_env_type
41 : USE preconditioner_apply, ONLY: apply_preconditioner_cfm_complex,&
42 : apply_preconditioner_dbcsr,&
43 : apply_preconditioner_dbcsr_complex,&
44 : apply_preconditioner_fm
45 : USE preconditioner_makes, ONLY: make_complex_full_all,&
46 : make_complex_full_kinetic,&
47 : make_complex_full_s_inverse,&
48 : make_complex_full_single,&
49 : make_complex_full_single_inverse,&
50 : make_preconditioner_matrix
51 : USE preconditioner_solvers, ONLY: solve_preconditioner,&
52 : transfer_dbcsr_to_fm,&
53 : transfer_fm_to_dbcsr
54 : USE preconditioner_types, ONLY: destroy_preconditioner,&
55 : init_preconditioner,&
56 : preconditioner_p_type,&
57 : preconditioner_type
58 : USE qs_environment_types, ONLY: get_qs_env,&
59 : qs_environment_type
60 : USE qs_mo_methods, ONLY: calculate_subspace_eigenvalues
61 : USE qs_mo_types, ONLY: get_mo_set,&
62 : mo_set_type,&
63 : set_mo_set
64 : #include "./base/base_uses.f90"
65 :
66 : IMPLICIT NONE
67 :
68 : PRIVATE
69 :
70 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'preconditioner'
71 :
72 : PUBLIC :: make_preconditioner, make_preconditioner_complex_full_all, &
73 : make_preconditioner_complex_full_kinetic, &
74 : make_preconditioner_complex_full_s_inverse, &
75 : make_preconditioner_complex_full_single, &
76 : make_preconditioner_complex_full_single_inverse, restart_preconditioner, &
77 : dbcsr_pair_to_cfm
78 : PUBLIC :: apply_preconditioner, prepare_preconditioner
79 :
80 : ! The public interface for apply preconditioner, the routines can be found in preconditioner_apply.F
81 : INTERFACE apply_preconditioner
82 : MODULE PROCEDURE apply_preconditioner_cfm_complex
83 : MODULE PROCEDURE apply_preconditioner_dbcsr
84 : MODULE PROCEDURE apply_preconditioner_dbcsr_complex
85 : MODULE PROCEDURE apply_preconditioner_fm
86 : END INTERFACE
87 :
88 : ! **************************************************************************************************
89 :
90 : CONTAINS
91 :
92 : ! **************************************************************************************************
93 :
94 : ! creates a preconditioner for the system (H-energy_homo S)
95 : ! this preconditioner is (must be) symmetric positive definite.
96 : ! currently uses a atom-block-diagonal form
97 : ! each block will be ....
98 : ! might overwrite matrix_h, matrix_t
99 :
100 : ! **************************************************************************************************
101 : !> \brief ...
102 : !> \param preconditioner_env ...
103 : !> \param precon_type ...
104 : !> \param solver_type ...
105 : !> \param matrix_h ...
106 : !> \param matrix_s ...
107 : !> \param matrix_t ...
108 : !> \param mo_set ...
109 : !> \param energy_gap ...
110 : !> \param convert_precond_to_dbcsr ...
111 : !> \param chol_type ...
112 : !> \par History
113 : !> 09.2014 removed some unused or unfinished methods
114 : !> removed sparse preconditioners and the
115 : !> sparse approximate inverse at rev 14341 [Florian Schiffmann]
116 : ! **************************************************************************************************
117 9592 : SUBROUTINE make_preconditioner(preconditioner_env, precon_type, solver_type, matrix_h, matrix_s, &
118 : matrix_t, mo_set, energy_gap, convert_precond_to_dbcsr, chol_type)
119 :
120 : TYPE(preconditioner_type) :: preconditioner_env
121 : INTEGER, INTENT(IN) :: precon_type, solver_type
122 : TYPE(dbcsr_type), POINTER :: matrix_h
123 : TYPE(dbcsr_type), OPTIONAL, POINTER :: matrix_s, matrix_t
124 : TYPE(mo_set_type), INTENT(IN) :: mo_set
125 : REAL(KIND=dp) :: energy_gap
126 : LOGICAL, INTENT(IN), OPTIONAL :: convert_precond_to_dbcsr
127 : INTEGER, INTENT(IN), OPTIONAL :: chol_type
128 :
129 : CHARACTER(len=*), PARAMETER :: routineN = 'make_preconditioner'
130 :
131 : INTEGER :: handle, k, my_solver_type, nao, nhomo
132 : LOGICAL :: my_convert_precond_to_dbcsr, &
133 : needs_full_spectrum, needs_homo, &
134 : use_mo_coeff_b
135 : REAL(KIND=dp) :: energy_homo
136 9592 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eigenvalues_ot
137 : TYPE(cp_fm_struct_type), POINTER :: fm_struct
138 : TYPE(cp_fm_type) :: mo_occ
139 : TYPE(cp_fm_type), POINTER :: mo_coeff
140 : TYPE(dbcsr_type), POINTER :: mo_coeff_b
141 :
142 9592 : CALL timeset(routineN, handle)
143 :
144 9592 : CALL get_mo_set(mo_set=mo_set, mo_coeff=mo_coeff, mo_coeff_b=mo_coeff_b, homo=nhomo)
145 9592 : use_mo_coeff_b = mo_set%use_mo_coeff_b
146 9592 : CALL cp_fm_get_info(mo_coeff, ncol_global=k, nrow_global=nao)
147 :
148 : ! Starting some matrix mess, check where to store the result in preconditioner_env, fm or dbcsr_matrix
149 9592 : my_convert_precond_to_dbcsr = .FALSE.
150 9592 : IF (PRESENT(convert_precond_to_dbcsr)) my_convert_precond_to_dbcsr = convert_precond_to_dbcsr
151 :
152 : ! Thanks to the mess with the matrices we need to make sure in this case that the
153 : ! Previous inverse is properly stored as a sparse matrix, fm gets deallocated here
154 : ! if it wasn't anyway
155 9592 : IF (preconditioner_env%solver == ot_precond_solver_update) THEN
156 4 : CALL transfer_fm_to_dbcsr(preconditioner_env%fm, preconditioner_env%dbcsr_matrix, matrix_h)
157 : END IF
158 :
159 9592 : needs_full_spectrum = .FALSE.
160 9592 : needs_homo = .FALSE.
161 :
162 13306 : SELECT CASE (precon_type)
163 : CASE (ot_precond_full_all)
164 3714 : needs_full_spectrum = .TRUE.
165 : ! both of them need the coefficients as fm's, more matrix mess
166 3714 : IF (use_mo_coeff_b) THEN
167 3410 : CALL copy_dbcsr_to_fm(mo_coeff_b, mo_coeff)
168 : END IF
169 : CASE (ot_precond_full_single)
170 38 : needs_homo = .TRUE.
171 : ! XXXX to be removed if homo estimate only is implemented
172 38 : needs_full_spectrum = .TRUE.
173 : CASE (ot_precond_full_kinetic, ot_precond_s_inverse, ot_precond_full_single_inverse)
174 : ! these should be happy without an estimate for the homo energy
175 : ! preconditioning can not depend on an absolute eigenvalue, only on eigenvalue differences
176 : CASE DEFAULT
177 9592 : CPABORT("The preconditioner is unknown ...")
178 : END SELECT
179 :
180 28592 : ALLOCATE (eigenvalues_ot(k))
181 9592 : energy_homo = 0.0_dp
182 9592 : IF (needs_full_spectrum) THEN
183 : ! XXXXXXXXXXXXXXXX do not touch the initial MOs, could be harmful for either
184 : ! the case of non-equivalent MOs but also for the derivate
185 : ! we could already have all eigenvalues e.g. full_all and we could skip this
186 : ! to be optimised later.
187 : ! one flaw is that not all SCF methods (i.e. that go over mo_derivs directly)
188 : ! have a 'valid' matrix_h... (we even don't know what evals are in that case)
189 3752 : IF (use_mo_coeff_b) THEN
190 : CALL calculate_subspace_eigenvalues(mo_coeff_b, matrix_h, &
191 : eigenvalues_ot, do_rotation=.FALSE., &
192 : para_env=mo_coeff%matrix_struct%para_env, &
193 3440 : blacs_env=mo_coeff%matrix_struct%context)
194 : ELSE
195 : CALL calculate_subspace_eigenvalues(mo_coeff, matrix_h, &
196 312 : eigenvalues_ot, do_rotation=.FALSE.)
197 : END IF
198 3752 : IF (k > 0) THEN
199 3642 : CPASSERT(nhomo > 0 .AND. nhomo <= k)
200 3642 : energy_homo = eigenvalues_ot(nhomo)
201 : END IF
202 : ELSE
203 5840 : IF (needs_homo) THEN
204 0 : CPABORT("Not yet implemented")
205 : END IF
206 : END IF
207 :
208 : ! After all bits and pieces of checking and initialization, here comes the
209 : ! part where the preconditioner matrix gets created and solved.
210 : ! This will give the matrices for later use
211 9592 : my_solver_type = solver_type
212 9592 : preconditioner_env%in_use = precon_type
213 9592 : preconditioner_env%cholesky_use = cholesky_reduce
214 9592 : IF (PRESENT(chol_type)) preconditioner_env%cholesky_use = chol_type
215 : preconditioner_env%in_use = precon_type
216 9592 : IF (nhomo == k) THEN
217 : CALL make_preconditioner_matrix(preconditioner_env, matrix_h, matrix_s, matrix_t, mo_coeff, &
218 9504 : energy_homo, eigenvalues_ot, energy_gap, my_solver_type)
219 : ELSE
220 : CALL cp_fm_struct_create(fm_struct, nrow_global=nao, ncol_global=nhomo, &
221 : context=preconditioner_env%ctxt, &
222 88 : para_env=preconditioner_env%para_env)
223 88 : CALL cp_fm_create(mo_occ, fm_struct)
224 88 : CALL cp_fm_to_fm(mo_coeff, mo_occ, nhomo)
225 88 : CALL cp_fm_struct_release(fm_struct)
226 : !
227 : CALL make_preconditioner_matrix(preconditioner_env, matrix_h, matrix_s, matrix_t, mo_occ, &
228 88 : energy_homo, eigenvalues_ot(1:nhomo), energy_gap, my_solver_type)
229 : !
230 88 : CALL cp_fm_release(mo_occ)
231 : END IF
232 :
233 9592 : CALL solve_preconditioner(my_solver_type, preconditioner_env, matrix_s, matrix_h)
234 :
235 : ! Here comes more matrix mess, make sure to output the correct matrix format,
236 : ! A bit pointless to convert the cholesky factorized version as it doesn't work in
237 : ! dbcsr form and will crash later,...
238 9592 : IF (my_convert_precond_to_dbcsr) THEN
239 7916 : CALL transfer_fm_to_dbcsr(preconditioner_env%fm, preconditioner_env%dbcsr_matrix, matrix_h)
240 : ELSE
241 : CALL transfer_dbcsr_to_fm(preconditioner_env%dbcsr_matrix, preconditioner_env%fm, &
242 1676 : preconditioner_env%para_env, preconditioner_env%ctxt)
243 : END IF
244 :
245 9592 : DEALLOCATE (eigenvalues_ot)
246 :
247 9592 : CALL timestop(handle)
248 :
249 9592 : END SUBROUTINE make_preconditioner
250 :
251 : ! **************************************************************************************************
252 : !> \brief Construct FULL_ALL directly from one complex H(k), S(k), and C(k) channel.
253 : !> \param preconditioner_env preconditioner storage
254 : !> \param matrix_c_re real part of the active OT reference orbitals
255 : !> \param matrix_c_im imaginary part of the active OT reference orbitals
256 : !> \param matrix_h_re real part of H(k)
257 : !> \param matrix_h_im imaginary part of H(k)
258 : !> \param matrix_s_re real part of S(k)
259 : !> \param matrix_s_im imaginary part of S(k)
260 : !> \param mo_set reference-orbital energy labels
261 : !> \param energy_gap denominator floor
262 : !> \param solver_type inverse application selected in the OT input
263 : ! **************************************************************************************************
264 520 : SUBROUTINE make_preconditioner_complex_full_all(preconditioner_env, &
265 : matrix_c_re, matrix_c_im, &
266 : matrix_h_re, matrix_h_im, &
267 : matrix_s_re, matrix_s_im, &
268 : mo_set, energy_gap, solver_type)
269 :
270 : TYPE(preconditioner_type) :: preconditioner_env
271 : TYPE(dbcsr_type), POINTER :: matrix_c_re, matrix_c_im, matrix_h_re, &
272 : matrix_h_im, matrix_s_re, matrix_s_im
273 : TYPE(mo_set_type), INTENT(IN) :: mo_set
274 : REAL(KIND=dp), INTENT(IN) :: energy_gap
275 : INTEGER, INTENT(IN) :: solver_type
276 :
277 : CHARACTER(len=*), PARAMETER :: routineN = 'make_preconditioner_complex_full_all'
278 :
279 : INTEGER :: handle, k
280 104 : REAL(KIND=dp), DIMENSION(:), POINTER :: eigenvalues
281 : TYPE(cp_cfm_type) :: matrix_c, matrix_h, matrix_s
282 :
283 104 : CALL timeset(routineN, handle)
284 :
285 104 : NULLIFY (eigenvalues)
286 104 : CPASSERT(ASSOCIATED(matrix_c_re))
287 104 : CPASSERT(ASSOCIATED(matrix_c_im))
288 104 : CPASSERT(ASSOCIATED(matrix_h_re))
289 104 : CPASSERT(ASSOCIATED(matrix_h_im))
290 104 : CPASSERT(ASSOCIATED(matrix_s_re))
291 104 : CPASSERT(ASSOCIATED(matrix_s_im))
292 104 : IF (solver_type /= ot_precond_solver_default) THEN
293 0 : CPABORT('Complex FULL_ALL supports only PRECOND_SOLVER DEFAULT')
294 : END IF
295 104 : CALL dbcsr_get_info(matrix_c_re, nfullcols_total=k)
296 104 : CPASSERT(k > 0)
297 104 : CALL get_mo_set(mo_set, eigenvalues=eigenvalues)
298 104 : CPASSERT(ASSOCIATED(eigenvalues))
299 104 : CPASSERT(SIZE(eigenvalues) >= k)
300 :
301 : CALL dbcsr_pair_to_cfm(matrix_c_re, matrix_c_im, preconditioner_env, &
302 104 : 'complex FULL_ALL C', matrix_c)
303 : CALL dbcsr_pair_to_cfm(matrix_h_re, matrix_h_im, preconditioner_env, &
304 104 : 'complex FULL_ALL H', matrix_h)
305 : CALL dbcsr_pair_to_cfm(matrix_s_re, matrix_s_im, preconditioner_env, &
306 104 : 'complex FULL_ALL S', matrix_s)
307 : CALL make_complex_full_all(preconditioner_env, matrix_c, matrix_h, matrix_s, &
308 104 : eigenvalues(1:k), energy_gap)
309 :
310 104 : CALL cp_cfm_release(matrix_c)
311 104 : CALL cp_cfm_release(matrix_s)
312 104 : CALL cp_cfm_release(matrix_h)
313 :
314 104 : CALL timestop(handle)
315 :
316 104 : END SUBROUTINE make_preconditioner_complex_full_all
317 :
318 : ! **************************************************************************************************
319 : !> \brief Construct a complex FULL_SINGLE preconditioner from H(k) and S(k).
320 : !> \param preconditioner_env preconditioner storage
321 : !> \param matrix_h_re real part of H(k)
322 : !> \param matrix_h_im imaginary part of H(k)
323 : !> \param matrix_s_re real part of S(k)
324 : !> \param matrix_s_im imaginary part of S(k)
325 : !> \param mo_set orbital energy labels defining the occupied edge
326 : !> \param energy_gap denominator floor
327 : !> \param solver_type requested inverse solver
328 : ! **************************************************************************************************
329 104 : SUBROUTINE make_preconditioner_complex_full_single(preconditioner_env, &
330 : matrix_h_re, matrix_h_im, &
331 : matrix_s_re, matrix_s_im, &
332 : mo_set, energy_gap, solver_type)
333 :
334 : TYPE(preconditioner_type) :: preconditioner_env
335 : TYPE(dbcsr_type), POINTER :: matrix_h_re, matrix_h_im, matrix_s_re, &
336 : matrix_s_im
337 : TYPE(mo_set_type), INTENT(IN) :: mo_set
338 : REAL(KIND=dp), INTENT(IN) :: energy_gap
339 : INTEGER, INTENT(IN) :: solver_type
340 :
341 : INTEGER :: homo
342 26 : REAL(KIND=dp), DIMENSION(:), POINTER :: eigenvalues
343 : TYPE(cp_cfm_type) :: matrix_h, matrix_s
344 :
345 26 : NULLIFY (eigenvalues)
346 0 : CPASSERT(ASSOCIATED(matrix_h_re))
347 26 : CPASSERT(ASSOCIATED(matrix_h_im))
348 26 : CPASSERT(ASSOCIATED(matrix_s_re))
349 26 : CPASSERT(ASSOCIATED(matrix_s_im))
350 26 : IF (solver_type /= ot_precond_solver_default) THEN
351 0 : CPABORT('Complex FULL_SINGLE supports only PRECOND_SOLVER DEFAULT')
352 : END IF
353 26 : CALL get_mo_set(mo_set, homo=homo, eigenvalues=eigenvalues)
354 26 : CPASSERT(ASSOCIATED(eigenvalues))
355 26 : CPASSERT(homo > 0 .AND. homo <= SIZE(eigenvalues))
356 :
357 : CALL dbcsr_pair_to_cfm(matrix_h_re, matrix_h_im, preconditioner_env, &
358 26 : 'complex FULL_SINGLE H', matrix_h)
359 : CALL dbcsr_pair_to_cfm(matrix_s_re, matrix_s_im, preconditioner_env, &
360 26 : 'complex FULL_SINGLE S', matrix_s)
361 : CALL make_complex_full_single(preconditioner_env, matrix_h, matrix_s, &
362 26 : eigenvalues(homo), energy_gap)
363 26 : CALL cp_cfm_release(matrix_s)
364 26 : CALL cp_cfm_release(matrix_h)
365 :
366 26 : END SUBROUTINE make_preconditioner_complex_full_single
367 :
368 : ! **************************************************************************************************
369 : !> \brief Construct a complex FULL_SINGLE_INVERSE preconditioner without discarding Im(H,S,C).
370 : !> \param preconditioner_env preconditioner storage
371 : !> \param matrix_c_re real part of the occupied reference orbitals
372 : !> \param matrix_c_im imaginary part of the occupied reference orbitals
373 : !> \param matrix_h_re real part of H(k)
374 : !> \param matrix_h_im imaginary part of H(k)
375 : !> \param matrix_s_re real part of S(k)
376 : !> \param matrix_s_im imaginary part of S(k)
377 : !> \param energy_gap lower spectral bound
378 : !> \param solver_type requested inverse solver
379 : ! **************************************************************************************************
380 484 : SUBROUTINE make_preconditioner_complex_full_single_inverse(preconditioner_env, &
381 : matrix_c_re, matrix_c_im, &
382 : matrix_h_re, matrix_h_im, &
383 : matrix_s_re, matrix_s_im, &
384 : energy_gap, solver_type)
385 :
386 : TYPE(preconditioner_type) :: preconditioner_env
387 : TYPE(dbcsr_type), POINTER :: matrix_c_re, matrix_c_im, matrix_h_re, &
388 : matrix_h_im, matrix_s_re, matrix_s_im
389 : REAL(KIND=dp), INTENT(IN) :: energy_gap
390 : INTEGER, INTENT(IN) :: solver_type
391 :
392 : INTEGER :: k, n
393 : TYPE(cp_cfm_type) :: matrix_c, matrix_h, matrix_s
394 : TYPE(cp_fm_struct_type), POINTER :: fm_struct_mo, fm_struct_square
395 : TYPE(cp_fm_type) :: matrix_c_im_fm, matrix_c_re_fm, &
396 : matrix_h_im_fm, matrix_h_re_fm, &
397 : matrix_s_im_fm, matrix_s_re_fm
398 :
399 44 : NULLIFY (fm_struct_mo, fm_struct_square)
400 0 : CPASSERT(ASSOCIATED(matrix_c_re))
401 44 : CPASSERT(ASSOCIATED(matrix_c_im))
402 44 : CPASSERT(ASSOCIATED(matrix_h_re))
403 44 : CPASSERT(ASSOCIATED(matrix_h_im))
404 44 : CPASSERT(ASSOCIATED(matrix_s_re))
405 44 : CPASSERT(ASSOCIATED(matrix_s_im))
406 44 : IF (solver_type /= ot_precond_solver_default .AND. &
407 : solver_type /= ot_precond_solver_inv_chol) THEN
408 0 : CPABORT("Complex FULL_SINGLE_INVERSE supports only PRECOND_SOLVER DEFAULT/INVERSE_CHOLESKY")
409 : END IF
410 44 : CALL dbcsr_get_info(matrix_c_re, nfullrows_total=n, nfullcols_total=k)
411 44 : CPASSERT(n > 0 .AND. k > 0 .AND. k <= n)
412 :
413 : CALL cp_fm_struct_create(fm_struct_mo, nrow_global=n, ncol_global=k, &
414 : context=preconditioner_env%ctxt, &
415 44 : para_env=preconditioner_env%para_env)
416 : CALL cp_fm_struct_create(fm_struct_square, nrow_global=n, ncol_global=n, &
417 : context=preconditioner_env%ctxt, &
418 44 : para_env=preconditioner_env%para_env)
419 44 : CALL cp_fm_create(matrix_c_re_fm, fm_struct_mo, name='complex preconditioner C real')
420 44 : CALL cp_fm_create(matrix_c_im_fm, fm_struct_mo, name='complex preconditioner C imaginary')
421 44 : CALL cp_fm_create(matrix_h_re_fm, fm_struct_square, name='complex preconditioner H real')
422 44 : CALL cp_fm_create(matrix_h_im_fm, fm_struct_square, name='complex preconditioner H imaginary')
423 44 : CALL cp_fm_create(matrix_s_re_fm, fm_struct_square, name='complex preconditioner S real')
424 44 : CALL cp_fm_create(matrix_s_im_fm, fm_struct_square, name='complex preconditioner S imaginary')
425 44 : CALL cp_fm_struct_release(fm_struct_mo)
426 44 : CALL cp_fm_struct_release(fm_struct_square)
427 44 : CALL copy_dbcsr_to_fm(matrix_c_re, matrix_c_re_fm)
428 44 : CALL copy_dbcsr_to_fm(matrix_c_im, matrix_c_im_fm)
429 44 : CALL copy_dbcsr_to_fm(matrix_h_re, matrix_h_re_fm)
430 44 : CALL copy_dbcsr_to_fm(matrix_h_im, matrix_h_im_fm)
431 44 : CALL copy_dbcsr_to_fm(matrix_s_re, matrix_s_re_fm)
432 44 : CALL copy_dbcsr_to_fm(matrix_s_im, matrix_s_im_fm)
433 :
434 44 : CALL cp_cfm_create(matrix_c, matrix_c_re_fm%matrix_struct, name='complex preconditioner C')
435 44 : CALL cp_cfm_create(matrix_h, matrix_h_re_fm%matrix_struct, name='complex preconditioner H')
436 44 : CALL cp_cfm_create(matrix_s, matrix_s_re_fm%matrix_struct, name='complex preconditioner S')
437 44 : CALL cp_fm_to_cfm(matrix_c_re_fm, matrix_c_im_fm, matrix_c)
438 44 : CALL cp_fm_to_cfm(matrix_h_re_fm, matrix_h_im_fm, matrix_h)
439 44 : CALL cp_fm_to_cfm(matrix_s_re_fm, matrix_s_im_fm, matrix_s)
440 : CALL make_complex_full_single_inverse(preconditioner_env, matrix_c, matrix_h, matrix_s, &
441 44 : energy_gap)
442 44 : preconditioner_env%solver = solver_type
443 :
444 44 : CALL cp_cfm_release(matrix_s)
445 44 : CALL cp_cfm_release(matrix_h)
446 44 : CALL cp_cfm_release(matrix_c)
447 44 : CALL cp_fm_release(matrix_s_im_fm)
448 44 : CALL cp_fm_release(matrix_s_re_fm)
449 44 : CALL cp_fm_release(matrix_h_im_fm)
450 44 : CALL cp_fm_release(matrix_h_re_fm)
451 44 : CALL cp_fm_release(matrix_c_im_fm)
452 44 : CALL cp_fm_release(matrix_c_re_fm)
453 :
454 44 : END SUBROUTINE make_preconditioner_complex_full_single_inverse
455 :
456 : ! **************************************************************************************************
457 : !> \brief Construct a complex FULL_S_INVERSE preconditioner.
458 : !> \param preconditioner_env preconditioner storage
459 : !> \param matrix_s_re real part of S(k)
460 : !> \param matrix_s_im imaginary part of S(k)
461 : !> \param solver_type requested inverse solver
462 : ! **************************************************************************************************
463 180 : SUBROUTINE make_preconditioner_complex_full_s_inverse(preconditioner_env, &
464 : matrix_s_re, matrix_s_im, solver_type)
465 :
466 : TYPE(preconditioner_type) :: preconditioner_env
467 : TYPE(dbcsr_type), POINTER :: matrix_s_re, matrix_s_im
468 : INTEGER, INTENT(IN) :: solver_type
469 :
470 : TYPE(cp_cfm_type) :: matrix_s
471 :
472 90 : IF (solver_type /= ot_precond_solver_default .AND. &
473 : solver_type /= ot_precond_solver_inv_chol) THEN
474 0 : CPABORT("Complex FULL_S_INVERSE supports only PRECOND_SOLVER DEFAULT/INVERSE_CHOLESKY")
475 : END IF
476 : CALL dbcsr_pair_to_cfm(matrix_s_re, matrix_s_im, preconditioner_env, &
477 90 : 'complex preconditioner S', matrix_s)
478 90 : CALL make_complex_full_s_inverse(preconditioner_env, matrix_s)
479 90 : preconditioner_env%solver = solver_type
480 90 : CALL cp_cfm_release(matrix_s)
481 :
482 90 : END SUBROUTINE make_preconditioner_complex_full_s_inverse
483 :
484 : ! **************************************************************************************************
485 : !> \brief Construct a complex FULL_KINETIC preconditioner.
486 : !> \param preconditioner_env preconditioner storage
487 : !> \param matrix_t_re real part of T(k)
488 : !> \param matrix_t_im imaginary part of T(k)
489 : !> \param matrix_s_re real part of S(k)
490 : !> \param matrix_s_im imaginary part of S(k)
491 : !> \param energy_gap non-negative overlap shift
492 : !> \param solver_type requested inverse solver
493 : ! **************************************************************************************************
494 360 : SUBROUTINE make_preconditioner_complex_full_kinetic(preconditioner_env, &
495 : matrix_t_re, matrix_t_im, &
496 : matrix_s_re, matrix_s_im, &
497 : energy_gap, solver_type)
498 :
499 : TYPE(preconditioner_type) :: preconditioner_env
500 : TYPE(dbcsr_type), POINTER :: matrix_t_re, matrix_t_im, matrix_s_re, &
501 : matrix_s_im
502 : REAL(KIND=dp), INTENT(IN) :: energy_gap
503 : INTEGER, INTENT(IN) :: solver_type
504 :
505 : TYPE(cp_cfm_type) :: matrix_s, matrix_t
506 :
507 120 : IF (solver_type /= ot_precond_solver_default .AND. &
508 : solver_type /= ot_precond_solver_inv_chol) THEN
509 0 : CPABORT("Complex FULL_KINETIC supports only PRECOND_SOLVER DEFAULT/INVERSE_CHOLESKY")
510 : END IF
511 : CALL dbcsr_pair_to_cfm(matrix_t_re, matrix_t_im, preconditioner_env, &
512 120 : 'complex preconditioner T', matrix_t)
513 : CALL dbcsr_pair_to_cfm(matrix_s_re, matrix_s_im, preconditioner_env, &
514 120 : 'complex preconditioner S', matrix_s)
515 120 : CALL make_complex_full_kinetic(preconditioner_env, matrix_t, matrix_s, energy_gap)
516 120 : preconditioner_env%solver = solver_type
517 120 : CALL cp_cfm_release(matrix_s)
518 120 : CALL cp_cfm_release(matrix_t)
519 :
520 120 : END SUBROUTINE make_preconditioner_complex_full_kinetic
521 :
522 : ! **************************************************************************************************
523 : !> \brief Copy a real/imaginary DBCSR pair to one distributed complex full matrix.
524 : !> \param matrix_re real matrix component
525 : !> \param matrix_im imaginary matrix component
526 : !> \param preconditioner_env source of the distribution context
527 : !> \param matrix_name matrix label
528 : !> \param matrix complex output matrix
529 : ! **************************************************************************************************
530 2442 : SUBROUTINE dbcsr_pair_to_cfm(matrix_re, matrix_im, preconditioner_env, matrix_name, matrix)
531 :
532 : TYPE(dbcsr_type), POINTER :: matrix_re, matrix_im
533 : TYPE(preconditioner_type) :: preconditioner_env
534 : CHARACTER(LEN=*), INTENT(IN) :: matrix_name
535 : TYPE(cp_cfm_type), INTENT(OUT) :: matrix
536 :
537 : INTEGER :: ncol, nrow
538 : TYPE(cp_fm_struct_type), POINTER :: fm_struct
539 : TYPE(cp_fm_type) :: matrix_im_fm, matrix_re_fm
540 :
541 814 : NULLIFY (fm_struct)
542 0 : CPASSERT(ASSOCIATED(matrix_re))
543 814 : CPASSERT(ASSOCIATED(matrix_im))
544 814 : CALL dbcsr_get_info(matrix_re, nfullrows_total=nrow, nfullcols_total=ncol)
545 : CALL cp_fm_struct_create(fm_struct, nrow_global=nrow, ncol_global=ncol, &
546 : context=preconditioner_env%ctxt, &
547 814 : para_env=preconditioner_env%para_env)
548 814 : CALL cp_fm_create(matrix_re_fm, fm_struct, name=TRIM(matrix_name)//' real')
549 814 : CALL cp_fm_create(matrix_im_fm, fm_struct, name=TRIM(matrix_name)//' imaginary')
550 814 : CALL cp_fm_struct_release(fm_struct)
551 814 : CALL copy_dbcsr_to_fm(matrix_re, matrix_re_fm)
552 814 : CALL copy_dbcsr_to_fm(matrix_im, matrix_im_fm)
553 814 : CALL cp_cfm_create(matrix, matrix_re_fm%matrix_struct, name=matrix_name)
554 814 : CALL cp_fm_to_cfm(matrix_re_fm, matrix_im_fm, matrix)
555 814 : CALL cp_fm_release(matrix_im_fm)
556 814 : CALL cp_fm_release(matrix_re_fm)
557 :
558 814 : END SUBROUTINE dbcsr_pair_to_cfm
559 :
560 : ! **************************************************************************************************
561 : !> \brief Allows for a restart of the preconditioner
562 : !> depending on the method it purges all arrays or keeps them
563 : !> \param qs_env ...
564 : !> \param preconditioner ...
565 : !> \param prec_type ...
566 : !> \param nspins ...
567 : ! **************************************************************************************************
568 7443 : SUBROUTINE restart_preconditioner(qs_env, preconditioner, prec_type, nspins)
569 :
570 : TYPE(qs_environment_type), POINTER :: qs_env
571 : TYPE(preconditioner_p_type), DIMENSION(:), POINTER :: preconditioner
572 : INTEGER, INTENT(IN) :: prec_type, nspins
573 :
574 : INTEGER :: ispin
575 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
576 : TYPE(mp_para_env_type), POINTER :: para_env
577 :
578 7443 : NULLIFY (para_env, blacs_env)
579 7443 : CALL get_qs_env(qs_env, para_env=para_env, blacs_env=blacs_env)
580 :
581 7443 : IF (ASSOCIATED(preconditioner)) THEN
582 6638 : SELECT CASE (prec_type)
583 : CASE (ot_precond_full_all, ot_precond_full_single) ! these depend on the ks matrix
584 3424 : DO ispin = 1, SIZE(preconditioner)
585 1908 : CALL destroy_preconditioner(preconditioner(ispin)%preconditioner)
586 3424 : DEALLOCATE (preconditioner(ispin)%preconditioner)
587 : END DO
588 1516 : DEALLOCATE (preconditioner)
589 : CASE (ot_precond_none, ot_precond_full_kinetic, ot_precond_s_inverse, &
590 : ot_precond_full_single_inverse) ! these are 'independent'
591 : ! do nothing
592 : CASE DEFAULT
593 5122 : CPABORT("Unknown preconditioner type")
594 : END SELECT
595 : END IF
596 :
597 : ! add an OT preconditioner if none is present
598 7443 : IF (.NOT. ASSOCIATED(preconditioner)) THEN
599 7033 : SELECT CASE (prec_type)
600 : CASE (ot_precond_full_all, ot_precond_full_single_inverse)
601 13683 : ALLOCATE (preconditioner(nspins))
602 : CASE DEFAULT
603 4478 : ALLOCATE (preconditioner(1))
604 : END SELECT
605 8573 : DO ispin = 1, SIZE(preconditioner)
606 4736 : ALLOCATE (preconditioner(ispin)%preconditioner)
607 : CALL init_preconditioner(preconditioner(ispin)%preconditioner, &
608 : para_env=para_env, &
609 8573 : blacs_env=blacs_env)
610 : END DO
611 : END IF
612 :
613 7443 : END SUBROUTINE restart_preconditioner
614 :
615 : ! **************************************************************************************************
616 : !> \brief ...
617 : !> \param qs_env ...
618 : !> \param mos ...
619 : !> \param matrix_ks ...
620 : !> \param matrix_s ...
621 : !> \param ot_preconditioner ...
622 : !> \param prec_type ...
623 : !> \param solver_type ...
624 : !> \param energy_gap ...
625 : !> \param nspins ...
626 : !> \param has_unit_metric ...
627 : !> \param convert_to_dbcsr ...
628 : !> \param chol_type ...
629 : !> \param full_mo_set ...
630 : !> \param chebyshev_degree ...
631 : ! **************************************************************************************************
632 7443 : SUBROUTINE prepare_preconditioner(qs_env, mos, matrix_ks, matrix_s, &
633 : ot_preconditioner, prec_type, solver_type, &
634 : energy_gap, nspins, has_unit_metric, &
635 : convert_to_dbcsr, chol_type, full_mo_set, chebyshev_degree)
636 :
637 : TYPE(qs_environment_type), POINTER :: qs_env
638 : TYPE(mo_set_type), DIMENSION(:), INTENT(INOUT) :: mos
639 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s
640 : TYPE(preconditioner_p_type), DIMENSION(:), POINTER :: ot_preconditioner
641 : INTEGER, INTENT(IN) :: prec_type, solver_type
642 : REAL(dp), INTENT(IN) :: energy_gap
643 : INTEGER, INTENT(IN) :: nspins
644 : LOGICAL, INTENT(IN), OPTIONAL :: has_unit_metric, convert_to_dbcsr
645 : INTEGER, INTENT(IN), OPTIONAL :: chol_type
646 : LOGICAL, INTENT(IN), OPTIONAL :: full_mo_set
647 : INTEGER, INTENT(IN), OPTIONAL :: chebyshev_degree
648 :
649 : CHARACTER(LEN=*), PARAMETER :: routineN = 'prepare_preconditioner'
650 :
651 : CHARACTER(LEN=default_string_length) :: msg
652 : INTEGER :: handle, icall, ispin, &
653 : my_chebyshev_degree, n_loops
654 : INTEGER, DIMENSION(5) :: nocc, norb
655 : LOGICAL :: do_co_rotate, my_convert_to_dbcsr, &
656 : my_full_mo_set, my_has_unit_metric, &
657 : use_mo_coeff_b
658 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
659 : TYPE(cp_fm_type), POINTER :: mo_coeff
660 7443 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: kinetic
661 : TYPE(dbcsr_type), POINTER :: matrix_t, mo_coeff_b
662 : TYPE(dft_control_type), POINTER :: dft_control
663 : TYPE(mp_para_env_type), POINTER :: para_env
664 :
665 7443 : CALL timeset(routineN, handle)
666 7443 : NULLIFY (matrix_t, mo_coeff_b, mo_coeff, kinetic, dft_control, para_env, blacs_env)
667 7443 : my_has_unit_metric = .FALSE.
668 7443 : IF (PRESENT(has_unit_metric)) my_has_unit_metric = has_unit_metric
669 7443 : my_convert_to_dbcsr = .TRUE.
670 7443 : IF (PRESENT(convert_to_dbcsr)) my_convert_to_dbcsr = convert_to_dbcsr
671 7443 : my_full_mo_set = .FALSE.
672 7443 : IF (PRESENT(full_mo_set)) my_full_mo_set = full_mo_set
673 7443 : my_chebyshev_degree = 8
674 7443 : IF (PRESENT(chebyshev_degree)) my_chebyshev_degree = chebyshev_degree
675 :
676 : CALL get_qs_env(qs_env, &
677 : dft_control=dft_control, &
678 : para_env=para_env, &
679 7443 : blacs_env=blacs_env)
680 :
681 7443 : IF (dft_control%qs_control%semi_empirical .OR. dft_control%qs_control%dftb .OR. &
682 : dft_control%qs_control%xtb) THEN
683 2148 : IF (prec_type == ot_precond_full_kinetic) THEN
684 0 : msg = "Full_kinetic not available for semi-empirical methods"
685 0 : CPABORT(TRIM(msg))
686 : END IF
687 2148 : matrix_t => matrix_s(1)%matrix
688 : ELSE
689 5295 : CPASSERT(.NOT. my_has_unit_metric)
690 5295 : CALL get_qs_env(qs_env, kinetic=kinetic)
691 5295 : matrix_t => kinetic(1)%matrix
692 : END IF
693 :
694 : ! use full set of MOs or just occupied MOs
695 7443 : nocc = 0
696 7443 : norb = 0
697 7443 : IF (my_full_mo_set) THEN
698 38 : DO ispin = 1, nspins
699 20 : CALL get_mo_set(mo_set=mos(ispin), homo=nocc(ispin), nmo=norb(ispin))
700 38 : CALL set_mo_set(mo_set=mos(ispin), homo=norb(ispin))
701 : END DO
702 : END IF
703 : !determines how often make preconditioner is called, spin dependent methods have to be called twice
704 7443 : n_loops = 1
705 7443 : IF (prec_type == ot_precond_full_single_inverse) n_loops = nspins
706 : ! check whether we need the ev and rotate the MOs
707 2524 : SELECT CASE (prec_type)
708 : CASE (ot_precond_full_all)
709 : ! if one of these preconditioners is used every spin needs to call make_preconditioner
710 2524 : n_loops = nspins
711 :
712 2524 : do_co_rotate = ASSOCIATED(qs_env%mo_derivs)
713 10745 : DO ispin = 1, nspins
714 3302 : CALL get_mo_set(mo_set=mos(ispin), mo_coeff_b=mo_coeff_b, mo_coeff=mo_coeff)
715 3302 : use_mo_coeff_b = mos(ispin)%use_mo_coeff_b
716 5826 : IF (use_mo_coeff_b .AND. do_co_rotate) THEN
717 : CALL calculate_subspace_eigenvalues(mo_coeff_b, matrix_ks(ispin)%matrix, &
718 : do_rotation=.TRUE., &
719 : co_rotate=qs_env%mo_derivs(ispin)%matrix, &
720 : para_env=para_env, &
721 3282 : blacs_env=blacs_env)
722 20 : ELSE IF (use_mo_coeff_b) THEN
723 : CALL calculate_subspace_eigenvalues(mo_coeff_b, matrix_ks(ispin)%matrix, &
724 : do_rotation=.TRUE., &
725 : para_env=para_env, &
726 20 : blacs_env=blacs_env)
727 : ELSE
728 : CALL calculate_subspace_eigenvalues(mo_coeff, matrix_ks(ispin)%matrix, &
729 0 : do_rotation=.TRUE.)
730 : END IF
731 : END DO
732 : CASE DEFAULT
733 : ! No need to rotate the MOs
734 : END SELECT
735 :
736 : ! check whether we have a preconditioner
737 724 : SELECT CASE (prec_type)
738 : CASE (ot_precond_none)
739 1448 : DO ispin = 1, SIZE(ot_preconditioner)
740 1448 : ot_preconditioner(ispin)%preconditioner%in_use = 0
741 : END DO
742 : CASE DEFAULT
743 22006 : DO icall = 1, n_loops
744 7844 : ot_preconditioner(icall)%preconditioner%polynomial_degree = my_chebyshev_degree
745 14563 : IF (my_has_unit_metric) THEN
746 : CALL make_preconditioner(ot_preconditioner(icall)%preconditioner, &
747 : prec_type, &
748 : solver_type, &
749 : matrix_h=matrix_ks(icall)%matrix, &
750 : mo_set=mos(icall), &
751 : energy_gap=energy_gap, &
752 474 : convert_precond_to_dbcsr=my_convert_to_dbcsr)
753 : ELSE
754 : CALL make_preconditioner(ot_preconditioner(icall)%preconditioner, &
755 : prec_type, &
756 : solver_type, &
757 : matrix_h=matrix_ks(icall)%matrix, &
758 : matrix_s=matrix_s(1)%matrix, &
759 : matrix_t=matrix_t, &
760 : mo_set=mos(icall), &
761 : energy_gap=energy_gap, &
762 7370 : convert_precond_to_dbcsr=my_convert_to_dbcsr, chol_type=chol_type)
763 : END IF
764 : END DO
765 : END SELECT
766 :
767 : ! reset homo values
768 7443 : IF (my_full_mo_set) THEN
769 38 : DO ispin = 1, nspins
770 38 : CALL set_mo_set(mo_set=mos(ispin), homo=nocc(ispin))
771 : END DO
772 : END IF
773 :
774 7443 : CALL timestop(handle)
775 :
776 7443 : END SUBROUTINE prepare_preconditioner
777 :
778 : END MODULE preconditioner
|
