(git:2d1258d)
Loading...
Searching...
No Matches
qs_ot_eigensolver.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 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! **************************************************************************************************
15 USE cp_cfm_diag, ONLY: cp_cfm_heevd
16 USE cp_cfm_types, ONLY: cp_cfm_create,&
21 USE cp_dbcsr_api, ONLY: &
23 dbcsr_scale, dbcsr_set, dbcsr_type, dbcsr_type_no_symmetry
34 USE cp_fm_types, ONLY: cp_fm_create,&
45 USE kinds, ONLY: dp
46 USE mathconstants, ONLY: z_one,&
47 z_zero
58 USE qs_mo_types, ONLY: mo_set_type
59 USE qs_ot, ONLY: qs_ot_get_orbitals,&
66 USE qs_ot_minimizer, ONLY: ot_mini
67 USE qs_ot_types, ONLY: qs_ot_allocate,&
74#include "./base/base_uses.f90"
75
76 IMPLICIT NONE
77 PRIVATE
78
79! *** Global parameters ***
80
81 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_ot_eigensolver'
82
83! *** Public subroutines ***
84
86
87CONTAINS
88
89! on input c contains the initial guess (should not be zero !)
90! on output c spans the subspace
91! **************************************************************************************************
92!> \brief ...
93!> \param matrix_h ...
94!> \param matrix_s ...
95!> \param matrix_orthogonal_space_fm ...
96!> \param matrix_c_fm ...
97!> \param preconditioner ...
98!> \param eps_gradient ...
99!> \param iter_max ...
100!> \param size_ortho_space ...
101!> \param silent ...
102!> \param ot_settings ...
103! **************************************************************************************************
104 SUBROUTINE ot_eigensolver(matrix_h, matrix_s, matrix_orthogonal_space_fm, &
105 matrix_c_fm, preconditioner, eps_gradient, &
106 iter_max, size_ortho_space, silent, ot_settings)
107
108 TYPE(dbcsr_type), POINTER :: matrix_h, matrix_s
109 TYPE(cp_fm_type), INTENT(IN), OPTIONAL :: matrix_orthogonal_space_fm
110 TYPE(cp_fm_type), INTENT(INOUT) :: matrix_c_fm
111 TYPE(preconditioner_type), OPTIONAL, POINTER :: preconditioner
112 REAL(kind=dp) :: eps_gradient
113 INTEGER, INTENT(IN) :: iter_max
114 INTEGER, INTENT(IN), OPTIONAL :: size_ortho_space
115 LOGICAL, INTENT(IN), OPTIONAL :: silent
116 TYPE(qs_ot_settings_type), INTENT(IN), OPTIONAL :: ot_settings
117
118 CHARACTER(len=*), PARAMETER :: routinen = 'ot_eigensolver'
119 INTEGER, PARAMETER :: max_iter_inner_loop = 40
120 REAL(kind=dp), PARAMETER :: rone = 1.0_dp, rzero = 0.0_dp
121
122 INTEGER :: handle, ieigensolver, iter_total, k, n, &
123 ortho_k, ortho_space_k, output_unit
124 LOGICAL :: energy_only, my_silent, ortho, &
125 ref_algorithm
126 REAL(kind=dp) :: delta, energy
127 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_hc
128 TYPE(dbcsr_type), POINTER :: matrix_buf1_ortho, matrix_buf2_ortho, &
129 matrix_c, matrix_orthogonal_space, &
130 matrix_os_ortho, matrix_s_ortho
131 TYPE(qs_ot_type), DIMENSION(:), POINTER :: qs_ot_env
132
133 CALL timeset(routinen, handle)
134
135 output_unit = cp_logger_get_default_io_unit()
136
137 IF (PRESENT(silent)) THEN
138 my_silent = silent
139 ELSE
140 my_silent = .false.
141 END IF
142
143 NULLIFY (matrix_c) ! fm->dbcsr
144
145 CALL cp_fm_get_info(matrix_c_fm, nrow_global=n, ncol_global=k) ! fm->dbcsr
146 ALLOCATE (matrix_c)
147 ref_algorithm = .false.
148 IF (PRESENT(ot_settings)) ref_algorithm = ot_settings%ot_algorithm == "REF"
149 IF (ref_algorithm) THEN
150 CALL dbcsr_init_p(matrix_c)
151 CALL cp_dbcsr_m_by_n_from_row_template(matrix_c, template=matrix_h, n=k, &
152 sym=dbcsr_type_no_symmetry)
153 CALL copy_fm_to_dbcsr(matrix_c_fm, matrix_c)
154 ELSE
155 CALL cp_fm_to_dbcsr_row_template(matrix_c, fm_in=matrix_c_fm, template=matrix_h)
156 END IF
157
158 iter_total = 0
159
160 outer_scf: DO
161
162 NULLIFY (qs_ot_env)
163
164 NULLIFY (matrix_s_ortho)
165 NULLIFY (matrix_os_ortho)
166 NULLIFY (matrix_buf1_ortho)
167 NULLIFY (matrix_buf2_ortho)
168 NULLIFY (matrix_orthogonal_space)
169
170 ALLOCATE (qs_ot_env(1))
171 ALLOCATE (matrix_hc(1))
172 NULLIFY (matrix_hc(1)%matrix)
173 CALL dbcsr_init_p(matrix_hc(1)%matrix)
174
175 ortho = .false.
176 IF (PRESENT(matrix_orthogonal_space_fm)) ortho = .true.
177
178 ! decide settings
179 IF (PRESENT(ot_settings)) THEN
180 qs_ot_env(1)%settings = ot_settings
181 ELSE
182 CALL qs_ot_settings_init(qs_ot_env(1)%settings)
183 ! overwrite defaults
184 qs_ot_env(1)%settings%ds_min = 0.10_dp
185 END IF
186
187 IF (ortho) THEN
188 ALLOCATE (matrix_orthogonal_space)
189 CALL cp_fm_to_dbcsr_row_template(matrix_orthogonal_space, fm_in=matrix_orthogonal_space_fm, template=matrix_h)
190 CALL cp_fm_get_info(matrix_orthogonal_space_fm, ncol_global=ortho_space_k)
191
192 IF (PRESENT(size_ortho_space)) ortho_space_k = size_ortho_space
193 ortho_k = ortho_space_k + k
194 ELSE
195 ortho_k = k
196 END IF
197
198 ! allocate
199 CALL qs_ot_allocate(qs_ot_env(1), matrix_s, matrix_c_fm%matrix_struct, ortho_k=ortho_k)
200 IF (ref_algorithm) THEN
201 CALL dbcsr_copy(matrix_hc(1)%matrix, qs_ot_env(1)%matrix_x, 'matrix_hc')
202 ELSE
203 CALL dbcsr_copy(matrix_hc(1)%matrix, matrix_c, 'matrix_hc')
204 END IF
205
206 IF (ortho) THEN
207 ! construct an initial guess that is orthogonal to matrix_orthogonal_space
208
209 CALL dbcsr_init_p(matrix_s_ortho)
210 CALL dbcsr_copy(matrix_s_ortho, matrix_orthogonal_space, name="matrix_s_ortho")
211
212 CALL dbcsr_init_p(matrix_os_ortho)
213 CALL cp_dbcsr_m_by_n_from_template(matrix_os_ortho, template=matrix_h, m=ortho_space_k, n=ortho_space_k, &
214 sym=dbcsr_type_no_symmetry)
215
216 CALL dbcsr_init_p(matrix_buf1_ortho)
217 CALL cp_dbcsr_m_by_n_from_template(matrix_buf1_ortho, template=matrix_h, m=ortho_space_k, n=k, &
218 sym=dbcsr_type_no_symmetry)
219
220 CALL dbcsr_init_p(matrix_buf2_ortho)
221 CALL cp_dbcsr_m_by_n_from_template(matrix_buf2_ortho, template=matrix_h, m=ortho_space_k, n=k, &
222 sym=dbcsr_type_no_symmetry)
223
224 CALL dbcsr_multiply('N', 'N', 1.0_dp, matrix_s, matrix_orthogonal_space, &
225 0.0_dp, matrix_s_ortho)
226 CALL dbcsr_multiply('T', 'N', rone, matrix_s_ortho, matrix_s_ortho, &
227 rzero, matrix_os_ortho)
228
229 CALL cp_dbcsr_cholesky_decompose(matrix_os_ortho, &
230 para_env=qs_ot_env(1)%para_env, blacs_env=qs_ot_env(1)%blacs_env)
231 CALL cp_dbcsr_cholesky_invert(matrix_os_ortho, &
232 para_env=qs_ot_env(1)%para_env, blacs_env=qs_ot_env(1)%blacs_env, &
233 uplo_to_full=.true.)
234
235 CALL dbcsr_multiply('T', 'N', rone, matrix_s_ortho, matrix_c, &
236 rzero, matrix_buf1_ortho)
237 CALL dbcsr_multiply('N', 'N', rone, matrix_os_ortho, matrix_buf1_ortho, &
238 rzero, matrix_buf2_ortho)
239 CALL dbcsr_multiply('N', 'N', -rone, matrix_s_ortho, matrix_buf2_ortho, &
240 rone, matrix_c)
241
242 ! make matrix_c0 an orthogonal basis, matrix_c contains sc0
243 CALL dbcsr_copy(qs_ot_env(1)%matrix_c0, matrix_c)
244 CALL dbcsr_multiply('N', 'N', 1.0_dp, matrix_s, qs_ot_env(1)%matrix_c0, &
245 0.0_dp, matrix_c)
246
247 CALL make_basis_sv(qs_ot_env(1)%matrix_c0, k, matrix_c, &
248 qs_ot_env(1)%para_env, qs_ot_env(1)%blacs_env)
249
250 ! copy sc0 and matrix_s_ortho in qs_ot_env(1)%matrix_sc0
251 !CALL dbcsr_copy_columns(qs_ot_env(1)%matrix_sc0,matrix_s_ortho,ortho_space_k,1,1)
252 CALL dbcsr_copy_columns_hack(qs_ot_env(1)%matrix_sc0, matrix_s_ortho, ortho_space_k, 1, 1, &
253 para_env=qs_ot_env(1)%para_env, blacs_env=qs_ot_env(1)%blacs_env)
254 !CALL dbcsr_copy_columns(qs_ot_env(1)%matrix_sc0,matrix_c,k,1,ortho_space_k+1)
255 CALL dbcsr_copy_columns_hack(qs_ot_env(1)%matrix_sc0, matrix_c, k, 1, ortho_space_k + 1, &
256 para_env=qs_ot_env(1)%para_env, blacs_env=qs_ot_env(1)%blacs_env)
257
258 CALL dbcsr_release_p(matrix_buf1_ortho)
259 CALL dbcsr_release_p(matrix_buf2_ortho)
260 CALL dbcsr_release_p(matrix_os_ortho)
261 CALL dbcsr_release_p(matrix_s_ortho)
262
263 ELSE
264
265 ! set c0,sc0
266 CALL dbcsr_copy(qs_ot_env(1)%matrix_c0, matrix_c)
267 CALL dbcsr_multiply('N', 'N', 1.0_dp, matrix_s, qs_ot_env(1)%matrix_c0, &
268 0.0_dp, qs_ot_env(1)%matrix_sc0)
269
270 CALL make_basis_sv(qs_ot_env(1)%matrix_c0, k, qs_ot_env(1)%matrix_sc0, &
271 qs_ot_env(1)%para_env, qs_ot_env(1)%blacs_env)
272 END IF
273
274 ! init
275 CALL qs_ot_init(qs_ot_env(1))
276 energy_only = qs_ot_env(1)%energy_only
277
278 SELECT CASE (qs_ot_env(1)%settings%ot_algorithm)
279 CASE ("TOD")
280 CALL dbcsr_set(qs_ot_env(1)%matrix_x, 0.0_dp)
281 CALL dbcsr_set(qs_ot_env(1)%matrix_sx, 0.0_dp)
282 CALL qs_ot_get_p(qs_ot_env(1)%matrix_x, qs_ot_env(1)%matrix_sx, qs_ot_env(1))
283 CALL qs_ot_get_orbitals(matrix_c, qs_ot_env(1)%matrix_x, qs_ot_env(1))
284 CASE ("REF")
285 CALL dbcsr_copy(qs_ot_env(1)%matrix_x, qs_ot_env(1)%matrix_c0)
286 CALL dbcsr_copy(qs_ot_env(1)%matrix_sx, qs_ot_env(1)%matrix_sc0)
287 CALL dbcsr_copy(matrix_c, qs_ot_env(1)%matrix_c0)
288 CASE DEFAULT
289 cpabort("OT eigensolver supports ALGORITHM STRICT or IRAC")
290 END SELECT
291
292 ! if present preconditioner, use it
293
294 IF (PRESENT(preconditioner)) THEN
295 IF (ASSOCIATED(preconditioner)) THEN
297 CALL qs_ot_new_preconditioner(qs_ot_env(1), preconditioner)
298 ELSE
299 ! we should presumably make one
300 END IF
301 END IF
302 END IF
303
304 ! *** Eigensolver loop ***
305 ieigensolver = 0
306 eigensolver_loop: DO
307
308 ieigensolver = ieigensolver + 1
309 iter_total = iter_total + 1
310
311 ! the energy is cHc, the gradient is 2*H*c
312 CALL dbcsr_multiply('N', 'N', 1.0_dp, matrix_h, matrix_c, &
313 0.0_dp, matrix_hc(1)%matrix)
314 CALL dbcsr_dot(matrix_c, matrix_hc(1)%matrix, energy)
315 IF (.NOT. energy_only) THEN
316 CALL dbcsr_scale(matrix_hc(1)%matrix, 2.0_dp)
317 END IF
318
319 qs_ot_env(1)%etotal = energy
320 CALL ot_mini(qs_ot_env, matrix_hc)
321 delta = qs_ot_env(1)%delta
322 energy_only = qs_ot_env(1)%energy_only
323
324 SELECT CASE (qs_ot_env(1)%settings%ot_algorithm)
325 CASE ("TOD")
326 CALL dbcsr_multiply('N', 'N', 1.0_dp, matrix_s, qs_ot_env(1)%matrix_x, &
327 0.0_dp, qs_ot_env(1)%matrix_sx)
328 CALL qs_ot_get_p(qs_ot_env(1)%matrix_x, qs_ot_env(1)%matrix_sx, qs_ot_env(1))
329 CALL qs_ot_get_orbitals(matrix_c, qs_ot_env(1)%matrix_x, qs_ot_env(1))
330 CASE ("REF")
331 CALL qs_ot_get_orbitals_ref(matrix_c, matrix_s, qs_ot_env(1)%matrix_x, &
332 qs_ot_env(1)%matrix_sx, qs_ot_env(1)%matrix_gx_old, &
333 qs_ot_env(1)%matrix_dx, qs_ot_env(1), qs_ot_env(1))
334 END SELECT
335
336 ! exit on convergence or if maximum of inner loop cycles is reached
337 IF (delta < eps_gradient .OR. ieigensolver >= max_iter_inner_loop) EXIT eigensolver_loop
338 ! exit if total number of steps is reached, but not during a line search step
339 IF (iter_total >= iter_max .AND. qs_ot_env(1)%OT_METHOD_FULL /= "OT LS") EXIT eigensolver_loop
340
341 END DO eigensolver_loop
342
343 CALL qs_ot_destroy(qs_ot_env(1))
344 DEALLOCATE (qs_ot_env)
345 CALL dbcsr_release_p(matrix_hc(1)%matrix)
346 DEALLOCATE (matrix_hc)
347 CALL dbcsr_release_p(matrix_orthogonal_space)
348
349 IF (delta < eps_gradient) THEN
350 IF ((output_unit > 0) .AND. .NOT. my_silent) THEN
351 WRITE (unit=output_unit, fmt="(T2,A,I0,A)") &
352 "OT| Eigensolver reached convergence in ", iter_total, " iterations"
353 END IF
354 EXIT outer_scf
355 END IF
356 IF (iter_total >= iter_max) THEN
357 IF (output_unit > 0) THEN
358 IF (my_silent) THEN
359 WRITE (output_unit, "(A,T60,E20.10)") " WARNING OT eigensolver did not converge: current gradient", delta
360 ELSE
361 WRITE (output_unit, *) "WARNING : did not converge in ot_eigensolver"
362 WRITE (output_unit, *) "number of iterations ", iter_total, " exceeded maximum"
363 WRITE (output_unit, *) "current gradient / target gradient", delta, " / ", eps_gradient
364 END IF
365 END IF
366 EXIT outer_scf
367 END IF
368
369 END DO outer_scf
370
371 CALL copy_dbcsr_to_fm(matrix_c, matrix_c_fm) ! fm->dbcsr
372 CALL dbcsr_release_p(matrix_c) ! fm->dbcsr
373
374 CALL timestop(handle)
375
376 END SUBROUTINE ot_eigensolver
377
378! **************************************************************************************************
379!> \brief solve a fixed complex Hermitian generalized eigenproblem by OT
380!> \param matrix_h real part of H(k)
381!> \param matrix_h_im imaginary part of H(k)
382!> \param matrix_s real part of S(k)
383!> \param matrix_s_im imaginary part of S(k)
384!> \param matrix_c_fm real part of the orbital coefficients
385!> \param matrix_c_fm_im imaginary part of the orbital coefficients
386!> \param preconditioner optional complex k-point preconditioner
387!> \param eps_gradient requested OT gradient accuracy
388!> \param iter_max maximum number of OT iterations
389!> \param eigenvalues ...
390!> \param silent suppress successful convergence output
391!> \param ot_settings OT algorithm and minimizer settings
392!> \param matrix_t ...
393!> \param matrix_t_im ...
394!> \param mo_set ...
395! **************************************************************************************************
396 SUBROUTINE ot_eigensolver_complex(matrix_h, matrix_h_im, matrix_s, matrix_s_im, &
397 matrix_c_fm, matrix_c_fm_im, preconditioner, &
398 eps_gradient, iter_max, eigenvalues, silent, ot_settings, &
399 matrix_t, matrix_t_im, mo_set)
400
401 TYPE(dbcsr_type), POINTER :: matrix_h, matrix_h_im, matrix_s, &
402 matrix_s_im
403 TYPE(cp_fm_type), INTENT(INOUT) :: matrix_c_fm, matrix_c_fm_im
404 TYPE(preconditioner_type), OPTIONAL, POINTER :: preconditioner
405 REAL(kind=dp), INTENT(IN) :: eps_gradient
406 INTEGER, INTENT(IN) :: iter_max
407 REAL(kind=dp), DIMENSION(:), INTENT(OUT), OPTIONAL :: eigenvalues
408 LOGICAL, INTENT(IN), OPTIONAL :: silent
409 TYPE(qs_ot_settings_type), INTENT(IN), OPTIONAL :: ot_settings
410 TYPE(dbcsr_type), OPTIONAL, POINTER :: matrix_t, matrix_t_im
411 TYPE(mo_set_type), INTENT(IN), OPTIONAL :: mo_set
412
413 CHARACTER(len=*), PARAMETER :: routinen = 'ot_eigensolver_complex'
414 INTEGER, PARAMETER :: max_iter_inner_loop = 40
415
416 INTEGER :: handle, ieigensolver, iter_total, k, n, &
417 output_unit
418 LOGICAL :: energy_only, my_silent, &
419 preconditioner_rejected
420 REAL(kind=dp) :: delta, energy_im, energy_re
421 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_hc, matrix_hc_im
422 TYPE(dbcsr_type), POINTER :: matrix_c, matrix_c_im, matrix_tmp
423 TYPE(preconditioner_type), POINTER :: local_preconditioner
424 TYPE(qs_ot_type), DIMENSION(:), POINTER :: qs_ot_env
425
426 CALL timeset(routinen, handle)
427 output_unit = cp_logger_get_default_io_unit()
428 my_silent = .false.
429 IF (PRESENT(silent)) my_silent = silent
430 NULLIFY (local_preconditioner)
431
432 CALL cp_fm_get_info(matrix_c_fm, nrow_global=n, ncol_global=k)
433 cpassert(n > 0 .AND. k > 0)
434 ALLOCATE (matrix_c, matrix_c_im, matrix_tmp)
435 CALL dbcsr_init_p(matrix_c)
436 CALL cp_dbcsr_m_by_n_from_row_template(matrix_c, template=matrix_h, n=k, &
437 sym=dbcsr_type_no_symmetry)
438 CALL dbcsr_init_p(matrix_c_im)
439 CALL cp_dbcsr_m_by_n_from_row_template(matrix_c_im, template=matrix_h, n=k, &
440 sym=dbcsr_type_no_symmetry)
441 CALL dbcsr_init_p(matrix_tmp)
442 CALL cp_dbcsr_m_by_n_from_row_template(matrix_tmp, template=matrix_h, n=k, &
443 sym=dbcsr_type_no_symmetry)
444 CALL copy_fm_to_dbcsr(matrix_c_fm, matrix_c)
445 CALL copy_fm_to_dbcsr(matrix_c_fm_im, matrix_c_im)
446
447 iter_total = 0
448 outer_scf: DO
449 ALLOCATE (qs_ot_env(1), matrix_hc(1), matrix_hc_im(1))
450 NULLIFY (matrix_hc(1)%matrix, matrix_hc_im(1)%matrix)
451
452 IF (PRESENT(ot_settings)) THEN
453 qs_ot_env(1)%settings = ot_settings
454 ELSE
455 CALL qs_ot_settings_init(qs_ot_env(1)%settings)
456 qs_ot_env(1)%settings%ds_min = 0.10_dp
457 END IF
458
459 CALL qs_ot_allocate(qs_ot_env(1), matrix_s, matrix_c_fm%matrix_struct)
460 CALL qs_ot_allocate_complex_state(qs_ot_env(1), matrix_s)
461 CALL dbcsr_init_p(matrix_hc(1)%matrix)
462 CALL dbcsr_copy(matrix_hc(1)%matrix, qs_ot_env(1)%matrix_x, 'matrix_hc')
463 CALL dbcsr_init_p(matrix_hc_im(1)%matrix)
464 CALL dbcsr_copy(matrix_hc_im(1)%matrix, qs_ot_env(1)%matrix_x_im, 'matrix_hc_im')
465
466 CALL dbcsr_copy(qs_ot_env(1)%matrix_c0, matrix_c)
467 CALL dbcsr_copy(qs_ot_env(1)%matrix_c0_im, matrix_c_im)
468 CALL dbcsr_multiply('N', 'N', 1.0_dp, matrix_s, matrix_c, &
469 0.0_dp, qs_ot_env(1)%matrix_sc0)
470 CALL dbcsr_multiply('N', 'N', 1.0_dp, matrix_s_im, matrix_c_im, &
471 0.0_dp, matrix_tmp)
472 CALL dbcsr_add(qs_ot_env(1)%matrix_sc0, matrix_tmp, &
473 alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
474 CALL dbcsr_multiply('N', 'N', 1.0_dp, matrix_s, matrix_c_im, &
475 0.0_dp, qs_ot_env(1)%matrix_sc0_im)
476 CALL dbcsr_multiply('N', 'N', 1.0_dp, matrix_s_im, matrix_c, &
477 0.0_dp, matrix_tmp)
478 CALL dbcsr_add(qs_ot_env(1)%matrix_sc0_im, matrix_tmp, &
479 alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
480
481 CALL qs_ot_init(qs_ot_env(1))
482 energy_only = qs_ot_env(1)%energy_only
483 SELECT CASE (qs_ot_env(1)%settings%ot_algorithm)
484 CASE ('TOD')
485 CALL dbcsr_set(qs_ot_env(1)%matrix_x, 0.0_dp)
486 CALL dbcsr_set(qs_ot_env(1)%matrix_x_im, 0.0_dp)
487 CALL dbcsr_set(qs_ot_env(1)%matrix_sx, 0.0_dp)
488 CALL dbcsr_set(qs_ot_env(1)%matrix_sx_im, 0.0_dp)
489 CALL qs_ot_get_orbitals_complex(matrix_c, matrix_c_im, matrix_s, matrix_s_im, &
490 qs_ot_env(1))
491 CASE ('REF')
492 CALL dbcsr_copy(qs_ot_env(1)%matrix_x, qs_ot_env(1)%matrix_c0)
493 CALL dbcsr_copy(qs_ot_env(1)%matrix_x_im, qs_ot_env(1)%matrix_c0_im)
494 CALL dbcsr_copy(qs_ot_env(1)%matrix_sx, qs_ot_env(1)%matrix_sc0)
495 CALL dbcsr_copy(qs_ot_env(1)%matrix_sx_im, qs_ot_env(1)%matrix_sc0_im)
496 CALL qs_ot_get_orbitals_ref_complex(matrix_c, matrix_c_im, matrix_s, matrix_s_im, &
497 qs_ot_env(1), qs_ot_env(1))
498 CASE DEFAULT
499 cpabort('Complex OT eigensolver supports ALGORITHM STRICT or IRAC')
500 END SELECT
501
502 IF (.NOT. PRESENT(preconditioner) .AND. .NOT. ASSOCIATED(local_preconditioner) .AND. &
503 qs_ot_env(1)%settings%preconditioner_type /= ot_precond_none) THEN
504 ALLOCATE (local_preconditioner)
505 CALL init_preconditioner(local_preconditioner, qs_ot_env(1)%para_env, qs_ot_env(1)%blacs_env)
506 SELECT CASE (qs_ot_env(1)%settings%preconditioner_type)
508 IF (qs_ot_env(1)%settings%ot_algorithm == 'TOD') THEN
509 ! The column-selective FULL_ALL model can lose its positive tangent metric in a
510 ! finite STRICT fixed-H chart. Retain a well-defined eigensolver with its S^-1 base.
512 local_preconditioner, matrix_s, matrix_s_im, &
513 qs_ot_env(1)%settings%precond_solver_type)
514 ELSE
515 IF (.NOT. PRESENT(mo_set)) THEN
516 cpabort('Complex FULL_ALL eigensolver preconditioning requires MO energy labels')
517 END IF
519 local_preconditioner, qs_ot_env(1)%matrix_c0, qs_ot_env(1)%matrix_c0_im, &
520 matrix_h, matrix_h_im, matrix_s, matrix_s_im, mo_set, &
521 qs_ot_env(1)%settings%energy_gap, qs_ot_env(1)%settings%precond_solver_type)
522 END IF
524 IF (.NOT. PRESENT(mo_set)) THEN
525 cpabort('Complex FULL_SINGLE eigensolver preconditioning requires MO energy labels')
526 END IF
528 local_preconditioner, matrix_h, matrix_h_im, matrix_s, matrix_s_im, mo_set, &
529 qs_ot_env(1)%settings%energy_gap, qs_ot_env(1)%settings%precond_solver_type)
532 local_preconditioner, qs_ot_env(1)%matrix_c0, qs_ot_env(1)%matrix_c0_im, &
533 matrix_h, matrix_h_im, matrix_s, matrix_s_im, &
534 qs_ot_env(1)%settings%energy_gap, qs_ot_env(1)%settings%precond_solver_type)
536 IF (.NOT. PRESENT(matrix_t) .OR. .NOT. PRESENT(matrix_t_im)) THEN
537 cpabort('Complex FULL_KINETIC eigensolver preconditioning requires T(k)')
538 END IF
539 IF (.NOT. ASSOCIATED(matrix_t) .OR. .NOT. ASSOCIATED(matrix_t_im)) THEN
540 cpabort('Complex FULL_KINETIC eigensolver preconditioning requires T(k)')
541 END IF
543 local_preconditioner, matrix_t, matrix_t_im, matrix_s, matrix_s_im, &
544 qs_ot_env(1)%settings%energy_gap, qs_ot_env(1)%settings%precond_solver_type)
547 local_preconditioner, matrix_s, matrix_s_im, &
548 qs_ot_env(1)%settings%precond_solver_type)
549 CASE DEFAULT
550 cpabort('Unsupported complex K-point OT eigensolver preconditioner')
551 END SELECT
552 END IF
553
554 IF (PRESENT(preconditioner)) THEN
555 IF (ASSOCIATED(preconditioner)) THEN
557 CALL qs_ot_new_preconditioner(qs_ot_env(1), preconditioner)
558 END IF
559 END IF
560 ELSE IF (ASSOCIATED(local_preconditioner)) THEN
561 CALL qs_ot_new_preconditioner(qs_ot_env(1), local_preconditioner)
562 END IF
563
564 IF (qs_ot_env(1)%settings%ot_algorithm == 'TOD') THEN
565 CALL qs_ot_prepare_complex_tangent_metric(qs_ot_env(1), preconditioner_rejected)
566 IF (preconditioner_rejected .AND. ASSOCIATED(local_preconditioner)) THEN
567 CALL destroy_preconditioner(local_preconditioner)
568 CALL init_preconditioner(local_preconditioner, &
569 qs_ot_env(1)%para_env, qs_ot_env(1)%blacs_env)
571 local_preconditioner, matrix_s, matrix_s_im, &
572 qs_ot_env(1)%settings%precond_solver_type)
573 CALL qs_ot_new_preconditioner(qs_ot_env(1), local_preconditioner)
574 CALL qs_ot_prepare_complex_tangent_metric(qs_ot_env(1), preconditioner_rejected)
575 cpassert(.NOT. preconditioner_rejected)
576 END IF
577 END IF
578
579 ieigensolver = 0
580 eigensolver_loop: DO
581 ieigensolver = ieigensolver + 1
582 iter_total = iter_total + 1
583
584 ! H*C = (Hre*Cre-Him*Cim) + i*(Hre*Cim+Him*Cre).
585 CALL dbcsr_multiply('N', 'N', 1.0_dp, matrix_h, matrix_c, &
586 0.0_dp, matrix_hc(1)%matrix)
587 CALL dbcsr_multiply('N', 'N', 1.0_dp, matrix_h_im, matrix_c_im, &
588 0.0_dp, matrix_tmp)
589 CALL dbcsr_add(matrix_hc(1)%matrix, matrix_tmp, &
590 alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
591 CALL dbcsr_multiply('N', 'N', 1.0_dp, matrix_h, matrix_c_im, &
592 0.0_dp, matrix_hc_im(1)%matrix)
593 CALL dbcsr_multiply('N', 'N', 1.0_dp, matrix_h_im, matrix_c, &
594 0.0_dp, matrix_tmp)
595 CALL dbcsr_add(matrix_hc_im(1)%matrix, matrix_tmp, &
596 alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
597
598 CALL dbcsr_dot(matrix_c, matrix_hc(1)%matrix, energy_re)
599 CALL dbcsr_dot(matrix_c_im, matrix_hc_im(1)%matrix, energy_im)
600 qs_ot_env(1)%etotal = energy_re + energy_im
601 IF (.NOT. energy_only) THEN
602 CALL dbcsr_scale(matrix_hc(1)%matrix, 2.0_dp)
603 CALL dbcsr_scale(matrix_hc_im(1)%matrix, 2.0_dp)
604 END IF
605 CALL ot_mini(qs_ot_env, matrix_hc, matrix_hc_im=matrix_hc_im)
606 delta = qs_ot_env(1)%delta
607 energy_only = qs_ot_env(1)%energy_only
608
609 SELECT CASE (qs_ot_env(1)%settings%ot_algorithm)
610 CASE ('TOD')
611 CALL qs_ot_get_orbitals_complex(matrix_c, matrix_c_im, matrix_s, matrix_s_im, &
612 qs_ot_env(1))
613 CASE ('REF')
614 CALL qs_ot_get_orbitals_ref_complex(matrix_c, matrix_c_im, matrix_s, matrix_s_im, &
615 qs_ot_env(1), qs_ot_env(1))
616 END SELECT
617
618 IF (delta < eps_gradient .OR. ieigensolver >= max_iter_inner_loop) EXIT eigensolver_loop
619 IF (iter_total >= iter_max .AND. qs_ot_env(1)%OT_METHOD_FULL /= 'OT LS') EXIT eigensolver_loop
620 END DO eigensolver_loop
621
622 CALL qs_ot_destroy(qs_ot_env(1))
623 DEALLOCATE (qs_ot_env)
624 CALL dbcsr_release_p(matrix_hc(1)%matrix)
625 CALL dbcsr_release_p(matrix_hc_im(1)%matrix)
626 DEALLOCATE (matrix_hc, matrix_hc_im)
627
628 IF (delta < eps_gradient) THEN
629 IF (output_unit > 0 .AND. .NOT. my_silent) THEN
630 WRITE (unit=output_unit, fmt='(T2,A,I0,A)') &
631 'OT| Complex eigensolver reached convergence in ', iter_total, ' iterations'
632 END IF
633 EXIT outer_scf
634 END IF
635 IF (iter_total >= iter_max) THEN
636 IF (output_unit > 0) THEN
637 WRITE (output_unit, '(A,T60,E20.10)') &
638 ' WARNING complex OT eigensolver did not converge: current gradient', delta
639 END IF
640 EXIT outer_scf
641 END IF
642 END DO outer_scf
643
644 CALL copy_dbcsr_to_fm(matrix_c, matrix_c_fm)
645 CALL copy_dbcsr_to_fm(matrix_c_im, matrix_c_fm_im)
646 IF (PRESENT(eigenvalues)) THEN
647 CALL canonicalize_complex_subspace(matrix_h, matrix_h_im, matrix_c_fm, &
648 matrix_c_fm_im, eigenvalues)
649 END IF
650 CALL dbcsr_release_p(matrix_c)
651 CALL dbcsr_release_p(matrix_c_im)
652 CALL dbcsr_release_p(matrix_tmp)
653 IF (ASSOCIATED(local_preconditioner)) THEN
654 CALL destroy_preconditioner(local_preconditioner)
655 DEALLOCATE (local_preconditioner)
656 END IF
657 CALL timestop(handle)
658
659 END SUBROUTINE ot_eigensolver_complex
660
661! **************************************************************************************************
662!> \brief diagonalize C^H H C and rotate a complex orthonormal orbital subspace
663!> \param matrix_h ...
664!> \param matrix_h_im ...
665!> \param coeff_re ...
666!> \param coeff_im ...
667!> \param eigenvalues ...
668! **************************************************************************************************
669 SUBROUTINE canonicalize_complex_subspace(matrix_h, matrix_h_im, coeff_re, coeff_im, eigenvalues)
670 TYPE(dbcsr_type), POINTER :: matrix_h, matrix_h_im
671 TYPE(cp_fm_type), INTENT(INOUT) :: coeff_re, coeff_im
672 REAL(kind=dp), DIMENSION(:), INTENT(OUT) :: eigenvalues
673
674 INTEGER :: nao, nmo
675 TYPE(cp_cfm_type) :: coeff, eigenvectors, hblock, hc, rotated
676 TYPE(cp_fm_type) :: hc_im, hc_re
677
678 CALL cp_fm_get_info(coeff_re, nrow_global=nao, ncol_global=nmo)
679 cpassert(SIZE(eigenvalues) >= nmo)
680 CALL cp_fm_create(hc_re, coeff_re%matrix_struct)
681 CALL cp_fm_create(hc_im, coeff_re%matrix_struct)
682 CALL cp_dbcsr_sm_fm_multiply(matrix_h, coeff_re, hc_re, nmo)
683 CALL cp_dbcsr_sm_fm_multiply(matrix_h_im, coeff_im, hc_re, nmo, &
684 alpha=-1.0_dp, beta=1.0_dp)
685 CALL cp_dbcsr_sm_fm_multiply(matrix_h, coeff_im, hc_im, nmo)
686 CALL cp_dbcsr_sm_fm_multiply(matrix_h_im, coeff_re, hc_im, nmo, &
687 alpha=1.0_dp, beta=1.0_dp)
688
689 CALL cp_cfm_create(coeff, coeff_re%matrix_struct)
690 CALL cp_cfm_create(hc, coeff_re%matrix_struct)
691 CALL cp_cfm_create(rotated, coeff_re%matrix_struct)
692 CALL cp_cfm_create(hblock, coeff_re%matrix_struct, nrow=nmo, ncol=nmo)
693 CALL cp_cfm_create(eigenvectors, coeff_re%matrix_struct, nrow=nmo, ncol=nmo)
694 CALL cp_fm_to_cfm(coeff_re, coeff_im, coeff)
695 CALL cp_fm_to_cfm(hc_re, hc_im, hc)
696 CALL cp_cfm_gemm('C', 'N', nmo, nmo, nao, z_one, coeff, hc, z_zero, hblock)
697 CALL cp_cfm_heevd(hblock, eigenvectors, eigenvalues(1:nmo))
698 CALL cp_cfm_gemm('N', 'N', nao, nmo, nmo, z_one, coeff, eigenvectors, z_zero, rotated)
699 CALL cp_cfm_to_fm(rotated, coeff_re, coeff_im)
700
701 CALL cp_cfm_release(eigenvectors)
702 CALL cp_cfm_release(hblock)
703 CALL cp_cfm_release(rotated)
704 CALL cp_cfm_release(hc)
705 CALL cp_cfm_release(coeff)
706 CALL cp_fm_release(hc_im)
707 CALL cp_fm_release(hc_re)
708 END SUBROUTINE canonicalize_complex_subspace
709
710END MODULE qs_ot_eigensolver
Basic linear algebra operations for complex full matrices.
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 ) + ...
used for collecting diagonalization schemes available for cp_cfm_type
Definition cp_cfm_diag.F:14
subroutine, public cp_cfm_heevd(matrix, eigenvectors, eigenvalues)
Perform a diagonalisation of a complex matrix.
Definition cp_cfm_diag.F:82
Represents a complex full matrix distributed on many processors.
subroutine, public cp_cfm_release(matrix)
Releases a full matrix.
subroutine, public cp_fm_to_cfm(msourcer, msourcei, mtarget)
Construct a complex full matrix by taking its real and imaginary parts from two separate real-value f...
subroutine, public cp_cfm_create(matrix, matrix_struct, name, nrow, ncol, set_zero)
Creates a new full matrix with the given structure.
subroutine, public cp_cfm_to_fm(msource, mtargetr, mtargeti)
Copy real and imaginary parts of a complex full matrix into separate real-value full matrices.
subroutine, public dbcsr_release_p(matrix)
...
subroutine, public dbcsr_scale(matrix, alpha_scalar)
...
subroutine, public dbcsr_copy(matrix_b, matrix_a, name, keep_sparsity, keep_imaginary)
...
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)
...
subroutine, public dbcsr_init_p(matrix)
...
subroutine, public dbcsr_set(matrix, alpha)
...
subroutine, public dbcsr_add(matrix_a, matrix_b, alpha_scalar, beta_scalar)
...
Interface to (sca)lapack for the Cholesky based procedures.
subroutine, public cp_dbcsr_cholesky_decompose(matrix, n, para_env, blacs_env)
used to replace a symmetric positive def. matrix M with its cholesky decomposition U: M = U^T * U,...
subroutine, public cp_dbcsr_cholesky_invert(matrix, n, para_env, blacs_env, uplo_to_full)
used to replace the cholesky decomposition by the inverse
subroutine, public dbcsr_dot(matrix_a, matrix_b, trace)
Computes the dot product of two matrices, also known as the trace of their matrix product.
DBCSR operations in CP2K.
subroutine, public cp_dbcsr_sm_fm_multiply(matrix, fm_in, fm_out, ncol, alpha, beta)
multiply a dbcsr with a fm matrix
subroutine, public copy_dbcsr_to_fm(matrix, fm)
Copy a DBCSR matrix to a BLACS matrix.
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...
subroutine, public cp_dbcsr_m_by_n_from_row_template(matrix, template, n, sym)
Utility function to create dbcsr matrix, m x n matrix (n arbitrary) with the same processor grid and ...
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...
subroutine, public dbcsr_copy_columns_hack(matrix_b, matrix_a, ncol, source_start, target_start, para_env, blacs_env)
hack for dbcsr_copy_columns
subroutine, public copy_fm_to_dbcsr(fm, matrix, keep_sparsity)
Copy a BLACS matrix to a dbcsr matrix.
represent a full matrix distributed on many processors
Definition cp_fm_types.F:15
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
subroutine, public cp_fm_create(matrix, matrix_struct, name, nrow, ncol, set_zero)
creates a new full matrix with the given structure
various routines to log and control the output. The idea is that decisions about where to log should ...
integer function, public cp_logger_get_default_io_unit(logger)
returns the unit nr for the ionode (-1 on all other processors) skips as well checks if the procs cal...
collects all constants needed in input so that they can be used without circular dependencies
integer, parameter, public ot_precond_full_kinetic
integer, parameter, public ot_precond_full_single
integer, parameter, public ot_precond_none
integer, parameter, public ot_precond_full_single_inverse
integer, parameter, public ot_precond_s_inverse
integer, parameter, public ot_precond_full_all
Defines the basic variable types.
Definition kinds.F:23
integer, parameter, public dp
Definition kinds.F:34
Definition of mathematical constants and functions.
complex(kind=dp), parameter, public z_one
complex(kind=dp), parameter, public z_zero
types of preconditioners
subroutine, public init_preconditioner(preconditioner_env, para_env, blacs_env)
...
subroutine, public destroy_preconditioner(preconditioner_env)
...
logical function, public preconditioner_in_use(preconditioner)
...
computes preconditioners, and implements methods to apply them currently used in qs_ot
subroutine, public make_preconditioner_complex_full_all(preconditioner_env, matrix_c_re, matrix_c_im, matrix_h_re, matrix_h_im, matrix_s_re, matrix_s_im, mo_set, energy_gap, solver_type)
Construct FULL_ALL directly from one complex H(k), S(k), and C(k) channel.
subroutine, public make_preconditioner_complex_full_s_inverse(preconditioner_env, matrix_s_re, matrix_s_im, solver_type)
Construct a complex FULL_S_INVERSE preconditioner.
subroutine, public make_preconditioner_complex_full_single(preconditioner_env, matrix_h_re, matrix_h_im, matrix_s_re, matrix_s_im, mo_set, energy_gap, solver_type)
Construct a complex FULL_SINGLE preconditioner from H(k) and S(k).
subroutine, public make_preconditioner_complex_full_single_inverse(preconditioner_env, matrix_c_re, matrix_c_im, matrix_h_re, matrix_h_im, matrix_s_re, matrix_s_im, energy_gap, solver_type)
Construct a complex FULL_SINGLE_INVERSE preconditioner without discarding Im(H,S,C).
subroutine, public make_preconditioner_complex_full_kinetic(preconditioner_env, matrix_t_re, matrix_t_im, matrix_s_re, matrix_s_im, energy_gap, solver_type)
Construct a complex FULL_KINETIC preconditioner.
collects routines that perform operations directly related to MOs
Definition and initialisation of the mo data type.
Definition qs_mo_types.F:22
an eigen-space solver for the generalised symmetric eigenvalue problem for sparse matrices,...
subroutine, public ot_eigensolver_complex(matrix_h, matrix_h_im, matrix_s, matrix_s_im, matrix_c_fm, matrix_c_fm_im, preconditioner, eps_gradient, iter_max, eigenvalues, silent, ot_settings, matrix_t, matrix_t_im, mo_set)
solve a fixed complex Hermitian generalized eigenproblem by OT
subroutine, public ot_eigensolver(matrix_h, matrix_s, matrix_orthogonal_space_fm, matrix_c_fm, preconditioner, eps_gradient, iter_max, size_ortho_space, silent, ot_settings)
...
orbital transformations
subroutine, public ot_mini(qs_ot_env, matrix_hc, matrix_hc_im, matrix_hc_physical, matrix_hc_physical_im, para_env_inter_kp, gradient_only, gradient_prepared)
...
orbital transformations
Definition qs_ot_types.F:15
subroutine, public qs_ot_init(qs_ot_env)
init matrices, needs c0 and sc0 so that c0*sc0=1
subroutine, public qs_ot_settings_init(settings)
sets default values for the settings type
subroutine, public qs_ot_allocate(qs_ot_env, matrix_s, fm_struct_ref, ortho_k, energy_dimension)
allocates the data in qs_ot_env, for a calculation with fm_struct_ref ortho_k allows for specifying a...
subroutine, public qs_ot_allocate_complex_state(qs_ot_env, matrix_s)
...
subroutine, public qs_ot_destroy(qs_ot_env)
deallocates data
orbital transformations
Definition qs_ot.F:15
subroutine, public qs_ot_get_p(matrix_x, matrix_sx, qs_ot_env)
computes p=x*S*x and the matrix functionals related matrices
Definition qs_ot.F:2594
subroutine, public qs_ot_get_orbitals_ref_complex(matrix_c, matrix_c_im, matrix_s, matrix_s_im, qs_ot_env, qs_ot_env1)
update complex REF k-point orbitals and their S(k)C(k) images
Definition qs_ot.F:1937
subroutine, public qs_ot_prepare_complex_tangent_metric(qs_ot_env, preconditioner_rejected)
Prepare the inverse metric used to project a complex STRICT gradient. An unusable preconditioner is d...
Definition qs_ot.F:3370
subroutine, public qs_ot_get_orbitals(matrix_c, matrix_x, qs_ot_env)
c=(c0*cos(p^0.5)+x*sin(p^0.5)*p^(-0.5)) x rot_mat_u this assumes that x is already ortho to S*C0,...
Definition qs_ot.F:3158
subroutine, public qs_ot_get_orbitals_ref(matrix_c, matrix_s, matrix_x, matrix_sx, matrix_gx_old, matrix_dx, qs_ot_env, qs_ot_env1)
...
Definition qs_ot.F:1832
subroutine, public qs_ot_new_preconditioner(qs_ot_env, preconditioner)
gets ready to use the preconditioner/ or renew the preconditioner only keeps a pointer to the precond...
Definition qs_ot.F:1322
subroutine, public qs_ot_get_orbitals_complex(matrix_c, matrix_c_im, matrix_s, matrix_s_im, qs_ot_env)
update complex K-point orbitals with the finite STRICT transformation
Definition qs_ot.F:3208
Represent a complex full matrix.
represent a full matrix
notice, this variable needs to be copyable, needed for spins as e.g. in qs_ot_scf
Definition qs_ot_types.F:78