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