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