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