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