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