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 basic functionality for using ot in the scf routines.
10 : !> \par History
11 : !> 01.2003 : Joost VandeVondele : adapted for LSD
12 : !> \author Joost VandeVondele (25.08.2002)
13 : ! **************************************************************************************************
14 : MODULE qs_ot_scf
15 : USE cp_array_utils, ONLY: cp_1d_r_p_type
16 : USE cp_dbcsr_api, ONLY: &
17 : dbcsr_copy, dbcsr_get_info, dbcsr_init_p, dbcsr_multiply, dbcsr_p_type, dbcsr_release, &
18 : dbcsr_set, dbcsr_type, dbcsr_type_no_symmetry
19 : USE cp_dbcsr_contrib, ONLY: dbcsr_add_on_diag,&
20 : dbcsr_dot,&
21 : dbcsr_get_diag,&
22 : dbcsr_scale_by_vector,&
23 : dbcsr_set_diag
24 : USE cp_dbcsr_operations, ONLY: copy_fm_to_dbcsr,&
25 : cp_dbcsr_m_by_n_from_row_template
26 : USE cp_fm_types, ONLY: cp_fm_type
27 : USE cp_log_handling, ONLY: cp_get_default_logger,&
28 : cp_logger_type
29 : USE cp_output_handling, ONLY: cp_print_key_finished_output,&
30 : cp_print_key_unit_nr
31 : USE input_section_types, ONLY: section_vals_get,&
32 : section_vals_get_subs_vals,&
33 : section_vals_type
34 : USE kinds, ONLY: dp
35 : USE qs_mo_occupation, ONLY: set_mo_occupation
36 : USE qs_mo_types, ONLY: get_mo_set,&
37 : mo_set_restrict,&
38 : mo_set_type
39 : USE qs_ot, ONLY: qs_ot_get_orbitals,&
40 : qs_ot_get_orbitals_ref,&
41 : qs_ot_get_p
42 : USE qs_ot_minimizer, ONLY: ot_mini
43 : USE qs_ot_types, ONLY: ot_readwrite_input,&
44 : qs_ot_allocate,&
45 : qs_ot_destroy,&
46 : qs_ot_init,&
47 : qs_ot_set_context,&
48 : qs_ot_settings_init,&
49 : qs_ot_type
50 : USE scf_control_types, ONLY: smear_type
51 : #include "./base/base_uses.f90"
52 :
53 : IMPLICIT NONE
54 :
55 : PRIVATE
56 :
57 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_ot_scf'
58 : ! *** Public subroutines ***
59 :
60 : PUBLIC :: ot_scf_init
61 : PUBLIC :: ot_scf_mini
62 : PUBLIC :: ot_scf_destroy
63 : PUBLIC :: ot_scf_read_input
64 :
65 : CONTAINS
66 :
67 : ! **************************************************************************************************
68 : !> \brief ...
69 : !> \param qs_ot_env ...
70 : !> \param scf_section ...
71 : !> \param do_kpoints ...
72 : ! **************************************************************************************************
73 15462 : SUBROUTINE ot_scf_read_input(qs_ot_env, scf_section, do_kpoints)
74 : TYPE(qs_ot_type), DIMENSION(:), POINTER :: qs_ot_env
75 : TYPE(section_vals_type), POINTER :: scf_section
76 : LOGICAL, INTENT(IN) :: do_kpoints
77 :
78 : CHARACTER(len=*), PARAMETER :: routineN = 'ot_scf_read_input'
79 :
80 : INTEGER :: handle, ispin, nspin, output_unit
81 : LOGICAL :: explicit
82 : TYPE(cp_logger_type), POINTER :: logger
83 : TYPE(section_vals_type), POINTER :: ot_section
84 :
85 7731 : CALL timeset(routineN, handle)
86 :
87 7731 : logger => cp_get_default_logger()
88 : output_unit = cp_print_key_unit_nr(logger, scf_section, "PRINT%PROGRAM_RUN_INFO", &
89 7731 : extension=".log")
90 :
91 : ! decide default settings
92 7731 : CALL qs_ot_settings_init(qs_ot_env(1)%settings)
93 :
94 : ! use ot input new style
95 7731 : ot_section => section_vals_get_subs_vals(scf_section, "OT")
96 7731 : CALL section_vals_get(ot_section, explicit=explicit)
97 :
98 : CALL ot_readwrite_input(qs_ot_env(1)%settings, ot_section, output_unit, &
99 7731 : complex_kpoints=do_kpoints)
100 :
101 : CALL cp_print_key_finished_output(output_unit, logger, scf_section, &
102 7731 : "PRINT%PROGRAM_RUN_INFO")
103 :
104 : ! copy the ot settings type so it is identical
105 7731 : nspin = SIZE(qs_ot_env)
106 9440 : DO ispin = 2, nspin
107 9440 : qs_ot_env(ispin)%settings = qs_ot_env(1)%settings
108 : END DO
109 :
110 7731 : CALL timestop(handle)
111 :
112 7731 : END SUBROUTINE ot_scf_read_input
113 :
114 : ! **************************************************************************************************
115 : !> \brief performs the actual minimisation, needs only limited info
116 : !> updated for restricted calculations
117 : !> matrix_dedc is the derivative of the energy with respect to the orbitals (except for a factor 2*fi)
118 : !> a null pointer for matrix_s implies that matrix_s is the unit matrix
119 : !> \param mo_array ...
120 : !> \param matrix_dedc ...
121 : !> \param smear ...
122 : !> \param matrix_s ...
123 : !> \param energy ...
124 : !> \param energy_only ...
125 : !> \param delta ...
126 : !> \param qs_ot_env ...
127 : ! **************************************************************************************************
128 81420 : SUBROUTINE ot_scf_mini(mo_array, matrix_dedc, smear, matrix_s, energy, &
129 : energy_only, delta, qs_ot_env)
130 :
131 : TYPE(mo_set_type), DIMENSION(:), INTENT(INOUT) :: mo_array
132 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_dedc
133 : TYPE(smear_type), POINTER :: smear
134 : TYPE(dbcsr_type), POINTER :: matrix_s
135 : REAL(KIND=dp) :: energy
136 : LOGICAL, INTENT(INOUT) :: energy_only
137 : REAL(KIND=dp) :: delta
138 : TYPE(qs_ot_type), DIMENSION(:), POINTER :: qs_ot_env
139 :
140 : CHARACTER(len=*), PARAMETER :: routineN = 'ot_scf_mini'
141 :
142 : INTEGER :: handle, ispin, k, n, nspin
143 : REAL(KIND=dp) :: ener_nondiag, trace
144 81420 : TYPE(cp_1d_r_p_type), ALLOCATABLE, DIMENSION(:) :: expectation_values, occupation_numbers, &
145 81420 : scaling_factor
146 : TYPE(cp_logger_type), POINTER :: logger
147 81420 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_dedc_physical, matrix_dedc_scaled
148 : TYPE(dbcsr_type), POINTER :: mo_coeff
149 :
150 81420 : CALL timeset(routineN, handle)
151 :
152 81420 : NULLIFY (logger)
153 81420 : logger => cp_get_default_logger()
154 :
155 81420 : nspin = SIZE(mo_array)
156 :
157 341577 : ALLOCATE (occupation_numbers(nspin))
158 260157 : ALLOCATE (scaling_factor(nspin))
159 :
160 81420 : IF (qs_ot_env(1)%settings%do_ener) THEN
161 0 : ALLOCATE (expectation_values(nspin))
162 : END IF
163 :
164 178737 : DO ispin = 1, nspin
165 97317 : CALL get_mo_set(mo_set=mo_array(ispin), occupation_numbers=occupation_numbers(ispin)%array)
166 289839 : ALLOCATE (scaling_factor(ispin)%array(SIZE(occupation_numbers(ispin)%array)))
167 1012270 : scaling_factor(ispin)%array = 2.0_dp*occupation_numbers(ispin)%array
168 178737 : IF (qs_ot_env(1)%settings%do_ener) THEN
169 0 : ALLOCATE (expectation_values(ispin)%array(SIZE(occupation_numbers(ispin)%array)))
170 : END IF
171 : END DO
172 :
173 : ! optimizing orbital energies somehow implies non-equivalent orbitals
174 81420 : IF (qs_ot_env(1)%settings%do_ener) THEN
175 0 : CPASSERT(qs_ot_env(1)%settings%do_rotation)
176 : END IF
177 : ! add_nondiag_energy requires do_ener
178 81420 : IF (qs_ot_env(1)%settings%add_nondiag_energy) THEN
179 0 : CPASSERT(qs_ot_env(1)%settings%do_ener)
180 : END IF
181 :
182 : ! get a rotational force
183 81420 : IF (.NOT. energy_only) THEN
184 65730 : IF (qs_ot_env(1)%settings%do_rotation) THEN
185 3652 : DO ispin = 1, SIZE(qs_ot_env)
186 1882 : CALL get_mo_set(mo_set=mo_array(ispin), mo_coeff_b=mo_coeff)
187 1882 : CALL dbcsr_get_info(mo_coeff, nfullrows_total=n, nfullcols_total=k)
188 : CALL dbcsr_multiply('T', 'N', 1.0_dp, mo_coeff, matrix_dedc(ispin)%matrix, &
189 1882 : 0.0_dp, qs_ot_env(ispin)%rot_mat_chc)
190 1882 : CALL dbcsr_copy(qs_ot_env(ispin)%matrix_buf1, qs_ot_env(ispin)%rot_mat_chc)
191 :
192 1882 : CALL dbcsr_scale_by_vector(qs_ot_env(ispin)%matrix_buf1, alpha=scaling_factor(ispin)%array, side='right')
193 : ! create the derivative of the energy wrt to rot_mat_u
194 : CALL dbcsr_multiply('N', 'N', 1.0_dp, qs_ot_env(ispin)%rot_mat_u, qs_ot_env(ispin)%matrix_buf1, &
195 5534 : 0.0_dp, qs_ot_env(ispin)%rot_mat_dedu)
196 : END DO
197 :
198 : ! here we construct the derivative of the free energy with respect to the evals
199 : ! (note that this requires the diagonal elements of chc)
200 : ! the mo occupations should in principle remain unaltered
201 1770 : IF (qs_ot_env(1)%settings%do_ener) THEN
202 0 : DO ispin = 1, SIZE(mo_array)
203 0 : CALL dbcsr_get_diag(qs_ot_env(ispin)%rot_mat_chc, expectation_values(ispin)%array)
204 0 : qs_ot_env(ispin)%ener_gx = expectation_values(ispin)%array
205 : CALL set_mo_occupation(mo_set=mo_array(ispin), &
206 0 : smear=smear, eval_deriv=qs_ot_env(ispin)%ener_gx)
207 : END DO
208 : END IF
209 :
210 : ! chc only needs to be stored in u independent form if we require add_nondiag_energy,
211 : ! which will use it in non-selfconsistent form for e.g. the linesearch
212 : ! transform C^T H C -> U C^T H C U ^ T
213 1770 : IF (qs_ot_env(1)%settings%add_nondiag_energy) THEN
214 0 : DO ispin = 1, SIZE(qs_ot_env)
215 0 : CALL dbcsr_get_info(qs_ot_env(ispin)%rot_mat_u, nfullcols_total=k)
216 : CALL dbcsr_multiply('N', 'N', 1.0_dp, qs_ot_env(ispin)%rot_mat_u, qs_ot_env(ispin)%rot_mat_chc, &
217 0 : 0.0_dp, qs_ot_env(ispin)%matrix_buf1)
218 : CALL dbcsr_multiply('N', 'T', 1.0_dp, qs_ot_env(ispin)%matrix_buf1, qs_ot_env(ispin)%rot_mat_u, &
219 0 : 0.0_dp, qs_ot_env(ispin)%rot_mat_chc)
220 : END DO
221 : END IF
222 : END IF
223 : END IF
224 :
225 : ! evaluate non-diagonal energy contribution
226 81420 : ener_nondiag = 0.0_dp
227 81420 : IF (qs_ot_env(1)%settings%add_nondiag_energy) THEN
228 0 : DO ispin = 1, SIZE(qs_ot_env)
229 : ! transform \tilde H to the current basis of C (assuming non-selfconsistent H)
230 0 : CALL dbcsr_get_info(qs_ot_env(ispin)%rot_mat_u, nfullcols_total=k)
231 : CALL dbcsr_multiply('T', 'N', 1.0_dp, qs_ot_env(ispin)%rot_mat_u, qs_ot_env(ispin)%rot_mat_chc, &
232 0 : 0.0_dp, qs_ot_env(ispin)%matrix_buf1)
233 : CALL dbcsr_multiply('N', 'N', 1.0_dp, qs_ot_env(ispin)%matrix_buf1, qs_ot_env(ispin)%rot_mat_u, &
234 0 : 0.0_dp, qs_ot_env(ispin)%matrix_buf2)
235 :
236 : ! subtract the current ener_x from the diagonal
237 0 : CALL dbcsr_get_diag(qs_ot_env(ispin)%matrix_buf2, expectation_values(ispin)%array)
238 0 : expectation_values(ispin)%array = expectation_values(ispin)%array - qs_ot_env(ispin)%ener_x
239 0 : CALL dbcsr_set_diag(qs_ot_env(ispin)%matrix_buf2, expectation_values(ispin)%array)
240 :
241 : ! get nondiag energy trace (D^T D)
242 0 : CALL dbcsr_dot(qs_ot_env(ispin)%matrix_buf2, qs_ot_env(ispin)%matrix_buf2, trace)
243 0 : ener_nondiag = ener_nondiag + 0.5_dp*qs_ot_env(1)%settings%nondiag_energy_strength*trace
244 :
245 : ! get gradient (again ignoring dependencies of H)
246 0 : IF (.NOT. energy_only) THEN
247 : ! first for the ener_x (-2*(diag(C^T H C)-ener_x))
248 : qs_ot_env(ispin)%ener_gx = qs_ot_env(ispin)%ener_gx - &
249 0 : qs_ot_env(1)%settings%nondiag_energy_strength*expectation_values(ispin)%array
250 :
251 : ! next for the rot_mat_u derivative (2 * k * \tilde H U D)
252 : CALL dbcsr_multiply('N', 'N', 1.0_dp, qs_ot_env(ispin)%rot_mat_chc, qs_ot_env(ispin)%rot_mat_u, &
253 0 : 0.0_dp, qs_ot_env(ispin)%matrix_buf1)
254 : CALL dbcsr_multiply('N', 'N', 2.0_dp*qs_ot_env(1)%settings%nondiag_energy_strength, &
255 : qs_ot_env(ispin)%matrix_buf1, qs_ot_env(ispin)%matrix_buf2, &
256 0 : 1.0_dp, qs_ot_env(ispin)%rot_mat_dedu)
257 : END IF
258 : END DO
259 : END IF
260 :
261 : ! this is kind of a hack so far (costly memory wise), we locally recreate the scaled matrix_hc, and
262 : ! use it in the following, eventually, as occupations numbers get more integrated, it should become possible
263 : ! to remove this.
264 340675 : ALLOCATE (matrix_dedc_scaled(SIZE(matrix_dedc)))
265 81420 : NULLIFY (matrix_dedc_physical)
266 81420 : IF (qs_ot_env(1)%settings%occupation_preconditioner) THEN
267 0 : ALLOCATE (matrix_dedc_physical(SIZE(matrix_dedc)))
268 : END IF
269 177835 : DO ispin = 1, SIZE(matrix_dedc)
270 96415 : ALLOCATE (matrix_dedc_scaled(ispin)%matrix)
271 96415 : CALL dbcsr_copy(matrix_dedc_scaled(ispin)%matrix, matrix_dedc(ispin)%matrix)
272 :
273 96415 : IF (qs_ot_env(1)%settings%occupation_preconditioner) THEN
274 0 : ALLOCATE (matrix_dedc_physical(ispin)%matrix)
275 0 : CALL dbcsr_copy(matrix_dedc_physical(ispin)%matrix, matrix_dedc(ispin)%matrix)
276 : CALL dbcsr_scale_by_vector(matrix_dedc_physical(ispin)%matrix, &
277 0 : alpha=2.0_dp*occupation_numbers(ispin)%array, side='right')
278 0 : scaling_factor(ispin)%array = 2.0_dp
279 : END IF
280 177835 : CALL dbcsr_scale_by_vector(matrix_dedc_scaled(ispin)%matrix, alpha=scaling_factor(ispin)%array, side='right')
281 : END DO
282 :
283 : ! notice we use qs_ot_env(1) for driving all output and the minimization in case of LSD
284 81420 : qs_ot_env(1)%etotal = energy + ener_nondiag
285 :
286 81420 : IF (qs_ot_env(1)%settings%occupation_preconditioner .AND. &
287 : (qs_ot_env(1)%settings%ot_method == "CG" .OR. &
288 : qs_ot_env(1)%settings%ot_method == "SD")) THEN
289 0 : CALL ot_mini(qs_ot_env, matrix_dedc_scaled, matrix_hc_physical=matrix_dedc_physical)
290 : ELSE
291 81420 : CALL ot_mini(qs_ot_env, matrix_dedc_scaled)
292 : END IF
293 :
294 81420 : delta = qs_ot_env(1)%delta
295 81420 : energy_only = qs_ot_env(1)%energy_only
296 :
297 : ! generate the orbitals using the new matrix_x
298 177835 : DO ispin = 1, SIZE(qs_ot_env)
299 96415 : CALL get_mo_set(mo_set=mo_array(ispin), mo_coeff_b=mo_coeff)
300 96415 : CALL dbcsr_get_info(mo_coeff, nfullrows_total=n, nfullcols_total=k)
301 177835 : SELECT CASE (qs_ot_env(1)%settings%ot_algorithm)
302 : CASE ("TOD")
303 91667 : IF (ASSOCIATED(matrix_s)) THEN
304 : CALL dbcsr_multiply('N', 'N', 1.0_dp, matrix_s, qs_ot_env(ispin)%matrix_x, &
305 72259 : 0.0_dp, qs_ot_env(ispin)%matrix_sx)
306 : ELSE
307 19408 : CALL dbcsr_copy(qs_ot_env(ispin)%matrix_sx, qs_ot_env(ispin)%matrix_x)
308 : END IF
309 91667 : CALL qs_ot_get_p(qs_ot_env(ispin)%matrix_x, qs_ot_env(ispin)%matrix_sx, qs_ot_env(ispin))
310 91667 : CALL qs_ot_get_orbitals(mo_coeff, qs_ot_env(ispin)%matrix_x, qs_ot_env(ispin))
311 : CASE ("REF")
312 : CALL qs_ot_get_orbitals_ref(mo_coeff, matrix_s, qs_ot_env(ispin)%matrix_x, &
313 : qs_ot_env(ispin)%matrix_sx, qs_ot_env(ispin)%matrix_gx_old, &
314 4748 : qs_ot_env(ispin)%matrix_dx, qs_ot_env(ispin), qs_ot_env(1))
315 : CASE DEFAULT
316 96415 : CPABORT("Algorithm not yet implemented")
317 : END SELECT
318 : END DO
319 :
320 81420 : IF (qs_ot_env(1)%restricted) THEN
321 902 : CALL mo_set_restrict(mo_array, convert_dbcsr=.TRUE.)
322 : END IF
323 : !
324 : ! obtain the new set of OT eigenvalues and set the occupations accordingly
325 : !
326 81420 : IF (qs_ot_env(1)%settings%do_ener) THEN
327 0 : DO ispin = 1, SIZE(mo_array)
328 0 : mo_array(ispin)%eigenvalues = qs_ot_env(ispin)%ener_x
329 : CALL set_mo_occupation(mo_set=mo_array(ispin), &
330 0 : smear=smear)
331 : END DO
332 : END IF
333 :
334 : ! cleanup
335 178737 : DO ispin = 1, SIZE(scaling_factor)
336 178737 : DEALLOCATE (scaling_factor(ispin)%array)
337 : END DO
338 81420 : DEALLOCATE (scaling_factor)
339 81420 : IF (qs_ot_env(1)%settings%do_ener) THEN
340 0 : DO ispin = 1, SIZE(expectation_values)
341 0 : DEALLOCATE (expectation_values(ispin)%array)
342 : END DO
343 0 : DEALLOCATE (expectation_values)
344 : END IF
345 81420 : DEALLOCATE (occupation_numbers)
346 177835 : DO ispin = 1, SIZE(matrix_dedc_scaled)
347 96415 : CALL dbcsr_release(matrix_dedc_scaled(ispin)%matrix)
348 177835 : DEALLOCATE (matrix_dedc_scaled(ispin)%matrix)
349 : END DO
350 81420 : DEALLOCATE (matrix_dedc_scaled)
351 81420 : IF (ASSOCIATED(matrix_dedc_physical)) THEN
352 0 : DO ispin = 1, SIZE(matrix_dedc_physical)
353 0 : CALL dbcsr_release(matrix_dedc_physical(ispin)%matrix)
354 0 : DEALLOCATE (matrix_dedc_physical(ispin)%matrix)
355 : END DO
356 0 : DEALLOCATE (matrix_dedc_physical)
357 : END IF
358 :
359 81420 : CALL timestop(handle)
360 :
361 162840 : END SUBROUTINE ot_scf_mini
362 :
363 : ! **************************************************************************************************
364 : !> \brief initialises qs_ot_env so that mo_coeff is the current point
365 : !> and that the minization can be started.
366 : !> \param mo_array ...
367 : !> \param matrix_s ...
368 : !> \param qs_ot_env ...
369 : !> \param matrix_ks ...
370 : !> \param broyden_adaptive_sigma ...
371 : ! **************************************************************************************************
372 7577 : SUBROUTINE ot_scf_init(mo_array, matrix_s, qs_ot_env, matrix_ks, broyden_adaptive_sigma)
373 :
374 : TYPE(mo_set_type), DIMENSION(:), INTENT(IN) :: mo_array
375 : TYPE(dbcsr_type), POINTER :: matrix_s
376 : TYPE(qs_ot_type), DIMENSION(:), POINTER :: qs_ot_env
377 : TYPE(dbcsr_type), POINTER :: matrix_ks
378 : REAL(KIND=dp) :: broyden_adaptive_sigma
379 :
380 : CHARACTER(len=*), PARAMETER :: routineN = 'ot_scf_init'
381 :
382 : INTEGER :: handle, ispin, k, n, nspin
383 : LOGICAL :: is_equal
384 : TYPE(cp_fm_type), POINTER :: mo_coeff_fm
385 : TYPE(dbcsr_type), POINTER :: mo_coeff
386 :
387 7577 : CALL timeset(routineN, handle)
388 :
389 16767 : DO ispin = 1, SIZE(mo_array)
390 16767 : IF (.NOT. ASSOCIATED(mo_array(ispin)%mo_coeff_b)) THEN
391 0 : CPABORT("Shouldn't get there")
392 : ! we do ot then copy fm to dbcsr
393 : ! allocate that somewhere else ! fm -> dbcsr
394 0 : CALL dbcsr_init_p(mo_array(ispin)%mo_coeff_b)
395 : CALL cp_dbcsr_m_by_n_from_row_template(mo_array(ispin)%mo_coeff_b, template=matrix_ks, &
396 : n=mo_array(ispin)%nmo, &
397 0 : sym=dbcsr_type_no_symmetry)
398 : END IF
399 : END DO
400 :
401 : ! *** set a history for broyden
402 16675 : DO ispin = 1, SIZE(qs_ot_env)
403 16675 : qs_ot_env(ispin)%broyden_adaptive_sigma = broyden_adaptive_sigma
404 : END DO
405 :
406 : ! **** SCP
407 : ! **** SCP
408 : ! adapted for work with the restricted keyword
409 7577 : nspin = SIZE(qs_ot_env)
410 :
411 16675 : DO ispin = 1, nspin
412 :
413 9098 : CALL qs_ot_set_context(qs_ot_env(ispin), spin_index=ispin)
414 :
415 9098 : NULLIFY (mo_coeff)
416 9098 : CALL get_mo_set(mo_set=mo_array(ispin), mo_coeff_b=mo_coeff, mo_coeff=mo_coeff_fm)
417 9098 : CALL copy_fm_to_dbcsr(mo_coeff_fm, mo_coeff) !fm -> dbcsr
418 :
419 9098 : CALL dbcsr_get_info(mo_coeff, nfullrows_total=n, nfullcols_total=k)
420 :
421 : ! allocate
422 9098 : CALL qs_ot_allocate(qs_ot_env(ispin), matrix_ks, mo_coeff_fm%matrix_struct)
423 :
424 : ! set c0,sc0
425 9098 : CALL dbcsr_copy(qs_ot_env(ispin)%matrix_c0, mo_coeff)
426 9098 : IF (ASSOCIATED(matrix_s)) THEN
427 : CALL dbcsr_multiply('N', 'N', 1.0_dp, matrix_s, qs_ot_env(ispin)%matrix_c0, &
428 7908 : 0.0_dp, qs_ot_env(ispin)%matrix_sc0)
429 : ELSE
430 1190 : CALL dbcsr_copy(qs_ot_env(ispin)%matrix_sc0, qs_ot_env(ispin)%matrix_c0)
431 : END IF
432 :
433 : ! init
434 9098 : CALL qs_ot_init(qs_ot_env(ispin))
435 :
436 : ! set x
437 9098 : CALL dbcsr_set(qs_ot_env(ispin)%matrix_x, 0.0_dp)
438 9098 : CALL dbcsr_set(qs_ot_env(ispin)%matrix_sx, 0.0_dp)
439 :
440 9098 : IF (qs_ot_env(ispin)%settings%do_rotation) THEN
441 234 : CALL dbcsr_set(qs_ot_env(ispin)%rot_mat_x, 0.0_dp)
442 234 : CALL dbcsr_set(qs_ot_env(ispin)%rot_mat_u, 0.0_dp)
443 234 : CALL dbcsr_add_on_diag(qs_ot_env(ispin)%rot_mat_u, 1.0_dp)
444 : END IF
445 :
446 9098 : IF (qs_ot_env(ispin)%settings%do_ener) THEN
447 0 : is_equal = SIZE(qs_ot_env(ispin)%ener_x) == SIZE(mo_array(ispin)%eigenvalues)
448 0 : CPASSERT(is_equal)
449 0 : qs_ot_env(ispin)%ener_x = mo_array(ispin)%eigenvalues
450 : END IF
451 :
452 16675 : SELECT CASE (qs_ot_env(1)%settings%ot_algorithm)
453 : CASE ("TOD")
454 : ! get c
455 8016 : CALL qs_ot_get_p(qs_ot_env(ispin)%matrix_x, qs_ot_env(ispin)%matrix_sx, qs_ot_env(ispin))
456 : CASE ("REF")
457 1082 : CALL dbcsr_copy(qs_ot_env(ispin)%matrix_x, qs_ot_env(ispin)%matrix_c0)
458 1082 : CALL dbcsr_copy(qs_ot_env(ispin)%matrix_sx, qs_ot_env(ispin)%matrix_sc0)
459 : CASE DEFAULT
460 9098 : CPABORT("Algorithm not yet implemented")
461 : END SELECT
462 :
463 : END DO
464 7577 : CALL timestop(handle)
465 7577 : END SUBROUTINE ot_scf_init
466 :
467 : ! **************************************************************************************************
468 : !> \brief ...
469 : !> \param qs_ot_env ...
470 : ! **************************************************************************************************
471 9440 : SUBROUTINE ot_scf_destroy(qs_ot_env)
472 :
473 : TYPE(qs_ot_type) :: qs_ot_env
474 :
475 9440 : CALL qs_ot_destroy(qs_ot_env)
476 :
477 9440 : END SUBROUTINE ot_scf_destroy
478 :
479 : END MODULE qs_ot_scf
|
