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 1060 : 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 530 : 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 530 : TYPE(qs_ot_type), DIMENSION(:), POINTER :: qs_ot_env
136 :
137 530 : CALL timeset(routineN, handle)
138 :
139 530 : output_unit = cp_logger_get_default_io_unit()
140 :
141 530 : IF (PRESENT(silent)) THEN
142 254 : my_silent = silent
143 : ELSE
144 : my_silent = .FALSE.
145 : END IF
146 :
147 530 : NULLIFY (matrix_c) ! fm->dbcsr
148 :
149 530 : CALL cp_fm_get_info(matrix_c_fm, nrow_global=n, ncol_global=k) ! fm->dbcsr
150 530 : ALLOCATE (matrix_c)
151 530 : ref_algorithm = .FALSE.
152 530 : 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 506 : CALL cp_fm_to_dbcsr_row_template(matrix_c, fm_in=matrix_c_fm, template=matrix_h)
160 : END IF
161 :
162 530 : iter_total = 0
163 :
164 : outer_scf: DO
165 :
166 : NULLIFY (qs_ot_env)
167 :
168 712 : NULLIFY (matrix_s_ortho)
169 712 : NULLIFY (matrix_os_ortho)
170 712 : NULLIFY (matrix_buf1_ortho)
171 712 : NULLIFY (matrix_buf2_ortho)
172 712 : NULLIFY (matrix_orthogonal_space)
173 :
174 117480 : ALLOCATE (qs_ot_env(1))
175 1424 : ALLOCATE (matrix_hc(1))
176 712 : NULLIFY (matrix_hc(1)%matrix)
177 712 : CALL dbcsr_init_p(matrix_hc(1)%matrix)
178 :
179 712 : ortho = .FALSE.
180 712 : IF (PRESENT(matrix_orthogonal_space_fm)) ortho = .TRUE.
181 :
182 : ! decide settings
183 712 : IF (PRESENT(ot_settings)) THEN
184 262 : qs_ot_env(1)%settings = ot_settings
185 : ELSE
186 450 : CALL qs_ot_settings_init(qs_ot_env(1)%settings)
187 : ! overwrite defaults
188 450 : qs_ot_env(1)%settings%ds_min = 0.10_dp
189 : END IF
190 :
191 712 : IF (ortho) THEN
192 450 : ALLOCATE (matrix_orthogonal_space)
193 450 : CALL cp_fm_to_dbcsr_row_template(matrix_orthogonal_space, fm_in=matrix_orthogonal_space_fm, template=matrix_h)
194 450 : CALL cp_fm_get_info(matrix_orthogonal_space_fm, ncol_global=ortho_space_k)
195 :
196 450 : IF (PRESENT(size_ortho_space)) ortho_space_k = size_ortho_space
197 450 : ortho_k = ortho_space_k + k
198 : ELSE
199 262 : ortho_k = k
200 : END IF
201 :
202 : ! allocate
203 712 : CALL qs_ot_allocate(qs_ot_env(1), matrix_s, matrix_c_fm%matrix_struct, ortho_k=ortho_k)
204 712 : IF (ref_algorithm) THEN
205 30 : CALL dbcsr_copy(matrix_hc(1)%matrix, qs_ot_env(1)%matrix_x, 'matrix_hc')
206 : ELSE
207 682 : CALL dbcsr_copy(matrix_hc(1)%matrix, matrix_c, 'matrix_hc')
208 : END IF
209 :
210 712 : IF (ortho) THEN
211 : ! construct an initial guess that is orthogonal to matrix_orthogonal_space
212 :
213 450 : CALL dbcsr_init_p(matrix_s_ortho)
214 450 : CALL dbcsr_copy(matrix_s_ortho, matrix_orthogonal_space, name="matrix_s_ortho")
215 :
216 450 : 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 450 : sym=dbcsr_type_no_symmetry)
219 :
220 450 : 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 450 : sym=dbcsr_type_no_symmetry)
223 :
224 450 : 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 450 : sym=dbcsr_type_no_symmetry)
227 :
228 : CALL dbcsr_multiply('N', 'N', 1.0_dp, matrix_s, matrix_orthogonal_space, &
229 450 : 0.0_dp, matrix_s_ortho)
230 : CALL dbcsr_multiply('T', 'N', rone, matrix_s_ortho, matrix_s_ortho, &
231 450 : rzero, matrix_os_ortho)
232 :
233 : CALL cp_dbcsr_cholesky_decompose(matrix_os_ortho, &
234 450 : 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 450 : uplo_to_full=.TRUE.)
238 :
239 : CALL dbcsr_multiply('T', 'N', rone, matrix_s_ortho, matrix_c, &
240 450 : rzero, matrix_buf1_ortho)
241 : CALL dbcsr_multiply('N', 'N', rone, matrix_os_ortho, matrix_buf1_ortho, &
242 450 : rzero, matrix_buf2_ortho)
243 : CALL dbcsr_multiply('N', 'N', -rone, matrix_s_ortho, matrix_buf2_ortho, &
244 450 : rone, matrix_c)
245 :
246 : ! make matrix_c0 an orthogonal basis, matrix_c contains sc0
247 450 : 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 450 : 0.0_dp, matrix_c)
250 :
251 : CALL make_basis_sv(qs_ot_env(1)%matrix_c0, k, matrix_c, &
252 450 : 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 450 : 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 450 : para_env=qs_ot_env(1)%para_env, blacs_env=qs_ot_env(1)%blacs_env)
261 :
262 450 : CALL dbcsr_release_p(matrix_buf1_ortho)
263 450 : CALL dbcsr_release_p(matrix_buf2_ortho)
264 450 : CALL dbcsr_release_p(matrix_os_ortho)
265 450 : CALL dbcsr_release_p(matrix_s_ortho)
266 :
267 : ELSE
268 :
269 : ! set c0,sc0
270 262 : 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 262 : 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 262 : qs_ot_env(1)%para_env, qs_ot_env(1)%blacs_env)
276 : END IF
277 :
278 : ! init
279 712 : CALL qs_ot_init(qs_ot_env(1))
280 712 : energy_only = qs_ot_env(1)%energy_only
281 :
282 682 : SELECT CASE (qs_ot_env(1)%settings%ot_algorithm)
283 : CASE ("TOD")
284 682 : CALL dbcsr_set(qs_ot_env(1)%matrix_x, 0.0_dp)
285 682 : CALL dbcsr_set(qs_ot_env(1)%matrix_sx, 0.0_dp)
286 682 : CALL qs_ot_get_p(qs_ot_env(1)%matrix_x, qs_ot_env(1)%matrix_sx, qs_ot_env(1))
287 682 : 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 712 : CPABORT("OT eigensolver supports ALGORITHM STRICT or IRAC")
294 : END SELECT
295 :
296 : ! if present preconditioner, use it
297 :
298 712 : IF (PRESENT(preconditioner)) THEN
299 712 : IF (ASSOCIATED(preconditioner)) THEN
300 462 : IF (preconditioner_in_use(preconditioner)) THEN
301 462 : 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 12586 : eigensolver_loop: DO
311 :
312 12586 : ieigensolver = ieigensolver + 1
313 12586 : 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 12586 : 0.0_dp, matrix_hc(1)%matrix)
318 12586 : CALL dbcsr_dot(matrix_c, matrix_hc(1)%matrix, energy)
319 12586 : IF (.NOT. energy_only) THEN
320 6684 : CALL dbcsr_scale(matrix_hc(1)%matrix, 2.0_dp)
321 : END IF
322 :
323 12586 : qs_ot_env(1)%etotal = energy
324 12586 : CALL ot_mini(qs_ot_env, matrix_hc)
325 12586 : delta = qs_ot_env(1)%delta
326 12586 : energy_only = qs_ot_env(1)%energy_only
327 :
328 12042 : 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 12042 : 0.0_dp, qs_ot_env(1)%matrix_sx)
332 12042 : CALL qs_ot_get_p(qs_ot_env(1)%matrix_x, qs_ot_env(1)%matrix_sx, qs_ot_env(1))
333 12042 : 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 12586 : 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 12586 : 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 12586 : 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 712 : CALL qs_ot_destroy(qs_ot_env(1))
348 712 : DEALLOCATE (qs_ot_env)
349 712 : CALL dbcsr_release_p(matrix_hc(1)%matrix)
350 712 : DEALLOCATE (matrix_hc)
351 712 : CALL dbcsr_release_p(matrix_orthogonal_space)
352 :
353 712 : IF (delta < eps_gradient) THEN
354 384 : IF ((output_unit > 0) .AND. .NOT. my_silent) THEN
355 : WRITE (UNIT=output_unit, FMT="(T2,A,I0,A)") &
356 144 : "OT| Eigensolver reached convergence in ", iter_total, " iterations"
357 : END IF
358 : EXIT outer_scf
359 : END IF
360 328 : IF (iter_total >= iter_max) THEN
361 146 : IF (output_unit > 0) THEN
362 73 : IF (my_silent) THEN
363 73 : 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 530 : CALL copy_dbcsr_to_fm(matrix_c, matrix_c_fm) ! fm->dbcsr
376 530 : CALL dbcsr_release_p(matrix_c) ! fm->dbcsr
377 :
378 530 : CALL timestop(handle)
379 :
380 530 : 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_cfm complex orbital coefficients
389 : !> \param preconditioner optional complex k-point preconditioner
390 : !> \param eps_gradient requested OT gradient accuracy
391 : !> \param iter_max maximum number of OT iterations
392 : !> \param eigenvalues ...
393 : !> \param silent suppress successful convergence output
394 : !> \param ot_settings OT algorithm and minimizer settings
395 : !> \param matrix_t ...
396 : !> \param matrix_t_im ...
397 : !> \param mo_set ...
398 : ! **************************************************************************************************
399 860 : SUBROUTINE ot_eigensolver_complex(matrix_h, matrix_h_im, matrix_s, matrix_s_im, &
400 : matrix_c_cfm, preconditioner, &
401 172 : eps_gradient, iter_max, eigenvalues, silent, ot_settings, &
402 : matrix_t, matrix_t_im, mo_set)
403 :
404 : TYPE(dbcsr_type), POINTER :: matrix_h, matrix_h_im, matrix_s, &
405 : matrix_s_im
406 : TYPE(cp_cfm_type), INTENT(INOUT) :: matrix_c_cfm
407 : TYPE(preconditioner_type), OPTIONAL, POINTER :: preconditioner
408 : REAL(KIND=dp), INTENT(IN) :: eps_gradient
409 : INTEGER, INTENT(IN) :: iter_max
410 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT), OPTIONAL :: eigenvalues
411 : LOGICAL, INTENT(IN), OPTIONAL :: silent
412 : TYPE(qs_ot_settings_type), INTENT(IN), OPTIONAL :: ot_settings
413 : TYPE(dbcsr_type), OPTIONAL, POINTER :: matrix_t, matrix_t_im
414 : TYPE(mo_set_type), INTENT(IN), OPTIONAL :: mo_set
415 :
416 : CHARACTER(len=*), PARAMETER :: routineN = 'ot_eigensolver_complex'
417 : INTEGER, PARAMETER :: max_iter_inner_loop = 40
418 :
419 : INTEGER :: handle, ieigensolver, iter_total, k, n, &
420 : output_unit
421 : LOGICAL :: energy_only, my_silent, &
422 : preconditioner_rejected
423 : REAL(KIND=dp) :: delta, energy_im, energy_re
424 : TYPE(cp_fm_type) :: matrix_c_fm, matrix_c_fm_im
425 172 : 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 172 : TYPE(qs_ot_type), DIMENSION(:), POINTER :: qs_ot_env
429 :
430 172 : CALL timeset(routineN, handle)
431 172 : output_unit = cp_logger_get_default_io_unit()
432 172 : my_silent = .FALSE.
433 172 : IF (PRESENT(silent)) my_silent = silent
434 172 : NULLIFY (local_preconditioner)
435 :
436 172 : CALL cp_fm_create(matrix_c_fm, matrix_c_cfm%matrix_struct)
437 172 : CALL cp_fm_create(matrix_c_fm_im, matrix_c_cfm%matrix_struct)
438 172 : CALL cp_cfm_to_fm(matrix_c_cfm, matrix_c_fm, matrix_c_fm_im)
439 172 : CALL cp_fm_get_info(matrix_c_fm, nrow_global=n, ncol_global=k)
440 172 : CPASSERT(n > 0 .AND. k > 0)
441 172 : ALLOCATE (matrix_c, matrix_c_im, matrix_tmp)
442 172 : CALL dbcsr_init_p(matrix_c)
443 : CALL cp_dbcsr_m_by_n_from_row_template(matrix_c, template=matrix_h, n=k, &
444 172 : sym=dbcsr_type_no_symmetry)
445 172 : CALL dbcsr_init_p(matrix_c_im)
446 : CALL cp_dbcsr_m_by_n_from_row_template(matrix_c_im, template=matrix_h, n=k, &
447 172 : sym=dbcsr_type_no_symmetry)
448 172 : CALL dbcsr_init_p(matrix_tmp)
449 : CALL cp_dbcsr_m_by_n_from_row_template(matrix_tmp, template=matrix_h, n=k, &
450 172 : sym=dbcsr_type_no_symmetry)
451 172 : CALL copy_fm_to_dbcsr(matrix_c_fm, matrix_c)
452 172 : CALL copy_fm_to_dbcsr(matrix_c_fm_im, matrix_c_im)
453 :
454 172 : iter_total = 0
455 : outer_scf: DO
456 29392 : ALLOCATE (qs_ot_env(1), matrix_hc(1), matrix_hc_im(1))
457 176 : NULLIFY (matrix_hc(1)%matrix, matrix_hc_im(1)%matrix)
458 :
459 176 : IF (PRESENT(ot_settings)) THEN
460 176 : qs_ot_env(1)%settings = ot_settings
461 : ELSE
462 0 : CALL qs_ot_settings_init(qs_ot_env(1)%settings)
463 0 : qs_ot_env(1)%settings%ds_min = 0.10_dp
464 : END IF
465 :
466 176 : CALL qs_ot_allocate(qs_ot_env(1), matrix_s, matrix_c_fm%matrix_struct)
467 176 : CALL qs_ot_allocate_complex_state(qs_ot_env(1), matrix_s)
468 176 : CALL dbcsr_init_p(matrix_hc(1)%matrix)
469 176 : CALL dbcsr_copy(matrix_hc(1)%matrix, qs_ot_env(1)%matrix_x, 'matrix_hc')
470 176 : CALL dbcsr_init_p(matrix_hc_im(1)%matrix)
471 176 : CALL dbcsr_copy(matrix_hc_im(1)%matrix, qs_ot_env(1)%matrix_x_im, 'matrix_hc_im')
472 :
473 176 : CALL dbcsr_copy(qs_ot_env(1)%matrix_c0, matrix_c)
474 176 : CALL dbcsr_copy(qs_ot_env(1)%matrix_c0_im, matrix_c_im)
475 : CALL dbcsr_multiply('N', 'N', 1.0_dp, matrix_s, matrix_c, &
476 176 : 0.0_dp, qs_ot_env(1)%matrix_sc0)
477 : CALL dbcsr_multiply('N', 'N', 1.0_dp, matrix_s_im, matrix_c_im, &
478 176 : 0.0_dp, matrix_tmp)
479 : CALL dbcsr_add(qs_ot_env(1)%matrix_sc0, matrix_tmp, &
480 176 : alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
481 : CALL dbcsr_multiply('N', 'N', 1.0_dp, matrix_s, matrix_c_im, &
482 176 : 0.0_dp, qs_ot_env(1)%matrix_sc0_im)
483 : CALL dbcsr_multiply('N', 'N', 1.0_dp, matrix_s_im, matrix_c, &
484 176 : 0.0_dp, matrix_tmp)
485 : CALL dbcsr_add(qs_ot_env(1)%matrix_sc0_im, matrix_tmp, &
486 176 : alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
487 :
488 176 : CALL qs_ot_init(qs_ot_env(1))
489 176 : energy_only = qs_ot_env(1)%energy_only
490 100 : SELECT CASE (qs_ot_env(1)%settings%ot_algorithm)
491 : CASE ('TOD')
492 100 : CALL dbcsr_set(qs_ot_env(1)%matrix_x, 0.0_dp)
493 100 : CALL dbcsr_set(qs_ot_env(1)%matrix_x_im, 0.0_dp)
494 100 : CALL dbcsr_set(qs_ot_env(1)%matrix_sx, 0.0_dp)
495 100 : CALL dbcsr_set(qs_ot_env(1)%matrix_sx_im, 0.0_dp)
496 : CALL qs_ot_get_orbitals_complex(matrix_c, matrix_c_im, matrix_s, matrix_s_im, &
497 100 : qs_ot_env(1))
498 : CASE ('REF')
499 76 : CALL dbcsr_copy(qs_ot_env(1)%matrix_x, qs_ot_env(1)%matrix_c0)
500 76 : CALL dbcsr_copy(qs_ot_env(1)%matrix_x_im, qs_ot_env(1)%matrix_c0_im)
501 76 : CALL dbcsr_copy(qs_ot_env(1)%matrix_sx, qs_ot_env(1)%matrix_sc0)
502 76 : CALL dbcsr_copy(qs_ot_env(1)%matrix_sx_im, qs_ot_env(1)%matrix_sc0_im)
503 : CALL qs_ot_get_orbitals_ref_complex(matrix_c, matrix_c_im, matrix_s, matrix_s_im, &
504 76 : qs_ot_env(1), qs_ot_env(1))
505 : CASE DEFAULT
506 176 : CPABORT('Complex OT eigensolver supports ALGORITHM STRICT or IRAC')
507 : END SELECT
508 :
509 176 : IF (.NOT. PRESENT(preconditioner) .AND. .NOT. ASSOCIATED(local_preconditioner) .AND. &
510 : qs_ot_env(1)%settings%preconditioner_type /= ot_precond_none) THEN
511 688 : ALLOCATE (local_preconditioner)
512 172 : CALL init_preconditioner(local_preconditioner, qs_ot_env(1)%para_env, qs_ot_env(1)%blacs_env)
513 212 : SELECT CASE (qs_ot_env(1)%settings%preconditioner_type)
514 : CASE (ot_precond_fermi_low_rank)
515 : CALL make_preconditioner_complex_fermi_low_rank( &
516 : local_preconditioner, qs_ot_env(1)%matrix_c0, qs_ot_env(1)%matrix_c0_im, &
517 : matrix_h, matrix_h_im, matrix_s, matrix_s_im, &
518 : qs_ot_env(1)%settings%energy_gap, &
519 : qs_ot_env(1)%settings%fermi_low_rank_max_rank, &
520 40 : qs_ot_env(1)%settings%precond_solver_type)
521 : CASE (ot_precond_full_all, ot_precond_full_all_covariant)
522 92 : IF (qs_ot_env(1)%settings%ot_algorithm == 'TOD') THEN
523 : ! State-selective FULL_ALL models can destabilize the finite STRICT fixed-H chart
524 : ! even when the tangent metric is positive. Retain a well-defined eigensolver with
525 : ! their S^-1 base; IRAC can use the complete state-dependent inverse below.
526 : CALL make_preconditioner_complex_full_s_inverse( &
527 : local_preconditioner, matrix_s, matrix_s_im, &
528 16 : qs_ot_env(1)%settings%precond_solver_type)
529 76 : ELSE IF (qs_ot_env(1)%settings%preconditioner_type == ot_precond_full_all) THEN
530 36 : IF (.NOT. PRESENT(mo_set)) THEN
531 0 : CPABORT('Complex FULL_ALL eigensolver preconditioning requires MO energy labels')
532 : END IF
533 : CALL make_preconditioner_complex_full_all( &
534 : local_preconditioner, qs_ot_env(1)%matrix_c0, qs_ot_env(1)%matrix_c0_im, &
535 : matrix_h, matrix_h_im, matrix_s, matrix_s_im, mo_set, &
536 36 : qs_ot_env(1)%settings%energy_gap, qs_ot_env(1)%settings%precond_solver_type)
537 : ELSE
538 : CALL make_preconditioner_complex_full_all_covariant( &
539 : local_preconditioner, qs_ot_env(1)%matrix_c0, qs_ot_env(1)%matrix_c0_im, &
540 : matrix_h, matrix_h_im, matrix_s, matrix_s_im, &
541 40 : qs_ot_env(1)%settings%energy_gap, qs_ot_env(1)%settings%precond_solver_type)
542 : END IF
543 : CASE (ot_precond_full_single)
544 0 : IF (.NOT. PRESENT(mo_set)) THEN
545 0 : CPABORT('Complex FULL_SINGLE eigensolver preconditioning requires MO energy labels')
546 : END IF
547 : CALL make_preconditioner_complex_full_single( &
548 : local_preconditioner, matrix_h, matrix_h_im, matrix_s, matrix_s_im, mo_set, &
549 0 : qs_ot_env(1)%settings%energy_gap, qs_ot_env(1)%settings%precond_solver_type)
550 : CASE (ot_precond_full_single_inverse)
551 : CALL make_preconditioner_complex_full_single_inverse( &
552 : local_preconditioner, qs_ot_env(1)%matrix_c0, qs_ot_env(1)%matrix_c0_im, &
553 : matrix_h, matrix_h_im, matrix_s, matrix_s_im, &
554 40 : qs_ot_env(1)%settings%energy_gap, qs_ot_env(1)%settings%precond_solver_type)
555 : CASE (ot_precond_full_kinetic)
556 0 : IF (.NOT. PRESENT(matrix_t) .OR. .NOT. PRESENT(matrix_t_im)) THEN
557 0 : CPABORT('Complex FULL_KINETIC eigensolver preconditioning requires T(k)')
558 : END IF
559 0 : IF (.NOT. ASSOCIATED(matrix_t) .OR. .NOT. ASSOCIATED(matrix_t_im)) THEN
560 0 : CPABORT('Complex FULL_KINETIC eigensolver preconditioning requires T(k)')
561 : END IF
562 : CALL make_preconditioner_complex_full_kinetic( &
563 : local_preconditioner, matrix_t, matrix_t_im, matrix_s, matrix_s_im, &
564 0 : qs_ot_env(1)%settings%energy_gap, qs_ot_env(1)%settings%precond_solver_type)
565 : CASE (ot_precond_s_inverse)
566 : CALL make_preconditioner_complex_full_s_inverse( &
567 : local_preconditioner, matrix_s, matrix_s_im, &
568 0 : qs_ot_env(1)%settings%precond_solver_type)
569 : CASE DEFAULT
570 172 : CPABORT('Unsupported complex K-point OT eigensolver preconditioner')
571 : END SELECT
572 : END IF
573 :
574 176 : IF (PRESENT(preconditioner)) THEN
575 0 : IF (ASSOCIATED(preconditioner)) THEN
576 0 : IF (preconditioner_in_use(preconditioner)) THEN
577 0 : CALL qs_ot_new_preconditioner(qs_ot_env(1), preconditioner)
578 : END IF
579 : END IF
580 176 : ELSE IF (ASSOCIATED(local_preconditioner)) THEN
581 176 : CALL qs_ot_new_preconditioner(qs_ot_env(1), local_preconditioner)
582 : END IF
583 :
584 176 : IF (qs_ot_env(1)%settings%ot_algorithm == 'TOD') THEN
585 100 : CALL qs_ot_prepare_complex_tangent_metric(qs_ot_env(1), preconditioner_rejected)
586 100 : IF (preconditioner_rejected .AND. ASSOCIATED(local_preconditioner)) THEN
587 0 : CALL destroy_preconditioner(local_preconditioner)
588 : CALL init_preconditioner(local_preconditioner, &
589 0 : qs_ot_env(1)%para_env, qs_ot_env(1)%blacs_env)
590 : CALL make_preconditioner_complex_full_s_inverse( &
591 : local_preconditioner, matrix_s, matrix_s_im, &
592 0 : qs_ot_env(1)%settings%precond_solver_type)
593 0 : CALL qs_ot_new_preconditioner(qs_ot_env(1), local_preconditioner)
594 0 : CALL qs_ot_prepare_complex_tangent_metric(qs_ot_env(1), preconditioner_rejected)
595 0 : CPASSERT(.NOT. preconditioner_rejected)
596 : END IF
597 : END IF
598 :
599 : ieigensolver = 0
600 1144 : eigensolver_loop: DO
601 572 : ieigensolver = ieigensolver + 1
602 572 : iter_total = iter_total + 1
603 :
604 : ! H*C = (Hre*Cre-Him*Cim) + i*(Hre*Cim+Him*Cre).
605 : CALL dbcsr_multiply('N', 'N', 1.0_dp, matrix_h, matrix_c, &
606 572 : 0.0_dp, matrix_hc(1)%matrix)
607 : CALL dbcsr_multiply('N', 'N', 1.0_dp, matrix_h_im, matrix_c_im, &
608 572 : 0.0_dp, matrix_tmp)
609 : CALL dbcsr_add(matrix_hc(1)%matrix, matrix_tmp, &
610 572 : alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
611 : CALL dbcsr_multiply('N', 'N', 1.0_dp, matrix_h, matrix_c_im, &
612 572 : 0.0_dp, matrix_hc_im(1)%matrix)
613 : CALL dbcsr_multiply('N', 'N', 1.0_dp, matrix_h_im, matrix_c, &
614 572 : 0.0_dp, matrix_tmp)
615 : CALL dbcsr_add(matrix_hc_im(1)%matrix, matrix_tmp, &
616 572 : alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
617 :
618 572 : CALL dbcsr_dot(matrix_c, matrix_hc(1)%matrix, energy_re)
619 572 : CALL dbcsr_dot(matrix_c_im, matrix_hc_im(1)%matrix, energy_im)
620 572 : qs_ot_env(1)%etotal = energy_re + energy_im
621 572 : IF (.NOT. energy_only) THEN
622 572 : CALL dbcsr_scale(matrix_hc(1)%matrix, 2.0_dp)
623 572 : CALL dbcsr_scale(matrix_hc_im(1)%matrix, 2.0_dp)
624 : END IF
625 572 : CALL ot_mini(qs_ot_env, matrix_hc, matrix_hc_im=matrix_hc_im)
626 572 : delta = qs_ot_env(1)%delta
627 572 : energy_only = qs_ot_env(1)%energy_only
628 :
629 426 : SELECT CASE (qs_ot_env(1)%settings%ot_algorithm)
630 : CASE ('TOD')
631 : CALL qs_ot_get_orbitals_complex(matrix_c, matrix_c_im, matrix_s, matrix_s_im, &
632 426 : qs_ot_env(1))
633 : CASE ('REF')
634 : CALL qs_ot_get_orbitals_ref_complex(matrix_c, matrix_c_im, matrix_s, matrix_s_im, &
635 572 : qs_ot_env(1), qs_ot_env(1))
636 : END SELECT
637 :
638 572 : IF (delta < eps_gradient .OR. ieigensolver >= max_iter_inner_loop) EXIT eigensolver_loop
639 572 : IF (iter_total >= iter_max .AND. qs_ot_env(1)%OT_METHOD_FULL /= 'OT LS') EXIT eigensolver_loop
640 : END DO eigensolver_loop
641 :
642 176 : CALL qs_ot_destroy(qs_ot_env(1))
643 176 : DEALLOCATE (qs_ot_env)
644 176 : CALL dbcsr_release_p(matrix_hc(1)%matrix)
645 176 : CALL dbcsr_release_p(matrix_hc_im(1)%matrix)
646 176 : DEALLOCATE (matrix_hc, matrix_hc_im)
647 :
648 176 : IF (delta < eps_gradient) THEN
649 172 : IF (output_unit > 0 .AND. .NOT. my_silent) THEN
650 : WRITE (UNIT=output_unit, FMT='(T2,A,I0,A)') &
651 0 : 'OT| Complex eigensolver reached convergence in ', iter_total, ' iterations'
652 : END IF
653 : EXIT outer_scf
654 : END IF
655 4 : IF (iter_total >= iter_max) THEN
656 0 : IF (output_unit > 0) THEN
657 : WRITE (output_unit, '(A,T60,E20.10)') &
658 0 : ' WARNING complex OT eigensolver did not converge: current gradient', delta
659 : END IF
660 : EXIT outer_scf
661 : END IF
662 : END DO outer_scf
663 :
664 172 : CALL copy_dbcsr_to_fm(matrix_c, matrix_c_fm)
665 172 : CALL copy_dbcsr_to_fm(matrix_c_im, matrix_c_fm_im)
666 172 : IF (PRESENT(eigenvalues)) THEN
667 : CALL canonicalize_complex_subspace(matrix_h, matrix_h_im, matrix_c_fm, &
668 172 : matrix_c_fm_im, eigenvalues)
669 : END IF
670 172 : CALL cp_fm_to_cfm(matrix_c_fm, matrix_c_fm_im, matrix_c_cfm)
671 172 : CALL cp_fm_release(matrix_c_fm)
672 172 : CALL cp_fm_release(matrix_c_fm_im)
673 172 : CALL dbcsr_release_p(matrix_c)
674 172 : CALL dbcsr_release_p(matrix_c_im)
675 172 : CALL dbcsr_release_p(matrix_tmp)
676 172 : IF (ASSOCIATED(local_preconditioner)) THEN
677 172 : CALL destroy_preconditioner(local_preconditioner)
678 172 : DEALLOCATE (local_preconditioner)
679 : END IF
680 172 : CALL timestop(handle)
681 :
682 172 : END SUBROUTINE ot_eigensolver_complex
683 :
684 : ! **************************************************************************************************
685 : !> \brief diagonalize C^H H C and rotate a complex orthonormal orbital subspace
686 : !> \param matrix_h ...
687 : !> \param matrix_h_im ...
688 : !> \param coeff_re ...
689 : !> \param coeff_im ...
690 : !> \param eigenvalues ...
691 : ! **************************************************************************************************
692 172 : SUBROUTINE canonicalize_complex_subspace(matrix_h, matrix_h_im, coeff_re, coeff_im, eigenvalues)
693 : TYPE(dbcsr_type), POINTER :: matrix_h, matrix_h_im
694 : TYPE(cp_fm_type), INTENT(INOUT) :: coeff_re, coeff_im
695 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: eigenvalues
696 :
697 : INTEGER :: nao, nmo
698 : TYPE(cp_cfm_type) :: coeff, eigenvectors, hblock, hc, rotated
699 : TYPE(cp_fm_type) :: hc_im, hc_re
700 :
701 172 : CALL cp_fm_get_info(coeff_re, nrow_global=nao, ncol_global=nmo)
702 172 : CPASSERT(SIZE(eigenvalues) >= nmo)
703 172 : CALL cp_fm_create(hc_re, coeff_re%matrix_struct)
704 172 : CALL cp_fm_create(hc_im, coeff_re%matrix_struct)
705 172 : CALL cp_dbcsr_sm_fm_multiply(matrix_h, coeff_re, hc_re, nmo)
706 : CALL cp_dbcsr_sm_fm_multiply(matrix_h_im, coeff_im, hc_re, nmo, &
707 172 : alpha=-1.0_dp, beta=1.0_dp)
708 172 : CALL cp_dbcsr_sm_fm_multiply(matrix_h, coeff_im, hc_im, nmo)
709 : CALL cp_dbcsr_sm_fm_multiply(matrix_h_im, coeff_re, hc_im, nmo, &
710 172 : alpha=1.0_dp, beta=1.0_dp)
711 :
712 172 : CALL cp_cfm_create(coeff, coeff_re%matrix_struct)
713 172 : CALL cp_cfm_create(hc, coeff_re%matrix_struct)
714 172 : CALL cp_cfm_create(rotated, coeff_re%matrix_struct)
715 172 : CALL cp_cfm_create(hblock, coeff_re%matrix_struct, nrow=nmo, ncol=nmo)
716 172 : CALL cp_cfm_create(eigenvectors, coeff_re%matrix_struct, nrow=nmo, ncol=nmo)
717 172 : CALL cp_fm_to_cfm(coeff_re, coeff_im, coeff)
718 172 : CALL cp_fm_to_cfm(hc_re, hc_im, hc)
719 172 : CALL cp_cfm_gemm('C', 'N', nmo, nmo, nao, z_one, coeff, hc, z_zero, hblock)
720 172 : CALL cp_cfm_heevd(hblock, eigenvectors, eigenvalues(1:nmo))
721 172 : CALL cp_cfm_gemm('N', 'N', nao, nmo, nmo, z_one, coeff, eigenvectors, z_zero, rotated)
722 172 : CALL cp_cfm_to_fm(rotated, coeff_re, coeff_im)
723 :
724 172 : CALL cp_cfm_release(eigenvectors)
725 172 : CALL cp_cfm_release(hblock)
726 172 : CALL cp_cfm_release(rotated)
727 172 : CALL cp_cfm_release(hc)
728 172 : CALL cp_cfm_release(coeff)
729 172 : CALL cp_fm_release(hc_im)
730 172 : CALL cp_fm_release(hc_re)
731 172 : END SUBROUTINE canonicalize_complex_subspace
732 :
733 : END MODULE qs_ot_eigensolver
|
