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