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 Utility routines for qs_scf
10 : ! **************************************************************************************************
11 : MODULE qs_scf_initialization
12 : USE atomic_kind_types, ONLY: atomic_kind_type
13 : USE cp_control_types, ONLY: dft_control_type
14 : USE cp_dbcsr_api, ONLY: dbcsr_create,&
15 : dbcsr_init_p,&
16 : dbcsr_p_type,&
17 : dbcsr_type,&
18 : dbcsr_type_no_symmetry
19 : USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
20 : copy_fm_to_dbcsr,&
21 : cp_dbcsr_m_by_n_from_row_template,&
22 : cp_dbcsr_sm_fm_multiply
23 : USE cp_dbcsr_output, ONLY: write_fm_with_basis_info
24 : USE cp_fm_basic_linalg, ONLY: cp_fm_column_scale,&
25 : cp_fm_row_scale,&
26 : cp_fm_transpose,&
27 : cp_fm_triangular_invert
28 : USE cp_fm_cholesky, ONLY: cp_fm_cholesky_decompose
29 : USE cp_fm_diag, ONLY: FM_DIAG_TYPE_CUSOLVER,&
30 : choose_eigv_solver,&
31 : cp_fm_power,&
32 : cusolver_generalized,&
33 : diag_type
34 : USE cp_fm_pool_types, ONLY: cp_fm_pool_p_type,&
35 : fm_pool_get_el_struct
36 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
37 : cp_fm_struct_get,&
38 : cp_fm_struct_release,&
39 : cp_fm_struct_type
40 : USE cp_fm_types, ONLY: cp_fm_create,&
41 : cp_fm_get_info,&
42 : cp_fm_release,&
43 : cp_fm_set_all,&
44 : cp_fm_to_fm,&
45 : cp_fm_to_fm_triangular,&
46 : cp_fm_type
47 : USE cp_log_handling, ONLY: cp_get_default_logger,&
48 : cp_logger_type,&
49 : cp_to_string
50 : USE cp_output_handling, ONLY: cp_p_file,&
51 : cp_print_key_finished_output,&
52 : cp_print_key_should_output,&
53 : cp_print_key_unit_nr
54 : USE hairy_probes, ONLY: AO_boundaries
55 : USE input_constants, ONLY: &
56 : broy_mix, cholesky_dbcsr, cholesky_inverse, cholesky_off, diag_block_davidson, &
57 : diag_block_krylov, diag_filter_matrix, diag_ot, diag_standard, direct_p_mix, kerker_mix, &
58 : multisec_mix, no_mix, ot2cdft, outer_scf_none, plus_u_lowdin, pulay_mix, &
59 : smeagol_runtype_emtransport, wfi_frozen_method_nr, wfi_ps_method_nr, &
60 : wfi_use_guess_method_nr
61 : USE input_section_types, ONLY: section_vals_get_subs_vals,&
62 : section_vals_type,&
63 : section_vals_val_get
64 : USE kinds, ONLY: dp
65 : USE kpoint_types, ONLY: kpoint_type
66 : USE message_passing, ONLY: mp_para_env_type
67 : USE parallel_gemm_api, ONLY: parallel_gemm
68 : USE particle_types, ONLY: particle_type
69 : USE pw_types, ONLY: pw_c1d_gs_type
70 : USE qmmm_image_charge, ONLY: conditional_calc_image_matrix
71 : USE qs_block_davidson_types, ONLY: block_davidson_allocate,&
72 : block_davidson_env_create
73 : USE qs_cdft_opt_types, ONLY: cdft_opt_type_copy
74 : USE qs_density_mixing_types, ONLY: direct_mixing_nr,&
75 : mixing_storage_create,&
76 : mixing_storage_release,&
77 : no_mixing_nr
78 : USE qs_environment_types, ONLY: get_qs_env,&
79 : qs_environment_type,&
80 : set_qs_env
81 : USE qs_fb_distribution_methods, ONLY: fb_distribution_build
82 : USE qs_fb_env_methods, ONLY: fb_env_build_atomic_halos,&
83 : fb_env_build_rcut_auto,&
84 : fb_env_read_input,&
85 : fb_env_write_info
86 : USE qs_fb_env_types, ONLY: fb_env_create,&
87 : fb_env_has_data
88 : USE qs_harris_types, ONLY: harris_type
89 : USE qs_harris_utils, ONLY: harris_density_update
90 : USE qs_initial_guess, ONLY: calculate_first_density_matrix
91 : USE qs_kind_types, ONLY: get_qs_kind,&
92 : qs_kind_type,&
93 : set_qs_kind
94 : USE qs_ks_types, ONLY: qs_ks_did_change
95 : USE qs_matrix_pools, ONLY: mpools_get
96 : USE qs_mixing_utils, ONLY: charge_mixing_init,&
97 : mixing_allocate,&
98 : mixing_init
99 : USE qs_mo_occupation, ONLY: set_mo_occupation
100 : USE qs_mo_types, ONLY: get_mo_set,&
101 : init_mo_set,&
102 : mo_set_type,&
103 : set_mo_set
104 : USE qs_outer_scf, ONLY: outer_loop_extrapolate,&
105 : outer_loop_switch,&
106 : outer_loop_variables_count
107 : USE qs_rho_atom_types, ONLY: rho_atom_type
108 : USE qs_rho_methods, ONLY: duplicate_rho_type,&
109 : qs_rho_update_rho
110 : USE qs_rho_types, ONLY: qs_rho_create,&
111 : qs_rho_get,&
112 : qs_rho_type
113 : USE qs_scf_diagonalization, ONLY: diag_kp_smat,&
114 : diag_subspace_allocate
115 : USE qs_scf_lanczos, ONLY: krylov_space_allocate
116 : USE qs_scf_output, ONLY: qs_scf_initial_info
117 : USE qs_scf_types, ONLY: &
118 : block_davidson_diag_method_nr, block_krylov_diag_method_nr, diag_subspace_env_create, &
119 : filter_matrix_diag_method_nr, general_diag_method_nr, krylov_space_create, &
120 : ot_diag_method_nr, ot_method_nr, qs_scf_env_type, scf_env_create, smeagol_method_nr, &
121 : special_diag_method_nr
122 : USE qs_wf_history_methods, ONLY: reorthogonalize_vectors,&
123 : wfi_extrapolate,&
124 : wfi_get_method_label,&
125 : wfi_update
126 : USE scf_control_types, ONLY: scf_control_type
127 : USE xas_env_types, ONLY: xas_environment_type
128 : USE xas_restart, ONLY: xas_initialize_rho
129 : #include "./base/base_uses.f90"
130 :
131 : IMPLICIT NONE
132 :
133 : PRIVATE
134 :
135 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_scf_initialization'
136 :
137 : PUBLIC:: qs_scf_env_initialize, qs_scf_env_init_basic
138 :
139 : CONTAINS
140 :
141 : ! **************************************************************************************************
142 : !> \brief initializes input parameters if needed or restores values from
143 : !> previous runs to fill scf_env with the values required for scf
144 : !> \param qs_env the qs_environment where to perform the scf procedure
145 : !> \param scf_env ...
146 : !> \param scf_control ...
147 : !> \param scf_section ...
148 : ! **************************************************************************************************
149 21249 : SUBROUTINE qs_scf_env_initialize(qs_env, scf_env, scf_control, scf_section)
150 : TYPE(qs_environment_type), POINTER :: qs_env
151 : TYPE(qs_scf_env_type), POINTER :: scf_env
152 : TYPE(scf_control_type), OPTIONAL, POINTER :: scf_control
153 : TYPE(section_vals_type), OPTIONAL, POINTER :: scf_section
154 :
155 : INTEGER :: ip, np
156 21249 : TYPE(atomic_kind_type), POINTER :: atomic_kind_set(:)
157 : TYPE(dft_control_type), POINTER :: dft_control
158 21249 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
159 21249 : TYPE(particle_type), POINTER :: particle_set(:)
160 21249 : TYPE(qs_kind_type), POINTER :: qs_kind_set(:)
161 : TYPE(scf_control_type), POINTER :: my_scf_control
162 : TYPE(section_vals_type), POINTER :: dft_section, input, my_scf_section
163 :
164 21249 : CALL get_qs_env(qs_env, input=input, dft_control=dft_control)
165 :
166 : !Initialize Hairy Probe calculation
167 21249 : IF (dft_control%hairy_probes .EQV. .TRUE.) THEN
168 : CALL get_qs_env(qs_env, &
169 : mos=mos, &
170 : atomic_kind_set=atomic_kind_set, &
171 : qs_kind_set=qs_kind_set, &
172 4 : particle_set=particle_set)
173 4 : np = SIZE(dft_control%probe)
174 12 : DO ip = 1, np
175 : CALL AO_boundaries(probe=dft_control%probe(ip), atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set, &
176 12 : particle_set=particle_set, nAO=mos(1)%nao) !FIX THIS!
177 : END DO
178 : END IF
179 :
180 21249 : IF (PRESENT(scf_control)) THEN
181 82 : my_scf_control => scf_control
182 : ELSE
183 21167 : CALL get_qs_env(qs_env, scf_control=my_scf_control)
184 : END IF
185 :
186 21249 : dft_section => section_vals_get_subs_vals(input, "DFT")
187 21249 : IF (PRESENT(scf_section)) THEN
188 82 : my_scf_section => scf_section
189 : ELSE
190 21167 : my_scf_section => section_vals_get_subs_vals(dft_section, "SCF")
191 : END IF
192 :
193 21249 : CALL qs_scf_ensure_scf_env(qs_env, scf_env)
194 :
195 21249 : CALL section_vals_val_get(my_scf_section, "CHOLESKY", i_val=scf_env%cholesky_method)
196 :
197 21249 : CALL qs_scf_ensure_mos(qs_env)
198 :
199 : ! set flags for diagonalization
200 : CALL qs_scf_ensure_diagonalization(scf_env, my_scf_section, qs_env, &
201 21249 : my_scf_control, qs_env%has_unit_metric)
202 : ! set parameters for mixing/DIIS during scf
203 21249 : CALL qs_scf_ensure_mixing(my_scf_control, my_scf_section, scf_env, dft_control)
204 :
205 21249 : CALL qs_scf_ensure_work_matrices(qs_env, scf_env)
206 :
207 21249 : CALL qs_scf_ensure_mixing_store(qs_env, scf_env)
208 :
209 : ! Initialize outer loop variables: handle CDFT and regular outer loop separately
210 21249 : IF (dft_control%qs_control%cdft) THEN
211 : CALL qs_scf_ensure_cdft_loop_vars(qs_env, scf_env, dft_control, &
212 326 : scf_control=my_scf_control)
213 : ELSE
214 20923 : CALL qs_scf_ensure_outer_loop_vars(scf_env, my_scf_control)
215 : END IF
216 :
217 21249 : CALL init_scf_run(scf_env, qs_env, my_scf_section, my_scf_control)
218 :
219 21249 : END SUBROUTINE qs_scf_env_initialize
220 :
221 : ! **************************************************************************************************
222 : !> \brief initializes input parameters if needed for non-scf calclulations using diagonalization
223 : !> \param qs_env the qs_environment where to perform the scf procedure
224 : !> \param scf_env ...
225 : ! **************************************************************************************************
226 2 : SUBROUTINE qs_scf_env_init_basic(qs_env, scf_env)
227 : TYPE(qs_environment_type), POINTER :: qs_env
228 : TYPE(qs_scf_env_type), POINTER :: scf_env
229 :
230 : TYPE(dft_control_type), POINTER :: dft_control
231 : TYPE(scf_control_type), POINTER :: scf_control
232 : TYPE(section_vals_type), POINTER :: dft_section, input, scf_section
233 :
234 2 : CALL get_qs_env(qs_env, input=input, dft_control=dft_control)
235 :
236 2 : CALL get_qs_env(qs_env, scf_control=scf_control)
237 2 : dft_section => section_vals_get_subs_vals(input, "DFT")
238 2 : scf_section => section_vals_get_subs_vals(dft_section, "SCF")
239 :
240 2 : CALL qs_scf_ensure_scf_env(qs_env, scf_env)
241 :
242 2 : CALL section_vals_val_get(scf_section, "CHOLESKY", i_val=scf_env%cholesky_method)
243 2 : scf_control%use_diag = .TRUE.
244 2 : scf_control%diagonalization%method = diag_standard
245 :
246 2 : CALL qs_scf_ensure_mos(qs_env)
247 :
248 : ! set flags for diagonalization
249 : CALL qs_scf_ensure_diagonalization(scf_env, scf_section, qs_env, &
250 2 : scf_control, qs_env%has_unit_metric)
251 2 : CALL qs_scf_ensure_work_matrices(qs_env, scf_env)
252 :
253 2 : CALL init_scf_run(scf_env, qs_env, scf_section, scf_control)
254 :
255 2 : END SUBROUTINE qs_scf_env_init_basic
256 :
257 : ! **************************************************************************************************
258 : !> \brief makes sure scf_env is allocated (might already be from before)
259 : !> in case it is present the g-space mixing storage is reset
260 : !> \param qs_env ...
261 : !> \param scf_env ...
262 : ! **************************************************************************************************
263 21251 : SUBROUTINE qs_scf_ensure_scf_env(qs_env, scf_env)
264 : TYPE(qs_environment_type), POINTER :: qs_env
265 : TYPE(qs_scf_env_type), POINTER :: scf_env
266 :
267 21251 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
268 : TYPE(qs_rho_type), POINTER :: rho
269 :
270 21251 : NULLIFY (rho_g)
271 :
272 27844 : IF (.NOT. ASSOCIATED(scf_env)) THEN ! i.e. for MD this is associated on the second step (it so seems)
273 6593 : ALLOCATE (scf_env)
274 6593 : CALL scf_env_create(scf_env)
275 : ELSE
276 : ! Reallocate mixing store, if the g space grid (cell) has changed
277 14746 : SELECT CASE (scf_env%mixing_method)
278 : CASE (kerker_mix, pulay_mix, broy_mix, multisec_mix)
279 14658 : IF (ASSOCIATED(scf_env%mixing_store)) THEN
280 : ! The current mixing_store data structure does not allow for an unique
281 : ! grid comparison, but the probability that the 1d lengths of the old and
282 : ! the new grid are accidentily equal is rather low
283 88 : CALL get_qs_env(qs_env, rho=rho)
284 88 : CALL qs_rho_get(rho, rho_g=rho_g)
285 88 : IF (ASSOCIATED(scf_env%mixing_store%rhoin)) THEN
286 52 : IF (SIZE(rho_g(1)%pw_grid%gsq) /= SIZE(scf_env%mixing_store%rhoin(1)%cc)) THEN
287 0 : CALL mixing_storage_release(scf_env%mixing_store)
288 0 : DEALLOCATE (scf_env%mixing_store)
289 : END IF
290 : END IF
291 : END IF
292 : END SELECT
293 : END IF
294 :
295 21251 : END SUBROUTINE qs_scf_ensure_scf_env
296 :
297 : ! **************************************************************************************************
298 : !> \brief performs allocation of outer SCF variables
299 : !> \param scf_env the SCF environment which contains the outer SCF variables
300 : !> \param scf_control control settings for the outer SCF loop
301 : !> \param nvar (optional) set number of outer SCF variables externally if CDFT SCF is active
302 : ! **************************************************************************************************
303 21249 : SUBROUTINE qs_scf_ensure_outer_loop_vars(scf_env, scf_control, nvar)
304 : TYPE(qs_scf_env_type), POINTER :: scf_env
305 : TYPE(scf_control_type), POINTER :: scf_control
306 : INTEGER, OPTIONAL :: nvar
307 :
308 : INTEGER :: nhistory, nvariables
309 :
310 21249 : IF (scf_control%outer_scf%have_scf) THEN
311 4139 : nhistory = scf_control%outer_scf%max_scf + 1
312 4139 : IF (PRESENT(nvar)) THEN
313 326 : IF (nvar > 0) THEN
314 : nvariables = nvar
315 : ELSE
316 0 : nvariables = outer_loop_variables_count(scf_control)
317 : END IF
318 : ELSE
319 3813 : nvariables = outer_loop_variables_count(scf_control)
320 : END IF
321 16556 : ALLOCATE (scf_env%outer_scf%variables(nvariables, nhistory))
322 12417 : ALLOCATE (scf_env%outer_scf%count(nhistory))
323 77593 : scf_env%outer_scf%count = 0
324 12417 : ALLOCATE (scf_env%outer_scf%gradient(nvariables, nhistory))
325 12417 : ALLOCATE (scf_env%outer_scf%energy(nhistory))
326 : END IF
327 :
328 21249 : END SUBROUTINE qs_scf_ensure_outer_loop_vars
329 :
330 : ! **************************************************************************************************
331 : !> \brief performs allocation of CDFT SCF variables
332 : !> \param qs_env the qs_env where to perform the allocation
333 : !> \param scf_env the currently active scf_env
334 : !> \param dft_control the dft_control that holds the cdft_control type
335 : !> \param scf_control the currently active scf_control
336 : ! **************************************************************************************************
337 326 : SUBROUTINE qs_scf_ensure_cdft_loop_vars(qs_env, scf_env, dft_control, scf_control)
338 : TYPE(qs_environment_type), POINTER :: qs_env
339 : TYPE(qs_scf_env_type), POINTER :: scf_env
340 : TYPE(dft_control_type), POINTER :: dft_control
341 : TYPE(scf_control_type), POINTER :: scf_control
342 :
343 : INTEGER :: nhistory, nvariables
344 : LOGICAL :: do_kpoints
345 326 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: gradient_history, outer_scf_history, &
346 326 : variable_history
347 :
348 326 : NULLIFY (outer_scf_history, gradient_history, variable_history)
349 326 : CALL get_qs_env(qs_env=qs_env, do_kpoints=do_kpoints)
350 : ! Test kpoints
351 326 : IF (do_kpoints) &
352 0 : CPABORT("CDFT calculation not possible with kpoints")
353 : ! Check that OUTER_SCF section in DFT&SCF is active
354 : ! This section must always be active to facilitate
355 : ! switching of the CDFT and SCF control parameters in outer_loop_switch
356 326 : IF (.NOT. scf_control%outer_scf%have_scf) &
357 0 : CPABORT("Section SCF&OUTER_SCF must be active for CDFT calculations.")
358 : ! Initialize CDFT and outer_loop variables (constraint settings active in scf_control)
359 326 : IF (dft_control%qs_control%cdft_control%constraint_control%have_scf) THEN
360 326 : nhistory = dft_control%qs_control%cdft_control%constraint_control%max_scf + 1
361 326 : IF (scf_control%outer_scf%type /= outer_scf_none) THEN
362 : nvariables = outer_loop_variables_count(scf_control, &
363 62 : dft_control%qs_control%cdft_control)
364 : ELSE
365 : ! First iteration: scf_control has not yet been updated
366 264 : nvariables = SIZE(dft_control%qs_control%cdft_control%target)
367 : END IF
368 1304 : ALLOCATE (dft_control%qs_control%cdft_control%constraint%variables(nvariables, nhistory))
369 978 : ALLOCATE (dft_control%qs_control%cdft_control%constraint%count(nhistory))
370 2246 : dft_control%qs_control%cdft_control%constraint%count = 0
371 978 : ALLOCATE (dft_control%qs_control%cdft_control%constraint%gradient(nvariables, nhistory))
372 978 : ALLOCATE (dft_control%qs_control%cdft_control%constraint%energy(nhistory))
373 326 : CALL qs_scf_ensure_outer_loop_vars(scf_env, scf_control, nvariables)
374 : END IF
375 : ! Executed only on first call (OT settings active in scf_control)
376 : ! Save OT settings and constraint initial values in CDFT control
377 : ! Then switch to constraint outer_scf settings for proper initialization of history
378 326 : IF (scf_control%outer_scf%have_scf) THEN
379 326 : IF (scf_control%outer_scf%type == outer_scf_none) THEN
380 264 : dft_control%qs_control%cdft_control%ot_control%have_scf = .TRUE.
381 264 : dft_control%qs_control%cdft_control%ot_control%max_scf = scf_control%outer_scf%max_scf
382 264 : dft_control%qs_control%cdft_control%ot_control%eps_scf = scf_control%outer_scf%eps_scf
383 264 : dft_control%qs_control%cdft_control%ot_control%step_size = scf_control%outer_scf%step_size
384 264 : dft_control%qs_control%cdft_control%ot_control%type = scf_control%outer_scf%type
385 264 : dft_control%qs_control%cdft_control%ot_control%optimizer = scf_control%outer_scf%optimizer
386 264 : dft_control%qs_control%cdft_control%ot_control%diis_buffer_length = scf_control%outer_scf%diis_buffer_length
387 264 : dft_control%qs_control%cdft_control%ot_control%bisect_trust_count = scf_control%outer_scf%bisect_trust_count
388 : CALL cdft_opt_type_copy(dft_control%qs_control%cdft_control%ot_control%cdft_opt_control, &
389 264 : scf_control%outer_scf%cdft_opt_control)
390 : ! In case constraint and OT extrapolation orders are different, make sure to use former
391 264 : nvariables = SIZE(dft_control%qs_control%cdft_control%target)
392 : IF (scf_control%outer_scf%extrapolation_order /= &
393 : dft_control%qs_control%cdft_control%constraint_control%extrapolation_order &
394 264 : .OR. nvariables /= 1) THEN
395 256 : DEALLOCATE (qs_env%outer_scf_history)
396 256 : DEALLOCATE (qs_env%gradient_history)
397 256 : DEALLOCATE (qs_env%variable_history)
398 256 : nhistory = dft_control%qs_control%cdft_control%constraint_control%extrapolation_order
399 1024 : ALLOCATE (outer_scf_history(nvariables, nhistory))
400 768 : ALLOCATE (gradient_history(nvariables, 2))
401 1324 : gradient_history = 0.0_dp
402 512 : ALLOCATE (variable_history(nvariables, 2))
403 1324 : variable_history = 0.0_dp
404 : CALL set_qs_env(qs_env, outer_scf_history=outer_scf_history, &
405 256 : gradient_history=gradient_history, variable_history=variable_history)
406 : END IF
407 264 : CALL outer_loop_switch(scf_env, scf_control, dft_control%qs_control%cdft_control, ot2cdft)
408 : END IF
409 : END IF
410 :
411 326 : END SUBROUTINE qs_scf_ensure_cdft_loop_vars
412 :
413 : ! **************************************************************************************************
414 : !> \brief performs allocation of the mixing storage
415 : !> \param qs_env ...
416 : !> \param scf_env ...
417 : ! **************************************************************************************************
418 21249 : SUBROUTINE qs_scf_ensure_mixing_store(qs_env, scf_env)
419 : TYPE(qs_environment_type), POINTER :: qs_env
420 : TYPE(qs_scf_env_type), POINTER :: scf_env
421 :
422 : TYPE(dft_control_type), POINTER :: dft_control
423 :
424 21249 : NULLIFY (dft_control)
425 21249 : CALL get_qs_env(qs_env=qs_env, dft_control=dft_control)
426 :
427 21249 : IF (scf_env%mixing_method > 0) THEN
428 : CALL mixing_allocate(qs_env, scf_env%mixing_method, scf_env%p_mix_new, &
429 : scf_env%p_delta, dft_control%nspins, &
430 15132 : scf_env%mixing_store)
431 : ELSE
432 6117 : NULLIFY (scf_env%p_mix_new)
433 : END IF
434 :
435 21249 : END SUBROUTINE qs_scf_ensure_mixing_store
436 :
437 : ! **************************************************************************************************
438 : !> \brief Performs allocation of the SCF work matrices
439 : !> In case of kpoints we probably don't need most of these matrices,
440 : !> maybe we have to initialize some matrices in the fm_pool in kpoints
441 : !> \param qs_env ...
442 : !> \param scf_env ...
443 : ! **************************************************************************************************
444 63753 : SUBROUTINE qs_scf_ensure_work_matrices(qs_env, scf_env)
445 :
446 : TYPE(qs_environment_type), POINTER :: qs_env
447 : TYPE(qs_scf_env_type), POINTER :: scf_env
448 :
449 : CHARACTER(LEN=*), PARAMETER :: routineN = 'qs_scf_ensure_work_matrices'
450 :
451 : INTEGER :: handle, is, nao, nrow_block, nw
452 : LOGICAL :: do_kpoints
453 21251 : TYPE(cp_fm_pool_p_type), DIMENSION(:), POINTER :: ao_mo_fm_pools
454 : TYPE(cp_fm_struct_type), POINTER :: ao_ao_fmstruct, ao_mo_fmstruct
455 21251 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s
456 : TYPE(dbcsr_type), POINTER :: ref_matrix
457 : TYPE(dft_control_type), POINTER :: dft_control
458 21251 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
459 : TYPE(scf_control_type), POINTER :: scf_control
460 :
461 21251 : CALL timeset(routineN, handle)
462 :
463 21251 : NULLIFY (ao_mo_fm_pools, ao_mo_fmstruct, ao_ao_fmstruct, dft_control, matrix_s, mos)
464 :
465 : CALL get_qs_env(qs_env=qs_env, &
466 : dft_control=dft_control, &
467 : matrix_s_kp=matrix_s, &
468 : mos=mos, &
469 : scf_control=scf_control, &
470 21251 : do_kpoints=do_kpoints)
471 21251 : CALL mpools_get(qs_env%mpools, ao_mo_fm_pools=ao_mo_fm_pools)
472 :
473 : ! create an ao_ao parallel matrix structure
474 21251 : ao_mo_fmstruct => fm_pool_get_el_struct(ao_mo_fm_pools(1)%pool)
475 21251 : CALL cp_fm_struct_get(ao_mo_fmstruct, nrow_block=nrow_block)
476 21251 : CALL get_mo_set(mos(1), nao=nao)
477 : CALL cp_fm_struct_create(fmstruct=ao_ao_fmstruct, &
478 : nrow_block=nrow_block, &
479 : ncol_block=nrow_block, &
480 : nrow_global=nao, &
481 : ncol_global=nao, &
482 21251 : template_fmstruct=ao_mo_fmstruct)
483 :
484 21251 : IF ((scf_env%method /= ot_method_nr) .AND. &
485 : (scf_env%method /= block_davidson_diag_method_nr)) THEN
486 15118 : IF (.NOT. ASSOCIATED(scf_env%scf_work1)) THEN
487 13670 : nw = dft_control%nspins
488 13670 : IF (do_kpoints) nw = 4
489 58244 : ALLOCATE (scf_env%scf_work1(nw))
490 30904 : DO is = 1, SIZE(scf_env%scf_work1)
491 : CALL cp_fm_create(scf_env%scf_work1(is), &
492 : matrix_struct=ao_ao_fmstruct, &
493 30904 : name="SCF-WORK_MATRIX-1-"//TRIM(ADJUSTL(cp_to_string(is))))
494 : END DO
495 : END IF
496 : IF ((.NOT. ASSOCIATED(scf_env%ortho)) .AND. &
497 15118 : (scf_env%method /= ot_diag_method_nr) .AND. &
498 : (scf_env%method /= special_diag_method_nr)) THEN
499 : ! Initialize fm matrix to store the Cholesky decomposition
500 11006 : ALLOCATE (scf_env%ortho)
501 : CALL cp_fm_create(scf_env%ortho, &
502 : matrix_struct=ao_ao_fmstruct, &
503 11006 : name="SCF-ORTHO_MATRIX")
504 : ! Initialize dbcsr matrix to store the Cholesky decomposition
505 11006 : IF (scf_env%cholesky_method == cholesky_dbcsr) THEN
506 58 : ref_matrix => matrix_s(1, 1)%matrix
507 58 : CALL dbcsr_init_p(scf_env%ortho_dbcsr)
508 : CALL dbcsr_create(scf_env%ortho_dbcsr, template=ref_matrix, &
509 58 : matrix_type=dbcsr_type_no_symmetry)
510 58 : CALL dbcsr_init_p(scf_env%buf1_dbcsr)
511 : CALL dbcsr_create(scf_env%buf1_dbcsr, template=ref_matrix, &
512 58 : matrix_type=dbcsr_type_no_symmetry)
513 58 : CALL dbcsr_init_p(scf_env%buf2_dbcsr)
514 : CALL dbcsr_create(scf_env%buf2_dbcsr, template=ref_matrix, &
515 58 : matrix_type=dbcsr_type_no_symmetry)
516 10948 : ELSE IF (scf_env%cholesky_method == cholesky_inverse .OR. &
517 : (scf_control%level_shift /= 0.0_dp .AND. &
518 : scf_env%cholesky_method == cholesky_off)) THEN
519 52 : ALLOCATE (scf_env%ortho_m1)
520 : CALL cp_fm_create(scf_env%ortho_m1, &
521 : matrix_struct=ao_ao_fmstruct, &
522 52 : name="SCF-ORTHO_MATRIX-1")
523 : END IF
524 : END IF
525 15118 : IF (.NOT. ASSOCIATED(scf_env%scf_work2)) THEN
526 13670 : ALLOCATE (scf_env%scf_work2)
527 : CALL cp_fm_create(scf_env%scf_work2, &
528 : matrix_struct=ao_ao_fmstruct, &
529 13670 : name="SCF-WORK_MATRIX-2")
530 : END IF
531 : END IF
532 :
533 21251 : IF (dft_control%dft_plus_u) THEN
534 92 : IF (dft_control%plus_u_method_id == plus_u_lowdin) THEN
535 14 : IF (.NOT. ASSOCIATED(scf_env%s_half)) THEN
536 10 : ALLOCATE (scf_env%s_half)
537 : CALL cp_fm_create(scf_env%s_half, &
538 : matrix_struct=ao_ao_fmstruct, &
539 10 : name="S**(1/2) MATRIX")
540 : END IF
541 : END IF
542 : END IF
543 :
544 21251 : IF (do_kpoints) THEN
545 974 : IF (.NOT. ASSOCIATED(scf_env%scf_work1)) THEN
546 0 : nw = 4
547 0 : ALLOCATE (scf_env%scf_work1(nw))
548 0 : DO is = 1, SIZE(scf_env%scf_work1)
549 : CALL cp_fm_create(scf_env%scf_work1(is), &
550 : matrix_struct=ao_ao_fmstruct, &
551 0 : name="SCF-WORK_MATRIX-1-"//TRIM(ADJUSTL(cp_to_string(is))))
552 : END DO
553 : END IF
554 : END IF
555 :
556 21251 : CALL cp_fm_struct_release(ao_ao_fmstruct)
557 :
558 21251 : CALL timestop(handle)
559 :
560 21251 : END SUBROUTINE qs_scf_ensure_work_matrices
561 :
562 : ! **************************************************************************************************
563 : !> \brief performs allocation of the MO matrices
564 : !> \param qs_env ...
565 : ! **************************************************************************************************
566 21251 : SUBROUTINE qs_scf_ensure_mos(qs_env)
567 : TYPE(qs_environment_type), POINTER :: qs_env
568 :
569 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_scf_ensure_mos'
570 :
571 : INTEGER :: handle, ic, ik, ikk, ispin, nmo, nmo_mat
572 21251 : TYPE(cp_fm_pool_p_type), DIMENSION(:), POINTER :: ao_mo_fm_pools
573 : TYPE(cp_fm_type), POINTER :: mo_coeff, mo_coeff_last
574 21251 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mo_derivs
575 21251 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s
576 : TYPE(dbcsr_type), POINTER :: mo_coeff_b
577 : TYPE(dft_control_type), POINTER :: dft_control
578 : TYPE(kpoint_type), POINTER :: kpoints
579 21251 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos, mos_last_converged
580 21251 : TYPE(mo_set_type), DIMENSION(:, :), POINTER :: mos_k
581 : TYPE(xas_environment_type), POINTER :: xas_env
582 :
583 21251 : CALL timeset(routineN, handle)
584 :
585 21251 : NULLIFY (ao_mo_fm_pools, dft_control, mos, xas_env, matrix_s, mos_last_converged, mo_coeff_last)
586 :
587 : CALL get_qs_env(qs_env=qs_env, &
588 : dft_control=dft_control, &
589 : mos=mos, &
590 : matrix_s_kp=matrix_s, &
591 21251 : xas_env=xas_env)
592 21251 : CALL mpools_get(qs_env%mpools, ao_mo_fm_pools=ao_mo_fm_pools)
593 21251 : IF (dft_control%switch_surf_dip) THEN
594 2 : CALL get_qs_env(qs_env, mos_last_converged=mos_last_converged)
595 : END IF
596 :
597 21251 : nmo_mat = dft_control%nspins
598 21251 : IF (dft_control%restricted) nmo_mat = 1 ! right now, there might be more mos than needed derivs
599 :
600 : ! Finish initialization of the MOs
601 21251 : CPASSERT(ASSOCIATED(mos))
602 45070 : DO ispin = 1, SIZE(mos)
603 23819 : CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff, mo_coeff_b=mo_coeff_b)
604 23819 : IF (.NOT. ASSOCIATED(mo_coeff)) THEN
605 : CALL init_mo_set(mos(ispin), &
606 : fm_pool=ao_mo_fm_pools(ispin)%pool, &
607 8084 : name="qs_env%mo"//TRIM(ADJUSTL(cp_to_string(ispin))))
608 : END IF
609 45070 : IF (.NOT. ASSOCIATED(mo_coeff_b)) THEN
610 8084 : CALL cp_fm_get_info(mos(ispin)%mo_coeff, ncol_global=nmo)
611 8084 : CALL dbcsr_init_p(mos(ispin)%mo_coeff_b)
612 : CALL cp_dbcsr_m_by_n_from_row_template(mos(ispin)%mo_coeff_b, template=matrix_s(1, 1)%matrix, n=nmo, &
613 8084 : sym=dbcsr_type_no_symmetry)
614 : END IF
615 : END DO
616 : ! Get the mo_derivs OK if needed
617 21251 : IF (qs_env%requires_mo_derivs) THEN
618 6123 : CALL get_qs_env(qs_env, mo_derivs=mo_derivs)
619 6123 : IF (.NOT. ASSOCIATED(mo_derivs)) THEN
620 9067 : ALLOCATE (mo_derivs(nmo_mat))
621 4833 : DO ispin = 1, nmo_mat
622 2716 : CALL get_mo_set(mos(ispin), mo_coeff_b=mo_coeff_b)
623 2716 : NULLIFY (mo_derivs(ispin)%matrix)
624 2716 : CALL dbcsr_init_p(mo_derivs(ispin)%matrix)
625 : CALL dbcsr_create(mo_derivs(ispin)%matrix, template=mo_coeff_b, &
626 4833 : name="mo_derivs", matrix_type=dbcsr_type_no_symmetry)
627 : END DO
628 2117 : CALL set_qs_env(qs_env, mo_derivs=mo_derivs)
629 : END IF
630 :
631 : ELSE
632 : ! nothing should be done
633 : END IF
634 :
635 : ! Finish initialization of the MOs for ADMM and derivs if needed ***
636 21251 : IF (dft_control%do_admm) THEN
637 910 : IF (dft_control%restricted) CPABORT("ROKS with ADMM is not implemented")
638 : END IF
639 :
640 : ! Finish initialization of mos_last_converged [SGh]
641 21251 : IF (dft_control%switch_surf_dip) THEN
642 2 : CPASSERT(ASSOCIATED(mos_last_converged))
643 4 : DO ispin = 1, SIZE(mos_last_converged)
644 2 : CALL get_mo_set(mos_last_converged(ispin), mo_coeff=mo_coeff_last)
645 4 : IF (.NOT. ASSOCIATED(mo_coeff_last)) THEN
646 : CALL init_mo_set(mos_last_converged(ispin), &
647 : fm_ref=mos(ispin)%mo_coeff, &
648 2 : name="qs_env%mos_last_converged"//TRIM(ADJUSTL(cp_to_string(ispin))))
649 : END IF
650 : END DO
651 : END IF
652 : ! kpoints: we have to initialize all the k-point MOs
653 21251 : CALL get_qs_env(qs_env=qs_env, kpoints=kpoints)
654 21251 : IF (kpoints%nkp /= 0) THEN
655 : ! check for some incompatible options
656 974 : IF (qs_env%requires_mo_derivs) THEN
657 2 : CPWARN("MO derivative methods flag has been switched off for kpoint calculation")
658 : ! we switch it off to make band structure calculations
659 : ! possible for OT gamma point calculations
660 2 : qs_env%requires_mo_derivs = .FALSE.
661 : END IF
662 974 : IF (dft_control%do_xas_calculation) &
663 0 : CPABORT("No XAS implemented with kpoints")
664 974 : IF (qs_env%do_rixs) &
665 0 : CPABORT("RIXS not implemented with kpoints")
666 4134 : DO ik = 1, SIZE(kpoints%kp_env)
667 3160 : CALL mpools_get(kpoints%mpools, ao_mo_fm_pools=ao_mo_fm_pools)
668 3160 : mos_k => kpoints%kp_env(ik)%kpoint_env%mos
669 3160 : ikk = kpoints%kp_range(1) + ik - 1
670 3160 : CPASSERT(ASSOCIATED(mos_k))
671 7776 : DO ispin = 1, SIZE(mos_k, 2)
672 14072 : DO ic = 1, SIZE(mos_k, 1)
673 7270 : CALL get_mo_set(mos_k(ic, ispin), mo_coeff=mo_coeff, mo_coeff_b=mo_coeff_b)
674 7270 : IF (.NOT. ASSOCIATED(mo_coeff)) THEN
675 : CALL init_mo_set(mos_k(ic, ispin), &
676 : fm_pool=ao_mo_fm_pools(ispin)%pool, &
677 : name="kpoints_"//TRIM(ADJUSTL(cp_to_string(ikk)))// &
678 3090 : "%mo"//TRIM(ADJUSTL(cp_to_string(ispin))))
679 : END IF
680 : ! no sparse matrix representation of kpoint MO vectors
681 10912 : CPASSERT(.NOT. ASSOCIATED(mo_coeff_b))
682 : END DO
683 : END DO
684 : END DO
685 : END IF
686 :
687 21251 : CALL timestop(handle)
688 :
689 21251 : END SUBROUTINE qs_scf_ensure_mos
690 :
691 : ! **************************************************************************************************
692 : !> \brief sets flag for mixing/DIIS during scf
693 : !> \param scf_control ...
694 : !> \param scf_section ...
695 : !> \param scf_env ...
696 : !> \param dft_control ...
697 : ! **************************************************************************************************
698 21249 : SUBROUTINE qs_scf_ensure_mixing(scf_control, scf_section, scf_env, dft_control)
699 : TYPE(scf_control_type), POINTER :: scf_control
700 : TYPE(section_vals_type), POINTER :: scf_section
701 : TYPE(qs_scf_env_type), POINTER :: scf_env
702 : TYPE(dft_control_type), POINTER :: dft_control
703 :
704 : TYPE(section_vals_type), POINTER :: mixing_section
705 :
706 21249 : SELECT CASE (scf_control%mixing_method)
707 : CASE (no_mix)
708 0 : scf_env%mixing_method = no_mixing_nr
709 0 : scf_env%p_mix_alpha = 1.0_dp
710 : CASE (direct_p_mix, kerker_mix, pulay_mix, broy_mix, multisec_mix)
711 21249 : scf_env%mixing_method = scf_control%mixing_method
712 21249 : mixing_section => section_vals_get_subs_vals(scf_section, "MIXING")
713 21249 : IF (.NOT. ASSOCIATED(scf_env%mixing_store)) THEN
714 26364 : ALLOCATE (scf_env%mixing_store)
715 : CALL mixing_storage_create(scf_env%mixing_store, mixing_section, scf_env%mixing_method, &
716 6591 : dft_control%qs_control%cutoff)
717 : END IF
718 : CASE DEFAULT
719 21249 : CPABORT("Unknown mixing method")
720 : END SELECT
721 :
722 : ! Disable DIIS for OT and g-space density mixing methods
723 21249 : IF (scf_env%method == ot_method_nr) THEN
724 : ! No mixing is used with OT
725 6117 : scf_env%mixing_method = no_mixing_nr
726 6117 : scf_env%p_mix_alpha = 1.0_dp
727 6117 : scf_env%skip_diis = .TRUE.
728 : END IF
729 :
730 21249 : IF (scf_control%use_diag .AND. scf_env%mixing_method == no_mixing_nr) THEN
731 0 : CPABORT("Diagonalization procedures without mixing are not recommendable")
732 : END IF
733 :
734 21249 : IF (scf_env%mixing_method > direct_mixing_nr) THEN
735 320 : scf_env%skip_diis = .TRUE.
736 320 : scf_env%p_mix_alpha = scf_env%mixing_store%alpha
737 320 : IF (scf_env%mixing_store%beta == 0.0_dp) THEN
738 0 : CPABORT("Mixing employing the Kerker damping factor needs BETA /= 0.0")
739 : END IF
740 : END IF
741 :
742 21249 : IF (scf_env%mixing_method == direct_mixing_nr) THEN
743 14812 : scf_env%p_mix_alpha = scf_env%mixing_store%alpha
744 14812 : IF (scf_control%eps_diis < scf_control%eps_scf) THEN
745 42 : scf_env%skip_diis = .TRUE.
746 42 : CPWARN("the DIIS scheme is disabled, since EPS_DIIS < EPS_SCF")
747 : END IF
748 : END IF
749 :
750 21249 : END SUBROUTINE qs_scf_ensure_mixing
751 :
752 : ! **************************************************************************************************
753 : !> \brief sets flags for diagonalization and ensure that everything is
754 : !> allocated
755 : !> \param scf_env ...
756 : !> \param scf_section ...
757 : !> \param qs_env ...
758 : !> \param scf_control ...
759 : !> \param has_unit_metric ...
760 : ! **************************************************************************************************
761 21251 : SUBROUTINE qs_scf_ensure_diagonalization(scf_env, scf_section, qs_env, &
762 : scf_control, has_unit_metric)
763 : TYPE(qs_scf_env_type), POINTER :: scf_env
764 : TYPE(section_vals_type), POINTER :: scf_section
765 : TYPE(qs_environment_type), POINTER :: qs_env
766 : TYPE(scf_control_type), POINTER :: scf_control
767 : LOGICAL :: has_unit_metric
768 :
769 : INTEGER :: ispin, nao, nmo
770 : LOGICAL :: do_kpoints, need_coeff_b, not_se_or_tb
771 : TYPE(cp_fm_type), POINTER :: mo_coeff
772 : TYPE(dft_control_type), POINTER :: dft_control
773 21251 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
774 :
775 21251 : CALL get_qs_env(qs_env=qs_env, do_kpoints=do_kpoints, dft_control=dft_control, mos=mos)
776 : not_se_or_tb = .NOT. (dft_control%qs_control%dftb .OR. dft_control%qs_control%xtb .OR. &
777 21251 : dft_control%qs_control%semi_empirical)
778 21251 : need_coeff_b = .FALSE.
779 21251 : scf_env%needs_ortho = .FALSE.
780 :
781 21251 : IF (dft_control%smeagol_control%smeagol_enabled .AND. &
782 : dft_control%smeagol_control%run_type == smeagol_runtype_emtransport) THEN
783 0 : scf_env%method = smeagol_method_nr
784 0 : scf_env%skip_diis = .TRUE.
785 0 : scf_control%use_diag = .FALSE.
786 :
787 0 : IF (.NOT. do_kpoints) THEN
788 0 : CPABORT("SMEAGOL requires kpoint calculations")
789 : END IF
790 0 : CPWARN_IF(scf_control%use_ot, "OT is irrelevant to NEGF method")
791 : END IF
792 :
793 21251 : IF (scf_control%use_diag) THEN
794 : ! sanity check whether combinations are allowed
795 15134 : IF (dft_control%restricted) &
796 0 : CPABORT("OT only for restricted (ROKS)")
797 15166 : SELECT CASE (scf_control%diagonalization%method)
798 : CASE (diag_ot, diag_block_krylov, diag_block_davidson)
799 32 : IF (.NOT. not_se_or_tb) &
800 15134 : CPABORT("TB and SE not possible with OT diagonalization")
801 : END SELECT
802 30226 : SELECT CASE (scf_control%diagonalization%method)
803 : ! Diagonalization: additional check whether we are in an orthonormal basis
804 : CASE (diag_standard)
805 15092 : scf_env%method = general_diag_method_nr
806 15092 : scf_env%needs_ortho = (.NOT. has_unit_metric) .AND. (.NOT. do_kpoints)
807 : IF (diag_type == FM_DIAG_TYPE_CUSOLVER .AND. &
808 : cusolver_generalized .AND. &
809 15092 : scf_control%level_shift == 0.0_dp .AND. &
810 : scf_env%cholesky_method /= cholesky_off) THEN
811 0 : scf_env%needs_ortho = .FALSE.
812 : END IF
813 15092 : IF (has_unit_metric) THEN
814 2656 : scf_env%method = special_diag_method_nr
815 : END IF
816 : ! OT Diagonalization: not possible with ROKS
817 : CASE (diag_ot)
818 8 : IF (dft_control%roks) &
819 0 : CPABORT("ROKS with OT diagonalization not possible")
820 8 : IF (do_kpoints) &
821 0 : CPABORT("OT diagonalization not possible with kpoint calculations")
822 8 : scf_env%method = ot_diag_method_nr
823 8 : need_coeff_b = .TRUE.
824 : ! Block Krylov diagonlization: not possible with ROKS,
825 : ! allocation of additional matrices is needed
826 : CASE (diag_block_krylov)
827 8 : IF (dft_control%roks) &
828 0 : CPABORT("ROKS with block PF diagonalization not possible")
829 8 : IF (do_kpoints) &
830 0 : CPABORT("Block Krylov diagonalization not possible with kpoint calculations")
831 8 : scf_env%method = block_krylov_diag_method_nr
832 8 : scf_env%needs_ortho = .TRUE.
833 8 : IF (.NOT. ASSOCIATED(scf_env%krylov_space)) &
834 4 : CALL krylov_space_create(scf_env%krylov_space, scf_section)
835 8 : CALL krylov_space_allocate(scf_env%krylov_space, scf_control, mos)
836 : ! Block davidson diagonlization: allocation of additional matrices is needed
837 : CASE (diag_block_davidson)
838 16 : IF (do_kpoints) &
839 0 : CPABORT("Block Davidson diagonalization not possible with kpoint calculations")
840 16 : scf_env%method = block_davidson_diag_method_nr
841 16 : IF (.NOT. ASSOCIATED(scf_env%block_davidson_env)) &
842 : CALL block_davidson_env_create(scf_env%block_davidson_env, dft_control%nspins, &
843 12 : scf_section)
844 34 : DO ispin = 1, dft_control%nspins
845 18 : CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, nao=nao, nmo=nmo)
846 34 : CALL block_davidson_allocate(scf_env%block_davidson_env(ispin), mo_coeff, nao, nmo)
847 : END DO
848 10 : need_coeff_b = .TRUE.
849 : ! Filter matrix diagonalisation method
850 : CASE (diag_filter_matrix)
851 10 : scf_env%method = filter_matrix_diag_method_nr
852 10 : IF (.NOT. fb_env_has_data(scf_env%filter_matrix_env)) THEN
853 10 : CALL fb_env_create(scf_env%filter_matrix_env)
854 : END IF
855 10 : CALL fb_env_read_input(scf_env%filter_matrix_env, scf_section)
856 10 : CALL fb_env_build_rcut_auto(scf_env%filter_matrix_env, qs_env)
857 10 : CALL fb_env_write_info(scf_env%filter_matrix_env, qs_env, scf_section)
858 10 : CALL fb_distribution_build(scf_env%filter_matrix_env, qs_env, scf_section)
859 10 : CALL fb_env_build_atomic_halos(scf_env%filter_matrix_env, qs_env, scf_section)
860 : CASE DEFAULT
861 15134 : CPABORT("Unknown diagonalization method")
862 : END SELECT
863 : ! Check if subspace diagonlization is requested: allocation of additional matrices is needed
864 15134 : IF (scf_control%do_diag_sub) THEN
865 2 : scf_env%needs_ortho = .TRUE.
866 2 : IF (.NOT. ASSOCIATED(scf_env%subspace_env)) &
867 : CALL diag_subspace_env_create(scf_env%subspace_env, scf_section, &
868 2 : dft_control%qs_control%cutoff)
869 2 : CALL diag_subspace_allocate(scf_env%subspace_env, qs_env, mos)
870 2 : IF (do_kpoints) &
871 0 : CPABORT("No subspace diagonlization with kpoint calculation")
872 : END IF
873 : ! OT: check if OT is used instead of diagonalization. Not possible with added MOS at the moment
874 6117 : ELSEIF (scf_control%use_ot) THEN
875 6117 : scf_env%method = ot_method_nr
876 6117 : need_coeff_b = .TRUE.
877 18351 : IF (SUM(ABS(scf_control%added_mos)) > 0) &
878 0 : CPABORT("OT with ADDED_MOS/=0 not implemented")
879 6117 : IF (dft_control%restricted .AND. dft_control%nspins /= 2) &
880 0 : CPABORT("nspin must be 2 for restricted (ROKS)")
881 6117 : IF (do_kpoints) &
882 0 : CPABORT("OT not possible with kpoint calculations")
883 0 : ELSEIF (scf_env%method /= smeagol_method_nr) THEN
884 0 : CPABORT("OT or DIAGONALIZATION have to be set")
885 : END IF
886 45070 : DO ispin = 1, dft_control%nspins
887 45070 : mos(ispin)%use_mo_coeff_b = need_coeff_b
888 : END DO
889 :
890 21251 : END SUBROUTINE qs_scf_ensure_diagonalization
891 :
892 : ! **************************************************************************************************
893 : !> \brief performs those initialisations that need to be done only once
894 : !> (e.g. that only depend on the atomic positions)
895 : !> this will be called in scf
896 : !> \param scf_env ...
897 : !> \param qs_env ...
898 : !> \param scf_section ...
899 : !> \param scf_control ...
900 : !> \par History
901 : !> 03.2006 created [Joost VandeVondele]
902 : ! **************************************************************************************************
903 21251 : SUBROUTINE init_scf_run(scf_env, qs_env, scf_section, scf_control)
904 :
905 : TYPE(qs_scf_env_type), POINTER :: scf_env
906 : TYPE(qs_environment_type), POINTER :: qs_env
907 : TYPE(section_vals_type), POINTER :: scf_section
908 : TYPE(scf_control_type), POINTER :: scf_control
909 :
910 : CHARACTER(LEN=*), PARAMETER :: routineN = 'init_scf_run'
911 :
912 : INTEGER :: after, handle, homo, ii, ikind, ispin, &
913 : iw, nao, ndep, needed_evals, nmo, &
914 : output_unit
915 : LOGICAL :: dft_plus_u_atom, do_kpoints, &
916 : init_u_ramping_each_scf, omit_headers, &
917 : s_minus_half_available
918 : REAL(KIND=dp) :: u_ramping
919 21251 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: evals
920 21251 : REAL(KIND=dp), DIMENSION(:), POINTER :: eigenvalues
921 : TYPE(cp_fm_struct_type), POINTER :: fm_struct
922 : TYPE(cp_fm_type) :: evecs, fm_w
923 : TYPE(cp_fm_type), POINTER :: mo_coeff
924 : TYPE(cp_logger_type), POINTER :: logger
925 21251 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
926 21251 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s_kp
927 : TYPE(dft_control_type), POINTER :: dft_control
928 : TYPE(kpoint_type), POINTER :: kpoints
929 21251 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
930 : TYPE(mp_para_env_type), POINTER :: para_env
931 21251 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
932 : TYPE(qs_kind_type), POINTER :: qs_kind
933 : TYPE(qs_rho_type), POINTER :: rho
934 : TYPE(xas_environment_type), POINTER :: xas_env
935 :
936 21251 : CALL timeset(routineN, handle)
937 :
938 21251 : NULLIFY (qs_kind_set, matrix_s, dft_control, mos, qs_kind, rho, xas_env, mo_coeff)
939 :
940 21251 : logger => cp_get_default_logger()
941 :
942 21251 : CPASSERT(ASSOCIATED(scf_env))
943 21251 : CPASSERT(ASSOCIATED(qs_env))
944 21251 : NULLIFY (para_env)
945 :
946 21251 : s_minus_half_available = .FALSE.
947 : CALL get_qs_env(qs_env, &
948 : dft_control=dft_control, &
949 : qs_kind_set=qs_kind_set, &
950 : mos=mos, &
951 : rho=rho, &
952 : nelectron_total=scf_env%nelectron, &
953 : do_kpoints=do_kpoints, &
954 : para_env=para_env, &
955 21251 : xas_env=xas_env)
956 :
957 : !Check restricted optimizers available for tblite library
958 21251 : IF (dft_control%qs_control%xtb_control%do_tblite) THEN
959 894 : IF (dft_control%lsd) THEN
960 0 : CPABORT("LSD option not compatible with tblite library.")
961 : END IF
962 894 : IF (scf_env%method == ot_method_nr) THEN
963 0 : CPABORT("OT SCF option not compatible with tblite library.")
964 : END IF
965 : END IF
966 :
967 : ! Calculate ortho matrix
968 21251 : ndep = 0
969 21251 : IF (scf_env%needs_ortho) THEN
970 11470 : CALL get_qs_env(qs_env, matrix_s=matrix_s)
971 11470 : CALL copy_dbcsr_to_fm(matrix_s(1)%matrix, scf_env%ortho)
972 11470 : IF (scf_env%cholesky_method > cholesky_off) THEN
973 11422 : CALL cp_fm_cholesky_decompose(scf_env%ortho)
974 11422 : IF (scf_env%cholesky_method == cholesky_dbcsr) THEN
975 58 : CALL cp_fm_triangular_invert(scf_env%ortho)
976 58 : CALL cp_fm_set_all(scf_env%scf_work2, 0.0_dp)
977 58 : CALL cp_fm_to_fm_triangular(scf_env%ortho, scf_env%scf_work2, "U")
978 58 : CALL copy_fm_to_dbcsr(scf_env%scf_work2, scf_env%ortho_dbcsr)
979 11364 : ELSE IF (scf_env%cholesky_method == cholesky_inverse) THEN
980 34 : CALL cp_fm_to_fm(scf_env%ortho, scf_env%ortho_m1)
981 34 : CALL cp_fm_triangular_invert(scf_env%ortho_m1)
982 : END IF
983 : ELSE
984 48 : CALL cp_fm_get_info(scf_env%ortho, ncol_global=nao)
985 144 : ALLOCATE (evals(nao))
986 1908 : evals = 0
987 :
988 48 : CALL cp_fm_create(evecs, scf_env%ortho%matrix_struct)
989 :
990 : ! Perform an EVD
991 48 : CALL choose_eigv_solver(scf_env%ortho, evecs, evals)
992 :
993 : ! Determine the number of neglectable eigenvalues assuming that the eigenvalues are in ascending order
994 : ! (Required by Lapack)
995 : ndep = 0
996 112 : DO ii = 1, nao
997 112 : IF (evals(ii) > scf_control%eps_eigval) THEN
998 48 : ndep = ii - 1
999 48 : EXIT
1000 : END IF
1001 : END DO
1002 48 : needed_evals = nao - ndep
1003 :
1004 : ! Set the eigenvalue of the eigenvectors belonging to the linear subspace to zero
1005 112 : evals(1:ndep) = 0.0_dp
1006 : ! Determine the eigenvalues of the inverse square root
1007 1844 : evals(ndep + 1:nao) = 1.0_dp/SQRT(evals(ndep + 1:nao))
1008 :
1009 : ! Create reduced matrices
1010 48 : NULLIFY (fm_struct)
1011 : CALL cp_fm_struct_create(fm_struct, template_fmstruct=scf_env%ortho%matrix_struct, &
1012 48 : nrow_global=nao, ncol_global=needed_evals)
1013 :
1014 48 : ALLOCATE (scf_env%ortho_red, scf_env%scf_work2_red)
1015 48 : CALL cp_fm_create(scf_env%ortho_red, fm_struct)
1016 48 : CALL cp_fm_create(scf_env%scf_work2_red, fm_struct)
1017 48 : CALL cp_fm_struct_release(fm_struct)
1018 :
1019 48 : IF (scf_control%level_shift /= 0.0_dp) THEN
1020 : CALL cp_fm_struct_create(fm_struct, template_fmstruct=scf_env%ortho%matrix_struct, &
1021 6 : nrow_global=needed_evals, ncol_global=nao)
1022 :
1023 6 : ALLOCATE (scf_env%ortho_m1_red)
1024 6 : CALL cp_fm_create(scf_env%ortho_m1_red, fm_struct)
1025 6 : CALL cp_fm_struct_release(fm_struct)
1026 : END IF
1027 :
1028 206 : ALLOCATE (scf_env%scf_work1_red(SIZE(scf_env%scf_work1)))
1029 110 : DO ispin = 1, SIZE(scf_env%scf_work1)
1030 : CALL cp_fm_struct_create(fm_struct, template_fmstruct=scf_env%ortho%matrix_struct, &
1031 62 : nrow_global=needed_evals, ncol_global=needed_evals)
1032 62 : CALL cp_fm_create(scf_env%scf_work1_red(ispin), fm_struct)
1033 110 : CALL cp_fm_struct_release(fm_struct)
1034 : END DO
1035 :
1036 : ! Scale the eigenvalues and copy them to
1037 48 : CALL cp_fm_to_fm(evecs, scf_env%ortho_red, needed_evals, ndep + 1, 1)
1038 :
1039 48 : IF (scf_control%level_shift /= 0.0_dp) THEN
1040 6 : CALL cp_fm_transpose(scf_env%ortho_red, scf_env%ortho_m1_red)
1041 : END IF
1042 :
1043 48 : CALL cp_fm_column_scale(scf_env%ortho_red, evals(ndep + 1:))
1044 :
1045 : ! Copy the linear dependent columns to the MO sets and set their orbital energies
1046 : ! to a very large value to reduce the probability of occupying them
1047 110 : DO ispin = 1, SIZE(mos)
1048 62 : CALL get_mo_set(mos(ispin), nmo=nmo, mo_coeff=mo_coeff, homo=homo, eigenvalues=eigenvalues)
1049 62 : IF (needed_evals < nmo) THEN
1050 2 : IF (needed_evals < homo) THEN
1051 : CALL cp_abort(__LOCATION__, &
1052 : "The numerical rank of the overlap matrix is lower than the "// &
1053 : "number of orbitals to be occupied! Check the geometry or increase "// &
1054 0 : "EPS_DEFAULT or EPS_PGF_ORB!")
1055 : END IF
1056 : CALL cp_warn(__LOCATION__, &
1057 : "The numerical rank of the overlap matrix is lower than the number of requested MOs! "// &
1058 : "Reduce the number of MOs to the number of available MOs. If necessary, "// &
1059 2 : "request a lower number of MOs or increase EPS_DEFAULT or EPS_PGF_ORB.")
1060 2 : CALL set_mo_set(mos(ispin), nmo=needed_evals)
1061 : END IF
1062 : ! Copy the last columns to mo_coeff if the container is large enough
1063 62 : CALL cp_fm_to_fm(evecs, mo_coeff, MIN(ndep, MAX(0, nmo - needed_evals)), 1, needed_evals + 1)
1064 : ! Set the corresponding eigenvalues to a large value
1065 : ! This prevents their occupation but still keeps the information on them
1066 182 : eigenvalues(needed_evals + 1:MIN(nao, nmo)) = 1.0_dp/scf_control%eps_eigval
1067 : END DO
1068 :
1069 : ! Obtain ortho from (P)DGEMM, skip the linear dependent columns
1070 : CALL parallel_gemm("N", "T", nao, nao, needed_evals, 1.0_dp, scf_env%ortho_red, evecs, &
1071 48 : 0.0_dp, scf_env%ortho, b_first_col=ndep + 1)
1072 :
1073 48 : IF (scf_control%level_shift /= 0.0_dp) THEN
1074 : ! We need SQRT(evals) of the eigenvalues of H, so 1/SQRT(evals) of ortho_red
1075 168 : evals(ndep + 1:nao) = 1.0_dp/evals(ndep + 1:nao)
1076 6 : CALL cp_fm_row_scale(scf_env%ortho_m1_red, evals(ndep + 1:))
1077 :
1078 : CALL parallel_gemm("T", "T", nao, nao, needed_evals, 1.0_dp, scf_env%ortho_m1_red, evecs, &
1079 6 : 0.0_dp, scf_env%ortho_m1, b_first_col=ndep + 1)
1080 : END IF
1081 :
1082 48 : CALL cp_fm_release(evecs)
1083 :
1084 144 : s_minus_half_available = .TRUE.
1085 : END IF
1086 :
1087 11470 : IF (BTEST(cp_print_key_should_output(logger%iter_info, &
1088 : qs_env%input, "DFT%PRINT%AO_MATRICES/ORTHO"), cp_p_file)) THEN
1089 : iw = cp_print_key_unit_nr(logger, qs_env%input, "DFT%PRINT%AO_MATRICES/ORTHO", &
1090 4 : extension=".Log")
1091 4 : CALL section_vals_val_get(qs_env%input, "DFT%PRINT%AO_MATRICES%NDIGITS", i_val=after)
1092 4 : CALL section_vals_val_get(qs_env%input, "DFT%PRINT%AO_MATRICES%OMIT_HEADERS", l_val=omit_headers)
1093 4 : after = MIN(MAX(after, 1), 16)
1094 : CALL write_fm_with_basis_info(scf_env%ortho, 4, after, qs_env, &
1095 4 : para_env, output_unit=iw, omit_headers=omit_headers)
1096 : CALL cp_print_key_finished_output(iw, logger, qs_env%input, &
1097 4 : "DFT%PRINT%AO_MATRICES/ORTHO")
1098 : END IF
1099 : END IF
1100 :
1101 21251 : CALL get_mo_set(mo_set=mos(1), nao=nao)
1102 :
1103 : ! DFT+U methods based on Lowdin charges need S^(1/2)
1104 21251 : IF (dft_control%dft_plus_u) THEN
1105 92 : IF (dft_control%plus_u_method_id == plus_u_lowdin) THEN
1106 14 : IF (do_kpoints) THEN
1107 0 : CALL get_qs_env(qs_env, kpoints=kpoints, matrix_s_kp=matrix_s_kp)
1108 0 : CALL diag_kp_smat(matrix_s_kp, kpoints, scf_env%scf_work1)
1109 : ELSE
1110 14 : CALL get_qs_env(qs_env, matrix_s=matrix_s)
1111 14 : IF (s_minus_half_available) THEN
1112 : CALL cp_dbcsr_sm_fm_multiply(matrix_s(1)%matrix, scf_env%ortho, &
1113 0 : scf_env%s_half, nao)
1114 : ELSE
1115 14 : CALL copy_dbcsr_to_fm(matrix_s(1)%matrix, scf_env%s_half)
1116 14 : CALL cp_fm_create(fm_w, scf_env%s_half%matrix_struct)
1117 14 : CALL cp_fm_power(scf_env%s_half, fm_w, 0.5_dp, scf_control%eps_eigval, ndep)
1118 14 : CALL cp_fm_release(fm_w)
1119 : END IF
1120 : END IF
1121 : END IF
1122 276 : DO ikind = 1, SIZE(qs_kind_set)
1123 184 : qs_kind => qs_kind_set(ikind)
1124 : CALL get_qs_kind(qs_kind=qs_kind, &
1125 : dft_plus_u_atom=dft_plus_u_atom, &
1126 : u_ramping=u_ramping, &
1127 184 : init_u_ramping_each_scf=init_u_ramping_each_scf)
1128 276 : IF (dft_plus_u_atom .AND. (u_ramping /= 0.0_dp)) THEN
1129 24 : IF (init_u_ramping_each_scf) THEN
1130 12 : CALL set_qs_kind(qs_kind=qs_kind, u_minus_j=0.0_dp)
1131 : END IF
1132 : END IF
1133 : END DO
1134 : END IF
1135 :
1136 : ! extrapolate outer loop variables
1137 21251 : IF (scf_control%outer_scf%have_scf) THEN
1138 4141 : CALL outer_loop_extrapolate(qs_env)
1139 : END IF
1140 :
1141 : ! initializes rho and the mos
1142 21251 : IF (ASSOCIATED(qs_env%xas_env)) THEN
1143 : ! if just optimized wfn, e.g. ground state
1144 : ! changes come from a perturbation, e.g., the occupation numbers
1145 : ! it could be generalized for other cases, at the moment used only for core level spectroscopy
1146 : ! initialize the density with the localized mos
1147 82 : CALL xas_initialize_rho(qs_env, scf_env, scf_control)
1148 : ELSE
1149 : CALL scf_env_initial_rho_setup(scf_env, qs_env=qs_env, &
1150 21169 : scf_section=scf_section, scf_control=scf_control)
1151 : END IF
1152 :
1153 : ! Frozen density approximation
1154 21251 : IF (ASSOCIATED(qs_env%wf_history)) THEN
1155 21251 : IF (qs_env%wf_history%interpolation_method_nr == wfi_frozen_method_nr) THEN
1156 12 : IF (.NOT. ASSOCIATED(qs_env%wf_history%past_states(1)%snapshot)) THEN
1157 4 : CALL wfi_update(qs_env%wf_history, qs_env=qs_env, dt=1.0_dp)
1158 4 : ALLOCATE (qs_env%wf_history%past_states(1)%snapshot%rho_frozen)
1159 4 : CALL qs_rho_create(qs_env%wf_history%past_states(1)%snapshot%rho_frozen)
1160 : CALL duplicate_rho_type(rho_input=rho, &
1161 : rho_output=qs_env%wf_history%past_states(1)%snapshot%rho_frozen, &
1162 4 : qs_env=qs_env)
1163 : END IF
1164 : END IF
1165 : END IF
1166 :
1167 : !image charge method, calculate image_matrix if required
1168 21251 : IF (qs_env%qmmm) THEN
1169 3802 : IF (qs_env%qmmm .AND. qs_env%qmmm_env_qm%image_charge) THEN
1170 : CALL conditional_calc_image_matrix(qs_env=qs_env, &
1171 20 : qmmm_env=qs_env%qmmm_env_qm)
1172 : END IF
1173 : END IF
1174 :
1175 : output_unit = cp_print_key_unit_nr(logger, scf_section, "PRINT%PROGRAM_RUN_INFO", &
1176 21251 : extension=".scfLog")
1177 21251 : CALL qs_scf_initial_info(output_unit, mos, dft_control, ndep)
1178 : CALL cp_print_key_finished_output(output_unit, logger, scf_section, &
1179 21251 : "PRINT%PROGRAM_RUN_INFO")
1180 :
1181 21251 : CALL timestop(handle)
1182 :
1183 42502 : END SUBROUTINE init_scf_run
1184 :
1185 : ! **************************************************************************************************
1186 : !> \brief Initializes rho and the mos, so that an scf cycle can start
1187 : !> \param scf_env the scf env in which to do the scf
1188 : !> \param qs_env the qs env the scf_env lives in
1189 : !> \param scf_section ...
1190 : !> \param scf_control ...
1191 : !> \par History
1192 : !> 02.2003 created [fawzi]
1193 : !> \author fawzi
1194 : ! **************************************************************************************************
1195 21169 : SUBROUTINE scf_env_initial_rho_setup(scf_env, qs_env, scf_section, scf_control)
1196 : TYPE(qs_scf_env_type), POINTER :: scf_env
1197 : TYPE(qs_environment_type), POINTER :: qs_env
1198 : TYPE(section_vals_type), POINTER :: scf_section
1199 : TYPE(scf_control_type), POINTER :: scf_control
1200 :
1201 : CHARACTER(len=*), PARAMETER :: routineN = 'scf_env_initial_rho_setup'
1202 :
1203 : INTEGER :: extrapolation_method_nr, handle, ispin, &
1204 : nmo, output_unit
1205 : LOGICAL :: do_harris, orthogonal_wf
1206 : TYPE(cp_fm_type), POINTER :: mo_coeff
1207 : TYPE(cp_logger_type), POINTER :: logger
1208 : TYPE(dft_control_type), POINTER :: dft_control
1209 : TYPE(harris_type), POINTER :: harris_env
1210 21169 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
1211 : TYPE(mp_para_env_type), POINTER :: para_env
1212 : TYPE(qs_rho_type), POINTER :: rho
1213 21169 : TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom
1214 :
1215 21169 : CALL timeset(routineN, handle)
1216 21169 : NULLIFY (mo_coeff, rho, dft_control, para_env, mos)
1217 21169 : logger => cp_get_default_logger()
1218 21169 : CPASSERT(ASSOCIATED(scf_env))
1219 21169 : CPASSERT(ASSOCIATED(qs_env))
1220 :
1221 : CALL get_qs_env(qs_env, &
1222 : rho=rho, &
1223 : mos=mos, &
1224 : dft_control=dft_control, &
1225 21169 : para_env=para_env)
1226 :
1227 21169 : do_harris = qs_env%harris_method
1228 :
1229 21169 : extrapolation_method_nr = wfi_use_guess_method_nr
1230 21169 : IF (ASSOCIATED(qs_env%wf_history)) THEN
1231 : CALL wfi_extrapolate(qs_env%wf_history, &
1232 : qs_env=qs_env, dt=1.0_dp, &
1233 : extrapolation_method_nr=extrapolation_method_nr, &
1234 21169 : orthogonal_wf=orthogonal_wf)
1235 : ! wfi_use_guess_method_nr the wavefunctions are not yet initialized
1236 : IF ((.NOT. orthogonal_wf) .AND. &
1237 21169 : (scf_env%method == ot_method_nr) .AND. &
1238 : (.NOT. (extrapolation_method_nr == wfi_use_guess_method_nr))) THEN
1239 0 : DO ispin = 1, SIZE(mos)
1240 0 : CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff, nmo=nmo)
1241 0 : CALL reorthogonalize_vectors(qs_env, v_matrix=mo_coeff, n_col=nmo)
1242 0 : IF (dft_control%hairy_probes .EQV. .TRUE.) THEN
1243 0 : scf_control%smear%do_smear = .FALSE.
1244 : CALL set_mo_occupation(mo_set=mos(ispin), &
1245 0 : smear=scf_control%smear, probe=dft_control%probe)
1246 : ELSE
1247 : CALL set_mo_occupation(mo_set=mos(ispin), &
1248 0 : smear=scf_control%smear)
1249 : END IF
1250 : END DO
1251 : END IF
1252 : END IF
1253 :
1254 21169 : IF (.NOT. do_harris) THEN
1255 : output_unit = cp_print_key_unit_nr(logger, scf_section, "PRINT%PROGRAM_RUN_INFO", &
1256 21153 : extension=".scfLog")
1257 21153 : IF (output_unit > 0) THEN
1258 : WRITE (UNIT=output_unit, FMT="(/,T2,A,I0)") &
1259 : "Extrapolation method: "// &
1260 10759 : TRIM(wfi_get_method_label(extrapolation_method_nr))
1261 10759 : IF (extrapolation_method_nr == wfi_ps_method_nr) THEN
1262 : WRITE (UNIT=output_unit, FMT="(T2,A,I0,A)") &
1263 188 : "Extrapolation order: ", &
1264 376 : MAX((MIN(qs_env%wf_history%memory_depth, qs_env%wf_history%snapshot_count) - 1), 0)
1265 : END IF
1266 : END IF
1267 : CALL cp_print_key_finished_output(output_unit, logger, scf_section, &
1268 21153 : "PRINT%PROGRAM_RUN_INFO")
1269 : END IF
1270 :
1271 : IF (do_harris) THEN
1272 16 : CALL get_qs_env(qs_env, harris_env=harris_env)
1273 16 : CALL harris_density_update(qs_env, harris_env)
1274 16 : CALL qs_rho_update_rho(rho, qs_env=qs_env)
1275 16 : CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.TRUE.)
1276 21153 : ELSE IF (extrapolation_method_nr == wfi_use_guess_method_nr) THEN
1277 6791 : CALL calculate_first_density_matrix(scf_env=scf_env, qs_env=qs_env)
1278 6791 : CALL qs_rho_update_rho(rho, qs_env=qs_env)
1279 6791 : CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.TRUE.)
1280 : END IF
1281 :
1282 : ! Some preparation for the mixing
1283 21169 : IF (scf_env%mixing_method > 1) THEN
1284 314 : IF (dft_control%qs_control%gapw) THEN
1285 40 : CALL get_qs_env(qs_env=qs_env, rho_atom_set=rho_atom)
1286 : CALL mixing_init(scf_env%mixing_method, rho, scf_env%mixing_store, &
1287 40 : para_env, rho_atom=rho_atom)
1288 274 : ELSEIF (dft_control%qs_control%dftb .OR. dft_control%qs_control%xtb) THEN
1289 54 : CALL charge_mixing_init(scf_env%mixing_store)
1290 220 : ELSEIF (dft_control%qs_control%semi_empirical) THEN
1291 0 : CPABORT('SE Code not possible')
1292 : ELSE
1293 : CALL mixing_init(scf_env%mixing_method, rho, scf_env%mixing_store, &
1294 220 : para_env)
1295 : END IF
1296 : END IF
1297 :
1298 44824 : DO ispin = 1, SIZE(mos) !fm->dbcsr
1299 44824 : IF (mos(ispin)%use_mo_coeff_b) THEN
1300 : CALL copy_fm_to_dbcsr(mos(ispin)%mo_coeff, &
1301 7127 : mos(ispin)%mo_coeff_b) !fm->dbcsr
1302 : END IF
1303 : END DO !fm->dbcsr
1304 :
1305 21169 : CALL timestop(handle)
1306 :
1307 21169 : END SUBROUTINE scf_env_initial_rho_setup
1308 :
1309 : END MODULE qs_scf_initialization
|
