Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 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_dbcsr_api, ONLY: &
15 : dbcsr_copy, dbcsr_init_p, dbcsr_multiply, dbcsr_p_type, dbcsr_release_p, dbcsr_scale, &
16 : dbcsr_set, dbcsr_type, dbcsr_type_no_symmetry
17 : USE cp_dbcsr_cholesky, ONLY: cp_dbcsr_cholesky_decompose,&
18 : cp_dbcsr_cholesky_invert
19 : USE cp_dbcsr_contrib, ONLY: dbcsr_dot
20 : USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
21 : cp_dbcsr_m_by_n_from_template,&
22 : cp_fm_to_dbcsr_row_template,&
23 : dbcsr_copy_columns_hack
24 : USE cp_fm_types, ONLY: cp_fm_get_info,&
25 : cp_fm_type
26 : USE cp_log_handling, ONLY: cp_logger_get_default_io_unit
27 : USE kinds, ONLY: dp
28 : USE preconditioner_types, ONLY: preconditioner_in_use,&
29 : preconditioner_type
30 : USE qs_mo_methods, ONLY: make_basis_sv
31 : USE qs_ot, ONLY: qs_ot_get_orbitals,&
32 : qs_ot_get_p,&
33 : qs_ot_new_preconditioner
34 : USE qs_ot_minimizer, ONLY: ot_mini
35 : USE qs_ot_types, ONLY: qs_ot_allocate,&
36 : qs_ot_destroy,&
37 : qs_ot_init,&
38 : qs_ot_settings_init,&
39 : qs_ot_settings_type,&
40 : qs_ot_type
41 : #include "./base/base_uses.f90"
42 :
43 : IMPLICIT NONE
44 : PRIVATE
45 :
46 : ! *** Global parameters ***
47 :
48 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_ot_eigensolver'
49 :
50 : ! *** Public subroutines ***
51 :
52 : PUBLIC :: ot_eigensolver
53 :
54 : CONTAINS
55 :
56 : ! on input c contains the initial guess (should not be zero !)
57 : ! on output c spans the subspace
58 : ! **************************************************************************************************
59 : !> \brief ...
60 : !> \param matrix_h ...
61 : !> \param matrix_s ...
62 : !> \param matrix_orthogonal_space_fm ...
63 : !> \param matrix_c_fm ...
64 : !> \param preconditioner ...
65 : !> \param eps_gradient ...
66 : !> \param iter_max ...
67 : !> \param size_ortho_space ...
68 : !> \param silent ...
69 : !> \param ot_settings ...
70 : ! **************************************************************************************************
71 704 : SUBROUTINE ot_eigensolver(matrix_h, matrix_s, matrix_orthogonal_space_fm, &
72 : matrix_c_fm, preconditioner, eps_gradient, &
73 : iter_max, size_ortho_space, silent, ot_settings)
74 :
75 : TYPE(dbcsr_type), POINTER :: matrix_h, matrix_s
76 : TYPE(cp_fm_type), INTENT(IN), OPTIONAL :: matrix_orthogonal_space_fm
77 : TYPE(cp_fm_type), INTENT(INOUT) :: matrix_c_fm
78 : TYPE(preconditioner_type), OPTIONAL, POINTER :: preconditioner
79 : REAL(KIND=dp) :: eps_gradient
80 : INTEGER, INTENT(IN) :: iter_max
81 : INTEGER, INTENT(IN), OPTIONAL :: size_ortho_space
82 : LOGICAL, INTENT(IN), OPTIONAL :: silent
83 : TYPE(qs_ot_settings_type), INTENT(IN), OPTIONAL :: ot_settings
84 :
85 : CHARACTER(len=*), PARAMETER :: routineN = 'ot_eigensolver'
86 : INTEGER, PARAMETER :: max_iter_inner_loop = 40
87 : REAL(KIND=dp), PARAMETER :: rone = 1.0_dp, rzero = 0.0_dp
88 :
89 : INTEGER :: handle, ieigensolver, iter_total, k, n, &
90 : ortho_k, ortho_space_k, output_unit
91 : LOGICAL :: energy_only, my_silent, ortho
92 : REAL(KIND=dp) :: delta, energy
93 352 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_hc
94 : TYPE(dbcsr_type), POINTER :: matrix_buf1_ortho, matrix_buf2_ortho, &
95 : matrix_c, matrix_orthogonal_space, &
96 : matrix_os_ortho, matrix_s_ortho
97 352 : TYPE(qs_ot_type), DIMENSION(:), POINTER :: qs_ot_env
98 :
99 352 : CALL timeset(routineN, handle)
100 :
101 352 : output_unit = cp_logger_get_default_io_unit()
102 :
103 352 : IF (PRESENT(silent)) THEN
104 116 : my_silent = silent
105 : ELSE
106 : my_silent = .FALSE.
107 : END IF
108 :
109 352 : NULLIFY (matrix_c) ! fm->dbcsr
110 :
111 352 : CALL cp_fm_get_info(matrix_c_fm, nrow_global=n, ncol_global=k) ! fm->dbcsr
112 352 : ALLOCATE (matrix_c)
113 352 : CALL cp_fm_to_dbcsr_row_template(matrix_c, fm_in=matrix_c_fm, template=matrix_h)
114 :
115 352 : iter_total = 0
116 :
117 : outer_scf: DO
118 :
119 : NULLIFY (qs_ot_env)
120 :
121 534 : NULLIFY (matrix_s_ortho)
122 534 : NULLIFY (matrix_os_ortho)
123 534 : NULLIFY (matrix_buf1_ortho)
124 534 : NULLIFY (matrix_buf2_ortho)
125 534 : NULLIFY (matrix_orthogonal_space)
126 :
127 85974 : ALLOCATE (qs_ot_env(1))
128 1068 : ALLOCATE (matrix_hc(1))
129 534 : NULLIFY (matrix_hc(1)%matrix)
130 534 : CALL dbcsr_init_p(matrix_hc(1)%matrix)
131 534 : CALL dbcsr_copy(matrix_hc(1)%matrix, matrix_c, 'matrix_hc')
132 :
133 534 : ortho = .FALSE.
134 534 : IF (PRESENT(matrix_orthogonal_space_fm)) ortho = .TRUE.
135 :
136 : ! decide settings
137 534 : IF (PRESENT(ot_settings)) THEN
138 122 : qs_ot_env(1)%settings = ot_settings
139 : ELSE
140 412 : CALL qs_ot_settings_init(qs_ot_env(1)%settings)
141 : ! overwrite defaults
142 412 : qs_ot_env(1)%settings%ds_min = 0.10_dp
143 : END IF
144 :
145 534 : IF (ortho) THEN
146 412 : ALLOCATE (matrix_orthogonal_space)
147 412 : CALL cp_fm_to_dbcsr_row_template(matrix_orthogonal_space, fm_in=matrix_orthogonal_space_fm, template=matrix_h)
148 412 : CALL cp_fm_get_info(matrix_orthogonal_space_fm, ncol_global=ortho_space_k)
149 :
150 412 : IF (PRESENT(size_ortho_space)) ortho_space_k = size_ortho_space
151 412 : ortho_k = ortho_space_k + k
152 : ELSE
153 122 : ortho_k = k
154 : END IF
155 :
156 : ! allocate
157 534 : CALL qs_ot_allocate(qs_ot_env(1), matrix_s, matrix_c_fm%matrix_struct, ortho_k=ortho_k)
158 :
159 534 : IF (ortho) THEN
160 : ! construct an initial guess that is orthogonal to matrix_orthogonal_space
161 :
162 412 : CALL dbcsr_init_p(matrix_s_ortho)
163 412 : CALL dbcsr_copy(matrix_s_ortho, matrix_orthogonal_space, name="matrix_s_ortho")
164 :
165 412 : CALL dbcsr_init_p(matrix_os_ortho)
166 : CALL cp_dbcsr_m_by_n_from_template(matrix_os_ortho, template=matrix_h, m=ortho_space_k, n=ortho_space_k, &
167 412 : sym=dbcsr_type_no_symmetry)
168 :
169 412 : CALL dbcsr_init_p(matrix_buf1_ortho)
170 : CALL cp_dbcsr_m_by_n_from_template(matrix_buf1_ortho, template=matrix_h, m=ortho_space_k, n=k, &
171 412 : sym=dbcsr_type_no_symmetry)
172 :
173 412 : CALL dbcsr_init_p(matrix_buf2_ortho)
174 : CALL cp_dbcsr_m_by_n_from_template(matrix_buf2_ortho, template=matrix_h, m=ortho_space_k, n=k, &
175 412 : sym=dbcsr_type_no_symmetry)
176 :
177 : CALL dbcsr_multiply('N', 'N', 1.0_dp, matrix_s, matrix_orthogonal_space, &
178 412 : 0.0_dp, matrix_s_ortho)
179 : CALL dbcsr_multiply('T', 'N', rone, matrix_s_ortho, matrix_s_ortho, &
180 412 : rzero, matrix_os_ortho)
181 :
182 : CALL cp_dbcsr_cholesky_decompose(matrix_os_ortho, &
183 412 : para_env=qs_ot_env(1)%para_env, blacs_env=qs_ot_env(1)%blacs_env)
184 : CALL cp_dbcsr_cholesky_invert(matrix_os_ortho, &
185 : para_env=qs_ot_env(1)%para_env, blacs_env=qs_ot_env(1)%blacs_env, &
186 412 : uplo_to_full=.TRUE.)
187 :
188 : CALL dbcsr_multiply('T', 'N', rone, matrix_s_ortho, matrix_c, &
189 412 : rzero, matrix_buf1_ortho)
190 : CALL dbcsr_multiply('N', 'N', rone, matrix_os_ortho, matrix_buf1_ortho, &
191 412 : rzero, matrix_buf2_ortho)
192 : CALL dbcsr_multiply('N', 'N', -rone, matrix_s_ortho, matrix_buf2_ortho, &
193 412 : rone, matrix_c)
194 :
195 : ! make matrix_c0 an orthogonal basis, matrix_c contains sc0
196 412 : CALL dbcsr_copy(qs_ot_env(1)%matrix_c0, matrix_c)
197 : CALL dbcsr_multiply('N', 'N', 1.0_dp, matrix_s, qs_ot_env(1)%matrix_c0, &
198 412 : 0.0_dp, matrix_c)
199 :
200 : CALL make_basis_sv(qs_ot_env(1)%matrix_c0, k, matrix_c, &
201 412 : qs_ot_env(1)%para_env, qs_ot_env(1)%blacs_env)
202 :
203 : ! copy sc0 and matrix_s_ortho in qs_ot_env(1)%matrix_sc0
204 : !CALL dbcsr_copy_columns(qs_ot_env(1)%matrix_sc0,matrix_s_ortho,ortho_space_k,1,1)
205 : CALL dbcsr_copy_columns_hack(qs_ot_env(1)%matrix_sc0, matrix_s_ortho, ortho_space_k, 1, 1, &
206 412 : para_env=qs_ot_env(1)%para_env, blacs_env=qs_ot_env(1)%blacs_env)
207 : !CALL dbcsr_copy_columns(qs_ot_env(1)%matrix_sc0,matrix_c,k,1,ortho_space_k+1)
208 : CALL dbcsr_copy_columns_hack(qs_ot_env(1)%matrix_sc0, matrix_c, k, 1, ortho_space_k + 1, &
209 412 : para_env=qs_ot_env(1)%para_env, blacs_env=qs_ot_env(1)%blacs_env)
210 :
211 412 : CALL dbcsr_release_p(matrix_buf1_ortho)
212 412 : CALL dbcsr_release_p(matrix_buf2_ortho)
213 412 : CALL dbcsr_release_p(matrix_os_ortho)
214 412 : CALL dbcsr_release_p(matrix_s_ortho)
215 :
216 : ELSE
217 :
218 : ! set c0,sc0
219 122 : CALL dbcsr_copy(qs_ot_env(1)%matrix_c0, matrix_c)
220 : CALL dbcsr_multiply('N', 'N', 1.0_dp, matrix_s, qs_ot_env(1)%matrix_c0, &
221 122 : 0.0_dp, qs_ot_env(1)%matrix_sc0)
222 :
223 : CALL make_basis_sv(qs_ot_env(1)%matrix_c0, k, qs_ot_env(1)%matrix_sc0, &
224 122 : qs_ot_env(1)%para_env, qs_ot_env(1)%blacs_env)
225 : END IF
226 :
227 : ! init
228 534 : CALL qs_ot_init(qs_ot_env(1))
229 534 : energy_only = qs_ot_env(1)%energy_only
230 :
231 : ! set x
232 534 : CALL dbcsr_set(qs_ot_env(1)%matrix_x, 0.0_dp)
233 534 : CALL dbcsr_set(qs_ot_env(1)%matrix_sx, 0.0_dp)
234 :
235 : ! get c
236 534 : CALL qs_ot_get_p(qs_ot_env(1)%matrix_x, qs_ot_env(1)%matrix_sx, qs_ot_env(1))
237 534 : CALL qs_ot_get_orbitals(matrix_c, qs_ot_env(1)%matrix_x, qs_ot_env(1))
238 :
239 : ! if present preconditioner, use it
240 :
241 534 : IF (PRESENT(preconditioner)) THEN
242 534 : IF (ASSOCIATED(preconditioner)) THEN
243 296 : IF (preconditioner_in_use(preconditioner)) THEN
244 296 : CALL qs_ot_new_preconditioner(qs_ot_env(1), preconditioner)
245 : ELSE
246 : ! we should presumably make one
247 : END IF
248 : END IF
249 : END IF
250 :
251 : ! *** Eigensolver loop ***
252 : ieigensolver = 0
253 11946 : eigensolver_loop: DO
254 :
255 11946 : ieigensolver = ieigensolver + 1
256 11946 : iter_total = iter_total + 1
257 :
258 : ! the energy is cHc, the gradient is 2*H*c
259 : CALL dbcsr_multiply('N', 'N', 1.0_dp, matrix_h, matrix_c, &
260 11946 : 0.0_dp, matrix_hc(1)%matrix)
261 11946 : CALL dbcsr_dot(matrix_c, matrix_hc(1)%matrix, energy)
262 11946 : IF (.NOT. energy_only) THEN
263 6240 : CALL dbcsr_scale(matrix_hc(1)%matrix, 2.0_dp)
264 : END IF
265 :
266 11946 : qs_ot_env(1)%etotal = energy
267 11946 : CALL ot_mini(qs_ot_env, matrix_hc)
268 11946 : delta = qs_ot_env(1)%delta
269 11946 : energy_only = qs_ot_env(1)%energy_only
270 :
271 : CALL dbcsr_multiply('N', 'N', 1.0_dp, matrix_s, qs_ot_env(1)%matrix_x, &
272 11946 : 0.0_dp, qs_ot_env(1)%matrix_sx)
273 :
274 11946 : CALL qs_ot_get_p(qs_ot_env(1)%matrix_x, qs_ot_env(1)%matrix_sx, qs_ot_env(1))
275 11946 : CALL qs_ot_get_orbitals(matrix_c, qs_ot_env(1)%matrix_x, qs_ot_env(1))
276 :
277 : ! exit on convergence or if maximum of inner loop cycles is reached
278 11946 : IF (delta < eps_gradient .OR. ieigensolver >= max_iter_inner_loop) EXIT eigensolver_loop
279 : ! exit if total number of steps is reached, but not during a line search step
280 11946 : IF (iter_total >= iter_max .AND. qs_ot_env(1)%OT_METHOD_FULL /= "OT LS") EXIT eigensolver_loop
281 :
282 : END DO eigensolver_loop
283 :
284 534 : CALL qs_ot_destroy(qs_ot_env(1))
285 534 : DEALLOCATE (qs_ot_env)
286 534 : CALL dbcsr_release_p(matrix_hc(1)%matrix)
287 534 : DEALLOCATE (matrix_hc)
288 534 : CALL dbcsr_release_p(matrix_orthogonal_space)
289 :
290 534 : IF (delta < eps_gradient) THEN
291 302 : IF ((output_unit > 0) .AND. .NOT. my_silent) THEN
292 : WRITE (UNIT=output_unit, FMT="(T2,A,I0,A)") &
293 122 : "OT| Eigensolver reached convergence in ", iter_total, " iterations"
294 : END IF
295 : EXIT outer_scf
296 : END IF
297 232 : IF (iter_total >= iter_max) THEN
298 50 : IF (output_unit > 0) THEN
299 25 : IF (my_silent) THEN
300 25 : WRITE (output_unit, "(A,T60,E20.10)") " WARNING OT eigensolver did not converge: current gradient", delta
301 : ELSE
302 0 : WRITE (output_unit, *) "WARNING : did not converge in ot_eigensolver"
303 0 : WRITE (output_unit, *) "number of iterations ", iter_total, " exceeded maximum"
304 0 : WRITE (output_unit, *) "current gradient / target gradient", delta, " / ", eps_gradient
305 : END IF
306 : END IF
307 : EXIT outer_scf
308 : END IF
309 :
310 : END DO outer_scf
311 :
312 352 : CALL copy_dbcsr_to_fm(matrix_c, matrix_c_fm) ! fm->dbcsr
313 352 : CALL dbcsr_release_p(matrix_c) ! fm->dbcsr
314 :
315 352 : CALL timestop(handle)
316 :
317 352 : END SUBROUTINE ot_eigensolver
318 :
319 : END MODULE qs_ot_eigensolver
|
