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