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_full_all,&
40 : ot_precond_full_kinetic,&
41 : ot_precond_full_single,&
42 : ot_precond_full_single_inverse,&
43 : ot_precond_none,&
44 : ot_precond_s_inverse
45 : USE kinds, ONLY: dp
46 : USE mathconstants, ONLY: z_one,&
47 : z_zero
48 : USE preconditioner, ONLY: make_preconditioner_complex_full_all,&
49 : make_preconditioner_complex_full_kinetic,&
50 : make_preconditioner_complex_full_s_inverse,&
51 : make_preconditioner_complex_full_single,&
52 : make_preconditioner_complex_full_single_inverse
53 : USE preconditioner_types, ONLY: destroy_preconditioner,&
54 : init_preconditioner,&
55 : preconditioner_in_use,&
56 : preconditioner_type
57 : USE qs_mo_methods, ONLY: make_basis_sv
58 : USE qs_mo_types, ONLY: mo_set_type
59 : USE qs_ot, ONLY: qs_ot_get_orbitals,&
60 : qs_ot_get_orbitals_complex,&
61 : qs_ot_get_orbitals_ref,&
62 : qs_ot_get_orbitals_ref_complex,&
63 : qs_ot_get_p,&
64 : qs_ot_new_preconditioner,&
65 : qs_ot_prepare_complex_tangent_metric
66 : USE qs_ot_minimizer, ONLY: ot_mini
67 : USE qs_ot_types, ONLY: qs_ot_allocate,&
68 : qs_ot_allocate_complex_state,&
69 : qs_ot_destroy,&
70 : qs_ot_init,&
71 : qs_ot_settings_init,&
72 : qs_ot_settings_type,&
73 : qs_ot_type
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 :
85 : PUBLIC :: ot_eigensolver, ot_eigensolver_complex
86 :
87 : CONTAINS
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 1008 : 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 504 : 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 504 : TYPE(qs_ot_type), DIMENSION(:), POINTER :: qs_ot_env
132 :
133 504 : CALL timeset(routineN, handle)
134 :
135 504 : output_unit = cp_logger_get_default_io_unit()
136 :
137 504 : IF (PRESENT(silent)) THEN
138 254 : my_silent = silent
139 : ELSE
140 : my_silent = .FALSE.
141 : END IF
142 :
143 504 : NULLIFY (matrix_c) ! fm->dbcsr
144 :
145 504 : CALL cp_fm_get_info(matrix_c_fm, nrow_global=n, ncol_global=k) ! fm->dbcsr
146 504 : ALLOCATE (matrix_c)
147 504 : ref_algorithm = .FALSE.
148 504 : IF (PRESENT(ot_settings)) ref_algorithm = ot_settings%ot_algorithm == "REF"
149 250 : IF (ref_algorithm) THEN
150 24 : CALL dbcsr_init_p(matrix_c)
151 : CALL cp_dbcsr_m_by_n_from_row_template(matrix_c, template=matrix_h, n=k, &
152 24 : sym=dbcsr_type_no_symmetry)
153 24 : CALL copy_fm_to_dbcsr(matrix_c_fm, matrix_c)
154 : ELSE
155 480 : CALL cp_fm_to_dbcsr_row_template(matrix_c, fm_in=matrix_c_fm, template=matrix_h)
156 : END IF
157 :
158 504 : iter_total = 0
159 :
160 : outer_scf: DO
161 :
162 : NULLIFY (qs_ot_env)
163 :
164 674 : NULLIFY (matrix_s_ortho)
165 674 : NULLIFY (matrix_os_ortho)
166 674 : NULLIFY (matrix_buf1_ortho)
167 674 : NULLIFY (matrix_buf2_ortho)
168 674 : NULLIFY (matrix_orthogonal_space)
169 :
170 108514 : ALLOCATE (qs_ot_env(1))
171 1348 : ALLOCATE (matrix_hc(1))
172 674 : NULLIFY (matrix_hc(1)%matrix)
173 674 : CALL dbcsr_init_p(matrix_hc(1)%matrix)
174 :
175 674 : ortho = .FALSE.
176 674 : IF (PRESENT(matrix_orthogonal_space_fm)) ortho = .TRUE.
177 :
178 : ! decide settings
179 674 : IF (PRESENT(ot_settings)) THEN
180 262 : qs_ot_env(1)%settings = ot_settings
181 : ELSE
182 412 : CALL qs_ot_settings_init(qs_ot_env(1)%settings)
183 : ! overwrite defaults
184 412 : qs_ot_env(1)%settings%ds_min = 0.10_dp
185 : END IF
186 :
187 674 : IF (ortho) THEN
188 412 : ALLOCATE (matrix_orthogonal_space)
189 412 : CALL cp_fm_to_dbcsr_row_template(matrix_orthogonal_space, fm_in=matrix_orthogonal_space_fm, template=matrix_h)
190 412 : CALL cp_fm_get_info(matrix_orthogonal_space_fm, ncol_global=ortho_space_k)
191 :
192 412 : IF (PRESENT(size_ortho_space)) ortho_space_k = size_ortho_space
193 412 : ortho_k = ortho_space_k + k
194 : ELSE
195 262 : ortho_k = k
196 : END IF
197 :
198 : ! allocate
199 674 : CALL qs_ot_allocate(qs_ot_env(1), matrix_s, matrix_c_fm%matrix_struct, ortho_k=ortho_k)
200 674 : IF (ref_algorithm) THEN
201 30 : CALL dbcsr_copy(matrix_hc(1)%matrix, qs_ot_env(1)%matrix_x, 'matrix_hc')
202 : ELSE
203 644 : CALL dbcsr_copy(matrix_hc(1)%matrix, matrix_c, 'matrix_hc')
204 : END IF
205 :
206 674 : IF (ortho) THEN
207 : ! construct an initial guess that is orthogonal to matrix_orthogonal_space
208 :
209 412 : CALL dbcsr_init_p(matrix_s_ortho)
210 412 : CALL dbcsr_copy(matrix_s_ortho, matrix_orthogonal_space, name="matrix_s_ortho")
211 :
212 412 : 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 412 : sym=dbcsr_type_no_symmetry)
215 :
216 412 : 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 412 : sym=dbcsr_type_no_symmetry)
219 :
220 412 : 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 412 : sym=dbcsr_type_no_symmetry)
223 :
224 : CALL dbcsr_multiply('N', 'N', 1.0_dp, matrix_s, matrix_orthogonal_space, &
225 412 : 0.0_dp, matrix_s_ortho)
226 : CALL dbcsr_multiply('T', 'N', rone, matrix_s_ortho, matrix_s_ortho, &
227 412 : rzero, matrix_os_ortho)
228 :
229 : CALL cp_dbcsr_cholesky_decompose(matrix_os_ortho, &
230 412 : 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 412 : uplo_to_full=.TRUE.)
234 :
235 : CALL dbcsr_multiply('T', 'N', rone, matrix_s_ortho, matrix_c, &
236 412 : rzero, matrix_buf1_ortho)
237 : CALL dbcsr_multiply('N', 'N', rone, matrix_os_ortho, matrix_buf1_ortho, &
238 412 : rzero, matrix_buf2_ortho)
239 : CALL dbcsr_multiply('N', 'N', -rone, matrix_s_ortho, matrix_buf2_ortho, &
240 412 : rone, matrix_c)
241 :
242 : ! make matrix_c0 an orthogonal basis, matrix_c contains sc0
243 412 : 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 412 : 0.0_dp, matrix_c)
246 :
247 : CALL make_basis_sv(qs_ot_env(1)%matrix_c0, k, matrix_c, &
248 412 : 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 412 : 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 412 : para_env=qs_ot_env(1)%para_env, blacs_env=qs_ot_env(1)%blacs_env)
257 :
258 412 : CALL dbcsr_release_p(matrix_buf1_ortho)
259 412 : CALL dbcsr_release_p(matrix_buf2_ortho)
260 412 : CALL dbcsr_release_p(matrix_os_ortho)
261 412 : CALL dbcsr_release_p(matrix_s_ortho)
262 :
263 : ELSE
264 :
265 : ! set c0,sc0
266 262 : 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 262 : 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 262 : qs_ot_env(1)%para_env, qs_ot_env(1)%blacs_env)
272 : END IF
273 :
274 : ! init
275 674 : CALL qs_ot_init(qs_ot_env(1))
276 674 : energy_only = qs_ot_env(1)%energy_only
277 :
278 644 : SELECT CASE (qs_ot_env(1)%settings%ot_algorithm)
279 : CASE ("TOD")
280 644 : CALL dbcsr_set(qs_ot_env(1)%matrix_x, 0.0_dp)
281 644 : CALL dbcsr_set(qs_ot_env(1)%matrix_sx, 0.0_dp)
282 644 : CALL qs_ot_get_p(qs_ot_env(1)%matrix_x, qs_ot_env(1)%matrix_sx, qs_ot_env(1))
283 644 : CALL qs_ot_get_orbitals(matrix_c, qs_ot_env(1)%matrix_x, qs_ot_env(1))
284 : CASE ("REF")
285 30 : CALL dbcsr_copy(qs_ot_env(1)%matrix_x, qs_ot_env(1)%matrix_c0)
286 30 : CALL dbcsr_copy(qs_ot_env(1)%matrix_sx, qs_ot_env(1)%matrix_sc0)
287 30 : CALL dbcsr_copy(matrix_c, qs_ot_env(1)%matrix_c0)
288 : CASE DEFAULT
289 674 : CPABORT("OT eigensolver supports ALGORITHM STRICT or IRAC")
290 : END SELECT
291 :
292 : ! if present preconditioner, use it
293 :
294 674 : IF (PRESENT(preconditioner)) THEN
295 674 : IF (ASSOCIATED(preconditioner)) THEN
296 446 : IF (preconditioner_in_use(preconditioner)) THEN
297 446 : 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 11886 : eigensolver_loop: DO
307 :
308 11886 : ieigensolver = ieigensolver + 1
309 11886 : 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 11886 : 0.0_dp, matrix_hc(1)%matrix)
314 11886 : CALL dbcsr_dot(matrix_c, matrix_hc(1)%matrix, energy)
315 11886 : IF (.NOT. energy_only) THEN
316 6322 : CALL dbcsr_scale(matrix_hc(1)%matrix, 2.0_dp)
317 : END IF
318 :
319 11886 : qs_ot_env(1)%etotal = energy
320 11886 : CALL ot_mini(qs_ot_env, matrix_hc)
321 11886 : delta = qs_ot_env(1)%delta
322 11886 : energy_only = qs_ot_env(1)%energy_only
323 :
324 11346 : 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 11346 : 0.0_dp, qs_ot_env(1)%matrix_sx)
328 11346 : CALL qs_ot_get_p(qs_ot_env(1)%matrix_x, qs_ot_env(1)%matrix_sx, qs_ot_env(1))
329 11346 : 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 11886 : 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 11886 : 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 11886 : 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 674 : CALL qs_ot_destroy(qs_ot_env(1))
344 674 : DEALLOCATE (qs_ot_env)
345 674 : CALL dbcsr_release_p(matrix_hc(1)%matrix)
346 674 : DEALLOCATE (matrix_hc)
347 674 : CALL dbcsr_release_p(matrix_orthogonal_space)
348 :
349 674 : IF (delta < eps_gradient) THEN
350 354 : IF ((output_unit > 0) .AND. .NOT. my_silent) THEN
351 : WRITE (UNIT=output_unit, FMT="(T2,A,I0,A)") &
352 131 : "OT| Eigensolver reached convergence in ", iter_total, " iterations"
353 : END IF
354 : EXIT outer_scf
355 : END IF
356 320 : IF (iter_total >= iter_max) THEN
357 150 : IF (output_unit > 0) THEN
358 75 : IF (my_silent) THEN
359 75 : WRITE (output_unit, "(A,T60,E20.10)") " WARNING OT eigensolver did not converge: current gradient", delta
360 : ELSE
361 0 : WRITE (output_unit, *) "WARNING : did not converge in ot_eigensolver"
362 0 : WRITE (output_unit, *) "number of iterations ", iter_total, " exceeded maximum"
363 0 : 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 504 : CALL copy_dbcsr_to_fm(matrix_c, matrix_c_fm) ! fm->dbcsr
372 504 : CALL dbcsr_release_p(matrix_c) ! fm->dbcsr
373 :
374 504 : CALL timestop(handle)
375 :
376 504 : 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 234 : SUBROUTINE ot_eigensolver_complex(matrix_h, matrix_h_im, matrix_s, matrix_s_im, &
397 : matrix_c_fm, matrix_c_fm_im, preconditioner, &
398 78 : 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 78 : 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 78 : TYPE(qs_ot_type), DIMENSION(:), POINTER :: qs_ot_env
425 :
426 78 : CALL timeset(routineN, handle)
427 78 : output_unit = cp_logger_get_default_io_unit()
428 78 : my_silent = .FALSE.
429 78 : IF (PRESENT(silent)) my_silent = silent
430 78 : NULLIFY (local_preconditioner)
431 :
432 78 : CALL cp_fm_get_info(matrix_c_fm, nrow_global=n, ncol_global=k)
433 78 : CPASSERT(n > 0 .AND. k > 0)
434 78 : ALLOCATE (matrix_c, matrix_c_im, matrix_tmp)
435 78 : CALL dbcsr_init_p(matrix_c)
436 : CALL cp_dbcsr_m_by_n_from_row_template(matrix_c, template=matrix_h, n=k, &
437 78 : sym=dbcsr_type_no_symmetry)
438 78 : 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 78 : sym=dbcsr_type_no_symmetry)
441 78 : CALL dbcsr_init_p(matrix_tmp)
442 : CALL cp_dbcsr_m_by_n_from_row_template(matrix_tmp, template=matrix_h, n=k, &
443 78 : sym=dbcsr_type_no_symmetry)
444 78 : CALL copy_fm_to_dbcsr(matrix_c_fm, matrix_c)
445 78 : CALL copy_fm_to_dbcsr(matrix_c_fm_im, matrix_c_im)
446 :
447 78 : iter_total = 0
448 : outer_scf: DO
449 14181 : ALLOCATE (qs_ot_env(1), matrix_hc(1), matrix_hc_im(1))
450 87 : NULLIFY (matrix_hc(1)%matrix, matrix_hc_im(1)%matrix)
451 :
452 87 : IF (PRESENT(ot_settings)) THEN
453 87 : qs_ot_env(1)%settings = ot_settings
454 : ELSE
455 0 : CALL qs_ot_settings_init(qs_ot_env(1)%settings)
456 0 : qs_ot_env(1)%settings%ds_min = 0.10_dp
457 : END IF
458 :
459 87 : CALL qs_ot_allocate(qs_ot_env(1), matrix_s, matrix_c_fm%matrix_struct)
460 87 : CALL qs_ot_allocate_complex_state(qs_ot_env(1), matrix_s)
461 87 : CALL dbcsr_init_p(matrix_hc(1)%matrix)
462 87 : CALL dbcsr_copy(matrix_hc(1)%matrix, qs_ot_env(1)%matrix_x, 'matrix_hc')
463 87 : CALL dbcsr_init_p(matrix_hc_im(1)%matrix)
464 87 : CALL dbcsr_copy(matrix_hc_im(1)%matrix, qs_ot_env(1)%matrix_x_im, 'matrix_hc_im')
465 :
466 87 : CALL dbcsr_copy(qs_ot_env(1)%matrix_c0, matrix_c)
467 87 : 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 87 : 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 87 : 0.0_dp, matrix_tmp)
472 : CALL dbcsr_add(qs_ot_env(1)%matrix_sc0, matrix_tmp, &
473 87 : 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 87 : 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 87 : 0.0_dp, matrix_tmp)
478 : CALL dbcsr_add(qs_ot_env(1)%matrix_sc0_im, matrix_tmp, &
479 87 : alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
480 :
481 87 : CALL qs_ot_init(qs_ot_env(1))
482 87 : energy_only = qs_ot_env(1)%energy_only
483 60 : SELECT CASE (qs_ot_env(1)%settings%ot_algorithm)
484 : CASE ('TOD')
485 60 : CALL dbcsr_set(qs_ot_env(1)%matrix_x, 0.0_dp)
486 60 : CALL dbcsr_set(qs_ot_env(1)%matrix_x_im, 0.0_dp)
487 60 : CALL dbcsr_set(qs_ot_env(1)%matrix_sx, 0.0_dp)
488 60 : 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 60 : qs_ot_env(1))
491 : CASE ('REF')
492 27 : CALL dbcsr_copy(qs_ot_env(1)%matrix_x, qs_ot_env(1)%matrix_c0)
493 27 : CALL dbcsr_copy(qs_ot_env(1)%matrix_x_im, qs_ot_env(1)%matrix_c0_im)
494 27 : CALL dbcsr_copy(qs_ot_env(1)%matrix_sx, qs_ot_env(1)%matrix_sc0)
495 27 : 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 27 : qs_ot_env(1), qs_ot_env(1))
498 : CASE DEFAULT
499 87 : CPABORT('Complex OT eigensolver supports ALGORITHM STRICT or IRAC')
500 : END SELECT
501 :
502 87 : IF (.NOT. PRESENT(preconditioner) .AND. .NOT. ASSOCIATED(local_preconditioner) .AND. &
503 : qs_ot_env(1)%settings%preconditioner_type /= ot_precond_none) THEN
504 78 : ALLOCATE (local_preconditioner)
505 78 : CALL init_preconditioner(local_preconditioner, qs_ot_env(1)%para_env, qs_ot_env(1)%blacs_env)
506 94 : SELECT CASE (qs_ot_env(1)%settings%preconditioner_type)
507 : CASE (ot_precond_full_all)
508 16 : 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.
511 : CALL make_preconditioner_complex_full_s_inverse( &
512 : local_preconditioner, matrix_s, matrix_s_im, &
513 16 : qs_ot_env(1)%settings%precond_solver_type)
514 : ELSE
515 0 : IF (.NOT. PRESENT(mo_set)) THEN
516 0 : CPABORT('Complex FULL_ALL eigensolver preconditioning requires MO energy labels')
517 : END IF
518 : CALL make_preconditioner_complex_full_all( &
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 0 : qs_ot_env(1)%settings%energy_gap, qs_ot_env(1)%settings%precond_solver_type)
522 : END IF
523 : CASE (ot_precond_full_single)
524 22 : IF (.NOT. PRESENT(mo_set)) THEN
525 0 : CPABORT('Complex FULL_SINGLE eigensolver preconditioning requires MO energy labels')
526 : END IF
527 : CALL make_preconditioner_complex_full_single( &
528 : local_preconditioner, matrix_h, matrix_h_im, matrix_s, matrix_s_im, mo_set, &
529 22 : qs_ot_env(1)%settings%energy_gap, qs_ot_env(1)%settings%precond_solver_type)
530 : CASE (ot_precond_full_single_inverse)
531 : CALL make_preconditioner_complex_full_single_inverse( &
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 40 : qs_ot_env(1)%settings%energy_gap, qs_ot_env(1)%settings%precond_solver_type)
535 : CASE (ot_precond_full_kinetic)
536 0 : IF (.NOT. PRESENT(matrix_t) .OR. .NOT. PRESENT(matrix_t_im)) THEN
537 0 : CPABORT('Complex FULL_KINETIC eigensolver preconditioning requires T(k)')
538 : END IF
539 0 : IF (.NOT. ASSOCIATED(matrix_t) .OR. .NOT. ASSOCIATED(matrix_t_im)) THEN
540 0 : CPABORT('Complex FULL_KINETIC eigensolver preconditioning requires T(k)')
541 : END IF
542 : CALL make_preconditioner_complex_full_kinetic( &
543 : local_preconditioner, matrix_t, matrix_t_im, matrix_s, matrix_s_im, &
544 0 : qs_ot_env(1)%settings%energy_gap, qs_ot_env(1)%settings%precond_solver_type)
545 : CASE (ot_precond_s_inverse)
546 : CALL make_preconditioner_complex_full_s_inverse( &
547 : local_preconditioner, matrix_s, matrix_s_im, &
548 0 : qs_ot_env(1)%settings%precond_solver_type)
549 : CASE DEFAULT
550 78 : CPABORT('Unsupported complex K-point OT eigensolver preconditioner')
551 : END SELECT
552 : END IF
553 :
554 87 : IF (PRESENT(preconditioner)) THEN
555 0 : IF (ASSOCIATED(preconditioner)) THEN
556 0 : IF (preconditioner_in_use(preconditioner)) THEN
557 0 : CALL qs_ot_new_preconditioner(qs_ot_env(1), preconditioner)
558 : END IF
559 : END IF
560 87 : ELSE IF (ASSOCIATED(local_preconditioner)) THEN
561 87 : CALL qs_ot_new_preconditioner(qs_ot_env(1), local_preconditioner)
562 : END IF
563 :
564 87 : IF (qs_ot_env(1)%settings%ot_algorithm == 'TOD') THEN
565 60 : CALL qs_ot_prepare_complex_tangent_metric(qs_ot_env(1), preconditioner_rejected)
566 60 : IF (preconditioner_rejected .AND. ASSOCIATED(local_preconditioner)) THEN
567 0 : CALL destroy_preconditioner(local_preconditioner)
568 : CALL init_preconditioner(local_preconditioner, &
569 0 : qs_ot_env(1)%para_env, qs_ot_env(1)%blacs_env)
570 : CALL make_preconditioner_complex_full_s_inverse( &
571 : local_preconditioner, matrix_s, matrix_s_im, &
572 0 : qs_ot_env(1)%settings%precond_solver_type)
573 0 : CALL qs_ot_new_preconditioner(qs_ot_env(1), local_preconditioner)
574 0 : CALL qs_ot_prepare_complex_tangent_metric(qs_ot_env(1), preconditioner_rejected)
575 0 : CPASSERT(.NOT. preconditioner_rejected)
576 : END IF
577 : END IF
578 :
579 : ieigensolver = 0
580 1658 : eigensolver_loop: DO
581 829 : ieigensolver = ieigensolver + 1
582 829 : 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 829 : 0.0_dp, matrix_hc(1)%matrix)
587 : CALL dbcsr_multiply('N', 'N', 1.0_dp, matrix_h_im, matrix_c_im, &
588 829 : 0.0_dp, matrix_tmp)
589 : CALL dbcsr_add(matrix_hc(1)%matrix, matrix_tmp, &
590 829 : 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 829 : 0.0_dp, matrix_hc_im(1)%matrix)
593 : CALL dbcsr_multiply('N', 'N', 1.0_dp, matrix_h_im, matrix_c, &
594 829 : 0.0_dp, matrix_tmp)
595 : CALL dbcsr_add(matrix_hc_im(1)%matrix, matrix_tmp, &
596 829 : alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
597 :
598 829 : CALL dbcsr_dot(matrix_c, matrix_hc(1)%matrix, energy_re)
599 829 : CALL dbcsr_dot(matrix_c_im, matrix_hc_im(1)%matrix, energy_im)
600 829 : qs_ot_env(1)%etotal = energy_re + energy_im
601 829 : IF (.NOT. energy_only) THEN
602 829 : CALL dbcsr_scale(matrix_hc(1)%matrix, 2.0_dp)
603 829 : CALL dbcsr_scale(matrix_hc_im(1)%matrix, 2.0_dp)
604 : END IF
605 829 : CALL ot_mini(qs_ot_env, matrix_hc, matrix_hc_im=matrix_hc_im)
606 829 : delta = qs_ot_env(1)%delta
607 829 : energy_only = qs_ot_env(1)%energy_only
608 :
609 386 : 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 386 : 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 829 : qs_ot_env(1), qs_ot_env(1))
616 : END SELECT
617 :
618 829 : IF (delta < eps_gradient .OR. ieigensolver >= max_iter_inner_loop) EXIT eigensolver_loop
619 829 : 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 87 : CALL qs_ot_destroy(qs_ot_env(1))
623 87 : DEALLOCATE (qs_ot_env)
624 87 : CALL dbcsr_release_p(matrix_hc(1)%matrix)
625 87 : CALL dbcsr_release_p(matrix_hc_im(1)%matrix)
626 87 : DEALLOCATE (matrix_hc, matrix_hc_im)
627 :
628 87 : IF (delta < eps_gradient) THEN
629 76 : IF (output_unit > 0 .AND. .NOT. my_silent) THEN
630 : WRITE (UNIT=output_unit, FMT='(T2,A,I0,A)') &
631 0 : 'OT| Complex eigensolver reached convergence in ', iter_total, ' iterations'
632 : END IF
633 : EXIT outer_scf
634 : END IF
635 11 : IF (iter_total >= iter_max) THEN
636 2 : IF (output_unit > 0) THEN
637 : WRITE (output_unit, '(A,T60,E20.10)') &
638 0 : ' 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 78 : CALL copy_dbcsr_to_fm(matrix_c, matrix_c_fm)
645 78 : CALL copy_dbcsr_to_fm(matrix_c_im, matrix_c_fm_im)
646 78 : IF (PRESENT(eigenvalues)) THEN
647 : CALL canonicalize_complex_subspace(matrix_h, matrix_h_im, matrix_c_fm, &
648 78 : matrix_c_fm_im, eigenvalues)
649 : END IF
650 78 : CALL dbcsr_release_p(matrix_c)
651 78 : CALL dbcsr_release_p(matrix_c_im)
652 78 : CALL dbcsr_release_p(matrix_tmp)
653 78 : IF (ASSOCIATED(local_preconditioner)) THEN
654 78 : CALL destroy_preconditioner(local_preconditioner)
655 78 : DEALLOCATE (local_preconditioner)
656 : END IF
657 78 : CALL timestop(handle)
658 :
659 78 : 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 78 : 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 78 : CALL cp_fm_get_info(coeff_re, nrow_global=nao, ncol_global=nmo)
679 78 : CPASSERT(SIZE(eigenvalues) >= nmo)
680 78 : CALL cp_fm_create(hc_re, coeff_re%matrix_struct)
681 78 : CALL cp_fm_create(hc_im, coeff_re%matrix_struct)
682 78 : 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 78 : alpha=-1.0_dp, beta=1.0_dp)
685 78 : 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 78 : alpha=1.0_dp, beta=1.0_dp)
688 :
689 78 : CALL cp_cfm_create(coeff, coeff_re%matrix_struct)
690 78 : CALL cp_cfm_create(hc, coeff_re%matrix_struct)
691 78 : CALL cp_cfm_create(rotated, coeff_re%matrix_struct)
692 78 : CALL cp_cfm_create(hblock, coeff_re%matrix_struct, nrow=nmo, ncol=nmo)
693 78 : CALL cp_cfm_create(eigenvectors, coeff_re%matrix_struct, nrow=nmo, ncol=nmo)
694 78 : CALL cp_fm_to_cfm(coeff_re, coeff_im, coeff)
695 78 : CALL cp_fm_to_cfm(hc_re, hc_im, hc)
696 78 : CALL cp_cfm_gemm('C', 'N', nmo, nmo, nao, z_one, coeff, hc, z_zero, hblock)
697 78 : CALL cp_cfm_heevd(hblock, eigenvectors, eigenvalues(1:nmo))
698 78 : CALL cp_cfm_gemm('N', 'N', nao, nmo, nmo, z_one, coeff, eigenvectors, z_zero, rotated)
699 78 : CALL cp_cfm_to_fm(rotated, coeff_re, coeff_im)
700 :
701 78 : CALL cp_cfm_release(eigenvectors)
702 78 : CALL cp_cfm_release(hblock)
703 78 : CALL cp_cfm_release(rotated)
704 78 : CALL cp_cfm_release(hc)
705 78 : CALL cp_cfm_release(coeff)
706 78 : CALL cp_fm_release(hc_im)
707 78 : CALL cp_fm_release(hc_re)
708 78 : END SUBROUTINE canonicalize_complex_subspace
709 :
710 : END MODULE qs_ot_eigensolver
|
