Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Adaptively Compressed Exchange (ACE) operator for HFX.
10 : !> Reference: Lin, J. Chem. Theory Comput. 2016, 12, 5, 2242-2249
11 : !>
12 : !> Algorithm (per spin):
13 : !>
14 : !> BUILD (first call, or every rebuild_frequency steps):
15 : !> 1. Full HFX: ks_matrix = K_HFX + H_core, energy%ex = E_x(full)
16 : !> 2. K_AO = ks_matrix - H_core (nao x nao, negative semidefinite)
17 : !> 3. C_occ = first nocc columns of mo_coeff, redistributed to a layout
18 : !> compatible with K_AO so that PDGEMM works correctly
19 : !> 4. xi = K_AO * C_occ (nao x nocc)
20 : !> 5. M = C_occ^T * xi (nocc x nocc, negative definite)
21 : !> 6. -M = U^T U via Cholesky (U upper triangular, stored in M_fm)
22 : !> 7. W = xi * U^{-1} (nao x nocc, the ACE projector)
23 : !> 8. Apply (see below) to update ks_matrix and energy%ex
24 : !>
25 : !> APPLY (all other steps):
26 : !> ks_matrix = H_core - W * W^T
27 : !> E_x = -0.5 * Tr[W^T * P * W]
28 : !>
29 : !> Diagnostics (controlled by DBG_STALE / DBG_EXACT_EX module flags):
30 : !>
31 : !> DIAG A projector staleness (cheap, runs every APPLY step)
32 : !> Computes ||W^T C_occ^current||_F / ||W^T C_occ^BUILD||_F.
33 : !> Ratio = 1 → W still accurate. Ratio -> 0 → W is stale.
34 : !>
35 : !> DIAG B exact vs ACE exchange energy (expensive: one full HFX per APPLY)
36 : !> Calls full HFX with just_energy=.TRUE. to get E_x^exact[P^k] and
37 : !> compares to E_x^ACE[P^k]. Growing |delta| confirms stale W.
38 : !> ACE ks_matrix and energy%ex are restored after the diagnostic.
39 : !>
40 : !> \author Ritama Kar
41 : ! **************************************************************************************************
42 :
43 : MODULE hfx_ace_methods
44 :
45 : USE admm_types, ONLY: admm_type,&
46 : get_admm_env
47 : USE bibliography, ONLY: Lin2016ACE,&
48 : cite_reference
49 : USE cp_blacs_env, ONLY: cp_blacs_env_type
50 : USE cp_control_types, ONLY: dft_control_type
51 : USE cp_dbcsr_api, ONLY: dbcsr_add,&
52 : dbcsr_copy,&
53 : dbcsr_create,&
54 : dbcsr_p_type,&
55 : dbcsr_release,&
56 : dbcsr_set,&
57 : dbcsr_type
58 : USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
59 : cp_dbcsr_plus_fm_fm_t
60 : USE cp_fm_basic_linalg, ONLY: cp_fm_scale,&
61 : cp_fm_trace,&
62 : cp_fm_triangular_multiply
63 : USE cp_fm_cholesky, ONLY: cp_fm_cholesky_decompose
64 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
65 : cp_fm_struct_release,&
66 : cp_fm_struct_type
67 : USE cp_fm_types, ONLY: cp_fm_create,&
68 : cp_fm_get_info,&
69 : cp_fm_release,&
70 : cp_fm_to_fm,&
71 : cp_fm_type
72 : USE cp_log_handling, ONLY: cp_get_default_logger,&
73 : cp_logger_get_default_io_unit,&
74 : cp_logger_type
75 : USE hfx_admm_utils, ONLY: hfx_ks_matrix
76 : USE hfx_types, ONLY: hfx_type
77 : USE input_section_types, ONLY: section_vals_type
78 : USE kinds, ONLY: dp
79 : USE message_passing, ONLY: mp_para_env_type
80 : USE parallel_gemm_api, ONLY: parallel_gemm
81 : USE pw_types, ONLY: pw_r3d_rs_type
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_mo_types, ONLY: get_mo_set,&
86 : mo_set_type
87 : USE qs_rho_types, ONLY: qs_rho_get,&
88 : qs_rho_type
89 : USE scf_control_types, ONLY: scf_control_type
90 : #include "./base/base_uses.f90"
91 :
92 : IMPLICIT NONE
93 : PRIVATE
94 :
95 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'hfx_ace_methods'
96 :
97 : ! -----------------------------------------------------------------------
98 : ! Module-level state: persists across SCF steps within one run.
99 : !
100 : ! ace_W(1,ispin) ACE projector W, shape nao x nocc
101 : ! ace_is_built .FALSE. until at least one successful BUILD
102 : ! ace_step_counter counts calls since last BUILD
103 : ! ace_W_ref_norm ||W^T C_occ^BUILD||_F stored at BUILD for DIAG A
104 : ! -----------------------------------------------------------------------
105 : TYPE(cp_fm_type), ALLOCATABLE, SAVE :: ace_W(:, :)
106 : LOGICAL, SAVE :: ace_is_built = .FALSE.
107 : INTEGER, SAVE :: ace_step_counter = 0
108 : REAL(dp), SAVE :: ace_W_ref_norm = 0.0_dp
109 : INTEGER, SAVE :: ace_geo_step = 0 ! NEW
110 :
111 : ! -----------------------------------------------------------------------
112 : ! Debug / diagnostic flags — set all .FALSE. for production.
113 : !
114 : ! DBG_ROUTING BUILD/APPLY/DEFER decisions, counters
115 : ! DBG_BUILD norms during BUILD
116 : ! DBG_ENERGY E_x(ACE) at every step; DIAG C on BUILD steps
117 : ! DBG_STALE DIAG A: staleness ratio at every APPLY step (cheap)
118 : ! DBG_EXACT_EX DIAG B: full HFX energy at every APPLY step (expensive)
119 : ! -----------------------------------------------------------------------
120 : LOGICAL, PARAMETER, PRIVATE :: DBG_ROUTING = .FALSE.
121 : LOGICAL, PARAMETER, PRIVATE :: DBG_BUILD = .FALSE.
122 : LOGICAL, PARAMETER, PRIVATE :: DBG_ENERGY = .FALSE.
123 : LOGICAL, PARAMETER, PRIVATE :: DBG_STALE = .FALSE.
124 : LOGICAL, PARAMETER, PRIVATE :: DBG_EXACT_EX = .FALSE.
125 :
126 : LOGICAL, SAVE :: ace_dynamic_mode = .FALSE.
127 : ! Set to .TRUE. by hfx_ace_set_dynamic_mode before geo_opt/MD starts.
128 : ! Stays .FALSE. for ENERGY/ENERGY_FORCE single-point runs.
129 :
130 : PUBLIC :: hfx_ace_ks_matrix, hfx_ace_release, hfx_ace_set_dynamic_mode
131 :
132 : CONTAINS
133 :
134 : ! **************************************************************************************************
135 : !> \brief Main ACE entry point, replacing hfx_ks_matrix in qs_ks_methods.
136 : !> \param qs_env ...
137 : !> \param ks_matrix ...
138 : !> \param rho ...
139 : !> \param energy ...
140 : !> \param calculate_forces ...
141 : !> \param just_energy ...
142 : !> \param v_rspace_new ...
143 : !> \param v_tau_rspace ...
144 : !> \param ace_rebuild_frequency ...
145 : !> \param ext_xc_section ...
146 : ! **************************************************************************************************
147 48 : SUBROUTINE hfx_ace_ks_matrix(qs_env, ks_matrix, rho, energy, &
148 : calculate_forces, just_energy, &
149 : v_rspace_new, v_tau_rspace, &
150 : ace_rebuild_frequency, ext_xc_section)
151 :
152 : TYPE(qs_environment_type), POINTER :: qs_env
153 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: ks_matrix
154 : TYPE(qs_rho_type), POINTER :: rho
155 : TYPE(qs_energy_type), POINTER :: energy
156 : LOGICAL, INTENT(IN) :: calculate_forces, just_energy
157 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: v_rspace_new, v_tau_rspace
158 : INTEGER, INTENT(IN) :: ace_rebuild_frequency
159 : TYPE(section_vals_type), OPTIONAL, POINTER :: ext_xc_section
160 :
161 : CHARACTER(LEN=*), PARAMETER :: routineN = 'hfx_ace_ks_matrix'
162 :
163 : INTEGER :: handle, iw, n_rep_hf, nspins, &
164 : rebuild_freq
165 : LOGICAL :: ace_built_now, rebuild_ace
166 : REAL(dp) :: ex_ace
167 : TYPE(cp_logger_type), POINTER :: logger
168 : TYPE(dft_control_type), POINTER :: dft_control
169 48 : TYPE(hfx_type), DIMENSION(:, :), POINTER :: x_data
170 : TYPE(scf_control_type), POINTER :: scf_control
171 :
172 48 : CALL timeset(routineN, handle)
173 :
174 48 : CALL cite_reference(Lin2016ACE)
175 48 : NULLIFY (logger, dft_control, x_data, scf_control)
176 :
177 48 : logger => cp_get_default_logger()
178 48 : iw = cp_logger_get_default_io_unit(logger)
179 :
180 48 : CALL get_qs_env(qs_env, x_data=x_data, dft_control=dft_control)
181 48 : n_rep_hf = SIZE(x_data, 1)
182 48 : nspins = dft_control%nspins
183 :
184 48 : IF (n_rep_hf /= 1) CPABORT("ACE: only one &HF section is supported.")
185 48 : IF (dft_control%nimages /= 1) THEN
186 0 : CPABORT("ACE: k-points / multiple images are not implemented.")
187 : END IF
188 :
189 : ! ACE requires explicit MO coefficients (C_occ) which are only available
190 : ! with diagonalization-based SCF. OT never constructs mo_coeff during
191 : ! the SCF, so the ACE projector build loop would silently get garbage.
192 48 : CALL get_qs_env(qs_env, scf_control=scf_control)
193 48 : IF (scf_control%use_ot) THEN
194 0 : CPABORT("ACE: OT doesn't work, use diagonalization-based SCF.")
195 : END IF
196 :
197 48 : rebuild_freq = MAX(1, ace_rebuild_frequency)
198 :
199 : ! ------------------------------------------------------------------
200 : ! Bypass A: energy-only call
201 : ! ------------------------------------------------------------------
202 48 : IF (just_energy) THEN
203 : IF (DBG_ROUTING .AND. iw > 0) THEN
204 : WRITE (iw, '(T2,A)') 'ACE | just_energy=T: full HFX (no matrix update)'
205 : END IF
206 : CALL hfx_call(qs_env, ks_matrix, rho, energy, &
207 : calculate_forces, just_energy, &
208 0 : v_rspace_new, v_tau_rspace, ext_xc_section)
209 0 : CALL timestop(handle)
210 8 : RETURN
211 : END IF
212 :
213 : ! ------------------------------------------------------------------
214 : ! Bypass B: ionic forces requested
215 : ! ------------------------------------------------------------------
216 48 : IF (calculate_forces) THEN
217 : IF (DBG_ROUTING .AND. iw > 0) THEN
218 : WRITE (iw, '(T2,A)') 'ACE | calculate_forces=T: full HFX for exact forces'
219 : END IF
220 : CALL hfx_call(qs_env, ks_matrix, rho, energy, &
221 : calculate_forces, just_energy, &
222 8 : v_rspace_new, v_tau_rspace, ext_xc_section)
223 8 : ace_is_built = .FALSE.
224 8 : ace_step_counter = 0
225 8 : ace_geo_step = ace_geo_step + 1 ! NEW: step 0 done, ACE active from now
226 8 : CALL timestop(handle)
227 8 : RETURN
228 : END IF
229 :
230 : ! ------------------------------------------------------------------
231 : ! Bypass C: first geometry step.
232 : !
233 : ! The ATOMIC initial guess gives C_occ far from self-consistency,
234 : ! which would produce an inaccurate W. Running full HFX for the
235 : ! entire first geometry step ensures that wavefunction extrapolation
236 : ! delivers a near-converged C_occ to geometry step 1, making the
237 : ! ACE BUILD there accurate from the first application.
238 : ! ------------------------------------------------------------------
239 40 : IF (ace_geo_step == 0 .AND. ace_dynamic_mode) THEN
240 0 : IF (iw > 0) WRITE (iw, '(T2,A)') &
241 0 : 'ACE | geo_step=0 (MD/GEO_OPT): full HFX for reference wavefunction'
242 : CALL hfx_call(qs_env, ks_matrix, rho, energy, &
243 : .FALSE., just_energy, &
244 0 : v_rspace_new, v_tau_rspace, ext_xc_section)
245 0 : CALL timestop(handle)
246 0 : RETURN
247 : END IF
248 :
249 : ! ------------------------------------------------------------------
250 : ! Rebuild decision
251 : ! ------------------------------------------------------------------
252 : rebuild_ace = (.NOT. ace_is_built) .OR. &
253 40 : (MOD(ace_step_counter, rebuild_freq) == 0)
254 :
255 : IF (DBG_ROUTING .AND. iw > 0) THEN
256 : WRITE (iw, '(/,T2,A)') REPEAT('-', 56)
257 : WRITE (iw, '(T2,A)') 'ACE | hfx_ace_ks_matrix'
258 : WRITE (iw, '(T4,A,L1)') 'ace_is_built = ', ace_is_built
259 : WRITE (iw, '(T4,A,L1)') 'rebuild_ace = ', rebuild_ace
260 : WRITE (iw, '(T4,A,I6)') 'step_counter = ', ace_step_counter
261 : WRITE (iw, '(T4,A,I6)') 'rebuild_freq = ', rebuild_freq
262 : WRITE (iw, '(T4,A,I4)') 'nspins = ', nspins
263 : WRITE (iw, '(T4,A)') MERGE('-> BUILD', '-> APPLY', rebuild_ace)
264 : WRITE (iw, '(T2,A)') REPEAT('-', 56)
265 : END IF
266 :
267 10 : IF (rebuild_ace) THEN
268 :
269 : ace_built_now = .FALSE.
270 : CALL hfx_ace_build_projector(qs_env, ks_matrix, rho, energy, &
271 : just_energy, &
272 : v_rspace_new, v_tau_rspace, &
273 : nspins, iw, ace_built_now, &
274 10 : ext_xc_section)
275 10 : IF (ace_built_now) THEN
276 8 : ace_is_built = .TRUE.
277 8 : ace_step_counter = 1
278 : IF (DBG_ROUTING .AND. iw > 0) THEN
279 : WRITE (iw, '(T4,A)') 'ACE | W built. Projector live from next step.'
280 : END IF
281 : ELSE
282 2 : ace_is_built = .FALSE.
283 2 : ace_step_counter = 0
284 : IF (DBG_ROUTING .AND. iw > 0) THEN
285 : WRITE (iw, '(T4,A)') 'ACE | Build deferred (C_occ=0). Full HFX in ks_matrix.'
286 : END IF
287 : END IF
288 :
289 : ELSE
290 :
291 30 : CALL hfx_ace_apply_projector(qs_env, ks_matrix, rho, energy, nspins, iw)
292 30 : ace_step_counter = ace_step_counter + 1
293 :
294 : ! ----------------------------------------------------------------
295 : ! DIAGNOSTIC B: compare E_x^ACE[P^k] with E_x^exact[P^k].
296 : !
297 : ! Calls full HFX (just_energy=.TRUE.) to get the exact exchange
298 : ! energy at the current ACE-converging density P^k. Compares
299 : ! with E_x^ACE[P^k] already stored in energy%ex.
300 : !
301 : ! Growing |delta| over the SCF confirms the root cause: W was
302 : ! built from C_occ^(step 1), which is far from the converged
303 : ! C_occ, so K_ACE = -WW^T no longer represents K_x accurately.
304 : !
305 : ! After the comparison, hfx_ace_apply_projector is called a
306 : ! second time to restore ks_matrix and energy%ex to ACE values
307 : ! so the SCF continues correctly. Cost: +2 full HFX per step.
308 : ! ----------------------------------------------------------------
309 : IF (DBG_EXACT_EX) THEN
310 : ex_ace = energy%ex
311 :
312 : CALL hfx_call(qs_env, ks_matrix, rho, energy, &
313 : .FALSE., .TRUE., &
314 : v_rspace_new, v_tau_rspace, ext_xc_section)
315 : IF (iw > 0) THEN
316 : WRITE (iw, '(/,T2,A)') REPEAT('-', 56)
317 : WRITE (iw, '(T2,A,I6)') 'ACE DIAG B | ace_step_counter = ', ace_step_counter
318 : WRITE (iw, '(T4,A,F20.10)') 'E_x(exact, P^k) = ', energy%ex
319 : WRITE (iw, '(T4,A,F20.10)') 'E_x(ACE, P^k) = ', ex_ace
320 : WRITE (iw, '(T4,A,ES12.4)') '|delta| = ', ABS(ex_ace - energy%ex)
321 : WRITE (iw, '(T4,A)') &
322 : '|delta|->0 on BUILD step; growth confirms stale projector'
323 : WRITE (iw, '(T2,A)') REPEAT('-', 56)
324 : END IF
325 :
326 : ! Restore ACE ks_matrix and energy%ex
327 : CALL hfx_ace_apply_projector(qs_env, ks_matrix, rho, energy, nspins, iw)
328 : END IF
329 :
330 : END IF
331 :
332 : IF (DBG_ROUTING .AND. iw > 0) THEN
333 : WRITE (iw, '(T4,A,F20.10)') 'energy%ex on exit = ', energy%ex
334 : WRITE (iw, '(T4,A,I6)') 'step_counter now = ', ace_step_counter
335 : END IF
336 :
337 40 : CALL timestop(handle)
338 :
339 48 : END SUBROUTINE hfx_ace_ks_matrix
340 :
341 : ! **************************************************************************************************
342 : !> \brief Build the ACE projector W.
343 : !> \param qs_env ...
344 : !> \param ks_matrix ...
345 : !> \param rho ...
346 : !> \param energy ...
347 : !> \param just_energy ...
348 : !> \param v_rspace_new ...
349 : !> \param v_tau_rspace ...
350 : !> \param nspins ...
351 : !> \param iw ...
352 : !> \param build_succeeded ...
353 : !> \param ext_xc_section ...
354 : ! **************************************************************************************************
355 10 : SUBROUTINE hfx_ace_build_projector(qs_env, ks_matrix, rho, energy, &
356 : just_energy, &
357 : v_rspace_new, v_tau_rspace, &
358 : nspins, iw, build_succeeded, &
359 : ext_xc_section)
360 :
361 : TYPE(qs_environment_type), POINTER :: qs_env
362 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: ks_matrix
363 : TYPE(qs_rho_type), POINTER :: rho
364 : TYPE(qs_energy_type), POINTER :: energy
365 : LOGICAL, INTENT(IN) :: just_energy
366 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: v_rspace_new, v_tau_rspace
367 : INTEGER, INTENT(IN) :: nspins, iw
368 : LOGICAL, INTENT(OUT) :: build_succeeded
369 : TYPE(section_vals_type), OPTIONAL, POINTER :: ext_xc_section
370 :
371 : CHARACTER(LEN=*), PARAMETER :: routineN = 'hfx_ace_build_projector'
372 :
373 : INTEGER :: handle, info_chol, ispin, nao, nmo, nocc
374 : LOGICAL :: do_admm
375 : REAL(dp) :: ehfx_full, frob
376 10 : REAL(dp), DIMENSION(:), POINTER :: occ_nums
377 : TYPE(admm_type), POINTER :: admm_env
378 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
379 : TYPE(cp_fm_struct_type), POINTER :: fmstruct
380 : TYPE(cp_fm_type) :: A_ref_fm, C_occ_fm, K_fm, M_fm, xi_fm
381 : TYPE(cp_fm_type), POINTER :: mo_coeff
382 10 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: ks_aux_fit
383 10 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_h
384 : TYPE(dbcsr_type) :: K_ao_dbcsr
385 : TYPE(dft_control_type), POINTER :: dft_control
386 10 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos, mos_for_ace
387 : TYPE(mp_para_env_type), POINTER :: para_env
388 :
389 : ! A_ref_fm: temporary nocc x nocc scratch for DIAG A reference norm
390 :
391 10 : CALL timeset(routineN, handle)
392 10 : NULLIFY (blacs_env, para_env, mos, mo_coeff, matrix_h, occ_nums, fmstruct)
393 10 : NULLIFY (dft_control, admm_env, mos_for_ace, ks_aux_fit)
394 :
395 10 : build_succeeded = .FALSE.
396 :
397 : CALL get_qs_env(qs_env, blacs_env=blacs_env, para_env=para_env, &
398 : mos=mos, matrix_h_kp=matrix_h, &
399 10 : dft_control=dft_control)
400 :
401 10 : do_admm = dft_control%do_admm
402 :
403 10 : IF (do_admm) THEN
404 6 : CALL get_qs_env(qs_env, admm_env=admm_env)
405 : CALL get_admm_env(admm_env, &
406 : matrix_ks_aux_fit=ks_aux_fit, &
407 6 : mos_aux_fit=mos_for_ace)
408 : ELSE
409 4 : mos_for_ace => mos
410 : END IF
411 :
412 : ! Step 1: full HFX
413 : CALL hfx_call(qs_env, ks_matrix, rho, energy, &
414 : .FALSE., just_energy, &
415 10 : v_rspace_new, v_tau_rspace, ext_xc_section)
416 10 : ehfx_full = energy%ex
417 :
418 : IF (DBG_BUILD .AND. iw > 0) THEN
419 : WRITE (iw, '(/,T2,A,F20.10)') 'ACE BUILD | E_x(full HFX) = ', ehfx_full
420 : END IF
421 :
422 : ! Allocate / resize module storage
423 10 : IF (ALLOCATED(ace_W)) THEN
424 2 : IF (SIZE(ace_W, 2) /= nspins) CALL hfx_ace_release()
425 : END IF
426 42 : IF (.NOT. ALLOCATED(ace_W)) ALLOCATE (ace_W(1, nspins))
427 :
428 : ! Reset reference norm; accumulated per spin in the loop below
429 10 : ace_W_ref_norm = 0.0_dp
430 :
431 : ! ----------------------------------------------------------------
432 : ! Per-spin build loop
433 : ! ----------------------------------------------------------------
434 18 : DO ispin = 1, nspins
435 :
436 10 : IF (mos_for_ace(ispin)%use_mo_coeff_b) THEN
437 : CALL copy_dbcsr_to_fm(mos_for_ace(ispin)%mo_coeff_b, &
438 0 : mos_for_ace(ispin)%mo_coeff)
439 : END IF
440 :
441 : CALL get_mo_set(mos_for_ace(ispin), mo_coeff=mo_coeff, &
442 : nao=nao, nmo=nmo, homo=nocc, &
443 10 : occupation_numbers=occ_nums)
444 :
445 10 : IF (nocc <= 0) CPABORT("ACE: homo <= 0.")
446 10 : IF (nocc > nmo) CPABORT("ACE: homo > nmo.")
447 10 : CPASSERT(ASSOCIATED(mo_coeff))
448 :
449 10 : CALL cp_fm_trace(mo_coeff, mo_coeff, frob)
450 :
451 : IF (DBG_BUILD .AND. iw > 0) THEN
452 : WRITE (iw, '(/,T2,A,I4)') 'ACE BUILD | ispin = ', ispin
453 : WRITE (iw, '(T4,A,I8)') 'nao = ', nao
454 : WRITE (iw, '(T4,A,I8)') 'nmo = ', nmo
455 : WRITE (iw, '(T4,A,I8)') 'nocc = ', nocc
456 : WRITE (iw, '(T4,A,L1)') 'use_mo_coeff_b = ', mos_for_ace(ispin)%use_mo_coeff_b
457 : WRITE (iw, '(T4,A,ES12.4)') '||mo_coeff||_F = ', SQRT(MAX(frob, 0.0_dp))
458 : END IF
459 :
460 10 : IF (frob < 1.0e-20_dp) THEN
461 : IF (DBG_BUILD .AND. iw > 0) THEN
462 : WRITE (iw, '(T4,A)') 'mo_coeff=0: build deferred to next step.'
463 : END IF
464 2 : CALL timestop(handle)
465 2 : RETURN
466 : END IF
467 :
468 : ! Step 2: K_AO
469 8 : IF (do_admm) THEN
470 : CALL dbcsr_create(K_ao_dbcsr, template=ks_aux_fit(ispin)%matrix, &
471 6 : name="K_ACE_aux")
472 6 : CALL dbcsr_copy(K_ao_dbcsr, ks_aux_fit(ispin)%matrix)
473 : ELSE
474 : CALL dbcsr_create(K_ao_dbcsr, template=ks_matrix(ispin, 1)%matrix, &
475 2 : name="K_AO")
476 2 : CALL dbcsr_copy(K_ao_dbcsr, ks_matrix(ispin, 1)%matrix)
477 2 : CALL dbcsr_add(K_ao_dbcsr, matrix_h(1, 1)%matrix, 1.0_dp, -1.0_dp)
478 : END IF
479 :
480 8 : NULLIFY (fmstruct)
481 : CALL cp_fm_struct_create(fmstruct, context=blacs_env, para_env=para_env, &
482 8 : nrow_global=nao, ncol_global=nao)
483 8 : CALL cp_fm_create(K_fm, fmstruct, name="K_dense")
484 8 : CALL cp_fm_struct_release(fmstruct)
485 8 : CALL copy_dbcsr_to_fm(K_ao_dbcsr, K_fm)
486 8 : CALL dbcsr_release(K_ao_dbcsr)
487 :
488 : IF (DBG_BUILD .AND. iw > 0) THEN
489 : CALL cp_fm_trace(K_fm, K_fm, frob)
490 : WRITE (iw, '(T4,A,ES12.4)') '||K_AO||_F = ', SQRT(MAX(frob, 0.0_dp))
491 : END IF
492 :
493 : ! Step 3: C_occ
494 8 : NULLIFY (fmstruct)
495 : CALL cp_fm_struct_create(fmstruct, context=blacs_env, para_env=para_env, &
496 8 : nrow_global=nao, ncol_global=nocc)
497 8 : CALL cp_fm_create(C_occ_fm, fmstruct, name="C_occ")
498 8 : CALL cp_fm_create(xi_fm, fmstruct, name="xi")
499 8 : CALL cp_fm_struct_release(fmstruct)
500 :
501 8 : CALL cp_fm_to_fm(mo_coeff, C_occ_fm)
502 :
503 : ! Step 4: xi = K_AO * C_occ
504 : CALL parallel_gemm('N', 'N', nao, nocc, nao, &
505 8 : 1.0_dp, K_fm, C_occ_fm, 0.0_dp, xi_fm)
506 8 : CALL cp_fm_release(K_fm)
507 :
508 : IF (DBG_BUILD .AND. iw > 0) THEN
509 : CALL cp_fm_trace(xi_fm, xi_fm, frob)
510 : WRITE (iw, '(T4,A,ES12.4)') '||xi||_F = ', SQRT(MAX(frob, 0.0_dp))
511 : END IF
512 :
513 : ! Step 5: M = C_occ^T * xi
514 8 : NULLIFY (fmstruct)
515 : CALL cp_fm_struct_create(fmstruct, context=blacs_env, para_env=para_env, &
516 8 : nrow_global=nocc, ncol_global=nocc)
517 8 : CALL cp_fm_create(M_fm, fmstruct, name="M")
518 8 : CALL cp_fm_struct_release(fmstruct)
519 :
520 : CALL parallel_gemm('T', 'N', nocc, nocc, nao, &
521 8 : 1.0_dp, C_occ_fm, xi_fm, 0.0_dp, M_fm)
522 8 : CALL cp_fm_release(C_occ_fm)
523 :
524 : IF (DBG_BUILD .AND. iw > 0) THEN
525 : CALL cp_fm_trace(M_fm, M_fm, frob)
526 : WRITE (iw, '(T4,A,ES12.4)') '||M||_F = ', SQRT(MAX(frob, 0.0_dp))
527 : END IF
528 :
529 : ! Step 6: Cholesky of -M = U^T U
530 8 : CALL cp_fm_scale(-1.0_dp, M_fm)
531 8 : CALL cp_fm_cholesky_decompose(M_fm, n=nocc, info_out=info_chol)
532 :
533 8 : IF (info_chol /= 0) THEN
534 0 : IF (iw > 0) THEN
535 0 : WRITE (iw, '(T4,A,I6)') 'ACE | Cholesky failed, info = ', info_chol
536 0 : WRITE (iw, '(T4,A,F20.10)') 'ACE | E_x(full) = ', ehfx_full
537 0 : WRITE (iw, '(T4,A,I8,A,I8)') 'ACE | nao=', nao, ' nocc=', nocc
538 : END IF
539 0 : CPABORT("ACE: Cholesky of -M failed (not positive definite).")
540 : END IF
541 :
542 : IF (DBG_BUILD .AND. iw > 0) THEN
543 : WRITE (iw, '(T4,A)') 'Cholesky OK (info=0).'
544 : END IF
545 :
546 : ! Step 7: W = xi * U^{-1}
547 8 : IF (ASSOCIATED(ace_W(1, ispin)%matrix_struct)) THEN
548 0 : CALL cp_fm_release(ace_W(1, ispin))
549 : END IF
550 :
551 8 : CALL cp_fm_create(ace_W(1, ispin), xi_fm%matrix_struct, name="W_ACE")
552 8 : CALL cp_fm_to_fm(xi_fm, ace_W(1, ispin))
553 :
554 : CALL cp_fm_triangular_multiply(M_fm, ace_W(1, ispin), &
555 : side='R', uplo_tr='U', &
556 : transpose_tr=.FALSE., &
557 : invert_tr=.TRUE., &
558 : n_rows=nao, n_cols=nocc, &
559 8 : alpha=1.0_dp)
560 :
561 8 : CALL cp_fm_release(xi_fm)
562 8 : CALL cp_fm_release(M_fm)
563 :
564 : IF (DBG_BUILD .AND. iw > 0) THEN
565 : CALL cp_fm_trace(ace_W(1, ispin), ace_W(1, ispin), frob)
566 : WRITE (iw, '(T4,A,I4,A,2I8,A,ES12.4)') &
567 : 'W spin=', ispin, ' shape=', nao, nocc, &
568 : ' ||W||_F=', SQRT(MAX(frob, 0.0_dp))
569 : END IF
570 :
571 : ! ----------------------------------------------------------------
572 : ! DIAG A reference: compute ||W^T C_occ^BUILD||_F for this spin.
573 : !
574 : ! C_occ_fm was released after step 5, but mo_coeff is still valid
575 : ! (it is a pointer into mos_for_ace, not allocated here).
576 : ! We re-create C_occ_fm from mo_coeff.
577 : !
578 : ! Theory: W^T C_occ^BUILD = U^{-T}(-U^T U) = -U → norm = ||U||_F
579 : ! (computed directly rather than storing U).
580 : ! ----------------------------------------------------------------
581 60 : IF (DBG_STALE) THEN
582 : NULLIFY (fmstruct)
583 : CALL cp_fm_struct_create(fmstruct, context=blacs_env, &
584 : para_env=para_env, &
585 : nrow_global=nao, ncol_global=nocc)
586 : CALL cp_fm_create(C_occ_fm, fmstruct, name="C_occ_ref_diag")
587 : CALL cp_fm_struct_release(fmstruct)
588 : CALL cp_fm_to_fm(mo_coeff, C_occ_fm)
589 :
590 : NULLIFY (fmstruct)
591 : CALL cp_fm_struct_create(fmstruct, context=blacs_env, &
592 : para_env=para_env, &
593 : nrow_global=nocc, ncol_global=nocc)
594 : CALL cp_fm_create(A_ref_fm, fmstruct, name="WtC_ref")
595 : CALL cp_fm_struct_release(fmstruct)
596 :
597 : CALL parallel_gemm('T', 'N', nocc, nocc, nao, &
598 : 1.0_dp, ace_W(1, ispin), C_occ_fm, 0.0_dp, A_ref_fm)
599 : CALL cp_fm_trace(A_ref_fm, A_ref_fm, frob)
600 : ace_W_ref_norm = ace_W_ref_norm + SQRT(MAX(frob, 0.0_dp))
601 :
602 : CALL cp_fm_release(C_occ_fm)
603 : CALL cp_fm_release(A_ref_fm)
604 : END IF
605 :
606 : END DO ! ispin
607 :
608 : IF (DBG_STALE .AND. iw > 0) THEN
609 : WRITE (iw, '(/,T2,A)') REPEAT('-', 56)
610 : WRITE (iw, '(T2,A)') 'ACE DIAG A | Reference norm stored at BUILD'
611 : WRITE (iw, '(T4,A,ES12.4)') &
612 : '||W^T C_occ^BUILD||_F (sum over spins) = ', ace_W_ref_norm
613 : WRITE (iw, '(T4,A)') &
614 : 'Staleness ratio = 1.0 at BUILD step; decreasing means W is becoming stale'
615 : WRITE (iw, '(T2,A)') REPEAT('-', 56)
616 : END IF
617 :
618 8 : build_succeeded = .TRUE.
619 :
620 : ! Step 8: apply immediately (DIAG C ratio printed via ehfx_full_ref)
621 : CALL hfx_ace_apply_projector(qs_env, ks_matrix, rho, energy, &
622 8 : nspins, iw, ehfx_full_ref=ehfx_full)
623 :
624 8 : CALL timestop(handle)
625 :
626 10 : END SUBROUTINE hfx_ace_build_projector
627 :
628 : ! **************************************************************************************************
629 : !> \brief Apply the stored ACE projector.
630 : !> \param qs_env ...
631 : !> \param ks_matrix ...
632 : !> \param rho ...
633 : !> \param energy ...
634 : !> \param nspins ...
635 : !> \param iw ...
636 : !> \param ehfx_full_ref ...
637 : ! **************************************************************************************************
638 38 : SUBROUTINE hfx_ace_apply_projector(qs_env, ks_matrix, rho, energy, &
639 : nspins, iw, ehfx_full_ref)
640 :
641 : TYPE(qs_environment_type), POINTER :: qs_env
642 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: ks_matrix
643 : TYPE(qs_rho_type), POINTER :: rho
644 : TYPE(qs_energy_type), POINTER :: energy
645 : INTEGER, INTENT(IN) :: nspins, iw
646 : REAL(dp), INTENT(IN), OPTIONAL :: ehfx_full_ref
647 :
648 : CHARACTER(LEN=*), PARAMETER :: routineN = 'hfx_ace_apply_projector'
649 :
650 : INTEGER :: handle, ispin, nao, nao_d, nmo_d, nocc, &
651 : nocc_d
652 : LOGICAL :: do_admm
653 : REAL(dp) :: ehfx_ace, frob_A, stale_norm, trace_val
654 : TYPE(admm_type), POINTER :: admm_env
655 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
656 : TYPE(cp_fm_struct_type), POINTER :: fmstruct, fmstruct_diag
657 : TYPE(cp_fm_type) :: A_diag_fm, C_occ_diag, P_fm, PW_fm
658 : TYPE(cp_fm_type), POINTER :: mo_coeff_diag
659 38 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: ks_aux_fit, ks_aux_fit_hfx
660 38 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_h, rho_ao
661 : TYPE(dft_control_type), POINTER :: dft_control
662 38 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos_aux_diag, mos_diag
663 : TYPE(mp_para_env_type), POINTER :: para_env
664 : TYPE(qs_rho_type), POINTER :: rho_aux_fit
665 :
666 : ! ------------------------------------------------------------------
667 : ! DIAG A local variables
668 : ! stale_norm ||W^T C_occ^current||_F summed over spins
669 : ! frob_A scratch for cp_fm_trace
670 : ! C_occ_diag current C_occ redistributed to W layout
671 : ! A_diag_fm nocc x nocc overlap W^T C_occ^current
672 : ! mos_diag primary mos (non-ADMM path)
673 : ! mos_aux_diag auxiliary mos (ADMM path)
674 : ! mo_coeff_diag pointer to the relevant mo_coeff
675 : ! nao_d, nmo_d, nocc_d dimensions from get_mo_set
676 : ! ------------------------------------------------------------------
677 :
678 38 : CALL timeset(routineN, handle)
679 38 : NULLIFY (blacs_env, para_env, matrix_h, rho_ao, fmstruct)
680 38 : NULLIFY (dft_control, admm_env, ks_aux_fit, ks_aux_fit_hfx, rho_aux_fit)
681 38 : NULLIFY (mos_diag, mos_aux_diag, mo_coeff_diag, fmstruct_diag)
682 38 : CPASSERT(ALLOCATED(ace_W))
683 :
684 : CALL get_qs_env(qs_env, blacs_env=blacs_env, para_env=para_env, &
685 38 : matrix_h_kp=matrix_h, dft_control=dft_control)
686 38 : do_admm = dft_control%do_admm
687 :
688 38 : IF (do_admm) THEN
689 30 : CALL get_qs_env(qs_env, admm_env=admm_env)
690 : CALL get_admm_env(admm_env, &
691 : matrix_ks_aux_fit=ks_aux_fit, &
692 : matrix_ks_aux_fit_hfx=ks_aux_fit_hfx, &
693 30 : rho_aux_fit=rho_aux_fit)
694 30 : CALL qs_rho_get(rho_aux_fit, rho_ao_kp=rho_ao)
695 : ELSE
696 8 : CALL qs_rho_get(rho, rho_ao_kp=rho_ao)
697 : END IF
698 :
699 76 : DO ispin = 1, nspins
700 38 : CALL cp_fm_get_info(ace_W(1, ispin), nrow_global=nao, ncol_global=nocc)
701 :
702 76 : IF (do_admm) THEN
703 30 : CALL dbcsr_set(ks_matrix(ispin, 1)%matrix, 0.0_dp)
704 : CALL dbcsr_add(ks_matrix(ispin, 1)%matrix, &
705 30 : matrix_h(1, 1)%matrix, 1.0_dp, 1.0_dp)
706 30 : CALL dbcsr_set(ks_aux_fit(ispin)%matrix, 0.0_dp)
707 : CALL cp_dbcsr_plus_fm_fm_t( &
708 : sparse_matrix=ks_aux_fit(ispin)%matrix, &
709 : matrix_v=ace_W(1, ispin), &
710 : ncol=nocc, &
711 : alpha=-1.0_dp, &
712 30 : keep_sparsity=.TRUE.)
713 : CALL dbcsr_add(ks_aux_fit_hfx(ispin)%matrix, &
714 30 : ks_aux_fit(ispin)%matrix, 0.0_dp, 1.0_dp)
715 : ELSE
716 8 : CALL dbcsr_set(ks_matrix(ispin, 1)%matrix, 0.0_dp)
717 : CALL cp_dbcsr_plus_fm_fm_t( &
718 : sparse_matrix=ks_matrix(ispin, 1)%matrix, &
719 : matrix_v=ace_W(1, ispin), &
720 : ncol=nocc, &
721 : alpha=-1.0_dp, &
722 8 : keep_sparsity=.TRUE.)
723 : CALL dbcsr_add(ks_matrix(ispin, 1)%matrix, &
724 8 : matrix_h(1, 1)%matrix, 1.0_dp, 1.0_dp)
725 : END IF
726 : END DO
727 :
728 : ! Exchange energy: E_x = -0.5 * sum_spin Tr[ W^T * P * W ]
729 : ehfx_ace = 0.0_dp
730 76 : DO ispin = 1, nspins
731 38 : CALL cp_fm_get_info(ace_W(1, ispin), nrow_global=nao, ncol_global=nocc)
732 :
733 38 : NULLIFY (fmstruct)
734 : CALL cp_fm_struct_create(fmstruct, context=blacs_env, para_env=para_env, &
735 38 : nrow_global=nao, ncol_global=nao)
736 38 : CALL cp_fm_create(P_fm, fmstruct, name="P_dense")
737 38 : CALL cp_fm_struct_release(fmstruct)
738 38 : CALL copy_dbcsr_to_fm(rho_ao(ispin, 1)%matrix, P_fm)
739 :
740 38 : CALL cp_fm_create(PW_fm, ace_W(1, ispin)%matrix_struct, name="PW")
741 : CALL parallel_gemm('N', 'N', nao, nocc, nao, &
742 38 : 1.0_dp, P_fm, ace_W(1, ispin), 0.0_dp, PW_fm)
743 38 : CALL cp_fm_trace(ace_W(1, ispin), PW_fm, trace_val)
744 38 : ehfx_ace = ehfx_ace - 0.5_dp*trace_val
745 :
746 38 : CALL cp_fm_release(P_fm)
747 38 : CALL cp_fm_release(PW_fm)
748 :
749 190 : IF (DBG_ENERGY .AND. iw > 0) THEN
750 : WRITE (iw, '(T4,A,I4,A,F20.10)') &
751 : 'ispin=', ispin, ' E_x(ACE) += ', -0.5_dp*trace_val
752 : END IF
753 : END DO
754 :
755 38 : energy%ex = ehfx_ace
756 :
757 : ! DIAG C: BUILD-step consistency check (printed when ehfx_full_ref present)
758 : IF (DBG_ENERGY .AND. iw > 0) THEN
759 : WRITE (iw, '(T2,A,F20.10)') 'ACE | E_x(ACE) = ', ehfx_ace
760 : IF (PRESENT(ehfx_full_ref)) THEN
761 : WRITE (iw, '(T2,A,F20.10)') 'ACE | E_x(full) = ', ehfx_full_ref
762 : WRITE (iw, '(T2,A,ES12.4)') 'ACE | |delta| = ', ABS(ehfx_ace - ehfx_full_ref)
763 : WRITE (iw, '(T2,A)') '(|delta| should be ~0 on BUILD steps; small is good)'
764 : END IF
765 : END IF
766 :
767 : ! ----------------------------------------------------------------
768 : ! DIAG A: projector staleness check.
769 : !
770 : ! Computes ||W^T C_occ^current||_F (summed over spins) and divides
771 : ! by ace_W_ref_norm = ||W^T C_occ^BUILD||_F stored at BUILD time.
772 : !
773 : ! staleness_ratio:
774 : ! 1.0 → C_occ hasn't changed since BUILD; projector is fresh
775 : ! < 1 → C_occ has rotated; how much depends on the SCF dynamics
776 : ! → 0 → C_occ is orthogonal to the BUILD-time span; W is useless
777 : !
778 : ! For non-ADMM: C_occ comes from primary mos (nao_orb x nocc).
779 : ! For ADMM: C_occ comes from mos_aux_fit (nao_aux x nocc_aux),
780 : ! consistent with ace_W dimensions.
781 : ! ----------------------------------------------------------------
782 : IF (DBG_STALE .AND. ace_W_ref_norm > 0.0_dp) THEN
783 : stale_norm = 0.0_dp
784 : CALL get_qs_env(qs_env, mos=mos_diag)
785 :
786 : DO ispin = 1, nspins
787 : CALL cp_fm_get_info(ace_W(1, ispin), nrow_global=nao_d, ncol_global=nocc_d)
788 :
789 : IF (do_admm) THEN
790 : CALL get_admm_env(admm_env, mos_aux_fit=mos_aux_diag)
791 : IF (mos_aux_diag(ispin)%use_mo_coeff_b) THEN
792 : CALL copy_dbcsr_to_fm(mos_aux_diag(ispin)%mo_coeff_b, &
793 : mos_aux_diag(ispin)%mo_coeff)
794 : END IF
795 : CALL get_mo_set(mos_aux_diag(ispin), mo_coeff=mo_coeff_diag, &
796 : nao=nao_d, nmo=nmo_d, homo=nocc_d)
797 : ELSE
798 : IF (mos_diag(ispin)%use_mo_coeff_b) THEN
799 : CALL copy_dbcsr_to_fm(mos_diag(ispin)%mo_coeff_b, &
800 : mos_diag(ispin)%mo_coeff)
801 : END IF
802 : CALL get_mo_set(mos_diag(ispin), mo_coeff=mo_coeff_diag, &
803 : nao=nao_d, nmo=nmo_d, homo=nocc_d)
804 : END IF
805 :
806 : NULLIFY (fmstruct_diag)
807 : CALL cp_fm_struct_create(fmstruct_diag, context=blacs_env, &
808 : para_env=para_env, &
809 : nrow_global=nao_d, ncol_global=nocc_d)
810 : CALL cp_fm_create(C_occ_diag, fmstruct_diag, name="C_stale")
811 : CALL cp_fm_struct_release(fmstruct_diag)
812 : CALL cp_fm_to_fm(mo_coeff_diag, C_occ_diag)
813 :
814 : NULLIFY (fmstruct_diag)
815 : CALL cp_fm_struct_create(fmstruct_diag, context=blacs_env, &
816 : para_env=para_env, &
817 : nrow_global=nocc_d, ncol_global=nocc_d)
818 : CALL cp_fm_create(A_diag_fm, fmstruct_diag, name="WtC_stale")
819 : CALL cp_fm_struct_release(fmstruct_diag)
820 :
821 : CALL parallel_gemm('T', 'N', nocc_d, nocc_d, nao_d, &
822 : 1.0_dp, ace_W(1, ispin), C_occ_diag, 0.0_dp, A_diag_fm)
823 : CALL cp_fm_trace(A_diag_fm, A_diag_fm, frob_A)
824 : stale_norm = stale_norm + SQRT(MAX(frob_A, 0.0_dp))
825 :
826 : CALL cp_fm_release(C_occ_diag)
827 : CALL cp_fm_release(A_diag_fm)
828 : END DO
829 : IF (iw > 0) THEN
830 : WRITE (iw, '(/,T2,A)') REPEAT('-', 56)
831 : WRITE (iw, '(T2,A,I6)') 'ACE DIAG A | ace_step_counter = ', ace_step_counter
832 : WRITE (iw, '(T4,A,ES12.4)') '||W^T C_occ^current||_F = ', stale_norm
833 : WRITE (iw, '(T4,A,ES12.4)') '||W^T C_occ^BUILD||_F (ref) = ', ace_W_ref_norm
834 : WRITE (iw, '(T4,A,F10.6)') 'staleness ratio (1=fresh, 0=stale) = ', &
835 : stale_norm/MAX(ace_W_ref_norm, 1.0e-30_dp)
836 : WRITE (iw, '(T2,A)') REPEAT('-', 56)
837 : END IF
838 : END IF
839 :
840 38 : CALL timestop(handle)
841 :
842 38 : END SUBROUTINE hfx_ace_apply_projector
843 :
844 : ! **************************************************************************************************
845 : !> \brief Release all ACE storage and reset state flags.
846 : !> \param iw_opt ...
847 : ! **************************************************************************************************
848 0 : SUBROUTINE hfx_ace_release(iw_opt)
849 :
850 : INTEGER, INTENT(IN), OPTIONAL :: iw_opt
851 :
852 : INTEGER :: i, iw, j
853 :
854 0 : iw = -1
855 0 : IF (PRESENT(iw_opt)) iw = iw_opt
856 :
857 0 : IF (ALLOCATED(ace_W)) THEN
858 0 : DO j = 1, SIZE(ace_W, 2)
859 0 : DO i = 1, SIZE(ace_W, 1)
860 0 : IF (ASSOCIATED(ace_W(i, j)%matrix_struct)) CALL cp_fm_release(ace_W(i, j))
861 : END DO
862 : END DO
863 0 : DEALLOCATE (ace_W)
864 : END IF
865 :
866 0 : ace_is_built = .FALSE.
867 0 : ace_step_counter = 0
868 0 : ace_W_ref_norm = 0.0_dp
869 0 : ace_geo_step = 0
870 0 : ace_dynamic_mode = .FALSE. ! ADDED: reset dynamic mode on release, so it must be explicitly re-enabled for GEO_OPT/MD runs
871 :
872 0 : IF (iw > 0) WRITE (iw, '(T2,A)') 'ACE | storage released, counters reset'
873 :
874 0 : END SUBROUTINE hfx_ace_release
875 :
876 : ! **************************************************************************************************
877 : !> \brief Mark this run as dynamic (GEO_OPT/MD) so Bypass C fires for geo step 0.
878 : !> Call this once from the geo_opt or MD driver before the first SCF.
879 : !> \param is_dynamic .TRUE. for GEO_OPT/MD, .FALSE. to reset.
880 : ! **************************************************************************************************
881 : ! **************************************************************************************************
882 : !> \brief ...
883 : !> \param is_dynamic ...
884 : ! **************************************************************************************************
885 2646 : SUBROUTINE hfx_ace_set_dynamic_mode(is_dynamic)
886 : LOGICAL, INTENT(IN) :: is_dynamic
887 :
888 2646 : ace_dynamic_mode = is_dynamic
889 2646 : END SUBROUTINE hfx_ace_set_dynamic_mode
890 :
891 : ! **************************************************************************************************
892 : !> \brief Private helper: call hfx_ks_matrix with or without ext_xc_section.
893 : !> \param qs_env ...
894 : !> \param ks_matrix ...
895 : !> \param rho ...
896 : !> \param energy ...
897 : !> \param calculate_forces ...
898 : !> \param just_energy ...
899 : !> \param v_rspace_new ...
900 : !> \param v_tau_rspace ...
901 : !> \param ext_xc_section ...
902 : ! **************************************************************************************************
903 : ! **************************************************************************************************
904 : !> \brief ...
905 : !> \param qs_env ...
906 : !> \param ks_matrix ...
907 : !> \param rho ...
908 : !> \param energy ...
909 : !> \param calculate_forces ...
910 : !> \param just_energy ...
911 : !> \param v_rspace_new ...
912 : !> \param v_tau_rspace ...
913 : !> \param ext_xc_section ...
914 : ! **************************************************************************************************
915 18 : SUBROUTINE hfx_call(qs_env, ks_matrix, rho, energy, &
916 : calculate_forces, just_energy, &
917 : v_rspace_new, v_tau_rspace, ext_xc_section)
918 :
919 : TYPE(qs_environment_type), POINTER :: qs_env
920 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: ks_matrix
921 : TYPE(qs_rho_type), POINTER :: rho
922 : TYPE(qs_energy_type), POINTER :: energy
923 : LOGICAL, INTENT(IN) :: calculate_forces, just_energy
924 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: v_rspace_new, v_tau_rspace
925 : TYPE(section_vals_type), OPTIONAL, POINTER :: ext_xc_section
926 :
927 18 : IF (PRESENT(ext_xc_section)) THEN
928 : CALL hfx_ks_matrix(qs_env, ks_matrix, rho, energy, &
929 : calculate_forces, just_energy, &
930 : v_rspace_new, v_tau_rspace, &
931 18 : ext_xc_section=ext_xc_section)
932 : ELSE
933 : CALL hfx_ks_matrix(qs_env, ks_matrix, rho, energy, &
934 : calculate_forces, just_energy, &
935 0 : v_rspace_new, v_tau_rspace)
936 : END IF
937 :
938 18 : END SUBROUTINE hfx_call
939 :
940 : END MODULE hfx_ace_methods
|
