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