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