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