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