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