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 an eigen-space solver for the generalised symmetric eigenvalue problem
10 : !> for sparse matrices, needing only multiplications
11 : !> \author Joost VandeVondele (25.08.2002)
12 : ! **************************************************************************************************
13 : MODULE qs_ot_eigensolver
14 : USE cp_cfm_basic_linalg, ONLY: cp_cfm_gemm
15 : USE cp_cfm_diag, ONLY: cp_cfm_heevd
16 : USE cp_cfm_types, ONLY: cp_cfm_create,&
17 : cp_cfm_release,&
18 : cp_cfm_to_fm,&
19 : cp_cfm_type,&
20 : cp_fm_to_cfm
21 : USE cp_dbcsr_api, ONLY: &
22 : dbcsr_add, dbcsr_copy, dbcsr_init_p, dbcsr_multiply, dbcsr_p_type, dbcsr_release_p, &
23 : dbcsr_scale, dbcsr_set, dbcsr_type, dbcsr_type_no_symmetry
24 : USE cp_dbcsr_cholesky, ONLY: cp_dbcsr_cholesky_decompose,&
25 : cp_dbcsr_cholesky_invert
26 : USE cp_dbcsr_contrib, ONLY: dbcsr_dot
27 : USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
28 : copy_fm_to_dbcsr,&
29 : cp_dbcsr_m_by_n_from_row_template,&
30 : cp_dbcsr_m_by_n_from_template,&
31 : cp_dbcsr_sm_fm_multiply,&
32 : cp_fm_to_dbcsr_row_template,&
33 : dbcsr_copy_columns_hack
34 : USE cp_fm_types, ONLY: cp_fm_create,&
35 : cp_fm_get_info,&
36 : cp_fm_release,&
37 : cp_fm_type
38 : USE cp_log_handling, ONLY: cp_logger_get_default_io_unit
39 : USE input_constants, ONLY: ot_precond_fermi_low_rank,&
40 : ot_precond_full_all,&
41 : ot_precond_full_all_covariant,&
42 : ot_precond_full_kinetic,&
43 : ot_precond_full_single,&
44 : ot_precond_full_single_inverse,&
45 : ot_precond_none,&
46 : ot_precond_s_inverse
47 : USE kinds, ONLY: dp
48 : USE mathconstants, ONLY: z_one,&
49 : z_zero
50 : USE preconditioner, ONLY: make_preconditioner_complex_fermi_low_rank,&
51 : make_preconditioner_complex_full_all,&
52 : make_preconditioner_complex_full_all_covariant,&
53 : make_preconditioner_complex_full_kinetic,&
54 : make_preconditioner_complex_full_s_inverse,&
55 : make_preconditioner_complex_full_single,&
56 : make_preconditioner_complex_full_single_inverse
57 : USE preconditioner_types, ONLY: destroy_preconditioner,&
58 : init_preconditioner,&
59 : preconditioner_in_use,&
60 : preconditioner_type
61 : USE qs_mo_methods, ONLY: make_basis_sv
62 : USE qs_mo_types, ONLY: mo_set_type
63 : USE qs_ot, ONLY: qs_ot_get_orbitals,&
64 : qs_ot_get_orbitals_complex,&
65 : qs_ot_get_orbitals_ref,&
66 : qs_ot_get_orbitals_ref_complex,&
67 : qs_ot_get_p,&
68 : qs_ot_new_preconditioner,&
69 : qs_ot_prepare_complex_tangent_metric
70 : USE qs_ot_minimizer, ONLY: ot_mini
71 : USE qs_ot_types, ONLY: qs_ot_allocate,&
72 : qs_ot_allocate_complex_state,&
73 : qs_ot_destroy,&
74 : qs_ot_init,&
75 : qs_ot_settings_init,&
76 : qs_ot_settings_type,&
77 : qs_ot_type
78 : #include "./base/base_uses.f90"
79 :
80 : IMPLICIT NONE
81 : PRIVATE
82 :
83 : ! *** Global parameters ***
84 :
85 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_ot_eigensolver'
86 :
87 : ! *** Public subroutines ***
88 :
89 : PUBLIC :: ot_eigensolver, ot_eigensolver_complex
90 :
91 : CONTAINS
92 :
93 : ! on input c contains the initial guess (should not be zero !)
94 : ! on output c spans the subspace
95 : ! **************************************************************************************************
96 : !> \brief ...
97 : !> \param matrix_h ...
98 : !> \param matrix_s ...
99 : !> \param matrix_orthogonal_space_fm ...
100 : !> \param matrix_c_fm ...
101 : !> \param preconditioner ...
102 : !> \param eps_gradient ...
103 : !> \param iter_max ...
104 : !> \param size_ortho_space ...
105 : !> \param silent ...
106 : !> \param ot_settings ...
107 : ! **************************************************************************************************
108 1008 : SUBROUTINE ot_eigensolver(matrix_h, matrix_s, matrix_orthogonal_space_fm, &
109 : matrix_c_fm, preconditioner, eps_gradient, &
110 : iter_max, size_ortho_space, silent, ot_settings)
111 :
112 : TYPE(dbcsr_type), POINTER :: matrix_h, matrix_s
113 : TYPE(cp_fm_type), INTENT(IN), OPTIONAL :: matrix_orthogonal_space_fm
114 : TYPE(cp_fm_type), INTENT(INOUT) :: matrix_c_fm
115 : TYPE(preconditioner_type), OPTIONAL, POINTER :: preconditioner
116 : REAL(KIND=dp) :: eps_gradient
117 : INTEGER, INTENT(IN) :: iter_max
118 : INTEGER, INTENT(IN), OPTIONAL :: size_ortho_space
119 : LOGICAL, INTENT(IN), OPTIONAL :: silent
120 : TYPE(qs_ot_settings_type), INTENT(IN), OPTIONAL :: ot_settings
121 :
122 : CHARACTER(len=*), PARAMETER :: routineN = 'ot_eigensolver'
123 : INTEGER, PARAMETER :: max_iter_inner_loop = 40
124 : REAL(KIND=dp), PARAMETER :: rone = 1.0_dp, rzero = 0.0_dp
125 :
126 : INTEGER :: handle, ieigensolver, iter_total, k, n, &
127 : ortho_k, ortho_space_k, output_unit
128 : LOGICAL :: energy_only, my_silent, ortho, &
129 : ref_algorithm
130 : REAL(KIND=dp) :: delta, energy
131 504 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_hc
132 : TYPE(dbcsr_type), POINTER :: matrix_buf1_ortho, matrix_buf2_ortho, &
133 : matrix_c, matrix_orthogonal_space, &
134 : matrix_os_ortho, matrix_s_ortho
135 504 : TYPE(qs_ot_type), DIMENSION(:), POINTER :: qs_ot_env
136 :
137 504 : CALL timeset(routineN, handle)
138 :
139 504 : output_unit = cp_logger_get_default_io_unit()
140 :
141 504 : IF (PRESENT(silent)) THEN
142 254 : my_silent = silent
143 : ELSE
144 : my_silent = .FALSE.
145 : END IF
146 :
147 504 : NULLIFY (matrix_c) ! fm->dbcsr
148 :
149 504 : CALL cp_fm_get_info(matrix_c_fm, nrow_global=n, ncol_global=k) ! fm->dbcsr
150 504 : ALLOCATE (matrix_c)
151 504 : ref_algorithm = .FALSE.
152 504 : IF (PRESENT(ot_settings)) ref_algorithm = ot_settings%ot_algorithm == "REF"
153 250 : IF (ref_algorithm) THEN
154 24 : CALL dbcsr_init_p(matrix_c)
155 : CALL cp_dbcsr_m_by_n_from_row_template(matrix_c, template=matrix_h, n=k, &
156 24 : sym=dbcsr_type_no_symmetry)
157 24 : CALL copy_fm_to_dbcsr(matrix_c_fm, matrix_c)
158 : ELSE
159 480 : CALL cp_fm_to_dbcsr_row_template(matrix_c, fm_in=matrix_c_fm, template=matrix_h)
160 : END IF
161 :
162 504 : iter_total = 0
163 :
164 : outer_scf: DO
165 :
166 : NULLIFY (qs_ot_env)
167 :
168 672 : NULLIFY (matrix_s_ortho)
169 672 : NULLIFY (matrix_os_ortho)
170 672 : NULLIFY (matrix_buf1_ortho)
171 672 : NULLIFY (matrix_buf2_ortho)
172 672 : NULLIFY (matrix_orthogonal_space)
173 :
174 110880 : ALLOCATE (qs_ot_env(1))
175 1344 : ALLOCATE (matrix_hc(1))
176 672 : NULLIFY (matrix_hc(1)%matrix)
177 672 : CALL dbcsr_init_p(matrix_hc(1)%matrix)
178 :
179 672 : ortho = .FALSE.
180 672 : IF (PRESENT(matrix_orthogonal_space_fm)) ortho = .TRUE.
181 :
182 : ! decide settings
183 672 : IF (PRESENT(ot_settings)) THEN
184 264 : qs_ot_env(1)%settings = ot_settings
185 : ELSE
186 408 : CALL qs_ot_settings_init(qs_ot_env(1)%settings)
187 : ! overwrite defaults
188 408 : qs_ot_env(1)%settings%ds_min = 0.10_dp
189 : END IF
190 :
191 672 : IF (ortho) THEN
192 408 : ALLOCATE (matrix_orthogonal_space)
193 408 : CALL cp_fm_to_dbcsr_row_template(matrix_orthogonal_space, fm_in=matrix_orthogonal_space_fm, template=matrix_h)
194 408 : CALL cp_fm_get_info(matrix_orthogonal_space_fm, ncol_global=ortho_space_k)
195 :
196 408 : IF (PRESENT(size_ortho_space)) ortho_space_k = size_ortho_space
197 408 : ortho_k = ortho_space_k + k
198 : ELSE
199 264 : ortho_k = k
200 : END IF
201 :
202 : ! allocate
203 672 : CALL qs_ot_allocate(qs_ot_env(1), matrix_s, matrix_c_fm%matrix_struct, ortho_k=ortho_k)
204 672 : IF (ref_algorithm) THEN
205 30 : CALL dbcsr_copy(matrix_hc(1)%matrix, qs_ot_env(1)%matrix_x, 'matrix_hc')
206 : ELSE
207 642 : CALL dbcsr_copy(matrix_hc(1)%matrix, matrix_c, 'matrix_hc')
208 : END IF
209 :
210 672 : IF (ortho) THEN
211 : ! construct an initial guess that is orthogonal to matrix_orthogonal_space
212 :
213 408 : CALL dbcsr_init_p(matrix_s_ortho)
214 408 : CALL dbcsr_copy(matrix_s_ortho, matrix_orthogonal_space, name="matrix_s_ortho")
215 :
216 408 : CALL dbcsr_init_p(matrix_os_ortho)
217 : CALL cp_dbcsr_m_by_n_from_template(matrix_os_ortho, template=matrix_h, m=ortho_space_k, n=ortho_space_k, &
218 408 : sym=dbcsr_type_no_symmetry)
219 :
220 408 : CALL dbcsr_init_p(matrix_buf1_ortho)
221 : CALL cp_dbcsr_m_by_n_from_template(matrix_buf1_ortho, template=matrix_h, m=ortho_space_k, n=k, &
222 408 : sym=dbcsr_type_no_symmetry)
223 :
224 408 : CALL dbcsr_init_p(matrix_buf2_ortho)
225 : CALL cp_dbcsr_m_by_n_from_template(matrix_buf2_ortho, template=matrix_h, m=ortho_space_k, n=k, &
226 408 : sym=dbcsr_type_no_symmetry)
227 :
228 : CALL dbcsr_multiply('N', 'N', 1.0_dp, matrix_s, matrix_orthogonal_space, &
229 408 : 0.0_dp, matrix_s_ortho)
230 : CALL dbcsr_multiply('T', 'N', rone, matrix_s_ortho, matrix_s_ortho, &
231 408 : rzero, matrix_os_ortho)
232 :
233 : CALL cp_dbcsr_cholesky_decompose(matrix_os_ortho, &
234 408 : para_env=qs_ot_env(1)%para_env, blacs_env=qs_ot_env(1)%blacs_env)
235 : CALL cp_dbcsr_cholesky_invert(matrix_os_ortho, &
236 : para_env=qs_ot_env(1)%para_env, blacs_env=qs_ot_env(1)%blacs_env, &
237 408 : uplo_to_full=.TRUE.)
238 :
239 : CALL dbcsr_multiply('T', 'N', rone, matrix_s_ortho, matrix_c, &
240 408 : rzero, matrix_buf1_ortho)
241 : CALL dbcsr_multiply('N', 'N', rone, matrix_os_ortho, matrix_buf1_ortho, &
242 408 : rzero, matrix_buf2_ortho)
243 : CALL dbcsr_multiply('N', 'N', -rone, matrix_s_ortho, matrix_buf2_ortho, &
244 408 : rone, matrix_c)
245 :
246 : ! make matrix_c0 an orthogonal basis, matrix_c contains sc0
247 408 : CALL dbcsr_copy(qs_ot_env(1)%matrix_c0, matrix_c)
248 : CALL dbcsr_multiply('N', 'N', 1.0_dp, matrix_s, qs_ot_env(1)%matrix_c0, &
249 408 : 0.0_dp, matrix_c)
250 :
251 : CALL make_basis_sv(qs_ot_env(1)%matrix_c0, k, matrix_c, &
252 408 : qs_ot_env(1)%para_env, qs_ot_env(1)%blacs_env)
253 :
254 : ! copy sc0 and matrix_s_ortho in qs_ot_env(1)%matrix_sc0
255 : !CALL dbcsr_copy_columns(qs_ot_env(1)%matrix_sc0,matrix_s_ortho,ortho_space_k,1,1)
256 : CALL dbcsr_copy_columns_hack(qs_ot_env(1)%matrix_sc0, matrix_s_ortho, ortho_space_k, 1, 1, &
257 408 : para_env=qs_ot_env(1)%para_env, blacs_env=qs_ot_env(1)%blacs_env)
258 : !CALL dbcsr_copy_columns(qs_ot_env(1)%matrix_sc0,matrix_c,k,1,ortho_space_k+1)
259 : CALL dbcsr_copy_columns_hack(qs_ot_env(1)%matrix_sc0, matrix_c, k, 1, ortho_space_k + 1, &
260 408 : para_env=qs_ot_env(1)%para_env, blacs_env=qs_ot_env(1)%blacs_env)
261 :
262 408 : CALL dbcsr_release_p(matrix_buf1_ortho)
263 408 : CALL dbcsr_release_p(matrix_buf2_ortho)
264 408 : CALL dbcsr_release_p(matrix_os_ortho)
265 408 : CALL dbcsr_release_p(matrix_s_ortho)
266 :
267 : ELSE
268 :
269 : ! set c0,sc0
270 264 : CALL dbcsr_copy(qs_ot_env(1)%matrix_c0, matrix_c)
271 : CALL dbcsr_multiply('N', 'N', 1.0_dp, matrix_s, qs_ot_env(1)%matrix_c0, &
272 264 : 0.0_dp, qs_ot_env(1)%matrix_sc0)
273 :
274 : CALL make_basis_sv(qs_ot_env(1)%matrix_c0, k, qs_ot_env(1)%matrix_sc0, &
275 264 : qs_ot_env(1)%para_env, qs_ot_env(1)%blacs_env)
276 : END IF
277 :
278 : ! init
279 672 : CALL qs_ot_init(qs_ot_env(1))
280 672 : energy_only = qs_ot_env(1)%energy_only
281 :
282 642 : SELECT CASE (qs_ot_env(1)%settings%ot_algorithm)
283 : CASE ("TOD")
284 642 : CALL dbcsr_set(qs_ot_env(1)%matrix_x, 0.0_dp)
285 642 : CALL dbcsr_set(qs_ot_env(1)%matrix_sx, 0.0_dp)
286 642 : CALL qs_ot_get_p(qs_ot_env(1)%matrix_x, qs_ot_env(1)%matrix_sx, qs_ot_env(1))
287 642 : CALL qs_ot_get_orbitals(matrix_c, qs_ot_env(1)%matrix_x, qs_ot_env(1))
288 : CASE ("REF")
289 30 : CALL dbcsr_copy(qs_ot_env(1)%matrix_x, qs_ot_env(1)%matrix_c0)
290 30 : CALL dbcsr_copy(qs_ot_env(1)%matrix_sx, qs_ot_env(1)%matrix_sc0)
291 30 : CALL dbcsr_copy(matrix_c, qs_ot_env(1)%matrix_c0)
292 : CASE DEFAULT
293 672 : CPABORT("OT eigensolver supports ALGORITHM STRICT or IRAC")
294 : END SELECT
295 :
296 : ! if present preconditioner, use it
297 :
298 672 : IF (PRESENT(preconditioner)) THEN
299 672 : IF (ASSOCIATED(preconditioner)) THEN
300 444 : IF (preconditioner_in_use(preconditioner)) THEN
301 444 : CALL qs_ot_new_preconditioner(qs_ot_env(1), preconditioner)
302 : ELSE
303 : ! we should presumably make one
304 : END IF
305 : END IF
306 : END IF
307 :
308 : ! *** Eigensolver loop ***
309 : ieigensolver = 0
310 11850 : eigensolver_loop: DO
311 :
312 11850 : ieigensolver = ieigensolver + 1
313 11850 : iter_total = iter_total + 1
314 :
315 : ! the energy is cHc, the gradient is 2*H*c
316 : CALL dbcsr_multiply('N', 'N', 1.0_dp, matrix_h, matrix_c, &
317 11850 : 0.0_dp, matrix_hc(1)%matrix)
318 11850 : CALL dbcsr_dot(matrix_c, matrix_hc(1)%matrix, energy)
319 11850 : IF (.NOT. energy_only) THEN
320 6304 : CALL dbcsr_scale(matrix_hc(1)%matrix, 2.0_dp)
321 : END IF
322 :
323 11850 : qs_ot_env(1)%etotal = energy
324 11850 : CALL ot_mini(qs_ot_env, matrix_hc)
325 11850 : delta = qs_ot_env(1)%delta
326 11850 : energy_only = qs_ot_env(1)%energy_only
327 :
328 11306 : SELECT CASE (qs_ot_env(1)%settings%ot_algorithm)
329 : CASE ("TOD")
330 : CALL dbcsr_multiply('N', 'N', 1.0_dp, matrix_s, qs_ot_env(1)%matrix_x, &
331 11306 : 0.0_dp, qs_ot_env(1)%matrix_sx)
332 11306 : CALL qs_ot_get_p(qs_ot_env(1)%matrix_x, qs_ot_env(1)%matrix_sx, qs_ot_env(1))
333 11306 : CALL qs_ot_get_orbitals(matrix_c, qs_ot_env(1)%matrix_x, qs_ot_env(1))
334 : CASE ("REF")
335 : CALL qs_ot_get_orbitals_ref(matrix_c, matrix_s, qs_ot_env(1)%matrix_x, &
336 : qs_ot_env(1)%matrix_sx, qs_ot_env(1)%matrix_gx_old, &
337 11850 : qs_ot_env(1)%matrix_dx, qs_ot_env(1), qs_ot_env(1))
338 : END SELECT
339 :
340 : ! exit on convergence or if maximum of inner loop cycles is reached
341 11850 : IF (delta < eps_gradient .OR. ieigensolver >= max_iter_inner_loop) EXIT eigensolver_loop
342 : ! exit if total number of steps is reached, but not during a line search step
343 11850 : IF (iter_total >= iter_max .AND. qs_ot_env(1)%OT_METHOD_FULL /= "OT LS") EXIT eigensolver_loop
344 :
345 : END DO eigensolver_loop
346 :
347 672 : CALL qs_ot_destroy(qs_ot_env(1))
348 672 : DEALLOCATE (qs_ot_env)
349 672 : CALL dbcsr_release_p(matrix_hc(1)%matrix)
350 672 : DEALLOCATE (matrix_hc)
351 672 : CALL dbcsr_release_p(matrix_orthogonal_space)
352 :
353 672 : IF (delta < eps_gradient) THEN
354 354 : IF ((output_unit > 0) .AND. .NOT. my_silent) THEN
355 : WRITE (UNIT=output_unit, FMT="(T2,A,I0,A)") &
356 131 : "OT| Eigensolver reached convergence in ", iter_total, " iterations"
357 : END IF
358 : EXIT outer_scf
359 : END IF
360 318 : IF (iter_total >= iter_max) THEN
361 150 : IF (output_unit > 0) THEN
362 75 : IF (my_silent) THEN
363 75 : WRITE (output_unit, "(A,T60,E20.10)") " WARNING OT eigensolver did not converge: current gradient", delta
364 : ELSE
365 0 : WRITE (output_unit, *) "WARNING : did not converge in ot_eigensolver"
366 0 : WRITE (output_unit, *) "number of iterations ", iter_total, " exceeded maximum"
367 0 : WRITE (output_unit, *) "current gradient / target gradient", delta, " / ", eps_gradient
368 : END IF
369 : END IF
370 : EXIT outer_scf
371 : END IF
372 :
373 : END DO outer_scf
374 :
375 504 : CALL copy_dbcsr_to_fm(matrix_c, matrix_c_fm) ! fm->dbcsr
376 504 : CALL dbcsr_release_p(matrix_c) ! fm->dbcsr
377 :
378 504 : CALL timestop(handle)
379 :
380 504 : END SUBROUTINE ot_eigensolver
381 :
382 : ! **************************************************************************************************
383 : !> \brief solve a fixed complex Hermitian generalized eigenproblem by OT
384 : !> \param matrix_h real part of H(k)
385 : !> \param matrix_h_im imaginary part of H(k)
386 : !> \param matrix_s real part of S(k)
387 : !> \param matrix_s_im imaginary part of S(k)
388 : !> \param matrix_c_fm real part of the orbital coefficients
389 : !> \param matrix_c_fm_im imaginary part of the orbital coefficients
390 : !> \param preconditioner optional complex k-point preconditioner
391 : !> \param eps_gradient requested OT gradient accuracy
392 : !> \param iter_max maximum number of OT iterations
393 : !> \param eigenvalues ...
394 : !> \param silent suppress successful convergence output
395 : !> \param ot_settings OT algorithm and minimizer settings
396 : !> \param matrix_t ...
397 : !> \param matrix_t_im ...
398 : !> \param mo_set ...
399 : ! **************************************************************************************************
400 474 : SUBROUTINE ot_eigensolver_complex(matrix_h, matrix_h_im, matrix_s, matrix_s_im, &
401 : matrix_c_fm, matrix_c_fm_im, preconditioner, &
402 158 : eps_gradient, iter_max, eigenvalues, silent, ot_settings, &
403 : matrix_t, matrix_t_im, mo_set)
404 :
405 : TYPE(dbcsr_type), POINTER :: matrix_h, matrix_h_im, matrix_s, &
406 : matrix_s_im
407 : TYPE(cp_fm_type), INTENT(INOUT) :: matrix_c_fm, matrix_c_fm_im
408 : TYPE(preconditioner_type), OPTIONAL, POINTER :: preconditioner
409 : REAL(KIND=dp), INTENT(IN) :: eps_gradient
410 : INTEGER, INTENT(IN) :: iter_max
411 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT), OPTIONAL :: eigenvalues
412 : LOGICAL, INTENT(IN), OPTIONAL :: silent
413 : TYPE(qs_ot_settings_type), INTENT(IN), OPTIONAL :: ot_settings
414 : TYPE(dbcsr_type), OPTIONAL, POINTER :: matrix_t, matrix_t_im
415 : TYPE(mo_set_type), INTENT(IN), OPTIONAL :: mo_set
416 :
417 : CHARACTER(len=*), PARAMETER :: routineN = 'ot_eigensolver_complex'
418 : INTEGER, PARAMETER :: max_iter_inner_loop = 40
419 :
420 : INTEGER :: handle, ieigensolver, iter_total, k, n, &
421 : output_unit
422 : LOGICAL :: energy_only, my_silent, &
423 : preconditioner_rejected
424 : REAL(KIND=dp) :: delta, energy_im, energy_re
425 158 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_hc, matrix_hc_im
426 : TYPE(dbcsr_type), POINTER :: matrix_c, matrix_c_im, matrix_tmp
427 : TYPE(preconditioner_type), POINTER :: local_preconditioner
428 158 : TYPE(qs_ot_type), DIMENSION(:), POINTER :: qs_ot_env
429 :
430 158 : CALL timeset(routineN, handle)
431 158 : output_unit = cp_logger_get_default_io_unit()
432 158 : my_silent = .FALSE.
433 158 : IF (PRESENT(silent)) my_silent = silent
434 158 : NULLIFY (local_preconditioner)
435 :
436 158 : CALL cp_fm_get_info(matrix_c_fm, nrow_global=n, ncol_global=k)
437 158 : CPASSERT(n > 0 .AND. k > 0)
438 158 : ALLOCATE (matrix_c, matrix_c_im, matrix_tmp)
439 158 : CALL dbcsr_init_p(matrix_c)
440 : CALL cp_dbcsr_m_by_n_from_row_template(matrix_c, template=matrix_h, n=k, &
441 158 : sym=dbcsr_type_no_symmetry)
442 158 : CALL dbcsr_init_p(matrix_c_im)
443 : CALL cp_dbcsr_m_by_n_from_row_template(matrix_c_im, template=matrix_h, n=k, &
444 158 : sym=dbcsr_type_no_symmetry)
445 158 : CALL dbcsr_init_p(matrix_tmp)
446 : CALL cp_dbcsr_m_by_n_from_row_template(matrix_tmp, template=matrix_h, n=k, &
447 158 : sym=dbcsr_type_no_symmetry)
448 158 : CALL copy_fm_to_dbcsr(matrix_c_fm, matrix_c)
449 158 : CALL copy_fm_to_dbcsr(matrix_c_fm_im, matrix_c_im)
450 :
451 158 : iter_total = 0
452 : outer_scf: DO
453 27889 : ALLOCATE (qs_ot_env(1), matrix_hc(1), matrix_hc_im(1))
454 167 : NULLIFY (matrix_hc(1)%matrix, matrix_hc_im(1)%matrix)
455 :
456 167 : IF (PRESENT(ot_settings)) THEN
457 167 : qs_ot_env(1)%settings = ot_settings
458 : ELSE
459 0 : CALL qs_ot_settings_init(qs_ot_env(1)%settings)
460 0 : qs_ot_env(1)%settings%ds_min = 0.10_dp
461 : END IF
462 :
463 167 : CALL qs_ot_allocate(qs_ot_env(1), matrix_s, matrix_c_fm%matrix_struct)
464 167 : CALL qs_ot_allocate_complex_state(qs_ot_env(1), matrix_s)
465 167 : CALL dbcsr_init_p(matrix_hc(1)%matrix)
466 167 : CALL dbcsr_copy(matrix_hc(1)%matrix, qs_ot_env(1)%matrix_x, 'matrix_hc')
467 167 : CALL dbcsr_init_p(matrix_hc_im(1)%matrix)
468 167 : CALL dbcsr_copy(matrix_hc_im(1)%matrix, qs_ot_env(1)%matrix_x_im, 'matrix_hc_im')
469 :
470 167 : CALL dbcsr_copy(qs_ot_env(1)%matrix_c0, matrix_c)
471 167 : CALL dbcsr_copy(qs_ot_env(1)%matrix_c0_im, matrix_c_im)
472 : CALL dbcsr_multiply('N', 'N', 1.0_dp, matrix_s, matrix_c, &
473 167 : 0.0_dp, qs_ot_env(1)%matrix_sc0)
474 : CALL dbcsr_multiply('N', 'N', 1.0_dp, matrix_s_im, matrix_c_im, &
475 167 : 0.0_dp, matrix_tmp)
476 : CALL dbcsr_add(qs_ot_env(1)%matrix_sc0, matrix_tmp, &
477 167 : alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
478 : CALL dbcsr_multiply('N', 'N', 1.0_dp, matrix_s, matrix_c_im, &
479 167 : 0.0_dp, qs_ot_env(1)%matrix_sc0_im)
480 : CALL dbcsr_multiply('N', 'N', 1.0_dp, matrix_s_im, matrix_c, &
481 167 : 0.0_dp, matrix_tmp)
482 : CALL dbcsr_add(qs_ot_env(1)%matrix_sc0_im, matrix_tmp, &
483 167 : alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
484 :
485 167 : CALL qs_ot_init(qs_ot_env(1))
486 167 : energy_only = qs_ot_env(1)%energy_only
487 100 : SELECT CASE (qs_ot_env(1)%settings%ot_algorithm)
488 : CASE ('TOD')
489 100 : CALL dbcsr_set(qs_ot_env(1)%matrix_x, 0.0_dp)
490 100 : CALL dbcsr_set(qs_ot_env(1)%matrix_x_im, 0.0_dp)
491 100 : CALL dbcsr_set(qs_ot_env(1)%matrix_sx, 0.0_dp)
492 100 : CALL dbcsr_set(qs_ot_env(1)%matrix_sx_im, 0.0_dp)
493 : CALL qs_ot_get_orbitals_complex(matrix_c, matrix_c_im, matrix_s, matrix_s_im, &
494 100 : qs_ot_env(1))
495 : CASE ('REF')
496 67 : CALL dbcsr_copy(qs_ot_env(1)%matrix_x, qs_ot_env(1)%matrix_c0)
497 67 : CALL dbcsr_copy(qs_ot_env(1)%matrix_x_im, qs_ot_env(1)%matrix_c0_im)
498 67 : CALL dbcsr_copy(qs_ot_env(1)%matrix_sx, qs_ot_env(1)%matrix_sc0)
499 67 : CALL dbcsr_copy(qs_ot_env(1)%matrix_sx_im, qs_ot_env(1)%matrix_sc0_im)
500 : CALL qs_ot_get_orbitals_ref_complex(matrix_c, matrix_c_im, matrix_s, matrix_s_im, &
501 67 : qs_ot_env(1), qs_ot_env(1))
502 : CASE DEFAULT
503 167 : CPABORT('Complex OT eigensolver supports ALGORITHM STRICT or IRAC')
504 : END SELECT
505 :
506 167 : IF (.NOT. PRESENT(preconditioner) .AND. .NOT. ASSOCIATED(local_preconditioner) .AND. &
507 : qs_ot_env(1)%settings%preconditioner_type /= ot_precond_none) THEN
508 632 : ALLOCATE (local_preconditioner)
509 158 : CALL init_preconditioner(local_preconditioner, qs_ot_env(1)%para_env, qs_ot_env(1)%blacs_env)
510 198 : SELECT CASE (qs_ot_env(1)%settings%preconditioner_type)
511 : CASE (ot_precond_fermi_low_rank)
512 : CALL make_preconditioner_complex_fermi_low_rank( &
513 : local_preconditioner, qs_ot_env(1)%matrix_c0, qs_ot_env(1)%matrix_c0_im, &
514 : matrix_h, matrix_h_im, matrix_s, matrix_s_im, &
515 : qs_ot_env(1)%settings%energy_gap, &
516 : qs_ot_env(1)%settings%fermi_low_rank_max_rank, &
517 40 : qs_ot_env(1)%settings%precond_solver_type)
518 : CASE (ot_precond_full_all, ot_precond_full_all_covariant)
519 56 : IF (qs_ot_env(1)%settings%ot_algorithm == 'TOD') THEN
520 : ! State-selective FULL_ALL models can destabilize the finite STRICT fixed-H chart
521 : ! even when the tangent metric is positive. Retain a well-defined eigensolver with
522 : ! their S^-1 base; IRAC can use the complete state-dependent inverse below.
523 : CALL make_preconditioner_complex_full_s_inverse( &
524 : local_preconditioner, matrix_s, matrix_s_im, &
525 16 : qs_ot_env(1)%settings%precond_solver_type)
526 40 : ELSE IF (qs_ot_env(1)%settings%preconditioner_type == ot_precond_full_all) THEN
527 0 : IF (.NOT. PRESENT(mo_set)) THEN
528 0 : CPABORT('Complex FULL_ALL eigensolver preconditioning requires MO energy labels')
529 : END IF
530 : CALL make_preconditioner_complex_full_all( &
531 : local_preconditioner, qs_ot_env(1)%matrix_c0, qs_ot_env(1)%matrix_c0_im, &
532 : matrix_h, matrix_h_im, matrix_s, matrix_s_im, mo_set, &
533 0 : qs_ot_env(1)%settings%energy_gap, qs_ot_env(1)%settings%precond_solver_type)
534 : ELSE
535 : CALL make_preconditioner_complex_full_all_covariant( &
536 : local_preconditioner, qs_ot_env(1)%matrix_c0, qs_ot_env(1)%matrix_c0_im, &
537 : matrix_h, matrix_h_im, matrix_s, matrix_s_im, &
538 40 : qs_ot_env(1)%settings%energy_gap, qs_ot_env(1)%settings%precond_solver_type)
539 : END IF
540 : CASE (ot_precond_full_single)
541 22 : IF (.NOT. PRESENT(mo_set)) THEN
542 0 : CPABORT('Complex FULL_SINGLE eigensolver preconditioning requires MO energy labels')
543 : END IF
544 : CALL make_preconditioner_complex_full_single( &
545 : local_preconditioner, matrix_h, matrix_h_im, matrix_s, matrix_s_im, mo_set, &
546 22 : qs_ot_env(1)%settings%energy_gap, qs_ot_env(1)%settings%precond_solver_type)
547 : CASE (ot_precond_full_single_inverse)
548 : CALL make_preconditioner_complex_full_single_inverse( &
549 : local_preconditioner, qs_ot_env(1)%matrix_c0, qs_ot_env(1)%matrix_c0_im, &
550 : matrix_h, matrix_h_im, matrix_s, matrix_s_im, &
551 40 : qs_ot_env(1)%settings%energy_gap, qs_ot_env(1)%settings%precond_solver_type)
552 : CASE (ot_precond_full_kinetic)
553 0 : IF (.NOT. PRESENT(matrix_t) .OR. .NOT. PRESENT(matrix_t_im)) THEN
554 0 : CPABORT('Complex FULL_KINETIC eigensolver preconditioning requires T(k)')
555 : END IF
556 0 : IF (.NOT. ASSOCIATED(matrix_t) .OR. .NOT. ASSOCIATED(matrix_t_im)) THEN
557 0 : CPABORT('Complex FULL_KINETIC eigensolver preconditioning requires T(k)')
558 : END IF
559 : CALL make_preconditioner_complex_full_kinetic( &
560 : local_preconditioner, matrix_t, matrix_t_im, matrix_s, matrix_s_im, &
561 0 : qs_ot_env(1)%settings%energy_gap, qs_ot_env(1)%settings%precond_solver_type)
562 : CASE (ot_precond_s_inverse)
563 : CALL make_preconditioner_complex_full_s_inverse( &
564 : local_preconditioner, matrix_s, matrix_s_im, &
565 0 : qs_ot_env(1)%settings%precond_solver_type)
566 : CASE DEFAULT
567 158 : CPABORT('Unsupported complex K-point OT eigensolver preconditioner')
568 : END SELECT
569 : END IF
570 :
571 167 : IF (PRESENT(preconditioner)) THEN
572 0 : IF (ASSOCIATED(preconditioner)) THEN
573 0 : IF (preconditioner_in_use(preconditioner)) THEN
574 0 : CALL qs_ot_new_preconditioner(qs_ot_env(1), preconditioner)
575 : END IF
576 : END IF
577 167 : ELSE IF (ASSOCIATED(local_preconditioner)) THEN
578 167 : CALL qs_ot_new_preconditioner(qs_ot_env(1), local_preconditioner)
579 : END IF
580 :
581 167 : IF (qs_ot_env(1)%settings%ot_algorithm == 'TOD') THEN
582 100 : CALL qs_ot_prepare_complex_tangent_metric(qs_ot_env(1), preconditioner_rejected)
583 100 : IF (preconditioner_rejected .AND. ASSOCIATED(local_preconditioner)) THEN
584 0 : CALL destroy_preconditioner(local_preconditioner)
585 : CALL init_preconditioner(local_preconditioner, &
586 0 : qs_ot_env(1)%para_env, qs_ot_env(1)%blacs_env)
587 : CALL make_preconditioner_complex_full_s_inverse( &
588 : local_preconditioner, matrix_s, matrix_s_im, &
589 0 : qs_ot_env(1)%settings%precond_solver_type)
590 0 : CALL qs_ot_new_preconditioner(qs_ot_env(1), local_preconditioner)
591 0 : CALL qs_ot_prepare_complex_tangent_metric(qs_ot_env(1), preconditioner_rejected)
592 0 : CPASSERT(.NOT. preconditioner_rejected)
593 : END IF
594 : END IF
595 :
596 : ieigensolver = 0
597 1802 : eigensolver_loop: DO
598 901 : ieigensolver = ieigensolver + 1
599 901 : iter_total = iter_total + 1
600 :
601 : ! H*C = (Hre*Cre-Him*Cim) + i*(Hre*Cim+Him*Cre).
602 : CALL dbcsr_multiply('N', 'N', 1.0_dp, matrix_h, matrix_c, &
603 901 : 0.0_dp, matrix_hc(1)%matrix)
604 : CALL dbcsr_multiply('N', 'N', 1.0_dp, matrix_h_im, matrix_c_im, &
605 901 : 0.0_dp, matrix_tmp)
606 : CALL dbcsr_add(matrix_hc(1)%matrix, matrix_tmp, &
607 901 : alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
608 : CALL dbcsr_multiply('N', 'N', 1.0_dp, matrix_h, matrix_c_im, &
609 901 : 0.0_dp, matrix_hc_im(1)%matrix)
610 : CALL dbcsr_multiply('N', 'N', 1.0_dp, matrix_h_im, matrix_c, &
611 901 : 0.0_dp, matrix_tmp)
612 : CALL dbcsr_add(matrix_hc_im(1)%matrix, matrix_tmp, &
613 901 : alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
614 :
615 901 : CALL dbcsr_dot(matrix_c, matrix_hc(1)%matrix, energy_re)
616 901 : CALL dbcsr_dot(matrix_c_im, matrix_hc_im(1)%matrix, energy_im)
617 901 : qs_ot_env(1)%etotal = energy_re + energy_im
618 901 : IF (.NOT. energy_only) THEN
619 901 : CALL dbcsr_scale(matrix_hc(1)%matrix, 2.0_dp)
620 901 : CALL dbcsr_scale(matrix_hc_im(1)%matrix, 2.0_dp)
621 : END IF
622 901 : CALL ot_mini(qs_ot_env, matrix_hc, matrix_hc_im=matrix_hc_im)
623 901 : delta = qs_ot_env(1)%delta
624 901 : energy_only = qs_ot_env(1)%energy_only
625 :
626 426 : SELECT CASE (qs_ot_env(1)%settings%ot_algorithm)
627 : CASE ('TOD')
628 : CALL qs_ot_get_orbitals_complex(matrix_c, matrix_c_im, matrix_s, matrix_s_im, &
629 426 : qs_ot_env(1))
630 : CASE ('REF')
631 : CALL qs_ot_get_orbitals_ref_complex(matrix_c, matrix_c_im, matrix_s, matrix_s_im, &
632 901 : qs_ot_env(1), qs_ot_env(1))
633 : END SELECT
634 :
635 901 : IF (delta < eps_gradient .OR. ieigensolver >= max_iter_inner_loop) EXIT eigensolver_loop
636 901 : IF (iter_total >= iter_max .AND. qs_ot_env(1)%OT_METHOD_FULL /= 'OT LS') EXIT eigensolver_loop
637 : END DO eigensolver_loop
638 :
639 167 : CALL qs_ot_destroy(qs_ot_env(1))
640 167 : DEALLOCATE (qs_ot_env)
641 167 : CALL dbcsr_release_p(matrix_hc(1)%matrix)
642 167 : CALL dbcsr_release_p(matrix_hc_im(1)%matrix)
643 167 : DEALLOCATE (matrix_hc, matrix_hc_im)
644 :
645 167 : IF (delta < eps_gradient) THEN
646 156 : IF (output_unit > 0 .AND. .NOT. my_silent) THEN
647 : WRITE (UNIT=output_unit, FMT='(T2,A,I0,A)') &
648 0 : 'OT| Complex eigensolver reached convergence in ', iter_total, ' iterations'
649 : END IF
650 : EXIT outer_scf
651 : END IF
652 11 : IF (iter_total >= iter_max) THEN
653 2 : IF (output_unit > 0) THEN
654 : WRITE (output_unit, '(A,T60,E20.10)') &
655 0 : ' WARNING complex OT eigensolver did not converge: current gradient', delta
656 : END IF
657 : EXIT outer_scf
658 : END IF
659 : END DO outer_scf
660 :
661 158 : CALL copy_dbcsr_to_fm(matrix_c, matrix_c_fm)
662 158 : CALL copy_dbcsr_to_fm(matrix_c_im, matrix_c_fm_im)
663 158 : IF (PRESENT(eigenvalues)) THEN
664 : CALL canonicalize_complex_subspace(matrix_h, matrix_h_im, matrix_c_fm, &
665 158 : matrix_c_fm_im, eigenvalues)
666 : END IF
667 158 : CALL dbcsr_release_p(matrix_c)
668 158 : CALL dbcsr_release_p(matrix_c_im)
669 158 : CALL dbcsr_release_p(matrix_tmp)
670 158 : IF (ASSOCIATED(local_preconditioner)) THEN
671 158 : CALL destroy_preconditioner(local_preconditioner)
672 158 : DEALLOCATE (local_preconditioner)
673 : END IF
674 158 : CALL timestop(handle)
675 :
676 158 : END SUBROUTINE ot_eigensolver_complex
677 :
678 : ! **************************************************************************************************
679 : !> \brief diagonalize C^H H C and rotate a complex orthonormal orbital subspace
680 : !> \param matrix_h ...
681 : !> \param matrix_h_im ...
682 : !> \param coeff_re ...
683 : !> \param coeff_im ...
684 : !> \param eigenvalues ...
685 : ! **************************************************************************************************
686 158 : SUBROUTINE canonicalize_complex_subspace(matrix_h, matrix_h_im, coeff_re, coeff_im, eigenvalues)
687 : TYPE(dbcsr_type), POINTER :: matrix_h, matrix_h_im
688 : TYPE(cp_fm_type), INTENT(INOUT) :: coeff_re, coeff_im
689 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: eigenvalues
690 :
691 : INTEGER :: nao, nmo
692 : TYPE(cp_cfm_type) :: coeff, eigenvectors, hblock, hc, rotated
693 : TYPE(cp_fm_type) :: hc_im, hc_re
694 :
695 158 : CALL cp_fm_get_info(coeff_re, nrow_global=nao, ncol_global=nmo)
696 158 : CPASSERT(SIZE(eigenvalues) >= nmo)
697 158 : CALL cp_fm_create(hc_re, coeff_re%matrix_struct)
698 158 : CALL cp_fm_create(hc_im, coeff_re%matrix_struct)
699 158 : CALL cp_dbcsr_sm_fm_multiply(matrix_h, coeff_re, hc_re, nmo)
700 : CALL cp_dbcsr_sm_fm_multiply(matrix_h_im, coeff_im, hc_re, nmo, &
701 158 : alpha=-1.0_dp, beta=1.0_dp)
702 158 : CALL cp_dbcsr_sm_fm_multiply(matrix_h, coeff_im, hc_im, nmo)
703 : CALL cp_dbcsr_sm_fm_multiply(matrix_h_im, coeff_re, hc_im, nmo, &
704 158 : alpha=1.0_dp, beta=1.0_dp)
705 :
706 158 : CALL cp_cfm_create(coeff, coeff_re%matrix_struct)
707 158 : CALL cp_cfm_create(hc, coeff_re%matrix_struct)
708 158 : CALL cp_cfm_create(rotated, coeff_re%matrix_struct)
709 158 : CALL cp_cfm_create(hblock, coeff_re%matrix_struct, nrow=nmo, ncol=nmo)
710 158 : CALL cp_cfm_create(eigenvectors, coeff_re%matrix_struct, nrow=nmo, ncol=nmo)
711 158 : CALL cp_fm_to_cfm(coeff_re, coeff_im, coeff)
712 158 : CALL cp_fm_to_cfm(hc_re, hc_im, hc)
713 158 : CALL cp_cfm_gemm('C', 'N', nmo, nmo, nao, z_one, coeff, hc, z_zero, hblock)
714 158 : CALL cp_cfm_heevd(hblock, eigenvectors, eigenvalues(1:nmo))
715 158 : CALL cp_cfm_gemm('N', 'N', nao, nmo, nmo, z_one, coeff, eigenvectors, z_zero, rotated)
716 158 : CALL cp_cfm_to_fm(rotated, coeff_re, coeff_im)
717 :
718 158 : CALL cp_cfm_release(eigenvectors)
719 158 : CALL cp_cfm_release(hblock)
720 158 : CALL cp_cfm_release(rotated)
721 158 : CALL cp_cfm_release(hc)
722 158 : CALL cp_cfm_release(coeff)
723 158 : CALL cp_fm_release(hc_im)
724 158 : CALL cp_fm_release(hc_re)
725 158 : END SUBROUTINE canonicalize_complex_subspace
726 :
727 : END MODULE qs_ot_eigensolver
|
