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_makes
16 : USE arnoldi_api, ONLY: arnoldi_env_type,&
17 : arnoldi_ev,&
18 : deallocate_arnoldi_env,&
19 : get_selected_ritz_val,&
20 : get_selected_ritz_vector,&
21 : set_arnoldi_initial_vector,&
22 : setup_arnoldi_env
23 : USE cp_cfm_basic_linalg, ONLY: cp_cfm_column_scale,&
24 : cp_cfm_gemm,&
25 : cp_cfm_scale_and_add,&
26 : cp_cfm_triangular_multiply,&
27 : cp_cfm_uplo_to_full
28 : USE cp_cfm_cholesky, ONLY: cp_cfm_cholesky_decompose,&
29 : cp_cfm_cholesky_invert
30 : USE cp_cfm_diag, ONLY: cp_cfm_geeig
31 : USE cp_cfm_types, ONLY: cp_cfm_create,&
32 : cp_cfm_get_info,&
33 : cp_cfm_release,&
34 : cp_cfm_set_element,&
35 : cp_cfm_to_cfm,&
36 : cp_cfm_type
37 : USE cp_dbcsr_api, ONLY: &
38 : dbcsr_add, dbcsr_copy, dbcsr_create, dbcsr_get_info, dbcsr_multiply, dbcsr_p_type, &
39 : dbcsr_release, dbcsr_type, dbcsr_type_symmetric
40 : USE cp_dbcsr_contrib, ONLY: dbcsr_add_on_diag
41 : USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
42 : cp_dbcsr_m_by_n_from_template,&
43 : cp_dbcsr_sm_fm_multiply,&
44 : cp_fm_to_dbcsr_row_template
45 : USE cp_fm_basic_linalg, ONLY: cp_fm_column_scale,&
46 : cp_fm_triangular_invert,&
47 : cp_fm_triangular_multiply,&
48 : cp_fm_uplo_to_full
49 : USE cp_fm_cholesky, ONLY: cp_fm_cholesky_decompose,&
50 : cp_fm_cholesky_reduce,&
51 : cp_fm_cholesky_restore
52 : USE cp_fm_diag, ONLY: choose_eigv_solver
53 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
54 : cp_fm_struct_release,&
55 : cp_fm_struct_type
56 : USE cp_fm_types, ONLY: cp_fm_create,&
57 : cp_fm_get_diag,&
58 : cp_fm_get_info,&
59 : cp_fm_release,&
60 : cp_fm_to_fm,&
61 : cp_fm_type
62 : USE input_constants, ONLY: &
63 : cholesky_inverse, cholesky_reduce, ot_precond_full_all, ot_precond_full_kinetic, &
64 : ot_precond_full_single, ot_precond_full_single_inverse, ot_precond_s_inverse, &
65 : ot_precond_solver_default, ot_precond_solver_inv_chol
66 : USE kinds, ONLY: dp
67 : USE mathconstants, ONLY: z_one,&
68 : z_zero
69 : USE parallel_gemm_api, ONLY: parallel_gemm
70 : USE preconditioner_types, ONLY: preconditioner_type
71 : #include "./base/base_uses.f90"
72 :
73 : IMPLICIT NONE
74 :
75 : PRIVATE
76 :
77 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'preconditioner_makes'
78 :
79 : PUBLIC :: make_complex_full_all, make_complex_full_kinetic, make_complex_full_s_inverse, &
80 : make_complex_full_single, make_complex_full_single_inverse, &
81 : make_preconditioner_matrix
82 :
83 : CONTAINS
84 :
85 : ! **************************************************************************************************
86 : !> \brief Build the state-selective FULL_ALL operator for a complex k-point channel.
87 : !> The occupied/reference subspace is retained in its current gauge, while the
88 : !> orthogonal complement carries the spectrum of H(k).
89 : !> \param preconditioner_env preconditioner storage
90 : !> \param matrix_c0 complex reference orbitals, C^H S C = I
91 : !> \param matrix_h complex Hermitian k-point Hamiltonian
92 : !> \param matrix_s complex Hermitian k-point overlap
93 : !> \param c0_evals reference-orbital energies in the current OT gauge
94 : !> \param energy_gap denominator floor
95 : ! **************************************************************************************************
96 152 : SUBROUTINE make_complex_full_all(preconditioner_env, matrix_c0, matrix_h, matrix_s, &
97 152 : c0_evals, energy_gap)
98 :
99 : TYPE(preconditioner_type) :: preconditioner_env
100 : TYPE(cp_cfm_type), INTENT(IN) :: matrix_c0, matrix_h, matrix_s
101 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: c0_evals
102 : REAL(KIND=dp), INTENT(IN) :: energy_gap
103 :
104 : CHARACTER(len=*), PARAMETER :: routineN = 'make_complex_full_all'
105 : REAL(KIND=dp), PARAMETER :: fudge_factor = 0.25_dp, &
106 : lambda_base = 10.0_dp
107 :
108 152 : COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:) :: shifted_evals
109 : COMPLEX(KIND=dp), CONTIGUOUS, DIMENSION(:, :), &
110 152 : POINTER :: local_data
111 : INTEGER :: handle, j, k, n, ncol_local, nrow_local
112 152 : INTEGER, DIMENSION(:), POINTER :: col_indices
113 : REAL(KIND=dp) :: error_estimate, lambda
114 152 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: residual_norm_sq
115 : TYPE(cp_cfm_type) :: matrix_chc, matrix_diag_work, matrix_hc0, matrix_s_copy, matrix_sc0, &
116 : matrix_sc_chc, matrix_shifted_sc, matrix_tmp
117 :
118 152 : CALL timeset(routineN, handle)
119 :
120 152 : CALL cp_cfm_get_info(matrix_c0, nrow_global=n, ncol_global=k)
121 152 : CPASSERT(n > 0)
122 152 : CPASSERT(k > 0 .AND. k <= n)
123 152 : CPASSERT(SIZE(c0_evals) >= k)
124 152 : CPASSERT(energy_gap > 0.0_dp)
125 :
126 152 : IF (ASSOCIATED(preconditioner_env%complex_fm)) THEN
127 0 : CALL cp_cfm_release(preconditioner_env%complex_fm)
128 0 : DEALLOCATE (preconditioner_env%complex_fm)
129 : END IF
130 152 : IF (ASSOCIATED(preconditioner_env%occ_evals)) DEALLOCATE (preconditioner_env%occ_evals)
131 152 : IF (ASSOCIATED(preconditioner_env%full_evals)) DEALLOCATE (preconditioner_env%full_evals)
132 :
133 152 : ALLOCATE (preconditioner_env%complex_fm)
134 : CALL cp_cfm_create(preconditioner_env%complex_fm, matrix_h%matrix_struct, &
135 152 : name='complex FULL_ALL eigenvectors')
136 760 : ALLOCATE (preconditioner_env%full_evals(n), preconditioner_env%occ_evals(k))
137 :
138 152 : CALL cp_cfm_create(matrix_hc0, matrix_c0%matrix_struct, name='complex FULL_ALL HC')
139 152 : CALL cp_cfm_create(matrix_sc0, matrix_c0%matrix_struct, name='complex FULL_ALL SC')
140 : CALL cp_cfm_create(matrix_chc, matrix_c0%matrix_struct, nrow=k, ncol=k, &
141 152 : name='complex FULL_ALL CHC')
142 : CALL cp_cfm_create(matrix_sc_chc, matrix_c0%matrix_struct, &
143 152 : name='complex FULL_ALL SC CHC')
144 : CALL cp_cfm_create(matrix_shifted_sc, matrix_c0%matrix_struct, &
145 152 : name='complex FULL_ALL shifted SC')
146 152 : CALL cp_cfm_create(matrix_tmp, matrix_h%matrix_struct, name='complex FULL_ALL projected H')
147 152 : CALL cp_cfm_create(matrix_s_copy, matrix_s%matrix_struct, name='complex FULL_ALL S copy')
148 : CALL cp_cfm_create(matrix_diag_work, matrix_h%matrix_struct, &
149 152 : name='complex FULL_ALL diagonalization work')
150 :
151 152 : CALL cp_cfm_gemm('N', 'N', n, k, n, z_one, matrix_h, matrix_c0, z_zero, matrix_hc0)
152 152 : CALL cp_cfm_gemm('N', 'N', n, k, n, z_one, matrix_s, matrix_c0, z_zero, matrix_sc0)
153 152 : CALL cp_cfm_gemm('C', 'N', k, k, n, z_one, matrix_c0, matrix_hc0, z_zero, matrix_chc)
154 :
155 : ! Estimate the S^-1 norm of the Ritz residual, R = H C - S C epsilon.
156 : ! As in the real FULL_ALL implementation, use it to prevent an inaccurate
157 : ! reference subspace from producing an overly aggressive preconditioner.
158 608 : ALLOCATE (shifted_evals(k), residual_norm_sq(k))
159 1562 : shifted_evals(:) = CMPLX(c0_evals(1:k), 0.0_dp, KIND=dp)
160 152 : CALL cp_cfm_to_cfm(matrix_sc0, matrix_sc_chc)
161 152 : CALL cp_cfm_column_scale(matrix_sc_chc, shifted_evals)
162 152 : CALL cp_cfm_to_cfm(matrix_hc0, matrix_shifted_sc)
163 152 : CALL cp_cfm_scale_and_add(z_one, matrix_shifted_sc, -z_one, matrix_sc_chc)
164 152 : CALL cp_cfm_to_cfm(matrix_s, matrix_s_copy)
165 152 : CALL cp_cfm_cholesky_decompose(matrix_s_copy)
166 : CALL cp_cfm_triangular_multiply(matrix_s_copy, matrix_shifted_sc, side='L', &
167 152 : transa_tr='C', invert_tr=.TRUE., uplo_tr='U')
168 :
169 152 : residual_norm_sq(:) = 0.0_dp
170 : CALL cp_cfm_get_info(matrix_shifted_sc, nrow_local=nrow_local, ncol_local=ncol_local, &
171 152 : col_indices=col_indices, local_data=local_data)
172 1562 : DO j = 1, ncol_local
173 : residual_norm_sq(col_indices(j)) = residual_norm_sq(col_indices(j)) + &
174 57430 : SUM(ABS(local_data(1:nrow_local, j))**2)
175 : END DO
176 152 : CALL preconditioner_env%para_env%sum(residual_norm_sq)
177 1562 : error_estimate = SQRT(MAXVAL(residual_norm_sq))
178 152 : preconditioner_env%energy_gap = MAX(energy_gap, error_estimate*fudge_factor)
179 152 : lambda = lambda_base + error_estimate
180 :
181 : ! Q^H H Q with Q = I - C C^H S. This removes occupied/complement cross terms
182 : ! without rotating the reference columns.
183 152 : CALL cp_cfm_to_cfm(matrix_h, matrix_tmp)
184 152 : CALL cp_cfm_gemm('N', 'C', n, n, k, -z_one, matrix_hc0, matrix_sc0, z_one, matrix_tmp)
185 152 : CALL cp_cfm_gemm('N', 'C', n, n, k, -z_one, matrix_sc0, matrix_hc0, z_one, matrix_tmp)
186 152 : CALL cp_cfm_gemm('N', 'N', n, k, k, z_one, matrix_sc0, matrix_chc, z_zero, matrix_sc_chc)
187 152 : CALL cp_cfm_gemm('N', 'C', n, n, k, z_one, matrix_sc_chc, matrix_sc0, z_one, matrix_tmp)
188 :
189 : ! Shift the retained reference subspace below the complementary spectrum, diagonalize,
190 : ! and then restore the original reference columns and their energy labels exactly.
191 1562 : shifted_evals(:) = CMPLX(c0_evals(1:k) - lambda, 0.0_dp, KIND=dp)
192 152 : CALL cp_cfm_to_cfm(matrix_sc0, matrix_shifted_sc)
193 152 : CALL cp_cfm_column_scale(matrix_shifted_sc, shifted_evals)
194 152 : CALL cp_cfm_gemm('N', 'C', n, n, k, z_one, matrix_shifted_sc, matrix_sc0, z_one, matrix_tmp)
195 :
196 152 : CALL cp_cfm_to_cfm(matrix_s, matrix_s_copy)
197 : CALL cp_cfm_geeig(matrix_tmp, matrix_s_copy, preconditioner_env%complex_fm, &
198 152 : preconditioner_env%full_evals, matrix_diag_work)
199 :
200 1562 : preconditioner_env%occ_evals(:) = c0_evals(1:k)
201 1562 : preconditioner_env%full_evals(1:k) = c0_evals(1:k)
202 152 : CALL cp_cfm_to_cfm(matrix_c0, preconditioner_env%complex_fm, k)
203 152 : preconditioner_env%in_use = ot_precond_full_all
204 152 : preconditioner_env%solver = ot_precond_solver_default
205 :
206 152 : DEALLOCATE (residual_norm_sq, shifted_evals)
207 152 : CALL cp_cfm_release(matrix_diag_work)
208 152 : CALL cp_cfm_release(matrix_s_copy)
209 152 : CALL cp_cfm_release(matrix_tmp)
210 152 : CALL cp_cfm_release(matrix_shifted_sc)
211 152 : CALL cp_cfm_release(matrix_sc_chc)
212 152 : CALL cp_cfm_release(matrix_chc)
213 152 : CALL cp_cfm_release(matrix_sc0)
214 152 : CALL cp_cfm_release(matrix_hc0)
215 :
216 152 : CALL timestop(handle)
217 :
218 456 : END SUBROUTINE make_complex_full_all
219 :
220 : ! **************************************************************************************************
221 : !> \brief Build the complex spectral FULL_SINGLE preconditioner.
222 : !> \param preconditioner_env preconditioner storage
223 : !> \param matrix_h complex Hermitian k-point Hamiltonian
224 : !> \param matrix_s complex Hermitian k-point overlap
225 : !> \param energy_homo occupied spectral edge
226 : !> \param energy_gap denominator floor
227 : ! **************************************************************************************************
228 28 : SUBROUTINE make_complex_full_single(preconditioner_env, matrix_h, matrix_s, &
229 : energy_homo, energy_gap)
230 :
231 : TYPE(preconditioner_type) :: preconditioner_env
232 : TYPE(cp_cfm_type), INTENT(IN) :: matrix_h, matrix_s
233 : REAL(KIND=dp), INTENT(IN) :: energy_homo, energy_gap
234 :
235 : CHARACTER(len=*), PARAMETER :: routineN = 'make_complex_full_single'
236 :
237 28 : COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:) :: scaling
238 : INTEGER :: handle, i, n
239 28 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eigenvalues
240 : TYPE(cp_cfm_type) :: matrix_diag_a, matrix_diag_b, &
241 : matrix_eigenvectors, &
242 : matrix_scaled_eigenvectors, matrix_work
243 :
244 28 : CALL timeset(routineN, handle)
245 :
246 28 : CALL cp_cfm_get_info(matrix_h, nrow_global=n)
247 28 : CPASSERT(n > 0)
248 28 : CPASSERT(energy_gap > 0.0_dp)
249 :
250 28 : IF (ASSOCIATED(preconditioner_env%complex_fm)) THEN
251 0 : CALL cp_cfm_release(preconditioner_env%complex_fm)
252 0 : DEALLOCATE (preconditioner_env%complex_fm)
253 : END IF
254 28 : IF (ASSOCIATED(preconditioner_env%occ_evals)) DEALLOCATE (preconditioner_env%occ_evals)
255 28 : IF (ASSOCIATED(preconditioner_env%full_evals)) DEALLOCATE (preconditioner_env%full_evals)
256 :
257 28 : ALLOCATE (preconditioner_env%complex_fm)
258 : CALL cp_cfm_create(preconditioner_env%complex_fm, matrix_h%matrix_struct, &
259 28 : name='complex FULL_SINGLE preconditioner')
260 : CALL cp_cfm_create(matrix_diag_a, matrix_h%matrix_struct, &
261 28 : name='complex FULL_SINGLE diagonalization A')
262 : CALL cp_cfm_create(matrix_diag_b, matrix_s%matrix_struct, &
263 28 : name='complex FULL_SINGLE diagonalization B')
264 : CALL cp_cfm_create(matrix_eigenvectors, matrix_h%matrix_struct, &
265 28 : name='complex FULL_SINGLE eigenvectors')
266 : CALL cp_cfm_create(matrix_scaled_eigenvectors, matrix_h%matrix_struct, &
267 28 : name='complex FULL_SINGLE scaled eigenvectors')
268 : CALL cp_cfm_create(matrix_work, matrix_h%matrix_struct, &
269 28 : name='complex FULL_SINGLE diagonalization work')
270 :
271 140 : ALLOCATE (eigenvalues(n), scaling(n))
272 28 : CALL cp_cfm_to_cfm(matrix_h, matrix_diag_a)
273 28 : CALL cp_cfm_to_cfm(matrix_s, matrix_diag_b)
274 : CALL cp_cfm_geeig(matrix_diag_a, matrix_diag_b, matrix_eigenvectors, &
275 28 : eigenvalues, matrix_work)
276 372 : DO i = 1, n
277 : scaling(i) = CMPLX(1.0_dp/MAX(eigenvalues(i) - energy_homo, energy_gap), &
278 372 : 0.0_dp, KIND=dp)
279 : END DO
280 28 : CALL cp_cfm_to_cfm(matrix_eigenvectors, matrix_scaled_eigenvectors)
281 28 : CALL cp_cfm_column_scale(matrix_scaled_eigenvectors, scaling)
282 : CALL cp_cfm_gemm('N', 'C', n, n, n, z_one, matrix_scaled_eigenvectors, &
283 28 : matrix_eigenvectors, z_zero, preconditioner_env%complex_fm)
284 :
285 28 : preconditioner_env%energy_gap = energy_gap
286 28 : preconditioner_env%in_use = ot_precond_full_single
287 28 : preconditioner_env%solver = ot_precond_solver_default
288 :
289 28 : DEALLOCATE (scaling, eigenvalues)
290 28 : CALL cp_cfm_release(matrix_work)
291 28 : CALL cp_cfm_release(matrix_scaled_eigenvectors)
292 28 : CALL cp_cfm_release(matrix_eigenvectors)
293 28 : CALL cp_cfm_release(matrix_diag_b)
294 28 : CALL cp_cfm_release(matrix_diag_a)
295 :
296 28 : CALL timestop(handle)
297 :
298 56 : END SUBROUTINE make_complex_full_single
299 :
300 : ! **************************************************************************************************
301 : !> \brief ...
302 : !> \param preconditioner_env ...
303 : !> \param matrix_h ...
304 : !> \param matrix_s ...
305 : !> \param matrix_t ...
306 : !> \param mo_coeff ...
307 : !> \param energy_homo ...
308 : !> \param eigenvalues_ot ...
309 : !> \param energy_gap ...
310 : !> \param my_solver_type ...
311 : ! **************************************************************************************************
312 9592 : SUBROUTINE make_preconditioner_matrix(preconditioner_env, matrix_h, matrix_s, matrix_t, mo_coeff, &
313 9592 : energy_homo, eigenvalues_ot, energy_gap, &
314 : my_solver_type)
315 : TYPE(preconditioner_type) :: preconditioner_env
316 : TYPE(dbcsr_type), POINTER :: matrix_h
317 : TYPE(dbcsr_type), OPTIONAL, POINTER :: matrix_s, matrix_t
318 : TYPE(cp_fm_type), INTENT(IN) :: mo_coeff
319 : REAL(KIND=dp) :: energy_homo
320 : REAL(KIND=dp), DIMENSION(:) :: eigenvalues_ot
321 : REAL(KIND=dp) :: energy_gap
322 : INTEGER :: my_solver_type
323 :
324 : INTEGER :: precon_type
325 :
326 9592 : precon_type = preconditioner_env%in_use
327 38 : SELECT CASE (precon_type)
328 : CASE (ot_precond_full_single)
329 38 : IF (my_solver_type /= ot_precond_solver_default) THEN
330 0 : CPABORT("Only PRECOND_SOLVER DEFAULT for the moment")
331 : END IF
332 38 : IF (PRESENT(matrix_s)) THEN
333 : CALL make_full_single(preconditioner_env, preconditioner_env%fm, &
334 32 : matrix_h, matrix_s, energy_homo, energy_gap)
335 : ELSE
336 : CALL make_full_single_ortho(preconditioner_env, preconditioner_env%fm, &
337 6 : matrix_h, energy_homo, energy_gap)
338 : END IF
339 :
340 : CASE (ot_precond_s_inverse)
341 78 : IF (my_solver_type == ot_precond_solver_default) my_solver_type = ot_precond_solver_inv_chol
342 78 : IF (.NOT. PRESENT(matrix_s)) THEN
343 0 : CPABORT("Type for S=1 not implemented")
344 : END IF
345 78 : CALL make_full_s_inverse(preconditioner_env, matrix_s)
346 :
347 : CASE (ot_precond_full_kinetic)
348 1321 : IF (my_solver_type == ot_precond_solver_default) my_solver_type = ot_precond_solver_inv_chol
349 1321 : IF (.NOT. (PRESENT(matrix_s) .AND. PRESENT(matrix_t))) THEN
350 0 : CPABORT("Type for S=1 not implemented")
351 : END IF
352 1321 : CALL make_full_kinetic(preconditioner_env, matrix_t, matrix_s, energy_gap)
353 : CASE (ot_precond_full_single_inverse)
354 4441 : IF (my_solver_type == ot_precond_solver_default) my_solver_type = ot_precond_solver_inv_chol
355 : CALL make_full_single_inverse(preconditioner_env, mo_coeff, matrix_h, energy_gap, &
356 4441 : matrix_s=matrix_s)
357 : CASE (ot_precond_full_all)
358 3714 : IF (my_solver_type /= ot_precond_solver_default) THEN
359 0 : CPABORT("Only PRECOND_SOLVER DEFAULT for the moment")
360 : END IF
361 3714 : IF (PRESENT(matrix_s)) THEN
362 : CALL make_full_all(preconditioner_env, mo_coeff, matrix_h, matrix_s, &
363 3634 : eigenvalues_ot, energy_gap)
364 : ELSE
365 : CALL make_full_all_ortho(preconditioner_env, mo_coeff, matrix_h, &
366 80 : eigenvalues_ot, energy_gap)
367 : END IF
368 :
369 : CASE DEFAULT
370 9592 : CPABORT("Type not implemented")
371 : END SELECT
372 :
373 9592 : END SUBROUTINE make_preconditioner_matrix
374 :
375 : ! **************************************************************************************************
376 : !> \brief Simply takes the overlap matrix as preconditioner
377 : !> \param preconditioner_env ...
378 : !> \param matrix_s ...
379 : ! **************************************************************************************************
380 78 : SUBROUTINE make_full_s_inverse(preconditioner_env, matrix_s)
381 : TYPE(preconditioner_type) :: preconditioner_env
382 : TYPE(dbcsr_type), POINTER :: matrix_s
383 :
384 : CHARACTER(len=*), PARAMETER :: routineN = 'make_full_s_inverse'
385 :
386 : INTEGER :: handle
387 :
388 78 : CALL timeset(routineN, handle)
389 :
390 78 : CPASSERT(ASSOCIATED(matrix_s))
391 :
392 78 : IF (.NOT. ASSOCIATED(preconditioner_env%sparse_matrix)) THEN
393 78 : ALLOCATE (preconditioner_env%sparse_matrix)
394 : END IF
395 78 : CALL dbcsr_copy(preconditioner_env%sparse_matrix, matrix_s, name="full_kinetic")
396 :
397 78 : CALL timestop(handle)
398 :
399 78 : END SUBROUTINE make_full_s_inverse
400 :
401 : ! **************************************************************************************************
402 : !> \brief kinetic matrix+shift*overlap as preconditioner. Cheap but could
403 : !> be better
404 : !> \param preconditioner_env ...
405 : !> \param matrix_t ...
406 : !> \param matrix_s ...
407 : !> \param energy_gap ...
408 : ! **************************************************************************************************
409 1321 : SUBROUTINE make_full_kinetic(preconditioner_env, matrix_t, matrix_s, &
410 : energy_gap)
411 : TYPE(preconditioner_type) :: preconditioner_env
412 : TYPE(dbcsr_type), POINTER :: matrix_t, matrix_s
413 : REAL(KIND=dp) :: energy_gap
414 :
415 : CHARACTER(len=*), PARAMETER :: routineN = 'make_full_kinetic'
416 :
417 : INTEGER :: handle
418 : REAL(KIND=dp) :: shift
419 :
420 1321 : CALL timeset(routineN, handle)
421 :
422 1321 : CPASSERT(ASSOCIATED(matrix_t))
423 1321 : CPASSERT(ASSOCIATED(matrix_s))
424 :
425 1321 : IF (.NOT. ASSOCIATED(preconditioner_env%sparse_matrix)) THEN
426 1319 : ALLOCATE (preconditioner_env%sparse_matrix)
427 : END IF
428 1321 : CALL dbcsr_copy(preconditioner_env%sparse_matrix, matrix_t, name="full_kinetic")
429 :
430 1321 : shift = MAX(0.0_dp, energy_gap)
431 :
432 : CALL dbcsr_add(preconditioner_env%sparse_matrix, matrix_s, &
433 1321 : alpha_scalar=1.0_dp, beta_scalar=shift)
434 :
435 1321 : CALL timestop(handle)
436 :
437 1321 : END SUBROUTINE make_full_kinetic
438 :
439 : ! **************************************************************************************************
440 : !> \brief full_single_preconditioner
441 : !> \param preconditioner_env ...
442 : !> \param fm ...
443 : !> \param matrix_h ...
444 : !> \param matrix_s ...
445 : !> \param energy_homo ...
446 : !> \param energy_gap ...
447 : ! **************************************************************************************************
448 32 : SUBROUTINE make_full_single(preconditioner_env, fm, matrix_h, matrix_s, &
449 : energy_homo, energy_gap)
450 : TYPE(preconditioner_type) :: preconditioner_env
451 : TYPE(cp_fm_type), POINTER :: fm
452 : TYPE(dbcsr_type), POINTER :: matrix_h, matrix_s
453 : REAL(KIND=dp) :: energy_homo, energy_gap
454 :
455 : CHARACTER(len=*), PARAMETER :: routineN = 'make_full_single'
456 :
457 : INTEGER :: handle, i, n
458 32 : REAL(KIND=dp), DIMENSION(:), POINTER :: evals
459 : TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp
460 : TYPE(cp_fm_type) :: fm_h, fm_s
461 :
462 32 : CALL timeset(routineN, handle)
463 :
464 32 : NULLIFY (fm_struct_tmp, evals)
465 :
466 32 : IF (ASSOCIATED(fm)) THEN
467 0 : CALL cp_fm_release(fm)
468 0 : DEALLOCATE (fm)
469 : NULLIFY (fm)
470 : END IF
471 32 : CALL dbcsr_get_info(matrix_h, nfullrows_total=n)
472 96 : ALLOCATE (evals(n))
473 :
474 : CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=n, ncol_global=n, &
475 : context=preconditioner_env%ctxt, &
476 32 : para_env=preconditioner_env%para_env)
477 32 : ALLOCATE (fm)
478 32 : CALL cp_fm_create(fm, fm_struct_tmp, name="preconditioner")
479 32 : CALL cp_fm_create(fm_h, fm_struct_tmp, name="fm_h")
480 32 : CALL cp_fm_create(fm_s, fm_struct_tmp, name="fm_s")
481 32 : CALL cp_fm_struct_release(fm_struct_tmp)
482 :
483 32 : CALL copy_dbcsr_to_fm(matrix_h, fm_h)
484 32 : CALL copy_dbcsr_to_fm(matrix_s, fm_s)
485 32 : CALL cp_fm_cholesky_decompose(fm_s)
486 :
487 32 : SELECT CASE (preconditioner_env%cholesky_use)
488 : CASE (cholesky_inverse)
489 : ! if cho inverse
490 0 : CALL cp_fm_triangular_invert(fm_s)
491 0 : CALL cp_fm_uplo_to_full(fm_h, fm)
492 :
493 : CALL cp_fm_triangular_multiply(fm_s, fm_h, side="R", transpose_tr=.FALSE., &
494 0 : invert_tr=.FALSE., uplo_tr="U", n_rows=n, n_cols=n, alpha=1.0_dp)
495 : CALL cp_fm_triangular_multiply(fm_s, fm_h, side="L", transpose_tr=.TRUE., &
496 0 : invert_tr=.FALSE., uplo_tr="U", n_rows=n, n_cols=n, alpha=1.0_dp)
497 : CASE (cholesky_reduce)
498 32 : CALL cp_fm_cholesky_reduce(fm_h, fm_s)
499 : CASE DEFAULT
500 32 : CPABORT("cholesky type not implemented")
501 : END SELECT
502 :
503 32 : CALL choose_eigv_solver(fm_h, fm, evals)
504 :
505 32 : SELECT CASE (preconditioner_env%cholesky_use)
506 : CASE (cholesky_inverse)
507 : CALL cp_fm_triangular_multiply(fm_s, fm, side="L", transpose_tr=.FALSE., &
508 0 : invert_tr=.FALSE., uplo_tr="U", n_rows=n, n_cols=n, alpha=1.0_dp)
509 0 : DO i = 1, n
510 0 : evals(i) = 1.0_dp/MAX(evals(i) - energy_homo, energy_gap)
511 : END DO
512 0 : CALL cp_fm_to_fm(fm, fm_h)
513 : CASE (cholesky_reduce)
514 32 : CALL cp_fm_cholesky_restore(fm, n, fm_s, fm_h, "SOLVE")
515 568 : DO i = 1, n
516 568 : evals(i) = 1.0_dp/MAX(evals(i) - energy_homo, energy_gap)
517 : END DO
518 64 : CALL cp_fm_to_fm(fm_h, fm)
519 : END SELECT
520 :
521 32 : CALL cp_fm_column_scale(fm, evals)
522 32 : CALL parallel_gemm('N', 'T', n, n, n, 1.0_dp, fm, fm_h, 0.0_dp, fm_s)
523 32 : CALL cp_fm_to_fm(fm_s, fm)
524 :
525 32 : DEALLOCATE (evals)
526 32 : CALL cp_fm_release(fm_h)
527 32 : CALL cp_fm_release(fm_s)
528 :
529 32 : CALL timestop(handle)
530 :
531 96 : END SUBROUTINE make_full_single
532 :
533 : ! **************************************************************************************************
534 : !> \brief full single in the orthonormal basis
535 : !> \param preconditioner_env ...
536 : !> \param fm ...
537 : !> \param matrix_h ...
538 : !> \param energy_homo ...
539 : !> \param energy_gap ...
540 : ! **************************************************************************************************
541 6 : SUBROUTINE make_full_single_ortho(preconditioner_env, fm, matrix_h, &
542 : energy_homo, energy_gap)
543 : TYPE(preconditioner_type) :: preconditioner_env
544 : TYPE(cp_fm_type), POINTER :: fm
545 : TYPE(dbcsr_type), POINTER :: matrix_h
546 : REAL(KIND=dp) :: energy_homo, energy_gap
547 :
548 : CHARACTER(len=*), PARAMETER :: routineN = 'make_full_single_ortho'
549 :
550 : INTEGER :: handle, i, n
551 6 : REAL(KIND=dp), DIMENSION(:), POINTER :: evals
552 : TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp
553 : TYPE(cp_fm_type) :: fm_h, fm_s
554 :
555 6 : CALL timeset(routineN, handle)
556 6 : NULLIFY (fm_struct_tmp, evals)
557 :
558 6 : IF (ASSOCIATED(fm)) THEN
559 0 : CALL cp_fm_release(fm)
560 0 : DEALLOCATE (fm)
561 : NULLIFY (fm)
562 : END IF
563 6 : CALL dbcsr_get_info(matrix_h, nfullrows_total=n)
564 18 : ALLOCATE (evals(n))
565 :
566 : CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=n, ncol_global=n, &
567 : context=preconditioner_env%ctxt, &
568 6 : para_env=preconditioner_env%para_env)
569 6 : ALLOCATE (fm)
570 6 : CALL cp_fm_create(fm, fm_struct_tmp, name="preconditioner")
571 6 : CALL cp_fm_create(fm_h, fm_struct_tmp, name="fm_h")
572 6 : CALL cp_fm_create(fm_s, fm_struct_tmp, name="fm_s")
573 6 : CALL cp_fm_struct_release(fm_struct_tmp)
574 :
575 6 : CALL copy_dbcsr_to_fm(matrix_h, fm_h)
576 :
577 6 : CALL choose_eigv_solver(fm_h, fm, evals)
578 282 : DO i = 1, n
579 282 : evals(i) = 1.0_dp/MAX(evals(i) - energy_homo, energy_gap)
580 : END DO
581 6 : CALL cp_fm_to_fm(fm, fm_h)
582 6 : CALL cp_fm_column_scale(fm, evals)
583 6 : CALL parallel_gemm('N', 'T', n, n, n, 1.0_dp, fm, fm_h, 0.0_dp, fm_s)
584 6 : CALL cp_fm_to_fm(fm_s, fm)
585 :
586 6 : DEALLOCATE (evals)
587 6 : CALL cp_fm_release(fm_h)
588 6 : CALL cp_fm_release(fm_s)
589 :
590 6 : CALL timestop(handle)
591 :
592 18 : END SUBROUTINE make_full_single_ortho
593 :
594 : ! **************************************************************************************************
595 : !> \brief generates a state by state preconditioner based on the full hamiltonian matrix
596 : !> \param preconditioner_env ...
597 : !> \param matrix_c0 ...
598 : !> \param matrix_h ...
599 : !> \param matrix_s ...
600 : !> \param c0_evals ...
601 : !> \param energy_gap should be a slight underestimate of the physical energy gap for almost all systems
602 : !> the c0 are already ritz states of (h,s)
603 : !> \par History
604 : !> 10.2006 made more stable [Joost VandeVondele]
605 : !> \note
606 : !> includes error estimate on the hamiltonian matrix to result in a stable preconditioner
607 : !> a preconditioner for each eigenstate i is generated by keeping the factorized form
608 : !> U diag( something i ) U^T. It is important to only precondition in the subspace orthogonal to c0.
609 : !> not only is it the only part that matters, it also simplifies the computation of
610 : !> the lagrangian multipliers in the OT minimization (i.e. if the c0 here is different
611 : !> from the c0 used in the OT setup, there will be a bug).
612 : ! **************************************************************************************************
613 3634 : SUBROUTINE make_full_all(preconditioner_env, matrix_c0, matrix_h, matrix_s, c0_evals, energy_gap)
614 : TYPE(preconditioner_type) :: preconditioner_env
615 : TYPE(cp_fm_type), INTENT(IN) :: matrix_c0
616 : TYPE(dbcsr_type), POINTER :: matrix_h, matrix_s
617 : REAL(KIND=dp), DIMENSION(:) :: c0_evals
618 : REAL(KIND=dp) :: energy_gap
619 :
620 : CHARACTER(len=*), PARAMETER :: routineN = 'make_full_all'
621 : REAL(KIND=dp), PARAMETER :: fudge_factor = 0.25_dp, &
622 : lambda_base = 10.0_dp
623 :
624 : INTEGER :: handle, k, n
625 : REAL(KIND=dp) :: error_estimate, lambda
626 3634 : REAL(KIND=dp), DIMENSION(:), POINTER :: diag, norms, shifted_evals
627 : TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp
628 : TYPE(cp_fm_type) :: matrix_hc0, matrix_left, matrix_s1, &
629 : matrix_s2, matrix_sc0, matrix_shc0, &
630 : matrix_tmp, ortho
631 : TYPE(cp_fm_type), POINTER :: matrix_pre
632 :
633 3634 : CALL timeset(routineN, handle)
634 :
635 3634 : IF (ASSOCIATED(preconditioner_env%fm)) THEN
636 0 : CALL cp_fm_release(preconditioner_env%fm)
637 0 : DEALLOCATE (preconditioner_env%fm)
638 : NULLIFY (preconditioner_env%fm)
639 : END IF
640 3634 : CALL cp_fm_get_info(matrix_c0, nrow_global=n, ncol_global=k)
641 : CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=n, ncol_global=n, &
642 : context=preconditioner_env%ctxt, &
643 3634 : para_env=preconditioner_env%para_env)
644 3634 : ALLOCATE (preconditioner_env%fm)
645 3634 : CALL cp_fm_create(preconditioner_env%fm, fm_struct_tmp, name="preconditioner_env%fm")
646 3634 : CALL cp_fm_create(ortho, fm_struct_tmp, name="ortho")
647 3634 : CALL cp_fm_create(matrix_tmp, fm_struct_tmp, name="matrix_tmp")
648 3634 : CALL cp_fm_struct_release(fm_struct_tmp)
649 10902 : ALLOCATE (preconditioner_env%full_evals(n))
650 10792 : ALLOCATE (preconditioner_env%occ_evals(k))
651 :
652 : ! 0) cholesky decompose the overlap matrix, if this fails the basis is singular,
653 : ! more than EPS_DEFAULT
654 3634 : CALL copy_dbcsr_to_fm(matrix_s, ortho)
655 3634 : CALL cp_fm_cholesky_decompose(ortho)
656 : ! if cho inverse
657 3634 : IF (preconditioner_env%cholesky_use == cholesky_inverse) THEN
658 0 : CALL cp_fm_triangular_invert(ortho)
659 : END IF
660 : ! 1) Construct a new H matrix, which has the current C0 as eigenvectors,
661 : ! possibly shifted by an amount lambda,
662 : ! and the same spectrum as the original H matrix in the space orthogonal to the C0
663 : ! with P=C0 C0 ^ T
664 : ! (1 - PS)^T H (1-PS) + (PS)^T (H - lambda S ) (PS)
665 : ! we exploit that the C0 are already the ritz states of H
666 3634 : CALL cp_fm_create(matrix_sc0, matrix_c0%matrix_struct, name="sc0")
667 3634 : CALL cp_dbcsr_sm_fm_multiply(matrix_s, matrix_c0, matrix_sc0, k)
668 3634 : CALL cp_fm_create(matrix_hc0, matrix_c0%matrix_struct, name="hc0")
669 3634 : CALL cp_dbcsr_sm_fm_multiply(matrix_h, matrix_c0, matrix_hc0, k)
670 :
671 : ! An aside, try to estimate the error on the ritz values, we'll need it later on
672 3634 : CALL cp_fm_create(matrix_shc0, matrix_c0%matrix_struct, name="shc0")
673 :
674 3634 : SELECT CASE (preconditioner_env%cholesky_use)
675 : CASE (cholesky_inverse)
676 : ! if cho inverse
677 0 : CALL cp_fm_to_fm(matrix_hc0, matrix_shc0)
678 : CALL cp_fm_triangular_multiply(ortho, matrix_shc0, side="L", transpose_tr=.TRUE., &
679 0 : invert_tr=.FALSE., uplo_tr="U", n_rows=n, n_cols=k, alpha=1.0_dp)
680 : CASE (cholesky_reduce)
681 3634 : CALL cp_fm_cholesky_restore(matrix_hc0, k, ortho, matrix_shc0, "SOLVE", transa="T")
682 : CASE DEFAULT
683 3634 : CPABORT("cholesky type not implemented")
684 : END SELECT
685 : CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=k, ncol_global=k, &
686 : context=preconditioner_env%ctxt, &
687 3634 : para_env=preconditioner_env%para_env)
688 3634 : CALL cp_fm_create(matrix_s1, fm_struct_tmp, name="matrix_s1")
689 3634 : CALL cp_fm_struct_release(fm_struct_tmp)
690 : ! since we only use diagonal elements this is a bit of a waste
691 3634 : CALL parallel_gemm('T', 'N', k, k, n, 1.0_dp, matrix_shc0, matrix_shc0, 0.0_dp, matrix_s1)
692 7158 : ALLOCATE (diag(k))
693 3634 : CALL cp_fm_get_diag(matrix_s1, diag)
694 20918 : error_estimate = MAXVAL(SQRT(ABS(diag - c0_evals**2)))
695 3634 : DEALLOCATE (diag)
696 3634 : CALL cp_fm_release(matrix_s1)
697 3634 : CALL cp_fm_release(matrix_shc0)
698 : ! we'll only use the energy gap, if our estimate of the error on the eigenvalues
699 : ! is small enough. A large error combined with a small energy gap would otherwise lead to
700 : ! an aggressive but bad preconditioner. Only when the error is small (MD), we can precondition
701 : ! aggressively
702 3634 : preconditioner_env%energy_gap = MAX(energy_gap, error_estimate*fudge_factor)
703 3634 : CALL copy_dbcsr_to_fm(matrix_h, matrix_tmp)
704 3634 : matrix_pre => preconditioner_env%fm
705 3634 : CALL cp_fm_uplo_to_full(matrix_tmp, matrix_pre)
706 : ! tmp = H ( 1 - PS )
707 3634 : CALL parallel_gemm('N', 'T', n, n, k, -1.0_dp, matrix_hc0, matrix_sc0, 1.0_dp, matrix_tmp)
708 :
709 : CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=k, ncol_global=n, &
710 : context=preconditioner_env%ctxt, &
711 3634 : para_env=preconditioner_env%para_env)
712 3634 : CALL cp_fm_create(matrix_left, fm_struct_tmp, name="matrix_left")
713 3634 : CALL cp_fm_struct_release(fm_struct_tmp)
714 3634 : CALL parallel_gemm('T', 'N', k, n, n, 1.0_dp, matrix_c0, matrix_tmp, 0.0_dp, matrix_left)
715 : ! tmp = (1 - PS)^T H (1-PS)
716 3634 : CALL parallel_gemm('N', 'N', n, n, k, -1.0_dp, matrix_sc0, matrix_left, 1.0_dp, matrix_tmp)
717 3634 : CALL cp_fm_release(matrix_left)
718 :
719 7158 : ALLOCATE (shifted_evals(k))
720 3634 : lambda = lambda_base + error_estimate
721 20808 : shifted_evals = c0_evals - lambda
722 3634 : CALL cp_fm_to_fm(matrix_sc0, matrix_hc0)
723 3634 : CALL cp_fm_column_scale(matrix_hc0, shifted_evals)
724 3634 : CALL parallel_gemm('N', 'T', n, n, k, 1.0_dp, matrix_hc0, matrix_sc0, 1.0_dp, matrix_tmp)
725 :
726 : ! 2) diagonalize this operator
727 3634 : SELECT CASE (preconditioner_env%cholesky_use)
728 : CASE (cholesky_inverse)
729 : CALL cp_fm_triangular_multiply(ortho, matrix_tmp, side="R", transpose_tr=.FALSE., &
730 0 : invert_tr=.FALSE., uplo_tr="U", n_rows=n, n_cols=n, alpha=1.0_dp)
731 : CALL cp_fm_triangular_multiply(ortho, matrix_tmp, side="L", transpose_tr=.TRUE., &
732 0 : invert_tr=.FALSE., uplo_tr="U", n_rows=n, n_cols=n, alpha=1.0_dp)
733 : CASE (cholesky_reduce)
734 3634 : CALL cp_fm_cholesky_reduce(matrix_tmp, ortho)
735 : END SELECT
736 3634 : CALL choose_eigv_solver(matrix_tmp, matrix_pre, preconditioner_env%full_evals)
737 3634 : SELECT CASE (preconditioner_env%cholesky_use)
738 : CASE (cholesky_inverse)
739 : CALL cp_fm_triangular_multiply(ortho, matrix_pre, side="L", transpose_tr=.FALSE., &
740 0 : invert_tr=.FALSE., uplo_tr="U", n_rows=n, n_cols=n, alpha=1.0_dp)
741 0 : CALL cp_fm_to_fm(matrix_pre, matrix_tmp)
742 : CASE (cholesky_reduce)
743 3634 : CALL cp_fm_cholesky_restore(matrix_pre, n, ortho, matrix_tmp, "SOLVE")
744 7268 : CALL cp_fm_to_fm(matrix_tmp, matrix_pre)
745 : END SELECT
746 :
747 : ! test that the subspace remained conserved
748 : IF (.FALSE.) THEN
749 : CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=k, ncol_global=k, &
750 : context=preconditioner_env%ctxt, &
751 : para_env=preconditioner_env%para_env)
752 : CALL cp_fm_create(matrix_s1, fm_struct_tmp, name="matrix_s1")
753 : CALL cp_fm_create(matrix_s2, fm_struct_tmp, name="matrix_s2")
754 : CALL cp_fm_struct_release(fm_struct_tmp)
755 : ALLOCATE (norms(k))
756 : CALL parallel_gemm('T', 'N', k, k, n, 1.0_dp, matrix_sc0, matrix_tmp, 0.0_dp, matrix_s1)
757 : CALL choose_eigv_solver(matrix_s1, matrix_s2, norms)
758 : WRITE (*, *) "matrix norm deviation (should be close to zero): ", MAXVAL(ABS(ABS(norms) - 1.0_dp))
759 : DEALLOCATE (norms)
760 : CALL cp_fm_release(matrix_s1)
761 : CALL cp_fm_release(matrix_s2)
762 : END IF
763 :
764 : ! 3) replace the lowest k evals and evecs with what they should be
765 20808 : preconditioner_env%occ_evals = c0_evals
766 : ! notice, this choice causes the preconditioner to be constant when applied to sc0 (see apply_full_all)
767 20808 : preconditioner_env%full_evals(1:k) = c0_evals
768 3634 : CALL cp_fm_to_fm(matrix_c0, matrix_pre, k, 1, 1)
769 :
770 3634 : CALL cp_fm_release(matrix_sc0)
771 3634 : CALL cp_fm_release(matrix_hc0)
772 3634 : CALL cp_fm_release(ortho)
773 3634 : CALL cp_fm_release(matrix_tmp)
774 3634 : DEALLOCATE (shifted_evals)
775 3634 : CALL timestop(handle)
776 :
777 29072 : END SUBROUTINE make_full_all
778 :
779 : ! **************************************************************************************************
780 : !> \brief full all in the orthonormal basis
781 : !> \param preconditioner_env ...
782 : !> \param matrix_c0 ...
783 : !> \param matrix_h ...
784 : !> \param c0_evals ...
785 : !> \param energy_gap ...
786 : ! **************************************************************************************************
787 80 : SUBROUTINE make_full_all_ortho(preconditioner_env, matrix_c0, matrix_h, c0_evals, energy_gap)
788 :
789 : TYPE(preconditioner_type) :: preconditioner_env
790 : TYPE(cp_fm_type), INTENT(IN) :: matrix_c0
791 : TYPE(dbcsr_type), POINTER :: matrix_h
792 : REAL(KIND=dp), DIMENSION(:) :: c0_evals
793 : REAL(KIND=dp) :: energy_gap
794 :
795 : CHARACTER(len=*), PARAMETER :: routineN = 'make_full_all_ortho'
796 : REAL(KIND=dp), PARAMETER :: fudge_factor = 0.25_dp, &
797 : lambda_base = 10.0_dp
798 :
799 : INTEGER :: handle, k, n
800 : REAL(KIND=dp) :: error_estimate, lambda
801 80 : REAL(KIND=dp), DIMENSION(:), POINTER :: diag, norms, shifted_evals
802 : TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp
803 : TYPE(cp_fm_type) :: matrix_hc0, matrix_left, matrix_s1, &
804 : matrix_s2, matrix_sc0, matrix_tmp
805 : TYPE(cp_fm_type), POINTER :: matrix_pre
806 :
807 80 : CALL timeset(routineN, handle)
808 :
809 80 : IF (ASSOCIATED(preconditioner_env%fm)) THEN
810 0 : CALL cp_fm_release(preconditioner_env%fm)
811 0 : DEALLOCATE (preconditioner_env%fm)
812 : NULLIFY (preconditioner_env%fm)
813 : END IF
814 80 : CALL cp_fm_get_info(matrix_c0, nrow_global=n, ncol_global=k)
815 : CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=n, ncol_global=n, &
816 : context=preconditioner_env%ctxt, &
817 80 : para_env=preconditioner_env%para_env)
818 80 : ALLOCATE (preconditioner_env%fm)
819 80 : CALL cp_fm_create(preconditioner_env%fm, fm_struct_tmp, name="preconditioner_env%fm")
820 80 : CALL cp_fm_create(matrix_tmp, fm_struct_tmp, name="matrix_tmp")
821 80 : CALL cp_fm_struct_release(fm_struct_tmp)
822 240 : ALLOCATE (preconditioner_env%full_evals(n))
823 240 : ALLOCATE (preconditioner_env%occ_evals(k))
824 :
825 : ! 1) Construct a new H matrix, which has the current C0 as eigenvectors,
826 : ! possibly shifted by an amount lambda,
827 : ! and the same spectrum as the original H matrix in the space orthogonal to the C0
828 : ! with P=C0 C0 ^ T
829 : ! (1 - PS)^T H (1-PS) + (PS)^T (H - lambda S ) (PS)
830 : ! we exploit that the C0 are already the ritz states of H
831 80 : CALL cp_fm_create(matrix_sc0, matrix_c0%matrix_struct, name="sc0")
832 80 : CALL cp_fm_to_fm(matrix_c0, matrix_sc0)
833 80 : CALL cp_fm_create(matrix_hc0, matrix_c0%matrix_struct, name="hc0")
834 80 : CALL cp_dbcsr_sm_fm_multiply(matrix_h, matrix_c0, matrix_hc0, k)
835 :
836 : ! An aside, try to estimate the error on the ritz values, we'll need it later on
837 : CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=k, ncol_global=k, &
838 : context=preconditioner_env%ctxt, &
839 80 : para_env=preconditioner_env%para_env)
840 80 : CALL cp_fm_create(matrix_s1, fm_struct_tmp, name="matrix_s1")
841 80 : CALL cp_fm_struct_release(fm_struct_tmp)
842 : ! since we only use diagonal elements this is a bit of a waste
843 80 : CALL parallel_gemm('T', 'N', k, k, n, 1.0_dp, matrix_hc0, matrix_hc0, 0.0_dp, matrix_s1)
844 160 : ALLOCATE (diag(k))
845 80 : CALL cp_fm_get_diag(matrix_s1, diag)
846 826 : error_estimate = MAXVAL(SQRT(ABS(diag - c0_evals**2)))
847 80 : DEALLOCATE (diag)
848 80 : CALL cp_fm_release(matrix_s1)
849 : ! we'll only use the energy gap, if our estimate of the error on the eigenvalues
850 : ! is small enough. A large error combined with a small energy gap would otherwise lead to
851 : ! an aggressive but bad preconditioner. Only when the error is small (MD), we can precondition
852 : ! aggressively
853 80 : preconditioner_env%energy_gap = MAX(energy_gap, error_estimate*fudge_factor)
854 :
855 80 : matrix_pre => preconditioner_env%fm
856 80 : CALL copy_dbcsr_to_fm(matrix_h, matrix_tmp)
857 80 : CALL cp_fm_uplo_to_full(matrix_tmp, matrix_pre)
858 : ! tmp = H ( 1 - PS )
859 80 : CALL parallel_gemm('N', 'T', n, n, k, -1.0_dp, matrix_hc0, matrix_sc0, 1.0_dp, matrix_tmp)
860 :
861 : CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=k, ncol_global=n, &
862 : context=preconditioner_env%ctxt, &
863 80 : para_env=preconditioner_env%para_env)
864 80 : CALL cp_fm_create(matrix_left, fm_struct_tmp, name="matrix_left")
865 80 : CALL cp_fm_struct_release(fm_struct_tmp)
866 80 : CALL parallel_gemm('T', 'N', k, n, n, 1.0_dp, matrix_c0, matrix_tmp, 0.0_dp, matrix_left)
867 : ! tmp = (1 - PS)^T H (1-PS)
868 80 : CALL parallel_gemm('N', 'N', n, n, k, -1.0_dp, matrix_sc0, matrix_left, 1.0_dp, matrix_tmp)
869 80 : CALL cp_fm_release(matrix_left)
870 :
871 160 : ALLOCATE (shifted_evals(k))
872 80 : lambda = lambda_base + error_estimate
873 826 : shifted_evals = c0_evals - lambda
874 80 : CALL cp_fm_to_fm(matrix_sc0, matrix_hc0)
875 80 : CALL cp_fm_column_scale(matrix_hc0, shifted_evals)
876 80 : CALL parallel_gemm('N', 'T', n, n, k, 1.0_dp, matrix_hc0, matrix_sc0, 1.0_dp, matrix_tmp)
877 :
878 : ! 2) diagonalize this operator
879 80 : CALL choose_eigv_solver(matrix_tmp, matrix_pre, preconditioner_env%full_evals)
880 :
881 : ! test that the subspace remained conserved
882 : IF (.FALSE.) THEN
883 : CALL cp_fm_to_fm(matrix_pre, matrix_tmp)
884 : CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=k, ncol_global=k, &
885 : context=preconditioner_env%ctxt, &
886 : para_env=preconditioner_env%para_env)
887 : CALL cp_fm_create(matrix_s1, fm_struct_tmp, name="matrix_s1")
888 : CALL cp_fm_create(matrix_s2, fm_struct_tmp, name="matrix_s2")
889 : CALL cp_fm_struct_release(fm_struct_tmp)
890 : ALLOCATE (norms(k))
891 : CALL parallel_gemm('T', 'N', k, k, n, 1.0_dp, matrix_sc0, matrix_tmp, 0.0_dp, matrix_s1)
892 : CALL choose_eigv_solver(matrix_s1, matrix_s2, norms)
893 :
894 : WRITE (*, *) "matrix norm deviation (should be close to zero): ", MAXVAL(ABS(ABS(norms) - 1.0_dp))
895 : DEALLOCATE (norms)
896 : CALL cp_fm_release(matrix_s1)
897 : CALL cp_fm_release(matrix_s2)
898 : END IF
899 :
900 : ! 3) replace the lowest k evals and evecs with what they should be
901 826 : preconditioner_env%occ_evals = c0_evals
902 : ! notice, this choice causes the preconditioner to be constant when applied to sc0 (see apply_full_all)
903 826 : preconditioner_env%full_evals(1:k) = c0_evals
904 80 : CALL cp_fm_to_fm(matrix_c0, matrix_pre, k, 1, 1)
905 :
906 80 : CALL cp_fm_release(matrix_sc0)
907 80 : CALL cp_fm_release(matrix_hc0)
908 80 : CALL cp_fm_release(matrix_tmp)
909 80 : DEALLOCATE (shifted_evals)
910 :
911 80 : CALL timestop(handle)
912 :
913 560 : END SUBROUTINE make_full_all_ortho
914 :
915 : ! **************************************************************************************************
916 : !> \brief generates a preconditioner matrix H-lambda S+(SC)(2.0*CT*H*C+delta)(SC)^T
917 : !> for later inversion.
918 : !> H is the Kohn Sham matrix
919 : !> lambda*S shifts the spectrum of the generalized form up by lambda
920 : !> the last term only shifts the occupied space (reversing them in energy order)
921 : !> This form is implicitly multiplied from both sides by S^0.5
922 : !> This ensures we precondition the correct quantity
923 : !> Before this reads S^-0.5 H S^-0.5 + lambda + (S^0.5 C)shifts(S^0.5 C)T
924 : !> which might be a bit more obvious
925 : !> Replaced the old full_single_inverse at revision 14616
926 : !> \param preconditioner_env the preconditioner env
927 : !> \param matrix_c0 the MO coefficient matrix (fm)
928 : !> \param matrix_h Kohn-Sham matrix (dbcsr)
929 : !> \param energy_gap an additional shift in lambda=-E_homo+energy_gap
930 : !> \param matrix_s the overlap matrix if not orthonormal (dbcsr, optional)
931 : ! **************************************************************************************************
932 4441 : SUBROUTINE make_full_single_inverse(preconditioner_env, matrix_c0, matrix_h, energy_gap, matrix_s)
933 : TYPE(preconditioner_type) :: preconditioner_env
934 : TYPE(cp_fm_type), INTENT(IN) :: matrix_c0
935 : TYPE(dbcsr_type), POINTER :: matrix_h
936 : REAL(KIND=dp) :: energy_gap
937 : TYPE(dbcsr_type), OPTIONAL, POINTER :: matrix_s
938 :
939 : CHARACTER(len=*), PARAMETER :: routineN = 'make_full_single_inverse'
940 :
941 : INTEGER :: handle, k, n
942 : REAL(KIND=dp) :: max_ev, min_ev, pre_shift
943 : TYPE(arnoldi_env_type) :: arnoldi_env
944 4441 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrices
945 : TYPE(dbcsr_type), TARGET :: dbcsr_cThc, dbcsr_hc, dbcsr_sc, mo_dbcsr
946 :
947 4441 : CALL timeset(routineN, handle)
948 :
949 : ! Allocate all working matrices needed
950 4441 : CALL cp_fm_get_info(matrix_c0, nrow_global=n, ncol_global=k)
951 : ! copy the fm MO's to a sparse matrix, can be solved better if the sparse version is already present
952 : ! but for the time beeing this will do
953 4441 : CALL cp_fm_to_dbcsr_row_template(mo_dbcsr, matrix_c0, matrix_h)
954 4441 : CALL dbcsr_create(dbcsr_sc, template=mo_dbcsr)
955 4441 : CALL dbcsr_create(dbcsr_hc, template=mo_dbcsr)
956 4441 : CALL cp_dbcsr_m_by_n_from_template(dbcsr_cThc, matrix_h, k, k, sym=dbcsr_type_symmetric)
957 :
958 : ! Check whether the output matrix was already created, if not do it now
959 4441 : IF (.NOT. ASSOCIATED(preconditioner_env%sparse_matrix)) THEN
960 4441 : ALLOCATE (preconditioner_env%sparse_matrix)
961 : END IF
962 :
963 : ! Put the first term of the preconditioner (H) into the output matrix
964 4441 : CALL dbcsr_copy(preconditioner_env%sparse_matrix, matrix_h)
965 :
966 : ! Precompute some matrices
967 : ! S*C, if orthonormal this will be simply C so a copy will do
968 4441 : IF (PRESENT(matrix_s)) THEN
969 4053 : CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_s, mo_dbcsr, 0.0_dp, dbcsr_sc)
970 : ELSE
971 388 : CALL dbcsr_copy(dbcsr_sc, mo_dbcsr)
972 : END IF
973 :
974 : !----------------------------compute the occupied subspace and shift it ------------------------------------
975 : ! cT*H*C which will be used to shift the occupied states to 0
976 4441 : CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_h, mo_dbcsr, 0.0_dp, dbcsr_hc)
977 4441 : CALL dbcsr_multiply("T", "N", 1.0_dp, mo_dbcsr, dbcsr_hc, 0.0_dp, dbcsr_cThc)
978 :
979 : ! Compute the Energy of the HOMO. We will use this as a reference energy
980 8882 : ALLOCATE (matrices(1))
981 4441 : matrices(1)%matrix => dbcsr_cThc
982 : CALL setup_arnoldi_env(arnoldi_env, matrices, max_iter=20, threshold=1.0E-3_dp, selection_crit=2, &
983 4441 : nval_request=1, nrestarts=8, generalized_ev=.FALSE., iram=.FALSE.)
984 4441 : IF (ASSOCIATED(preconditioner_env%max_ev_vector)) THEN
985 2440 : CALL set_arnoldi_initial_vector(arnoldi_env, preconditioner_env%max_ev_vector)
986 : END IF
987 4441 : CALL arnoldi_ev(matrices, arnoldi_env)
988 4441 : max_ev = REAL(get_selected_ritz_val(arnoldi_env, 1), dp)
989 :
990 : ! save the ev as guess for the next time
991 4441 : IF (.NOT. ASSOCIATED(preconditioner_env%max_ev_vector)) ALLOCATE (preconditioner_env%max_ev_vector)
992 4441 : CALL get_selected_ritz_vector(arnoldi_env, 1, matrices(1)%matrix, preconditioner_env%max_ev_vector)
993 4441 : CALL deallocate_arnoldi_env(arnoldi_env)
994 4441 : DEALLOCATE (matrices)
995 :
996 : ! Lets shift the occupied states a bit further up, -1.0 because we gonna subtract it from H
997 4441 : CALL dbcsr_add_on_diag(dbcsr_cThc, -0.5_dp)
998 : ! Get the AO representation of the shift (see above why S is needed), W-matrix like object
999 4441 : CALL dbcsr_multiply("N", "N", 2.0_dp, dbcsr_sc, dbcsr_cThc, 0.0_dp, dbcsr_hc)
1000 4441 : CALL dbcsr_multiply("N", "T", -1.0_dp, dbcsr_hc, dbcsr_sc, 1.0_dp, preconditioner_env%sparse_matrix)
1001 :
1002 : !-------------------------------------compute eigenvalues of H ----------------------------------------------
1003 : ! Setup the arnoldi procedure to compute the lowest ev. if S is present this has to be the generalized ev
1004 4441 : IF (PRESENT(matrix_s)) THEN
1005 12159 : ALLOCATE (matrices(2))
1006 4053 : matrices(1)%matrix => preconditioner_env%sparse_matrix
1007 4053 : matrices(2)%matrix => matrix_s
1008 : CALL setup_arnoldi_env(arnoldi_env, matrices, max_iter=20, threshold=2.0E-2_dp, selection_crit=3, &
1009 4053 : nval_request=1, nrestarts=21, generalized_ev=.TRUE., iram=.FALSE.)
1010 : ELSE
1011 776 : ALLOCATE (matrices(1))
1012 388 : matrices(1)%matrix => preconditioner_env%sparse_matrix
1013 : CALL setup_arnoldi_env(arnoldi_env, matrices, max_iter=20, threshold=2.0E-2_dp, selection_crit=3, &
1014 388 : nval_request=1, nrestarts=8, generalized_ev=.FALSE., iram=.FALSE.)
1015 : END IF
1016 4441 : IF (ASSOCIATED(preconditioner_env%min_ev_vector)) THEN
1017 2440 : CALL set_arnoldi_initial_vector(arnoldi_env, preconditioner_env%min_ev_vector)
1018 : END IF
1019 :
1020 : ! compute the LUMO energy
1021 4441 : CALL arnoldi_ev(matrices, arnoldi_env)
1022 4441 : min_eV = REAL(get_selected_ritz_val(arnoldi_env, 1), dp)
1023 :
1024 : ! save the lumo vector for restarting in the next step
1025 4441 : IF (.NOT. ASSOCIATED(preconditioner_env%min_ev_vector)) ALLOCATE (preconditioner_env%min_ev_vector)
1026 4441 : CALL get_selected_ritz_vector(arnoldi_env, 1, matrices(1)%matrix, preconditioner_env%min_ev_vector)
1027 4441 : CALL deallocate_arnoldi_env(arnoldi_env)
1028 4441 : DEALLOCATE (matrices)
1029 :
1030 : !-------------------------------------compute eigenvalues of H ----------------------------------------------
1031 : ! Shift the Lumo to the 1.5*the computed energy_gap or the external energy gap value
1032 : ! The factor 1.5 is determined by trying. If the LUMO is positive, enough, just leave it alone
1033 4441 : pre_shift = MAX(1.5_dp*(min_ev - max_ev), energy_gap)
1034 4441 : IF (min_ev < pre_shift) THEN
1035 4423 : pre_shift = pre_shift - min_ev
1036 : ELSE
1037 18 : pre_shift = 0.0_dp
1038 : END IF
1039 4441 : IF (PRESENT(matrix_s)) THEN
1040 4053 : CALL dbcsr_add(preconditioner_env%sparse_matrix, matrix_s, 1.0_dp, pre_shift)
1041 : ELSE
1042 388 : CALL dbcsr_add_on_diag(preconditioner_env%sparse_matrix, pre_shift)
1043 : END IF
1044 :
1045 4441 : CALL dbcsr_release(mo_dbcsr)
1046 4441 : CALL dbcsr_release(dbcsr_hc)
1047 4441 : CALL dbcsr_release(dbcsr_sc)
1048 4441 : CALL dbcsr_release(dbcsr_cThc)
1049 :
1050 4441 : CALL timestop(handle)
1051 :
1052 4441 : END SUBROUTINE make_full_single_inverse
1053 :
1054 : ! **************************************************************************************************
1055 : !> \brief Build a gauge-covariant FULL_SINGLE_INVERSE operator for a complex k-point channel.
1056 : !> \param preconditioner_env preconditioner storage
1057 : !> \param matrix_c0 complex reference orbitals, C^H S C = I
1058 : !> \param matrix_h complex Hermitian k-point Hamiltonian
1059 : !> \param matrix_s complex Hermitian k-point overlap
1060 : !> \param energy_gap lower spectral bound of the positive operator
1061 : ! **************************************************************************************************
1062 64 : SUBROUTINE make_complex_full_single_inverse(preconditioner_env, matrix_c0, matrix_h, matrix_s, &
1063 : energy_gap)
1064 :
1065 : TYPE(preconditioner_type) :: preconditioner_env
1066 : TYPE(cp_cfm_type), INTENT(IN) :: matrix_c0, matrix_h, matrix_s
1067 : REAL(KIND=dp), INTENT(IN) :: energy_gap
1068 :
1069 : CHARACTER(len=*), PARAMETER :: routineN = 'make_complex_full_single_inverse'
1070 :
1071 : INTEGER :: handle, i, k, n
1072 : REAL(KIND=dp) :: max_ev, min_ev, pre_shift, target_edge
1073 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: occupied_evals, operator_evals
1074 : TYPE(cp_cfm_type) :: matrix_chc, matrix_diag_a, matrix_diag_b, matrix_diag_evec, &
1075 : matrix_diag_work, matrix_hc0, matrix_occ_diag_a, matrix_occ_diag_b, matrix_occ_diag_evec, &
1076 : matrix_occ_diag_work, matrix_operator, matrix_sc0, matrix_sc_chc
1077 :
1078 64 : CALL timeset(routineN, handle)
1079 :
1080 64 : CALL cp_cfm_get_info(matrix_c0, nrow_global=n, ncol_global=k)
1081 64 : CPASSERT(n > 0)
1082 64 : CPASSERT(k > 0 .AND. k <= n)
1083 64 : CPASSERT(energy_gap > 0.0_dp)
1084 :
1085 64 : CALL cp_cfm_create(matrix_hc0, matrix_c0%matrix_struct, name='complex FULL_SINGLE HC')
1086 64 : CALL cp_cfm_create(matrix_sc0, matrix_c0%matrix_struct, name='complex FULL_SINGLE SC')
1087 : CALL cp_cfm_create(matrix_chc, matrix_c0%matrix_struct, nrow=k, ncol=k, &
1088 64 : name='complex FULL_SINGLE CHC')
1089 : CALL cp_cfm_create(matrix_sc_chc, matrix_c0%matrix_struct, &
1090 64 : name='complex FULL_SINGLE SC CHC')
1091 : CALL cp_cfm_create(matrix_operator, matrix_h%matrix_struct, &
1092 64 : name='complex FULL_SINGLE operator')
1093 : CALL cp_cfm_create(matrix_diag_a, matrix_h%matrix_struct, &
1094 64 : name='complex FULL_SINGLE diagonalization A')
1095 : CALL cp_cfm_create(matrix_diag_b, matrix_s%matrix_struct, &
1096 64 : name='complex FULL_SINGLE diagonalization B')
1097 : CALL cp_cfm_create(matrix_diag_evec, matrix_h%matrix_struct, nrow=n, ncol=1, &
1098 64 : name='complex FULL_SINGLE eigenvectors')
1099 : CALL cp_cfm_create(matrix_diag_work, matrix_h%matrix_struct, &
1100 64 : name='complex FULL_SINGLE diagonalization work')
1101 : CALL cp_cfm_create(matrix_occ_diag_a, matrix_c0%matrix_struct, nrow=k, ncol=k, &
1102 64 : name='complex FULL_SINGLE occupied diagonalization A')
1103 : CALL cp_cfm_create(matrix_occ_diag_b, matrix_c0%matrix_struct, nrow=k, ncol=k, &
1104 64 : name='complex FULL_SINGLE occupied diagonalization B', set_zero=.TRUE.)
1105 : CALL cp_cfm_create(matrix_occ_diag_evec, matrix_c0%matrix_struct, nrow=k, ncol=k, &
1106 64 : name='complex FULL_SINGLE occupied eigenvectors')
1107 : CALL cp_cfm_create(matrix_occ_diag_work, matrix_c0%matrix_struct, nrow=k, ncol=k, &
1108 64 : name='complex FULL_SINGLE occupied diagonalization work')
1109 :
1110 64 : CALL cp_cfm_gemm('N', 'N', n, k, n, z_one, matrix_h, matrix_c0, z_zero, matrix_hc0)
1111 64 : CALL cp_cfm_gemm('N', 'N', n, k, n, z_one, matrix_s, matrix_c0, z_zero, matrix_sc0)
1112 64 : CALL cp_cfm_gemm('C', 'N', k, k, n, z_one, matrix_c0, matrix_hc0, z_zero, matrix_chc)
1113 :
1114 : ! A = H - 2*S*C*(C^H*H*C)*C^H*S + S*C*C^H*S is invariant under C -> C*U.
1115 64 : CALL cp_cfm_to_cfm(matrix_h, matrix_operator)
1116 : CALL cp_cfm_gemm('N', 'N', n, k, k, z_one, matrix_sc0, matrix_chc, &
1117 64 : z_zero, matrix_sc_chc)
1118 : CALL cp_cfm_gemm('N', 'C', n, n, k, -2.0_dp*z_one, matrix_sc_chc, matrix_sc0, &
1119 64 : z_one, matrix_operator)
1120 : CALL cp_cfm_gemm('N', 'C', n, n, k, z_one, matrix_sc0, matrix_sc0, &
1121 64 : z_one, matrix_operator)
1122 :
1123 192 : ALLOCATE (occupied_evals(k), operator_evals(1))
1124 64 : CALL cp_cfm_to_cfm(matrix_chc, matrix_occ_diag_a)
1125 484 : DO i = 1, k
1126 484 : CALL cp_cfm_set_element(matrix_occ_diag_b, i, i, z_one)
1127 : END DO
1128 : CALL cp_cfm_geeig(matrix_occ_diag_a, matrix_occ_diag_b, matrix_occ_diag_evec, &
1129 64 : occupied_evals, matrix_occ_diag_work)
1130 484 : max_ev = MAXVAL(occupied_evals)
1131 64 : CALL cp_cfm_to_cfm(matrix_operator, matrix_diag_a)
1132 64 : CALL cp_cfm_to_cfm(matrix_s, matrix_diag_b)
1133 : CALL cp_cfm_geeig(matrix_diag_a, matrix_diag_b, matrix_diag_evec, &
1134 64 : operator_evals, matrix_diag_work, lowest_subset=.TRUE.)
1135 64 : min_ev = operator_evals(1)
1136 64 : target_edge = MAX(1.5_dp*(min_ev - max_ev), energy_gap)
1137 64 : pre_shift = MAX(0.0_dp, target_edge - min_ev)
1138 64 : IF (pre_shift > 0.0_dp) THEN
1139 : CALL cp_cfm_scale_and_add(z_one, matrix_operator, &
1140 64 : CMPLX(pre_shift, 0.0_dp, KIND=dp), matrix_s)
1141 : END IF
1142 :
1143 : CALL store_complex_inverse(preconditioner_env, matrix_operator, &
1144 64 : ot_precond_full_single_inverse, energy_gap)
1145 :
1146 64 : DEALLOCATE (operator_evals, occupied_evals)
1147 64 : CALL cp_cfm_release(matrix_occ_diag_work)
1148 64 : CALL cp_cfm_release(matrix_occ_diag_evec)
1149 64 : CALL cp_cfm_release(matrix_occ_diag_b)
1150 64 : CALL cp_cfm_release(matrix_occ_diag_a)
1151 64 : CALL cp_cfm_release(matrix_diag_work)
1152 64 : CALL cp_cfm_release(matrix_diag_evec)
1153 64 : CALL cp_cfm_release(matrix_diag_b)
1154 64 : CALL cp_cfm_release(matrix_diag_a)
1155 64 : CALL cp_cfm_release(matrix_operator)
1156 64 : CALL cp_cfm_release(matrix_sc_chc)
1157 64 : CALL cp_cfm_release(matrix_chc)
1158 64 : CALL cp_cfm_release(matrix_sc0)
1159 64 : CALL cp_cfm_release(matrix_hc0)
1160 :
1161 64 : CALL timestop(handle)
1162 :
1163 256 : END SUBROUTINE make_complex_full_single_inverse
1164 :
1165 : ! **************************************************************************************************
1166 : !> \brief Build the inverse complex overlap preconditioner.
1167 : !> \param preconditioner_env preconditioner storage
1168 : !> \param matrix_s complex Hermitian k-point overlap
1169 : ! **************************************************************************************************
1170 92 : SUBROUTINE make_complex_full_s_inverse(preconditioner_env, matrix_s)
1171 :
1172 : TYPE(preconditioner_type) :: preconditioner_env
1173 : TYPE(cp_cfm_type), INTENT(IN) :: matrix_s
1174 :
1175 92 : CALL store_complex_inverse(preconditioner_env, matrix_s, ot_precond_s_inverse, 0.0_dp)
1176 :
1177 92 : END SUBROUTINE make_complex_full_s_inverse
1178 :
1179 : ! **************************************************************************************************
1180 : !> \brief Build the inverse complex kinetic-plus-overlap preconditioner.
1181 : !> \param preconditioner_env preconditioner storage
1182 : !> \param matrix_t complex Hermitian k-point kinetic operator
1183 : !> \param matrix_s complex Hermitian k-point overlap
1184 : !> \param energy_gap non-negative overlap shift
1185 : ! **************************************************************************************************
1186 122 : SUBROUTINE make_complex_full_kinetic(preconditioner_env, matrix_t, matrix_s, energy_gap)
1187 :
1188 : TYPE(preconditioner_type) :: preconditioner_env
1189 : TYPE(cp_cfm_type), INTENT(IN) :: matrix_t, matrix_s
1190 : REAL(KIND=dp), INTENT(IN) :: energy_gap
1191 :
1192 : REAL(KIND=dp) :: shift
1193 : TYPE(cp_cfm_type) :: matrix_operator
1194 :
1195 122 : shift = MAX(0.0_dp, energy_gap)
1196 : CALL cp_cfm_create(matrix_operator, matrix_t%matrix_struct, &
1197 122 : name='complex FULL_KINETIC operator')
1198 122 : CALL cp_cfm_to_cfm(matrix_t, matrix_operator)
1199 122 : CALL cp_cfm_scale_and_add(z_one, matrix_operator, CMPLX(shift, 0.0_dp, KIND=dp), matrix_s)
1200 : CALL store_complex_inverse(preconditioner_env, matrix_operator, &
1201 122 : ot_precond_full_kinetic, energy_gap)
1202 122 : CALL cp_cfm_release(matrix_operator)
1203 :
1204 122 : END SUBROUTINE make_complex_full_kinetic
1205 :
1206 : ! **************************************************************************************************
1207 : !> \brief Store an explicitly inverted positive complex Hermitian operator.
1208 : !> \param preconditioner_env preconditioner storage
1209 : !> \param matrix_operator positive complex Hermitian operator
1210 : !> \param preconditioner_kind selected OT preconditioner
1211 : !> \param energy_gap configured spectral shift
1212 : ! **************************************************************************************************
1213 834 : SUBROUTINE store_complex_inverse(preconditioner_env, matrix_operator, &
1214 : preconditioner_kind, energy_gap)
1215 :
1216 : TYPE(preconditioner_type) :: preconditioner_env
1217 : TYPE(cp_cfm_type), INTENT(IN) :: matrix_operator
1218 : INTEGER, INTENT(IN) :: preconditioner_kind
1219 : REAL(KIND=dp), INTENT(IN) :: energy_gap
1220 :
1221 : INTEGER :: info
1222 :
1223 278 : IF (ASSOCIATED(preconditioner_env%complex_fm)) THEN
1224 0 : CALL cp_cfm_release(preconditioner_env%complex_fm)
1225 0 : DEALLOCATE (preconditioner_env%complex_fm)
1226 : END IF
1227 278 : ALLOCATE (preconditioner_env%complex_fm)
1228 : CALL cp_cfm_create(preconditioner_env%complex_fm, matrix_operator%matrix_struct, &
1229 278 : name='complex inverse preconditioner')
1230 278 : CALL cp_cfm_to_cfm(matrix_operator, preconditioner_env%complex_fm)
1231 278 : CALL cp_cfm_cholesky_decompose(preconditioner_env%complex_fm, info_out=info)
1232 278 : CPASSERT(info == 0)
1233 278 : CALL cp_cfm_cholesky_invert(preconditioner_env%complex_fm, info_out=info)
1234 278 : CPASSERT(info == 0)
1235 278 : CALL cp_cfm_uplo_to_full(preconditioner_env%complex_fm)
1236 :
1237 278 : preconditioner_env%energy_gap = energy_gap
1238 278 : preconditioner_env%in_use = preconditioner_kind
1239 278 : preconditioner_env%solver = ot_precond_solver_default
1240 :
1241 278 : END SUBROUTINE store_complex_inverse
1242 :
1243 : END MODULE preconditioner_makes
|
