(git:d3d49ac)
Loading...
Searching...
No Matches
mtlr_u_j_methods.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 Driver for self-consistent minimum tracking linear response U and J calculations.
9!> \author Ziwei Chai
10!> \date 29.07.2026
11!> \version 1.0
12! **************************************************************************************************
13MODULE mtlr_u_j_methods
14
15 USE atomic_kind_types, ONLY: atomic_kind_type,&
16 get_atomic_kind
17 USE cp_control_types, ONLY: dft_control_type
18 USE cp_log_handling, ONLY: cp_get_default_logger,&
19 cp_logger_get_default_io_unit,&
20 cp_logger_type
21 USE force_env_methods, ONLY: force_env_calc_energy_force
22 USE force_env_types, ONLY: force_env_get,&
23 force_env_type
24 USE input_constants, ONLY: atomic_guess,&
25 restart_guess
26 USE input_section_types, ONLY: section_vals_get_subs_vals,&
27 section_vals_type,&
28 section_vals_val_get
29 USE kinds, ONLY: dp
30 USE physcon, ONLY: evolt
31 USE qs_environment_types, ONLY: get_qs_env
32 USE qs_kind_types, ONLY: qs_kind_type
33 USE scf_control_types, ONLY: scf_control_type
34#include "./base/base_uses.f90"
35
36 IMPLICIT NONE
37
38 PRIVATE
39 PUBLIC :: do_mtlr_u_j
40
41CONTAINS
42! **************************************************************************************************
43!> \brief Driver for self-consistent MTLR U/J iteration.
44!> Each outer iteration performs:
45!> 1) one standard ENERGY SCF
46!> 2) one MTLR evaluation of U and J
47!> 3) one update of the Hubbard parameters
48!> until U and J are converged.
49!> using a method based on Lowdin charges
50!> \f[Q = S^{1/2} P S^{1/2}\f]
51!> where \b P and \b S are the density and the
52!> overlap matrix, respectively.
53!> \param[in,out] force_env ...
54!> \date 29.07.2026
55!> \author Ziwei Chai
56!> \version 1.0
57! **************************************************************************************************
58 SUBROUTINE do_mtlr_u_j(force_env)
59
60 TYPE(force_env_type), INTENT(INOUT), POINTER :: force_env
61
62 CHARACTER(LEN=*), PARAMETER :: routinen = 'do_mtlr_u_j'
63
64 INTEGER :: handle, ikind, k, max_mtlr_iter, n, &
65 nkind, output_unit, p_iter, u_iter
66 INTEGER, DIMENSION(:), POINTER :: atom_list
67 LOGICAL :: any_dft_plus_u, any_mtlr_kind, &
68 converged, do_reference_scf, &
69 wfn_restart_file_explicit
70 LOGICAL, ALLOCATABLE, DIMENSION(:) :: mtlr_kind
71 REAL(kind=dp) :: delta_j, delta_u, denominator_minus, denominator_plus, eps_u_j_loop, &
72 fhxc_minus, fhxc_plus, intercept_minus, intercept_plus, l, max_delta_j, max_delta_u, &
73 std_j, std_minus, std_plus, std_u, sum_trq_minus, sum_trq_minus_x_trq_minus, &
74 sum_trq_minus_x_vhxc_minus, sum_trq_plus, sum_trq_plus_x_trq_plus, &
75 sum_trq_plus_x_vhxc_plus, sum_vhxc_minus, sum_vhxc_plus
76 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: j_new, j_old, perturbation_strength, &
77 trq_minus, trq_plus, u_new, u_old, &
78 vhxc_minus, vhxc_plus
79 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
80 TYPE(cp_logger_type), POINTER :: logger
81 TYPE(dft_control_type), POINTER :: dft_control
82 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
83 TYPE(scf_control_type), POINTER :: scf_control
84 TYPE(section_vals_type), POINTER :: dft_section, force_env_section
85
86 CALL timeset(routinen, handle)
87
88 NULLIFY (atom_list, qs_kind_set, dft_control, logger, dft_section, force_env_section)
89
90 logger => cp_get_default_logger()
91 output_unit = cp_logger_get_default_io_unit(logger)
92
93 cpassert(ASSOCIATED(force_env))
94 cpassert(ASSOCIATED(force_env%qs_env))
95
96 CALL get_qs_env(force_env%qs_env, &
97 qs_kind_set=qs_kind_set, &
98 dft_control=dft_control, &
99 scf_control=scf_control, &
100 atomic_kind_set=atomic_kind_set)
101
102 cpassert(ASSOCIATED(atomic_kind_set))
103 cpassert(ASSOCIATED(dft_control))
104 cpassert(ASSOCIATED(qs_kind_set))
105 cpassert(ASSOCIATED(scf_control))
106
107 nkind = SIZE(atomic_kind_set)
108 IF (SIZE(qs_kind_set) /= nkind) THEN
109 cpabort("The atomic-kind and Quickstep-kind arrays have inconsistent sizes.")
110 END IF
111 ALLOCATE (mtlr_kind(nkind))
112 mtlr_kind(:) = .false.
113 any_dft_plus_u = .false.
114 any_mtlr_kind = .false.
115 DO ikind = 1, nkind
116 IF (.NOT. ASSOCIATED(qs_kind_set(ikind)%dft_plus_u)) cycle
117 any_dft_plus_u = .true.
118 IF (qs_kind_set(ikind)%dft_plus_u%do_mtlr) THEN
119 mtlr_kind(ikind) = .true.
120 any_mtlr_kind = .true.
121 END IF
122 END DO
123 IF (.NOT. any_dft_plus_u) THEN
124 CALL cp_abort(__location__, "RUN_TYPE MTLR requires at least one active "// &
125 "DFT_PLUS_U section.")
126 END IF
127 IF (.NOT. any_mtlr_kind) THEN
128 CALL cp_abort(__location__, "RUN_TYPE MTLR requires an active "// &
129 "MINIMUM_TRACKING_LINEAR_RESPONSE subsection.")
130 END IF
131
132 CALL force_env_get(force_env, force_env_section=force_env_section)
133 dft_section => section_vals_get_subs_vals(force_env_section, "DFT")
134 CALL section_vals_val_get(dft_section, "WFN_RESTART_FILE_NAME", &
135 explicit=wfn_restart_file_explicit)
136 IF (wfn_restart_file_explicit) THEN
137 CALL cp_abort(__location__, "MTLR does not allow an explicit WFN_RESTART_FILE_NAME. "// &
138 "Remove this keyword; the reference WFN is managed internally.")
139 END IF
140
141 SELECT CASE (scf_control%density_guess)
142 CASE (restart_guess)
143 do_reference_scf = .true.
144 CASE (atomic_guess)
145 do_reference_scf = .false.
146 CASE DEFAULT
147 cpabort("MTLR requires SCF_GUESS RESTART or SCF_GUESS ATOMIC.")
148 END SELECT
149
150 IF (output_unit > 0) THEN
151 WRITE (unit=output_unit, fmt="(/,T2,78('='))")
152 WRITE (unit=output_unit, fmt="(T2,A)") &
153 "MTLR| SCF initialization settings"
154 WRITE (unit=output_unit, fmt="(T2,78('-'))")
155 IF (scf_control%density_guess == restart_guess) THEN
156 WRITE (unit=output_unit, fmt="(T2,A,T68,A12)") &
157 "MTLR| SCF initial guess:", "RESTART"
158 WRITE (unit=output_unit, fmt="(T2,A,T68,A12)") &
159 "MTLR| QS extrapolation:", "USE_GUESS"
160 WRITE (unit=output_unit, fmt="(T2,A)") &
161 "MTLR| A reference SCF will precede each U/J iteration."
162 WRITE (unit=output_unit, fmt="(T2,A)") &
163 "MTLR| Every perturbation SCF will restart from the reference WFN."
164 ELSE
165 WRITE (unit=output_unit, fmt="(T2,A,T68,A12)") &
166 "MTLR| SCF initial guess:", "ATOMIC"
167 WRITE (unit=output_unit, fmt="(T2,A,T68,A12)") &
168 "MTLR| QS extrapolation:", "USE_GUESS"
169 WRITE (unit=output_unit, fmt="(T2,A)") &
170 "MTLR| No separate reference SCF will be performed."
171 WRITE (unit=output_unit, fmt="(T2,A)") &
172 "MTLR| Every perturbation SCF will start from an atomic guess."
173 END IF
174 WRITE (unit=output_unit, fmt="(T2,78('='))")
175 END IF
176
177 ALLOCATE (u_new(nkind))
178 ALLOCATE (j_new(nkind))
179 ALLOCATE (u_old(nkind))
180 ALLOCATE (j_old(nkind))
181 u_new(:) = 0.0_dp
182 j_new(:) = 0.0_dp
183 u_old(:) = 0.0_dp
184 j_old(:) = 0.0_dp
185 converged = .false.
186 eps_u_j_loop = dft_control%eps_u_j_loop
187 max_mtlr_iter = dft_control%max_mtlr_iter
188 IF (dft_control%nspins /= 2) THEN
189 cpabort("Unrestricted KS has to be used (the number of spin channels should be 2).")
190 END IF
191 IF (max_mtlr_iter < 1) THEN
192 cpabort("MAX_MTLR_LOOP must be at least one.")
193 END IF
194 IF (eps_u_j_loop <= 0.0_dp) THEN
195 cpabort("EPS_U_J_LOOP must be positive.")
196 END IF
197
198 ! Ensure that the DFT+U+J machinery remains active during the
199 ! initial MTLR calculation, even when a parameter starts from zero.
200 DO ikind = 1, nkind
201 IF (.NOT. mtlr_kind(ikind)) cycle
202 IF (qs_kind_set(ikind)%dft_plus_u%u_minus_j == 0.0_dp) THEN
203 qs_kind_set(ikind)%dft_plus_u%u_minus_j = epsilon(1.0_dp)
204 END IF
205 IF (qs_kind_set(ikind)%dft_plus_u%hund_j == 0.0_dp) THEN
206 qs_kind_set(ikind)%dft_plus_u%hund_j = epsilon(1.0_dp)
207 END IF
208 j_old(ikind) = qs_kind_set(ikind)%dft_plus_u%hund_j
209 u_old(ikind) = qs_kind_set(ikind)%dft_plus_u%u_minus_j + &
210 qs_kind_set(ikind)%dft_plus_u%hund_j
211 END DO
212
213 DO u_iter = 1, max_mtlr_iter
214
215 dft_control%mtlr_dft_with_perturbation = .false.
216
217 IF (do_reference_scf) THEN
218 IF (output_unit > 0) THEN
219 WRITE (unit=output_unit, fmt="(/,T2,78('='))")
220 WRITE (unit=output_unit, fmt="(T2,A)") &
221 "MTLR| Starting the unperturbed reference SCF."
222 WRITE (unit=output_unit, fmt="(T2,A,T72,I8)") &
223 "MTLR| U/J iteration:", u_iter
224 WRITE (unit=output_unit, fmt="(T2,78('='))")
225 END IF
226 CALL force_env_calc_energy_force(force_env, calc_force=.false.)
227 END IF
228
229 DO ikind = 1, nkind
230
231 IF (.NOT. mtlr_kind(ikind)) cycle
232 CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom_list)
233 IF (.NOT. any(atom_list == qs_kind_set(ikind)%dft_plus_u%lr_atom)) THEN
234 cpabort("INDEX_PERTURBED_ATOM does not belong to the KIND containing the MTLR section.")
235 END IF
236 IF (.NOT. ALLOCATED( &
237 qs_kind_set(ikind)%dft_plus_u%perturbation_strength)) THEN
238 cpabort("MTLR target does not contain perturbation strengths.")
239 END IF
240
241 dft_control%mtlr_ikind = ikind
242
243 n = SIZE(qs_kind_set(ikind)%dft_plus_u%perturbation_strength)
244 IF (n < 3) THEN
245 cpabort("MTLR linear regression requires at least three perturbation strengths.")
246 END IF
247 l = real(n, dp)
248 ALLOCATE (perturbation_strength(n))
249 ALLOCATE (trq_plus(n))
250 ALLOCATE (vhxc_plus(n))
251 ALLOCATE (trq_minus(n))
252 ALLOCATE (vhxc_minus(n))
253 trq_plus = 0.0_dp
254 trq_minus = 0.0_dp
255 vhxc_plus = 0.0_dp
256 vhxc_minus = 0.0_dp
257 sum_trq_plus = 0.0_dp
258 sum_vhxc_plus = 0.0_dp
259 sum_trq_plus_x_trq_plus = 0.0_dp
260 sum_trq_plus_x_vhxc_plus = 0.0_dp
261 sum_trq_minus = 0.0_dp
262 sum_vhxc_minus = 0.0_dp
263 sum_trq_minus_x_trq_minus = 0.0_dp
264 sum_trq_minus_x_vhxc_minus = 0.0_dp
265 std_plus = 0.0_dp
266 std_minus = 0.0_dp
267 perturbation_strength(:) = qs_kind_set(ikind)%dft_plus_u%perturbation_strength(:)
268
269 DO p_iter = 1, n
270
271 IF (output_unit > 0) THEN
272 WRITE (unit=output_unit, fmt="(/,T2,78('='))")
273 WRITE (unit=output_unit, fmt="(T2,A,T72,I8)") &
274 "MTLR| U/J iteration:", u_iter
275 WRITE (unit=output_unit, fmt="(T2,A,T74,I1,A4,I1)") &
276 "MTLR| Perturbation SCF:", p_iter, " of ", n
277 WRITE (unit=output_unit, fmt="(T2,A,T72,I8)") &
278 "MTLR| Target KIND index:", ikind
279 WRITE (unit=output_unit, fmt="(T2,A,T72,I8)") &
280 "MTLR| Target atom index:", &
281 qs_kind_set(ikind)%dft_plus_u%lr_atom
282 WRITE (unit=output_unit, fmt="(T2,A,T63,F14.8,A3)") &
283 "MTLR| Perturbation strength:", &
284 perturbation_strength(p_iter)*evolt, " eV"
285 WRITE (unit=output_unit, fmt="(T2,78('='))")
286 END IF
287
288 dft_control%perturbation_strength = perturbation_strength(p_iter)
289 dft_control%mtlr_dft_with_perturbation = .true.
290
291 CALL force_env_calc_energy_force(force_env, calc_force=.false.)
292
293 trq_plus(p_iter) = dft_control%trq(1) + dft_control%trq(2)
294 trq_minus(p_iter) = dft_control%trq(1) - dft_control%trq(2)
295 vhxc_plus(p_iter) = dft_control%vhxc(1) + dft_control%vhxc(2)
296 vhxc_minus(p_iter) = dft_control%vhxc(1) - dft_control%vhxc(2)
297
298 END DO
299
300 !Calculate all the number about vhxc and trq.
301 DO k = 1, n
302 sum_trq_plus_x_vhxc_plus = sum_trq_plus_x_vhxc_plus + trq_plus(k)*vhxc_plus(k)
303 sum_trq_plus = sum_trq_plus + trq_plus(k)
304 sum_vhxc_plus = sum_vhxc_plus + vhxc_plus(k)
305 sum_trq_plus_x_trq_plus = sum_trq_plus_x_trq_plus + trq_plus(k)*trq_plus(k)
306 END DO
307
308 DO k = 1, n
309 sum_trq_minus_x_vhxc_minus = sum_trq_minus_x_vhxc_minus + trq_minus(k)*vhxc_minus(k)
310 sum_trq_minus = sum_trq_minus + trq_minus(k)
311 sum_vhxc_minus = sum_vhxc_minus + vhxc_minus(k)
312 sum_trq_minus_x_trq_minus = sum_trq_minus_x_trq_minus + trq_minus(k)*trq_minus(k)
313 END DO
314
315 denominator_plus = l*sum_trq_plus_x_trq_plus - sum_trq_plus**2.0_dp
316 denominator_minus = l*sum_trq_minus_x_trq_minus - sum_trq_minus**2.0_dp
317 IF (abs(denominator_plus) < 100.0_dp*epsilon(1.0_dp)) THEN
318 cpabort("MTLR regression is singular.")
319 END IF
320 IF (abs(denominator_minus) < 100.0_dp*epsilon(1.0_dp)) THEN
321 cpabort("MTLR regression is singular.")
322 END IF
323 fhxc_minus = (l*sum_trq_minus_x_vhxc_minus - sum_trq_minus*sum_vhxc_minus) &
324 /denominator_minus
325 fhxc_plus = (l*sum_trq_plus_x_vhxc_plus - sum_trq_plus*sum_vhxc_plus) &
326 /denominator_plus
327 intercept_minus = (sum_vhxc_minus - fhxc_minus*sum_trq_minus)/l
328 intercept_plus = (sum_vhxc_plus - fhxc_plus*sum_trq_plus)/l
329 DO k = 1, n
330 std_minus = std_minus + (vhxc_minus(k) - (trq_minus(k)*fhxc_minus + intercept_minus))**2/(l - 2)
331 END DO
332 DO k = 1, n
333 std_plus = std_plus + (vhxc_plus(k) - (trq_plus(k)*fhxc_plus + intercept_plus))**2/(l - 2)
334 END DO
335 std_minus = sqrt(std_minus)
336 std_plus = sqrt(std_plus)
337
338 u_new(ikind) = 0.5_dp*fhxc_plus
339 j_new(ikind) = -0.5_dp*fhxc_minus
340 delta_u = u_new(ikind) - u_old(ikind)
341 delta_j = j_new(ikind) - j_old(ikind)
342 std_u = 0.5_dp*std_plus
343 std_j = 0.5_dp*std_minus
344
345 IF (output_unit > 0) THEN
346 WRITE (unit=output_unit, fmt="(/,T2,78('='))")
347 WRITE (unit=output_unit, fmt="(T2,A,I0,A,I0,A)") &
348 "MTLR| U/J iteration ", u_iter, &
349 " results for KIND ", ikind, "."
350 WRITE (unit=output_unit, fmt="(T2,78('-'))")
351 WRITE (unit=output_unit, fmt="(T2,A)") &
352 "MTLR| Linear-regression data for Hubbard U"
353 WRITE (unit=output_unit, &
354 fmt="(T2,A3,T8,A10,T22,A12,T37,A13,T52,A13,T67,A13)") &
355 "Pt", &
356 "Pert.[eV]", &
357 "trq(+)", &
358 "vhxc(+)[eV]", &
359 "vhxc(+)-fit", &
360 "vhxc(+)-resid"
361 DO k = 1, n
362 WRITE (unit=output_unit, &
363 fmt="(T2,I3,T8,F10.4,T22,ES12.4,T37,ES13.5,"// &
364 " T52,ES13.5,T67,ES13.4)") &
365 k, &
366 perturbation_strength(k)*evolt, &
367 trq_plus(k), &
368 vhxc_plus(k)*evolt, &
369 (fhxc_plus*trq_plus(k) + intercept_plus)*evolt, &
370 (vhxc_plus(k) - &
371 (fhxc_plus*trq_plus(k) + intercept_plus))*evolt
372 END DO
373 WRITE (unit=output_unit, fmt="(T2,78('-'))")
374 WRITE (unit=output_unit, fmt="(T2,A)") &
375 "MTLR| Linear-regression data for Hund J"
376 WRITE (unit=output_unit, &
377 fmt="(T2,A3,T8,A10,T22,A12,T37,A13,T52,A13,T67,A13)") &
378 "Pt", &
379 "Pert.[eV]", &
380 "trq(-)", &
381 "vhxc(-)[eV]", &
382 "vhxc(-)-fit", &
383 "vhxc(-)-resid"
384 DO k = 1, n
385 WRITE (unit=output_unit, &
386 fmt="(T2,I3,T8,F10.4,T22,ES12.4,T37,ES13.5,"// &
387 " T52,ES13.5,T67,ES13.4)") &
388 k, &
389 perturbation_strength(k)*evolt, &
390 trq_minus(k), &
391 vhxc_minus(k)*evolt, &
392 (fhxc_minus*trq_minus(k) + intercept_minus)*evolt, &
393 (vhxc_minus(k) - &
394 (fhxc_minus*trq_minus(k) + intercept_minus))*evolt
395 END DO
396 WRITE (unit=output_unit, fmt="(T2,78('-'))")
397 WRITE (unit=output_unit, &
398 fmt="(T2,A,T22,A16,T43,A16,T64,A16)") &
399 "MTLR| Parameter", &
400 "Old [eV]", &
401 "New [eV]", &
402 "Change [eV]"
403 WRITE (unit=output_unit, fmt="(T2,78('-'))")
404 WRITE (unit=output_unit, &
405 fmt="(T2,A,T22,ES16.8,T43,ES16.8,T64,ES16.8)") &
406 "MTLR| Hubbard U", &
407 u_old(ikind)*evolt, &
408 u_new(ikind)*evolt, &
409 (u_new(ikind) - u_old(ikind))*evolt
410 WRITE (unit=output_unit, &
411 fmt="(T2,A,T22,ES16.8,T43,ES16.8,T64,ES16.8)") &
412 "MTLR| Hund J", &
413 j_old(ikind)*evolt, &
414 j_new(ikind)*evolt, &
415 (j_new(ikind) - j_old(ikind))*evolt
416 WRITE (unit=output_unit, fmt="(T2,78('-'))")
417 WRITE (unit=output_unit, fmt="(T2,A,T61,ES16.8,A3)") &
418 "MTLR| Absolute change in U:", &
419 abs(u_new(ikind) - u_old(ikind))*evolt, " eV"
420 WRITE (unit=output_unit, fmt="(T2,A,T61,ES16.8,A3)") &
421 "MTLR| Absolute change in J:", &
422 abs(j_new(ikind) - j_old(ikind))*evolt, " eV"
423 WRITE (unit=output_unit, fmt="(T2,A,T61,ES16.8,A3)") &
424 "MTLR| U fit residual:", &
425 std_u*evolt, " eV"
426 WRITE (unit=output_unit, fmt="(T2,A,T61,ES16.8,A3)") &
427 "MTLR| J fit residual:", &
428 std_j*evolt, " eV"
429 WRITE (unit=output_unit, fmt="(T2,78('='))")
430 END IF
431
432 DEALLOCATE (perturbation_strength)
433 DEALLOCATE (trq_plus)
434 DEALLOCATE (vhxc_plus)
435 DEALLOCATE (trq_minus)
436 DEALLOCATE (vhxc_minus)
437
438 END DO
439
440 DO ikind = 1, nkind
441 IF (.NOT. mtlr_kind(ikind)) cycle
442 qs_kind_set(ikind)%dft_plus_u%u_minus_j = u_new(ikind) - j_new(ikind)
443 qs_kind_set(ikind)%dft_plus_u%hund_j = j_new(ikind)
444 END DO
445
446 max_delta_u = maxval(abs(pack(u_new(:) - u_old(:), mtlr_kind)))
447 max_delta_j = maxval(abs(pack(j_new(:) - j_old(:), mtlr_kind)))
448 IF (max_delta_u < eps_u_j_loop .AND. max_delta_j < eps_u_j_loop) THEN
449 converged = .true.
450 END IF
451
452 IF (output_unit > 0) THEN
453 WRITE (unit=output_unit, fmt="(/,T2,78('='))")
454
455 WRITE (unit=output_unit, fmt="(T2,A,I0,A)") &
456 "MTLR| Iteration ", u_iter, &
457 " completed for all atomic kinds."
458
459 WRITE (unit=output_unit, fmt="(T2,78('-'))")
460
461 WRITE (unit=output_unit, fmt="(T2,A,ES16.8,A)") &
462 "MTLR| Maximum change in U: ", &
463 max_delta_u*evolt, " eV"
464
465 WRITE (unit=output_unit, fmt="(T2,A,ES16.8,A)") &
466 "MTLR| Maximum change in J: ", &
467 max_delta_j*evolt, " eV"
468
469 WRITE (unit=output_unit, fmt="(T2,A,ES16.8,A)") &
470 "MTLR| Convergence threshold: ", &
471 eps_u_j_loop*evolt, " eV"
472
473 WRITE (unit=output_unit, fmt="(T2,78('-'))")
474
475 IF (converged) THEN
476 WRITE (unit=output_unit, fmt="(T2,A)") &
477 "MTLR| U and J have converged."
478
479 WRITE (unit=output_unit, fmt="(T2,A)") &
480 "MTLR| Both maximum changes are below the convergence threshold."
481
482 ELSE IF (u_iter < max_mtlr_iter) THEN
483 WRITE (unit=output_unit, fmt="(T2,A)") &
484 "MTLR| U and J have not yet converged."
485
486 WRITE (unit=output_unit, fmt="(T2,A)") &
487 "MTLR| Proceeding to the next linear-response U/J iteration."
488
489 ELSE
490 WRITE (unit=output_unit, fmt="(T2,A)") &
491 "MTLR| U and J have not converged."
492
493 WRITE (unit=output_unit, fmt="(T2,A)") &
494 "MTLR| The maximum number of MTLR iterations has been reached."
495 END IF
496
497 WRITE (unit=output_unit, fmt="(T2,78('='))")
498 END IF
499
500 IF (converged) THEN
501 IF (output_unit > 0) THEN
502 WRITE (unit=output_unit, fmt="(/,T2,78('*'))")
503 WRITE (unit=output_unit, fmt="(T2,A,I0,A)") &
504 "MTLR| U and J converged after ", &
505 u_iter, " iterations."
506 WRITE (unit=output_unit, fmt="(T2,78('*'),/)")
507 END IF
508 ELSE IF (u_iter == max_mtlr_iter) THEN
509 IF (output_unit > 0) THEN
510 WRITE (unit=output_unit, fmt="(/,T2,78('*'))")
511 WRITE (unit=output_unit, fmt="(T2,A,I0,A)") &
512 "MTLR| U and J did not converge within the maximum of ", &
513 max_mtlr_iter, " iterations."
514 WRITE (unit=output_unit, fmt="(T2,A)") &
515 "MTLR| Results from the final iteration will be reported."
516 WRITE (unit=output_unit, fmt="(T2,78('*'),/)")
517 END IF
518 END IF
519
520 IF (converged .OR. u_iter == max_mtlr_iter) THEN
521 IF (output_unit > 0) THEN
522 WRITE (unit=output_unit, fmt="(/,T2,78('*'))")
523 WRITE (unit=output_unit, fmt="(T2,A,ES16.8,A)") &
524 "MTLR| Maximum change in U: ", &
525 max_delta_u*evolt, " eV"
526 WRITE (unit=output_unit, fmt="(T2,A,ES16.8,A)") &
527 "MTLR| Maximum change in J: ", &
528 max_delta_j*evolt, " eV"
529 WRITE (unit=output_unit, fmt="(T2,A,ES16.8,A)") &
530 "MTLR| Maximum change over U and J: ", &
531 max(max_delta_u, max_delta_j)*evolt, " eV"
532 WRITE (unit=output_unit, fmt="(T2,78('*'))")
533 DO ikind = 1, nkind
534 IF (.NOT. mtlr_kind(ikind)) cycle
535 WRITE (unit=output_unit, fmt="(/,T2,A,I0,A)") &
536 "MTLR| Final parameters for KIND ", ikind, ":"
537 WRITE (unit=output_unit, fmt="(T2,A,F16.8,A)") &
538 "MTLR| Calculated Hubbard U: ", &
539 u_new(ikind)*evolt, " eV"
540 WRITE (unit=output_unit, fmt="(T2,A,F16.8,A)") &
541 "MTLR| Calculated Hund J: ", &
542 j_new(ikind)*evolt, " eV"
543 WRITE (unit=output_unit, fmt="(T2,A)") &
544 "MTLR| Recommended CP2K input parameters:"
545 WRITE (unit=output_unit, fmt="(T2,A,F16.8)") &
546 "MTLR| U_MINUS_J [eV] ", &
547 (u_new(ikind) - j_new(ikind))*evolt
548 WRITE (unit=output_unit, fmt="(T2,A,F16.8)") &
549 "MTLR| J [eV] ", &
550 j_new(ikind)*evolt
551 END DO
552 WRITE (unit=output_unit, fmt="(/,T2,78('*'),/)")
553 END IF
554 END IF
555
556 IF (converged) THEN
557 EXIT
558 END IF
559
560 u_old(:) = u_new
561 j_old(:) = j_new
562
563 END DO
564
565 DO ikind = 1, nkind
566 IF (.NOT. mtlr_kind(ikind)) cycle
567 qs_kind_set(ikind)%dft_plus_u%u_minus_j = u_new(ikind) - j_new(ikind)
568 qs_kind_set(ikind)%dft_plus_u%hund_j = j_new(ikind)
569 END DO
570 dft_control%mtlr_dft_with_perturbation = .false.
571 dft_control%perturbation_strength = 0.0_dp
572
573 CALL force_env_calc_energy_force(force_env, calc_force=.false.)
574
575 DEALLOCATE (u_new)
576 DEALLOCATE (j_new)
577 DEALLOCATE (u_old)
578 DEALLOCATE (j_old)
579 DEALLOCATE (mtlr_kind)
580
581 CALL timestop(handle)
582
583 END SUBROUTINE do_mtlr_u_j
584
585END MODULE mtlr_u_j_methods
atomic_kind_types
Define the atomic kind types and their sub types.
Definition atomic_kind_types.F:23
atomic_kind_types::get_atomic_kind
subroutine, public get_atomic_kind(atomic_kind, fist_potential, element_symbol, name, mass, kind_number, natom, atom_list, rcov, rvdw, z, qeff, apol, cpol, mm_radius, shell, shell_active, damping)
Get attributes of an atomic kind.
Definition atomic_kind_types.F:138
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_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
cp_log_handling::cp_logger_get_default_io_unit
integer function, public cp_logger_get_default_io_unit(logger)
returns the unit nr for the ionode (-1 on all other processors) skips as well checks if the procs cal...
Definition cp_log_handling.F:523
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition cp_log_handling.F:234
force_env_methods
Interface for the force calculations.
Definition force_env_methods.F:17
force_env_methods::force_env_calc_energy_force
recursive subroutine, public force_env_calc_energy_force(force_env, calc_force, consistent_energies, skip_external_control, eval_energy_forces, require_consistent_energy_force, linres, calc_stress_tensor)
Interface routine for force and energy calculations.
Definition force_env_methods.F:218
force_env_types
Interface for the force calculations.
Definition force_env_types.F:18
force_env_types::force_env_get
recursive subroutine, public force_env_get(force_env, in_use, fist_env, qs_env, meta_env, fp_env, subsys, para_env, potential_energy, additional_potential, kinetic_energy, harmonic_shell, kinetic_shell, cell, sub_force_env, qmmm_env, qmmmx_env, eip_env, pwdft_env, globenv, input, force_env_section, method_name_id, root_section, mixed_env, nnp_env, embed_env, ipi_env)
returns various attributes about the force environment
Definition force_env_types.F:344
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:18
input_constants::atomic_guess
integer, parameter, public atomic_guess
Definition input_constants.F:212
input_constants::restart_guess
integer, parameter, public restart_guess
Definition input_constants.F:212
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_val_get
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Definition input_section_types.F:1052
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
mtlr_u_j_methods
Driver for self-consistent minimum tracking linear response U and J calculations.
Definition mtlr_u_j_methods.F:13
mtlr_u_j_methods::do_mtlr_u_j
subroutine, public do_mtlr_u_j(force_env)
Driver for self-consistent MTLR U/J iteration. Each outer iteration performs: 1) one standard ENERGY ...
Definition mtlr_u_j_methods.F:59
physcon
Definition of physical constants:
Definition physcon.F:68
physcon::evolt
real(kind=dp), parameter, public evolt
Definition physcon.F:183
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:534
qs_kind_types
Define the quickstep kind type and their sub types.
Definition qs_kind_types.F:23
scf_control_types
parameters that control an scf iteration
Definition scf_control_types.F:17
atomic_kind_types::atomic_kind_type
Provides all information about an atomic kind.
Definition atomic_kind_types.F:45
cp_control_types::dft_control_type
Definition cp_control_types.F:757
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
force_env_types::force_env_type
wrapper to abstract the force evaluation of the various methods
Definition force_env_types.F:145
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:127
qs_kind_types::qs_kind_type
Provides all information about a quickstep kind.
Definition qs_kind_types.F:190
scf_control_types::scf_control_type
Definition scf_control_types.F:126