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_fermi_low_rank, ot_precond_full_all, &
64 : ot_precond_full_kinetic, ot_precond_full_single, ot_precond_full_single_inverse, &
65 : ot_precond_s_inverse, 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 1610 : 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 1610 : DO j = 1, ncol_local
173 : residual_norm_sq(col_indices(j)) = residual_norm_sq(col_indices(j)) + &
174 62054 : SUM(ABS(local_data(1:nrow_local, j))**2)
175 : END DO
176 152 : CALL preconditioner_env%para_env%sum(residual_norm_sq)
177 1610 : 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 1610 : 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 1610 : preconditioner_env%occ_evals(:) = c0_evals(1:k)
201 1610 : 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 9704 : SUBROUTINE make_preconditioner_matrix(preconditioner_env, matrix_h, matrix_s, matrix_t, mo_coeff, &
313 9704 : 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 9704 : 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 4501 : 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 4501 : matrix_s=matrix_s)
357 : CASE (ot_precond_full_all)
358 3762 : IF (my_solver_type /= ot_precond_solver_default) THEN
359 0 : CPABORT("Only PRECOND_SOLVER DEFAULT for the moment")
360 : END IF
361 3762 : IF (PRESENT(matrix_s)) THEN
362 : CALL make_full_all(preconditioner_env, mo_coeff, matrix_h, matrix_s, &
363 3682 : 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 (ot_precond_fermi_low_rank)
370 4 : IF (my_solver_type /= ot_precond_solver_default) THEN
371 0 : CPABORT("Only PRECOND_SOLVER DEFAULT for the moment")
372 : END IF
373 4 : IF (PRESENT(matrix_s)) THEN
374 : CALL make_full_all(preconditioner_env, mo_coeff, matrix_h, matrix_s, &
375 4 : eigenvalues_ot, energy_gap, common_reference=energy_homo)
376 : ELSE
377 : CALL make_full_all_ortho(preconditioner_env, mo_coeff, matrix_h, &
378 0 : eigenvalues_ot, energy_gap, common_reference=energy_homo)
379 : END IF
380 :
381 : CASE DEFAULT
382 9704 : CPABORT("Type not implemented")
383 : END SELECT
384 :
385 9704 : END SUBROUTINE make_preconditioner_matrix
386 :
387 : ! **************************************************************************************************
388 : !> \brief Simply takes the overlap matrix as preconditioner
389 : !> \param preconditioner_env ...
390 : !> \param matrix_s ...
391 : ! **************************************************************************************************
392 78 : SUBROUTINE make_full_s_inverse(preconditioner_env, matrix_s)
393 : TYPE(preconditioner_type) :: preconditioner_env
394 : TYPE(dbcsr_type), POINTER :: matrix_s
395 :
396 : CHARACTER(len=*), PARAMETER :: routineN = 'make_full_s_inverse'
397 :
398 : INTEGER :: handle
399 :
400 78 : CALL timeset(routineN, handle)
401 :
402 78 : CPASSERT(ASSOCIATED(matrix_s))
403 :
404 78 : IF (.NOT. ASSOCIATED(preconditioner_env%sparse_matrix)) THEN
405 78 : ALLOCATE (preconditioner_env%sparse_matrix)
406 : END IF
407 78 : CALL dbcsr_copy(preconditioner_env%sparse_matrix, matrix_s, name="full_kinetic")
408 :
409 78 : CALL timestop(handle)
410 :
411 78 : END SUBROUTINE make_full_s_inverse
412 :
413 : ! **************************************************************************************************
414 : !> \brief kinetic matrix+shift*overlap as preconditioner. Cheap but could
415 : !> be better
416 : !> \param preconditioner_env ...
417 : !> \param matrix_t ...
418 : !> \param matrix_s ...
419 : !> \param energy_gap ...
420 : ! **************************************************************************************************
421 1321 : SUBROUTINE make_full_kinetic(preconditioner_env, matrix_t, matrix_s, &
422 : energy_gap)
423 : TYPE(preconditioner_type) :: preconditioner_env
424 : TYPE(dbcsr_type), POINTER :: matrix_t, matrix_s
425 : REAL(KIND=dp) :: energy_gap
426 :
427 : CHARACTER(len=*), PARAMETER :: routineN = 'make_full_kinetic'
428 :
429 : INTEGER :: handle
430 : REAL(KIND=dp) :: shift
431 :
432 1321 : CALL timeset(routineN, handle)
433 :
434 1321 : CPASSERT(ASSOCIATED(matrix_t))
435 1321 : CPASSERT(ASSOCIATED(matrix_s))
436 :
437 1321 : IF (.NOT. ASSOCIATED(preconditioner_env%sparse_matrix)) THEN
438 1319 : ALLOCATE (preconditioner_env%sparse_matrix)
439 : END IF
440 1321 : CALL dbcsr_copy(preconditioner_env%sparse_matrix, matrix_t, name="full_kinetic")
441 :
442 1321 : shift = MAX(0.0_dp, energy_gap)
443 :
444 : CALL dbcsr_add(preconditioner_env%sparse_matrix, matrix_s, &
445 1321 : alpha_scalar=1.0_dp, beta_scalar=shift)
446 :
447 1321 : CALL timestop(handle)
448 :
449 1321 : END SUBROUTINE make_full_kinetic
450 :
451 : ! **************************************************************************************************
452 : !> \brief full_single_preconditioner
453 : !> \param preconditioner_env ...
454 : !> \param fm ...
455 : !> \param matrix_h ...
456 : !> \param matrix_s ...
457 : !> \param energy_homo ...
458 : !> \param energy_gap ...
459 : ! **************************************************************************************************
460 32 : SUBROUTINE make_full_single(preconditioner_env, fm, matrix_h, matrix_s, &
461 : energy_homo, energy_gap)
462 : TYPE(preconditioner_type) :: preconditioner_env
463 : TYPE(cp_fm_type), POINTER :: fm
464 : TYPE(dbcsr_type), POINTER :: matrix_h, matrix_s
465 : REAL(KIND=dp) :: energy_homo, energy_gap
466 :
467 : CHARACTER(len=*), PARAMETER :: routineN = 'make_full_single'
468 :
469 : INTEGER :: handle, i, n
470 32 : REAL(KIND=dp), DIMENSION(:), POINTER :: evals
471 : TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp
472 : TYPE(cp_fm_type) :: fm_h, fm_s
473 :
474 32 : CALL timeset(routineN, handle)
475 :
476 32 : NULLIFY (fm_struct_tmp, evals)
477 :
478 32 : IF (ASSOCIATED(fm)) THEN
479 0 : CALL cp_fm_release(fm)
480 0 : DEALLOCATE (fm)
481 : NULLIFY (fm)
482 : END IF
483 32 : CALL dbcsr_get_info(matrix_h, nfullrows_total=n)
484 96 : ALLOCATE (evals(n))
485 :
486 : CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=n, ncol_global=n, &
487 : context=preconditioner_env%ctxt, &
488 32 : para_env=preconditioner_env%para_env)
489 32 : ALLOCATE (fm)
490 32 : CALL cp_fm_create(fm, fm_struct_tmp, name="preconditioner")
491 32 : CALL cp_fm_create(fm_h, fm_struct_tmp, name="fm_h")
492 32 : CALL cp_fm_create(fm_s, fm_struct_tmp, name="fm_s")
493 32 : CALL cp_fm_struct_release(fm_struct_tmp)
494 :
495 32 : CALL copy_dbcsr_to_fm(matrix_h, fm_h)
496 32 : CALL copy_dbcsr_to_fm(matrix_s, fm_s)
497 32 : CALL cp_fm_cholesky_decompose(fm_s)
498 :
499 32 : SELECT CASE (preconditioner_env%cholesky_use)
500 : CASE (cholesky_inverse)
501 : ! if cho inverse
502 0 : CALL cp_fm_triangular_invert(fm_s)
503 0 : CALL cp_fm_uplo_to_full(fm_h, fm)
504 :
505 : CALL cp_fm_triangular_multiply(fm_s, fm_h, side="R", transpose_tr=.FALSE., &
506 0 : invert_tr=.FALSE., uplo_tr="U", n_rows=n, n_cols=n, alpha=1.0_dp)
507 : CALL cp_fm_triangular_multiply(fm_s, fm_h, side="L", transpose_tr=.TRUE., &
508 0 : invert_tr=.FALSE., uplo_tr="U", n_rows=n, n_cols=n, alpha=1.0_dp)
509 : CASE (cholesky_reduce)
510 32 : CALL cp_fm_cholesky_reduce(fm_h, fm_s)
511 : CASE DEFAULT
512 32 : CPABORT("cholesky type not implemented")
513 : END SELECT
514 :
515 32 : CALL choose_eigv_solver(fm_h, fm, evals)
516 :
517 32 : SELECT CASE (preconditioner_env%cholesky_use)
518 : CASE (cholesky_inverse)
519 : CALL cp_fm_triangular_multiply(fm_s, fm, side="L", transpose_tr=.FALSE., &
520 0 : invert_tr=.FALSE., uplo_tr="U", n_rows=n, n_cols=n, alpha=1.0_dp)
521 0 : DO i = 1, n
522 0 : evals(i) = 1.0_dp/MAX(evals(i) - energy_homo, energy_gap)
523 : END DO
524 0 : CALL cp_fm_to_fm(fm, fm_h)
525 : CASE (cholesky_reduce)
526 32 : CALL cp_fm_cholesky_restore(fm, n, fm_s, fm_h, "SOLVE")
527 568 : DO i = 1, n
528 568 : evals(i) = 1.0_dp/MAX(evals(i) - energy_homo, energy_gap)
529 : END DO
530 64 : CALL cp_fm_to_fm(fm_h, fm)
531 : END SELECT
532 :
533 32 : CALL cp_fm_column_scale(fm, evals)
534 32 : CALL parallel_gemm('N', 'T', n, n, n, 1.0_dp, fm, fm_h, 0.0_dp, fm_s)
535 32 : CALL cp_fm_to_fm(fm_s, fm)
536 :
537 32 : DEALLOCATE (evals)
538 32 : CALL cp_fm_release(fm_h)
539 32 : CALL cp_fm_release(fm_s)
540 :
541 32 : CALL timestop(handle)
542 :
543 96 : END SUBROUTINE make_full_single
544 :
545 : ! **************************************************************************************************
546 : !> \brief full single in the orthonormal basis
547 : !> \param preconditioner_env ...
548 : !> \param fm ...
549 : !> \param matrix_h ...
550 : !> \param energy_homo ...
551 : !> \param energy_gap ...
552 : ! **************************************************************************************************
553 6 : SUBROUTINE make_full_single_ortho(preconditioner_env, fm, matrix_h, &
554 : energy_homo, energy_gap)
555 : TYPE(preconditioner_type) :: preconditioner_env
556 : TYPE(cp_fm_type), POINTER :: fm
557 : TYPE(dbcsr_type), POINTER :: matrix_h
558 : REAL(KIND=dp) :: energy_homo, energy_gap
559 :
560 : CHARACTER(len=*), PARAMETER :: routineN = 'make_full_single_ortho'
561 :
562 : INTEGER :: handle, i, n
563 6 : REAL(KIND=dp), DIMENSION(:), POINTER :: evals
564 : TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp
565 : TYPE(cp_fm_type) :: fm_h, fm_s
566 :
567 6 : CALL timeset(routineN, handle)
568 6 : NULLIFY (fm_struct_tmp, evals)
569 :
570 6 : IF (ASSOCIATED(fm)) THEN
571 0 : CALL cp_fm_release(fm)
572 0 : DEALLOCATE (fm)
573 : NULLIFY (fm)
574 : END IF
575 6 : CALL dbcsr_get_info(matrix_h, nfullrows_total=n)
576 18 : ALLOCATE (evals(n))
577 :
578 : CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=n, ncol_global=n, &
579 : context=preconditioner_env%ctxt, &
580 6 : para_env=preconditioner_env%para_env)
581 6 : ALLOCATE (fm)
582 6 : CALL cp_fm_create(fm, fm_struct_tmp, name="preconditioner")
583 6 : CALL cp_fm_create(fm_h, fm_struct_tmp, name="fm_h")
584 6 : CALL cp_fm_create(fm_s, fm_struct_tmp, name="fm_s")
585 6 : CALL cp_fm_struct_release(fm_struct_tmp)
586 :
587 6 : CALL copy_dbcsr_to_fm(matrix_h, fm_h)
588 :
589 6 : CALL choose_eigv_solver(fm_h, fm, evals)
590 282 : DO i = 1, n
591 282 : evals(i) = 1.0_dp/MAX(evals(i) - energy_homo, energy_gap)
592 : END DO
593 6 : CALL cp_fm_to_fm(fm, fm_h)
594 6 : CALL cp_fm_column_scale(fm, evals)
595 6 : CALL parallel_gemm('N', 'T', n, n, n, 1.0_dp, fm, fm_h, 0.0_dp, fm_s)
596 6 : CALL cp_fm_to_fm(fm_s, fm)
597 :
598 6 : DEALLOCATE (evals)
599 6 : CALL cp_fm_release(fm_h)
600 6 : CALL cp_fm_release(fm_s)
601 :
602 6 : CALL timestop(handle)
603 :
604 18 : END SUBROUTINE make_full_single_ortho
605 :
606 : ! **************************************************************************************************
607 : !> \brief generates a state by state preconditioner based on the full hamiltonian matrix
608 : !> \param preconditioner_env ...
609 : !> \param matrix_c0 ...
610 : !> \param matrix_h ...
611 : !> \param matrix_s ...
612 : !> \param c0_evals ...
613 : !> \param energy_gap should be a slight underestimate of the physical energy gap for almost all systems
614 : !> the c0 are already ritz states of (h,s)
615 : !> \param common_reference optional scalar occupied-space level used by FERMI_LOW_RANK
616 : !> \par History
617 : !> 10.2006 made more stable [Joost VandeVondele]
618 : !> \note
619 : !> includes error estimate on the hamiltonian matrix to result in a stable preconditioner
620 : !> a preconditioner for each eigenstate i is generated by keeping the factorized form
621 : !> U diag( something i ) U^T. It is important to only precondition in the subspace orthogonal to c0.
622 : !> not only is it the only part that matters, it also simplifies the computation of
623 : !> the lagrangian multipliers in the OT minimization (i.e. if the c0 here is different
624 : !> from the c0 used in the OT setup, there will be a bug).
625 : ! **************************************************************************************************
626 3686 : SUBROUTINE make_full_all(preconditioner_env, matrix_c0, matrix_h, matrix_s, c0_evals, energy_gap, &
627 : common_reference)
628 : TYPE(preconditioner_type) :: preconditioner_env
629 : TYPE(cp_fm_type), INTENT(IN) :: matrix_c0
630 : TYPE(dbcsr_type), POINTER :: matrix_h, matrix_s
631 : REAL(KIND=dp), DIMENSION(:) :: c0_evals
632 : REAL(KIND=dp) :: energy_gap
633 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: common_reference
634 :
635 : CHARACTER(len=*), PARAMETER :: routineN = 'make_full_all'
636 : REAL(KIND=dp), PARAMETER :: fudge_factor = 0.25_dp, &
637 : lambda_base = 10.0_dp
638 :
639 : INTEGER :: handle, k, n
640 : REAL(KIND=dp) :: error_estimate, lambda
641 3686 : REAL(KIND=dp), DIMENSION(:), POINTER :: diag, norms, shifted_evals
642 : TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp
643 : TYPE(cp_fm_type) :: matrix_hc0, matrix_left, matrix_s1, &
644 : matrix_s2, matrix_sc0, matrix_shc0, &
645 : matrix_tmp, ortho
646 : TYPE(cp_fm_type), POINTER :: matrix_pre
647 :
648 3686 : CALL timeset(routineN, handle)
649 :
650 3686 : IF (ASSOCIATED(preconditioner_env%fm)) THEN
651 0 : CALL cp_fm_release(preconditioner_env%fm)
652 0 : DEALLOCATE (preconditioner_env%fm)
653 : NULLIFY (preconditioner_env%fm)
654 : END IF
655 3686 : CALL cp_fm_get_info(matrix_c0, nrow_global=n, ncol_global=k)
656 : CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=n, ncol_global=n, &
657 : context=preconditioner_env%ctxt, &
658 3686 : para_env=preconditioner_env%para_env)
659 3686 : ALLOCATE (preconditioner_env%fm)
660 3686 : CALL cp_fm_create(preconditioner_env%fm, fm_struct_tmp, name="preconditioner_env%fm")
661 3686 : CALL cp_fm_create(ortho, fm_struct_tmp, name="ortho")
662 3686 : CALL cp_fm_create(matrix_tmp, fm_struct_tmp, name="matrix_tmp")
663 3686 : CALL cp_fm_struct_release(fm_struct_tmp)
664 11058 : ALLOCATE (preconditioner_env%full_evals(n))
665 10948 : ALLOCATE (preconditioner_env%occ_evals(k))
666 :
667 : ! 0) cholesky decompose the overlap matrix, if this fails the basis is singular,
668 : ! more than EPS_DEFAULT
669 3686 : CALL copy_dbcsr_to_fm(matrix_s, ortho)
670 3686 : CALL cp_fm_cholesky_decompose(ortho)
671 : ! if cho inverse
672 3686 : IF (preconditioner_env%cholesky_use == cholesky_inverse) THEN
673 0 : CALL cp_fm_triangular_invert(ortho)
674 : END IF
675 : ! 1) Construct a new H matrix, which has the current C0 as eigenvectors,
676 : ! possibly shifted by an amount lambda,
677 : ! and the same spectrum as the original H matrix in the space orthogonal to the C0
678 : ! with P=C0 C0 ^ T
679 : ! (1 - PS)^T H (1-PS) + (PS)^T (H - lambda S ) (PS)
680 : ! we exploit that the C0 are already the ritz states of H
681 3686 : CALL cp_fm_create(matrix_sc0, matrix_c0%matrix_struct, name="sc0")
682 3686 : CALL cp_dbcsr_sm_fm_multiply(matrix_s, matrix_c0, matrix_sc0, k)
683 3686 : CALL cp_fm_create(matrix_hc0, matrix_c0%matrix_struct, name="hc0")
684 3686 : CALL cp_dbcsr_sm_fm_multiply(matrix_h, matrix_c0, matrix_hc0, k)
685 :
686 : ! An aside, try to estimate the error on the ritz values, we'll need it later on
687 3686 : CALL cp_fm_create(matrix_shc0, matrix_c0%matrix_struct, name="shc0")
688 :
689 3686 : SELECT CASE (preconditioner_env%cholesky_use)
690 : CASE (cholesky_inverse)
691 : ! if cho inverse
692 0 : CALL cp_fm_to_fm(matrix_hc0, matrix_shc0)
693 : CALL cp_fm_triangular_multiply(ortho, matrix_shc0, side="L", transpose_tr=.TRUE., &
694 0 : invert_tr=.FALSE., uplo_tr="U", n_rows=n, n_cols=k, alpha=1.0_dp)
695 : CASE (cholesky_reduce)
696 3686 : CALL cp_fm_cholesky_restore(matrix_hc0, k, ortho, matrix_shc0, "SOLVE", transa="T")
697 : CASE DEFAULT
698 3686 : CPABORT("cholesky type not implemented")
699 : END SELECT
700 : CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=k, ncol_global=k, &
701 : context=preconditioner_env%ctxt, &
702 3686 : para_env=preconditioner_env%para_env)
703 3686 : CALL cp_fm_create(matrix_s1, fm_struct_tmp, name="matrix_s1")
704 3686 : CALL cp_fm_struct_release(fm_struct_tmp)
705 : ! since we only use diagonal elements this is a bit of a waste
706 3686 : CALL parallel_gemm('T', 'N', k, k, n, 1.0_dp, matrix_shc0, matrix_shc0, 0.0_dp, matrix_s1)
707 7262 : ALLOCATE (diag(k))
708 3686 : CALL cp_fm_get_diag(matrix_s1, diag)
709 3686 : IF (PRESENT(common_reference)) THEN
710 : ! ||(1-CC^T S)S^-1 H C||_F is invariant under C -> C U.
711 68 : error_estimate = SQRT(MAX(0.0_dp, SUM(diag) - SUM(c0_evals**2)))
712 : ELSE
713 21142 : error_estimate = MAXVAL(SQRT(ABS(diag - c0_evals**2)))
714 : END IF
715 3686 : DEALLOCATE (diag)
716 3686 : CALL cp_fm_release(matrix_s1)
717 3686 : CALL cp_fm_release(matrix_shc0)
718 : ! we'll only use the energy gap, if our estimate of the error on the eigenvalues
719 : ! is small enough. A large error combined with a small energy gap would otherwise lead to
720 : ! an aggressive but bad preconditioner. Only when the error is small (MD), we can precondition
721 : ! aggressively
722 3686 : preconditioner_env%energy_gap = MAX(energy_gap, error_estimate*fudge_factor)
723 3686 : CALL copy_dbcsr_to_fm(matrix_h, matrix_tmp)
724 3686 : matrix_pre => preconditioner_env%fm
725 3686 : CALL cp_fm_uplo_to_full(matrix_tmp, matrix_pre)
726 : ! tmp = H ( 1 - PS )
727 3686 : CALL parallel_gemm('N', 'T', n, n, k, -1.0_dp, matrix_hc0, matrix_sc0, 1.0_dp, matrix_tmp)
728 :
729 : CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=k, ncol_global=n, &
730 : context=preconditioner_env%ctxt, &
731 3686 : para_env=preconditioner_env%para_env)
732 3686 : CALL cp_fm_create(matrix_left, fm_struct_tmp, name="matrix_left")
733 3686 : CALL cp_fm_struct_release(fm_struct_tmp)
734 3686 : CALL parallel_gemm('T', 'N', k, n, n, 1.0_dp, matrix_c0, matrix_tmp, 0.0_dp, matrix_left)
735 : ! tmp = (1 - PS)^T H (1-PS)
736 3686 : CALL parallel_gemm('N', 'N', n, n, k, -1.0_dp, matrix_sc0, matrix_left, 1.0_dp, matrix_tmp)
737 3686 : CALL cp_fm_release(matrix_left)
738 :
739 7262 : ALLOCATE (shifted_evals(k))
740 3686 : lambda = lambda_base + error_estimate
741 3686 : IF (PRESENT(common_reference)) THEN
742 36 : shifted_evals = common_reference - lambda
743 : ELSE
744 21032 : shifted_evals = c0_evals - lambda
745 : END IF
746 3686 : CALL cp_fm_to_fm(matrix_sc0, matrix_hc0)
747 3686 : CALL cp_fm_column_scale(matrix_hc0, shifted_evals)
748 3686 : CALL parallel_gemm('N', 'T', n, n, k, 1.0_dp, matrix_hc0, matrix_sc0, 1.0_dp, matrix_tmp)
749 :
750 : ! 2) diagonalize this operator
751 3686 : SELECT CASE (preconditioner_env%cholesky_use)
752 : CASE (cholesky_inverse)
753 : CALL cp_fm_triangular_multiply(ortho, matrix_tmp, side="R", transpose_tr=.FALSE., &
754 0 : invert_tr=.FALSE., uplo_tr="U", n_rows=n, n_cols=n, alpha=1.0_dp)
755 : CALL cp_fm_triangular_multiply(ortho, matrix_tmp, side="L", transpose_tr=.TRUE., &
756 0 : invert_tr=.FALSE., uplo_tr="U", n_rows=n, n_cols=n, alpha=1.0_dp)
757 : CASE (cholesky_reduce)
758 3686 : CALL cp_fm_cholesky_reduce(matrix_tmp, ortho)
759 : END SELECT
760 3686 : CALL choose_eigv_solver(matrix_tmp, matrix_pre, preconditioner_env%full_evals)
761 3686 : SELECT CASE (preconditioner_env%cholesky_use)
762 : CASE (cholesky_inverse)
763 : CALL cp_fm_triangular_multiply(ortho, matrix_pre, side="L", transpose_tr=.FALSE., &
764 0 : invert_tr=.FALSE., uplo_tr="U", n_rows=n, n_cols=n, alpha=1.0_dp)
765 0 : CALL cp_fm_to_fm(matrix_pre, matrix_tmp)
766 : CASE (cholesky_reduce)
767 3686 : CALL cp_fm_cholesky_restore(matrix_pre, n, ortho, matrix_tmp, "SOLVE")
768 7372 : CALL cp_fm_to_fm(matrix_tmp, matrix_pre)
769 : END SELECT
770 :
771 : ! test that the subspace remained conserved
772 : IF (.FALSE.) THEN
773 : CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=k, ncol_global=k, &
774 : context=preconditioner_env%ctxt, &
775 : para_env=preconditioner_env%para_env)
776 : CALL cp_fm_create(matrix_s1, fm_struct_tmp, name="matrix_s1")
777 : CALL cp_fm_create(matrix_s2, fm_struct_tmp, name="matrix_s2")
778 : CALL cp_fm_struct_release(fm_struct_tmp)
779 : ALLOCATE (norms(k))
780 : CALL parallel_gemm('T', 'N', k, k, n, 1.0_dp, matrix_sc0, matrix_tmp, 0.0_dp, matrix_s1)
781 : CALL choose_eigv_solver(matrix_s1, matrix_s2, norms)
782 : WRITE (*, *) "matrix norm deviation (should be close to zero): ", MAXVAL(ABS(ABS(norms) - 1.0_dp))
783 : DEALLOCATE (norms)
784 : CALL cp_fm_release(matrix_s1)
785 : CALL cp_fm_release(matrix_s2)
786 : END IF
787 :
788 : ! 3) replace the lowest k evals and evecs with what they should be
789 3686 : IF (PRESENT(common_reference)) THEN
790 36 : preconditioner_env%occ_evals = common_reference
791 36 : preconditioner_env%full_evals(1:k) = common_reference
792 : ELSE
793 21032 : preconditioner_env%occ_evals = c0_evals
794 : ! This choice makes FULL_ALL constant when applied to sc0 (see apply_full_all).
795 21032 : preconditioner_env%full_evals(1:k) = c0_evals
796 : END IF
797 3686 : CALL cp_fm_to_fm(matrix_c0, matrix_pre, k, 1, 1)
798 :
799 3686 : CALL cp_fm_release(matrix_sc0)
800 3686 : CALL cp_fm_release(matrix_hc0)
801 3686 : CALL cp_fm_release(ortho)
802 3686 : CALL cp_fm_release(matrix_tmp)
803 3686 : DEALLOCATE (shifted_evals)
804 3686 : CALL timestop(handle)
805 :
806 29488 : END SUBROUTINE make_full_all
807 :
808 : ! **************************************************************************************************
809 : !> \brief full all in the orthonormal basis
810 : !> \param preconditioner_env ...
811 : !> \param matrix_c0 ...
812 : !> \param matrix_h ...
813 : !> \param c0_evals ...
814 : !> \param energy_gap ...
815 : !> \param common_reference optional scalar occupied-space level used by FERMI_LOW_RANK
816 : ! **************************************************************************************************
817 80 : SUBROUTINE make_full_all_ortho(preconditioner_env, matrix_c0, matrix_h, c0_evals, energy_gap, &
818 : common_reference)
819 :
820 : TYPE(preconditioner_type) :: preconditioner_env
821 : TYPE(cp_fm_type), INTENT(IN) :: matrix_c0
822 : TYPE(dbcsr_type), POINTER :: matrix_h
823 : REAL(KIND=dp), DIMENSION(:) :: c0_evals
824 : REAL(KIND=dp) :: energy_gap
825 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: common_reference
826 :
827 : CHARACTER(len=*), PARAMETER :: routineN = 'make_full_all_ortho'
828 : REAL(KIND=dp), PARAMETER :: fudge_factor = 0.25_dp, &
829 : lambda_base = 10.0_dp
830 :
831 : INTEGER :: handle, k, n
832 : REAL(KIND=dp) :: error_estimate, lambda
833 80 : REAL(KIND=dp), DIMENSION(:), POINTER :: diag, norms, shifted_evals
834 : TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp
835 : TYPE(cp_fm_type) :: matrix_hc0, matrix_left, matrix_s1, &
836 : matrix_s2, matrix_sc0, matrix_tmp
837 : TYPE(cp_fm_type), POINTER :: matrix_pre
838 :
839 80 : CALL timeset(routineN, handle)
840 :
841 80 : IF (ASSOCIATED(preconditioner_env%fm)) THEN
842 0 : CALL cp_fm_release(preconditioner_env%fm)
843 0 : DEALLOCATE (preconditioner_env%fm)
844 : NULLIFY (preconditioner_env%fm)
845 : END IF
846 80 : CALL cp_fm_get_info(matrix_c0, nrow_global=n, ncol_global=k)
847 : CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=n, ncol_global=n, &
848 : context=preconditioner_env%ctxt, &
849 80 : para_env=preconditioner_env%para_env)
850 80 : ALLOCATE (preconditioner_env%fm)
851 80 : CALL cp_fm_create(preconditioner_env%fm, fm_struct_tmp, name="preconditioner_env%fm")
852 80 : CALL cp_fm_create(matrix_tmp, fm_struct_tmp, name="matrix_tmp")
853 80 : CALL cp_fm_struct_release(fm_struct_tmp)
854 240 : ALLOCATE (preconditioner_env%full_evals(n))
855 240 : ALLOCATE (preconditioner_env%occ_evals(k))
856 :
857 : ! 1) Construct a new H matrix, which has the current C0 as eigenvectors,
858 : ! possibly shifted by an amount lambda,
859 : ! and the same spectrum as the original H matrix in the space orthogonal to the C0
860 : ! with P=C0 C0 ^ T
861 : ! (1 - PS)^T H (1-PS) + (PS)^T (H - lambda S ) (PS)
862 : ! we exploit that the C0 are already the ritz states of H
863 80 : CALL cp_fm_create(matrix_sc0, matrix_c0%matrix_struct, name="sc0")
864 80 : CALL cp_fm_to_fm(matrix_c0, matrix_sc0)
865 80 : CALL cp_fm_create(matrix_hc0, matrix_c0%matrix_struct, name="hc0")
866 80 : CALL cp_dbcsr_sm_fm_multiply(matrix_h, matrix_c0, matrix_hc0, k)
867 :
868 : ! An aside, try to estimate the error on the ritz values, we'll need it later on
869 : CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=k, ncol_global=k, &
870 : context=preconditioner_env%ctxt, &
871 80 : para_env=preconditioner_env%para_env)
872 80 : CALL cp_fm_create(matrix_s1, fm_struct_tmp, name="matrix_s1")
873 80 : CALL cp_fm_struct_release(fm_struct_tmp)
874 : ! since we only use diagonal elements this is a bit of a waste
875 80 : CALL parallel_gemm('T', 'N', k, k, n, 1.0_dp, matrix_hc0, matrix_hc0, 0.0_dp, matrix_s1)
876 160 : ALLOCATE (diag(k))
877 80 : CALL cp_fm_get_diag(matrix_s1, diag)
878 80 : IF (PRESENT(common_reference)) THEN
879 : ! ||(1-CC^T)H C||_F is invariant under C -> C U.
880 0 : error_estimate = SQRT(MAX(0.0_dp, SUM(diag) - SUM(c0_evals**2)))
881 : ELSE
882 826 : error_estimate = MAXVAL(SQRT(ABS(diag - c0_evals**2)))
883 : END IF
884 80 : DEALLOCATE (diag)
885 80 : CALL cp_fm_release(matrix_s1)
886 : ! we'll only use the energy gap, if our estimate of the error on the eigenvalues
887 : ! is small enough. A large error combined with a small energy gap would otherwise lead to
888 : ! an aggressive but bad preconditioner. Only when the error is small (MD), we can precondition
889 : ! aggressively
890 80 : preconditioner_env%energy_gap = MAX(energy_gap, error_estimate*fudge_factor)
891 :
892 80 : matrix_pre => preconditioner_env%fm
893 80 : CALL copy_dbcsr_to_fm(matrix_h, matrix_tmp)
894 80 : CALL cp_fm_uplo_to_full(matrix_tmp, matrix_pre)
895 : ! tmp = H ( 1 - PS )
896 80 : CALL parallel_gemm('N', 'T', n, n, k, -1.0_dp, matrix_hc0, matrix_sc0, 1.0_dp, matrix_tmp)
897 :
898 : CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=k, ncol_global=n, &
899 : context=preconditioner_env%ctxt, &
900 80 : para_env=preconditioner_env%para_env)
901 80 : CALL cp_fm_create(matrix_left, fm_struct_tmp, name="matrix_left")
902 80 : CALL cp_fm_struct_release(fm_struct_tmp)
903 80 : CALL parallel_gemm('T', 'N', k, n, n, 1.0_dp, matrix_c0, matrix_tmp, 0.0_dp, matrix_left)
904 : ! tmp = (1 - PS)^T H (1-PS)
905 80 : CALL parallel_gemm('N', 'N', n, n, k, -1.0_dp, matrix_sc0, matrix_left, 1.0_dp, matrix_tmp)
906 80 : CALL cp_fm_release(matrix_left)
907 :
908 160 : ALLOCATE (shifted_evals(k))
909 80 : lambda = lambda_base + error_estimate
910 80 : IF (PRESENT(common_reference)) THEN
911 0 : shifted_evals = common_reference - lambda
912 : ELSE
913 826 : shifted_evals = c0_evals - lambda
914 : END IF
915 80 : CALL cp_fm_to_fm(matrix_sc0, matrix_hc0)
916 80 : CALL cp_fm_column_scale(matrix_hc0, shifted_evals)
917 80 : CALL parallel_gemm('N', 'T', n, n, k, 1.0_dp, matrix_hc0, matrix_sc0, 1.0_dp, matrix_tmp)
918 :
919 : ! 2) diagonalize this operator
920 80 : CALL choose_eigv_solver(matrix_tmp, matrix_pre, preconditioner_env%full_evals)
921 :
922 : ! test that the subspace remained conserved
923 : IF (.FALSE.) THEN
924 : CALL cp_fm_to_fm(matrix_pre, matrix_tmp)
925 : CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=k, ncol_global=k, &
926 : context=preconditioner_env%ctxt, &
927 : para_env=preconditioner_env%para_env)
928 : CALL cp_fm_create(matrix_s1, fm_struct_tmp, name="matrix_s1")
929 : CALL cp_fm_create(matrix_s2, fm_struct_tmp, name="matrix_s2")
930 : CALL cp_fm_struct_release(fm_struct_tmp)
931 : ALLOCATE (norms(k))
932 : CALL parallel_gemm('T', 'N', k, k, n, 1.0_dp, matrix_sc0, matrix_tmp, 0.0_dp, matrix_s1)
933 : CALL choose_eigv_solver(matrix_s1, matrix_s2, norms)
934 :
935 : WRITE (*, *) "matrix norm deviation (should be close to zero): ", MAXVAL(ABS(ABS(norms) - 1.0_dp))
936 : DEALLOCATE (norms)
937 : CALL cp_fm_release(matrix_s1)
938 : CALL cp_fm_release(matrix_s2)
939 : END IF
940 :
941 : ! 3) replace the lowest k evals and evecs with what they should be
942 80 : IF (PRESENT(common_reference)) THEN
943 0 : preconditioner_env%occ_evals = common_reference
944 0 : preconditioner_env%full_evals(1:k) = common_reference
945 : ELSE
946 826 : preconditioner_env%occ_evals = c0_evals
947 : ! This choice makes FULL_ALL constant when applied to sc0 (see apply_full_all).
948 826 : preconditioner_env%full_evals(1:k) = c0_evals
949 : END IF
950 80 : CALL cp_fm_to_fm(matrix_c0, matrix_pre, k, 1, 1)
951 :
952 80 : CALL cp_fm_release(matrix_sc0)
953 80 : CALL cp_fm_release(matrix_hc0)
954 80 : CALL cp_fm_release(matrix_tmp)
955 80 : DEALLOCATE (shifted_evals)
956 :
957 80 : CALL timestop(handle)
958 :
959 560 : END SUBROUTINE make_full_all_ortho
960 :
961 : ! **************************************************************************************************
962 : !> \brief generates a preconditioner matrix H-lambda S+(SC)(2.0*CT*H*C+delta)(SC)^T
963 : !> for later inversion.
964 : !> H is the Kohn Sham matrix
965 : !> lambda*S shifts the spectrum of the generalized form up by lambda
966 : !> the last term only shifts the occupied space (reversing them in energy order)
967 : !> This form is implicitly multiplied from both sides by S^0.5
968 : !> This ensures we precondition the correct quantity
969 : !> Before this reads S^-0.5 H S^-0.5 + lambda + (S^0.5 C)shifts(S^0.5 C)T
970 : !> which might be a bit more obvious
971 : !> Replaced the old full_single_inverse at revision 14616
972 : !> \param preconditioner_env the preconditioner env
973 : !> \param matrix_c0 the MO coefficient matrix (fm)
974 : !> \param matrix_h Kohn-Sham matrix (dbcsr)
975 : !> \param energy_gap an additional shift in lambda=-E_homo+energy_gap
976 : !> \param matrix_s the overlap matrix if not orthonormal (dbcsr, optional)
977 : ! **************************************************************************************************
978 4501 : SUBROUTINE make_full_single_inverse(preconditioner_env, matrix_c0, matrix_h, energy_gap, matrix_s)
979 : TYPE(preconditioner_type) :: preconditioner_env
980 : TYPE(cp_fm_type), INTENT(IN) :: matrix_c0
981 : TYPE(dbcsr_type), POINTER :: matrix_h
982 : REAL(KIND=dp) :: energy_gap
983 : TYPE(dbcsr_type), OPTIONAL, POINTER :: matrix_s
984 :
985 : CHARACTER(len=*), PARAMETER :: routineN = 'make_full_single_inverse'
986 :
987 : INTEGER :: handle, k, n
988 : REAL(KIND=dp) :: max_ev, min_ev, pre_shift
989 : TYPE(arnoldi_env_type) :: arnoldi_env
990 4501 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrices
991 : TYPE(dbcsr_type), TARGET :: dbcsr_cThc, dbcsr_hc, dbcsr_sc, mo_dbcsr
992 :
993 4501 : CALL timeset(routineN, handle)
994 :
995 : ! Allocate all working matrices needed
996 4501 : CALL cp_fm_get_info(matrix_c0, nrow_global=n, ncol_global=k)
997 : ! copy the fm MO's to a sparse matrix, can be solved better if the sparse version is already present
998 : ! but for the time beeing this will do
999 4501 : CALL cp_fm_to_dbcsr_row_template(mo_dbcsr, matrix_c0, matrix_h)
1000 4501 : CALL dbcsr_create(dbcsr_sc, template=mo_dbcsr)
1001 4501 : CALL dbcsr_create(dbcsr_hc, template=mo_dbcsr)
1002 4501 : CALL cp_dbcsr_m_by_n_from_template(dbcsr_cThc, matrix_h, k, k, sym=dbcsr_type_symmetric)
1003 :
1004 : ! Check whether the output matrix was already created, if not do it now
1005 4501 : IF (.NOT. ASSOCIATED(preconditioner_env%sparse_matrix)) THEN
1006 4501 : ALLOCATE (preconditioner_env%sparse_matrix)
1007 : END IF
1008 :
1009 : ! Put the first term of the preconditioner (H) into the output matrix
1010 4501 : CALL dbcsr_copy(preconditioner_env%sparse_matrix, matrix_h)
1011 :
1012 : ! Precompute some matrices
1013 : ! S*C, if orthonormal this will be simply C so a copy will do
1014 4501 : IF (PRESENT(matrix_s)) THEN
1015 4113 : CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_s, mo_dbcsr, 0.0_dp, dbcsr_sc)
1016 : ELSE
1017 388 : CALL dbcsr_copy(dbcsr_sc, mo_dbcsr)
1018 : END IF
1019 :
1020 : !----------------------------compute the occupied subspace and shift it ------------------------------------
1021 : ! cT*H*C which will be used to shift the occupied states to 0
1022 4501 : CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_h, mo_dbcsr, 0.0_dp, dbcsr_hc)
1023 4501 : CALL dbcsr_multiply("T", "N", 1.0_dp, mo_dbcsr, dbcsr_hc, 0.0_dp, dbcsr_cThc)
1024 :
1025 : ! Compute the Energy of the HOMO. We will use this as a reference energy
1026 9002 : ALLOCATE (matrices(1))
1027 4501 : matrices(1)%matrix => dbcsr_cThc
1028 : CALL setup_arnoldi_env(arnoldi_env, matrices, max_iter=20, threshold=1.0E-3_dp, selection_crit=2, &
1029 4501 : nval_request=1, nrestarts=8, generalized_ev=.FALSE., iram=.FALSE.)
1030 4501 : IF (ASSOCIATED(preconditioner_env%max_ev_vector)) THEN
1031 2444 : CALL set_arnoldi_initial_vector(arnoldi_env, preconditioner_env%max_ev_vector)
1032 : END IF
1033 4501 : CALL arnoldi_ev(matrices, arnoldi_env)
1034 4501 : max_ev = REAL(get_selected_ritz_val(arnoldi_env, 1), dp)
1035 :
1036 : ! save the ev as guess for the next time
1037 4501 : IF (.NOT. ASSOCIATED(preconditioner_env%max_ev_vector)) ALLOCATE (preconditioner_env%max_ev_vector)
1038 4501 : CALL get_selected_ritz_vector(arnoldi_env, 1, matrices(1)%matrix, preconditioner_env%max_ev_vector)
1039 4501 : CALL deallocate_arnoldi_env(arnoldi_env)
1040 4501 : DEALLOCATE (matrices)
1041 :
1042 : ! Lets shift the occupied states a bit further up, -1.0 because we gonna subtract it from H
1043 4501 : CALL dbcsr_add_on_diag(dbcsr_cThc, -0.5_dp)
1044 : ! Get the AO representation of the shift (see above why S is needed), W-matrix like object
1045 4501 : CALL dbcsr_multiply("N", "N", 2.0_dp, dbcsr_sc, dbcsr_cThc, 0.0_dp, dbcsr_hc)
1046 4501 : CALL dbcsr_multiply("N", "T", -1.0_dp, dbcsr_hc, dbcsr_sc, 1.0_dp, preconditioner_env%sparse_matrix)
1047 :
1048 : !-------------------------------------compute eigenvalues of H ----------------------------------------------
1049 : ! Setup the arnoldi procedure to compute the lowest ev. if S is present this has to be the generalized ev
1050 4501 : IF (PRESENT(matrix_s)) THEN
1051 12339 : ALLOCATE (matrices(2))
1052 4113 : matrices(1)%matrix => preconditioner_env%sparse_matrix
1053 4113 : matrices(2)%matrix => matrix_s
1054 : CALL setup_arnoldi_env(arnoldi_env, matrices, max_iter=20, threshold=2.0E-2_dp, selection_crit=3, &
1055 4113 : nval_request=1, nrestarts=21, generalized_ev=.TRUE., iram=.FALSE.)
1056 : ELSE
1057 776 : ALLOCATE (matrices(1))
1058 388 : matrices(1)%matrix => preconditioner_env%sparse_matrix
1059 : CALL setup_arnoldi_env(arnoldi_env, matrices, max_iter=20, threshold=2.0E-2_dp, selection_crit=3, &
1060 388 : nval_request=1, nrestarts=8, generalized_ev=.FALSE., iram=.FALSE.)
1061 : END IF
1062 4501 : IF (ASSOCIATED(preconditioner_env%min_ev_vector)) THEN
1063 2444 : CALL set_arnoldi_initial_vector(arnoldi_env, preconditioner_env%min_ev_vector)
1064 : END IF
1065 :
1066 : ! compute the LUMO energy
1067 4501 : CALL arnoldi_ev(matrices, arnoldi_env)
1068 4501 : min_eV = REAL(get_selected_ritz_val(arnoldi_env, 1), dp)
1069 :
1070 : ! save the lumo vector for restarting in the next step
1071 4501 : IF (.NOT. ASSOCIATED(preconditioner_env%min_ev_vector)) ALLOCATE (preconditioner_env%min_ev_vector)
1072 4501 : CALL get_selected_ritz_vector(arnoldi_env, 1, matrices(1)%matrix, preconditioner_env%min_ev_vector)
1073 4501 : CALL deallocate_arnoldi_env(arnoldi_env)
1074 4501 : DEALLOCATE (matrices)
1075 :
1076 : !-------------------------------------compute eigenvalues of H ----------------------------------------------
1077 : ! Shift the Lumo to the 1.5*the computed energy_gap or the external energy gap value
1078 : ! The factor 1.5 is determined by trying. If the LUMO is positive, enough, just leave it alone
1079 4501 : pre_shift = MAX(1.5_dp*(min_ev - max_ev), energy_gap)
1080 4501 : IF (min_ev < pre_shift) THEN
1081 4483 : pre_shift = pre_shift - min_ev
1082 : ELSE
1083 18 : pre_shift = 0.0_dp
1084 : END IF
1085 4501 : IF (PRESENT(matrix_s)) THEN
1086 4113 : CALL dbcsr_add(preconditioner_env%sparse_matrix, matrix_s, 1.0_dp, pre_shift)
1087 : ELSE
1088 388 : CALL dbcsr_add_on_diag(preconditioner_env%sparse_matrix, pre_shift)
1089 : END IF
1090 :
1091 4501 : CALL dbcsr_release(mo_dbcsr)
1092 4501 : CALL dbcsr_release(dbcsr_hc)
1093 4501 : CALL dbcsr_release(dbcsr_sc)
1094 4501 : CALL dbcsr_release(dbcsr_cThc)
1095 :
1096 4501 : CALL timestop(handle)
1097 :
1098 4501 : END SUBROUTINE make_full_single_inverse
1099 :
1100 : ! **************************************************************************************************
1101 : !> \brief Build a gauge-covariant FULL_SINGLE_INVERSE operator for a complex k-point channel.
1102 : !> \param preconditioner_env preconditioner storage
1103 : !> \param matrix_c0 complex reference orbitals, C^H S C = I
1104 : !> \param matrix_h complex Hermitian k-point Hamiltonian
1105 : !> \param matrix_s complex Hermitian k-point overlap
1106 : !> \param energy_gap lower spectral bound of the positive operator
1107 : ! **************************************************************************************************
1108 64 : SUBROUTINE make_complex_full_single_inverse(preconditioner_env, matrix_c0, matrix_h, matrix_s, &
1109 : energy_gap)
1110 :
1111 : TYPE(preconditioner_type) :: preconditioner_env
1112 : TYPE(cp_cfm_type), INTENT(IN) :: matrix_c0, matrix_h, matrix_s
1113 : REAL(KIND=dp), INTENT(IN) :: energy_gap
1114 :
1115 : CHARACTER(len=*), PARAMETER :: routineN = 'make_complex_full_single_inverse'
1116 :
1117 : INTEGER :: handle, i, k, n
1118 : REAL(KIND=dp) :: max_ev, min_ev, pre_shift, target_edge
1119 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: occupied_evals, operator_evals
1120 : TYPE(cp_cfm_type) :: matrix_chc, matrix_diag_a, matrix_diag_b, matrix_diag_evec, &
1121 : matrix_diag_work, matrix_hc0, matrix_occ_diag_a, matrix_occ_diag_b, matrix_occ_diag_evec, &
1122 : matrix_occ_diag_work, matrix_operator, matrix_sc0, matrix_sc_chc
1123 :
1124 64 : CALL timeset(routineN, handle)
1125 :
1126 64 : CALL cp_cfm_get_info(matrix_c0, nrow_global=n, ncol_global=k)
1127 64 : CPASSERT(n > 0)
1128 64 : CPASSERT(k > 0 .AND. k <= n)
1129 64 : CPASSERT(energy_gap > 0.0_dp)
1130 :
1131 64 : CALL cp_cfm_create(matrix_hc0, matrix_c0%matrix_struct, name='complex FULL_SINGLE HC')
1132 64 : CALL cp_cfm_create(matrix_sc0, matrix_c0%matrix_struct, name='complex FULL_SINGLE SC')
1133 : CALL cp_cfm_create(matrix_chc, matrix_c0%matrix_struct, nrow=k, ncol=k, &
1134 64 : name='complex FULL_SINGLE CHC')
1135 : CALL cp_cfm_create(matrix_sc_chc, matrix_c0%matrix_struct, &
1136 64 : name='complex FULL_SINGLE SC CHC')
1137 : CALL cp_cfm_create(matrix_operator, matrix_h%matrix_struct, &
1138 64 : name='complex FULL_SINGLE operator')
1139 : CALL cp_cfm_create(matrix_diag_a, matrix_h%matrix_struct, &
1140 64 : name='complex FULL_SINGLE diagonalization A')
1141 : CALL cp_cfm_create(matrix_diag_b, matrix_s%matrix_struct, &
1142 64 : name='complex FULL_SINGLE diagonalization B')
1143 : CALL cp_cfm_create(matrix_diag_evec, matrix_h%matrix_struct, nrow=n, ncol=1, &
1144 64 : name='complex FULL_SINGLE eigenvectors')
1145 : CALL cp_cfm_create(matrix_diag_work, matrix_h%matrix_struct, &
1146 64 : name='complex FULL_SINGLE diagonalization work')
1147 : CALL cp_cfm_create(matrix_occ_diag_a, matrix_c0%matrix_struct, nrow=k, ncol=k, &
1148 64 : name='complex FULL_SINGLE occupied diagonalization A')
1149 : CALL cp_cfm_create(matrix_occ_diag_b, matrix_c0%matrix_struct, nrow=k, ncol=k, &
1150 64 : name='complex FULL_SINGLE occupied diagonalization B', set_zero=.TRUE.)
1151 : CALL cp_cfm_create(matrix_occ_diag_evec, matrix_c0%matrix_struct, nrow=k, ncol=k, &
1152 64 : name='complex FULL_SINGLE occupied eigenvectors')
1153 : CALL cp_cfm_create(matrix_occ_diag_work, matrix_c0%matrix_struct, nrow=k, ncol=k, &
1154 64 : name='complex FULL_SINGLE occupied diagonalization work')
1155 :
1156 64 : CALL cp_cfm_gemm('N', 'N', n, k, n, z_one, matrix_h, matrix_c0, z_zero, matrix_hc0)
1157 64 : CALL cp_cfm_gemm('N', 'N', n, k, n, z_one, matrix_s, matrix_c0, z_zero, matrix_sc0)
1158 64 : CALL cp_cfm_gemm('C', 'N', k, k, n, z_one, matrix_c0, matrix_hc0, z_zero, matrix_chc)
1159 :
1160 : ! A = H - 2*S*C*(C^H*H*C)*C^H*S + S*C*C^H*S is invariant under C -> C*U.
1161 64 : CALL cp_cfm_to_cfm(matrix_h, matrix_operator)
1162 : CALL cp_cfm_gemm('N', 'N', n, k, k, z_one, matrix_sc0, matrix_chc, &
1163 64 : z_zero, matrix_sc_chc)
1164 : CALL cp_cfm_gemm('N', 'C', n, n, k, -2.0_dp*z_one, matrix_sc_chc, matrix_sc0, &
1165 64 : z_one, matrix_operator)
1166 : CALL cp_cfm_gemm('N', 'C', n, n, k, z_one, matrix_sc0, matrix_sc0, &
1167 64 : z_one, matrix_operator)
1168 :
1169 192 : ALLOCATE (occupied_evals(k), operator_evals(1))
1170 64 : CALL cp_cfm_to_cfm(matrix_chc, matrix_occ_diag_a)
1171 484 : DO i = 1, k
1172 484 : CALL cp_cfm_set_element(matrix_occ_diag_b, i, i, z_one)
1173 : END DO
1174 : CALL cp_cfm_geeig(matrix_occ_diag_a, matrix_occ_diag_b, matrix_occ_diag_evec, &
1175 64 : occupied_evals, matrix_occ_diag_work)
1176 484 : max_ev = MAXVAL(occupied_evals)
1177 64 : CALL cp_cfm_to_cfm(matrix_operator, matrix_diag_a)
1178 64 : CALL cp_cfm_to_cfm(matrix_s, matrix_diag_b)
1179 : CALL cp_cfm_geeig(matrix_diag_a, matrix_diag_b, matrix_diag_evec, &
1180 64 : operator_evals, matrix_diag_work, lowest_subset=.TRUE.)
1181 64 : min_ev = operator_evals(1)
1182 64 : target_edge = MAX(1.5_dp*(min_ev - max_ev), energy_gap)
1183 64 : pre_shift = MAX(0.0_dp, target_edge - min_ev)
1184 64 : IF (pre_shift > 0.0_dp) THEN
1185 : CALL cp_cfm_scale_and_add(z_one, matrix_operator, &
1186 64 : CMPLX(pre_shift, 0.0_dp, KIND=dp), matrix_s)
1187 : END IF
1188 :
1189 : CALL store_complex_inverse(preconditioner_env, matrix_operator, &
1190 64 : ot_precond_full_single_inverse, energy_gap)
1191 :
1192 64 : DEALLOCATE (operator_evals, occupied_evals)
1193 64 : CALL cp_cfm_release(matrix_occ_diag_work)
1194 64 : CALL cp_cfm_release(matrix_occ_diag_evec)
1195 64 : CALL cp_cfm_release(matrix_occ_diag_b)
1196 64 : CALL cp_cfm_release(matrix_occ_diag_a)
1197 64 : CALL cp_cfm_release(matrix_diag_work)
1198 64 : CALL cp_cfm_release(matrix_diag_evec)
1199 64 : CALL cp_cfm_release(matrix_diag_b)
1200 64 : CALL cp_cfm_release(matrix_diag_a)
1201 64 : CALL cp_cfm_release(matrix_operator)
1202 64 : CALL cp_cfm_release(matrix_sc_chc)
1203 64 : CALL cp_cfm_release(matrix_chc)
1204 64 : CALL cp_cfm_release(matrix_sc0)
1205 64 : CALL cp_cfm_release(matrix_hc0)
1206 :
1207 64 : CALL timestop(handle)
1208 :
1209 256 : END SUBROUTINE make_complex_full_single_inverse
1210 :
1211 : ! **************************************************************************************************
1212 : !> \brief Build the inverse complex overlap preconditioner.
1213 : !> \param preconditioner_env preconditioner storage
1214 : !> \param matrix_s complex Hermitian k-point overlap
1215 : ! **************************************************************************************************
1216 92 : SUBROUTINE make_complex_full_s_inverse(preconditioner_env, matrix_s)
1217 :
1218 : TYPE(preconditioner_type) :: preconditioner_env
1219 : TYPE(cp_cfm_type), INTENT(IN) :: matrix_s
1220 :
1221 92 : CALL store_complex_inverse(preconditioner_env, matrix_s, ot_precond_s_inverse, 0.0_dp)
1222 :
1223 92 : END SUBROUTINE make_complex_full_s_inverse
1224 :
1225 : ! **************************************************************************************************
1226 : !> \brief Build the inverse complex kinetic-plus-overlap preconditioner.
1227 : !> \param preconditioner_env preconditioner storage
1228 : !> \param matrix_t complex Hermitian k-point kinetic operator
1229 : !> \param matrix_s complex Hermitian k-point overlap
1230 : !> \param energy_gap non-negative overlap shift
1231 : ! **************************************************************************************************
1232 124 : SUBROUTINE make_complex_full_kinetic(preconditioner_env, matrix_t, matrix_s, energy_gap)
1233 :
1234 : TYPE(preconditioner_type) :: preconditioner_env
1235 : TYPE(cp_cfm_type), INTENT(IN) :: matrix_t, matrix_s
1236 : REAL(KIND=dp), INTENT(IN) :: energy_gap
1237 :
1238 : REAL(KIND=dp) :: shift
1239 : TYPE(cp_cfm_type) :: matrix_operator
1240 :
1241 124 : shift = MAX(0.0_dp, energy_gap)
1242 : CALL cp_cfm_create(matrix_operator, matrix_t%matrix_struct, &
1243 124 : name='complex FULL_KINETIC operator')
1244 124 : CALL cp_cfm_to_cfm(matrix_t, matrix_operator)
1245 124 : CALL cp_cfm_scale_and_add(z_one, matrix_operator, CMPLX(shift, 0.0_dp, KIND=dp), matrix_s)
1246 : CALL store_complex_inverse(preconditioner_env, matrix_operator, &
1247 124 : ot_precond_full_kinetic, energy_gap)
1248 124 : CALL cp_cfm_release(matrix_operator)
1249 :
1250 124 : END SUBROUTINE make_complex_full_kinetic
1251 :
1252 : ! **************************************************************************************************
1253 : !> \brief Store an explicitly inverted positive complex Hermitian operator.
1254 : !> \param preconditioner_env preconditioner storage
1255 : !> \param matrix_operator positive complex Hermitian operator
1256 : !> \param preconditioner_kind selected OT preconditioner
1257 : !> \param energy_gap configured spectral shift
1258 : ! **************************************************************************************************
1259 840 : SUBROUTINE store_complex_inverse(preconditioner_env, matrix_operator, &
1260 : preconditioner_kind, energy_gap)
1261 :
1262 : TYPE(preconditioner_type) :: preconditioner_env
1263 : TYPE(cp_cfm_type), INTENT(IN) :: matrix_operator
1264 : INTEGER, INTENT(IN) :: preconditioner_kind
1265 : REAL(KIND=dp), INTENT(IN) :: energy_gap
1266 :
1267 : INTEGER :: info
1268 :
1269 280 : IF (ASSOCIATED(preconditioner_env%complex_fm)) THEN
1270 0 : CALL cp_cfm_release(preconditioner_env%complex_fm)
1271 0 : DEALLOCATE (preconditioner_env%complex_fm)
1272 : END IF
1273 280 : ALLOCATE (preconditioner_env%complex_fm)
1274 : CALL cp_cfm_create(preconditioner_env%complex_fm, matrix_operator%matrix_struct, &
1275 280 : name='complex inverse preconditioner')
1276 280 : CALL cp_cfm_to_cfm(matrix_operator, preconditioner_env%complex_fm)
1277 280 : CALL cp_cfm_cholesky_decompose(preconditioner_env%complex_fm, info_out=info)
1278 280 : CPASSERT(info == 0)
1279 280 : CALL cp_cfm_cholesky_invert(preconditioner_env%complex_fm, info_out=info)
1280 280 : CPASSERT(info == 0)
1281 280 : CALL cp_cfm_uplo_to_full(preconditioner_env%complex_fm)
1282 :
1283 280 : preconditioner_env%energy_gap = energy_gap
1284 280 : preconditioner_env%in_use = preconditioner_kind
1285 280 : preconditioner_env%solver = ot_precond_solver_default
1286 :
1287 280 : END SUBROUTINE store_complex_inverse
1288 : END MODULE preconditioner_makes
|
