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