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