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