(git:f2099e5)
Loading...
Searching...
No Matches
preconditioner_makes.F
Go to the documentation of this file.
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! **************************************************************************************************
15MODULE 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
89CONTAINS
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 SUBROUTINE make_complex_full_all(preconditioner_env, matrix_c0, matrix_h, matrix_s, &
103 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 COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:) :: shifted_evals
115 COMPLEX(KIND=dp), CONTIGUOUS, DIMENSION(:, :), &
116 POINTER :: local_data
117 INTEGER :: handle, j, k, n, ncol_local, nrow_local
118 INTEGER, DIMENSION(:), POINTER :: col_indices
119 REAL(kind=dp) :: error_estimate, lambda
120 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 CALL timeset(routinen, handle)
125
126 CALL cp_cfm_get_info(matrix_c0, nrow_global=n, ncol_global=k)
127 cpassert(n > 0)
128 cpassert(k > 0 .AND. k <= n)
129 cpassert(SIZE(c0_evals) >= k)
130 cpassert(energy_gap > 0.0_dp)
131
132 IF (ASSOCIATED(preconditioner_env%complex_fm)) THEN
133 CALL cp_cfm_release(preconditioner_env%complex_fm)
134 DEALLOCATE (preconditioner_env%complex_fm)
135 END IF
136 IF (ASSOCIATED(preconditioner_env%occ_evals)) DEALLOCATE (preconditioner_env%occ_evals)
137 IF (ASSOCIATED(preconditioner_env%full_evals)) DEALLOCATE (preconditioner_env%full_evals)
138
139 ALLOCATE (preconditioner_env%complex_fm)
140 CALL cp_cfm_create(preconditioner_env%complex_fm, matrix_h%matrix_struct, &
141 name='complex FULL_ALL eigenvectors')
142 ALLOCATE (preconditioner_env%full_evals(n), preconditioner_env%occ_evals(k))
143
144 CALL cp_cfm_create(matrix_hc0, matrix_c0%matrix_struct, name='complex FULL_ALL HC')
145 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 name='complex FULL_ALL CHC')
148 CALL cp_cfm_create(matrix_sc_chc, matrix_c0%matrix_struct, &
149 name='complex FULL_ALL SC CHC')
150 CALL cp_cfm_create(matrix_shifted_sc, matrix_c0%matrix_struct, &
151 name='complex FULL_ALL shifted SC')
152 CALL cp_cfm_create(matrix_tmp, matrix_h%matrix_struct, name='complex FULL_ALL projected H')
153 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 name='complex FULL_ALL diagonalization work')
156
157 CALL cp_cfm_gemm('N', 'N', n, k, n, z_one, matrix_h, matrix_c0, z_zero, matrix_hc0)
158 CALL cp_cfm_gemm('N', 'N', n, k, n, z_one, matrix_s, matrix_c0, z_zero, matrix_sc0)
159 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 ALLOCATE (shifted_evals(k), residual_norm_sq(k))
165 shifted_evals(:) = cmplx(c0_evals(1:k), 0.0_dp, kind=dp)
166 CALL cp_cfm_to_cfm(matrix_sc0, matrix_sc_chc)
167 CALL cp_cfm_column_scale(matrix_sc_chc, shifted_evals)
168 CALL cp_cfm_to_cfm(matrix_hc0, matrix_shifted_sc)
169 CALL cp_cfm_scale_and_add(z_one, matrix_shifted_sc, -z_one, matrix_sc_chc)
170 CALL cp_cfm_to_cfm(matrix_s, matrix_s_copy)
171 CALL cp_cfm_cholesky_decompose(matrix_s_copy)
172 CALL cp_cfm_triangular_multiply(matrix_s_copy, matrix_shifted_sc, side='L', &
173 transa_tr='C', invert_tr=.true., uplo_tr='U')
174
175 residual_norm_sq(:) = 0.0_dp
176 CALL cp_cfm_get_info(matrix_shifted_sc, nrow_local=nrow_local, ncol_local=ncol_local, &
177 col_indices=col_indices, local_data=local_data)
178 DO j = 1, ncol_local
179 residual_norm_sq(col_indices(j)) = residual_norm_sq(col_indices(j)) + &
180 sum(abs(local_data(1:nrow_local, j))**2)
181 END DO
182 CALL preconditioner_env%para_env%sum(residual_norm_sq)
183 error_estimate = sqrt(maxval(residual_norm_sq))
184 preconditioner_env%energy_gap = max(energy_gap, error_estimate*fudge_factor)
185 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 CALL cp_cfm_to_cfm(matrix_h, matrix_tmp)
190 CALL cp_cfm_gemm('N', 'C', n, n, k, -z_one, matrix_hc0, matrix_sc0, z_one, matrix_tmp)
191 CALL cp_cfm_gemm('N', 'C', n, n, k, -z_one, matrix_sc0, matrix_hc0, z_one, matrix_tmp)
192 CALL cp_cfm_gemm('N', 'N', n, k, k, z_one, matrix_sc0, matrix_chc, z_zero, matrix_sc_chc)
193 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 shifted_evals(:) = cmplx(c0_evals(1:k) - lambda, 0.0_dp, kind=dp)
198 CALL cp_cfm_to_cfm(matrix_sc0, matrix_shifted_sc)
199 CALL cp_cfm_column_scale(matrix_shifted_sc, shifted_evals)
200 CALL cp_cfm_gemm('N', 'C', n, n, k, z_one, matrix_shifted_sc, matrix_sc0, z_one, matrix_tmp)
201
202 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 preconditioner_env%full_evals, matrix_diag_work)
205
206 preconditioner_env%occ_evals(:) = c0_evals(1:k)
207 preconditioner_env%full_evals(1:k) = c0_evals(1:k)
208 CALL cp_cfm_to_cfm(matrix_c0, preconditioner_env%complex_fm, k)
209 preconditioner_env%in_use = ot_precond_full_all
210 preconditioner_env%solver = ot_precond_solver_default
211
212 DEALLOCATE (residual_norm_sq, shifted_evals)
213 CALL cp_cfm_release(matrix_diag_work)
214 CALL cp_cfm_release(matrix_s_copy)
215 CALL cp_cfm_release(matrix_tmp)
216 CALL cp_cfm_release(matrix_shifted_sc)
217 CALL cp_cfm_release(matrix_sc_chc)
218 CALL cp_cfm_release(matrix_chc)
219 CALL cp_cfm_release(matrix_sc0)
220 CALL cp_cfm_release(matrix_hc0)
221
222 CALL timestop(handle)
223
224 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 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 COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:) :: scaling
244 INTEGER :: handle, i, n
245 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 CALL timeset(routinen, handle)
251
252 CALL cp_cfm_get_info(matrix_h, nrow_global=n)
253 cpassert(n > 0)
254 cpassert(energy_gap > 0.0_dp)
255
256 IF (ASSOCIATED(preconditioner_env%complex_fm)) THEN
257 CALL cp_cfm_release(preconditioner_env%complex_fm)
258 DEALLOCATE (preconditioner_env%complex_fm)
259 END IF
260 IF (ASSOCIATED(preconditioner_env%occ_evals)) DEALLOCATE (preconditioner_env%occ_evals)
261 IF (ASSOCIATED(preconditioner_env%full_evals)) DEALLOCATE (preconditioner_env%full_evals)
262
263 ALLOCATE (preconditioner_env%complex_fm)
264 CALL cp_cfm_create(preconditioner_env%complex_fm, matrix_h%matrix_struct, &
265 name='complex FULL_SINGLE preconditioner')
266 CALL cp_cfm_create(matrix_diag_a, matrix_h%matrix_struct, &
267 name='complex FULL_SINGLE diagonalization A')
268 CALL cp_cfm_create(matrix_diag_b, matrix_s%matrix_struct, &
269 name='complex FULL_SINGLE diagonalization B')
270 CALL cp_cfm_create(matrix_eigenvectors, matrix_h%matrix_struct, &
271 name='complex FULL_SINGLE eigenvectors')
272 CALL cp_cfm_create(matrix_scaled_eigenvectors, matrix_h%matrix_struct, &
273 name='complex FULL_SINGLE scaled eigenvectors')
274 CALL cp_cfm_create(matrix_work, matrix_h%matrix_struct, &
275 name='complex FULL_SINGLE diagonalization work')
276
277 ALLOCATE (eigenvalues(n), scaling(n))
278 CALL cp_cfm_to_cfm(matrix_h, matrix_diag_a)
279 CALL cp_cfm_to_cfm(matrix_s, matrix_diag_b)
280 CALL cp_cfm_geeig(matrix_diag_a, matrix_diag_b, matrix_eigenvectors, &
281 eigenvalues, matrix_work)
282 DO i = 1, n
283 scaling(i) = cmplx(1.0_dp/max(eigenvalues(i) - energy_homo, energy_gap), &
284 0.0_dp, kind=dp)
285 END DO
286 CALL cp_cfm_to_cfm(matrix_eigenvectors, matrix_scaled_eigenvectors)
287 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 matrix_eigenvectors, z_zero, preconditioner_env%complex_fm)
290
291 preconditioner_env%energy_gap = energy_gap
292 preconditioner_env%in_use = ot_precond_full_single
293 preconditioner_env%solver = ot_precond_solver_default
294
295 DEALLOCATE (scaling, eigenvalues)
296 CALL cp_cfm_release(matrix_work)
297 CALL cp_cfm_release(matrix_scaled_eigenvectors)
298 CALL cp_cfm_release(matrix_eigenvectors)
299 CALL cp_cfm_release(matrix_diag_b)
300 CALL cp_cfm_release(matrix_diag_a)
301
302 CALL timestop(handle)
303
304 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 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 CALL cp_cfm_get_info(matrix_c0, nrow_global=n, ncol_global=k)
332 cpassert(k > 0 .AND. k <= n)
333 cpassert(max_rank == -1 .OR. max_rank > 0)
334 cpassert(spectral_window > 0.0_dp)
335 cpassert(degeneracy_tolerance >= 0.0_dp)
336
337 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 name='complex FERMI_LOW_RANK CHC')
340 CALL cp_cfm_create(matrix_vectors, matrix_chc%matrix_struct, &
341 name='complex FERMI_LOW_RANK active vectors')
342 ALLOCATE (active_evals(k), common_evals(k))
343 CALL cp_cfm_gemm('N', 'N', n, k, n, z_one, matrix_h, matrix_c0, z_zero, matrix_hc0)
344 CALL cp_cfm_gemm('C', 'N', k, k, n, z_one, matrix_c0, matrix_hc0, z_zero, matrix_chc)
345 CALL cp_cfm_heevd(matrix_chc, matrix_vectors, active_evals)
346 preconditioner_env%spectral_reference = maxval(active_evals)
347 common_evals(:) = preconditioner_env%spectral_reference
348
349 CALL make_complex_full_all(preconditioner_env, matrix_c0, matrix_h, matrix_s, &
350 common_evals, energy_gap)
351 preconditioner_env%spectral_window = spectral_window
352 preconditioner_env%spectral_base_scale = 1.0_dp/spectral_window
353 preconditioner_env%spectral_max_rank = max_rank
354 max_rank_eff = max_rank
355 IF (max_rank_eff == -1) max_rank_eff = n - k
356 preconditioner_env%spectral_rank = low_rank_select_rank( &
357 preconditioner_env%full_evals, k, max_rank_eff, degeneracy_tolerance)
358 preconditioner_env%low_rank_base = ot_low_rank_base_overlap
359 preconditioner_env%lattice_fft_active = .false.
360 preconditioner_env%in_use = ot_precond_fermi_low_rank
361 preconditioner_env%solver = ot_precond_solver_default
362
363 DEALLOCATE (common_evals, active_evals)
364 CALL cp_cfm_release(matrix_vectors)
365 CALL cp_cfm_release(matrix_chc)
366 CALL cp_cfm_release(matrix_hc0)
367
368 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 SUBROUTINE make_complex_full_all_covariant( &
384 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 POINTER :: chc_local
394 INTEGER :: group, i, j, k, n, ncol_local, ngroup, &
395 nrow_local, output_column
396 INTEGER, ALLOCATABLE, DIMENSION(:) :: group_representative, occupation_group
397 INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
398 REAL(kind=dp) :: occupation_tolerance, rayleigh_bound, &
399 shift_stride
400 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 CALL cp_cfm_get_info(matrix_c0, nrow_global=n, ncol_global=k)
405 cpassert(k > 0 .AND. k <= n)
406
407 CALL cp_cfm_create(matrix_hc0, matrix_c0%matrix_struct, &
408 name='complex covariant FULL_ALL HC')
409 CALL cp_cfm_create(matrix_chc, matrix_c0%matrix_struct, nrow=k, ncol=k, &
410 name='complex covariant FULL_ALL CHC')
411 CALL cp_cfm_create(matrix_occ_vectors, matrix_chc%matrix_struct, &
412 name='complex covariant FULL_ALL occupied rotation')
413 CALL cp_cfm_create(matrix_c_canonical, matrix_c0%matrix_struct, &
414 name='complex covariant FULL_ALL canonical C')
415 ALLOCATE (occupied_evals(k))
416
417 CALL cp_cfm_gemm('N', 'N', n, k, n, z_one, matrix_h, matrix_c0, z_zero, matrix_hc0)
418 CALL cp_cfm_gemm('C', 'N', k, k, n, z_one, matrix_c0, matrix_hc0, z_zero, matrix_chc)
419
420 IF (PRESENT(occupation_signature)) THEN
421 cpassert(SIZE(occupation_signature, 1) == k)
422 cpassert(SIZE(occupation_signature, 2) > 0)
423 ALLOCATE (group_representative(k), occupation_group(k))
424 occupation_tolerance = 64.0_dp*epsilon(1.0_dp)* &
425 max(1.0_dp, maxval(abs(occupation_signature)))
426 ngroup = 0
427 DO i = 1, k
428 occupation_group(i) = 0
429 DO group = 1, ngroup
430 IF (maxval(abs(occupation_signature(i, :) - &
431 occupation_signature(group_representative(group), :))) <= &
432 occupation_tolerance) THEN
433 occupation_group(i) = group
434 EXIT
435 END IF
436 END DO
437 IF (occupation_group(i) == 0) THEN
438 ngroup = ngroup + 1
439 group_representative(ngroup) = i
440 occupation_group(i) = ngroup
441 END IF
442 END DO
443 DEALLOCATE (group_representative)
444
445 IF (ngroup > 1) THEN
446 rayleigh_bound = cp_cfm_norm(matrix_chc, 'I')
447 shift_stride = 4.0_dp*max(1.0_dp, rayleigh_bound)
448 ALLOCATE (group_shift(ngroup))
449 DO group = 1, ngroup
450 group_shift(group) = shift_stride* &
451 (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 local_data=chc_local)
456 DO j = 1, ncol_local
457 DO i = 1, nrow_local
458 IF (occupation_group(row_indices(i)) /= occupation_group(col_indices(j))) THEN
459 chc_local(i, j) = z_zero
460 ELSE IF (row_indices(i) == col_indices(j)) THEN
461 group = occupation_group(row_indices(i))
462 chc_local(i, j) = chc_local(i, j) + &
463 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 CALL cp_cfm_heevd(matrix_chc, matrix_occ_vectors, occupied_evals)
471 IF (PRESENT(occupation_signature)) THEN
472 IF (ngroup > 1) THEN
473 output_column = 0
474 DO group = 1, ngroup
475 DO i = 1, count(occupation_group == group)
476 output_column = output_column + 1
477 occupied_evals(output_column) = occupied_evals(output_column) - group_shift(group)
478 END DO
479 END DO
480 cpassert(output_column == k)
481 DEALLOCATE (group_shift)
482 END IF
483 DEALLOCATE (occupation_group)
484 END IF
485 CALL cp_cfm_gemm('N', 'N', n, k, k, z_one, matrix_c0, matrix_occ_vectors, &
486 z_zero, matrix_c_canonical)
487
488 CALL make_complex_full_all(preconditioner_env, matrix_c_canonical, matrix_h, matrix_s, &
489 occupied_evals, energy_gap)
490 IF (ASSOCIATED(preconditioner_env%occ_rotation_complex)) THEN
491 CALL cp_cfm_release(preconditioner_env%occ_rotation_complex)
492 DEALLOCATE (preconditioner_env%occ_rotation_complex)
493 END IF
494 ALLOCATE (preconditioner_env%occ_rotation_complex)
495 CALL cp_cfm_create(preconditioner_env%occ_rotation_complex, matrix_occ_vectors%matrix_struct, &
496 name='complex covariant FULL_ALL occupied rotation')
497 CALL cp_cfm_to_cfm(matrix_occ_vectors, preconditioner_env%occ_rotation_complex)
498 preconditioner_env%in_use = ot_precond_full_all_covariant
499
500 DEALLOCATE (occupied_evals)
501 CALL cp_cfm_release(matrix_c_canonical)
502 CALL cp_cfm_release(matrix_occ_vectors)
503 CALL cp_cfm_release(matrix_chc)
504 CALL cp_cfm_release(matrix_hc0)
505
506 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 SUBROUTINE make_preconditioner_matrix(preconditioner_env, matrix_h, matrix_s, matrix_t, mo_coeff, &
521 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 precon_type = preconditioner_env%in_use
535 SELECT CASE (precon_type)
536 CASE (ot_precond_full_single)
537 IF (my_solver_type /= ot_precond_solver_default) THEN
538 cpabort("Only PRECOND_SOLVER DEFAULT for the moment")
539 END IF
540 IF (PRESENT(matrix_s)) THEN
541 CALL make_full_single(preconditioner_env, preconditioner_env%fm, &
542 matrix_h, matrix_s, energy_homo, energy_gap)
543 ELSE
544 CALL make_full_single_ortho(preconditioner_env, preconditioner_env%fm, &
545 matrix_h, energy_homo, energy_gap)
546 END IF
547
548 CASE (ot_precond_s_inverse)
549 IF (my_solver_type == ot_precond_solver_default) my_solver_type = ot_precond_solver_inv_chol
550 IF (.NOT. PRESENT(matrix_s)) THEN
551 cpabort("Type for S=1 not implemented")
552 END IF
553 CALL make_full_s_inverse(preconditioner_env, matrix_s)
554
555 CASE (ot_precond_full_kinetic)
556 IF (my_solver_type == ot_precond_solver_default) my_solver_type = ot_precond_solver_inv_chol
557 IF (.NOT. (PRESENT(matrix_s) .AND. PRESENT(matrix_t))) THEN
558 cpabort("Type for S=1 not implemented")
559 END IF
560 CALL make_full_kinetic(preconditioner_env, matrix_t, matrix_s, energy_gap)
561 CASE (ot_precond_full_single_inverse)
562 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 matrix_s=matrix_s)
565 CASE (ot_precond_full_all, ot_precond_full_all_covariant)
566 IF (my_solver_type /= ot_precond_solver_default) THEN
567 cpabort("Only PRECOND_SOLVER DEFAULT for the moment")
568 END IF
569 IF (PRESENT(matrix_s)) THEN
570 CALL make_full_all(preconditioner_env, mo_coeff, matrix_h, matrix_s, &
571 eigenvalues_ot, energy_gap)
572 ELSE
573 CALL make_full_all_ortho(preconditioner_env, mo_coeff, matrix_h, &
574 eigenvalues_ot, energy_gap)
575 END IF
576
577 CASE (ot_precond_fermi_low_rank)
578 IF (my_solver_type /= ot_precond_solver_default) THEN
579 cpabort("Only PRECOND_SOLVER DEFAULT for the moment")
580 END IF
581 IF (PRESENT(matrix_s)) THEN
582 CALL make_full_all(preconditioner_env, mo_coeff, matrix_h, matrix_s, &
583 eigenvalues_ot, energy_gap, common_reference=energy_homo)
584 ELSE
585 CALL make_full_all_ortho(preconditioner_env, mo_coeff, matrix_h, &
586 eigenvalues_ot, energy_gap, common_reference=energy_homo)
587 END IF
588
589 CASE DEFAULT
590 cpabort("Type not implemented")
591 END SELECT
592
593 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 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 CALL timeset(routinen, handle)
609
610 cpassert(ASSOCIATED(matrix_s))
611
612 IF (.NOT. ASSOCIATED(preconditioner_env%sparse_matrix)) THEN
613 ALLOCATE (preconditioner_env%sparse_matrix)
614 END IF
615 CALL dbcsr_copy(preconditioner_env%sparse_matrix, matrix_s, name="full_kinetic")
616
617 CALL timestop(handle)
618
619 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 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 CALL timeset(routinen, handle)
641
642 cpassert(ASSOCIATED(matrix_t))
643 cpassert(ASSOCIATED(matrix_s))
644
645 IF (.NOT. ASSOCIATED(preconditioner_env%sparse_matrix)) THEN
646 ALLOCATE (preconditioner_env%sparse_matrix)
647 END IF
648 CALL dbcsr_copy(preconditioner_env%sparse_matrix, matrix_t, name="full_kinetic")
649
650 shift = max(0.0_dp, energy_gap)
651
652 CALL dbcsr_add(preconditioner_env%sparse_matrix, matrix_s, &
653 alpha_scalar=1.0_dp, beta_scalar=shift)
654
655 CALL timestop(handle)
656
657 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 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 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 CALL timeset(routinen, handle)
683
684 NULLIFY (fm_struct_tmp, evals)
685
686 IF (ASSOCIATED(fm)) THEN
687 CALL cp_fm_release(fm)
688 DEALLOCATE (fm)
689 NULLIFY (fm)
690 END IF
691 CALL dbcsr_get_info(matrix_h, nfullrows_total=n)
692 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 para_env=preconditioner_env%para_env)
697 ALLOCATE (fm)
698 CALL cp_fm_create(fm, fm_struct_tmp, name="preconditioner")
699 CALL cp_fm_create(fm_h, fm_struct_tmp, name="fm_h")
700 CALL cp_fm_create(fm_s, fm_struct_tmp, name="fm_s")
701 CALL cp_fm_struct_release(fm_struct_tmp)
702
703 CALL copy_dbcsr_to_fm(matrix_h, fm_h)
704 CALL copy_dbcsr_to_fm(matrix_s, fm_s)
705 CALL cp_fm_cholesky_decompose(fm_s)
706
707 SELECT CASE (preconditioner_env%cholesky_use)
708 CASE (cholesky_inverse)
709! if cho inverse
710 CALL cp_fm_triangular_invert(fm_s)
711 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 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 invert_tr=.false., uplo_tr="U", n_rows=n, n_cols=n, alpha=1.0_dp)
717 CASE (cholesky_reduce)
718 CALL cp_fm_cholesky_reduce(fm_h, fm_s)
719 CASE DEFAULT
720 cpabort("cholesky type not implemented")
721 END SELECT
722
723 CALL choose_eigv_solver(fm_h, fm, evals)
724
725 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 invert_tr=.false., uplo_tr="U", n_rows=n, n_cols=n, alpha=1.0_dp)
729 DO i = 1, n
730 evals(i) = 1.0_dp/max(evals(i) - energy_homo, energy_gap)
731 END DO
732 CALL cp_fm_to_fm(fm, fm_h)
733 CASE (cholesky_reduce)
734 CALL cp_fm_cholesky_restore(fm, n, fm_s, fm_h, "SOLVE")
735 DO i = 1, n
736 evals(i) = 1.0_dp/max(evals(i) - energy_homo, energy_gap)
737 END DO
738 CALL cp_fm_to_fm(fm_h, fm)
739 END SELECT
740
741 CALL cp_fm_column_scale(fm, evals)
742 CALL parallel_gemm('N', 'T', n, n, n, 1.0_dp, fm, fm_h, 0.0_dp, fm_s)
743 CALL cp_fm_to_fm(fm_s, fm)
744
745 DEALLOCATE (evals)
746 CALL cp_fm_release(fm_h)
747 CALL cp_fm_release(fm_s)
748
749 CALL timestop(handle)
750
751 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 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 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 CALL timeset(routinen, handle)
776 NULLIFY (fm_struct_tmp, evals)
777
778 IF (ASSOCIATED(fm)) THEN
779 CALL cp_fm_release(fm)
780 DEALLOCATE (fm)
781 NULLIFY (fm)
782 END IF
783 CALL dbcsr_get_info(matrix_h, nfullrows_total=n)
784 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 para_env=preconditioner_env%para_env)
789 ALLOCATE (fm)
790 CALL cp_fm_create(fm, fm_struct_tmp, name="preconditioner")
791 CALL cp_fm_create(fm_h, fm_struct_tmp, name="fm_h")
792 CALL cp_fm_create(fm_s, fm_struct_tmp, name="fm_s")
793 CALL cp_fm_struct_release(fm_struct_tmp)
794
795 CALL copy_dbcsr_to_fm(matrix_h, fm_h)
796
797 CALL choose_eigv_solver(fm_h, fm, evals)
798 DO i = 1, n
799 evals(i) = 1.0_dp/max(evals(i) - energy_homo, energy_gap)
800 END DO
801 CALL cp_fm_to_fm(fm, fm_h)
802 CALL cp_fm_column_scale(fm, evals)
803 CALL parallel_gemm('N', 'T', n, n, n, 1.0_dp, fm, fm_h, 0.0_dp, fm_s)
804 CALL cp_fm_to_fm(fm_s, fm)
805
806 DEALLOCATE (evals)
807 CALL cp_fm_release(fm_h)
808 CALL cp_fm_release(fm_s)
809
810 CALL timestop(handle)
811
812 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 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 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 CALL timeset(routinen, handle)
857
858 IF (ASSOCIATED(preconditioner_env%fm)) THEN
859 CALL cp_fm_release(preconditioner_env%fm)
860 DEALLOCATE (preconditioner_env%fm)
861 NULLIFY (preconditioner_env%fm)
862 END IF
863 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 para_env=preconditioner_env%para_env)
867 ALLOCATE (preconditioner_env%fm)
868 CALL cp_fm_create(preconditioner_env%fm, fm_struct_tmp, name="preconditioner_env%fm")
869 CALL cp_fm_create(ortho, fm_struct_tmp, name="ortho")
870 CALL cp_fm_create(matrix_tmp, fm_struct_tmp, name="matrix_tmp")
871 CALL cp_fm_struct_release(fm_struct_tmp)
872 ALLOCATE (preconditioner_env%full_evals(n))
873 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 CALL copy_dbcsr_to_fm(matrix_s, ortho)
878 CALL cp_fm_cholesky_decompose(ortho)
879! if cho inverse
880 IF (preconditioner_env%cholesky_use == cholesky_inverse) THEN
881 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 CALL cp_fm_create(matrix_sc0, matrix_c0%matrix_struct, name="sc0")
890 CALL cp_dbcsr_sm_fm_multiply(matrix_s, matrix_c0, matrix_sc0, k)
891 CALL cp_fm_create(matrix_hc0, matrix_c0%matrix_struct, name="hc0")
892 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 CALL cp_fm_create(matrix_shc0, matrix_c0%matrix_struct, name="shc0")
896
897 SELECT CASE (preconditioner_env%cholesky_use)
898 CASE (cholesky_inverse)
899! if cho inverse
900 CALL cp_fm_to_fm(matrix_hc0, matrix_shc0)
901 CALL cp_fm_triangular_multiply(ortho, matrix_shc0, side="L", transpose_tr=.true., &
902 invert_tr=.false., uplo_tr="U", n_rows=n, n_cols=k, alpha=1.0_dp)
903 CASE (cholesky_reduce)
904 CALL cp_fm_cholesky_restore(matrix_hc0, k, ortho, matrix_shc0, "SOLVE", transa="T")
905 CASE DEFAULT
906 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 para_env=preconditioner_env%para_env)
911 CALL cp_fm_create(matrix_s1, fm_struct_tmp, name="matrix_s1")
912 CALL cp_fm_struct_release(fm_struct_tmp)
913 ! since we only use diagonal elements this is a bit of a waste
914 CALL parallel_gemm('T', 'N', k, k, n, 1.0_dp, matrix_shc0, matrix_shc0, 0.0_dp, matrix_s1)
915 ALLOCATE (diag(k))
916 CALL cp_fm_get_diag(matrix_s1, diag)
917 IF (PRESENT(common_reference)) THEN
918 ! ||(1-CC^T S)S^-1 H C||_F is invariant under C -> C U.
919 error_estimate = sqrt(max(0.0_dp, sum(diag) - sum(c0_evals**2)))
920 ELSE
921 error_estimate = maxval(sqrt(abs(diag - c0_evals**2)))
922 END IF
923 DEALLOCATE (diag)
924 CALL cp_fm_release(matrix_s1)
925 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 preconditioner_env%energy_gap = max(energy_gap, error_estimate*fudge_factor)
931 CALL copy_dbcsr_to_fm(matrix_h, matrix_tmp)
932 matrix_pre => preconditioner_env%fm
933 CALL cp_fm_uplo_to_full(matrix_tmp, matrix_pre)
934 ! tmp = H ( 1 - PS )
935 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 para_env=preconditioner_env%para_env)
940 CALL cp_fm_create(matrix_left, fm_struct_tmp, name="matrix_left")
941 CALL cp_fm_struct_release(fm_struct_tmp)
942 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 CALL parallel_gemm('N', 'N', n, n, k, -1.0_dp, matrix_sc0, matrix_left, 1.0_dp, matrix_tmp)
945 CALL cp_fm_release(matrix_left)
946
947 ALLOCATE (shifted_evals(k))
948 lambda = lambda_base + error_estimate
949 IF (PRESENT(common_reference)) THEN
950 shifted_evals = common_reference - lambda
951 ELSE
952 shifted_evals = c0_evals - lambda
953 END IF
954 CALL cp_fm_to_fm(matrix_sc0, matrix_hc0)
955 CALL cp_fm_column_scale(matrix_hc0, shifted_evals)
956 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 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 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 invert_tr=.false., uplo_tr="U", n_rows=n, n_cols=n, alpha=1.0_dp)
965 CASE (cholesky_reduce)
966 CALL cp_fm_cholesky_reduce(matrix_tmp, ortho)
967 END SELECT
968 CALL choose_eigv_solver(matrix_tmp, matrix_pre, preconditioner_env%full_evals)
969 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 invert_tr=.false., uplo_tr="U", n_rows=n, n_cols=n, alpha=1.0_dp)
973 CALL cp_fm_to_fm(matrix_pre, matrix_tmp)
974 CASE (cholesky_reduce)
975 CALL cp_fm_cholesky_restore(matrix_pre, n, ortho, matrix_tmp, "SOLVE")
976 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 IF (PRESENT(common_reference)) THEN
998 preconditioner_env%occ_evals = common_reference
999 preconditioner_env%full_evals(1:k) = common_reference
1000 ELSE
1001 preconditioner_env%occ_evals = c0_evals
1002 ! This choice makes FULL_ALL constant when applied to sc0 (see apply_full_all).
1003 preconditioner_env%full_evals(1:k) = c0_evals
1004 END IF
1005 CALL cp_fm_to_fm(matrix_c0, matrix_pre, k, 1, 1)
1006
1007 CALL cp_fm_release(matrix_sc0)
1008 CALL cp_fm_release(matrix_hc0)
1009 CALL cp_fm_release(ortho)
1010 CALL cp_fm_release(matrix_tmp)
1011 DEALLOCATE (shifted_evals)
1012 CALL timestop(handle)
1013
1014 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 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 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 CALL timeset(routinen, handle)
1048
1049 IF (ASSOCIATED(preconditioner_env%fm)) THEN
1050 CALL cp_fm_release(preconditioner_env%fm)
1051 DEALLOCATE (preconditioner_env%fm)
1052 NULLIFY (preconditioner_env%fm)
1053 END IF
1054 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 para_env=preconditioner_env%para_env)
1058 ALLOCATE (preconditioner_env%fm)
1059 CALL cp_fm_create(preconditioner_env%fm, fm_struct_tmp, name="preconditioner_env%fm")
1060 CALL cp_fm_create(matrix_tmp, fm_struct_tmp, name="matrix_tmp")
1061 CALL cp_fm_struct_release(fm_struct_tmp)
1062 ALLOCATE (preconditioner_env%full_evals(n))
1063 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 CALL cp_fm_create(matrix_sc0, matrix_c0%matrix_struct, name="sc0")
1072 CALL cp_fm_to_fm(matrix_c0, matrix_sc0)
1073 CALL cp_fm_create(matrix_hc0, matrix_c0%matrix_struct, name="hc0")
1074 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 para_env=preconditioner_env%para_env)
1080 CALL cp_fm_create(matrix_s1, fm_struct_tmp, name="matrix_s1")
1081 CALL cp_fm_struct_release(fm_struct_tmp)
1082 ! since we only use diagonal elements this is a bit of a waste
1083 CALL parallel_gemm('T', 'N', k, k, n, 1.0_dp, matrix_hc0, matrix_hc0, 0.0_dp, matrix_s1)
1084 ALLOCATE (diag(k))
1085 CALL cp_fm_get_diag(matrix_s1, diag)
1086 IF (PRESENT(common_reference)) THEN
1087 ! ||(1-CC^T)H C||_F is invariant under C -> C U.
1088 error_estimate = sqrt(max(0.0_dp, sum(diag) - sum(c0_evals**2)))
1089 ELSE
1090 error_estimate = maxval(sqrt(abs(diag - c0_evals**2)))
1091 END IF
1092 DEALLOCATE (diag)
1093 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 preconditioner_env%energy_gap = max(energy_gap, error_estimate*fudge_factor)
1099
1100 matrix_pre => preconditioner_env%fm
1101 CALL copy_dbcsr_to_fm(matrix_h, matrix_tmp)
1102 CALL cp_fm_uplo_to_full(matrix_tmp, matrix_pre)
1103 ! tmp = H ( 1 - PS )
1104 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 para_env=preconditioner_env%para_env)
1109 CALL cp_fm_create(matrix_left, fm_struct_tmp, name="matrix_left")
1110 CALL cp_fm_struct_release(fm_struct_tmp)
1111 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 CALL parallel_gemm('N', 'N', n, n, k, -1.0_dp, matrix_sc0, matrix_left, 1.0_dp, matrix_tmp)
1114 CALL cp_fm_release(matrix_left)
1115
1116 ALLOCATE (shifted_evals(k))
1117 lambda = lambda_base + error_estimate
1118 IF (PRESENT(common_reference)) THEN
1119 shifted_evals = common_reference - lambda
1120 ELSE
1121 shifted_evals = c0_evals - lambda
1122 END IF
1123 CALL cp_fm_to_fm(matrix_sc0, matrix_hc0)
1124 CALL cp_fm_column_scale(matrix_hc0, shifted_evals)
1125 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 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 IF (PRESENT(common_reference)) THEN
1151 preconditioner_env%occ_evals = common_reference
1152 preconditioner_env%full_evals(1:k) = common_reference
1153 ELSE
1154 preconditioner_env%occ_evals = c0_evals
1155 ! This choice makes FULL_ALL constant when applied to sc0 (see apply_full_all).
1156 preconditioner_env%full_evals(1:k) = c0_evals
1157 END IF
1158 CALL cp_fm_to_fm(matrix_c0, matrix_pre, k, 1, 1)
1159
1160 CALL cp_fm_release(matrix_sc0)
1161 CALL cp_fm_release(matrix_hc0)
1162 CALL cp_fm_release(matrix_tmp)
1163 DEALLOCATE (shifted_evals)
1164
1165 CALL timestop(handle)
1166
1167 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 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 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrices
1199 TYPE(dbcsr_type), TARGET :: dbcsr_cthc, dbcsr_hc, dbcsr_sc, mo_dbcsr
1200
1201 CALL timeset(routinen, handle)
1202
1203 ! Allocate all working matrices needed
1204 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 CALL cp_fm_to_dbcsr_row_template(mo_dbcsr, matrix_c0, matrix_h)
1208 CALL dbcsr_create(dbcsr_sc, template=mo_dbcsr)
1209 CALL dbcsr_create(dbcsr_hc, template=mo_dbcsr)
1210 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 IF (.NOT. ASSOCIATED(preconditioner_env%sparse_matrix)) THEN
1214 ALLOCATE (preconditioner_env%sparse_matrix)
1215 END IF
1216
1217 ! Put the first term of the preconditioner (H) into the output matrix
1218 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 IF (PRESENT(matrix_s)) THEN
1223 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_s, mo_dbcsr, 0.0_dp, dbcsr_sc)
1224 ELSE
1225 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 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_h, mo_dbcsr, 0.0_dp, dbcsr_hc)
1231 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 ALLOCATE (matrices(1))
1235 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 nval_request=1, nrestarts=8, generalized_ev=.false., iram=.false.)
1238 IF (ASSOCIATED(preconditioner_env%max_ev_vector)) THEN
1239 CALL set_arnoldi_initial_vector(arnoldi_env, preconditioner_env%max_ev_vector)
1240 END IF
1241 CALL arnoldi_ev(matrices, arnoldi_env)
1242 max_ev = real(get_selected_ritz_val(arnoldi_env, 1), dp)
1243
1244 ! save the ev as guess for the next time
1245 IF (.NOT. ASSOCIATED(preconditioner_env%max_ev_vector)) ALLOCATE (preconditioner_env%max_ev_vector)
1246 CALL get_selected_ritz_vector(arnoldi_env, 1, matrices(1)%matrix, preconditioner_env%max_ev_vector)
1247 CALL deallocate_arnoldi_env(arnoldi_env)
1248 DEALLOCATE (matrices)
1249
1250 ! Lets shift the occupied states a bit further up, -1.0 because we gonna subtract it from H
1251 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 CALL dbcsr_multiply("N", "N", 2.0_dp, dbcsr_sc, dbcsr_cthc, 0.0_dp, dbcsr_hc)
1254 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 IF (PRESENT(matrix_s)) THEN
1259 ALLOCATE (matrices(2))
1260 matrices(1)%matrix => preconditioner_env%sparse_matrix
1261 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 nval_request=1, nrestarts=21, generalized_ev=.true., iram=.false.)
1264 ELSE
1265 ALLOCATE (matrices(1))
1266 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 nval_request=1, nrestarts=8, generalized_ev=.false., iram=.false.)
1269 END IF
1270 IF (ASSOCIATED(preconditioner_env%min_ev_vector)) THEN
1271 CALL set_arnoldi_initial_vector(arnoldi_env, preconditioner_env%min_ev_vector)
1272 END IF
1273
1274 ! compute the LUMO energy
1275 CALL arnoldi_ev(matrices, arnoldi_env)
1276 min_ev = real(get_selected_ritz_val(arnoldi_env, 1), dp)
1277
1278 ! save the lumo vector for restarting in the next step
1279 IF (.NOT. ASSOCIATED(preconditioner_env%min_ev_vector)) ALLOCATE (preconditioner_env%min_ev_vector)
1280 CALL get_selected_ritz_vector(arnoldi_env, 1, matrices(1)%matrix, preconditioner_env%min_ev_vector)
1281 CALL deallocate_arnoldi_env(arnoldi_env)
1282 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 pre_shift = max(1.5_dp*(min_ev - max_ev), energy_gap)
1288 IF (min_ev < pre_shift) THEN
1289 pre_shift = pre_shift - min_ev
1290 ELSE
1291 pre_shift = 0.0_dp
1292 END IF
1293 IF (PRESENT(matrix_s)) THEN
1294 CALL dbcsr_add(preconditioner_env%sparse_matrix, matrix_s, 1.0_dp, pre_shift)
1295 ELSE
1296 CALL dbcsr_add_on_diag(preconditioner_env%sparse_matrix, pre_shift)
1297 END IF
1298
1299 CALL dbcsr_release(mo_dbcsr)
1300 CALL dbcsr_release(dbcsr_hc)
1301 CALL dbcsr_release(dbcsr_sc)
1302 CALL dbcsr_release(dbcsr_cthc)
1303
1304 CALL timestop(handle)
1305
1306 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 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 CALL timeset(routinen, handle)
1333
1334 CALL cp_cfm_get_info(matrix_c0, nrow_global=n, ncol_global=k)
1335 cpassert(n > 0)
1336 cpassert(k > 0 .AND. k <= n)
1337 cpassert(energy_gap > 0.0_dp)
1338
1339 CALL cp_cfm_create(matrix_hc0, matrix_c0%matrix_struct, name='complex FULL_SINGLE HC')
1340 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 name='complex FULL_SINGLE CHC')
1343 CALL cp_cfm_create(matrix_sc_chc, matrix_c0%matrix_struct, &
1344 name='complex FULL_SINGLE SC CHC')
1345 CALL cp_cfm_create(matrix_operator, matrix_h%matrix_struct, &
1346 name='complex FULL_SINGLE operator')
1347 CALL cp_cfm_create(matrix_diag_a, matrix_h%matrix_struct, &
1348 name='complex FULL_SINGLE diagonalization A')
1349 CALL cp_cfm_create(matrix_diag_b, matrix_s%matrix_struct, &
1350 name='complex FULL_SINGLE diagonalization B')
1351 CALL cp_cfm_create(matrix_diag_evec, matrix_h%matrix_struct, nrow=n, ncol=1, &
1352 name='complex FULL_SINGLE eigenvectors')
1353 CALL cp_cfm_create(matrix_diag_work, matrix_h%matrix_struct, &
1354 name='complex FULL_SINGLE diagonalization work')
1355 CALL cp_cfm_create(matrix_occ_diag_a, matrix_c0%matrix_struct, nrow=k, ncol=k, &
1356 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 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 name='complex FULL_SINGLE occupied eigenvectors')
1361 CALL cp_cfm_create(matrix_occ_diag_work, matrix_c0%matrix_struct, nrow=k, ncol=k, &
1362 name='complex FULL_SINGLE occupied diagonalization work')
1363
1364 CALL cp_cfm_gemm('N', 'N', n, k, n, z_one, matrix_h, matrix_c0, z_zero, matrix_hc0)
1365 CALL cp_cfm_gemm('N', 'N', n, k, n, z_one, matrix_s, matrix_c0, z_zero, matrix_sc0)
1366 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 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 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 z_one, matrix_operator)
1374 CALL cp_cfm_gemm('N', 'C', n, n, k, z_one, matrix_sc0, matrix_sc0, &
1375 z_one, matrix_operator)
1376
1377 ALLOCATE (occupied_evals(k), operator_evals(1))
1378 CALL cp_cfm_to_cfm(matrix_chc, matrix_occ_diag_a)
1379 DO i = 1, k
1380 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 occupied_evals, matrix_occ_diag_work)
1384 max_ev = maxval(occupied_evals)
1385 CALL cp_cfm_to_cfm(matrix_operator, matrix_diag_a)
1386 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 operator_evals, matrix_diag_work, lowest_subset=.true.)
1389 min_ev = operator_evals(1)
1390 target_edge = max(1.5_dp*(min_ev - max_ev), energy_gap)
1391 pre_shift = max(0.0_dp, target_edge - min_ev)
1392 IF (pre_shift > 0.0_dp) THEN
1393 CALL cp_cfm_scale_and_add(z_one, matrix_operator, &
1394 cmplx(pre_shift, 0.0_dp, kind=dp), matrix_s)
1395 END IF
1396
1397 CALL store_complex_inverse(preconditioner_env, matrix_operator, &
1398 ot_precond_full_single_inverse, energy_gap)
1399
1400 DEALLOCATE (operator_evals, occupied_evals)
1401 CALL cp_cfm_release(matrix_occ_diag_work)
1402 CALL cp_cfm_release(matrix_occ_diag_evec)
1403 CALL cp_cfm_release(matrix_occ_diag_b)
1404 CALL cp_cfm_release(matrix_occ_diag_a)
1405 CALL cp_cfm_release(matrix_diag_work)
1406 CALL cp_cfm_release(matrix_diag_evec)
1407 CALL cp_cfm_release(matrix_diag_b)
1408 CALL cp_cfm_release(matrix_diag_a)
1409 CALL cp_cfm_release(matrix_operator)
1410 CALL cp_cfm_release(matrix_sc_chc)
1411 CALL cp_cfm_release(matrix_chc)
1412 CALL cp_cfm_release(matrix_sc0)
1413 CALL cp_cfm_release(matrix_hc0)
1414
1415 CALL timestop(handle)
1416
1417 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 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 CALL store_complex_inverse(preconditioner_env, matrix_s, ot_precond_s_inverse, 0.0_dp)
1430
1431 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 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 shift = max(0.0_dp, energy_gap)
1450 CALL cp_cfm_create(matrix_operator, matrix_t%matrix_struct, &
1451 name='complex FULL_KINETIC operator')
1452 CALL cp_cfm_to_cfm(matrix_t, matrix_operator)
1453 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 ot_precond_full_kinetic, energy_gap)
1456 CALL cp_cfm_release(matrix_operator)
1457
1458 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 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 IF (ASSOCIATED(preconditioner_env%complex_fm)) THEN
1478 CALL cp_cfm_release(preconditioner_env%complex_fm)
1479 DEALLOCATE (preconditioner_env%complex_fm)
1480 END IF
1481 ALLOCATE (preconditioner_env%complex_fm)
1482 CALL cp_cfm_create(preconditioner_env%complex_fm, matrix_operator%matrix_struct, &
1483 name='complex inverse preconditioner')
1484 CALL cp_cfm_to_cfm(matrix_operator, preconditioner_env%complex_fm)
1485 CALL cp_cfm_cholesky_decompose(preconditioner_env%complex_fm, info_out=info)
1486 cpassert(info == 0)
1487 CALL cp_cfm_cholesky_invert(preconditioner_env%complex_fm, info_out=info)
1488 cpassert(info == 0)
1489 CALL cp_cfm_uplo_to_full(preconditioner_env%complex_fm)
1490
1491 preconditioner_env%energy_gap = energy_gap
1492 preconditioner_env%in_use = preconditioner_kind
1493 preconditioner_env%solver = ot_precond_solver_default
1494
1495 END SUBROUTINE store_complex_inverse
1496END MODULE preconditioner_makes
cp_cfm_types::cp_cfm_to_cfm
Definition cp_cfm_types.F:57
cp_dbcsr_api::dbcsr_create
Definition cp_dbcsr_api.F:194
cp_fm_types::cp_fm_release
Definition cp_fm_types.F:88
cp_fm_types::cp_fm_to_fm
Definition cp_fm_types.F:83
parallel_gemm_api::parallel_gemm
Definition parallel_gemm_api.F:39
arnoldi_api
arnoldi iteration using dbcsr
Definition arnoldi_api.F:16
arnoldi_api::arnoldi_ev
subroutine, public arnoldi_ev(matrix, arnoldi_env)
Driver routine for different arnoldi eigenvalue methods the selection which one is to be taken is mad...
Definition arnoldi_api.F:69
cp_cfm_basic_linalg
Basic linear algebra operations for complex full matrices.
Definition cp_cfm_basic_linalg.F:16
cp_cfm_basic_linalg::cp_cfm_scale_and_add
subroutine, public cp_cfm_scale_and_add(alpha, matrix_a, beta, matrix_b)
Scale and add two BLACS matrices (a = alpha*a + beta*b).
Definition cp_cfm_basic_linalg.F:243
cp_cfm_basic_linalg::cp_cfm_norm
real(kind=dp) function, public cp_cfm_norm(matrix, mode)
Norm of matrix using (p)zlange.
Definition cp_cfm_basic_linalg.F:1185
cp_cfm_basic_linalg::cp_cfm_gemm
subroutine, public cp_cfm_gemm(transa, transb, m, n, k, alpha, matrix_a, matrix_b, beta, matrix_c, a_first_col, a_first_row, b_first_col, b_first_row, c_first_col, c_first_row)
Performs one of the matrix-matrix operations: matrix_c = alpha * op1( matrix_a ) * op2( matrix_b ) + ...
Definition cp_cfm_basic_linalg.F:562
cp_cfm_basic_linalg::cp_cfm_triangular_multiply
subroutine, public cp_cfm_triangular_multiply(triangular_matrix, matrix_b, side, transa_tr, invert_tr, uplo_tr, unit_diag_tr, n_rows, n_cols, alpha)
Multiplies in place by a triangular matrix: matrix_b = alpha op(triangular_matrix) matrix_b or (if si...
Definition cp_cfm_basic_linalg.F:978
cp_cfm_basic_linalg::cp_cfm_uplo_to_full
subroutine, public cp_cfm_uplo_to_full(matrix, workspace, uplo)
...
Definition cp_cfm_basic_linalg.F:1350
cp_cfm_basic_linalg::cp_cfm_column_scale
subroutine, public cp_cfm_column_scale(matrix_a, scaling)
Scales columns of the full matrix by corresponding factors.
Definition cp_cfm_basic_linalg.F:632
cp_cfm_cholesky
various cholesky decomposition related routines
Definition cp_cfm_cholesky.F:13
cp_cfm_cholesky::cp_cfm_cholesky_decompose
subroutine, public cp_cfm_cholesky_decompose(matrix, n, info_out)
Used to replace a symmetric positive definite matrix M with its Cholesky decomposition U: M = U^T * U...
Definition cp_cfm_cholesky.F:53
cp_cfm_cholesky::cp_cfm_cholesky_invert
subroutine, public cp_cfm_cholesky_invert(matrix, n, info_out)
Used to replace Cholesky decomposition by the inverse.
Definition cp_cfm_cholesky.F:124
cp_cfm_diag
used for collecting diagonalization schemes available for cp_cfm_type
Definition cp_cfm_diag.F:14
cp_cfm_diag::cp_cfm_heevd
subroutine, public cp_cfm_heevd(matrix, eigenvectors, eigenvalues)
Perform a diagonalisation of a complex matrix.
Definition cp_cfm_diag.F:92
cp_cfm_diag::cp_cfm_geeig
subroutine, public cp_cfm_geeig(amatrix, bmatrix, eigenvectors, eigenvalues, work, lowest_subset)
General Eigenvalue Problem AX = BXE Single option version: Cholesky decomposition of B.
Definition cp_cfm_diag.F:340
cp_cfm_types
Represents a complex full matrix distributed on many processors.
Definition cp_cfm_types.F:12
cp_cfm_types::cp_cfm_release
subroutine, public cp_cfm_release(matrix)
Releases a full matrix.
Definition cp_cfm_types.F:183
cp_cfm_types::cp_cfm_set_element
subroutine, public cp_cfm_set_element(matrix, irow_global, icol_global, alpha)
Set the matrix element (irow_global,icol_global) of the full matrix to alpha.
Definition cp_cfm_types.F:303
cp_cfm_types::cp_cfm_create
subroutine, public cp_cfm_create(matrix, matrix_struct, name, nrow, ncol, set_zero)
Creates a new full matrix with the given structure.
Definition cp_cfm_types.F:126
cp_cfm_types::cp_cfm_get_info
subroutine, public cp_cfm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, matrix_struct, para_env)
Returns information about a full matrix.
Definition cp_cfm_types.F:632
cp_dbcsr_api
Definition cp_dbcsr_api.F:8
cp_dbcsr_api::dbcsr_copy
subroutine, public dbcsr_copy(matrix_b, matrix_a, name, keep_sparsity, keep_imaginary)
...
Definition cp_dbcsr_api.F:371
cp_dbcsr_api::dbcsr_multiply
subroutine, public dbcsr_multiply(transa, transb, alpha, matrix_a, matrix_b, beta, matrix_c, first_row, last_row, first_column, last_column, first_k, last_k, retain_sparsity, filter_eps, flop)
...
Definition cp_dbcsr_api.F:1096
cp_dbcsr_api::dbcsr_get_info
subroutine, public dbcsr_get_info(matrix, nblkrows_total, nblkcols_total, nfullrows_total, nfullcols_total, nblkrows_local, nblkcols_local, nfullrows_local, nfullcols_local, my_prow, my_pcol, local_rows, local_cols, proc_row_dist, proc_col_dist, row_blk_size, col_blk_size, row_blk_offset, col_blk_offset, distribution, name, matrix_type, group)
...
Definition cp_dbcsr_api.F:808
cp_dbcsr_api::dbcsr_release
subroutine, public dbcsr_release(matrix)
...
Definition cp_dbcsr_api.F:1142
cp_dbcsr_api::dbcsr_add
subroutine, public dbcsr_add(matrix_a, matrix_b, alpha_scalar, beta_scalar)
...
Definition cp_dbcsr_api.F:254
cp_dbcsr_contrib
Definition cp_dbcsr_contrib.F:8
cp_dbcsr_contrib::dbcsr_add_on_diag
subroutine, public dbcsr_add_on_diag(matrix, alpha)
Adds the given scalar to the diagonal of the matrix. Reserves any missing diagonal blocks.
Definition cp_dbcsr_contrib.F:331
cp_dbcsr_operations
DBCSR operations in CP2K.
Definition cp_dbcsr_operations.F:18
cp_dbcsr_operations::cp_dbcsr_sm_fm_multiply
subroutine, public cp_dbcsr_sm_fm_multiply(matrix, fm_in, fm_out, ncol, alpha, beta)
multiply a dbcsr with a fm matrix
Definition cp_dbcsr_operations.F:799
cp_dbcsr_operations::copy_dbcsr_to_fm
subroutine, public copy_dbcsr_to_fm(matrix, fm, plan)
Copy a DBCSR matrix to a BLACS matrix.
Definition cp_dbcsr_operations.F:232
cp_dbcsr_operations::cp_fm_to_dbcsr_row_template
subroutine, public cp_fm_to_dbcsr_row_template(matrix, fm_in, template)
Utility function to copy a specially shaped fm to dbcsr_matrix The result matrix will be the matrix i...
Definition cp_dbcsr_operations.F:1109
cp_dbcsr_operations::cp_dbcsr_m_by_n_from_template
subroutine, public cp_dbcsr_m_by_n_from_template(matrix, template, m, n, sym)
Utility function to create an arbitrary shaped dbcsr matrix with the same processor grid as the templ...
Definition cp_dbcsr_operations.F:1143
cp_fm_basic_linalg
Basic linear algebra operations for full matrices.
Definition cp_fm_basic_linalg.F:14
cp_fm_basic_linalg::cp_fm_column_scale
subroutine, public cp_fm_column_scale(matrixa, scaling)
scales column i of matrix a with scaling(i)
Definition cp_fm_basic_linalg.F:1690
cp_fm_basic_linalg::cp_fm_uplo_to_full
subroutine, public cp_fm_uplo_to_full(matrix, work, uplo)
given a triangular matrix according to uplo, computes the corresponding full matrix
Definition cp_fm_basic_linalg.F:1584
cp_fm_basic_linalg::cp_fm_triangular_invert
subroutine, public cp_fm_triangular_invert(matrix_a, uplo_tr)
inverts a triangular matrix
Definition cp_fm_basic_linalg.F:1968
cp_fm_basic_linalg::cp_fm_triangular_multiply
subroutine, public cp_fm_triangular_multiply(triangular_matrix, matrix_b, side, transpose_tr, invert_tr, uplo_tr, unit_diag_tr, n_rows, n_cols, alpha)
multiplies in place by a triangular matrix: matrix_b = alpha op(triangular_matrix) matrix_b or (if si...
Definition cp_fm_basic_linalg.F:1413
cp_fm_cholesky
various cholesky decomposition related routines
Definition cp_fm_cholesky.F:14
cp_fm_cholesky::cp_fm_cholesky_restore
subroutine, public cp_fm_cholesky_restore(fm_matrix, neig, fm_matrixb, fm_matrixout, op, pos, transa)
apply Cholesky decomposition op can be "SOLVE" (out = U^-1 * in) or "MULTIPLY" (out = U * in) pos can...
Definition cp_fm_cholesky.F:313
cp_fm_cholesky::cp_fm_cholesky_decompose
subroutine, public cp_fm_cholesky_decompose(matrix, n, info_out)
used to replace a symmetric positive def. matrix M with its cholesky decomposition U: M = U^T * U,...
Definition cp_fm_cholesky.F:61
cp_fm_cholesky::cp_fm_cholesky_reduce
subroutine, public cp_fm_cholesky_reduce(matrix, matrixb, itype)
reduce a matrix pencil A,B to normal form B has to be cholesky decomposed with cp_fm_cholesky_decompo...
Definition cp_fm_cholesky.F:255
cp_fm_diag
used for collecting some of the diagonalization schemes available for cp_fm_type. cp_fm_power also mo...
Definition cp_fm_diag.F:17
cp_fm_diag::choose_eigv_solver
subroutine, public choose_eigv_solver(matrix, eigenvectors, eigenvalues, info)
Choose the Eigensolver depending on which library is available ELPA seems to be unstable for small sy...
Definition cp_fm_diag.F:262
cp_fm_struct
represent the structure of a full matrix
Definition cp_fm_struct.F:14
cp_fm_struct::cp_fm_struct_create
subroutine, public cp_fm_struct_create(fmstruct, para_env, context, nrow_global, ncol_global, nrow_block, ncol_block, descriptor, first_p_pos, local_leading_dimension, template_fmstruct, square_blocks, force_block)
allocates and initializes a full matrix structure
Definition cp_fm_struct.F:132
cp_fm_struct::cp_fm_struct_release
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
Definition cp_fm_struct.F:355
cp_fm_types
represent a full matrix distributed on many processors
Definition cp_fm_types.F:15
cp_fm_types::cp_fm_get_diag
subroutine, public cp_fm_get_diag(matrix, diag)
returns the diagonal elements of a fm
Definition cp_fm_types.F:564
cp_fm_types::cp_fm_get_info
subroutine, public cp_fm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, nrow_locals, ncol_locals, matrix_struct, para_env)
returns all kind of information about the full matrix
Definition cp_fm_types.F:1038
cp_fm_types::cp_fm_create
subroutine, public cp_fm_create(matrix, matrix_struct, name, nrow, ncol, set_zero)
creates a new full matrix with the given structure
Definition cp_fm_types.F:165
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:18
input_constants::ot_low_rank_base_overlap
integer, parameter, public ot_low_rank_base_overlap
Definition input_constants.F:548
input_constants::ot_precond_full_all_covariant
integer, parameter, public ot_precond_full_all_covariant
Definition input_constants.F:539
input_constants::ot_precond_full_kinetic
integer, parameter, public ot_precond_full_kinetic
Definition input_constants.F:539
input_constants::cholesky_reduce
integer, parameter, public cholesky_reduce
Definition input_constants.F:206
input_constants::cholesky_inverse
integer, parameter, public cholesky_inverse
Definition input_constants.F:206
input_constants::ot_precond_solver_default
integer, parameter, public ot_precond_solver_default
Definition input_constants.F:553
input_constants::ot_precond_full_single
integer, parameter, public ot_precond_full_single
Definition input_constants.F:539
input_constants::ot_precond_solver_inv_chol
integer, parameter, public ot_precond_solver_inv_chol
Definition input_constants.F:553
input_constants::ot_precond_full_single_inverse
integer, parameter, public ot_precond_full_single_inverse
Definition input_constants.F:539
input_constants::ot_precond_fermi_low_rank
integer, parameter, public ot_precond_fermi_low_rank
Definition input_constants.F:539
input_constants::ot_precond_s_inverse
integer, parameter, public ot_precond_s_inverse
Definition input_constants.F:539
input_constants::ot_precond_full_all
integer, parameter, public ot_precond_full_all
Definition input_constants.F:539
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
low_rank_preconditioner_model
Small algebraic helpers for the rotationally covariant low-rank OT preconditioner.
Definition low_rank_preconditioner_model.F:11
low_rank_preconditioner_model::low_rank_select_rank
pure integer function, public low_rank_select_rank(eigenvalues, nocc, max_rank, degeneracy_tolerance)
Select a bounded spectral rank without cutting a degenerate boundary manifold.
Definition low_rank_preconditioner_model.F:50
mathconstants
Definition of mathematical constants and functions.
Definition mathconstants.F:16
mathconstants::z_one
complex(kind=dp), parameter, public z_one
Definition mathconstants.F:144
mathconstants::z_zero
complex(kind=dp), parameter, public z_zero
Definition mathconstants.F:144
parallel_gemm_api
basic linear algebra operations for full matrixes
Definition parallel_gemm_api.F:14
preconditioner_makes
computes preconditioners, and implements methods to apply them currently used in qs_ot
Definition preconditioner_makes.F:15
preconditioner_makes::make_complex_full_single_inverse
subroutine, public make_complex_full_single_inverse(preconditioner_env, matrix_c0, matrix_h, matrix_s, energy_gap)
Build a gauge-covariant FULL_SINGLE_INVERSE operator for a complex k-point channel.
Definition preconditioner_makes.F:1318
preconditioner_makes::make_complex_full_kinetic
subroutine, public make_complex_full_kinetic(preconditioner_env, matrix_t, matrix_s, energy_gap)
Build the inverse complex kinetic-plus-overlap preconditioner.
Definition preconditioner_makes.F:1441
preconditioner_makes::make_complex_full_s_inverse
subroutine, public make_complex_full_s_inverse(preconditioner_env, matrix_s)
Build the inverse complex overlap preconditioner.
Definition preconditioner_makes.F:1425
preconditioner_makes::make_complex_full_all_covariant
subroutine, public make_complex_full_all_covariant(preconditioner_env, matrix_c0, matrix_h, matrix_s, energy_gap, occupation_signature)
Build the gauge-covariant FULL_ALL inverse for a complex k-point channel.
Definition preconditioner_makes.F:385
preconditioner_makes::make_complex_full_single
subroutine, public make_complex_full_single(preconditioner_env, matrix_h, matrix_s, energy_homo, energy_gap)
Build the complex spectral FULL_SINGLE preconditioner.
Definition preconditioner_makes.F:236
preconditioner_makes::make_complex_fermi_low_rank
subroutine, public make_complex_fermi_low_rank(preconditioner_env, matrix_c0, matrix_h, matrix_s, energy_gap, max_rank, spectral_window, degeneracy_tolerance)
Build the bounded low-rank spectral inverse for a complex k-point channel.
Definition preconditioner_makes.F:320
preconditioner_makes::make_preconditioner_matrix
subroutine, public make_preconditioner_matrix(preconditioner_env, matrix_h, matrix_s, matrix_t, mo_coeff, energy_homo, eigenvalues_ot, energy_gap, my_solver_type)
...
Definition preconditioner_makes.F:523
preconditioner_makes::make_complex_full_all
subroutine, public make_complex_full_all(preconditioner_env, matrix_c0, matrix_h, matrix_s, c0_evals, energy_gap)
Build the state-selective FULL_ALL operator for a complex k-point channel. The occupied/reference sub...
Definition preconditioner_makes.F:104
preconditioner_types
types of preconditioners
Definition preconditioner_types.F:14
cp_cfm_types::cp_cfm_type
Represent a complex full matrix.
Definition cp_cfm_types.F:70
cp_dbcsr_api::dbcsr_p_type
Definition cp_dbcsr_api.F:172
cp_dbcsr_api::dbcsr_type
Definition cp_dbcsr_api.F:176
cp_fm_struct::cp_fm_struct_type
keeps the information about the structure of a full matrix
Definition cp_fm_struct.F:82
cp_fm_types::cp_fm_type
represent a full matrix
Definition cp_fm_types.F:114
preconditioner_types::preconditioner_type
Definition preconditioner_types.F:60