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