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