(git:58e3e09)
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-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 ! **************************************************************************************************
20  USE cp_fm_types, ONLY: cp_fm_get_info,&
21  cp_fm_type
23  USE dbcsr_api, ONLY: &
24  dbcsr_copy, dbcsr_dot, dbcsr_init_p, dbcsr_multiply, dbcsr_p_type, dbcsr_release_p, &
25  dbcsr_scale, dbcsr_set, dbcsr_type, dbcsr_type_no_symmetry
26  USE kinds, ONLY: dp
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,&
33  USE qs_ot_minimizer, ONLY: ot_mini
34  USE qs_ot_types, ONLY: qs_ot_allocate,&
36  qs_ot_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  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  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  TYPE(qs_ot_type), DIMENSION(:), POINTER :: qs_ot_env
97 
98  CALL timeset(routinen, handle)
99 
100  output_unit = cp_logger_get_default_io_unit()
101 
102  IF (PRESENT(silent)) THEN
103  my_silent = silent
104  ELSE
105  my_silent = .false.
106  END IF
107 
108  NULLIFY (matrix_c) ! fm->dbcsr
109 
110  CALL cp_fm_get_info(matrix_c_fm, nrow_global=n, ncol_global=k) ! fm->dbcsr
111  ALLOCATE (matrix_c)
112  CALL cp_fm_to_dbcsr_row_template(matrix_c, fm_in=matrix_c_fm, template=matrix_h)
113 
114  iter_total = 0
115 
116  outer_scf: DO
117 
118  NULLIFY (qs_ot_env)
119 
120  NULLIFY (matrix_s_ortho)
121  NULLIFY (matrix_os_ortho)
122  NULLIFY (matrix_buf1_ortho)
123  NULLIFY (matrix_buf2_ortho)
124  NULLIFY (matrix_orthogonal_space)
125 
126  ALLOCATE (qs_ot_env(1))
127  ALLOCATE (matrix_hc(1))
128  NULLIFY (matrix_hc(1)%matrix)
129  CALL dbcsr_init_p(matrix_hc(1)%matrix)
130  CALL dbcsr_copy(matrix_hc(1)%matrix, matrix_c, 'matrix_hc')
131 
132  ortho = .false.
133  IF (PRESENT(matrix_orthogonal_space_fm)) ortho = .true.
134 
135  ! decide settings
136  IF (PRESENT(ot_settings)) THEN
137  qs_ot_env(1)%settings = ot_settings
138  ELSE
139  CALL qs_ot_settings_init(qs_ot_env(1)%settings)
140  ! overwrite defaults
141  qs_ot_env(1)%settings%ds_min = 0.10_dp
142  END IF
143 
144  IF (ortho) THEN
145  ALLOCATE (matrix_orthogonal_space)
146  CALL cp_fm_to_dbcsr_row_template(matrix_orthogonal_space, fm_in=matrix_orthogonal_space_fm, template=matrix_h)
147  CALL cp_fm_get_info(matrix_orthogonal_space_fm, ncol_global=ortho_space_k)
148 
149  IF (PRESENT(size_ortho_space)) ortho_space_k = size_ortho_space
150  ortho_k = ortho_space_k + k
151  ELSE
152  ortho_k = k
153  END IF
154 
155  ! allocate
156  CALL qs_ot_allocate(qs_ot_env(1), matrix_s, matrix_c_fm%matrix_struct, ortho_k=ortho_k)
157 
158  IF (ortho) THEN
159  ! construct an initial guess that is orthogonal to matrix_orthogonal_space
160 
161  CALL dbcsr_init_p(matrix_s_ortho)
162  CALL dbcsr_copy(matrix_s_ortho, matrix_orthogonal_space, name="matrix_s_ortho")
163 
164  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  sym=dbcsr_type_no_symmetry)
167 
168  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  sym=dbcsr_type_no_symmetry)
171 
172  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  sym=dbcsr_type_no_symmetry)
175 
176  CALL dbcsr_multiply('N', 'N', 1.0_dp, matrix_s, matrix_orthogonal_space, &
177  0.0_dp, matrix_s_ortho)
178  CALL dbcsr_multiply('T', 'N', rone, matrix_s_ortho, matrix_s_ortho, &
179  rzero, matrix_os_ortho)
180 
181  CALL cp_dbcsr_cholesky_decompose(matrix_os_ortho, &
182  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  upper_to_full=.true.)
186 
187  CALL dbcsr_multiply('T', 'N', rone, matrix_s_ortho, matrix_c, &
188  rzero, matrix_buf1_ortho)
189  CALL dbcsr_multiply('N', 'N', rone, matrix_os_ortho, matrix_buf1_ortho, &
190  rzero, matrix_buf2_ortho)
191  CALL dbcsr_multiply('N', 'N', -rone, matrix_s_ortho, matrix_buf2_ortho, &
192  rone, matrix_c)
193 
194  ! make matrix_c0 an orthogonal basis, matrix_c contains sc0
195  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  0.0_dp, matrix_c)
198 
199  CALL make_basis_sv(qs_ot_env(1)%matrix_c0, k, matrix_c, &
200  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  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  para_env=qs_ot_env(1)%para_env, blacs_env=qs_ot_env(1)%blacs_env)
209 
210  CALL dbcsr_release_p(matrix_buf1_ortho)
211  CALL dbcsr_release_p(matrix_buf2_ortho)
212  CALL dbcsr_release_p(matrix_os_ortho)
213  CALL dbcsr_release_p(matrix_s_ortho)
214 
215  ELSE
216 
217  ! set c0,sc0
218  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  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  qs_ot_env(1)%para_env, qs_ot_env(1)%blacs_env)
224  END IF
225 
226  ! init
227  CALL qs_ot_init(qs_ot_env(1))
228  energy_only = qs_ot_env(1)%energy_only
229 
230  ! set x
231  CALL dbcsr_set(qs_ot_env(1)%matrix_x, 0.0_dp)
232  CALL dbcsr_set(qs_ot_env(1)%matrix_sx, 0.0_dp)
233 
234  ! get c
235  CALL qs_ot_get_p(qs_ot_env(1)%matrix_x, qs_ot_env(1)%matrix_sx, qs_ot_env(1))
236  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  IF (PRESENT(preconditioner)) THEN
241  IF (ASSOCIATED(preconditioner)) THEN
243  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  ieigensolver = ieigensolver + 1
255  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  0.0_dp, matrix_hc(1)%matrix)
260  CALL dbcsr_dot(matrix_c, matrix_hc(1)%matrix, energy)
261  IF (.NOT. energy_only) THEN
262  CALL dbcsr_scale(matrix_hc(1)%matrix, 2.0_dp)
263  END IF
264 
265  qs_ot_env(1)%etotal = energy
266  CALL ot_mini(qs_ot_env, matrix_hc)
267  delta = qs_ot_env(1)%delta
268  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  0.0_dp, qs_ot_env(1)%matrix_sx)
272 
273  CALL qs_ot_get_p(qs_ot_env(1)%matrix_x, qs_ot_env(1)%matrix_sx, qs_ot_env(1))
274  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  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  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  CALL qs_ot_destroy(qs_ot_env(1))
284  DEALLOCATE (qs_ot_env)
285  CALL dbcsr_release_p(matrix_hc(1)%matrix)
286  DEALLOCATE (matrix_hc)
287  CALL dbcsr_release_p(matrix_orthogonal_space)
288 
289  IF (delta < eps_gradient) THEN
290  IF ((output_unit > 0) .AND. .NOT. my_silent) THEN
291  WRITE (unit=output_unit, fmt="(T2,A,I0,A)") &
292  "OT| Eigensolver reached convergence in ", iter_total, " iterations"
293  END IF
294  EXIT outer_scf
295  END IF
296  IF (iter_total >= iter_max) THEN
297  IF (output_unit > 0) THEN
298  IF (my_silent) THEN
299  WRITE (output_unit, "(A,T60,E20.10)") " WARNING OT eigensolver did not converge: current gradient", delta
300  ELSE
301  WRITE (output_unit, *) "WARNING : did not converge in ot_eigensolver"
302  WRITE (output_unit, *) "number of iterations ", iter_total, " exceeded maximum"
303  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  CALL copy_dbcsr_to_fm(matrix_c, matrix_c_fm) ! fm->dbcsr
312  CALL dbcsr_release_p(matrix_c) ! fm->dbcsr
313 
314  CALL timestop(handle)
315 
316  END SUBROUTINE ot_eigensolver
317 
318 END MODULE qs_ot_eigensolver
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, upper_to_full)
used to replace the cholesky decomposition by the inverse
DBCSR operations in CP2K.
subroutine, public cp_dbcsr_m_by_n_from_template(matrix, template, m, n, sym, data_type)
Utility function to create an arbitrary shaped dbcsr matrix with the same processor grid as the templ...
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 dbcsr_copy_columns_hack(matrix_b, matrix_a, ncol, source_start, target_start, para_env, blacs_env)
hack for dbcsr_copy_columns
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
Definition: cp_fm_types.F:1016
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...
Defines the basic variable types.
Definition: kinds.F:23
integer, parameter, public dp
Definition: kinds.F:34
types of preconditioners
logical function, public preconditioner_in_use(preconditioner)
...
computes preconditioners, and implements methods to apply them currently used in qs_ot
collects routines that perform operations directly related to MOs
Definition: qs_mo_methods.F:14
an eigen-space solver for the generalised symmetric eigenvalue problem for sparse matrices,...
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)
...
orbital transformations
Definition: qs_ot_types.F:15
subroutine, public qs_ot_init(qs_ot_env)
...
Definition: qs_ot_types.F:254
subroutine, public qs_ot_allocate(qs_ot_env, matrix_s, fm_struct_ref, ortho_k)
...
Definition: qs_ot_types.F:303
subroutine, public qs_ot_settings_init(settings)
sets default values for the settings type
Definition: qs_ot_types.F:215
subroutine, public qs_ot_destroy(qs_ot_env)
...
Definition: qs_ot_types.F:636
orbital transformations
Definition: qs_ot.F:15
subroutine, public qs_ot_get_p(matrix_x, matrix_sx, qs_ot_env)
...
Definition: qs_ot.F:752
subroutine, public qs_ot_get_orbitals(matrix_c, matrix_x, qs_ot_env)
...
Definition: qs_ot.F:1018
subroutine, public qs_ot_new_preconditioner(qs_ot_env, preconditioner)
...
Definition: qs_ot.F:70