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