Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 : ! **************************************************************************************************
8 : !> \brief Utility routines for qs_scf
9 : ! **************************************************************************************************
10 : MODULE qs_scf_loop_utils
11 : USE cp_control_types, ONLY: dft_control_type,&
12 : hairy_probes_type
13 : USE cp_dbcsr_api, ONLY: dbcsr_copy,&
14 : dbcsr_get_info,&
15 : dbcsr_p_type,&
16 : dbcsr_release,&
17 : dbcsr_release_p,&
18 : dbcsr_scale,&
19 : dbcsr_set,&
20 : dbcsr_type
21 : USE cp_dbcsr_contrib, ONLY: dbcsr_dot,&
22 : dbcsr_get_diag,&
23 : dbcsr_scale_by_vector
24 : USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
25 : copy_fm_to_dbcsr,&
26 : cp_dbcsr_sm_fm_multiply,&
27 : dbcsr_allocate_matrix_set,&
28 : dbcsr_deallocate_matrix_set
29 : USE cp_external_control, ONLY: external_control
30 : USE cp_fm_basic_linalg, ONLY: cp_complex_fm_gemm,&
31 : cp_fm_column_scale,&
32 : cp_fm_scale_and_add
33 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
34 : cp_fm_struct_release,&
35 : cp_fm_struct_type
36 : USE cp_fm_types, ONLY: cp_fm_create,&
37 : cp_fm_get_info,&
38 : cp_fm_get_submatrix,&
39 : cp_fm_release,&
40 : cp_fm_set_submatrix,&
41 : cp_fm_to_fm,&
42 : cp_fm_type
43 : USE cp_log_handling, ONLY: cp_to_string
44 : USE ieee_arithmetic, ONLY: ieee_is_finite
45 : USE input_constants, ONLY: diag_update_method_adiis,&
46 : smear_fermi_dirac,&
47 : smear_gaussian,&
48 : smear_mp,&
49 : smear_mv
50 : USE input_section_types, ONLY: section_vals_get,&
51 : section_vals_get_subs_vals,&
52 : section_vals_type
53 : USE kinds, ONLY: default_string_length,&
54 : dp
55 : USE kpoint_methods, ONLY: kpoint_density_matrices,&
56 : kpoint_density_transform,&
57 : kpoint_set_mo_occupation
58 : USE kpoint_types, ONLY: get_kpoint_info,&
59 : kpoint_env_type,&
60 : kpoint_type
61 : USE message_passing, ONLY: mp_para_env_type
62 : USE qs_density_matrices, ONLY: calculate_density_matrix
63 : USE qs_density_mixing_types, ONLY: broyden_mixing_nr,&
64 : direct_mixing_nr,&
65 : gspace_mixing_nr,&
66 : modified_broyden_mixing_nr,&
67 : multisecant_mixing_nr,&
68 : new_pulay_mixing_nr,&
69 : no_mixing_nr,&
70 : pulay_mixing_nr
71 : USE qs_energy_types, ONLY: qs_energy_type
72 : USE qs_environment_types, ONLY: get_qs_env,&
73 : qs_environment_type
74 : USE qs_fb_env_methods, ONLY: fb_env_do_diag
75 : USE qs_gspace_mixing, ONLY: gspace_mixing
76 : USE qs_kpoint_operators, ONLY: kpoint_operator_context_create,&
77 : kpoint_operator_context_release,&
78 : kpoint_operator_context_type,&
79 : kpoint_operator_get,&
80 : kpoint_operator_get_local,&
81 : kpoint_operator_prepare,&
82 : kpoint_operator_store,&
83 : kpoint_slot_ks
84 : USE qs_kpp1_env_types, ONLY: kpp1_release,&
85 : qs_kpp1_env_type
86 : USE qs_ks_types, ONLY: qs_ks_did_change,&
87 : qs_ks_env_type
88 : USE qs_linres_kernel, ONLY: apply_hxc_kernel_kp
89 : USE qs_mixing_utils, ONLY: self_consistency_check
90 : USE qs_mo_occupation, ONLY: set_mo_occupation
91 : USE qs_mo_types, ONLY: get_mo_set,&
92 : mo_set_type
93 : USE qs_mom_methods, ONLY: do_mom_diag
94 : USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
95 : USE qs_ot, ONLY: &
96 : qs_ot_antihermitian_spectral_norm, qs_ot_density_secant_orbital_overlaps, &
97 : qs_ot_density_secant_projected_hessian, qs_ot_density_tangent, &
98 : qs_ot_finite_rotation_response, qs_ot_fixed_n_energy_gradient, &
99 : qs_ot_fixed_n_multigroup_schur_block, qs_ot_get_orbitals_complex, &
100 : qs_ot_get_orbitals_ref_complex, qs_ot_projected_response_update, qs_ot_symmetric_abs_solve
101 : USE qs_ot_minimizer, ONLY: ot_mermin_response_preparation_needed,&
102 : ot_mermin_secant_curvature,&
103 : ot_mini,&
104 : ot_mini_prepare_gradient
105 : USE qs_ot_scf, ONLY: ot_scf_destroy,&
106 : ot_scf_mini
107 : USE qs_ot_types, ONLY: qs_ot_channel_index,&
108 : qs_ot_check_channel_context,&
109 : qs_ot_physical_secant_type,&
110 : qs_ot_type
111 : USE qs_outer_scf, ONLY: outer_loop_gradient
112 : USE qs_rho_methods, ONLY: qs_rho_update_rho
113 : USE qs_rho_types, ONLY: qs_rho_get,&
114 : qs_rho_type
115 : USE qs_scf_diagonalization, ONLY: &
116 : do_block_davidson_diag, do_block_davidson_diag_kp, do_block_krylov_diag, do_general_diag, &
117 : do_general_diag_kp, do_ot_diag, do_ot_diag_kp, do_roks_diag, do_scf_diag_subspace, &
118 : do_special_diag
119 : USE qs_scf_methods, ONLY: cp_sm_mix,&
120 : scf_env_density_mixing
121 : USE qs_scf_output, ONLY: qs_scf_print_summary
122 : USE qs_scf_types, ONLY: &
123 : block_davidson_diag_method_nr, block_krylov_diag_method_nr, filter_matrix_diag_method_nr, &
124 : general_diag_method_nr, ot_diag_method_nr, ot_method_nr, qs_scf_env_type, &
125 : smeagol_method_nr, special_diag_method_nr
126 : USE scf_control_types, ONLY: scf_control_type,&
127 : smear_type
128 : USE smeagol_interface, ONLY: run_smeagol_emtrans
129 : USE smearing_utils, ONLY: Smearkp,&
130 : Smearkp2,&
131 : smearing_response_weight
132 : USE tblite_interface, ONLY: tb_native_scc_mixer_active
133 : #include "./base/base_uses.f90"
134 :
135 : IMPLICIT NONE
136 :
137 : PRIVATE
138 :
139 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_scf_loop_utils'
140 :
141 : PUBLIC :: qs_scf_set_loop_flags, &
142 : qs_scf_new_mos, qs_scf_new_mos_kp, &
143 : qs_scf_kp_ref_refresh_ready, qs_scf_kp_search_endpoint, &
144 : qs_scf_candidate_density_delta, qs_scf_commit_density_candidate, &
145 : qs_scf_density_mixing, qs_scf_check_inner_exit, &
146 : qs_scf_check_outer_exit, qs_scf_inner_finalize, qs_scf_rho_update
147 :
148 : CONTAINS
149 :
150 : ! **************************************************************************************************
151 : !> \brief identify an accepted OT search endpoint from its iteration label
152 : !> \param method OT iteration label
153 : !> \return true for the primary endpoint of every supported minimizer
154 : ! **************************************************************************************************
155 2692 : PURE FUNCTION qs_scf_kp_search_endpoint(method) RESULT(endpoint)
156 :
157 : CHARACTER(LEN=*), INTENT(IN) :: method
158 : LOGICAL :: endpoint
159 :
160 2692 : SELECT CASE (TRIM(method))
161 : CASE ("OT CG", "OT SD", "OT DIIS", "OT BROY", "OT LBFGS")
162 1228 : endpoint = .TRUE.
163 : CASE DEFAULT
164 2692 : endpoint = .FALSE.
165 : END SELECT
166 :
167 2692 : END FUNCTION qs_scf_kp_search_endpoint
168 :
169 : ! **************************************************************************************************
170 : !> \brief decide whether a mature variable-occupation REF chart should be rebuilt
171 : !> \param accepted_searches accepted search directions in the current REF chart
172 : !> \param residual current OT residual
173 : !> \param eps_scf requested SCF threshold
174 : !> \param base_state whether the current point can start a new search direction
175 : !> \param allow_exit_window whether a prior physical REF rebuild permits an exit-only rebuild
176 : !> \return true inside the asymptotic residual window after sufficient accepted history
177 : ! **************************************************************************************************
178 1486 : PURE FUNCTION qs_scf_kp_ref_refresh_ready( &
179 : accepted_searches, residual, eps_scf, base_state, allow_exit_window) &
180 : RESULT(refresh)
181 :
182 : INTEGER, INTENT(IN) :: accepted_searches
183 : REAL(KIND=dp), INTENT(IN) :: residual, eps_scf
184 : LOGICAL, INTENT(IN) :: base_state
185 : LOGICAL, INTENT(IN), OPTIONAL :: allow_exit_window
186 : LOGICAL :: refresh
187 :
188 : INTEGER, PARAMETER :: min_accepted_searches = 12
189 :
190 : LOGICAL :: exit_window
191 : REAL(KIND=dp) :: eps
192 :
193 1486 : eps = MAX(eps_scf, EPSILON(1.0_dp))
194 1486 : exit_window = .FALSE.
195 1486 : IF (PRESENT(allow_exit_window)) exit_window = allow_exit_window
196 : refresh = base_state .AND. accepted_searches >= min_accepted_searches .AND. &
197 : ((residual > 10.0_dp*eps .AND. residual <= 0.2_dp*SQRT(eps)) .OR. &
198 1486 : (exit_window .AND. residual > 0.25_dp*eps .AND. residual <= 10.0_dp*eps))
199 :
200 1486 : END FUNCTION qs_scf_kp_ref_refresh_ready
201 :
202 : ! **************************************************************************************************
203 : !> \brief computes properties for a given hamiltonian using the current wfn
204 : !> \param scf_env ...
205 : !> \param diis_step ...
206 : !> \param energy_only ...
207 : !> \param just_energy ...
208 : !> \param exit_inner_loop ...
209 : ! **************************************************************************************************
210 25439 : SUBROUTINE qs_scf_set_loop_flags(scf_env, diis_step, &
211 : energy_only, just_energy, exit_inner_loop)
212 :
213 : TYPE(qs_scf_env_type), POINTER :: scf_env
214 : LOGICAL :: diis_step, energy_only, just_energy, &
215 : exit_inner_loop
216 :
217 : ! Some flags needed to be set at the beginning of the loop
218 :
219 25439 : diis_step = .FALSE.
220 25439 : energy_only = .FALSE.
221 25439 : just_energy = .FALSE.
222 :
223 : ! SCF loop, optimisation of the wfn coefficients
224 : ! qs_env%rho%rho_r and qs_env%rho%rho_g should be up to date here
225 :
226 25439 : scf_env%iter_count = 0
227 25439 : scf_env%adiis_check_next = .FALSE.
228 25439 : scf_env%adiis_validated = .FALSE.
229 25439 : scf_env%step_norm = 0.0_dp
230 25439 : scf_env%raw_map_delta = 0.0_dp
231 25439 : scf_env%raw_map_delta_valid = .FALSE.
232 25439 : exit_inner_loop = .FALSE.
233 :
234 25439 : END SUBROUTINE qs_scf_set_loop_flags
235 :
236 : ! **************************************************************************************************
237 : !> \brief takes known energy and derivatives and produces new wfns
238 : !> and or density matrix
239 : !> \param qs_env ...
240 : !> \param scf_env ...
241 : !> \param scf_control ...
242 : !> \param scf_section ...
243 : !> \param diis_step ...
244 : !> \param energy_only ...
245 : !> \param probe ...
246 : ! **************************************************************************************************
247 201249 : SUBROUTINE qs_scf_new_mos(qs_env, scf_env, scf_control, scf_section, diis_step, &
248 : energy_only, probe)
249 : TYPE(qs_environment_type), POINTER :: qs_env
250 : TYPE(qs_scf_env_type), POINTER :: scf_env
251 : TYPE(scf_control_type), POINTER :: scf_control
252 : TYPE(section_vals_type), POINTER :: scf_section
253 : LOGICAL :: diis_step, energy_only
254 : TYPE(hairy_probes_type), DIMENSION(:), OPTIONAL, &
255 : POINTER :: probe
256 :
257 : CHARACTER(LEN=*), PARAMETER :: routineN = 'qs_scf_new_mos'
258 :
259 : INTEGER :: handle, ispin
260 : LOGICAL :: disable_diis, has_unit_metric, &
261 : skip_diag_sub
262 : REAL(KIND=dp) :: saved_eps_diis
263 201249 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s
264 : TYPE(dft_control_type), POINTER :: dft_control
265 201249 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
266 : TYPE(qs_energy_type), POINTER :: energy
267 : TYPE(qs_ks_env_type), POINTER :: ks_env
268 : TYPE(qs_rho_type), POINTER :: rho
269 :
270 201249 : CALL timeset(routineN, handle)
271 :
272 201249 : NULLIFY (energy, ks_env, matrix_ks, matrix_s, rho, mos, dft_control)
273 :
274 : CALL get_qs_env(qs_env=qs_env, &
275 : matrix_s=matrix_s, energy=energy, &
276 : ks_env=ks_env, &
277 : matrix_ks=matrix_ks, rho=rho, mos=mos, &
278 : dft_control=dft_control, &
279 201249 : has_unit_metric=has_unit_metric)
280 201249 : scf_env%iter_param = 0.0_dp
281 :
282 : disable_diis = dft_control%qs_control%xtb_control%do_tblite .AND. &
283 201249 : tb_native_scc_mixer_active(dft_control)
284 : IF (disable_diis) THEN
285 12156 : saved_eps_diis = scf_control%eps_diis
286 12156 : scf_control%eps_diis = 0.0_dp
287 : END IF
288 :
289 : ! transfer total_zeff_corr from qs_env to scf_env only if
290 : ! correct_el_density_dip is switched on [SGh]
291 201249 : IF (dft_control%correct_el_density_dip) THEN
292 40 : scf_env%sum_zeff_corr = qs_env%total_zeff_corr
293 40 : IF (ABS(qs_env%total_zeff_corr) > 0.0_dp) THEN
294 40 : IF (scf_env%method /= general_diag_method_nr) THEN
295 : CALL cp_abort(__LOCATION__, &
296 : "Please use ALGORITHM STANDARD in "// &
297 : "SCF%DIAGONALIZATION if "// &
298 : "CORE_CORRECTION /= 0.0 and "// &
299 0 : "SURFACE_DIPOLE_CORRECTION TRUE ")
300 40 : ELSE IF (dft_control%roks) THEN
301 : CALL cp_abort(__LOCATION__, &
302 : "Combination of "// &
303 : "CORE_CORRECTION /= 0.0 and "// &
304 : "SURFACE_DIPOLE_CORRECTION TRUE "// &
305 0 : "is not implemented with ROKS")
306 40 : ELSE IF (scf_control%diagonalization%mom) THEN
307 : CALL cp_abort(__LOCATION__, &
308 : "Combination of "// &
309 : "CORE_CORRECTION /= 0.0 and "// &
310 : "SURFACE_DIPOLE_CORRECTION TRUE "// &
311 0 : "is not implemented with SCF%MOM")
312 : END IF
313 : END IF
314 : END IF
315 :
316 201249 : SELECT CASE (scf_env%method)
317 : CASE DEFAULT
318 : CALL cp_abort(__LOCATION__, &
319 : "unknown scf method: "// &
320 0 : cp_to_string(scf_env%method))
321 :
322 : ! *************************************************************************
323 : ! Filter matrix diagonalisation: ugly implementation at this point of time
324 : ! *************************************************************************
325 : CASE (filter_matrix_diag_method_nr)
326 :
327 80 : IF (ABS(qs_env%total_zeff_corr) > 0.0_dp) THEN
328 : CALL cp_abort(__LOCATION__, &
329 : "CORE_CORRECTION /= 0.0 plus SURFACE_DIPOLE_CORRECTION TRUE "// &
330 0 : "requires SCF%DIAGONALIZATION: ALGORITHM STANDARD")
331 : END IF
332 : CALL fb_env_do_diag(scf_env%filter_matrix_env, qs_env, &
333 80 : matrix_ks, matrix_s, scf_section, diis_step)
334 :
335 : ! Diagonlization in non orthonormal case
336 : CASE (general_diag_method_nr)
337 100743 : IF (dft_control%roks) THEN
338 : CALL do_roks_diag(scf_env, mos, matrix_ks, matrix_s, &
339 : scf_control, scf_section, diis_step, &
340 610 : has_unit_metric)
341 : ELSE
342 100133 : IF (scf_control%diagonalization%mom) THEN
343 : CALL do_mom_diag(scf_env, mos, matrix_ks, &
344 : matrix_s, scf_control, scf_section, &
345 324 : diis_step)
346 : ELSE
347 99809 : IF (dft_control%hairy_probes .EQV. .TRUE.) THEN
348 : CALL do_general_diag(scf_env, mos, matrix_ks, &
349 : matrix_s, scf_control, scf_section, &
350 : diis_step, &
351 14 : probe)
352 : ELSE
353 : CALL do_general_diag(scf_env, mos, matrix_ks, &
354 : matrix_s, scf_control, scf_section, &
355 99795 : diis_step)
356 : END IF
357 : END IF
358 100133 : IF (scf_control%do_diag_sub) THEN
359 : skip_diag_sub = (scf_env%subspace_env%eps_diag_sub > 0.0_dp) .AND. &
360 10 : (scf_env%iter_count == 1 .OR. scf_env%iter_delta > scf_env%subspace_env%eps_diag_sub)
361 : IF (.NOT. skip_diag_sub) THEN
362 : CALL do_scf_diag_subspace(qs_env, scf_env, scf_env%subspace_env, mos, rho, &
363 10 : ks_env, scf_section, scf_control)
364 : END IF
365 : END IF
366 : END IF
367 : ! Diagonlization in orthonormal case
368 : CASE (special_diag_method_nr)
369 18398 : IF (dft_control%roks) THEN
370 : CALL do_roks_diag(scf_env, mos, matrix_ks, matrix_s, &
371 : scf_control, scf_section, diis_step, &
372 522 : has_unit_metric)
373 : ELSE
374 : CALL do_special_diag(scf_env, mos, matrix_ks, &
375 : scf_control, scf_section, &
376 17876 : diis_step)
377 : END IF
378 : ! OT diagonalization
379 : CASE (ot_diag_method_nr)
380 : CALL do_ot_diag(scf_env, mos, matrix_ks, matrix_s, &
381 202 : scf_control, scf_section, diis_step)
382 : ! Block Krylov diagonlization
383 : CASE (block_krylov_diag_method_nr)
384 114 : IF ((scf_env%krylov_space%eps_std_diag > 0.0_dp) .AND. &
385 : (scf_env%iter_count == 1 .OR. scf_env%iter_delta > scf_env%krylov_space%eps_std_diag)) THEN
386 2 : IF (scf_env%krylov_space%always_check_conv) THEN
387 : CALL do_block_krylov_diag(scf_env, mos, matrix_ks, &
388 0 : scf_control, scf_section, check_moconv_only=.TRUE.)
389 : END IF
390 : CALL do_general_diag(scf_env, mos, matrix_ks, &
391 2 : matrix_s, scf_control, scf_section, diis_step)
392 : ELSE
393 : CALL do_block_krylov_diag(scf_env, mos, matrix_ks, &
394 112 : scf_control, scf_section)
395 : END IF
396 114 : IF (scf_control%do_diag_sub) THEN
397 : skip_diag_sub = (scf_env%subspace_env%eps_diag_sub > 0.0_dp) .AND. &
398 0 : (scf_env%iter_count == 1 .OR. scf_env%iter_delta > scf_env%subspace_env%eps_diag_sub)
399 : IF (.NOT. skip_diag_sub) THEN
400 : CALL do_scf_diag_subspace(qs_env, scf_env, scf_env%subspace_env, mos, rho, &
401 0 : ks_env, scf_section, scf_control)
402 : END IF
403 : END IF
404 : ! Block Davidson diagonlization
405 : CASE (block_davidson_diag_method_nr)
406 : CALL do_block_davidson_diag(qs_env, scf_env, mos, matrix_ks, matrix_s, scf_control, &
407 94 : scf_section, .FALSE.)
408 : ! OT without diagonlization. Needs special treatment for SCP runs
409 : CASE (ot_method_nr)
410 : CALL qs_scf_loop_do_ot(qs_env, scf_env, scf_control%smear, mos, rho, &
411 : qs_env%mo_derivs, energy%total, &
412 201249 : matrix_s, energy_only=energy_only, has_unit_metric=has_unit_metric)
413 : END SELECT
414 201249 : IF (disable_diis) scf_control%eps_diis = saved_eps_diis
415 :
416 201249 : energy%kTS = 0.0_dp
417 201249 : energy%efermi = 0.0_dp
418 201249 : CALL get_qs_env(qs_env, mos=mos)
419 437903 : DO ispin = 1, SIZE(mos)
420 236654 : energy%kTS = energy%kTS + mos(ispin)%kTS
421 437903 : energy%efermi = energy%efermi + mos(ispin)%mu
422 : END DO
423 201249 : energy%efermi = energy%efermi/REAL(SIZE(mos), KIND=dp)
424 :
425 201249 : CALL timestop(handle)
426 :
427 201249 : END SUBROUTINE qs_scf_new_mos
428 :
429 : ! **************************************************************************************************
430 : !> \brief Updates MOs and density matrix using diagonalization
431 : !> Kpoint code
432 : !> \param qs_env ...
433 : !> \param scf_env ...
434 : !> \param scf_control ...
435 : !> \param diis_step ...
436 : !> \param probe ...
437 : !> \param ot_kp_subspace_refresh ...
438 : !> \param allow_ot_kp_subspace_refresh ...
439 : !> \param allow_ot_kp_exit_refresh ...
440 : !> \param accepted_ot_kp_searches ...
441 : !> \param added_mos_auto_grow ...
442 : !> \param energy_only ...
443 : ! **************************************************************************************************
444 37306 : SUBROUTINE qs_scf_new_mos_kp(qs_env, scf_env, scf_control, diis_step, probe, &
445 : ot_kp_subspace_refresh, allow_ot_kp_subspace_refresh, &
446 : allow_ot_kp_exit_refresh, accepted_ot_kp_searches, &
447 : added_mos_auto_grow, energy_only)
448 : TYPE(qs_environment_type), POINTER :: qs_env
449 : TYPE(qs_scf_env_type), POINTER :: scf_env
450 : TYPE(scf_control_type), POINTER :: scf_control
451 : LOGICAL :: diis_step
452 : TYPE(hairy_probes_type), DIMENSION(:), OPTIONAL, &
453 : POINTER :: probe
454 : LOGICAL, INTENT(OUT), OPTIONAL :: ot_kp_subspace_refresh
455 : LOGICAL, INTENT(IN), OPTIONAL :: allow_ot_kp_subspace_refresh, &
456 : allow_ot_kp_exit_refresh
457 : INTEGER, INTENT(IN), OPTIONAL :: accepted_ot_kp_searches
458 : LOGICAL, INTENT(OUT), OPTIONAL :: added_mos_auto_grow, energy_only
459 :
460 : CHARACTER(LEN=*), PARAMETER :: routineN = 'qs_scf_new_mos_kp'
461 :
462 : INTEGER :: accepted_searches, handle, ispin, &
463 : refresh_count
464 : LOGICAL :: allow_exit_refresh, allow_refresh, base_state, disable_diis, has_unit_metric, &
465 : my_added_mos_auto_grow, ot_energy_only, refresh
466 : REAL(dp) :: diis_error, residual, saved_eps_diis
467 37306 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks, matrix_s, matrix_t
468 : TYPE(dft_control_type), POINTER :: dft_control
469 : TYPE(kpoint_type), POINTER :: kpoints
470 37306 : TYPE(mo_set_type), DIMENSION(:, :), POINTER :: mos
471 : TYPE(mp_para_env_type), POINTER :: para_env
472 : TYPE(qs_energy_type), POINTER :: energy
473 :
474 37306 : CALL timeset(routineN, handle)
475 :
476 37306 : allow_refresh = .FALSE.
477 37306 : allow_exit_refresh = .FALSE.
478 37306 : accepted_searches = 0
479 37306 : my_added_mos_auto_grow = .FALSE.
480 37306 : refresh = .FALSE.
481 37306 : IF (PRESENT(allow_ot_kp_subspace_refresh)) allow_refresh = allow_ot_kp_subspace_refresh
482 37306 : IF (PRESENT(allow_ot_kp_exit_refresh)) allow_exit_refresh = allow_ot_kp_exit_refresh
483 37306 : IF (PRESENT(accepted_ot_kp_searches)) accepted_searches = accepted_ot_kp_searches
484 37306 : IF (PRESENT(ot_kp_subspace_refresh)) ot_kp_subspace_refresh = .FALSE.
485 37306 : IF (PRESENT(added_mos_auto_grow)) added_mos_auto_grow = .FALSE.
486 37306 : IF (PRESENT(energy_only)) energy_only = .FALSE.
487 :
488 37306 : NULLIFY (dft_control, energy, kpoints, matrix_ks, matrix_s, matrix_t, para_env)
489 :
490 37306 : CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, energy=energy, kpoints=kpoints)
491 37306 : CALL get_kpoint_info(kpoints, para_env=para_env)
492 37306 : CPASSERT(ASSOCIATED(para_env))
493 37306 : scf_env%iter_param = 0.0_dp
494 : disable_diis = dft_control%qs_control%xtb_control%do_tblite .AND. &
495 37306 : tb_native_scc_mixer_active(dft_control)
496 : IF (disable_diis) THEN
497 10474 : saved_eps_diis = scf_control%eps_diis
498 10474 : scf_control%eps_diis = 0.0_dp
499 : END IF
500 :
501 37306 : IF (dft_control%roks .AND. scf_env%method /= ot_method_nr) THEN
502 0 : CPABORT("KP code: ROKS method not available: ")
503 : END IF
504 :
505 37306 : SELECT CASE (scf_env%method)
506 : CASE DEFAULT
507 : CALL cp_abort(__LOCATION__, &
508 : "KP code: Unknown scf method: "// &
509 0 : cp_to_string(scf_env%method))
510 : CASE (general_diag_method_nr)
511 : ! Diagonlization in non orthonormal case
512 33578 : CALL get_qs_env(qs_env, matrix_ks_kp=matrix_ks, matrix_s_kp=matrix_s)
513 33578 : IF (scf_control%diagonalization%update_method == diag_update_method_adiis) THEN
514 80 : CPASSERT(ASSOCIATED(scf_env%scf_subspace_buffer))
515 80 : diis_error = scf_env%scf_subspace_buffer%last_old_fock_weight
516 : END IF
517 33578 : IF (dft_control%hairy_probes .EQV. .TRUE.) THEN
518 0 : scf_control%smear%do_smear = .FALSE.
519 : CALL do_general_diag_kp(matrix_ks, matrix_s, kpoints, scf_env, scf_control, .TRUE., &
520 0 : diis_step, diis_error, qs_env, probe)
521 : ELSE
522 : CALL do_general_diag_kp(matrix_ks, matrix_s, kpoints, scf_env, scf_control, .TRUE., &
523 : diis_step, diis_error, qs_env, &
524 33578 : added_mos_auto_grow=my_added_mos_auto_grow)
525 : END IF
526 33578 : IF (my_added_mos_auto_grow) THEN
527 0 : IF (PRESENT(added_mos_auto_grow)) added_mos_auto_grow = .TRUE.
528 0 : IF (disable_diis) scf_control%eps_diis = saved_eps_diis
529 0 : CALL timestop(handle)
530 10 : RETURN
531 : END IF
532 33578 : IF (scf_control%diagonalization%update_method == diag_update_method_adiis .AND. &
533 0 : .NOT. diis_step) THEN
534 48 : scf_env%iter_param = diis_error
535 48 : IF (scf_env%scf_subspace_buffer%last_restart) THEN
536 0 : scf_env%iter_method = "ADIIS/Rst."
537 : ELSE
538 48 : scf_env%iter_method = "ADIIS/Diag."
539 : END IF
540 33530 : ELSE IF (diis_step) THEN
541 10688 : scf_env%iter_param = diis_error
542 10688 : scf_env%iter_method = "DIIS/Diag."
543 : ELSE
544 22842 : IF (scf_env%mixing_method == 0) THEN
545 0 : scf_env%iter_method = "NoMix/Diag."
546 22842 : ELSE IF (scf_env%mixing_method == 1) THEN
547 20096 : scf_env%iter_param = scf_env%p_mix_alpha
548 20096 : scf_env%iter_method = "P_Mix/Diag."
549 2746 : ELSE IF (scf_env%mixing_method > 1) THEN
550 2746 : scf_env%iter_param = scf_env%mixing_store%alpha
551 2746 : scf_env%iter_method = TRIM(scf_env%mixing_store%iter_method)//"/Diag."
552 : END IF
553 : END IF
554 : CASE (special_diag_method_nr)
555 0 : CALL get_qs_env(qs_env=qs_env, has_unit_metric=has_unit_metric)
556 0 : CPASSERT(has_unit_metric)
557 : ! Diagonlization in orthonormal case
558 : CALL cp_abort(__LOCATION__, &
559 : "KP code: Scf method not available: "// &
560 0 : cp_to_string(scf_env%method))
561 : CASE (ot_diag_method_nr)
562 98 : CALL get_qs_env(qs_env, matrix_ks_kp=matrix_ks, matrix_s_kp=matrix_s, kinetic_kp=matrix_t)
563 : CALL do_ot_diag_kp(matrix_ks, matrix_s, matrix_t, kpoints, scf_env, scf_control, &
564 98 : diis_step, my_added_mos_auto_grow)
565 98 : IF (my_added_mos_auto_grow) THEN
566 0 : IF (PRESENT(added_mos_auto_grow)) added_mos_auto_grow = .TRUE.
567 0 : IF (disable_diis) scf_control%eps_diis = saved_eps_diis
568 0 : CALL timestop(handle)
569 0 : RETURN
570 : END IF
571 : CASE (block_krylov_diag_method_nr)
572 : CALL cp_abort(__LOCATION__, &
573 : "KP code: Scf method not available: "// &
574 0 : cp_to_string(scf_env%method))
575 : CASE (block_davidson_diag_method_nr)
576 936 : CALL get_qs_env(qs_env, matrix_ks_kp=matrix_ks, matrix_s_kp=matrix_s)
577 : CALL do_block_davidson_diag_kp(qs_env, scf_env, matrix_ks, matrix_s, kpoints, &
578 936 : scf_control, added_mos_auto_grow=my_added_mos_auto_grow)
579 936 : IF (my_added_mos_auto_grow) THEN
580 4 : IF (PRESENT(added_mos_auto_grow)) added_mos_auto_grow = .TRUE.
581 4 : IF (disable_diis) scf_control%eps_diis = saved_eps_diis
582 4 : CALL timestop(handle)
583 4 : RETURN
584 : END IF
585 : CASE (ot_method_nr)
586 2694 : CALL get_qs_env(qs_env, matrix_ks_kp=matrix_ks, matrix_s_kp=matrix_s)
587 : base_state = .NOT. scf_env%qs_ot_env(1)%energy_only .AND. &
588 : (scf_env%qs_ot_env(1)%line_search_count == 0 .OR. &
589 2694 : scf_env%qs_ot_env(1)%line_search_might_be_done)
590 2694 : residual = MAX(scf_env%iter_delta, scf_env%qs_ot_env(1)%delta)
591 : refresh = (allow_refresh .OR. allow_exit_refresh) .AND. &
592 : scf_env%qs_ot_env(1)%settings%do_ener .AND. &
593 : scf_env%qs_ot_env(1)%settings%occupation_preconditioner .AND. &
594 : (dft_control%smear .OR. scf_control%smear%do_smear) .AND. &
595 : qs_scf_kp_ref_refresh_ready(accepted_searches, residual, &
596 2694 : scf_control%eps_scf, base_state, allow_exit_refresh)
597 2694 : refresh_count = MERGE(1, 0, refresh)
598 2694 : CALL para_env%sum(refresh_count)
599 2694 : refresh = refresh_count > 0
600 2694 : IF (refresh) THEN
601 24 : IF (PRESENT(ot_kp_subspace_refresh)) ot_kp_subspace_refresh = .TRUE.
602 : ELSE
603 : CALL qs_scf_loop_do_ot_kp(qs_env, scf_env, matrix_ks, matrix_s, &
604 2670 : my_added_mos_auto_grow, ot_energy_only)
605 2670 : IF (PRESENT(energy_only)) energy_only = ot_energy_only
606 2670 : IF (my_added_mos_auto_grow) THEN
607 6 : IF (PRESENT(added_mos_auto_grow)) added_mos_auto_grow = .TRUE.
608 6 : IF (disable_diis) scf_control%eps_diis = saved_eps_diis
609 6 : CALL timestop(handle)
610 6 : RETURN
611 : END IF
612 2664 : IF (allow_exit_refresh .AND. scf_env%iter_delta <= scf_control%eps_scf) THEN
613 14 : base_state = qs_scf_kp_search_endpoint(scf_env%iter_method)
614 : refresh = qs_scf_kp_ref_refresh_ready( &
615 : accepted_searches, scf_env%iter_delta, scf_control%eps_scf, &
616 14 : base_state, allow_exit_window=.TRUE.)
617 14 : IF (refresh .AND. PRESENT(ot_kp_subspace_refresh)) THEN
618 0 : ot_kp_subspace_refresh = .TRUE.
619 : END IF
620 : END IF
621 2664 : refresh_count = MERGE(1, 0, refresh)
622 2664 : CALL para_env%sum(refresh_count)
623 2664 : refresh = refresh_count > 0
624 2664 : IF (refresh .AND. PRESENT(ot_kp_subspace_refresh)) THEN
625 0 : ot_kp_subspace_refresh = .TRUE.
626 : END IF
627 : END IF
628 : CASE (smeagol_method_nr)
629 : ! SMEAGOL interface
630 0 : diis_step = .FALSE.
631 0 : IF (scf_env%mixing_method == 0) THEN
632 0 : scf_env%iter_method = "NoMix/SMGL"
633 0 : ELSE IF (scf_env%mixing_method == 1) THEN
634 0 : scf_env%iter_param = scf_env%p_mix_alpha
635 0 : scf_env%iter_method = "P_Mix/SMGL"
636 0 : ELSE IF (scf_env%mixing_method > 1) THEN
637 0 : scf_env%iter_param = scf_env%mixing_store%alpha
638 0 : scf_env%iter_method = TRIM(scf_env%mixing_store%iter_method)//"/SMGL"
639 : END IF
640 38340 : CALL run_smeagol_emtrans(qs_env, last=.FALSE., iter=scf_env%iter_count, rho_ao_kp=scf_env%p_mix_new)
641 : END SELECT
642 37296 : IF (disable_diis) scf_control%eps_diis = saved_eps_diis
643 :
644 37296 : CALL get_qs_env(qs_env=qs_env, energy=energy)
645 37296 : energy%kTS = 0.0_dp
646 37296 : energy%efermi = 0.0_dp
647 37296 : mos => kpoints%kp_env(1)%kpoint_env%mos
648 76856 : DO ispin = 1, SIZE(mos, 2)
649 39560 : energy%kTS = energy%kTS + mos(1, ispin)%kTS
650 76856 : energy%efermi = energy%efermi + mos(1, ispin)%mu
651 : END DO
652 37296 : energy%efermi = energy%efermi/REAL(SIZE(mos, 2), KIND=dp)
653 :
654 37296 : CALL timestop(handle)
655 :
656 37306 : END SUBROUTINE qs_scf_new_mos_kp
657 :
658 : ! **************************************************************************************************
659 : !> \brief K-point OT driver for local complex spin/k-point channels
660 : !> \param qs_env ...
661 : !> \param scf_env ...
662 : !> \param matrix_ks ...
663 : !> \param matrix_s ...
664 : !> \param added_mos_auto_grow ...
665 : !> \param energy_only ...
666 : ! **************************************************************************************************
667 2670 : SUBROUTINE qs_scf_loop_do_ot_kp(qs_env, scf_env, matrix_ks, matrix_s, added_mos_auto_grow, energy_only)
668 :
669 : TYPE(qs_environment_type), POINTER :: qs_env
670 : TYPE(qs_scf_env_type), POINTER :: scf_env
671 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks, matrix_s
672 : LOGICAL, INTENT(OUT) :: added_mos_auto_grow, energy_only
673 :
674 : CHARACTER(LEN=*), PARAMETER :: routineN = 'qs_scf_loop_do_ot_kp'
675 :
676 : INTEGER :: energy_spin, energy_start, first_channel, handle, homo, icomponent, ikpoint, &
677 : ispin, local_channel, local_kpoint, nao, nkp_groups, nkpoint, nlocal_channels, nmo, &
678 : nspin_energy, nspin_ot
679 : INTEGER, DIMENSION(2) :: kp_range
680 : LOGICAL :: ot_gradient_prepared, &
681 : physical_base_state, restricted_roks, &
682 : use_real_wfn
683 : REAL(KIND=dp) :: trial_kTS
684 2670 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: orbital_scaling_factor, roks_occupation, &
685 2670 : trial_mu
686 2670 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: trial_eigenvalues, trial_occupation
687 2670 : REAL(KIND=dp), DIMENSION(:), POINTER :: eigenvalues, occupation_numbers, wkp
688 : TYPE(cp_fm_struct_type), POINTER :: active_mo_struct, chc_struct
689 : TYPE(cp_fm_type) :: active_mo_coeff, active_mo_coeff_im, &
690 : chc_im, chc_re, hc_im, hc_re, &
691 : hc_roks_im, hc_roks_re
692 2670 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: fmwork
693 : TYPE(cp_fm_type), POINTER :: mo_coeff, mo_coeff_im, mo_coeff_target
694 2670 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_hc, matrix_hc_im, &
695 2670 : matrix_hc_physical, &
696 2670 : matrix_hc_physical_im
697 2670 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: density_target, rho_ao_kp
698 : TYPE(dbcsr_type), POINTER :: matrix_k_im, matrix_k_re, matrix_s_im, &
699 : matrix_s_re
700 : TYPE(dft_control_type), POINTER :: dft_control
701 : TYPE(kpoint_env_type), POINTER :: kp
702 : TYPE(kpoint_type), POINTER :: kpoints
703 : TYPE(mp_para_env_type), POINTER :: para_env_inter_kp
704 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
705 2670 : POINTER :: sab_nl
706 : TYPE(qs_energy_type), POINTER :: energy
707 2670 : TYPE(qs_ot_type), DIMENSION(:), POINTER :: local_ot_env
708 : TYPE(qs_rho_type), POINTER :: rho
709 : TYPE(scf_control_type), POINTER :: scf_control
710 :
711 2670 : CALL timeset(routineN, handle)
712 2670 : added_mos_auto_grow = .FALSE.
713 2670 : ot_gradient_prepared = .FALSE.
714 :
715 2670 : NULLIFY (active_mo_struct, chc_struct, density_target, dft_control, energy, fmwork, kpoints, &
716 2670 : local_ot_env, matrix_hc, matrix_hc_im, matrix_k_im, matrix_k_re, matrix_s_im, &
717 2670 : matrix_hc_physical, matrix_hc_physical_im, matrix_s_re, mo_coeff, mo_coeff_im, &
718 2670 : occupation_numbers, para_env_inter_kp, rho, &
719 2670 : rho_ao_kp, sab_nl, scf_control, wkp)
720 :
721 2670 : CPASSERT(ASSOCIATED(scf_env%qs_ot_env))
722 : CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, energy=energy, kpoints=kpoints, rho=rho, &
723 2670 : scf_control=scf_control)
724 : CALL get_kpoint_info(kpoints, nkp=nkpoint, wkp=wkp, use_real_wfn=use_real_wfn, &
725 : kp_range=kp_range, nkp_groups=nkp_groups, para_env_inter_kp=para_env_inter_kp, &
726 2670 : sab_nl=sab_nl)
727 2670 : fmwork => scf_env%scf_work1
728 2670 : CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp)
729 2670 : CPASSERT(ASSOCIATED(sab_nl))
730 2670 : CPASSERT(ASSOCIATED(matrix_ks))
731 2670 : CPASSERT(ASSOCIATED(matrix_s))
732 :
733 2670 : IF (use_real_wfn) THEN
734 : CALL cp_abort(__LOCATION__, &
735 0 : "K-point OT driver currently supports complex k-point wavefunctions only.")
736 : END IF
737 2670 : IF (scf_env%qs_ot_env(1)%settings%ot_algorithm /= "TOD" .AND. &
738 : scf_env%qs_ot_env(1)%settings%ot_algorithm /= "REF") THEN
739 : CALL cp_abort(__LOCATION__, &
740 0 : "K-point OT driver supports OT%ALGORITHM STRICT or IRAC.")
741 : END IF
742 2670 : kp => kpoints%kp_env(1)%kpoint_env
743 2670 : restricted_roks = dft_control%restricted
744 2670 : nspin_ot = MERGE(1, SIZE(kp%mos, 2), restricted_roks)
745 2670 : nspin_energy = MERGE(SIZE(kp%mos, 2), nspin_ot, restricted_roks)
746 : CALL qs_ot_check_channel_context(scf_env%qs_ot_env, nspin_ot, &
747 : nkpoint=nkpoint, &
748 : restricted=ANY(scf_env%qs_ot_env(:)%restricted), &
749 : require_kpoint=.TRUE., &
750 : kp_range=kp_range, &
751 : wkp=wkp, &
752 : require_local_state=.TRUE., &
753 8172 : require_complex_state=.TRUE.)
754 :
755 2670 : IF (nkp_groups > 1) THEN
756 942 : CALL kpoint_operator_store(kpoints, fmwork(1), matrix_ks, matrix_s)
757 : END IF
758 :
759 2670 : first_channel = qs_ot_channel_index(1, kp_range(1), nspin_ot)
760 2670 : nlocal_channels = SIZE(kpoints%kp_env)*nspin_ot
761 2670 : local_ot_env => scf_env%qs_ot_env(first_channel:first_channel + nlocal_channels - 1)
762 2670 : local_ot_env(1)%etotal = energy%total
763 : physical_base_state = local_ot_env(1)%settings%do_ener .AND. &
764 : .NOT. local_ot_env(1)%energy_only .AND. &
765 : (local_ot_env(1)%line_search_count == 0 .OR. &
766 2670 : local_ot_env(1)%line_search_might_be_done)
767 :
768 12289 : ALLOCATE (matrix_hc(nlocal_channels))
769 9619 : ALLOCATE (matrix_hc_im(nlocal_channels))
770 2670 : IF (local_ot_env(1)%settings%occupation_preconditioner) THEN
771 4808 : ALLOCATE (matrix_hc_physical(nlocal_channels))
772 4808 : ALLOCATE (matrix_hc_physical_im(nlocal_channels))
773 : END IF
774 6949 : DO local_channel = 1, nlocal_channels
775 4279 : NULLIFY (matrix_hc(local_channel)%matrix)
776 4279 : NULLIFY (matrix_hc_im(local_channel)%matrix)
777 4279 : ALLOCATE (matrix_hc(local_channel)%matrix)
778 4279 : ALLOCATE (matrix_hc_im(local_channel)%matrix)
779 4279 : CALL dbcsr_copy(matrix_hc(local_channel)%matrix, local_ot_env(local_channel)%matrix_x, name="matrix_hc")
780 : CALL dbcsr_copy(matrix_hc_im(local_channel)%matrix, local_ot_env(local_channel)%matrix_x_im, &
781 4279 : name="matrix_hc_im")
782 4279 : CALL dbcsr_set(matrix_hc(local_channel)%matrix, 0.0_dp)
783 4279 : CALL dbcsr_set(matrix_hc_im(local_channel)%matrix, 0.0_dp)
784 6949 : IF (local_ot_env(1)%settings%occupation_preconditioner) THEN
785 1948 : NULLIFY (matrix_hc_physical(local_channel)%matrix)
786 1948 : NULLIFY (matrix_hc_physical_im(local_channel)%matrix)
787 1948 : ALLOCATE (matrix_hc_physical(local_channel)%matrix)
788 1948 : ALLOCATE (matrix_hc_physical_im(local_channel)%matrix)
789 : CALL dbcsr_copy(matrix_hc_physical(local_channel)%matrix, &
790 1948 : local_ot_env(local_channel)%matrix_x, name="matrix_hc_physical")
791 : CALL dbcsr_copy(matrix_hc_physical_im(local_channel)%matrix, &
792 1948 : local_ot_env(local_channel)%matrix_x_im, name="matrix_hc_physical_im")
793 1948 : CALL dbcsr_set(matrix_hc_physical(local_channel)%matrix, 0.0_dp)
794 1948 : CALL dbcsr_set(matrix_hc_physical_im(local_channel)%matrix, 0.0_dp)
795 : END IF
796 : END DO
797 :
798 6643 : DO local_kpoint = 1, SIZE(kpoints%kp_env)
799 3973 : kp => kpoints%kp_env(local_kpoint)%kpoint_env
800 3973 : ikpoint = kp%nkpoint
801 3973 : CPASSERT(ASSOCIATED(kp%ot_hmat))
802 3973 : CPASSERT(SIZE(kp%ot_hmat, 1) >= 2)
803 10922 : DO ispin = 1, nspin_ot
804 4279 : local_channel = qs_ot_channel_index(ispin, ikpoint, nspin_ot) - first_channel + 1
805 : CALL kpoint_operator_get_local( &
806 : matrix_ks, kpoints, kp, MIN(ispin, SIZE(matrix_ks, 1)), &
807 : kp%ot_hmat(1, MIN(ispin, SIZE(kp%ot_hmat, 2))), &
808 4279 : kp%ot_hmat(2, MIN(ispin, SIZE(kp%ot_hmat, 2))), matrix_k_re, matrix_k_im)
809 4279 : IF (physical_base_state .AND. nkp_groups == 1) THEN
810 : ! The direct single-group transform bypasses the full-matrix cache. Synchronize it
811 : ! at accepted endpoints so the physical Mermin history sees the H(k) used for H(k)C(k).
812 : CALL copy_dbcsr_to_fm( &
813 888 : matrix_k_re, kp%ot_hmat(1, MIN(ispin, SIZE(kp%ot_hmat, 2))))
814 : CALL copy_dbcsr_to_fm( &
815 888 : matrix_k_im, kp%ot_hmat(2, MIN(ispin, SIZE(kp%ot_hmat, 2))))
816 : END IF
817 :
818 : CALL get_mo_set(kp%mos(1, ispin), homo=homo, mo_coeff=mo_coeff, nao=nao, nmo=nmo, &
819 4279 : occupation_numbers=occupation_numbers)
820 4279 : CALL get_mo_set(kp%mos(2, ispin), mo_coeff=mo_coeff_im)
821 4279 : IF (local_ot_env(local_channel)%settings%do_ener) homo = nmo
822 4279 : CPASSERT(homo >= 1 .AND. homo <= nmo)
823 4279 : IF (.NOT. local_ot_env(local_channel)%settings%do_ener) THEN
824 18639 : IF (MAXVAL(ABS(occupation_numbers(1:homo) - occupation_numbers(1))) > &
825 : 100.0_dp*EPSILON(1.0_dp)) THEN
826 : CALL cp_abort(__LOCATION__, &
827 0 : "Fixed-occupation K-point OT requires uniform occupations within each channel.")
828 : END IF
829 2179 : IF (occupation_numbers(1) <= 100.0_dp*EPSILON(1.0_dp)) THEN
830 : CALL cp_abort(__LOCATION__, &
831 0 : "Fixed-occupation K-point OT encountered an empty active orbital block.")
832 : END IF
833 2179 : IF (homo < nmo) THEN
834 332 : IF (MAXVAL(ABS(occupation_numbers(homo + 1:nmo))) > 100.0_dp*EPSILON(1.0_dp)) THEN
835 : CALL cp_abort(__LOCATION__, &
836 0 : "Fixed-occupation K-point OT cannot optimize occupied padded MO columns.")
837 : END IF
838 : END IF
839 : END IF
840 : CALL cp_fm_struct_create(active_mo_struct, template_fmstruct=mo_coeff%matrix_struct, &
841 4279 : ncol_global=homo)
842 4279 : CALL cp_fm_create(active_mo_coeff, active_mo_struct)
843 4279 : CALL cp_fm_create(active_mo_coeff_im, active_mo_struct)
844 4279 : CALL cp_fm_to_fm(mo_coeff, active_mo_coeff, homo)
845 4279 : CALL cp_fm_to_fm(mo_coeff_im, active_mo_coeff_im, homo)
846 4279 : CALL cp_fm_create(hc_re, active_mo_struct)
847 4279 : CALL cp_fm_create(hc_im, active_mo_struct)
848 :
849 4279 : CALL cp_dbcsr_sm_fm_multiply(matrix_k_re, active_mo_coeff, hc_re, homo)
850 : CALL cp_dbcsr_sm_fm_multiply(matrix_k_im, active_mo_coeff_im, hc_re, homo, &
851 4279 : alpha=-1.0_dp, beta=1.0_dp)
852 4279 : CALL cp_dbcsr_sm_fm_multiply(matrix_k_re, active_mo_coeff_im, hc_im, homo)
853 : CALL cp_dbcsr_sm_fm_multiply(matrix_k_im, active_mo_coeff, hc_im, homo, &
854 4279 : alpha=1.0_dp, beta=1.0_dp)
855 :
856 4279 : IF (restricted_roks) THEN
857 : ! ROKS has one complex orbital variable but two spin densities. Form the exact
858 : ! shared-orbital derivative before handing it to the one-channel OT minimizer.
859 137 : IF (local_ot_env(local_channel)%settings%do_ener) THEN
860 : CALL cp_fm_struct_create(chc_struct, template_fmstruct=active_mo_struct, &
861 50 : nrow_global=homo, ncol_global=homo)
862 50 : CALL cp_fm_create(chc_re, chc_struct)
863 50 : CALL cp_fm_create(chc_im, chc_struct)
864 : CALL multiply_complex_fm_adjoint( &
865 : active_mo_coeff, active_mo_coeff_im, hc_re, hc_im, &
866 50 : chc_re, chc_im, 1.0_dp, 0.0_dp)
867 50 : CALL copy_fm_to_dbcsr(chc_re, local_ot_env(local_channel)%rot_mat_chc)
868 : CALL dbcsr_get_diag(local_ot_env(local_channel)%rot_mat_chc, &
869 50 : local_ot_env(local_channel)%ener_rayleigh(1:homo))
870 : CALL local_ot_env(local_channel)%para_env%sum( &
871 1018 : local_ot_env(local_channel)%ener_rayleigh(1:homo))
872 50 : CALL cp_fm_release(chc_im)
873 50 : CALL cp_fm_release(chc_re)
874 50 : CALL cp_fm_struct_release(chc_struct)
875 : END IF
876 137 : CALL cp_fm_column_scale(hc_re, occupation_numbers(1:homo))
877 137 : CALL cp_fm_column_scale(hc_im, occupation_numbers(1:homo))
878 137 : CALL dbcsr_release_p(matrix_k_re)
879 137 : CALL dbcsr_release_p(matrix_k_im)
880 :
881 : CALL kpoint_operator_get_local( &
882 : matrix_ks, kpoints, kp, 2, kp%ot_hmat(1, 2), kp%ot_hmat(2, 2), &
883 137 : matrix_k_re, matrix_k_im)
884 137 : CALL get_mo_set(kp%mos(1, 2), occupation_numbers=occupation_numbers, nmo=nmo)
885 411 : ALLOCATE (roks_occupation(homo))
886 137 : roks_occupation(:) = 0.0_dp
887 795 : roks_occupation(1:MIN(homo, nmo)) = occupation_numbers(1:MIN(homo, nmo))
888 :
889 137 : CALL cp_fm_create(hc_roks_re, active_mo_struct)
890 137 : CALL cp_fm_create(hc_roks_im, active_mo_struct)
891 137 : CALL cp_dbcsr_sm_fm_multiply(matrix_k_re, active_mo_coeff, hc_roks_re, homo)
892 : CALL cp_dbcsr_sm_fm_multiply(matrix_k_im, active_mo_coeff_im, hc_roks_re, homo, &
893 137 : alpha=-1.0_dp, beta=1.0_dp)
894 137 : CALL cp_dbcsr_sm_fm_multiply(matrix_k_re, active_mo_coeff_im, hc_roks_im, homo)
895 : CALL cp_dbcsr_sm_fm_multiply(matrix_k_im, active_mo_coeff, hc_roks_im, homo, &
896 137 : alpha=1.0_dp, beta=1.0_dp)
897 137 : IF (local_ot_env(local_channel)%settings%do_ener) THEN
898 : CALL cp_fm_struct_create(chc_struct, template_fmstruct=active_mo_struct, &
899 50 : nrow_global=homo, ncol_global=homo)
900 50 : CALL cp_fm_create(chc_re, chc_struct)
901 50 : CALL cp_fm_create(chc_im, chc_struct)
902 : CALL multiply_complex_fm_adjoint( &
903 : active_mo_coeff, active_mo_coeff_im, hc_roks_re, hc_roks_im, &
904 50 : chc_re, chc_im, 1.0_dp, 0.0_dp)
905 50 : CALL copy_fm_to_dbcsr(chc_re, local_ot_env(local_channel)%rot_mat_chc)
906 : CALL dbcsr_get_diag(local_ot_env(local_channel)%rot_mat_chc, &
907 50 : local_ot_env(local_channel)%ener_rayleigh(homo + 1:2*homo))
908 : CALL local_ot_env(local_channel)%para_env%sum( &
909 1018 : local_ot_env(local_channel)%ener_rayleigh(homo + 1:2*homo))
910 50 : CALL cp_fm_release(chc_im)
911 50 : CALL cp_fm_release(chc_re)
912 50 : CALL cp_fm_struct_release(chc_struct)
913 : END IF
914 137 : CALL cp_fm_column_scale(hc_roks_re, roks_occupation)
915 137 : CALL cp_fm_column_scale(hc_roks_im, roks_occupation)
916 137 : CALL cp_fm_scale_and_add(1.0_dp, hc_re, 1.0_dp, hc_roks_re)
917 137 : CALL cp_fm_scale_and_add(1.0_dp, hc_im, 1.0_dp, hc_roks_im)
918 137 : CALL cp_fm_release(hc_roks_re)
919 137 : CALL cp_fm_release(hc_roks_im)
920 411 : DEALLOCATE (roks_occupation)
921 : END IF
922 :
923 4279 : IF (local_ot_env(local_channel)%settings%do_ener) THEN
924 : CALL cp_fm_struct_create(chc_struct, template_fmstruct=active_mo_struct, &
925 2100 : nrow_global=homo, ncol_global=homo)
926 2100 : CALL cp_fm_create(chc_re, chc_struct)
927 2100 : CALL cp_fm_create(chc_im, chc_struct)
928 : CALL multiply_complex_fm_adjoint( &
929 : active_mo_coeff, active_mo_coeff_im, hc_re, hc_im, &
930 2100 : chc_re, chc_im, 1.0_dp, 0.0_dp)
931 2100 : CALL copy_fm_to_dbcsr(chc_re, local_ot_env(local_channel)%rot_mat_chc)
932 2100 : CPASSERT(ASSOCIATED(local_ot_env(local_channel)%rot_mat_chc_im))
933 2100 : CALL copy_fm_to_dbcsr(chc_im, local_ot_env(local_channel)%rot_mat_chc_im)
934 2100 : IF (.NOT. restricted_roks) THEN
935 : CALL dbcsr_get_diag(local_ot_env(local_channel)%rot_mat_chc, &
936 2050 : local_ot_env(local_channel)%ener_rayleigh)
937 2050 : CPASSERT(ASSOCIATED(local_ot_env(local_channel)%para_env))
938 : CALL local_ot_env(local_channel)%para_env%sum( &
939 29962 : local_ot_env(local_channel)%ener_rayleigh)
940 : END IF
941 2100 : CALL cp_fm_release(chc_im)
942 2100 : CALL cp_fm_release(chc_re)
943 2100 : CALL cp_fm_struct_release(chc_struct)
944 : END IF
945 :
946 4279 : CALL copy_fm_to_dbcsr(hc_re, matrix_hc(local_channel)%matrix)
947 4279 : CALL copy_fm_to_dbcsr(hc_im, matrix_hc_im(local_channel)%matrix)
948 :
949 4279 : CALL cp_fm_release(active_mo_coeff)
950 4279 : CALL cp_fm_release(active_mo_coeff_im)
951 4279 : CALL cp_fm_release(hc_re)
952 4279 : CALL cp_fm_release(hc_im)
953 4279 : CALL cp_fm_struct_release(active_mo_struct)
954 4279 : CALL dbcsr_release_p(matrix_k_re)
955 21089 : CALL dbcsr_release_p(matrix_k_im)
956 : END DO
957 : END DO
958 :
959 2670 : IF (local_ot_env(1)%settings%do_ener) THEN
960 : CALL qs_scf_loop_collect_ot_kp_trial_smearing( &
961 : kpoints, local_ot_env, first_channel, nspin_ot, nspin_energy, restricted_roks, &
962 : kp_range, wkp, para_env_inter_kp, scf_control%smear, &
963 1582 : trial_eigenvalues, trial_occupation, trial_mu, trial_kTS)
964 : END IF
965 :
966 6643 : DO local_kpoint = 1, SIZE(kpoints%kp_env)
967 3973 : kp => kpoints%kp_env(local_kpoint)%kpoint_env
968 3973 : ikpoint = kp%nkpoint
969 10922 : DO ispin = 1, nspin_ot
970 4279 : local_channel = qs_ot_channel_index(ispin, ikpoint, nspin_ot) - first_channel + 1
971 : CALL get_mo_set(kp%mos(1, ispin), homo=homo, mo_coeff=mo_coeff, nmo=nmo, &
972 4279 : occupation_numbers=occupation_numbers)
973 4279 : CALL get_mo_set(kp%mos(2, ispin), mo_coeff=mo_coeff_im)
974 4279 : IF (local_ot_env(local_channel)%settings%do_ener) homo = nmo
975 4279 : IF (physical_base_state .AND. &
976 : local_ot_env(local_channel)%settings%occupation_preconditioner) THEN
977 976 : IF (restricted_roks) THEN
978 22 : IF (.NOT. (nspin_energy <= SIZE(local_ot_env(local_channel)%mermin_physical))) THEN
979 0 : CALL cp_abort(__LOCATION__, "Missing Mermin physical spin state")
980 : END IF
981 66 : DO energy_spin = 1, nspin_energy
982 : CALL qs_scf_kp_roll_physical_secant( &
983 : local_ot_env(local_channel)%mermin_physical(energy_spin), &
984 : mo_coeff, mo_coeff_im, kp%ot_hmat(1, energy_spin), &
985 : kp%ot_hmat(2, energy_spin), &
986 66 : trial_occupation(:, ikpoint, energy_spin))
987 : END DO
988 : ELSE
989 : CALL qs_scf_kp_roll_physical_secant( &
990 : local_ot_env(local_channel)%mermin_physical(1), mo_coeff, mo_coeff_im, &
991 : kp%ot_hmat(1, MIN(ispin, SIZE(kp%ot_hmat, 2))), &
992 : kp%ot_hmat(2, MIN(ispin, SIZE(kp%ot_hmat, 2))), &
993 954 : trial_occupation(:, ikpoint, ispin))
994 : END IF
995 : END IF
996 4279 : IF (local_ot_env(local_channel)%settings%occupation_preconditioner) THEN
997 : CALL dbcsr_copy(matrix_hc_physical(local_channel)%matrix, &
998 1948 : matrix_hc(local_channel)%matrix)
999 : CALL dbcsr_copy(matrix_hc_physical_im(local_channel)%matrix, &
1000 1948 : matrix_hc_im(local_channel)%matrix)
1001 5844 : ALLOCATE (orbital_scaling_factor(homo))
1002 1948 : IF (restricted_roks) THEN
1003 534 : orbital_scaling_factor(:) = 2.0_dp*wkp(ikpoint)
1004 : ELSE
1005 : orbital_scaling_factor(:) = 2.0_dp*wkp(ikpoint)* &
1006 12510 : trial_occupation(1:homo, ikpoint, ispin)
1007 : END IF
1008 : CALL dbcsr_scale_by_vector(matrix_hc_physical(local_channel)%matrix, &
1009 1948 : alpha=orbital_scaling_factor, side='right')
1010 : CALL dbcsr_scale_by_vector(matrix_hc_physical_im(local_channel)%matrix, &
1011 1948 : alpha=orbital_scaling_factor, side='right')
1012 1948 : DEALLOCATE (orbital_scaling_factor)
1013 : END IF
1014 12837 : ALLOCATE (orbital_scaling_factor(homo))
1015 4279 : IF (local_ot_env(local_channel)%settings%occupation_preconditioner) THEN
1016 13044 : orbital_scaling_factor(:) = 2.0_dp*wkp(ikpoint)
1017 2331 : ELSE IF (restricted_roks) THEN
1018 261 : orbital_scaling_factor(:) = 2.0_dp*wkp(ikpoint)
1019 2244 : ELSE IF (local_ot_env(local_channel)%settings%do_ener) THEN
1020 : orbital_scaling_factor(:) = 2.0_dp*wkp(ikpoint)* &
1021 3496 : trial_occupation(1:homo, ikpoint, ispin)
1022 : ELSE
1023 18378 : orbital_scaling_factor(:) = 2.0_dp*wkp(ikpoint)*occupation_numbers(1:homo)
1024 : END IF
1025 : CALL dbcsr_scale_by_vector(matrix_hc(local_channel)%matrix, &
1026 4279 : alpha=orbital_scaling_factor, side='right')
1027 : CALL dbcsr_scale_by_vector(matrix_hc_im(local_channel)%matrix, &
1028 4279 : alpha=orbital_scaling_factor, side='right')
1029 12531 : DEALLOCATE (orbital_scaling_factor)
1030 : END DO
1031 : END DO
1032 :
1033 2670 : IF (local_ot_env(1)%settings%do_ener) THEN
1034 1582 : local_ot_env(1)%etotal = energy%total - energy%kTS + trial_kTS
1035 : CALL qs_scf_loop_prepare_ot_kp_energy_gradient( &
1036 : kpoints, local_ot_env, first_channel, nspin_ot, nspin_energy, restricted_roks, &
1037 : wkp, para_env_inter_kp, scf_control%smear, &
1038 1582 : trial_eigenvalues, trial_occupation, trial_mu)
1039 1582 : IF (local_ot_env(1)%settings%occupation_preconditioner) THEN
1040 : CALL ot_mini_prepare_gradient( &
1041 : local_ot_env, matrix_hc, matrix_hc_im=matrix_hc_im, &
1042 : matrix_hc_physical=matrix_hc_physical, &
1043 : matrix_hc_physical_im=matrix_hc_physical_im, &
1044 1430 : para_env_inter_kp=para_env_inter_kp)
1045 : ELSE
1046 : CALL ot_mini_prepare_gradient( &
1047 : local_ot_env, matrix_hc, matrix_hc_im=matrix_hc_im, &
1048 152 : para_env_inter_kp=para_env_inter_kp)
1049 : END IF
1050 1582 : ot_gradient_prepared = .TRUE.
1051 : ! A coupled response direction is consumed only when the minimizer starts a new
1052 : ! search. Line-search trial points still need their physical gradient, but rebuilding
1053 : ! the dense finite-rotation Hessian there cannot affect the accepted direction.
1054 1582 : IF (physical_base_state) THEN
1055 : CALL qs_scf_loop_prepare_ot_kp_rotation_response( &
1056 : qs_env, kpoints, local_ot_env, first_channel, nspin_ot, restricted_roks, &
1057 : wkp, para_env_inter_kp, sab_nl, fmwork, rho_ao_kp, &
1058 828 : scf_control%smear, trial_eigenvalues, trial_occupation, trial_mu)
1059 : END IF
1060 : IF (physical_base_state .AND. &
1061 828 : local_ot_env(1)%settings%occupation_preconditioner .AND. &
1062 : (local_ot_env(1)%settings%ot_method == "CG" .OR. &
1063 : local_ot_env(1)%settings%ot_method == "LBFG")) THEN
1064 856 : DO local_channel = 1, nlocal_channels
1065 1328 : CALL qs_scf_kp_store_mermin_gradient(local_ot_env(local_channel))
1066 : END DO
1067 : END IF
1068 : END IF
1069 2670 : IF (local_ot_env(1)%settings%occupation_preconditioner) THEN
1070 : CALL ot_mini(local_ot_env, matrix_hc, matrix_hc_im=matrix_hc_im, &
1071 : matrix_hc_physical=matrix_hc_physical, &
1072 : matrix_hc_physical_im=matrix_hc_physical_im, &
1073 : para_env_inter_kp=para_env_inter_kp, &
1074 1430 : gradient_prepared=ot_gradient_prepared)
1075 : ELSE
1076 : CALL ot_mini(local_ot_env, matrix_hc, matrix_hc_im=matrix_hc_im, &
1077 : para_env_inter_kp=para_env_inter_kp, &
1078 1240 : gradient_prepared=ot_gradient_prepared)
1079 : END IF
1080 2670 : scf_env%iter_delta = local_ot_env(1)%delta
1081 2670 : energy_only = local_ot_env(1)%energy_only
1082 :
1083 6643 : DO local_kpoint = 1, SIZE(kpoints%kp_env)
1084 3973 : kp => kpoints%kp_env(local_kpoint)%kpoint_env
1085 3973 : ikpoint = kp%nkpoint
1086 3973 : CPASSERT(ASSOCIATED(kp%ot_smat))
1087 3973 : CPASSERT(SIZE(kp%ot_smat) >= 2)
1088 : CALL kpoint_operator_get_local(matrix_s, kpoints, kp, 1, &
1089 : kp%ot_smat(1), kp%ot_smat(2), &
1090 3973 : matrix_s_re, matrix_s_im)
1091 8252 : DO ispin = 1, nspin_ot
1092 4279 : local_channel = qs_ot_channel_index(ispin, ikpoint, nspin_ot) - first_channel + 1
1093 : CALL get_mo_set(kp%mos(1, ispin), homo=homo, mo_coeff=mo_coeff, eigenvalues=eigenvalues, &
1094 4279 : nao=nao, nmo=nmo)
1095 4279 : CALL get_mo_set(kp%mos(2, ispin), mo_coeff=mo_coeff_im)
1096 4279 : IF (local_ot_env(local_channel)%settings%do_ener) homo = nmo
1097 1382 : SELECT CASE (local_ot_env(local_channel)%settings%ot_algorithm)
1098 : CASE ("TOD")
1099 : CALL qs_ot_get_orbitals_complex(matrix_hc(local_channel)%matrix, &
1100 : matrix_hc_im(local_channel)%matrix, &
1101 : matrix_s_re, matrix_s_im, &
1102 1382 : local_ot_env(local_channel))
1103 : CASE ("REF")
1104 : CALL qs_ot_get_orbitals_ref_complex(matrix_hc(local_channel)%matrix, &
1105 : matrix_hc_im(local_channel)%matrix, &
1106 : matrix_s_re, matrix_s_im, &
1107 2897 : local_ot_env(local_channel), local_ot_env(1))
1108 : CASE DEFAULT
1109 4279 : CPABORT("K-point OT orbital update requires ALGORITHM STRICT or IRAC")
1110 : END SELECT
1111 : CALL cp_fm_struct_create(active_mo_struct, template_fmstruct=mo_coeff%matrix_struct, &
1112 4279 : ncol_global=homo)
1113 4279 : CALL cp_fm_create(active_mo_coeff, active_mo_struct)
1114 4279 : CALL cp_fm_create(active_mo_coeff_im, active_mo_struct)
1115 4279 : CALL copy_dbcsr_to_fm(matrix_hc(local_channel)%matrix, active_mo_coeff)
1116 4279 : CALL copy_dbcsr_to_fm(matrix_hc_im(local_channel)%matrix, active_mo_coeff_im)
1117 4279 : CALL cp_fm_to_fm(active_mo_coeff, mo_coeff, homo)
1118 4279 : CALL cp_fm_to_fm(active_mo_coeff_im, mo_coeff_im, homo)
1119 4279 : IF (local_ot_env(local_channel)%settings%do_ener .AND. .NOT. restricted_roks) THEN
1120 16006 : eigenvalues(1:homo) = local_ot_env(local_channel)%ener_x(:)
1121 29962 : kp%mos(2, ispin)%eigenvalues(1:homo) = local_ot_env(local_channel)%ener_x(:)
1122 : END IF
1123 4279 : CALL cp_fm_release(active_mo_coeff)
1124 4279 : CALL cp_fm_release(active_mo_coeff_im)
1125 16810 : CALL cp_fm_struct_release(active_mo_struct)
1126 : END DO
1127 3973 : IF (restricted_roks) THEN
1128 137 : CPASSERT(SIZE(kp%mos, 2) == 2)
1129 411 : DO icomponent = 1, 2
1130 274 : CALL get_mo_set(kp%mos(icomponent, 1), mo_coeff=mo_coeff)
1131 274 : CALL get_mo_set(kp%mos(icomponent, 2), mo_coeff=mo_coeff_target)
1132 411 : CALL cp_fm_to_fm(mo_coeff, mo_coeff_target)
1133 : END DO
1134 137 : IF (local_ot_env(local_channel)%settings%do_ener) THEN
1135 150 : DO energy_spin = 1, nspin_energy
1136 100 : energy_start = (energy_spin - 1)*homo + 1
1137 : kp%mos(1, energy_spin)%eigenvalues(1:homo) = &
1138 2036 : local_ot_env(local_channel)%ener_x(energy_start:energy_start + homo - 1)
1139 : kp%mos(2, energy_spin)%eigenvalues(1:homo) = &
1140 2086 : local_ot_env(local_channel)%ener_x(energy_start:energy_start + homo - 1)
1141 : END DO
1142 : END IF
1143 : END IF
1144 3973 : CALL dbcsr_release_p(matrix_s_re)
1145 6643 : CALL dbcsr_release_p(matrix_s_im)
1146 : END DO
1147 :
1148 2670 : IF (local_ot_env(1)%settings%do_ener) THEN
1149 : CALL kpoint_set_mo_occupation(kpoints, scf_control%smear, &
1150 : added_mos_auto=scf_control%added_mos_auto, &
1151 : added_mos_auto_grow=added_mos_auto_grow, &
1152 1582 : separate_spin_occupations=restricted_roks)
1153 : END IF
1154 :
1155 2670 : IF (ALLOCATED(trial_eigenvalues)) DEALLOCATE (trial_eigenvalues)
1156 2670 : IF (ALLOCATED(trial_occupation)) DEALLOCATE (trial_occupation)
1157 2670 : IF (ALLOCATED(trial_mu)) DEALLOCATE (trial_mu)
1158 :
1159 2670 : IF (.NOT. added_mos_auto_grow) THEN
1160 2664 : fmwork => scf_env%scf_work1
1161 2664 : IF (ASSOCIATED(scf_env%p_mix_new)) THEN
1162 0 : density_target => scf_env%p_mix_new
1163 : ELSE
1164 2664 : CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp)
1165 2664 : density_target => rho_ao_kp
1166 : END IF
1167 2664 : CALL kpoint_density_matrices(kpoints)
1168 : CALL kpoint_density_transform(kpoints, density_target, .FALSE., &
1169 : matrix_s(1, 1)%matrix, sab_nl, fmwork, &
1170 2664 : overlap_rs=matrix_s)
1171 : END IF
1172 2670 : scf_env%iter_method = scf_env%qs_ot_env(first_channel)%OT_METHOD_FULL
1173 2670 : scf_env%iter_param = scf_env%qs_ot_env(first_channel)%ds_min
1174 2670 : qs_env%broyden_adaptive_sigma = scf_env%qs_ot_env(first_channel)%broyden_adaptive_sigma
1175 :
1176 6949 : DO local_channel = 1, nlocal_channels
1177 4279 : CALL dbcsr_release(matrix_hc(local_channel)%matrix)
1178 4279 : CALL dbcsr_release(matrix_hc_im(local_channel)%matrix)
1179 4279 : DEALLOCATE (matrix_hc(local_channel)%matrix)
1180 4279 : DEALLOCATE (matrix_hc_im(local_channel)%matrix)
1181 6949 : IF (local_ot_env(1)%settings%occupation_preconditioner) THEN
1182 1948 : CALL dbcsr_release(matrix_hc_physical(local_channel)%matrix)
1183 1948 : CALL dbcsr_release(matrix_hc_physical_im(local_channel)%matrix)
1184 1948 : DEALLOCATE (matrix_hc_physical(local_channel)%matrix)
1185 1948 : DEALLOCATE (matrix_hc_physical_im(local_channel)%matrix)
1186 : END IF
1187 : END DO
1188 2670 : DEALLOCATE (matrix_hc, matrix_hc_im)
1189 2670 : IF (local_ot_env(1)%settings%occupation_preconditioner) THEN
1190 1430 : DEALLOCATE (matrix_hc_physical, matrix_hc_physical_im)
1191 : END IF
1192 :
1193 2670 : CALL timestop(handle)
1194 :
1195 5340 : END SUBROUTINE qs_scf_loop_do_ot_kp
1196 :
1197 : ! **************************************************************************************************
1198 : !> \brief roll one accepted physical K-point Mermin endpoint into a bounded response history
1199 : !> \param history physical history for one spin density
1200 : !> \param c_re current real orbital coefficients
1201 : !> \param c_im current imaginary orbital coefficients
1202 : !> \param h_re current real self-consistent Hamiltonian
1203 : !> \param h_im current imaginary self-consistent Hamiltonian
1204 : !> \param occupation current fixed-N occupations
1205 : ! **************************************************************************************************
1206 998 : SUBROUTINE qs_scf_kp_roll_physical_secant(history, c_re, c_im, h_re, h_im, occupation)
1207 : TYPE(qs_ot_physical_secant_type), INTENT(INOUT) :: history
1208 : TYPE(cp_fm_type), INTENT(IN), TARGET :: c_re, c_im, h_re, h_im
1209 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: occupation
1210 :
1211 : INTEGER :: nmo
1212 : LOGICAL :: previous_valid
1213 :
1214 998 : nmo = SIZE(occupation)
1215 998 : CPASSERT(nmo > 0)
1216 : previous_valid = history%reference_valid .AND. &
1217 : ASSOCIATED(history%c0) .AND. ASSOCIATED(history%c0_im) .AND. &
1218 : ASSOCIATED(history%h0) .AND. ASSOCIATED(history%h0_im) .AND. &
1219 998 : ALLOCATED(history%occupation0)
1220 874 : IF (previous_valid) previous_valid = SIZE(history%occupation0) == nmo
1221 124 : history%secant_valid = .FALSE.
1222 874 : IF (previous_valid) THEN
1223 874 : CALL store_physical_fm_matrix(history%c0, history%c_previous, "mermin_previous_c")
1224 874 : CALL store_physical_fm_matrix(history%c0_im, history%c_previous_im, "mermin_previous_c_im")
1225 874 : CALL store_physical_fm_matrix(h_re, history%y_previous, "mermin_previous_y")
1226 874 : CALL store_physical_fm_matrix(h_im, history%y_previous_im, "mermin_previous_y_im")
1227 : CALL cp_fm_scale_and_add( &
1228 874 : 1.0_dp, history%y_previous, -1.0_dp, history%h0)
1229 : CALL cp_fm_scale_and_add( &
1230 874 : 1.0_dp, history%y_previous_im, -1.0_dp, history%h0_im)
1231 874 : IF (ALLOCATED(history%occupation_previous)) THEN
1232 772 : IF (SIZE(history%occupation_previous) /= nmo) THEN
1233 0 : DEALLOCATE (history%occupation_previous)
1234 : END IF
1235 : END IF
1236 874 : IF (.NOT. ALLOCATED(history%occupation_previous)) THEN
1237 306 : ALLOCATE (history%occupation_previous(nmo))
1238 : END IF
1239 6392 : history%occupation_previous(:) = history%occupation0(:)
1240 874 : history%secant_valid = .TRUE.
1241 : END IF
1242 :
1243 998 : CALL store_physical_fm_matrix(c_re, history%c0, "mermin_c0")
1244 998 : CALL store_physical_fm_matrix(c_im, history%c0_im, "mermin_c0_im")
1245 998 : CALL store_physical_fm_matrix(h_re, history%h0, "mermin_h0")
1246 998 : CALL store_physical_fm_matrix(h_im, history%h0_im, "mermin_h0_im")
1247 998 : IF (ALLOCATED(history%occupation0)) THEN
1248 874 : IF (SIZE(history%occupation0) /= nmo) DEALLOCATE (history%occupation0)
1249 : END IF
1250 1246 : IF (.NOT. ALLOCATED(history%occupation0)) ALLOCATE (history%occupation0(nmo))
1251 7294 : history%occupation0(:) = occupation(:)
1252 998 : history%reference_valid = .TRUE.
1253 :
1254 998 : END SUBROUTINE qs_scf_kp_roll_physical_secant
1255 :
1256 : ! **************************************************************************************************
1257 : !> \brief store the raw physical Mermin gradient at one accepted K-point endpoint
1258 : !> \param qs_ot_env local spin/K-point OT channel
1259 : ! **************************************************************************************************
1260 500 : SUBROUTINE qs_scf_kp_store_mermin_gradient(qs_ot_env)
1261 : TYPE(qs_ot_type) :: qs_ot_env
1262 :
1263 500 : CPASSERT(ASSOCIATED(qs_ot_env%matrix_mermin_g0))
1264 500 : CPASSERT(ASSOCIATED(qs_ot_env%matrix_mermin_g0_im))
1265 500 : CPASSERT(ASSOCIATED(qs_ot_env%rot_mat_mermin_g0))
1266 500 : CPASSERT(ASSOCIATED(qs_ot_env%rot_mat_mermin_g0_im))
1267 500 : CPASSERT(ALLOCATED(qs_ot_env%ener_mermin_g0))
1268 500 : CPASSERT(SIZE(qs_ot_env%ener_mermin_g0) == SIZE(qs_ot_env%ener_gx))
1269 :
1270 500 : CALL dbcsr_copy(qs_ot_env%matrix_mermin_g0, qs_ot_env%matrix_gx)
1271 500 : CALL dbcsr_copy(qs_ot_env%matrix_mermin_g0_im, qs_ot_env%matrix_gx_im)
1272 500 : CALL dbcsr_copy(qs_ot_env%rot_mat_mermin_g0, qs_ot_env%rot_mat_gx)
1273 500 : CALL dbcsr_copy(qs_ot_env%rot_mat_mermin_g0_im, qs_ot_env%rot_mat_gx_im)
1274 3452 : qs_ot_env%ener_mermin_g0(:) = qs_ot_env%ener_gx(:)
1275 500 : qs_ot_env%mermin_gradient_ref_valid = .TRUE.
1276 :
1277 500 : END SUBROUTINE qs_scf_kp_store_mermin_gradient
1278 :
1279 : ! **************************************************************************************************
1280 : !> \brief copy an owned physical full-matrix snapshot, reallocating on shape changes
1281 : !> \param source transient or persistent source matrix
1282 : !> \param snapshot owned response-history matrix
1283 : !> \param name matrix name
1284 : ! **************************************************************************************************
1285 7488 : SUBROUTINE store_physical_fm_matrix(source, snapshot, name)
1286 : TYPE(cp_fm_type), INTENT(IN), TARGET :: source
1287 : TYPE(cp_fm_type), INTENT(INOUT), POINTER :: snapshot
1288 : CHARACTER(LEN=*), INTENT(IN) :: name
1289 :
1290 : INTEGER :: ncol_snapshot, ncol_source, &
1291 : nrow_snapshot, nrow_source
1292 : LOGICAL :: shape_changed
1293 :
1294 7488 : CALL cp_fm_get_info(source, nrow_global=nrow_source, ncol_global=ncol_source)
1295 7488 : shape_changed = .FALSE.
1296 7488 : IF (ASSOCIATED(snapshot)) THEN
1297 6584 : CALL cp_fm_get_info(snapshot, nrow_global=nrow_snapshot, ncol_global=ncol_snapshot)
1298 6584 : shape_changed = nrow_snapshot /= nrow_source .OR. ncol_snapshot /= ncol_source
1299 : IF (shape_changed) THEN
1300 0 : CALL cp_fm_release(snapshot)
1301 0 : DEALLOCATE (snapshot)
1302 : END IF
1303 : END IF
1304 7488 : IF (.NOT. ASSOCIATED(snapshot)) THEN
1305 904 : ALLOCATE (snapshot)
1306 904 : CALL cp_fm_create(snapshot, source%matrix_struct, name=name)
1307 : END IF
1308 7488 : CALL cp_fm_to_fm(source, snapshot)
1309 :
1310 7488 : END SUBROUTINE store_physical_fm_matrix
1311 :
1312 : ! **************************************************************************************************
1313 : !> \brief project one accepted physical Hxc secant onto current-orbital density tangents
1314 : !> \param history accepted endpoint snapshots for one physical spin density
1315 : !> \param overlap_re real part of the AO overlap matrix
1316 : !> \param overlap_im imaginary part of the AO overlap matrix
1317 : !> \param density_modes current-orbital density tangents including the K-point weight
1318 : !> \param kpoint_weight irreducible K-point weight
1319 : !> \param density_norm_sq physical density-step norm
1320 : !> \param response_work density/Hamiltonian response work
1321 : !> \param density_overlap projections of the density step onto density_modes
1322 : !> \param response_overlap projections of the Hamiltonian step onto density_modes
1323 : !> \param valid whether a finite accepted secant was available
1324 : ! **************************************************************************************************
1325 256 : SUBROUTINE qs_scf_kp_physical_secant_overlaps( &
1326 256 : history, overlap_re, overlap_im, density_modes, kpoint_weight, density_norm_sq, &
1327 256 : response_work, density_overlap, response_overlap, valid)
1328 :
1329 : TYPE(qs_ot_physical_secant_type), INTENT(IN) :: history
1330 : TYPE(cp_fm_type), INTENT(IN) :: overlap_re, overlap_im
1331 : COMPLEX(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: density_modes
1332 : REAL(KIND=dp), INTENT(IN) :: kpoint_weight
1333 : REAL(KIND=dp), INTENT(OUT) :: density_norm_sq, response_work
1334 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: density_overlap, response_overlap
1335 : LOGICAL, INTENT(OUT) :: valid
1336 :
1337 256 : COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: hamiltonian_step_current, &
1338 256 : hamiltonian_step_previous, &
1339 256 : overlap_previous_current
1340 : INTEGER :: nmo
1341 256 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: dense_im, dense_re
1342 : TYPE(cp_fm_type) :: action_im, action_re, cross_im, cross_re, &
1343 : projection_im, projection_re
1344 :
1345 256 : nmo = SIZE(density_modes, 1)
1346 0 : CPASSERT(nmo > 0)
1347 256 : CPASSERT(SIZE(density_modes, 2) == nmo)
1348 256 : CPASSERT(SIZE(density_overlap) == SIZE(density_modes, 3))
1349 256 : CPASSERT(SIZE(response_overlap) == SIZE(density_modes, 3))
1350 : valid = history%secant_valid .AND. &
1351 : ASSOCIATED(history%c0) .AND. ASSOCIATED(history%c0_im) .AND. &
1352 : ASSOCIATED(history%c_previous) .AND. ASSOCIATED(history%c_previous_im) .AND. &
1353 : ASSOCIATED(history%y_previous) .AND. ASSOCIATED(history%y_previous_im) .AND. &
1354 256 : ALLOCATED(history%occupation0) .AND. ALLOCATED(history%occupation_previous)
1355 256 : density_norm_sq = 0.0_dp
1356 256 : response_work = 0.0_dp
1357 768 : density_overlap(:) = 0.0_dp
1358 768 : response_overlap(:) = 0.0_dp
1359 256 : IF (.NOT. valid) RETURN
1360 256 : CPASSERT(SIZE(history%occupation0) == nmo)
1361 256 : CPASSERT(SIZE(history%occupation_previous) == nmo)
1362 :
1363 : CALL cp_fm_create(action_re, history%c0%matrix_struct, &
1364 256 : name="mermin_hxc_action_re", set_zero=.TRUE.)
1365 : CALL cp_fm_create(action_im, history%c0%matrix_struct, &
1366 256 : name="mermin_hxc_action_im", set_zero=.TRUE.)
1367 : CALL cp_fm_create(cross_re, history%c0%matrix_struct, &
1368 256 : name="mermin_hxc_cross_re", nrow=nmo, ncol=nmo, set_zero=.TRUE.)
1369 : CALL cp_fm_create(cross_im, history%c0%matrix_struct, &
1370 256 : name="mermin_hxc_cross_im", nrow=nmo, ncol=nmo, set_zero=.TRUE.)
1371 : CALL cp_fm_create(projection_re, history%c0%matrix_struct, &
1372 256 : name="mermin_hxc_projection_re", nrow=nmo, ncol=nmo, set_zero=.TRUE.)
1373 : CALL cp_fm_create(projection_im, history%c0%matrix_struct, &
1374 256 : name="mermin_hxc_projection_im", nrow=nmo, ncol=nmo, set_zero=.TRUE.)
1375 : ALLOCATE (dense_re(nmo, nmo), dense_im(nmo, nmo), &
1376 : hamiltonian_step_current(nmo, nmo), hamiltonian_step_previous(nmo, nmo), &
1377 3328 : overlap_previous_current(nmo, nmo))
1378 :
1379 : CALL multiply_complex_fm_nn( &
1380 : overlap_re, overlap_im, history%c0, history%c0_im, &
1381 256 : action_re, action_im, 1.0_dp, 0.0_dp)
1382 : CALL multiply_complex_fm_adjoint( &
1383 : history%c_previous, history%c_previous_im, &
1384 256 : action_re, action_im, cross_re, cross_im, 1.0_dp, 0.0_dp)
1385 256 : CALL cp_fm_get_submatrix(cross_re, dense_re)
1386 256 : CALL cp_fm_get_submatrix(cross_im, dense_im)
1387 13680 : overlap_previous_current(:, :) = CMPLX(dense_re, dense_im, KIND=dp)
1388 :
1389 : CALL project_physical_response( &
1390 : history, history%c0, history%c0_im, &
1391 256 : action_re, action_im, projection_re, projection_im, hamiltonian_step_current)
1392 : CALL project_physical_response( &
1393 : history, history%c_previous, history%c_previous_im, &
1394 256 : action_re, action_im, projection_re, projection_im, hamiltonian_step_previous)
1395 : CALL qs_ot_density_secant_orbital_overlaps( &
1396 : overlap_previous_current, history%occupation_previous, history%occupation0, &
1397 : hamiltonian_step_previous, hamiltonian_step_current, &
1398 : density_modes, kpoint_weight, density_norm_sq, response_work, density_overlap, &
1399 256 : response_overlap, valid)
1400 :
1401 0 : DEALLOCATE (dense_re, dense_im, hamiltonian_step_current, &
1402 256 : hamiltonian_step_previous, overlap_previous_current)
1403 256 : CALL cp_fm_release(action_re)
1404 256 : CALL cp_fm_release(action_im)
1405 256 : CALL cp_fm_release(cross_re)
1406 256 : CALL cp_fm_release(cross_im)
1407 256 : CALL cp_fm_release(projection_re)
1408 512 : CALL cp_fm_release(projection_im)
1409 :
1410 : CONTAINS
1411 :
1412 : ! **************************************************************************************************
1413 : !> \brief ...
1414 : !> \param env ...
1415 : !> \param c_re ...
1416 : !> \param c_im ...
1417 : !> \param action_re ...
1418 : !> \param action_im ...
1419 : !> \param projected_re ...
1420 : !> \param projected_im ...
1421 : !> \param projected ...
1422 : ! **************************************************************************************************
1423 512 : SUBROUTINE project_physical_response( &
1424 512 : env, c_re, c_im, action_re, action_im, projected_re, projected_im, projected)
1425 : TYPE(qs_ot_physical_secant_type), INTENT(IN) :: env
1426 : TYPE(cp_fm_type), INTENT(IN) :: c_re, c_im, action_re, action_im, &
1427 : projected_re, projected_im
1428 : COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: projected
1429 :
1430 : CALL multiply_complex_fm_nn( &
1431 : env%y_previous, env%y_previous_im, c_re, c_im, &
1432 512 : action_re, action_im, 1.0_dp, 0.0_dp)
1433 : CALL multiply_complex_fm_adjoint( &
1434 512 : c_re, c_im, action_re, action_im, projected_re, projected_im, 1.0_dp, 0.0_dp)
1435 512 : CALL cp_fm_get_submatrix(projected_re, dense_re)
1436 512 : CALL cp_fm_get_submatrix(projected_im, dense_im)
1437 27360 : projected(:, :) = CMPLX(dense_re, dense_im, KIND=dp)
1438 54720 : projected(:, :) = 0.5_dp*(projected + CONJG(TRANSPOSE(projected)))
1439 512 : END SUBROUTINE project_physical_response
1440 :
1441 : END SUBROUTINE qs_scf_kp_physical_secant_overlaps
1442 :
1443 : ! **************************************************************************************************
1444 : !> \brief accumulate C=beta*C+alpha*A*B for split-complex distributed full matrices
1445 : !> \param a_re ...
1446 : !> \param a_im ...
1447 : !> \param b_re ...
1448 : !> \param b_im ...
1449 : !> \param c_re ...
1450 : !> \param c_im ...
1451 : !> \param alpha ...
1452 : !> \param beta ...
1453 : ! **************************************************************************************************
1454 3102 : SUBROUTINE multiply_complex_fm_nn(a_re, a_im, b_re, b_im, c_re, c_im, alpha, beta)
1455 : TYPE(cp_fm_type), INTENT(IN) :: a_re, a_im, b_re, b_im, c_re, c_im
1456 : REAL(KIND=dp), INTENT(IN) :: alpha, beta
1457 :
1458 : INTEGER :: k, k_b, m, m_c, n, n_c
1459 :
1460 1034 : CALL cp_fm_get_info(a_re, nrow_global=m, ncol_global=k)
1461 1034 : CALL cp_fm_get_info(b_re, nrow_global=k_b, ncol_global=n)
1462 1034 : CALL cp_fm_get_info(c_re, nrow_global=m_c, ncol_global=n_c)
1463 1034 : CPASSERT(k_b == k .AND. m_c == m .AND. n_c == n)
1464 1034 : CALL cp_complex_fm_gemm('N', 'N', m, n, k, alpha, a_re, a_im, b_re, b_im, beta, c_re, c_im)
1465 :
1466 1034 : END SUBROUTINE multiply_complex_fm_nn
1467 :
1468 : ! **************************************************************************************************
1469 : !> \brief accumulate C=beta*C+alpha*A^H*B for split-complex distributed full matrices
1470 : !> \param a_re ...
1471 : !> \param a_im ...
1472 : !> \param b_re ...
1473 : !> \param b_im ...
1474 : !> \param c_re ...
1475 : !> \param c_im ...
1476 : !> \param alpha ...
1477 : !> \param beta ...
1478 : ! **************************************************************************************************
1479 9702 : SUBROUTINE multiply_complex_fm_adjoint(a_re, a_im, b_re, b_im, c_re, c_im, alpha, beta)
1480 : TYPE(cp_fm_type), INTENT(IN) :: a_re, a_im, b_re, b_im, c_re, c_im
1481 : REAL(KIND=dp), INTENT(IN) :: alpha, beta
1482 :
1483 : INTEGER :: k, k_b, m, m_c, n, n_c
1484 :
1485 3234 : CALL cp_fm_get_info(a_re, nrow_global=k, ncol_global=m)
1486 3234 : CALL cp_fm_get_info(b_re, nrow_global=k_b, ncol_global=n)
1487 3234 : CALL cp_fm_get_info(c_re, nrow_global=m_c, ncol_global=n_c)
1488 3234 : CPASSERT(k_b == k .AND. m_c == m .AND. n_c == n)
1489 3234 : CALL cp_complex_fm_gemm('C', 'N', m, n, k, alpha, a_re, a_im, b_re, b_im, beta, c_re, c_im)
1490 :
1491 3234 : END SUBROUTINE multiply_complex_fm_adjoint
1492 :
1493 : ! **************************************************************************************************
1494 : !> \brief accumulate C=beta*C+alpha*A*B^H for split-complex distributed full matrices
1495 : !> \param a_re real part of A
1496 : !> \param a_im imaginary part of A
1497 : !> \param b_re real part of B
1498 : !> \param b_im imaginary part of B
1499 : !> \param c_re real part of C
1500 : !> \param c_im imaginary part of C
1501 : !> \param alpha product scaling factor
1502 : !> \param beta initial C scaling factor
1503 : ! **************************************************************************************************
1504 714 : SUBROUTINE multiply_complex_fm_right_adjoint(a_re, a_im, b_re, b_im, c_re, c_im, alpha, beta)
1505 : TYPE(cp_fm_type), INTENT(IN) :: a_re, a_im, b_re, b_im, c_re, c_im
1506 : REAL(KIND=dp), INTENT(IN) :: alpha, beta
1507 :
1508 : INTEGER :: k, k_b, m, m_c, n, n_c
1509 :
1510 238 : CALL cp_fm_get_info(a_re, nrow_global=m, ncol_global=k)
1511 238 : CALL cp_fm_get_info(b_re, nrow_global=n, ncol_global=k_b)
1512 238 : CALL cp_fm_get_info(c_re, nrow_global=m_c, ncol_global=n_c)
1513 238 : CPASSERT(k_b == k .AND. m_c == m .AND. n_c == n)
1514 238 : CALL cp_complex_fm_gemm('N', 'C', m, n, k, alpha, a_re, a_im, b_re, b_im, beta, c_re, c_im)
1515 :
1516 238 : END SUBROUTINE multiply_complex_fm_right_adjoint
1517 :
1518 : ! **************************************************************************************************
1519 : !> \brief evaluate the exact GPW Hartree-XC curvature of a complex K-point density mode
1520 : !> \param qs_env ground-state QS environment defining the current density and XC kernel
1521 : !> \param kpoints K-point environment and group distribution
1522 : !> \param local_ot_env local spin/K-point OT channels
1523 : !> \param first_channel global index of the first local channel
1524 : !> \param nspin_ot number of physical spin channels
1525 : !> \param wkp irreducible K-point weights
1526 : !> \param para_env_inter_kp communicator between distributed K-point groups
1527 : !> \param sab_nl real-space neighbor lists used by the K-point density transform
1528 : !> \param fmwork full-matrix workspaces on the global communicator
1529 : !> \param rho_ao_kp template AO density matrices for the real-space images
1530 : !> \param channel_density_modes weighted current-orbital density modes per local channel
1531 : !> \param number_of_modes number of populated modes to apply
1532 : !> \param projected_hxc projected Hartree-XC Hessian in the supplied density-mode basis
1533 : !> \param valid whether the projected block is finite and Hermitian within tolerance
1534 : ! **************************************************************************************************
1535 114 : SUBROUTINE qs_scf_kp_exact_hxc_projected( &
1536 : qs_env, kpoints, local_ot_env, first_channel, nspin_ot, wkp, para_env_inter_kp, &
1537 228 : sab_nl, fmwork, rho_ao_kp, channel_density_modes, number_of_modes, &
1538 114 : projected_hxc, valid)
1539 :
1540 : TYPE(qs_environment_type), POINTER :: qs_env
1541 : TYPE(kpoint_type), POINTER :: kpoints
1542 : TYPE(qs_ot_type), DIMENSION(:), POINTER :: local_ot_env
1543 : INTEGER, INTENT(IN) :: first_channel, nspin_ot
1544 : REAL(KIND=dp), DIMENSION(:), POINTER :: wkp
1545 : TYPE(mp_para_env_type), POINTER :: para_env_inter_kp
1546 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1547 : POINTER :: sab_nl
1548 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: fmwork
1549 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rho_ao_kp
1550 : COMPLEX(KIND=dp), DIMENSION(:, :, :, :), &
1551 : INTENT(IN) :: channel_density_modes
1552 : INTEGER, INTENT(IN) :: number_of_modes
1553 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: projected_hxc
1554 : LOGICAL, INTENT(OUT) :: valid
1555 :
1556 114 : COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: projected
1557 : INTEGER :: ikpoint, ispin, local_channel, local_kpoint, nao, nimages, nkp_groups, &
1558 : nlocal_kpoints, nmo, nmode, response_mode, test_mode
1559 : LOGICAL :: do_hfx
1560 : REAL(KIND=dp) :: symmetry_error, symmetry_scale
1561 114 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: dense_im, dense_re
1562 : TYPE(cp_fm_type) :: action_im, action_re, density_im, &
1563 : density_re, inactive_im, inactive_re, &
1564 : projected_im, projected_re
1565 114 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: operator_im, operator_re
1566 114 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :, :) :: density_kp
1567 114 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: density_rs, potential_rs
1568 : TYPE(dbcsr_type), POINTER :: matrix_k_im, matrix_k_re, template
1569 : TYPE(dft_control_type), POINTER :: dft_control
1570 : TYPE(kpoint_env_type), POINTER :: kp
1571 114 : TYPE(kpoint_operator_context_type) :: hxc_op_ctx
1572 : TYPE(qs_kpp1_env_type) :: kernel_env
1573 : TYPE(section_vals_type), POINTER :: hfx_section, input
1574 :
1575 798 : projected_hxc(:, :) = 0.0_dp
1576 114 : symmetry_error = 0.0_dp
1577 114 : valid = .FALSE.
1578 114 : NULLIFY (density_rs, dft_control, hfx_section, input, matrix_k_im, matrix_k_re, &
1579 114 : potential_rs, template)
1580 114 : CALL get_qs_env(qs_env, dft_control=dft_control, input=input)
1581 : IF (dft_control%qs_control%semi_empirical .OR. dft_control%qs_control%dftb .OR. &
1582 : dft_control%qs_control%xtb .OR. dft_control%qs_control%gapw .OR. &
1583 : dft_control%qs_control%gapw_xc .OR. dft_control%qs_control%lrigpw .OR. &
1584 114 : dft_control%qs_control%rigpw .OR. dft_control%do_admm) RETURN
1585 98 : hfx_section => section_vals_get_subs_vals(input, "DFT%XC%HF")
1586 98 : CALL section_vals_get(hfx_section, explicit=do_hfx)
1587 98 : IF (do_hfx) RETURN
1588 98 : CALL get_kpoint_info(kpoints, nkp_groups=nkp_groups)
1589 :
1590 98 : nlocal_kpoints = SIZE(kpoints%kp_env)
1591 98 : nmo = SIZE(channel_density_modes, 1)
1592 98 : nmode = number_of_modes
1593 98 : nimages = SIZE(rho_ao_kp, 2)
1594 98 : CPASSERT(nmo > 0)
1595 98 : CPASSERT(nmode > 0 .AND. nmode <= SIZE(channel_density_modes, 3))
1596 98 : CPASSERT(SIZE(channel_density_modes, 2) == nmo)
1597 98 : CPASSERT(SIZE(channel_density_modes, 4) == SIZE(local_ot_env))
1598 98 : IF (.NOT. (SIZE(projected_hxc, 1) >= nmode .AND. SIZE(projected_hxc, 2) >= nmode)) THEN
1599 0 : CALL cp_abort(__LOCATION__, "Projected HXC matrix is too small")
1600 : END IF
1601 98 : CPASSERT(SIZE(rho_ao_kp, 1) == nspin_ot)
1602 98 : CPASSERT(ASSOCIATED(fmwork))
1603 :
1604 98 : template => rho_ao_kp(1, 1)%matrix
1605 98 : CALL cp_fm_get_info(local_ot_env(1)%mermin_physical(1)%c0, nrow_global=nao)
1606 1064 : ALLOCATE (operator_re(nlocal_kpoints, nspin_ot), operator_im(nlocal_kpoints, nspin_ot))
1607 238 : DO local_kpoint = 1, nlocal_kpoints
1608 378 : DO ispin = 1, nspin_ot
1609 : local_channel = qs_ot_channel_index( &
1610 : ispin, kpoints%kp_env(local_kpoint)%kpoint_env%nkpoint, nspin_ot) - &
1611 140 : first_channel + 1
1612 : CALL cp_fm_create(operator_re(local_kpoint, ispin), &
1613 : local_ot_env(local_channel)%mermin_physical(1)%c0%matrix_struct, &
1614 140 : name="hxc_operator_k_re", nrow=nao, ncol=nao, set_zero=.TRUE.)
1615 : CALL cp_fm_create(operator_im(local_kpoint, ispin), &
1616 : local_ot_env(local_channel)%mermin_physical(1)%c0%matrix_struct, &
1617 280 : name="hxc_operator_k_im", nrow=nao, ncol=nao, set_zero=.TRUE.)
1618 : END DO
1619 : END DO
1620 882 : ALLOCATE (dense_re(nmo, nmo), dense_im(nmo, nmo), projected(nmo, nmo))
1621 : ! The group communicators do not own the complete real-space operator;
1622 : ! redistribute it to every group through the service. The fm pair lives
1623 : ! on the OT channel distribution; the service transfers into it
1624 : ! directly, whatever that distribution is. The context spans every
1625 : ! response mode: nothing inside the loop changes the host k-point
1626 : ! environment, and every mode rebuilds potential_rs from the same
1627 : ! rho_ao_kp template, so the work trio distribution stays valid.
1628 98 : IF (nkp_groups > 1) THEN
1629 32 : CALL kpoint_operator_context_create(hxc_op_ctx, kpoints, fmwork(1), nspin_ot)
1630 : END IF
1631 266 : DO response_mode = 1, nmode
1632 1652 : ALLOCATE (density_kp(nlocal_kpoints, 2, nspin_ot))
1633 406 : DO local_kpoint = 1, nlocal_kpoints
1634 238 : kp => kpoints%kp_env(local_kpoint)%kpoint_env
1635 238 : ikpoint = kp%nkpoint
1636 644 : DO ispin = 1, nspin_ot
1637 238 : local_channel = qs_ot_channel_index(ispin, ikpoint, nspin_ot) - first_channel + 1
1638 238 : IF (.NOT. (ASSOCIATED(local_ot_env(local_channel)%mermin_physical(1)%c0))) THEN
1639 0 : CALL cp_abort(__LOCATION__, "Missing physical OT coefficients")
1640 : END IF
1641 238 : IF (.NOT. (ASSOCIATED(local_ot_env(local_channel)%mermin_physical(1)%c0_im))) THEN
1642 0 : CALL cp_abort(__LOCATION__, "Missing imaginary physical OT coefficients")
1643 : END IF
1644 : CALL cp_fm_create( &
1645 : density_re, local_ot_env(local_channel)%mermin_physical(1)%c0%matrix_struct, &
1646 238 : name="hxc_density_mode_re", nrow=nmo, ncol=nmo)
1647 : CALL cp_fm_create( &
1648 : density_im, local_ot_env(local_channel)%mermin_physical(1)%c0%matrix_struct, &
1649 238 : name="hxc_density_mode_im", nrow=nmo, ncol=nmo)
1650 : CALL cp_fm_set_submatrix( &
1651 : density_re, REAL(channel_density_modes(:, :, response_mode, local_channel), &
1652 8954 : KIND=dp)/wkp(ikpoint))
1653 : CALL cp_fm_set_submatrix( &
1654 : density_im, AIMAG(channel_density_modes(:, :, response_mode, local_channel))/ &
1655 8954 : wkp(ikpoint))
1656 : CALL cp_fm_create( &
1657 : action_re, local_ot_env(local_channel)%mermin_physical(1)%c0%matrix_struct, &
1658 238 : name="hxc_density_action_re", set_zero=.TRUE.)
1659 : CALL cp_fm_create( &
1660 : action_im, local_ot_env(local_channel)%mermin_physical(1)%c0%matrix_struct, &
1661 238 : name="hxc_density_action_im", set_zero=.TRUE.)
1662 : CALL multiply_complex_fm_nn( &
1663 : local_ot_env(local_channel)%mermin_physical(1)%c0, &
1664 : local_ot_env(local_channel)%mermin_physical(1)%c0_im, density_re, density_im, &
1665 238 : action_re, action_im, 1.0_dp, 0.0_dp)
1666 : CALL cp_fm_get_info( &
1667 238 : local_ot_env(local_channel)%mermin_physical(1)%c0, nrow_global=nao)
1668 : CALL cp_fm_create(density_kp(local_kpoint, 1, ispin), &
1669 : local_ot_env(local_channel)%mermin_physical(1)%c0%matrix_struct, &
1670 238 : name="hxc_density_k_re", nrow=nao, ncol=nao, set_zero=.TRUE.)
1671 : CALL cp_fm_create(density_kp(local_kpoint, 2, ispin), &
1672 : local_ot_env(local_channel)%mermin_physical(1)%c0%matrix_struct, &
1673 238 : name="hxc_density_k_im", nrow=nao, ncol=nao, set_zero=.TRUE.)
1674 : CALL multiply_complex_fm_right_adjoint( &
1675 : action_re, action_im, local_ot_env(local_channel)%mermin_physical(1)%c0, &
1676 : local_ot_env(local_channel)%mermin_physical(1)%c0_im, &
1677 : density_kp(local_kpoint, 1, ispin), density_kp(local_kpoint, 2, ispin), &
1678 238 : 1.0_dp, 0.0_dp)
1679 238 : CALL cp_fm_release(action_re)
1680 238 : CALL cp_fm_release(action_im)
1681 238 : CALL cp_fm_release(density_re)
1682 952 : CALL cp_fm_release(density_im)
1683 : END DO
1684 : END DO
1685 :
1686 168 : CALL dbcsr_allocate_matrix_set(density_rs, nspin_ot, nimages)
1687 336 : DO ispin = 1, nspin_ot
1688 37578 : DO local_kpoint = 1, nimages
1689 37242 : ALLOCATE (density_rs(ispin, local_kpoint)%matrix)
1690 : CALL dbcsr_copy(density_rs(ispin, local_kpoint)%matrix, &
1691 37242 : rho_ao_kp(ispin, local_kpoint)%matrix, name="Hxc density response")
1692 37410 : CALL dbcsr_set(density_rs(ispin, local_kpoint)%matrix, 0.0_dp)
1693 : END DO
1694 : END DO
1695 : CALL kpoint_density_transform(kpoints, density_rs, .FALSE., template, sab_nl, fmwork, &
1696 168 : pmat_ext=density_kp)
1697 406 : DO local_kpoint = 1, nlocal_kpoints
1698 644 : DO ispin = 1, nspin_ot
1699 238 : CALL cp_fm_release(density_kp(local_kpoint, 1, ispin))
1700 476 : CALL cp_fm_release(density_kp(local_kpoint, 2, ispin))
1701 : END DO
1702 : END DO
1703 168 : DEALLOCATE (density_kp)
1704 :
1705 168 : CALL apply_hxc_kernel_kp(qs_env, kernel_env, density_rs, potential_rs)
1706 168 : IF (nkp_groups > 1) THEN
1707 : ! Refresh every response mode. The OT channel can use a different
1708 : ! group-local BLACS layout from the default AO/AO pool.
1709 54 : CALL kpoint_operator_prepare(hxc_op_ctx, potential_rs, kpoint_slot_ks, TARGET=operator_re(1, 1))
1710 216 : DO local_kpoint = 1, MAXVAL(kpoints%kp_dist(2, :) - kpoints%kp_dist(1, :) + 1)
1711 162 : DO ispin = 1, nspin_ot
1712 108 : IF (local_kpoint <= nlocal_kpoints) THEN
1713 : CALL kpoint_operator_get(hxc_op_ctx, local_kpoint, ispin, potential_rs, &
1714 : fm_re=operator_re(local_kpoint, ispin), &
1715 54 : fm_im=operator_im(local_kpoint, ispin), matrix_row=ispin, slot=kpoint_slot_ks)
1716 : ELSE
1717 : ! Participate in the global transfers without a local output in this round.
1718 : CALL kpoint_operator_get(hxc_op_ctx, local_kpoint, ispin, potential_rs, &
1719 0 : fm_re=inactive_re, fm_im=inactive_im, matrix_row=ispin, slot=kpoint_slot_ks)
1720 : END IF
1721 : END DO
1722 : END DO
1723 : END IF
1724 406 : DO local_kpoint = 1, nlocal_kpoints
1725 238 : kp => kpoints%kp_env(local_kpoint)%kpoint_env
1726 238 : ikpoint = kp%nkpoint
1727 644 : DO ispin = 1, nspin_ot
1728 238 : local_channel = qs_ot_channel_index(ispin, ikpoint, nspin_ot) - first_channel + 1
1729 : CALL cp_fm_create( &
1730 : action_re, local_ot_env(local_channel)%mermin_physical(1)%c0%matrix_struct, &
1731 238 : name="hxc_potential_action_re", set_zero=.TRUE.)
1732 : CALL cp_fm_create( &
1733 : action_im, local_ot_env(local_channel)%mermin_physical(1)%c0%matrix_struct, &
1734 238 : name="hxc_potential_action_im", set_zero=.TRUE.)
1735 : CALL kpoint_operator_get_local( &
1736 : potential_rs, kpoints, kp, ispin, operator_re(local_kpoint, ispin), &
1737 238 : operator_im(local_kpoint, ispin), matrix_k_re, matrix_k_im)
1738 : CALL cp_dbcsr_sm_fm_multiply( &
1739 238 : matrix_k_re, local_ot_env(local_channel)%mermin_physical(1)%c0, action_re, nmo)
1740 : CALL cp_dbcsr_sm_fm_multiply( &
1741 : matrix_k_im, local_ot_env(local_channel)%mermin_physical(1)%c0_im, action_re, nmo, &
1742 238 : alpha=-1.0_dp, beta=1.0_dp)
1743 : CALL cp_dbcsr_sm_fm_multiply( &
1744 238 : matrix_k_re, local_ot_env(local_channel)%mermin_physical(1)%c0_im, action_im, nmo)
1745 : CALL cp_dbcsr_sm_fm_multiply( &
1746 : matrix_k_im, local_ot_env(local_channel)%mermin_physical(1)%c0, action_im, nmo, &
1747 238 : alpha=1.0_dp, beta=1.0_dp)
1748 : CALL cp_fm_create( &
1749 : projected_re, local_ot_env(local_channel)%mermin_physical(1)%c0%matrix_struct, &
1750 238 : name="hxc_projected_re", nrow=nmo, ncol=nmo, set_zero=.TRUE.)
1751 : CALL cp_fm_create( &
1752 : projected_im, local_ot_env(local_channel)%mermin_physical(1)%c0%matrix_struct, &
1753 238 : name="hxc_projected_im", nrow=nmo, ncol=nmo, set_zero=.TRUE.)
1754 : CALL multiply_complex_fm_adjoint( &
1755 : local_ot_env(local_channel)%mermin_physical(1)%c0, &
1756 : local_ot_env(local_channel)%mermin_physical(1)%c0_im, action_re, action_im, &
1757 238 : projected_re, projected_im, 1.0_dp, 0.0_dp)
1758 238 : CALL cp_fm_get_submatrix(projected_re, dense_re)
1759 238 : CALL cp_fm_get_submatrix(projected_im, dense_im)
1760 8954 : projected(:, :) = CMPLX(dense_re, dense_im, KIND=dp)
1761 672 : DO test_mode = 1, nmode
1762 : projected_hxc(test_mode, response_mode) = &
1763 : projected_hxc(test_mode, response_mode) + REAL(SUM( &
1764 16596 : CONJG(channel_density_modes(:, :, test_mode, local_channel))*projected), KIND=dp)
1765 : END DO
1766 238 : CALL cp_fm_release(projected_re)
1767 238 : CALL cp_fm_release(projected_im)
1768 238 : CALL cp_fm_release(action_re)
1769 238 : CALL cp_fm_release(action_im)
1770 238 : CALL dbcsr_release_p(matrix_k_re)
1771 476 : CALL dbcsr_release_p(matrix_k_im)
1772 : END DO
1773 : END DO
1774 168 : CALL dbcsr_deallocate_matrix_set(density_rs)
1775 168 : CALL dbcsr_deallocate_matrix_set(potential_rs)
1776 266 : NULLIFY (density_rs, potential_rs)
1777 : END DO
1778 98 : IF (nkp_groups > 1) THEN
1779 32 : CALL kpoint_operator_context_release(hxc_op_ctx)
1780 : END IF
1781 :
1782 1274 : CALL para_env_inter_kp%sum(projected_hxc)
1783 : symmetry_error = MAXVAL(ABS( &
1784 : projected_hxc(1:nmode, 1:nmode) - &
1785 574 : TRANSPOSE(projected_hxc(1:nmode, 1:nmode))))
1786 : symmetry_scale = MAX(SQRT(EPSILON(1.0_dp)), &
1787 574 : MAXVAL(ABS(projected_hxc(1:nmode, 1:nmode))))
1788 : projected_hxc(1:nmode, 1:nmode) = 0.5_dp*( &
1789 : projected_hxc(1:nmode, 1:nmode) + &
1790 1050 : TRANSPOSE(projected_hxc(1:nmode, 1:nmode)))
1791 : valid = ALL(ieee_is_finite(projected_hxc(1:nmode, 1:nmode))) .AND. &
1792 : ieee_is_finite(symmetry_error) .AND. &
1793 574 : symmetry_error <= 1.0E-10_dp*symmetry_scale
1794 :
1795 98 : DEALLOCATE (dense_re, dense_im, projected)
1796 238 : DO local_kpoint = 1, nlocal_kpoints
1797 378 : DO ispin = 1, nspin_ot
1798 140 : CALL cp_fm_release(operator_re(local_kpoint, ispin))
1799 280 : CALL cp_fm_release(operator_im(local_kpoint, ispin))
1800 : END DO
1801 : END DO
1802 98 : DEALLOCATE (operator_re, operator_im)
1803 98 : CALL kpp1_release(kernel_env)
1804 :
1805 424 : END SUBROUTINE qs_scf_kp_exact_hxc_projected
1806 :
1807 : ! **************************************************************************************************
1808 : !> \brief project a split-complex Hamiltonian into a complex orbital basis
1809 : !> \param c_re real orbital coefficients
1810 : !> \param c_im imaginary orbital coefficients
1811 : !> \param h_re real Hamiltonian
1812 : !> \param h_im imaginary Hamiltonian
1813 : !> \param projected dense Hermitian projected Hamiltonian
1814 : ! **************************************************************************************************
1815 28 : SUBROUTINE qs_scf_kp_project_hamiltonian(c_re, c_im, h_re, h_im, projected)
1816 : TYPE(cp_fm_type), INTENT(IN) :: c_re, c_im, h_re, h_im
1817 : COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: projected
1818 :
1819 : INTEGER :: nmo
1820 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: dense_im, dense_re
1821 : TYPE(cp_fm_type) :: action_im, action_re, projection_im, &
1822 : projection_re
1823 :
1824 28 : nmo = SIZE(projected, 1)
1825 28 : CPASSERT(nmo > 0 .AND. SIZE(projected, 2) == nmo)
1826 28 : CALL cp_fm_create(action_re, c_re%matrix_struct, name="ot_response_hc_re", set_zero=.TRUE.)
1827 28 : CALL cp_fm_create(action_im, c_re%matrix_struct, name="ot_response_hc_im", set_zero=.TRUE.)
1828 : CALL cp_fm_create(projection_re, c_re%matrix_struct, name="ot_response_chc_re", &
1829 28 : nrow=nmo, ncol=nmo, set_zero=.TRUE.)
1830 : CALL cp_fm_create(projection_im, c_re%matrix_struct, name="ot_response_chc_im", &
1831 28 : nrow=nmo, ncol=nmo, set_zero=.TRUE.)
1832 168 : ALLOCATE (dense_re(nmo, nmo), dense_im(nmo, nmo))
1833 :
1834 : CALL multiply_complex_fm_nn( &
1835 28 : h_re, h_im, c_re, c_im, action_re, action_im, 1.0_dp, 0.0_dp)
1836 : CALL multiply_complex_fm_adjoint( &
1837 28 : c_re, c_im, action_re, action_im, projection_re, projection_im, 1.0_dp, 0.0_dp)
1838 28 : CALL cp_fm_get_submatrix(projection_re, dense_re)
1839 28 : CALL cp_fm_get_submatrix(projection_im, dense_im)
1840 3108 : projected(:, :) = CMPLX(dense_re, dense_im, KIND=dp)
1841 6216 : projected(:, :) = 0.5_dp*(projected + CONJG(TRANSPOSE(projected)))
1842 :
1843 28 : DEALLOCATE (dense_re, dense_im)
1844 28 : CALL cp_fm_release(action_re)
1845 28 : CALL cp_fm_release(action_im)
1846 28 : CALL cp_fm_release(projection_re)
1847 28 : CALL cp_fm_release(projection_im)
1848 :
1849 28 : END SUBROUTINE qs_scf_kp_project_hamiltonian
1850 :
1851 : ! **************************************************************************************************
1852 : !> \brief collect the global fixed-N trial smearing state used by k-point OT
1853 : !> \param kpoints ...
1854 : !> \param local_ot_env ...
1855 : !> \param first_channel ...
1856 : !> \param nspin_ot ...
1857 : !> \param nspin_energy number of physical spin-resolved auxiliary-energy blocks
1858 : !> \param restricted_roks whether both energy blocks share one orbital channel
1859 : !> \param kp_range ...
1860 : !> \param wkp ...
1861 : !> \param para_env_inter_kp ...
1862 : !> \param smear ...
1863 : !> \param trial_eigenvalues ...
1864 : !> \param trial_occupation ...
1865 : !> \param trial_mu chemical potential for each spin channel
1866 : !> \param trial_kTS ...
1867 : ! **************************************************************************************************
1868 1582 : SUBROUTINE qs_scf_loop_collect_ot_kp_trial_smearing( &
1869 : kpoints, local_ot_env, first_channel, nspin_ot, nspin_energy, restricted_roks, &
1870 : kp_range, wkp, para_env_inter_kp, smear, trial_eigenvalues, trial_occupation, trial_mu, trial_kTS)
1871 :
1872 : TYPE(kpoint_type), POINTER :: kpoints
1873 : TYPE(qs_ot_type), DIMENSION(:), POINTER :: local_ot_env
1874 : INTEGER, INTENT(IN) :: first_channel, nspin_ot, nspin_energy
1875 : LOGICAL, INTENT(IN) :: restricted_roks
1876 : INTEGER, DIMENSION(2), INTENT(IN) :: kp_range
1877 : REAL(KIND=dp), DIMENSION(:), POINTER :: wkp
1878 : TYPE(mp_para_env_type), POINTER :: para_env_inter_kp
1879 : TYPE(smear_type), INTENT(IN) :: smear
1880 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :), &
1881 : INTENT(OUT) :: trial_eigenvalues, trial_occupation
1882 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), &
1883 : INTENT(OUT) :: trial_mu
1884 : REAL(KIND=dp), INTENT(OUT) :: trial_kTS
1885 :
1886 : INTEGER :: energy_start, ikpoint, ispin, &
1887 : local_channel, local_kpoint, ne_a, &
1888 : ne_b, nelectron, nkpoint, nmo
1889 : REAL(KIND=dp) :: kTS, mu, nel, sigma
1890 : REAL(KIND=dp), DIMENSION(2) :: mus
1891 : TYPE(kpoint_env_type), POINTER :: kp
1892 :
1893 1582 : CPASSERT(smear%do_smear)
1894 3014 : SELECT CASE (smear%method)
1895 : CASE (smear_fermi_dirac)
1896 1432 : sigma = smear%electronic_temperature
1897 : CASE (smear_gaussian, smear_mp, smear_mv)
1898 150 : sigma = smear%smearing_width
1899 : CASE DEFAULT
1900 1582 : CPABORT("K-point Mermin OT does not support the selected smearing method")
1901 : END SELECT
1902 1582 : CPASSERT(sigma > EPSILON(sigma))
1903 1582 : CPASSERT(ASSOCIATED(para_env_inter_kp))
1904 1582 : CALL get_kpoint_info(kpoints, nkp=nkpoint)
1905 1582 : kp => kpoints%kp_env(1)%kpoint_env
1906 1582 : CALL get_mo_set(kp%mos(1, 1), nelectron=nelectron, nmo=nmo)
1907 1582 : ne_a = nelectron
1908 1582 : ne_b = 0
1909 1582 : IF (nspin_energy == 2) CALL get_mo_set(kp%mos(1, 2), nelectron=ne_b)
1910 :
1911 7910 : ALLOCATE (trial_eigenvalues(nmo, nkpoint, nspin_energy))
1912 6328 : ALLOCATE (trial_occupation(nmo, nkpoint, nspin_energy))
1913 4746 : ALLOCATE (trial_mu(nspin_energy))
1914 1582 : trial_eigenvalues(:, :, :) = 0.0_dp
1915 1582 : trial_occupation(:, :, :) = 0.0_dp
1916 1582 : mus(:) = 0.0_dp
1917 :
1918 3542 : DO local_kpoint = 1, SIZE(kpoints%kp_env)
1919 1960 : kp => kpoints%kp_env(local_kpoint)%kpoint_env
1920 1960 : ikpoint = kp%nkpoint
1921 1960 : CPASSERT(ikpoint >= kp_range(1) .AND. ikpoint <= kp_range(2))
1922 5692 : DO ispin = 1, nspin_energy
1923 : local_channel = qs_ot_channel_index(MERGE(1, ispin, restricted_roks), &
1924 4200 : ikpoint, nspin_ot) - first_channel + 1
1925 2150 : energy_start = 1
1926 2150 : IF (restricted_roks) energy_start = (ispin - 1)*nmo + 1
1927 : trial_eigenvalues(1:nmo, ikpoint, ispin) = &
1928 19034 : local_ot_env(local_channel)%ener_x(energy_start:energy_start + nmo - 1)
1929 : END DO
1930 : END DO
1931 1582 : CALL para_env_inter_kp%sum(trial_eigenvalues)
1932 :
1933 1582 : IF (nspin_energy == 1) THEN
1934 1392 : nel = REAL(nelectron, KIND=dp)
1935 : CALL Smearkp(trial_occupation(:, :, 1), mus(1), trial_kTS, &
1936 1392 : trial_eigenvalues(:, :, 1), nel, wkp, sigma, 2.0_dp, smear%method)
1937 190 : ELSE IF (restricted_roks .OR. smear%fixed_mag_mom > 0.0_dp) THEN
1938 94 : nel = REAL(ne_a, KIND=dp)
1939 : CALL Smearkp(trial_occupation(:, :, 1), mus(1), kTS, &
1940 94 : trial_eigenvalues(:, :, 1), nel, wkp, sigma, 1.0_dp, smear%method)
1941 94 : trial_kTS = kTS
1942 94 : nel = REAL(ne_b, KIND=dp)
1943 : CALL Smearkp(trial_occupation(:, :, 2), mus(2), kTS, &
1944 94 : trial_eigenvalues(:, :, 2), nel, wkp, sigma, 1.0_dp, smear%method)
1945 94 : trial_kTS = trial_kTS + kTS
1946 : ELSE
1947 96 : nel = REAL(ne_a + ne_b, KIND=dp)
1948 : CALL Smearkp2(trial_occupation, mu, trial_kTS, trial_eigenvalues, nel, wkp, &
1949 96 : sigma, smear%method)
1950 288 : mus(1:2) = mu
1951 : END IF
1952 3354 : trial_mu(:) = mus(1:nspin_energy)
1953 1582 : END SUBROUTINE qs_scf_loop_collect_ot_kp_trial_smearing
1954 :
1955 : ! **************************************************************************************************
1956 : !> \brief construct auxiliary-energy gradients on the global fixed-N Mermin surface
1957 : !> \param kpoints ...
1958 : !> \param local_ot_env ...
1959 : !> \param first_channel ...
1960 : !> \param nspin_ot ...
1961 : !> \param nspin_energy ...
1962 : !> \param restricted_roks ...
1963 : !> \param wkp ...
1964 : !> \param para_env_inter_kp ...
1965 : !> \param smear ...
1966 : !> \param trial_eigenvalues ...
1967 : !> \param trial_occupation ...
1968 : !> \param trial_mu ...
1969 : ! **************************************************************************************************
1970 1582 : SUBROUTINE qs_scf_loop_prepare_ot_kp_energy_gradient( &
1971 : kpoints, local_ot_env, first_channel, nspin_ot, nspin_energy, restricted_roks, &
1972 1582 : wkp, para_env_inter_kp, smear, trial_eigenvalues, trial_occupation, trial_mu)
1973 :
1974 : TYPE(kpoint_type), POINTER :: kpoints
1975 : TYPE(qs_ot_type), DIMENSION(:), POINTER :: local_ot_env
1976 : INTEGER, INTENT(IN) :: first_channel, nspin_ot, nspin_energy
1977 : LOGICAL, INTENT(IN) :: restricted_roks
1978 : REAL(KIND=dp), DIMENSION(:), POINTER :: wkp
1979 : TYPE(mp_para_env_type), POINTER :: para_env_inter_kp
1980 : TYPE(smear_type), INTENT(IN) :: smear
1981 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: trial_eigenvalues, trial_occupation
1982 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: trial_mu
1983 :
1984 : INTEGER :: energy_start, group, ikpoint, imo, &
1985 : ispin, local_channel, local_kpoint, &
1986 : ngroups, nmo
1987 : REAL(KIND=dp) :: fixed_n_mean, maxocc, sigma
1988 1582 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: fixed_n_residual, fixed_n_weight, &
1989 1582 : preconditioner_count, &
1990 1582 : preconditioner_shift, response_weight
1991 : TYPE(kpoint_env_type), POINTER :: kp
1992 :
1993 1582 : IF (.NOT. local_ot_env(1)%settings%do_ener) RETURN
1994 :
1995 1582 : CPASSERT(smear%do_smear)
1996 1582 : CPASSERT(ASSOCIATED(para_env_inter_kp))
1997 3014 : SELECT CASE (smear%method)
1998 : CASE (smear_fermi_dirac)
1999 1432 : sigma = smear%electronic_temperature
2000 : CASE (smear_gaussian, smear_mp, smear_mv)
2001 150 : sigma = smear%smearing_width
2002 : CASE DEFAULT
2003 1582 : CPABORT("K-point Mermin OT does not support the selected smearing method")
2004 : END SELECT
2005 1582 : CPASSERT(sigma > EPSILON(sigma))
2006 :
2007 1582 : ngroups = 1
2008 1582 : IF (restricted_roks .OR. &
2009 94 : (nspin_energy == 2 .AND. smear%fixed_mag_mom > 0.0_dp)) ngroups = 2
2010 1582 : nmo = SIZE(trial_eigenvalues, 1)
2011 : ALLOCATE (fixed_n_residual(ngroups), fixed_n_weight(ngroups), &
2012 : preconditioner_count(ngroups), preconditioner_shift(ngroups), &
2013 11074 : response_weight(nmo))
2014 1582 : fixed_n_residual(:) = 0.0_dp
2015 1582 : fixed_n_weight(:) = 0.0_dp
2016 1582 : preconditioner_count(:) = 0.0_dp
2017 1582 : preconditioner_shift(:) = 0.0_dp
2018 :
2019 3542 : DO local_kpoint = 1, SIZE(kpoints%kp_env)
2020 1960 : kp => kpoints%kp_env(local_kpoint)%kpoint_env
2021 1960 : ikpoint = kp%nkpoint
2022 5692 : DO ispin = 1, nspin_energy
2023 2150 : group = MERGE(ispin, 1, ngroups == nspin_energy)
2024 : local_channel = qs_ot_channel_index(MERGE(1, ispin, restricted_roks), &
2025 4200 : ikpoint, nspin_ot) - first_channel + 1
2026 2150 : energy_start = 1
2027 2150 : IF (restricted_roks) energy_start = (ispin - 1)*nmo + 1
2028 2150 : CALL get_mo_set(kp%mos(1, ispin), maxocc=maxocc, nmo=nmo)
2029 2150 : CPASSERT(maxocc > EPSILON(maxocc))
2030 2150 : IF (.NOT. (energy_start + nmo - 1 <= SIZE(local_ot_env(local_channel)%ener_x))) THEN
2031 0 : CALL cp_abort(__LOCATION__, "OT energy slice exceeds channel storage")
2032 : END IF
2033 : CALL smearing_response_weight( &
2034 : response_weight, trial_occupation(:, ikpoint, ispin), &
2035 : trial_eigenvalues(:, ikpoint, ispin), trial_mu(ispin), sigma, maxocc, nmo, &
2036 2150 : smear%method)
2037 17074 : response_weight(:) = wkp(ikpoint)*response_weight
2038 21184 : DO imo = 1, nmo
2039 : fixed_n_weight(group) = fixed_n_weight(group) + &
2040 14924 : response_weight(imo)
2041 : fixed_n_residual(group) = fixed_n_residual(group) + &
2042 : response_weight(imo)* &
2043 : (local_ot_env(local_channel)%ener_rayleigh(energy_start + imo - 1) - &
2044 17074 : trial_eigenvalues(imo, ikpoint, ispin))
2045 : END DO
2046 : END DO
2047 : END DO
2048 1582 : CALL para_env_inter_kp%sum(fixed_n_weight)
2049 1582 : CALL para_env_inter_kp%sum(fixed_n_residual)
2050 :
2051 1582 : IF (local_ot_env(1)%settings%occupation_preconditioner) THEN
2052 3238 : DO local_kpoint = 1, SIZE(kpoints%kp_env)
2053 1808 : kp => kpoints%kp_env(local_kpoint)%kpoint_env
2054 1808 : ikpoint = kp%nkpoint
2055 5236 : DO ispin = 1, nspin_energy
2056 1998 : group = MERGE(ispin, 1, ngroups == nspin_energy)
2057 : local_channel = qs_ot_channel_index(MERGE(1, ispin, restricted_roks), &
2058 3896 : ikpoint, nspin_ot) - first_channel + 1
2059 1998 : CALL get_mo_set(kp%mos(1, ispin), nmo=nmo)
2060 1998 : energy_start = 1
2061 1998 : IF (restricted_roks) energy_start = (ispin - 1)*nmo + 1
2062 1998 : fixed_n_mean = 0.0_dp
2063 1998 : IF (ABS(fixed_n_weight(group)) > EPSILON(fixed_n_weight(group))) THEN
2064 1954 : fixed_n_mean = fixed_n_residual(group)/fixed_n_weight(group)
2065 : END IF
2066 : local_ot_env(local_channel)%ener_preconditioned_gx( &
2067 : energy_start:energy_start + nmo - 1) = fixed_n_mean - &
2068 : (local_ot_env(local_channel)%ener_rayleigh( &
2069 : energy_start:energy_start + nmo - 1) - &
2070 13578 : trial_eigenvalues(:, ikpoint, ispin))
2071 : preconditioner_shift(group) = preconditioner_shift(group) + &
2072 : SUM(local_ot_env(local_channel)%ener_preconditioned_gx( &
2073 13578 : energy_start:energy_start + nmo - 1))
2074 5804 : preconditioner_count(group) = preconditioner_count(group) + REAL(nmo, KIND=dp)
2075 : END DO
2076 : END DO
2077 1430 : CALL para_env_inter_kp%sum(preconditioner_shift)
2078 1430 : CALL para_env_inter_kp%sum(preconditioner_count)
2079 : END IF
2080 :
2081 3542 : DO local_kpoint = 1, SIZE(kpoints%kp_env)
2082 1960 : kp => kpoints%kp_env(local_kpoint)%kpoint_env
2083 1960 : ikpoint = kp%nkpoint
2084 5692 : DO ispin = 1, nspin_energy
2085 2150 : group = MERGE(ispin, 1, ngroups == nspin_energy)
2086 : local_channel = qs_ot_channel_index(MERGE(1, ispin, restricted_roks), &
2087 4200 : ikpoint, nspin_ot) - first_channel + 1
2088 2150 : CALL get_mo_set(kp%mos(1, ispin), maxocc=maxocc, nmo=nmo)
2089 2150 : energy_start = 1
2090 2150 : IF (restricted_roks) energy_start = (ispin - 1)*nmo + 1
2091 2150 : IF (local_ot_env(local_channel)%settings%occupation_preconditioner .AND. &
2092 : preconditioner_count(group) > 0.0_dp) THEN
2093 : local_ot_env(local_channel)%ener_preconditioned_gx( &
2094 : energy_start:energy_start + nmo - 1) = &
2095 : local_ot_env(local_channel)%ener_preconditioned_gx( &
2096 : energy_start:energy_start + nmo - 1) - &
2097 13578 : preconditioner_shift(group)/preconditioner_count(group)
2098 : END IF
2099 : CALL smearing_response_weight( &
2100 : response_weight, trial_occupation(:, ikpoint, ispin), &
2101 : trial_eigenvalues(:, ikpoint, ispin), trial_mu(ispin), sigma, maxocc, nmo, &
2102 2150 : smear%method)
2103 17074 : response_weight(:) = wkp(ikpoint)*response_weight
2104 : CALL qs_ot_fixed_n_energy_gradient( &
2105 : local_ot_env(local_channel)%ener_rayleigh(energy_start:energy_start + nmo - 1), &
2106 : trial_eigenvalues(:, ikpoint, ispin), &
2107 : response_weight, fixed_n_weight(group), fixed_n_residual(group), &
2108 6260 : local_ot_env(local_channel)%ener_gx(energy_start:energy_start + nmo - 1))
2109 : END DO
2110 : END DO
2111 :
2112 0 : DEALLOCATE (fixed_n_residual, fixed_n_weight, preconditioner_count, preconditioner_shift, &
2113 1582 : response_weight)
2114 :
2115 1582 : END SUBROUTINE qs_scf_loop_prepare_ot_kp_energy_gradient
2116 :
2117 : ! **************************************************************************************************
2118 : !> \brief build a coupled finite-complex rotation/energy preconditioner for K-point Mermin OT
2119 : !> \param qs_env ground-state QS environment defining the current density and XC kernel
2120 : !> \param kpoints K-point environment
2121 : !> \param local_ot_env local spin/K-point OT channels
2122 : !> \param first_channel global index of the first local channel
2123 : !> \param nspin_ot number of spin channels
2124 : !> \param restricted_roks whether two physical spin densities share each orbital channel
2125 : !> \param wkp irreducible K-point weights
2126 : !> \param para_env_inter_kp communicator between distributed K-point groups
2127 : !> \param sab_nl real-space neighbor lists used by the K-point density transform
2128 : !> \param fmwork full-matrix workspaces on the global communicator
2129 : !> \param rho_ao_kp template AO density matrices for the real-space images
2130 : !> \param smear smearing settings
2131 : !> \param trial_eigenvalues current auxiliary band energies
2132 : !> \param trial_occupation current fixed-N occupations
2133 : !> \param trial_mu chemical potential for each spin channel
2134 : ! **************************************************************************************************
2135 828 : SUBROUTINE qs_scf_loop_prepare_ot_kp_rotation_response( &
2136 : qs_env, kpoints, local_ot_env, first_channel, nspin_ot, restricted_roks, &
2137 : wkp, para_env_inter_kp, sab_nl, fmwork, rho_ao_kp, smear, &
2138 828 : trial_eigenvalues, trial_occupation, trial_mu)
2139 : TYPE(qs_environment_type), POINTER :: qs_env
2140 : TYPE(kpoint_type), POINTER :: kpoints
2141 : TYPE(qs_ot_type), DIMENSION(:), POINTER :: local_ot_env
2142 : INTEGER, INTENT(IN) :: first_channel, nspin_ot
2143 : LOGICAL, INTENT(IN) :: restricted_roks
2144 : REAL(KIND=dp), DIMENSION(:), POINTER :: wkp
2145 : TYPE(mp_para_env_type), POINTER :: para_env_inter_kp
2146 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
2147 : POINTER :: sab_nl
2148 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: fmwork
2149 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rho_ao_kp
2150 : TYPE(smear_type), INTENT(IN) :: smear
2151 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: trial_eigenvalues, trial_occupation
2152 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: trial_mu
2153 :
2154 : INTEGER, PARAMETER :: finite_response_max_rotation = 800
2155 : REAL(KIND=dp), PARAMETER :: rotation_trust_rms = 0.25_dp, rotation_trust_spectral = 0.25_dp
2156 :
2157 828 : COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: chc, generator
2158 828 : COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: channel_generator, density_modes
2159 : COMPLEX(KIND=dp), ALLOCATABLE, &
2160 828 : DIMENSION(:, :, :, :) :: channel_density_modes
2161 : INTEGER :: energy_spin, energy_start, energy_stop, failure, gradient_secant_failure, group, &
2162 : i, ikpoint, ispin, j, local_channel, local_kpoint, nenergy, ngroups, nlocal_channels, &
2163 : nmo, nrotation, nrotation_groups, nspin_energy, physical_spin, r, rotation_group
2164 828 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: channel_response_group
2165 : LOGICAL :: accepted_cross_valid, accepted_curvature_valid, exact_hxc_valid, &
2166 : hxc_channel_valid, hxc_valid, projected_update_valid, shadow_pending, solve_valid
2167 : REAL(KIND=dp) :: accepted_cross, accepted_curvature, denominator, hxc_density_norm_sq, &
2168 : hxc_density_norm_sq_channel, hxc_response_work, hxc_response_work_channel, maxocc, &
2169 : response_cross, response_curvature_total, response_scale, shadow_curvature_total, sigma, &
2170 : tmp, unscaled_rotation_norm
2171 828 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: curvature_sum, energy_count, energy_mean, &
2172 828 : energy_metric_step, energy_shift, energy_weighted_sum, rayleigh_step, response_weight, &
2173 828 : rotation_count, rotation_norm, rotation_scale, schur_rhs, shadow_energy_mean, &
2174 828 : shadow_energy_metric_step, shadow_energy_weighted_sum, shadow_rayleigh_step, &
2175 828 : spin_rotation_gradient, vy
2176 828 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: channel_response_weight, &
2177 828 : channel_rotation_gradient, chc_im, chc_re, coupling, dense_im, dense_re, fixed_n_matrix, &
2178 828 : fixed_n_rhs, fixed_n_solution, generator_im, generator_re, previous_rotation_step, &
2179 828 : rayleigh_response, rhs, rotation_gradient, rotation_hessian, rotation_step, rotation_y, &
2180 828 : schur_block, solution, spin_rayleigh_response, spin_rotation_hessian, vz
2181 828 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: channel_rayleigh_response, &
2182 828 : channel_rotation_hessian, rotation_z
2183 : REAL(KIND=dp), DIMENSION(2) :: hxc_density_overlap, hxc_density_overlap_channel, &
2184 : hxc_response_overlap, hxc_response_overlap_channel, projected_coefficients, &
2185 : projected_gradient
2186 : REAL(KIND=dp), DIMENSION(2, 2) :: exact_hxc_correction, hxc_correction, &
2187 : projected_hessian
2188 : TYPE(cp_fm_type), POINTER :: mo_coeff, mo_coeff_im
2189 : TYPE(kpoint_env_type), POINTER :: kp
2190 :
2191 828 : nlocal_channels = SIZE(local_ot_env)
2192 1882 : DO local_channel = 1, nlocal_channels
2193 1054 : local_ot_env(local_channel)%rotation_response_valid = .FALSE.
2194 1054 : IF (ASSOCIATED(local_ot_env(local_channel)%rot_mat_response_gx)) THEN
2195 976 : CALL dbcsr_set(local_ot_env(local_channel)%rot_mat_response_gx, 0.0_dp)
2196 : END IF
2197 1054 : IF (ASSOCIATED(local_ot_env(local_channel)%rot_mat_response_gx_im)) THEN
2198 976 : CALL dbcsr_set(local_ot_env(local_channel)%rot_mat_response_gx_im, 0.0_dp)
2199 : END IF
2200 1054 : IF (ASSOCIATED(local_ot_env(local_channel)%ener_response_gx)) THEN
2201 7272 : local_ot_env(local_channel)%ener_response_gx(:) = 0.0_dp
2202 : END IF
2203 1054 : IF (ASSOCIATED(local_ot_env(local_channel)%matrix_response_gx)) THEN
2204 500 : CALL dbcsr_set(local_ot_env(local_channel)%matrix_response_gx, 0.0_dp)
2205 : END IF
2206 1882 : IF (ASSOCIATED(local_ot_env(local_channel)%matrix_response_gx_im)) THEN
2207 500 : CALL dbcsr_set(local_ot_env(local_channel)%matrix_response_gx_im, 0.0_dp)
2208 : END IF
2209 : END DO
2210 828 : local_ot_env(1)%response_model_curvature = 0.0_dp
2211 828 : local_ot_env(1)%response_shadow_curvature = 0.0_dp
2212 828 : local_ot_env(1)%response_hxc_direction_valid = .FALSE.
2213 : IF (.NOT. local_ot_env(1)%settings%occupation_preconditioner .OR. &
2214 828 : .NOT. local_ot_env(1)%settings%do_rotation .OR. &
2215 : .NOT. local_ot_env(1)%settings%do_ener) RETURN
2216 : ! The dense candidate is accepted through a physical line search. DIIS and Broyden use the
2217 : ! cheaper fixed-N response in their history residuals, but have no comparable accepted-step
2218 : ! model and therefore must not pay for an unused finite-response construction.
2219 750 : IF (local_ot_env(1)%settings%ot_method /= "CG" .AND. &
2220 : local_ot_env(1)%settings%ot_method /= "LBFG") RETURN
2221 : ! ROKS has one orbital generator but two physical spin-density and fixed-N energy blocks.
2222 : ! Keep those blocks separate until their response metrics have been accumulated.
2223 356 : IF (.NOT. ot_mermin_response_preparation_needed( &
2224 : local_ot_env(1)%delta, local_ot_env(1)%response_candidate_directions, &
2225 : local_ot_env(1)%response_shadow_good_samples, &
2226 : local_ot_env(1)%response_candidate_good_samples, &
2227 : local_ot_env(1)%response_candidate_cooldown, &
2228 : local_ot_env(1)%response_shadow_pending)) RETURN
2229 214 : CPASSERT(smear%do_smear)
2230 214 : CPASSERT(ASSOCIATED(para_env_inter_kp))
2231 392 : SELECT CASE (smear%method)
2232 : CASE (smear_fermi_dirac)
2233 178 : sigma = smear%electronic_temperature
2234 : CASE (smear_gaussian, smear_mp, smear_mv)
2235 36 : sigma = smear%smearing_width
2236 : CASE DEFAULT
2237 214 : RETURN
2238 : END SELECT
2239 214 : CPASSERT(sigma > EPSILON(sigma))
2240 :
2241 214 : nmo = SIZE(trial_eigenvalues, 1)
2242 214 : nspin_energy = SIZE(trial_eigenvalues, 3)
2243 214 : nenergy = SIZE(local_ot_env(1)%ener_x)
2244 214 : CPASSERT(SIZE(trial_occupation, 1) == nmo)
2245 214 : CPASSERT(SIZE(trial_occupation, 3) == nspin_energy)
2246 214 : CPASSERT(nenergy == nmo*nspin_energy/nspin_ot)
2247 214 : nrotation = nmo*(nmo - 1)
2248 214 : IF (nrotation <= 0 .OR. nrotation > finite_response_max_rotation) RETURN
2249 532 : DO local_channel = 1, nlocal_channels
2250 318 : CPASSERT(local_ot_env(local_channel)%has_complex_kpoint_state)
2251 318 : CPASSERT(SIZE(local_ot_env(local_channel)%ener_x) == nenergy)
2252 318 : CPASSERT(ASSOCIATED(local_ot_env(local_channel)%rot_mat_response_gx))
2253 318 : CPASSERT(ASSOCIATED(local_ot_env(local_channel)%rot_mat_response_gx_im))
2254 532 : CPASSERT(ASSOCIATED(local_ot_env(local_channel)%ener_response_gx))
2255 : END DO
2256 :
2257 214 : ngroups = 1
2258 214 : IF (restricted_roks) THEN
2259 : ngroups = nspin_energy
2260 200 : ELSE IF (nspin_ot == 2 .AND. smear%fixed_mag_mom > 0.0_dp) THEN
2261 14 : ngroups = 2
2262 : END IF
2263 214 : nrotation_groups = ngroups
2264 214 : IF (restricted_roks) nrotation_groups = 1
2265 214 : shadow_pending = local_ot_env(1)%response_shadow_pending
2266 : ALLOCATE (coupling(nrotation, ngroups), curvature_sum(ngroups), energy_count(ngroups), &
2267 : energy_mean(ngroups), energy_metric_step(nenergy), energy_shift(ngroups), &
2268 : energy_weighted_sum(ngroups), rayleigh_step(nenergy), response_weight(nmo), &
2269 : rotation_count(nrotation_groups), rotation_norm(nrotation_groups), &
2270 : rotation_scale(nrotation_groups), schur_rhs(nrotation), &
2271 : shadow_energy_mean(ngroups), shadow_energy_metric_step(nenergy), &
2272 : shadow_energy_weighted_sum(ngroups), shadow_rayleigh_step(nenergy), &
2273 : vy(ngroups), vz(ngroups, ngroups), rayleigh_response(nenergy, nrotation), &
2274 : rhs(nrotation, ngroups + 1), rotation_gradient(nrotation, 1), &
2275 : rotation_hessian(nrotation, nrotation), &
2276 : spin_rotation_gradient(nrotation), &
2277 : spin_rayleigh_response(nmo, nrotation), &
2278 : spin_rotation_hessian(nrotation, nrotation), &
2279 : previous_rotation_step(nrotation, nlocal_channels), &
2280 : rotation_step(nrotation, nlocal_channels), &
2281 : rotation_y(nrotation, nlocal_channels), rotation_z(nrotation, nlocal_channels, ngroups), &
2282 : schur_block(nrotation, nrotation), solution(nrotation, ngroups + 1), &
2283 : fixed_n_matrix(ngroups, ngroups), fixed_n_rhs(ngroups, 1), &
2284 : fixed_n_solution(ngroups, 1), &
2285 : channel_rayleigh_response(nenergy, nrotation, nlocal_channels), &
2286 : channel_rotation_hessian(nrotation, nrotation, nlocal_channels), &
2287 : channel_rotation_gradient(nrotation, nlocal_channels), &
2288 : channel_response_weight(nenergy, nlocal_channels), &
2289 : channel_response_group(nenergy, nlocal_channels), &
2290 18832 : channel_generator(nmo, nmo, nlocal_channels), density_modes(nmo, nmo, 2))
2291 1070 : ALLOCATE (channel_density_modes(nmo, nmo, 2, nlocal_channels))
2292 214 : curvature_sum(:) = 0.0_dp
2293 214 : vy(:) = 0.0_dp
2294 214 : vz(:, :) = 0.0_dp
2295 214 : rotation_y(:, :) = 0.0_dp
2296 214 : rotation_z(:, :, :) = 0.0_dp
2297 214 : previous_rotation_step(:, :) = 0.0_dp
2298 214 : channel_rayleigh_response(:, :, :) = 0.0_dp
2299 214 : channel_rotation_hessian(:, :, :) = 0.0_dp
2300 214 : channel_rotation_gradient(:, :) = 0.0_dp
2301 214 : channel_response_weight(:, :) = 0.0_dp
2302 214 : channel_response_group(:, :) = 0
2303 214 : channel_generator(:, :, :) = CMPLX(0.0_dp, 0.0_dp, KIND=dp)
2304 214 : channel_density_modes(:, :, :, :) = CMPLX(0.0_dp, 0.0_dp, KIND=dp)
2305 214 : failure = 0
2306 214 : response_cross = 0.0_dp
2307 :
2308 512 : DO local_kpoint = 1, SIZE(kpoints%kp_env)
2309 298 : kp => kpoints%kp_env(local_kpoint)%kpoint_env
2310 298 : ikpoint = kp%nkpoint
2311 830 : DO ispin = 1, nspin_ot
2312 318 : group = MERGE(ispin, 1, ngroups == nspin_ot)
2313 318 : local_channel = qs_ot_channel_index(ispin, ikpoint, nspin_ot) - first_channel + 1
2314 : CALL qs_scf_loop_ot_matrix_to_dense( &
2315 318 : local_ot_env(local_channel)%rot_mat_x, local_ot_env(local_channel), generator_re)
2316 : CALL qs_scf_loop_ot_matrix_to_dense( &
2317 318 : local_ot_env(local_channel)%rot_mat_x_im, local_ot_env(local_channel), generator_im)
2318 1272 : ALLOCATE (generator(nmo, nmo))
2319 13266 : generator(:, :) = CMPLX(generator_re, generator_im, KIND=dp)
2320 13266 : channel_generator(:, :, local_channel) = generator
2321 9762 : rotation_gradient(:, 1) = 0.0_dp
2322 318 : rotation_hessian(:, :) = 0.0_dp
2323 318 : rayleigh_response(:, :) = 0.0_dp
2324 :
2325 318 : IF (restricted_roks) THEN
2326 14 : IF (.NOT. (ASSOCIATED(kp%ot_hmat) .AND. SIZE(kp%ot_hmat, 2) >= nspin_energy)) THEN
2327 0 : CALL cp_abort(__LOCATION__, "Missing OT Hamiltonian spin cache")
2328 : END IF
2329 14 : CALL get_mo_set(kp%mos(1, 1), mo_coeff=mo_coeff)
2330 14 : CALL get_mo_set(kp%mos(2, 1), mo_coeff=mo_coeff_im)
2331 42 : DO energy_spin = 1, nspin_energy
2332 28 : energy_start = (energy_spin - 1)*nmo + 1
2333 28 : energy_stop = energy_start + nmo - 1
2334 28 : group = energy_spin
2335 28 : CALL get_mo_set(kp%mos(1, energy_spin), maxocc=maxocc)
2336 28 : CPASSERT(maxocc > EPSILON(maxocc))
2337 : CALL smearing_response_weight( &
2338 : response_weight, trial_occupation(:, ikpoint, energy_spin), &
2339 : trial_eigenvalues(:, ikpoint, energy_spin), trial_mu(energy_spin), &
2340 28 : sigma, maxocc, nmo, smear%method)
2341 308 : response_weight(:) = wkp(ikpoint)*response_weight
2342 308 : curvature_sum(group) = curvature_sum(group) + SUM(response_weight)
2343 308 : channel_response_weight(energy_start:energy_stop, local_channel) = response_weight
2344 308 : channel_response_group(energy_start:energy_stop, local_channel) = group
2345 :
2346 84 : ALLOCATE (chc(nmo, nmo))
2347 : CALL qs_scf_kp_project_hamiltonian( &
2348 : mo_coeff, mo_coeff_im, kp%ot_hmat(1, energy_spin), &
2349 28 : kp%ot_hmat(2, energy_spin), chc)
2350 : CALL qs_ot_finite_rotation_response( &
2351 : chc, generator, trial_occupation(:, ikpoint, energy_spin), wkp(ikpoint), &
2352 28 : spin_rotation_gradient, spin_rotation_hessian, spin_rayleigh_response)
2353 2548 : rotation_gradient(:, 1) = rotation_gradient(:, 1) + spin_rotation_gradient
2354 229348 : rotation_hessian(:, :) = rotation_hessian + spin_rotation_hessian
2355 27748 : rayleigh_response(energy_start:energy_stop, :) = spin_rayleigh_response
2356 70 : DEALLOCATE (chc)
2357 : END DO
2358 : ELSE
2359 304 : CALL get_mo_set(kp%mos(1, ispin), maxocc=maxocc)
2360 304 : CPASSERT(maxocc > EPSILON(maxocc))
2361 : CALL smearing_response_weight( &
2362 : response_weight, trial_occupation(:, ikpoint, ispin), &
2363 : trial_eigenvalues(:, ikpoint, ispin), trial_mu(ispin), sigma, maxocc, nmo, &
2364 304 : smear%method)
2365 1916 : response_weight(:) = wkp(ikpoint)*response_weight
2366 1916 : curvature_sum(group) = curvature_sum(group) + SUM(response_weight)
2367 1916 : channel_response_weight(:, local_channel) = response_weight
2368 1916 : channel_response_group(:, local_channel) = group
2369 : CALL qs_scf_loop_ot_matrix_to_dense( &
2370 304 : local_ot_env(local_channel)%rot_mat_chc, local_ot_env(local_channel), chc_re)
2371 : CALL qs_scf_loop_ot_matrix_to_dense( &
2372 304 : local_ot_env(local_channel)%rot_mat_chc_im, local_ot_env(local_channel), chc_im)
2373 912 : ALLOCATE (chc(nmo, nmo))
2374 11712 : chc(:, :) = CMPLX(chc_re, chc_im, KIND=dp)
2375 : CALL qs_ot_finite_rotation_response( &
2376 : chc, generator, trial_occupation(:, ikpoint, ispin), wkp(ikpoint), &
2377 304 : rotation_gradient(:, 1), rotation_hessian, rayleigh_response)
2378 608 : DEALLOCATE (chc, chc_im, chc_re)
2379 : END IF
2380 9762 : channel_rotation_gradient(:, local_channel) = rotation_gradient(:, 1)
2381 502794 : channel_rotation_hessian(:, :, local_channel) = rotation_hessian
2382 : CALL qs_ot_fixed_n_multigroup_schur_block( &
2383 : rotation_hessian, rayleigh_response, channel_response_weight(:, local_channel), &
2384 : channel_response_group(:, local_channel), rotation_gradient(:, 1), &
2385 318 : local_ot_env(local_channel)%ener_gx, schur_block, coupling, schur_rhs)
2386 9762 : rhs(:, 1) = schur_rhs
2387 678 : DO group = 1, ngroups
2388 11718 : rhs(:, group + 1) = coupling(:, group)
2389 : END DO
2390 : CALL qs_ot_symmetric_abs_solve( &
2391 318 : schur_block, rhs, solution, solve_valid, relative_floor=1.0E-4_dp)
2392 318 : IF (solve_valid) THEN
2393 9618 : rotation_y(:, local_channel) = solution(:, 1)
2394 630 : DO group = 1, ngroups
2395 11550 : rotation_z(:, local_channel, group) = solution(:, group + 1)
2396 : END DO
2397 91662 : channel_rayleigh_response(:, :, local_channel) = rayleigh_response
2398 630 : DO group = 1, ngroups
2399 11256 : vy(group) = vy(group) + DOT_PRODUCT(coupling(:, group), solution(:, 1))
2400 1050 : DO j = 1, ngroups
2401 : vz(group, j) = vz(group, j) + &
2402 14868 : DOT_PRODUCT(coupling(:, group), solution(:, j + 1))
2403 : END DO
2404 : END DO
2405 : ELSE
2406 24 : failure = 1
2407 : END IF
2408 318 : IF (shadow_pending .AND. local_ot_env(local_channel)%use_dx) THEN
2409 : CALL qs_scf_loop_ot_matrix_to_dense( &
2410 176 : local_ot_env(local_channel)%rot_mat_dx, local_ot_env(local_channel), dense_re)
2411 : CALL qs_scf_loop_ot_matrix_to_dense( &
2412 176 : local_ot_env(local_channel)%rot_mat_dx_im, local_ot_env(local_channel), dense_im)
2413 176 : r = 0
2414 1116 : DO i = 1, nmo - 1
2415 4452 : DO j = i + 1, nmo
2416 3336 : r = r + 1
2417 3336 : previous_rotation_step(r, local_channel) = dense_re(i, j)
2418 3336 : r = r + 1
2419 4276 : previous_rotation_step(r, local_channel) = dense_im(i, j)
2420 : END DO
2421 : END DO
2422 176 : CPASSERT(r == nrotation)
2423 176 : DEALLOCATE (dense_im, dense_re)
2424 : END IF
2425 616 : DEALLOCATE (generator, generator_im, generator_re)
2426 : END DO
2427 : END DO
2428 :
2429 214 : CALL para_env_inter_kp%sum(curvature_sum)
2430 214 : CALL para_env_inter_kp%sum(vy)
2431 214 : CALL para_env_inter_kp%sum(vz)
2432 214 : CALL para_env_inter_kp%sum(failure)
2433 214 : IF (failure /= 0) THEN
2434 0 : DEALLOCATE (coupling, curvature_sum, energy_count, energy_mean, energy_metric_step, &
2435 0 : energy_shift, energy_weighted_sum, rayleigh_step, response_weight, &
2436 0 : shadow_energy_mean, &
2437 0 : shadow_energy_metric_step, shadow_energy_weighted_sum, shadow_rayleigh_step, &
2438 0 : rotation_count, rotation_norm, rotation_scale, schur_rhs, vy, vz, &
2439 0 : rayleigh_response, rhs, rotation_gradient, rotation_hessian, &
2440 0 : spin_rotation_gradient, spin_rotation_hessian, spin_rayleigh_response, &
2441 0 : previous_rotation_step, rotation_step, &
2442 0 : rotation_y, rotation_z, schur_block, solution, channel_rayleigh_response, &
2443 0 : channel_rotation_hessian, channel_rotation_gradient, channel_response_weight, &
2444 0 : channel_density_modes, channel_response_group, channel_generator, &
2445 0 : fixed_n_matrix, fixed_n_rhs, &
2446 0 : fixed_n_solution, &
2447 30 : density_modes)
2448 30 : RETURN
2449 : END IF
2450 :
2451 664 : fixed_n_matrix(:, :) = 0.5_dp*(vz + TRANSPOSE(vz))
2452 396 : DO group = 1, ngroups
2453 396 : fixed_n_matrix(group, group) = fixed_n_matrix(group, group) + curvature_sum(group)
2454 : END DO
2455 396 : fixed_n_rhs(:, 1) = vy
2456 396 : fixed_n_solution(:, 1) = 0.0_dp
2457 184 : solve_valid = .FALSE.
2458 184 : IF (ngroups == 1) THEN
2459 156 : denominator = fixed_n_matrix(1, 1)
2460 156 : solve_valid = ABS(denominator) > EPSILON(denominator)
2461 156 : IF (solve_valid) fixed_n_solution(1, 1) = fixed_n_rhs(1, 1)/denominator
2462 28 : ELSE IF (ngroups == 2) THEN
2463 : denominator = fixed_n_matrix(1, 1)*fixed_n_matrix(2, 2) - &
2464 28 : fixed_n_matrix(1, 2)*fixed_n_matrix(2, 1)
2465 : solve_valid = ABS(denominator) > EPSILON(denominator)* &
2466 196 : MAX(1.0_dp, MAXVAL(ABS(fixed_n_matrix))**2)
2467 28 : IF (solve_valid) THEN
2468 : fixed_n_solution(1, 1) = (fixed_n_matrix(2, 2)*fixed_n_rhs(1, 1) - &
2469 14 : fixed_n_matrix(1, 2)*fixed_n_rhs(2, 1))/denominator
2470 : fixed_n_solution(2, 1) = (fixed_n_matrix(1, 1)*fixed_n_rhs(2, 1) - &
2471 14 : fixed_n_matrix(2, 1)*fixed_n_rhs(1, 1))/denominator
2472 : END IF
2473 : END IF
2474 580 : IF (.NOT. solve_valid .OR. .NOT. ALL(fixed_n_solution == fixed_n_solution)) THEN
2475 0 : DEALLOCATE (coupling, curvature_sum, energy_count, energy_mean, energy_metric_step, &
2476 0 : energy_shift, energy_weighted_sum, rayleigh_step, response_weight, &
2477 0 : shadow_energy_mean, shadow_energy_metric_step, shadow_energy_weighted_sum, &
2478 0 : shadow_rayleigh_step, rotation_count, rotation_norm, rotation_scale, &
2479 0 : schur_rhs, vy, vz, rayleigh_response, rhs, rotation_gradient, rotation_hessian, &
2480 0 : spin_rotation_gradient, spin_rotation_hessian, spin_rayleigh_response, &
2481 0 : previous_rotation_step, rotation_step, rotation_y, rotation_z, schur_block, &
2482 0 : solution, channel_rayleigh_response, channel_rotation_hessian, &
2483 0 : channel_rotation_gradient, channel_response_weight, channel_response_group, &
2484 14 : channel_generator, fixed_n_matrix, fixed_n_rhs, fixed_n_solution, density_modes)
2485 14 : RETURN
2486 : END IF
2487 :
2488 170 : rotation_norm(:) = 0.0_dp
2489 170 : rotation_count(:) = 0.0_dp
2490 412 : DO local_kpoint = 1, SIZE(kpoints%kp_env)
2491 242 : kp => kpoints%kp_env(local_kpoint)%kpoint_env
2492 242 : ikpoint = kp%nkpoint
2493 654 : DO ispin = 1, nspin_ot
2494 242 : rotation_group = MERGE(ispin, 1, nrotation_groups == nspin_ot)
2495 242 : local_channel = qs_ot_channel_index(ispin, ikpoint, nspin_ot) - first_channel + 1
2496 9110 : rotation_step(:, local_channel) = rotation_y(:, local_channel)
2497 498 : DO group = 1, ngroups
2498 : rotation_step(:, local_channel) = rotation_step(:, local_channel) - &
2499 : rotation_z(:, local_channel, group)* &
2500 10626 : fixed_n_solution(group, 1)
2501 : END DO
2502 : rotation_norm(rotation_group) = rotation_norm(rotation_group) + &
2503 : DOT_PRODUCT(rotation_step(:, local_channel), &
2504 9110 : rotation_step(:, local_channel))
2505 484 : rotation_count(rotation_group) = rotation_count(rotation_group) + REAL(nrotation, KIND=dp)
2506 : END DO
2507 : END DO
2508 170 : CALL para_env_inter_kp%sum(rotation_norm)
2509 170 : CALL para_env_inter_kp%sum(rotation_count)
2510 340 : rotation_scale(:) = 1.0_dp
2511 340 : DO group = 1, nrotation_groups
2512 340 : IF (rotation_norm(group) > EPSILON(rotation_norm(group))) THEN
2513 : rotation_scale(group) = MIN( &
2514 170 : 1.0_dp, rotation_trust_rms*SQRT(rotation_count(group)/rotation_norm(group)))
2515 : END IF
2516 : END DO
2517 :
2518 170 : energy_shift(:) = 0.0_dp
2519 170 : energy_count(:) = 0.0_dp
2520 412 : DO local_kpoint = 1, SIZE(kpoints%kp_env)
2521 242 : kp => kpoints%kp_env(local_kpoint)%kpoint_env
2522 242 : ikpoint = kp%nkpoint
2523 654 : DO ispin = 1, nspin_ot
2524 242 : rotation_group = MERGE(ispin, 1, nrotation_groups == nspin_ot)
2525 242 : local_channel = qs_ot_channel_index(ispin, ikpoint, nspin_ot) - first_channel + 1
2526 : rotation_step(:, local_channel) = &
2527 9110 : rotation_scale(rotation_group)*rotation_step(:, local_channel)
2528 : local_ot_env(local_channel)%ener_response_gx(:) = &
2529 : local_ot_env(local_channel)%ener_preconditioned_gx + &
2530 : MATMUL(channel_rayleigh_response(:, :, local_channel), &
2531 92758 : rotation_step(:, local_channel))
2532 1890 : DO i = 1, nenergy
2533 1648 : group = channel_response_group(i, local_channel)
2534 : energy_shift(group) = energy_shift(group) + &
2535 1648 : local_ot_env(local_channel)%ener_response_gx(i)
2536 1890 : energy_count(group) = energy_count(group) + 1.0_dp
2537 : END DO
2538 :
2539 1452 : ALLOCATE (dense_re(nmo, nmo), dense_im(nmo, nmo))
2540 242 : dense_re(:, :) = 0.0_dp
2541 242 : dense_im(:, :) = 0.0_dp
2542 242 : r = 0
2543 1508 : DO i = 1, nmo - 1
2544 5942 : DO j = i + 1, nmo
2545 4434 : r = r + 1
2546 4434 : dense_re(i, j) = rotation_step(r, local_channel)
2547 4434 : dense_re(j, i) = -rotation_step(r, local_channel)
2548 4434 : r = r + 1
2549 4434 : dense_im(i, j) = rotation_step(r, local_channel)
2550 5700 : dense_im(j, i) = rotation_step(r, local_channel)
2551 : END DO
2552 : END DO
2553 242 : CPASSERT(r == nrotation)
2554 : CALL qs_scf_loop_ot_dense_to_matrix( &
2555 : dense_re, local_ot_env(local_channel), &
2556 242 : local_ot_env(local_channel)%rot_mat_response_gx)
2557 : CALL qs_scf_loop_ot_dense_to_matrix( &
2558 : dense_im, local_ot_env(local_channel), &
2559 242 : local_ot_env(local_channel)%rot_mat_response_gx_im)
2560 242 : local_ot_env(local_channel)%rotation_response_valid = .TRUE.
2561 484 : DEALLOCATE (dense_im, dense_re)
2562 : END DO
2563 : END DO
2564 170 : CALL para_env_inter_kp%sum(energy_shift)
2565 170 : CALL para_env_inter_kp%sum(energy_count)
2566 412 : DO local_kpoint = 1, SIZE(kpoints%kp_env)
2567 242 : kp => kpoints%kp_env(local_kpoint)%kpoint_env
2568 242 : ikpoint = kp%nkpoint
2569 654 : DO ispin = 1, nspin_ot
2570 242 : local_channel = qs_ot_channel_index(ispin, ikpoint, nspin_ot) - first_channel + 1
2571 2132 : DO i = 1, nenergy
2572 1648 : group = channel_response_group(i, local_channel)
2573 1890 : IF (energy_count(group) > 0.0_dp) THEN
2574 : local_ot_env(local_channel)%ener_response_gx(i) = &
2575 : local_ot_env(local_channel)%ener_response_gx(i) - &
2576 1648 : energy_shift(group)/energy_count(group)
2577 : END IF
2578 : END DO
2579 : END DO
2580 : END DO
2581 :
2582 : ! The accepted gradient difference supplies the missing total cross curvature between the
2583 : ! finite rotation/energy candidate and the preceding accepted product direction. Keep the
2584 : ! dedicated endpoint history separate from CG's preconditioned-gradient history.
2585 170 : accepted_cross = 0.0_dp
2586 170 : accepted_cross_valid = .FALSE.
2587 170 : gradient_secant_failure = 0
2588 170 : IF (shadow_pending .AND. ABS(local_ot_env(1)%ds_min) > SQRT(EPSILON(1.0_dp))) THEN
2589 300 : DO local_kpoint = 1, SIZE(kpoints%kp_env)
2590 176 : kp => kpoints%kp_env(local_kpoint)%kpoint_env
2591 176 : ikpoint = kp%nkpoint
2592 476 : DO ispin = 1, nspin_ot
2593 176 : local_channel = qs_ot_channel_index(ispin, ikpoint, nspin_ot) - first_channel + 1
2594 : IF (local_ot_env(local_channel)%mermin_gradient_ref_valid .AND. &
2595 : ASSOCIATED(local_ot_env(local_channel)%rot_mat_mermin_g0) .AND. &
2596 176 : ASSOCIATED(local_ot_env(local_channel)%rot_mat_mermin_g0_im) .AND. &
2597 176 : ALLOCATED(local_ot_env(local_channel)%ener_mermin_g0)) THEN
2598 : CALL dbcsr_dot(local_ot_env(local_channel)%rot_mat_gx, &
2599 176 : local_ot_env(local_channel)%rot_mat_response_gx, tmp)
2600 176 : accepted_cross = accepted_cross + 0.5_dp*tmp
2601 : CALL dbcsr_dot(local_ot_env(local_channel)%rot_mat_mermin_g0, &
2602 176 : local_ot_env(local_channel)%rot_mat_response_gx, tmp)
2603 176 : accepted_cross = accepted_cross - 0.5_dp*tmp
2604 : CALL dbcsr_dot(local_ot_env(local_channel)%rot_mat_gx_im, &
2605 176 : local_ot_env(local_channel)%rot_mat_response_gx_im, tmp)
2606 176 : accepted_cross = accepted_cross + 0.5_dp*tmp
2607 : CALL dbcsr_dot(local_ot_env(local_channel)%rot_mat_mermin_g0_im, &
2608 176 : local_ot_env(local_channel)%rot_mat_response_gx_im, tmp)
2609 176 : accepted_cross = accepted_cross - 0.5_dp*tmp
2610 : accepted_cross = accepted_cross + &
2611 : DOT_PRODUCT( &
2612 : local_ot_env(local_channel)%ener_gx - &
2613 : local_ot_env(local_channel)%ener_mermin_g0, &
2614 1392 : local_ot_env(local_channel)%ener_response_gx)
2615 : ELSE
2616 0 : gradient_secant_failure = 1
2617 : END IF
2618 : END DO
2619 : END DO
2620 124 : CALL para_env_inter_kp%sum(accepted_cross)
2621 124 : CALL para_env_inter_kp%sum(gradient_secant_failure)
2622 124 : IF (gradient_secant_failure == 0) THEN
2623 124 : accepted_cross = accepted_cross/local_ot_env(1)%ds_min
2624 124 : accepted_cross_valid = ieee_is_finite(accepted_cross)
2625 : END IF
2626 : END IF
2627 : ! Evaluate the signed finite response in the same fixed-N tangent used to build the Schur
2628 : ! step. This is a prediction only; the minimizer calibrates it against accepted Mermin drops
2629 : ! before the candidate can be reused.
2630 170 : energy_weighted_sum(:) = 0.0_dp
2631 412 : DO local_kpoint = 1, SIZE(kpoints%kp_env)
2632 242 : kp => kpoints%kp_env(local_kpoint)%kpoint_env
2633 242 : ikpoint = kp%nkpoint
2634 654 : DO ispin = 1, nspin_ot
2635 242 : local_channel = qs_ot_channel_index(ispin, ikpoint, nspin_ot) - first_channel + 1
2636 2132 : DO i = 1, nenergy
2637 1648 : group = channel_response_group(i, local_channel)
2638 : energy_weighted_sum(group) = energy_weighted_sum(group) + &
2639 : channel_response_weight(i, local_channel)* &
2640 1890 : local_ot_env(local_channel)%ener_response_gx(i)
2641 : END DO
2642 : END DO
2643 : END DO
2644 170 : CALL para_env_inter_kp%sum(energy_weighted_sum)
2645 170 : energy_mean(:) = 0.0_dp
2646 354 : DO group = 1, ngroups
2647 354 : IF (ABS(curvature_sum(group)) > EPSILON(curvature_sum(group))) THEN
2648 184 : energy_mean(group) = energy_weighted_sum(group)/curvature_sum(group)
2649 : END IF
2650 : END DO
2651 :
2652 170 : response_curvature_total = 0.0_dp
2653 412 : DO local_kpoint = 1, SIZE(kpoints%kp_env)
2654 242 : kp => kpoints%kp_env(local_kpoint)%kpoint_env
2655 242 : ikpoint = kp%nkpoint
2656 654 : DO ispin = 1, nspin_ot
2657 242 : local_channel = qs_ot_channel_index(ispin, ikpoint, nspin_ot) - first_channel + 1
2658 1890 : DO i = 1, nenergy
2659 1648 : group = channel_response_group(i, local_channel)
2660 : energy_metric_step(i) = channel_response_weight(i, local_channel)* &
2661 : (local_ot_env(local_channel)%ener_response_gx(i) - &
2662 1890 : energy_mean(group))
2663 : END DO
2664 : rayleigh_step(:) = MATMUL(channel_rayleigh_response(:, :, local_channel), &
2665 89462 : rotation_step(:, local_channel))
2666 : response_curvature_total = response_curvature_total + &
2667 : DOT_PRODUCT(rotation_step(:, local_channel), &
2668 726 : MATMUL(channel_rotation_hessian(:, :, local_channel), &
2669 : rotation_step(:, local_channel))) - &
2670 : 2.0_dp*DOT_PRODUCT(rayleigh_step, energy_metric_step) + &
2671 509172 : DOT_PRODUCT(local_ot_env(local_channel)%ener_response_gx, energy_metric_step)
2672 : END DO
2673 : END DO
2674 170 : CALL para_env_inter_kp%sum(response_curvature_total)
2675 170 : local_ot_env(1)%response_model_curvature = response_curvature_total
2676 :
2677 170 : IF (shadow_pending) THEN
2678 124 : shadow_energy_weighted_sum(:) = 0.0_dp
2679 300 : DO local_kpoint = 1, SIZE(kpoints%kp_env)
2680 176 : kp => kpoints%kp_env(local_kpoint)%kpoint_env
2681 176 : ikpoint = kp%nkpoint
2682 476 : DO ispin = 1, nspin_ot
2683 176 : local_channel = qs_ot_channel_index(ispin, ikpoint, nspin_ot) - first_channel + 1
2684 1568 : DO i = 1, nenergy
2685 1216 : group = channel_response_group(i, local_channel)
2686 : shadow_energy_weighted_sum(group) = shadow_energy_weighted_sum(group) + &
2687 : channel_response_weight(i, local_channel)* &
2688 1392 : local_ot_env(local_channel)%ener_dx(i)
2689 : END DO
2690 : END DO
2691 : END DO
2692 124 : CALL para_env_inter_kp%sum(shadow_energy_weighted_sum)
2693 124 : shadow_energy_mean(:) = 0.0_dp
2694 258 : DO group = 1, ngroups
2695 258 : IF (ABS(curvature_sum(group)) > EPSILON(curvature_sum(group))) THEN
2696 134 : shadow_energy_mean(group) = shadow_energy_weighted_sum(group)/curvature_sum(group)
2697 : END IF
2698 : END DO
2699 :
2700 124 : shadow_curvature_total = 0.0_dp
2701 124 : projected_gradient(:) = 0.0_dp
2702 300 : DO local_kpoint = 1, SIZE(kpoints%kp_env)
2703 176 : kp => kpoints%kp_env(local_kpoint)%kpoint_env
2704 176 : ikpoint = kp%nkpoint
2705 476 : DO ispin = 1, nspin_ot
2706 176 : local_channel = qs_ot_channel_index(ispin, ikpoint, nspin_ot) - first_channel + 1
2707 1392 : DO i = 1, nenergy
2708 1216 : group = channel_response_group(i, local_channel)
2709 : energy_metric_step(i) = channel_response_weight(i, local_channel)* &
2710 : (local_ot_env(local_channel)%ener_response_gx(i) - &
2711 1216 : energy_mean(group))
2712 : shadow_energy_metric_step(i) = channel_response_weight(i, local_channel)* &
2713 : (local_ot_env(local_channel)%ener_dx(i) - &
2714 1392 : shadow_energy_mean(group))
2715 : END DO
2716 : rayleigh_step(:) = MATMUL(channel_rayleigh_response(:, :, local_channel), &
2717 67256 : rotation_step(:, local_channel))
2718 : shadow_rayleigh_step(:) = MATMUL( &
2719 : channel_rayleigh_response(:, :, local_channel), &
2720 67256 : previous_rotation_step(:, local_channel))
2721 : shadow_curvature_total = shadow_curvature_total + &
2722 : DOT_PRODUCT(previous_rotation_step(:, local_channel), &
2723 528 : MATMUL(channel_rotation_hessian(:, :, local_channel), &
2724 : previous_rotation_step(:, local_channel))) - &
2725 : 2.0_dp*DOT_PRODUCT(shadow_rayleigh_step, &
2726 : shadow_energy_metric_step) + &
2727 : DOT_PRODUCT(local_ot_env(local_channel)%ener_dx, &
2728 387616 : shadow_energy_metric_step)
2729 : response_cross = response_cross + &
2730 : DOT_PRODUCT(rotation_step(:, local_channel), &
2731 528 : MATMUL(channel_rotation_hessian(:, :, local_channel), &
2732 : previous_rotation_step(:, local_channel))) - &
2733 : DOT_PRODUCT(rayleigh_step, shadow_energy_metric_step) - &
2734 : DOT_PRODUCT(shadow_rayleigh_step, energy_metric_step) + &
2735 : 0.5_dp*( &
2736 : DOT_PRODUCT(local_ot_env(local_channel)%ener_response_gx, &
2737 : shadow_energy_metric_step) + &
2738 390048 : DOT_PRODUCT(local_ot_env(local_channel)%ener_dx, energy_metric_step))
2739 : projected_gradient(1) = projected_gradient(1) + &
2740 : DOT_PRODUCT(channel_rotation_gradient(:, local_channel), &
2741 : rotation_step(:, local_channel)) + &
2742 : DOT_PRODUCT(local_ot_env(local_channel)%ener_gx, &
2743 8064 : local_ot_env(local_channel)%ener_response_gx)
2744 : CALL dbcsr_dot(local_ot_env(local_channel)%matrix_gx, &
2745 176 : local_ot_env(local_channel)%matrix_dx, tmp)
2746 176 : projected_gradient(2) = projected_gradient(2) + tmp
2747 : CALL dbcsr_dot(local_ot_env(local_channel)%matrix_gx_im, &
2748 176 : local_ot_env(local_channel)%matrix_dx_im, tmp)
2749 176 : projected_gradient(2) = projected_gradient(2) + tmp
2750 : projected_gradient(2) = projected_gradient(2) + &
2751 : DOT_PRODUCT(channel_rotation_gradient(:, local_channel), &
2752 : previous_rotation_step(:, local_channel)) + &
2753 : DOT_PRODUCT(local_ot_env(local_channel)%ener_gx, &
2754 8416 : local_ot_env(local_channel)%ener_dx)
2755 : END DO
2756 : END DO
2757 124 : CALL para_env_inter_kp%sum(shadow_curvature_total)
2758 124 : CALL para_env_inter_kp%sum(response_cross)
2759 124 : CALL para_env_inter_kp%sum(projected_gradient)
2760 124 : local_ot_env(1)%response_shadow_curvature = shadow_curvature_total
2761 : END IF
2762 :
2763 : ! The frozen-H finite response above omits the self-consistent Hxc density response. Project
2764 : ! the most recent accepted physical density/Hamiltonian secant onto both the candidate and
2765 : ! conventional accepted directions. All spin/K-point overlaps are accumulated before the
2766 : ! minimum-Frobenius secant is formed: this is the direct-sum spin response, including the Hxc
2767 : ! cross-spin coupling, and not a sum of independently fitted spin Hessians. The resulting
2768 : ! sign-aware rank-two response is calibrated against accepted Mermin drops before it can replace
2769 : ! the conventional search direction.
2770 170 : hxc_density_norm_sq = 0.0_dp
2771 170 : hxc_response_work = 0.0_dp
2772 170 : hxc_density_overlap(:) = 0.0_dp
2773 170 : hxc_response_overlap(:) = 0.0_dp
2774 170 : failure = 0
2775 412 : DO local_kpoint = 1, SIZE(kpoints%kp_env)
2776 242 : kp => kpoints%kp_env(local_kpoint)%kpoint_env
2777 242 : ikpoint = kp%nkpoint
2778 242 : CPASSERT(ASSOCIATED(kp%ot_smat))
2779 242 : CPASSERT(SIZE(kp%ot_smat) >= 2)
2780 654 : DO ispin = 1, nspin_ot
2781 242 : local_channel = qs_ot_channel_index(ispin, ikpoint, nspin_ot) - first_channel + 1
2782 968 : DO energy_spin = 1, MERGE(nspin_energy, 1, restricted_roks)
2783 256 : physical_spin = MERGE(energy_spin, ispin, restricted_roks)
2784 256 : energy_start = MERGE((energy_spin - 1)*nmo + 1, 1, restricted_roks)
2785 256 : energy_stop = energy_start + nmo - 1
2786 256 : group = channel_response_group(energy_start, local_channel)
2787 : energy_metric_step(energy_start:energy_stop) = &
2788 : -channel_response_weight(energy_start:energy_stop, local_channel)* &
2789 : (local_ot_env(local_channel)%ener_response_gx(energy_start:energy_stop) - &
2790 1904 : energy_mean(group))
2791 : CALL qs_ot_density_tangent( &
2792 : channel_generator(:, :, local_channel), &
2793 : trial_occupation(:, ikpoint, physical_spin), wkp(ikpoint), &
2794 : rotation_step(:, local_channel), &
2795 256 : energy_metric_step(energy_start:energy_stop), density_modes(:, :, 1))
2796 256 : IF (.NOT. restricted_roks) THEN
2797 10572 : channel_density_modes(:, :, 1, local_channel) = density_modes(:, :, 1)
2798 : END IF
2799 13680 : density_modes(:, :, 2) = CMPLX(0.0_dp, 0.0_dp, KIND=dp)
2800 256 : IF (shadow_pending) THEN
2801 : energy_metric_step(energy_start:energy_stop) = &
2802 : -channel_response_weight(energy_start:energy_stop, local_channel)* &
2803 : (local_ot_env(local_channel)%ener_dx(energy_start:energy_stop) - &
2804 1402 : shadow_energy_mean(group))
2805 : CALL qs_ot_density_tangent( &
2806 : channel_generator(:, :, local_channel), &
2807 : trial_occupation(:, ikpoint, physical_spin), wkp(ikpoint), &
2808 : previous_rotation_step(:, local_channel), &
2809 186 : energy_metric_step(energy_start:energy_stop), density_modes(:, :, 2))
2810 186 : IF (.NOT. restricted_roks) THEN
2811 7970 : channel_density_modes(:, :, 2, local_channel) = density_modes(:, :, 2)
2812 : END IF
2813 : END IF
2814 : CALL qs_scf_kp_physical_secant_overlaps( &
2815 : local_ot_env(local_channel)%mermin_physical(energy_spin), &
2816 : kp%ot_smat(1), kp%ot_smat(2), density_modes, wkp(ikpoint), &
2817 : hxc_density_norm_sq_channel, hxc_response_work_channel, &
2818 256 : hxc_density_overlap_channel, hxc_response_overlap_channel, hxc_channel_valid)
2819 498 : IF (hxc_channel_valid) THEN
2820 242 : hxc_density_norm_sq = hxc_density_norm_sq + hxc_density_norm_sq_channel
2821 242 : hxc_response_work = hxc_response_work + hxc_response_work_channel
2822 726 : hxc_density_overlap(:) = hxc_density_overlap + hxc_density_overlap_channel
2823 726 : hxc_response_overlap(:) = hxc_response_overlap + hxc_response_overlap_channel
2824 : ELSE
2825 14 : failure = 1
2826 : END IF
2827 : END DO
2828 : END DO
2829 : END DO
2830 170 : exact_hxc_correction(:, :) = 0.0_dp
2831 170 : exact_hxc_valid = .FALSE.
2832 170 : IF (local_ot_env(1)%settings%ot_method == "CG" .AND. .NOT. restricted_roks) THEN
2833 : CALL qs_scf_kp_exact_hxc_projected( &
2834 : qs_env, kpoints, local_ot_env, first_channel, nspin_ot, wkp, para_env_inter_kp, &
2835 : sab_nl, fmwork, rho_ao_kp, channel_density_modes, MERGE(2, 1, shadow_pending), &
2836 144 : exact_hxc_correction, exact_hxc_valid)
2837 : END IF
2838 170 : CALL para_env_inter_kp%sum(failure)
2839 170 : CALL para_env_inter_kp%sum(hxc_density_norm_sq)
2840 170 : CALL para_env_inter_kp%sum(hxc_response_work)
2841 170 : CALL para_env_inter_kp%sum(hxc_density_overlap)
2842 170 : CALL para_env_inter_kp%sum(hxc_response_overlap)
2843 170 : hxc_correction(:, :) = 0.0_dp
2844 170 : hxc_valid = .FALSE.
2845 170 : IF (failure == 0) THEN
2846 156 : IF (shadow_pending .AND. &
2847 : ABS(local_ot_env(1)%ds_min) > SQRT(EPSILON(1.0_dp))) THEN
2848 : CALL qs_ot_density_secant_projected_hessian( &
2849 : hxc_density_norm_sq, hxc_response_work, hxc_density_overlap, &
2850 : hxc_response_overlap, hxc_correction, hxc_valid, &
2851 112 : secant_mode=2, secant_position=local_ot_env(1)%ds_min)
2852 : ELSE
2853 : CALL qs_ot_density_secant_projected_hessian( &
2854 : hxc_density_norm_sq, hxc_response_work, hxc_density_overlap, &
2855 44 : hxc_response_overlap, hxc_correction, hxc_valid)
2856 : END IF
2857 : END IF
2858 170 : IF (exact_hxc_valid) THEN
2859 98 : hxc_correction(:, :) = exact_hxc_correction
2860 98 : hxc_valid = .TRUE.
2861 : END IF
2862 170 : IF (hxc_valid) THEN
2863 : local_ot_env(1)%response_model_curvature = &
2864 156 : local_ot_env(1)%response_model_curvature + hxc_correction(1, 1)
2865 156 : IF (shadow_pending) THEN
2866 : local_ot_env(1)%response_shadow_curvature = &
2867 112 : local_ot_env(1)%response_shadow_curvature + hxc_correction(2, 2)
2868 : END IF
2869 112 : IF (shadow_pending) THEN
2870 : CALL ot_mermin_secant_curvature( &
2871 : local_ot_env(1)%response_reference_energy, local_ot_env(1)%etotal, &
2872 : local_ot_env(1)%response_predicted_slope, local_ot_env(1)%ds_min, &
2873 112 : accepted_curvature, accepted_curvature_valid)
2874 112 : projected_hessian(1, 1) = response_curvature_total
2875 112 : projected_hessian(1, 2) = response_cross
2876 112 : projected_hessian(2, 1) = response_cross
2877 112 : projected_hessian(2, 2) = shadow_curvature_total
2878 112 : IF (accepted_cross_valid) THEN
2879 112 : projected_hessian(1, 2) = accepted_cross - hxc_correction(1, 2)
2880 112 : projected_hessian(2, 1) = projected_hessian(1, 2)
2881 : END IF
2882 112 : IF (accepted_curvature_valid) THEN
2883 112 : projected_hessian(2, 2) = accepted_curvature - hxc_correction(2, 2)
2884 112 : local_ot_env(1)%response_shadow_curvature = accepted_curvature
2885 : END IF
2886 : CALL qs_ot_projected_response_update( &
2887 : projected_hessian, hxc_correction, projected_coefficients, projected_update_valid, &
2888 112 : projected_gradient=projected_gradient)
2889 112 : IF (projected_update_valid) THEN
2890 112 : rotation_norm(:) = 0.0_dp
2891 112 : rotation_count(:) = 0.0_dp
2892 264 : DO local_kpoint = 1, SIZE(kpoints%kp_env)
2893 152 : kp => kpoints%kp_env(local_kpoint)%kpoint_env
2894 152 : ikpoint = kp%nkpoint
2895 416 : DO ispin = 1, nspin_ot
2896 152 : rotation_group = MERGE(ispin, 1, nrotation_groups == nspin_ot)
2897 : local_channel = &
2898 152 : qs_ot_channel_index(ispin, ikpoint, nspin_ot) - first_channel + 1
2899 : rotation_step(:, local_channel) = &
2900 : projected_coefficients(1)*rotation_step(:, local_channel) + &
2901 5096 : projected_coefficients(2)*previous_rotation_step(:, local_channel)
2902 : local_ot_env(local_channel)%ener_response_gx(:) = &
2903 : projected_coefficients(1)* &
2904 : local_ot_env(local_channel)%ener_response_gx + &
2905 1152 : projected_coefficients(2)*local_ot_env(local_channel)%ener_dx
2906 : rotation_norm(rotation_group) = rotation_norm(rotation_group) + &
2907 : DOT_PRODUCT(rotation_step(:, local_channel), &
2908 5096 : rotation_step(:, local_channel))
2909 : rotation_count(rotation_group) = &
2910 304 : rotation_count(rotation_group) + REAL(nrotation, KIND=dp)
2911 : END DO
2912 : END DO
2913 112 : CALL para_env_inter_kp%sum(rotation_norm)
2914 112 : CALL para_env_inter_kp%sum(rotation_count)
2915 112 : response_scale = 1.0_dp
2916 224 : IF (SUM(rotation_norm) > EPSILON(SUM(rotation_norm))) THEN
2917 : response_scale = MIN( &
2918 : 1.0_dp, rotation_trust_rms* &
2919 336 : SQRT(SUM(rotation_count)/SUM(rotation_norm)))
2920 : END IF
2921 112 : unscaled_rotation_norm = 0.0_dp
2922 264 : DO local_channel = 1, nlocal_channels
2923 912 : ALLOCATE (dense_re(nmo, nmo), dense_im(nmo, nmo))
2924 152 : dense_re(:, :) = 0.0_dp
2925 152 : dense_im(:, :) = 0.0_dp
2926 152 : r = 0
2927 900 : DO i = 1, nmo - 1
2928 3372 : DO j = i + 1, nmo
2929 2472 : r = r + 1
2930 2472 : dense_re(i, j) = rotation_step(r, local_channel)
2931 2472 : dense_re(j, i) = -rotation_step(r, local_channel)
2932 2472 : r = r + 1
2933 2472 : dense_im(i, j) = rotation_step(r, local_channel)
2934 3220 : dense_im(j, i) = rotation_step(r, local_channel)
2935 : END DO
2936 : END DO
2937 152 : CPASSERT(r == nrotation)
2938 : unscaled_rotation_norm = MAX( &
2939 : unscaled_rotation_norm, qs_ot_antihermitian_spectral_norm( &
2940 6896 : CMPLX(dense_re, dense_im, KIND=dp)))
2941 264 : DEALLOCATE (dense_im, dense_re)
2942 : END DO
2943 112 : CALL para_env_inter_kp%max(unscaled_rotation_norm)
2944 112 : IF (unscaled_rotation_norm > EPSILON(unscaled_rotation_norm)) THEN
2945 : response_scale = MIN( &
2946 : response_scale, &
2947 112 : rotation_trust_spectral/unscaled_rotation_norm)
2948 : END IF
2949 264 : DO local_kpoint = 1, SIZE(kpoints%kp_env)
2950 152 : kp => kpoints%kp_env(local_kpoint)%kpoint_env
2951 152 : ikpoint = kp%nkpoint
2952 416 : DO ispin = 1, nspin_ot
2953 : local_channel = &
2954 152 : qs_ot_channel_index(ispin, ikpoint, nspin_ot) - first_channel + 1
2955 : rotation_step(:, local_channel) = &
2956 5096 : response_scale*rotation_step(:, local_channel)
2957 : local_ot_env(local_channel)%ener_response_gx(:) = &
2958 1152 : response_scale*local_ot_env(local_channel)%ener_response_gx
2959 152 : IF (.NOT. (ASSOCIATED(local_ot_env(local_channel)%matrix_response_gx))) THEN
2960 0 : CALL cp_abort(__LOCATION__, "Missing OT response gradient")
2961 : END IF
2962 152 : IF (.NOT. (ASSOCIATED(local_ot_env(local_channel)%matrix_response_gx_im))) THEN
2963 0 : CALL cp_abort(__LOCATION__, "Missing imaginary OT response gradient")
2964 : END IF
2965 : CALL dbcsr_copy(local_ot_env(local_channel)%matrix_response_gx, &
2966 152 : local_ot_env(local_channel)%matrix_dx)
2967 : CALL dbcsr_scale(local_ot_env(local_channel)%matrix_response_gx, &
2968 152 : response_scale*projected_coefficients(2))
2969 : CALL dbcsr_copy(local_ot_env(local_channel)%matrix_response_gx_im, &
2970 152 : local_ot_env(local_channel)%matrix_dx_im)
2971 : CALL dbcsr_scale(local_ot_env(local_channel)%matrix_response_gx_im, &
2972 152 : response_scale*projected_coefficients(2))
2973 912 : ALLOCATE (dense_re(nmo, nmo), dense_im(nmo, nmo))
2974 152 : dense_re(:, :) = 0.0_dp
2975 152 : dense_im(:, :) = 0.0_dp
2976 152 : r = 0
2977 900 : DO i = 1, nmo - 1
2978 3372 : DO j = i + 1, nmo
2979 2472 : r = r + 1
2980 2472 : dense_re(i, j) = rotation_step(r, local_channel)
2981 2472 : dense_re(j, i) = -rotation_step(r, local_channel)
2982 2472 : r = r + 1
2983 2472 : dense_im(i, j) = rotation_step(r, local_channel)
2984 3220 : dense_im(j, i) = rotation_step(r, local_channel)
2985 : END DO
2986 : END DO
2987 152 : CPASSERT(r == nrotation)
2988 : CALL qs_scf_loop_ot_dense_to_matrix( &
2989 : dense_re, local_ot_env(local_channel), &
2990 152 : local_ot_env(local_channel)%rot_mat_response_gx)
2991 : CALL qs_scf_loop_ot_dense_to_matrix( &
2992 : dense_im, local_ot_env(local_channel), &
2993 152 : local_ot_env(local_channel)%rot_mat_response_gx_im)
2994 304 : DEALLOCATE (dense_im, dense_re)
2995 : END DO
2996 : END DO
2997 112 : local_ot_env(1)%response_hxc_direction_valid = .TRUE.
2998 : local_ot_env(1)%response_model_curvature = response_scale**2* &
2999 : DOT_PRODUCT(projected_coefficients, &
3000 : MATMUL(projected_hessian + hxc_correction, &
3001 1680 : projected_coefficients))
3002 : END IF
3003 : END IF
3004 : END IF
3005 :
3006 0 : DEALLOCATE (coupling, curvature_sum, energy_count, energy_mean, energy_metric_step, &
3007 0 : energy_shift, energy_weighted_sum, rayleigh_step, response_weight, shadow_energy_mean, &
3008 0 : shadow_energy_metric_step, shadow_energy_weighted_sum, shadow_rayleigh_step, &
3009 0 : rotation_count, rotation_norm, rotation_scale, schur_rhs, vy, vz, &
3010 0 : rayleigh_response, rhs, rotation_gradient, rotation_hessian, &
3011 0 : spin_rotation_gradient, spin_rotation_hessian, spin_rayleigh_response, &
3012 0 : previous_rotation_step, rotation_step, &
3013 0 : rotation_y, rotation_z, schur_block, solution, channel_rayleigh_response, &
3014 0 : channel_rotation_hessian, channel_rotation_gradient, channel_response_weight, &
3015 0 : channel_density_modes, channel_response_group, channel_generator, &
3016 0 : fixed_n_matrix, fixed_n_rhs, &
3017 0 : fixed_n_solution, &
3018 170 : density_modes)
3019 :
3020 828 : END SUBROUTINE qs_scf_loop_prepare_ot_kp_rotation_response
3021 :
3022 : ! **************************************************************************************************
3023 : !> \brief copy a distributed OT matrix to a replicated dense matrix
3024 : !> \param matrix distributed matrix
3025 : !> \param ot_env channel owning the matrix distribution
3026 : !> \param dense replicated result
3027 : ! **************************************************************************************************
3028 1596 : SUBROUTINE qs_scf_loop_ot_matrix_to_dense(matrix, ot_env, dense)
3029 : TYPE(dbcsr_type), INTENT(IN) :: matrix
3030 : TYPE(qs_ot_type), INTENT(IN) :: ot_env
3031 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :), &
3032 : INTENT(OUT) :: dense
3033 :
3034 : INTEGER :: ncol, nrow
3035 : TYPE(cp_fm_struct_type), POINTER :: fm_struct
3036 : TYPE(cp_fm_type) :: fm
3037 :
3038 1596 : NULLIFY (fm_struct)
3039 0 : CPASSERT(ASSOCIATED(ot_env%blacs_env))
3040 1596 : CPASSERT(ASSOCIATED(ot_env%para_env))
3041 1596 : CALL dbcsr_get_info(matrix, nfullrows_total=nrow, nfullcols_total=ncol)
3042 6384 : ALLOCATE (dense(nrow, ncol))
3043 : CALL cp_fm_struct_create(fm_struct, context=ot_env%blacs_env, para_env=ot_env%para_env, &
3044 1596 : nrow_global=nrow, ncol_global=ncol)
3045 1596 : CALL cp_fm_create(fm, fm_struct)
3046 1596 : CALL copy_dbcsr_to_fm(matrix, fm)
3047 1596 : CALL cp_fm_get_submatrix(fm, dense, n_rows=nrow, n_cols=ncol)
3048 1596 : CALL cp_fm_release(fm)
3049 1596 : CALL cp_fm_struct_release(fm_struct)
3050 :
3051 4788 : END SUBROUTINE qs_scf_loop_ot_matrix_to_dense
3052 :
3053 : ! **************************************************************************************************
3054 : !> \brief copy a replicated dense matrix into a distributed OT matrix
3055 : !> \param dense replicated source
3056 : !> \param ot_env channel owning the matrix distribution
3057 : !> \param matrix distributed result
3058 : ! **************************************************************************************************
3059 788 : SUBROUTINE qs_scf_loop_ot_dense_to_matrix(dense, ot_env, matrix)
3060 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: dense
3061 : TYPE(qs_ot_type), INTENT(IN) :: ot_env
3062 : TYPE(dbcsr_type), INTENT(INOUT) :: matrix
3063 :
3064 : TYPE(cp_fm_struct_type), POINTER :: fm_struct
3065 : TYPE(cp_fm_type) :: fm
3066 :
3067 788 : NULLIFY (fm_struct)
3068 0 : CPASSERT(ASSOCIATED(ot_env%blacs_env))
3069 788 : CPASSERT(ASSOCIATED(ot_env%para_env))
3070 : CALL cp_fm_struct_create(fm_struct, context=ot_env%blacs_env, para_env=ot_env%para_env, &
3071 788 : nrow_global=SIZE(dense, 1), ncol_global=SIZE(dense, 2))
3072 788 : CALL cp_fm_create(fm, fm_struct)
3073 788 : CALL cp_fm_set_submatrix(fm, dense)
3074 788 : CALL copy_fm_to_dbcsr(fm, matrix)
3075 788 : CALL cp_fm_release(fm)
3076 788 : CALL cp_fm_struct_release(fm_struct)
3077 :
3078 788 : END SUBROUTINE qs_scf_loop_ot_dense_to_matrix
3079 :
3080 : ! **************************************************************************************************
3081 : !> \brief the inner loop of scf, specific to using to the orbital transformation method
3082 : !> basically, in goes the ks matrix out goes a new p matrix
3083 : !> \param qs_env ...
3084 : !> \param scf_env ...
3085 : !> \param smear ...
3086 : !> \param mos ...
3087 : !> \param rho ...
3088 : !> \param mo_derivs ...
3089 : !> \param total_energy ...
3090 : !> \param matrix_s ...
3091 : !> \param energy_only ...
3092 : !> \param has_unit_metric ...
3093 : !> \par History
3094 : !> 03.2006 created [Joost VandeVondele]
3095 : !> 2013 moved from qs_scf [Florian Schiffmann]
3096 : ! **************************************************************************************************
3097 81618 : SUBROUTINE qs_scf_loop_do_ot(qs_env, scf_env, smear, mos, rho, mo_derivs, total_energy, &
3098 : matrix_s, energy_only, has_unit_metric)
3099 :
3100 : TYPE(qs_environment_type), POINTER :: qs_env
3101 : TYPE(qs_scf_env_type), POINTER :: scf_env
3102 : TYPE(smear_type), POINTER :: smear
3103 : TYPE(mo_set_type), DIMENSION(:), INTENT(INOUT) :: mos
3104 : TYPE(qs_rho_type), POINTER :: rho
3105 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mo_derivs
3106 : REAL(KIND=dp), INTENT(IN) :: total_energy
3107 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
3108 : LOGICAL, INTENT(INOUT) :: energy_only
3109 : LOGICAL, INTENT(IN) :: has_unit_metric
3110 :
3111 : CHARACTER(LEN=*), PARAMETER :: routineN = 'qs_scf_loop_do_ot'
3112 :
3113 : INTEGER :: handle, ispin
3114 81618 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_ao
3115 : TYPE(dbcsr_type), POINTER :: orthogonality_metric
3116 :
3117 81618 : CALL timeset(routineN, handle)
3118 81618 : NULLIFY (rho_ao)
3119 :
3120 81618 : CALL qs_rho_get(rho, rho_ao=rho_ao)
3121 :
3122 81618 : IF (has_unit_metric) THEN
3123 18446 : NULLIFY (orthogonality_metric)
3124 : ELSE
3125 63172 : orthogonality_metric => matrix_s(1)%matrix
3126 : END IF
3127 :
3128 : ! in case of LSD the first spin qs_ot_env will drive the minimization
3129 : ! in the case of a restricted calculation, it will make sure the spin orbitals are equal
3130 :
3131 : CALL ot_scf_mini(mos, mo_derivs, smear, orthogonality_metric, &
3132 : total_energy, energy_only, scf_env%iter_delta, &
3133 81618 : scf_env%qs_ot_env)
3134 :
3135 179145 : DO ispin = 1, SIZE(mos)
3136 179145 : CALL set_mo_occupation(mo_set=mos(ispin), smear=smear)
3137 : END DO
3138 :
3139 179145 : DO ispin = 1, SIZE(mos)
3140 : CALL calculate_density_matrix(mos(ispin), &
3141 : rho_ao(ispin)%matrix, &
3142 179145 : use_dbcsr=.TRUE.)
3143 : END DO
3144 :
3145 81618 : scf_env%iter_method = scf_env%qs_ot_env(1)%OT_METHOD_FULL
3146 81618 : scf_env%iter_param = scf_env%qs_ot_env(1)%ds_min
3147 81618 : qs_env%broyden_adaptive_sigma = scf_env%qs_ot_env(1)%broyden_adaptive_sigma
3148 :
3149 81618 : CALL timestop(handle)
3150 :
3151 81618 : END SUBROUTINE qs_scf_loop_do_ot
3152 :
3153 : ! **************************************************************************************************
3154 : !> \brief Measure the distance between the diagonalized candidate density and the current density.
3155 : !> \param scf_env Holds SCF environment information.
3156 : !> \param rho Current electron density.
3157 : !> \param para_env Parallel environment.
3158 : !> \param delta Maximum matrix-element difference.
3159 : ! **************************************************************************************************
3160 148 : SUBROUTINE qs_scf_candidate_density_delta(scf_env, rho, para_env, delta)
3161 : TYPE(qs_scf_env_type), POINTER :: scf_env
3162 : TYPE(qs_rho_type), POINTER :: rho
3163 : TYPE(mp_para_env_type), POINTER :: para_env
3164 : REAL(KIND=dp), INTENT(OUT) :: delta
3165 :
3166 : INTEGER :: ic, ispin
3167 : REAL(KIND=dp) :: tmp
3168 148 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rho_ao_kp
3169 :
3170 148 : NULLIFY (rho_ao_kp)
3171 148 : CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp)
3172 :
3173 148 : CPASSERT(ASSOCIATED(scf_env%p_mix_new))
3174 148 : delta = 0.0_dp
3175 7910 : DO ic = 1, SIZE(scf_env%p_mix_new, 2)
3176 22642 : DO ispin = 1, SIZE(scf_env%p_mix_new, 1)
3177 : ! p_mix=1 leaves the candidate unchanged while cp_sm_mix evaluates
3178 : ! the same max-element density difference used by direct mixing.
3179 : CALL cp_sm_mix(m1=scf_env%p_mix_new(ispin, ic)%matrix, &
3180 : m2=rho_ao_kp(ispin, ic)%matrix, p_mix=1.0_dp, &
3181 14732 : delta=tmp, para_env=para_env)
3182 22494 : delta = MAX(delta, tmp)
3183 : END DO
3184 : END DO
3185 :
3186 148 : END SUBROUTINE qs_scf_candidate_density_delta
3187 :
3188 : ! **************************************************************************************************
3189 : !> \brief Commit the diagonalized candidate density without numerical mixing.
3190 : !> \param scf_env Holds SCF environment information.
3191 : !> \param rho Electron density to update.
3192 : ! **************************************************************************************************
3193 148 : SUBROUTINE qs_scf_commit_density_candidate(scf_env, rho)
3194 : TYPE(qs_scf_env_type), POINTER :: scf_env
3195 : TYPE(qs_rho_type), POINTER :: rho
3196 :
3197 : CHARACTER(len=default_string_length) :: name
3198 : INTEGER :: ic, ispin
3199 148 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rho_ao_kp
3200 :
3201 148 : NULLIFY (rho_ao_kp)
3202 148 : CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp)
3203 :
3204 148 : CPASSERT(ASSOCIATED(scf_env%p_mix_new))
3205 7910 : DO ic = 1, SIZE(scf_env%p_mix_new, 2)
3206 22642 : DO ispin = 1, SIZE(scf_env%p_mix_new, 1)
3207 14732 : CALL dbcsr_get_info(rho_ao_kp(ispin, ic)%matrix, name=name)
3208 22494 : CALL dbcsr_copy(rho_ao_kp(ispin, ic)%matrix, scf_env%p_mix_new(ispin, ic)%matrix, name=name)
3209 : END DO
3210 : END DO
3211 :
3212 148 : END SUBROUTINE qs_scf_commit_density_candidate
3213 :
3214 : ! **************************************************************************************************
3215 : !> \brief Performs the requested density mixing if any needed
3216 : !> \param scf_env Holds SCF environment information
3217 : !> \param rho All data for the electron density
3218 : !> \param para_env Parallel environment
3219 : !> \param diis_step Did we do a DIIS step?
3220 : ! **************************************************************************************************
3221 235537 : SUBROUTINE qs_scf_density_mixing(scf_env, rho, para_env, diis_step)
3222 : TYPE(qs_scf_env_type), POINTER :: scf_env
3223 : TYPE(qs_rho_type), POINTER :: rho
3224 : TYPE(mp_para_env_type), POINTER :: para_env
3225 : LOGICAL :: diis_step
3226 :
3227 235537 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rho_ao_kp
3228 :
3229 235537 : NULLIFY (rho_ao_kp)
3230 :
3231 235537 : CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp)
3232 :
3233 378506 : SELECT CASE (scf_env%mixing_method)
3234 : CASE (direct_mixing_nr)
3235 : CALL scf_env_density_mixing(scf_env%p_mix_new, &
3236 : scf_env%mixing_store, rho_ao_kp, para_env, scf_env%iter_delta, scf_env%iter_count, &
3237 142969 : diis=diis_step)
3238 : CASE (gspace_mixing_nr, pulay_mixing_nr, broyden_mixing_nr, modified_broyden_mixing_nr, &
3239 : multisecant_mixing_nr, new_pulay_mixing_nr)
3240 : ! Compute the difference p_out-p_in
3241 : CALL self_consistency_check(rho_ao_kp, scf_env%p_delta, para_env, scf_env%p_mix_new, &
3242 8286 : delta=scf_env%iter_delta)
3243 : CASE (no_mixing_nr)
3244 : CASE DEFAULT
3245 : CALL cp_abort(__LOCATION__, &
3246 : "unknown scf mixing method: "// &
3247 235537 : cp_to_string(scf_env%mixing_method))
3248 : END SELECT
3249 :
3250 235537 : END SUBROUTINE qs_scf_density_mixing
3251 :
3252 : ! **************************************************************************************************
3253 : !> \brief checks whether exit conditions for outer loop are satisfied
3254 : !> \param qs_env ...
3255 : !> \param scf_env ...
3256 : !> \param scf_control ...
3257 : !> \param should_stop ...
3258 : !> \param outer_loop_converged ...
3259 : !> \param exit_outer_loop ...
3260 : ! **************************************************************************************************
3261 26001 : SUBROUTINE qs_scf_check_outer_exit(qs_env, scf_env, scf_control, should_stop, &
3262 : outer_loop_converged, exit_outer_loop)
3263 : TYPE(qs_environment_type), POINTER :: qs_env
3264 : TYPE(qs_scf_env_type), POINTER :: scf_env
3265 : TYPE(scf_control_type), POINTER :: scf_control
3266 : LOGICAL :: should_stop, outer_loop_converged, &
3267 : exit_outer_loop
3268 :
3269 : REAL(KIND=dp) :: outer_loop_eps
3270 :
3271 26001 : outer_loop_converged = .TRUE.
3272 26001 : IF (scf_control%outer_scf%have_scf) THEN
3273 : ! We have an outer SCF loop...
3274 6197 : scf_env%outer_scf%iter_count = scf_env%outer_scf%iter_count + 1
3275 6197 : outer_loop_converged = .FALSE.
3276 :
3277 6197 : CALL outer_loop_gradient(qs_env, scf_env)
3278 : ! Multiple constraints: get largest deviation
3279 12482 : outer_loop_eps = SQRT(MAXVAL(scf_env%outer_scf%gradient(:, scf_env%outer_scf%iter_count)**2))
3280 :
3281 6197 : IF (outer_loop_eps < scf_control%outer_scf%eps_scf) outer_loop_converged = .TRUE.
3282 : END IF
3283 :
3284 : exit_outer_loop = should_stop .OR. outer_loop_converged .OR. &
3285 26001 : scf_env%outer_scf%iter_count > scf_control%outer_scf%max_scf
3286 :
3287 26001 : END SUBROUTINE qs_scf_check_outer_exit
3288 :
3289 : ! **************************************************************************************************
3290 : !> \brief checks whether exit conditions for inner loop are satisfied
3291 : !> \param qs_env ...
3292 : !> \param scf_env ...
3293 : !> \param scf_control ...
3294 : !> \param should_stop ...
3295 : !> \param just_energy ...
3296 : !> \param exit_inner_loop ...
3297 : !> \param inner_loop_converged ...
3298 : !> \param output_unit ...
3299 : ! **************************************************************************************************
3300 471370 : SUBROUTINE qs_scf_check_inner_exit(qs_env, scf_env, scf_control, should_stop, just_energy, &
3301 : exit_inner_loop, inner_loop_converged, output_unit)
3302 : TYPE(qs_environment_type), POINTER :: qs_env
3303 : TYPE(qs_scf_env_type), POINTER :: scf_env
3304 : TYPE(scf_control_type), POINTER :: scf_control
3305 : LOGICAL :: should_stop, just_energy, &
3306 : exit_inner_loop, inner_loop_converged
3307 : INTEGER :: output_unit
3308 :
3309 : LOGICAL :: davidson_tightening, update_validated
3310 :
3311 235685 : inner_loop_converged = .FALSE.
3312 235685 : exit_inner_loop = .FALSE.
3313 235685 : update_validated = scf_control%diagonalization%update_method /= diag_update_method_adiis
3314 235685 : IF (.NOT. update_validated) THEN
3315 148 : CPASSERT(ASSOCIATED(scf_env%scf_subspace_buffer))
3316 : update_validated = scf_env%adiis_validated .OR. &
3317 148 : scf_env%scf_subspace_buffer%diis_weight >= 1.0_dp
3318 : END IF
3319 :
3320 : CALL external_control(should_stop, "SCF", target_time=qs_env%target_time, &
3321 235685 : start_time=qs_env%start_time)
3322 235685 : IF (scf_env%iter_delta < scf_control%eps_scf .AND. update_validated) THEN
3323 : ! nested IFs on purpose: a .AND. chain leaves the evaluation order
3324 : ! unspecified. On non-Davidson runs the env array has size zero
3325 : ! while the tightening test below indexes its first channel
3326 22199 : davidson_tightening = .FALSE.
3327 22199 : IF (scf_env%method == block_davidson_diag_method_nr) THEN
3328 36 : IF (SIZE(scf_env%block_davidson_env) > 0) THEN
3329 : davidson_tightening = scf_env%block_davidson_env(1)%eps_iter_used > &
3330 : scf_env%block_davidson_env(1)%eps_iter .AND. &
3331 36 : scf_env%iter_delta < scf_env%block_davidson_env(1)%eps_iter_used
3332 : END IF
3333 : END IF
3334 : IF (davidson_tightening) THEN
3335 : ! the adaptive inner tolerance of the Davidson solve that produced
3336 : ! this mixing delta was looser than the delta itself. The apparent
3337 : ! convergence may be nothing but the noise floor of the loose
3338 : ! solve. Tighten the ratchet one decade below what that solve used
3339 : ! and keep iterating until the delta survives a tighter solve.
3340 0 : IF (output_unit > 0) THEN
3341 : WRITE (UNIT=output_unit, FMT="(T3,A,ES9.2,A)") &
3342 0 : "*** SCF delta below the loose Davidson tolerance ", &
3343 0 : scf_env%block_davidson_env(1)%eps_iter_used, " - tightening and continuing ***"
3344 : END IF
3345 : scf_env%block_davidson_env(1)%eps_iter_prev = &
3346 : MIN(scf_env%block_davidson_env(1)%eps_iter_prev, &
3347 : MAX(scf_env%block_davidson_env(1)%eps_iter, &
3348 0 : 0.1_dp*scf_env%block_davidson_env(1)%eps_iter_used))
3349 : ELSE
3350 22199 : IF (output_unit > 0) THEN
3351 : WRITE (UNIT=output_unit, FMT="(/,T3,A,I5,A/)") &
3352 11269 : "*** SCF run converged in ", scf_env%iter_count, " steps ***"
3353 : END IF
3354 22199 : inner_loop_converged = .TRUE.
3355 22199 : exit_inner_loop = .TRUE.
3356 : END IF
3357 213486 : ELSE IF (should_stop .OR. scf_env%iter_count >= scf_control%max_scf) THEN
3358 4118 : inner_loop_converged = .FALSE.
3359 4118 : IF (just_energy) THEN
3360 912 : exit_inner_loop = .FALSE.
3361 : ELSE
3362 3206 : exit_inner_loop = .TRUE.
3363 3206 : IF (output_unit > 0) THEN
3364 : WRITE (UNIT=output_unit, FMT="(/,T3,A,I5,A/)") &
3365 1604 : "Leaving inner SCF loop after reaching ", scf_env%iter_count, " steps."
3366 : END IF
3367 : END IF
3368 : END IF
3369 :
3370 235685 : END SUBROUTINE qs_scf_check_inner_exit
3371 :
3372 : ! **************************************************************************************************
3373 : !> \brief undoing density mixing. Important upon convergence
3374 : !> \param scf_env ...
3375 : !> \param rho ...
3376 : !> \param dft_control ...
3377 : !> \param para_env ...
3378 : !> \param diis_step ...
3379 : ! **************************************************************************************************
3380 25405 : SUBROUTINE qs_scf_undo_mixing(scf_env, rho, dft_control, para_env, diis_step)
3381 : TYPE(qs_scf_env_type), POINTER :: scf_env
3382 : TYPE(qs_rho_type), POINTER :: rho
3383 : TYPE(dft_control_type), POINTER :: dft_control
3384 : TYPE(mp_para_env_type), POINTER :: para_env
3385 : LOGICAL :: diis_step
3386 :
3387 : CHARACTER(len=default_string_length) :: name
3388 : INTEGER :: ic, ispin, nc
3389 25405 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rho_ao_kp
3390 :
3391 25405 : NULLIFY (rho_ao_kp)
3392 :
3393 25405 : IF (scf_env%mixing_method > 0) THEN
3394 17690 : CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp)
3395 17690 : nc = SIZE(scf_env%p_mix_new, 2)
3396 34544 : SELECT CASE (scf_env%mixing_method)
3397 : CASE (direct_mixing_nr)
3398 : CALL scf_env_density_mixing(scf_env%p_mix_new, scf_env%mixing_store, &
3399 : rho_ao_kp, para_env, scf_env%iter_delta, &
3400 : scf_env%iter_count, diis=diis_step, &
3401 16854 : invert=.TRUE.)
3402 162566 : DO ic = 1, nc
3403 325684 : DO ispin = 1, dft_control%nspins
3404 163118 : CALL dbcsr_get_info(rho_ao_kp(ispin, ic)%matrix, name=name) ! keep the name
3405 308830 : CALL dbcsr_copy(rho_ao_kp(ispin, ic)%matrix, scf_env%p_mix_new(ispin, ic)%matrix, name=name)
3406 : END DO
3407 : END DO
3408 : CASE (gspace_mixing_nr, pulay_mixing_nr, broyden_mixing_nr, modified_broyden_mixing_nr, &
3409 : multisecant_mixing_nr, new_pulay_mixing_nr)
3410 66542 : DO ic = 1, nc
3411 97572 : DO ispin = 1, dft_control%nspins
3412 48720 : CALL dbcsr_get_info(rho_ao_kp(ispin, ic)%matrix, name=name) ! keep the name
3413 96736 : CALL dbcsr_copy(rho_ao_kp(ispin, ic)%matrix, scf_env%p_mix_new(ispin, ic)%matrix, name=name)
3414 : END DO
3415 : END DO
3416 : END SELECT
3417 : END IF
3418 25405 : END SUBROUTINE qs_scf_undo_mixing
3419 :
3420 : ! **************************************************************************************************
3421 : !> \brief Performs the updates rho (takes care of mixing as well)
3422 : !> \param rho ...
3423 : !> \param qs_env ...
3424 : !> \param scf_env ...
3425 : !> \param ks_env ...
3426 : !> \param mix_rho ...
3427 : ! **************************************************************************************************
3428 235685 : SUBROUTINE qs_scf_rho_update(rho, qs_env, scf_env, ks_env, mix_rho)
3429 : TYPE(qs_rho_type), POINTER :: rho
3430 : TYPE(qs_environment_type), POINTER :: qs_env
3431 : TYPE(qs_scf_env_type), POINTER :: scf_env
3432 : TYPE(qs_ks_env_type), POINTER :: ks_env
3433 : LOGICAL, INTENT(IN) :: mix_rho
3434 :
3435 : TYPE(mp_para_env_type), POINTER :: para_env
3436 :
3437 235685 : NULLIFY (para_env)
3438 235685 : CALL get_qs_env(qs_env, para_env=para_env)
3439 : ! ** update qs_env%rho
3440 235685 : CALL qs_rho_update_rho(rho, qs_env=qs_env)
3441 : ! ** Density mixing through density matrix or on the reciprocal space grid (exclusive)
3442 235685 : IF (mix_rho) THEN
3443 : CALL gspace_mixing(qs_env, scf_env%mixing_method, scf_env%mixing_store, rho, &
3444 7450 : para_env, scf_env%iter_count)
3445 :
3446 : END IF
3447 235685 : CALL qs_ks_did_change(ks_env, rho_changed=.TRUE.)
3448 :
3449 235685 : END SUBROUTINE qs_scf_rho_update
3450 :
3451 : ! **************************************************************************************************
3452 : !> \brief Performs the necessary steps before leaving innner scf loop
3453 : !> \param scf_env ...
3454 : !> \param qs_env ...
3455 : !> \param diis_step ...
3456 : !> \param output_unit ...
3457 : ! **************************************************************************************************
3458 25405 : SUBROUTINE qs_scf_inner_finalize(scf_env, qs_env, diis_step, output_unit)
3459 : TYPE(qs_scf_env_type), POINTER :: scf_env
3460 : TYPE(qs_environment_type), POINTER :: qs_env
3461 : LOGICAL :: diis_step
3462 : INTEGER, INTENT(IN) :: output_unit
3463 :
3464 : LOGICAL :: do_kpoints
3465 : TYPE(dft_control_type), POINTER :: dft_control
3466 : TYPE(mp_para_env_type), POINTER :: para_env
3467 : TYPE(qs_energy_type), POINTER :: energy
3468 : TYPE(qs_ks_env_type), POINTER :: ks_env
3469 : TYPE(qs_rho_type), POINTER :: rho
3470 :
3471 25405 : NULLIFY (energy, rho, dft_control, ks_env)
3472 :
3473 : CALL get_qs_env(qs_env=qs_env, energy=energy, ks_env=ks_env, &
3474 : rho=rho, dft_control=dft_control, para_env=para_env, &
3475 25405 : do_kpoints=do_kpoints)
3476 :
3477 25405 : CALL cleanup_scf_loop(scf_env)
3478 :
3479 : ! now, print out energies and charges corresponding to the obtained wfn
3480 : ! (this actually is not 100% consistent at this point)!
3481 25405 : CALL qs_scf_print_summary(output_unit, qs_env)
3482 :
3483 25405 : CALL qs_scf_undo_mixing(scf_env, rho, dft_control, para_env, diis_step)
3484 :
3485 : ! *** update rspace rho since the mo changed
3486 : ! *** this might not always be needed (i.e. no post calculation / no forces )
3487 : ! *** but guarantees that rho and wfn are consistent at this point
3488 25405 : CALL qs_scf_rho_update(rho, qs_env, scf_env, ks_env, mix_rho=.FALSE.)
3489 :
3490 25405 : END SUBROUTINE qs_scf_inner_finalize
3491 :
3492 : ! **************************************************************************************************
3493 : !> \brief perform cleanup operations at the end of an scf loop
3494 : !> \param scf_env ...
3495 : !> \par History
3496 : !> 03.2006 created [Joost VandeVondele]
3497 : ! **************************************************************************************************
3498 25405 : SUBROUTINE cleanup_scf_loop(scf_env)
3499 : TYPE(qs_scf_env_type), INTENT(INOUT) :: scf_env
3500 :
3501 : CHARACTER(len=*), PARAMETER :: routineN = 'cleanup_scf_loop'
3502 :
3503 : INTEGER :: handle, ispin
3504 :
3505 25405 : CALL timeset(routineN, handle)
3506 :
3507 33108 : SELECT CASE (scf_env%method)
3508 : CASE (ot_method_nr)
3509 17107 : DO ispin = 1, SIZE(scf_env%qs_ot_env)
3510 17107 : CALL ot_scf_destroy(scf_env%qs_ot_env(ispin))
3511 : END DO
3512 7703 : DEALLOCATE (scf_env%qs_ot_env)
3513 : CASE (ot_diag_method_nr)
3514 : !
3515 : CASE (general_diag_method_nr)
3516 : !
3517 : CASE (special_diag_method_nr)
3518 : !
3519 : CASE (block_krylov_diag_method_nr, block_davidson_diag_method_nr)
3520 : !
3521 : CASE (filter_matrix_diag_method_nr)
3522 : !
3523 : CASE (smeagol_method_nr)
3524 : !
3525 : CASE DEFAULT
3526 : CALL cp_abort(__LOCATION__, &
3527 : "unknown scf method method:"// &
3528 25405 : cp_to_string(scf_env%method))
3529 : END SELECT
3530 :
3531 25405 : CALL timestop(handle)
3532 :
3533 25405 : END SUBROUTINE cleanup_scf_loop
3534 :
3535 594 : END MODULE qs_scf_loop_utils
|
