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