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