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