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 : !> \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 : ! **************************************************************************************************
13 : MODULE 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 :
41 : CONTAINS
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 4 : 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 4 : 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 4 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: j_new, j_old, perturbation_strength, &
77 4 : trq_minus, trq_plus, u_new, u_old, &
78 4 : vhxc_minus, vhxc_plus
79 4 : 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 4 : 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 4 : CALL timeset(routineN, handle)
87 :
88 4 : NULLIFY (atom_list, qs_kind_set, dft_control, logger, dft_section, force_env_section)
89 :
90 4 : logger => cp_get_default_logger()
91 4 : output_unit = cp_logger_get_default_io_unit(logger)
92 :
93 4 : CPASSERT(ASSOCIATED(force_env))
94 4 : 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 4 : atomic_kind_set=atomic_kind_set)
101 :
102 4 : CPASSERT(ASSOCIATED(atomic_kind_set))
103 4 : CPASSERT(ASSOCIATED(dft_control))
104 4 : CPASSERT(ASSOCIATED(qs_kind_set))
105 4 : CPASSERT(ASSOCIATED(scf_control))
106 :
107 4 : nkind = SIZE(atomic_kind_set)
108 4 : IF (SIZE(qs_kind_set) /= nkind) THEN
109 0 : CPABORT("The atomic-kind and Quickstep-kind arrays have inconsistent sizes.")
110 : END IF
111 12 : ALLOCATE (mtlr_kind(nkind))
112 4 : mtlr_kind(:) = .FALSE.
113 4 : any_dft_plus_u = .FALSE.
114 4 : any_mtlr_kind = .FALSE.
115 8 : DO ikind = 1, nkind
116 4 : IF (.NOT. ASSOCIATED(qs_kind_set(ikind)%dft_plus_u)) CYCLE
117 4 : any_dft_plus_u = .TRUE.
118 8 : IF (qs_kind_set(ikind)%dft_plus_u%do_mtlr) THEN
119 4 : mtlr_kind(ikind) = .TRUE.
120 4 : any_mtlr_kind = .TRUE.
121 : END IF
122 : END DO
123 4 : IF (.NOT. any_dft_plus_u) THEN
124 : CALL cp_abort(__LOCATION__, "RUN_TYPE MTLR requires at least one active "// &
125 0 : "DFT_PLUS_U section.")
126 : END IF
127 4 : IF (.NOT. any_mtlr_kind) THEN
128 : CALL cp_abort(__LOCATION__, "RUN_TYPE MTLR requires an active "// &
129 0 : "MINIMUM_TRACKING_LINEAR_RESPONSE subsection.")
130 : END IF
131 :
132 4 : CALL force_env_get(force_env, force_env_section=force_env_section)
133 4 : dft_section => section_vals_get_subs_vals(force_env_section, "DFT")
134 : CALL section_vals_val_get(dft_section, "WFN_RESTART_FILE_NAME", &
135 4 : explicit=wfn_restart_file_explicit)
136 4 : IF (wfn_restart_file_explicit) THEN
137 : CALL cp_abort(__LOCATION__, "MTLR does not allow an explicit WFN_RESTART_FILE_NAME. "// &
138 0 : "Remove this keyword; the reference WFN is managed internally.")
139 : END IF
140 :
141 4 : SELECT CASE (scf_control%density_guess)
142 : CASE (restart_guess)
143 0 : do_reference_scf = .TRUE.
144 : CASE (atomic_guess)
145 0 : do_reference_scf = .FALSE.
146 : CASE DEFAULT
147 4 : CPABORT("MTLR requires SCF_GUESS RESTART or SCF_GUESS ATOMIC.")
148 : END SELECT
149 :
150 4 : IF (output_unit > 0) THEN
151 2 : WRITE (UNIT=output_unit, FMT="(/,T2,78('='))")
152 : WRITE (UNIT=output_unit, FMT="(T2,A)") &
153 2 : "MTLR| SCF initialization settings"
154 2 : WRITE (UNIT=output_unit, FMT="(T2,78('-'))")
155 2 : IF (scf_control%density_guess == restart_guess) THEN
156 : WRITE (UNIT=output_unit, FMT="(T2,A,T68,A12)") &
157 2 : "MTLR| SCF initial guess:", "RESTART"
158 : WRITE (UNIT=output_unit, FMT="(T2,A,T68,A12)") &
159 2 : "MTLR| QS extrapolation:", "USE_GUESS"
160 : WRITE (UNIT=output_unit, FMT="(T2,A)") &
161 2 : "MTLR| A reference SCF will precede each U/J iteration."
162 : WRITE (UNIT=output_unit, FMT="(T2,A)") &
163 2 : "MTLR| Every perturbation SCF will restart from the reference WFN."
164 : ELSE
165 : WRITE (UNIT=output_unit, FMT="(T2,A,T68,A12)") &
166 0 : "MTLR| SCF initial guess:", "ATOMIC"
167 : WRITE (UNIT=output_unit, FMT="(T2,A,T68,A12)") &
168 0 : "MTLR| QS extrapolation:", "USE_GUESS"
169 : WRITE (UNIT=output_unit, FMT="(T2,A)") &
170 0 : "MTLR| No separate reference SCF will be performed."
171 : WRITE (UNIT=output_unit, FMT="(T2,A)") &
172 0 : "MTLR| Every perturbation SCF will start from an atomic guess."
173 : END IF
174 2 : WRITE (UNIT=output_unit, FMT="(T2,78('='))")
175 : END IF
176 :
177 12 : ALLOCATE (u_new(nkind))
178 8 : ALLOCATE (j_new(nkind))
179 8 : ALLOCATE (u_old(nkind))
180 8 : ALLOCATE (j_old(nkind))
181 4 : u_new(:) = 0.0_dp
182 4 : j_new(:) = 0.0_dp
183 4 : u_old(:) = 0.0_dp
184 4 : j_old(:) = 0.0_dp
185 4 : converged = .FALSE.
186 4 : eps_u_j_loop = dft_control%eps_u_j_loop
187 4 : max_mtlr_iter = dft_control%max_mtlr_iter
188 4 : IF (dft_control%nspins /= 2) THEN
189 0 : CPABORT("Unrestricted KS has to be used (the number of spin channels should be 2).")
190 : END IF
191 4 : IF (max_mtlr_iter < 1) THEN
192 0 : CPABORT("MAX_MTLR_LOOP must be at least one.")
193 : END IF
194 4 : IF (eps_u_j_loop <= 0.0_dp) THEN
195 0 : 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 8 : DO ikind = 1, nkind
201 4 : IF (.NOT. mtlr_kind(ikind)) CYCLE
202 4 : IF (qs_kind_set(ikind)%dft_plus_u%u_minus_j == 0.0_dp) THEN
203 0 : qs_kind_set(ikind)%dft_plus_u%u_minus_j = EPSILON(1.0_dp)
204 : END IF
205 4 : IF (qs_kind_set(ikind)%dft_plus_u%hund_j == 0.0_dp) THEN
206 0 : qs_kind_set(ikind)%dft_plus_u%hund_j = EPSILON(1.0_dp)
207 : END IF
208 4 : 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 8 : qs_kind_set(ikind)%dft_plus_u%hund_j
211 : END DO
212 :
213 10 : DO u_iter = 1, max_mtlr_iter
214 :
215 10 : dft_control%mtlr_dft_with_perturbation = .FALSE.
216 :
217 10 : IF (do_reference_scf) THEN
218 10 : IF (output_unit > 0) THEN
219 5 : WRITE (UNIT=output_unit, FMT="(/,T2,78('='))")
220 : WRITE (UNIT=output_unit, FMT="(T2,A)") &
221 5 : "MTLR| Starting the unperturbed reference SCF."
222 : WRITE (UNIT=output_unit, FMT="(T2,A,T72,I8)") &
223 5 : "MTLR| U/J iteration:", u_iter
224 5 : WRITE (UNIT=output_unit, FMT="(T2,78('='))")
225 : END IF
226 10 : CALL force_env_calc_energy_force(force_env, calc_force=.FALSE.)
227 : END IF
228 :
229 20 : DO ikind = 1, nkind
230 :
231 10 : IF (.NOT. mtlr_kind(ikind)) CYCLE
232 10 : CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom_list)
233 10 : IF (.NOT. ANY(atom_list == qs_kind_set(ikind)%dft_plus_u%lr_atom)) THEN
234 0 : CPABORT("INDEX_PERTURBED_ATOM does not belong to the KIND containing the MTLR section.")
235 : END IF
236 10 : IF (.NOT. ALLOCATED( &
237 : qs_kind_set(ikind)%dft_plus_u%perturbation_strength)) THEN
238 0 : CPABORT("MTLR target does not contain perturbation strengths.")
239 : END IF
240 :
241 10 : dft_control%mtlr_ikind = ikind
242 :
243 10 : n = SIZE(qs_kind_set(ikind)%dft_plus_u%perturbation_strength)
244 10 : IF (n < 3) THEN
245 0 : CPABORT("MTLR linear regression requires at least three perturbation strengths.")
246 : END IF
247 10 : l = REAL(n, dp)
248 30 : ALLOCATE (perturbation_strength(n))
249 20 : ALLOCATE (trq_plus(n))
250 20 : ALLOCATE (vhxc_plus(n))
251 20 : ALLOCATE (trq_minus(n))
252 20 : ALLOCATE (vhxc_minus(n))
253 10 : trq_plus = 0.0_dp
254 10 : trq_minus = 0.0_dp
255 10 : vhxc_plus = 0.0_dp
256 10 : vhxc_minus = 0.0_dp
257 10 : sum_trq_plus = 0.0_dp
258 10 : sum_vhxc_plus = 0.0_dp
259 10 : sum_trq_plus_x_trq_plus = 0.0_dp
260 10 : sum_trq_plus_x_vhxc_plus = 0.0_dp
261 10 : sum_trq_minus = 0.0_dp
262 10 : sum_vhxc_minus = 0.0_dp
263 10 : sum_trq_minus_x_trq_minus = 0.0_dp
264 10 : sum_trq_minus_x_vhxc_minus = 0.0_dp
265 10 : std_plus = 0.0_dp
266 10 : std_minus = 0.0_dp
267 60 : perturbation_strength(:) = qs_kind_set(ikind)%dft_plus_u%perturbation_strength(:)
268 :
269 60 : DO p_iter = 1, n
270 :
271 50 : IF (output_unit > 0) THEN
272 25 : WRITE (UNIT=output_unit, FMT="(/,T2,78('='))")
273 : WRITE (UNIT=output_unit, FMT="(T2,A,T72,I8)") &
274 25 : "MTLR| U/J iteration:", u_iter
275 : WRITE (UNIT=output_unit, FMT="(T2,A,T74,I1,A4,I1)") &
276 25 : "MTLR| Perturbation SCF:", p_iter, " of ", n
277 : WRITE (UNIT=output_unit, FMT="(T2,A,T72,I8)") &
278 25 : "MTLR| Target KIND index:", ikind
279 : WRITE (UNIT=output_unit, FMT="(T2,A,T72,I8)") &
280 25 : "MTLR| Target atom index:", &
281 50 : qs_kind_set(ikind)%dft_plus_u%lr_atom
282 : WRITE (UNIT=output_unit, FMT="(T2,A,T63,F14.8,A3)") &
283 25 : "MTLR| Perturbation strength:", &
284 50 : perturbation_strength(p_iter)*evolt, " eV"
285 25 : WRITE (UNIT=output_unit, FMT="(T2,78('='))")
286 : END IF
287 :
288 50 : dft_control%perturbation_strength = perturbation_strength(p_iter)
289 50 : dft_control%mtlr_dft_with_perturbation = .TRUE.
290 :
291 50 : CALL force_env_calc_energy_force(force_env, calc_force=.FALSE.)
292 :
293 50 : trq_plus(p_iter) = dft_control%trq(1) + dft_control%trq(2)
294 50 : trq_minus(p_iter) = dft_control%trq(1) - dft_control%trq(2)
295 50 : vhxc_plus(p_iter) = dft_control%vhxc(1) + dft_control%vhxc(2)
296 60 : 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 60 : DO k = 1, n
302 50 : sum_trq_plus_x_vhxc_plus = sum_trq_plus_x_vhxc_plus + trq_plus(k)*vhxc_plus(k)
303 50 : sum_trq_plus = sum_trq_plus + trq_plus(k)
304 50 : sum_vhxc_plus = sum_vhxc_plus + vhxc_plus(k)
305 60 : sum_trq_plus_x_trq_plus = sum_trq_plus_x_trq_plus + trq_plus(k)*trq_plus(k)
306 : END DO
307 :
308 60 : DO k = 1, n
309 50 : sum_trq_minus_x_vhxc_minus = sum_trq_minus_x_vhxc_minus + trq_minus(k)*vhxc_minus(k)
310 50 : sum_trq_minus = sum_trq_minus + trq_minus(k)
311 50 : sum_vhxc_minus = sum_vhxc_minus + vhxc_minus(k)
312 60 : sum_trq_minus_x_trq_minus = sum_trq_minus_x_trq_minus + trq_minus(k)*trq_minus(k)
313 : END DO
314 :
315 10 : denominator_plus = l*sum_trq_plus_x_trq_plus - sum_trq_plus**2.0_dp
316 10 : denominator_minus = l*sum_trq_minus_x_trq_minus - sum_trq_minus**2.0_dp
317 10 : IF (ABS(denominator_plus) < 100.0_dp*EPSILON(1.0_dp)) THEN
318 0 : CPABORT("MTLR regression is singular.")
319 : END IF
320 10 : IF (ABS(denominator_minus) < 100.0_dp*EPSILON(1.0_dp)) THEN
321 0 : 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 10 : /denominator_minus
325 : fhxc_plus = (l*sum_trq_plus_x_vhxc_plus - sum_trq_plus*sum_vhxc_plus) &
326 10 : /denominator_plus
327 10 : intercept_minus = (sum_vhxc_minus - fhxc_minus*sum_trq_minus)/l
328 10 : intercept_plus = (sum_vhxc_plus - fhxc_plus*sum_trq_plus)/l
329 60 : DO k = 1, n
330 60 : std_minus = std_minus + (vhxc_minus(k) - (trq_minus(k)*fhxc_minus + intercept_minus))**2/(l - 2)
331 : END DO
332 60 : DO k = 1, n
333 60 : std_plus = std_plus + (vhxc_plus(k) - (trq_plus(k)*fhxc_plus + intercept_plus))**2/(l - 2)
334 : END DO
335 10 : std_minus = SQRT(std_minus)
336 10 : std_plus = SQRT(std_plus)
337 :
338 10 : u_new(ikind) = 0.5_dp*fhxc_plus
339 10 : j_new(ikind) = -0.5_dp*fhxc_minus
340 10 : delta_u = u_new(ikind) - u_old(ikind)
341 10 : delta_j = j_new(ikind) - j_old(ikind)
342 10 : std_u = 0.5_dp*std_plus
343 10 : std_j = 0.5_dp*std_minus
344 :
345 10 : IF (output_unit > 0) THEN
346 5 : WRITE (UNIT=output_unit, FMT="(/,T2,78('='))")
347 : WRITE (UNIT=output_unit, FMT="(T2,A,I0,A,I0,A)") &
348 5 : "MTLR| U/J iteration ", u_iter, &
349 10 : " results for KIND ", ikind, "."
350 5 : WRITE (UNIT=output_unit, FMT="(T2,78('-'))")
351 : WRITE (UNIT=output_unit, FMT="(T2,A)") &
352 5 : "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 5 : "Pt", &
356 5 : "Pert.[eV]", &
357 5 : "trq(+)", &
358 5 : "vhxc(+)[eV]", &
359 5 : "vhxc(+)-fit", &
360 10 : "vhxc(+)-resid"
361 30 : 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 25 : k, &
366 25 : perturbation_strength(k)*evolt, &
367 25 : trq_plus(k), &
368 25 : vhxc_plus(k)*evolt, &
369 25 : (fhxc_plus*trq_plus(k) + intercept_plus)*evolt, &
370 : (vhxc_plus(k) - &
371 55 : (fhxc_plus*trq_plus(k) + intercept_plus))*evolt
372 : END DO
373 5 : WRITE (UNIT=output_unit, FMT="(T2,78('-'))")
374 : WRITE (UNIT=output_unit, FMT="(T2,A)") &
375 5 : "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 5 : "Pt", &
379 5 : "Pert.[eV]", &
380 5 : "trq(-)", &
381 5 : "vhxc(-)[eV]", &
382 5 : "vhxc(-)-fit", &
383 10 : "vhxc(-)-resid"
384 30 : 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 25 : k, &
389 25 : perturbation_strength(k)*evolt, &
390 25 : trq_minus(k), &
391 25 : vhxc_minus(k)*evolt, &
392 25 : (fhxc_minus*trq_minus(k) + intercept_minus)*evolt, &
393 : (vhxc_minus(k) - &
394 55 : (fhxc_minus*trq_minus(k) + intercept_minus))*evolt
395 : END DO
396 5 : WRITE (UNIT=output_unit, FMT="(T2,78('-'))")
397 : WRITE (UNIT=output_unit, &
398 : FMT="(T2,A,T22,A16,T43,A16,T64,A16)") &
399 5 : "MTLR| Parameter", &
400 5 : "Old [eV]", &
401 5 : "New [eV]", &
402 10 : "Change [eV]"
403 5 : 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 5 : "MTLR| Hubbard U", &
407 5 : u_old(ikind)*evolt, &
408 5 : u_new(ikind)*evolt, &
409 10 : (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 5 : "MTLR| Hund J", &
413 5 : j_old(ikind)*evolt, &
414 5 : j_new(ikind)*evolt, &
415 10 : (j_new(ikind) - j_old(ikind))*evolt
416 5 : WRITE (UNIT=output_unit, FMT="(T2,78('-'))")
417 : WRITE (UNIT=output_unit, FMT="(T2,A,T61,ES16.8,A3)") &
418 5 : "MTLR| Absolute change in U:", &
419 10 : ABS(u_new(ikind) - u_old(ikind))*evolt, " eV"
420 : WRITE (UNIT=output_unit, FMT="(T2,A,T61,ES16.8,A3)") &
421 5 : "MTLR| Absolute change in J:", &
422 10 : ABS(j_new(ikind) - j_old(ikind))*evolt, " eV"
423 : WRITE (UNIT=output_unit, FMT="(T2,A,T61,ES16.8,A3)") &
424 5 : "MTLR| U fit residual:", &
425 10 : std_u*evolt, " eV"
426 : WRITE (UNIT=output_unit, FMT="(T2,A,T61,ES16.8,A3)") &
427 5 : "MTLR| J fit residual:", &
428 10 : std_j*evolt, " eV"
429 5 : WRITE (UNIT=output_unit, FMT="(T2,78('='))")
430 : END IF
431 :
432 10 : DEALLOCATE (perturbation_strength)
433 10 : DEALLOCATE (trq_plus)
434 10 : DEALLOCATE (vhxc_plus)
435 10 : DEALLOCATE (trq_minus)
436 20 : DEALLOCATE (vhxc_minus)
437 :
438 : END DO
439 :
440 20 : DO ikind = 1, nkind
441 10 : IF (.NOT. mtlr_kind(ikind)) CYCLE
442 10 : qs_kind_set(ikind)%dft_plus_u%u_minus_j = u_new(ikind) - j_new(ikind)
443 20 : qs_kind_set(ikind)%dft_plus_u%hund_j = j_new(ikind)
444 : END DO
445 :
446 30 : max_delta_u = MAXVAL(ABS(PACK(u_new(:) - u_old(:), mtlr_kind)))
447 30 : max_delta_j = MAXVAL(ABS(PACK(j_new(:) - j_old(:), mtlr_kind)))
448 10 : IF (max_delta_u < eps_u_j_loop .AND. max_delta_j < eps_u_j_loop) THEN
449 4 : converged = .TRUE.
450 : END IF
451 :
452 10 : IF (output_unit > 0) THEN
453 5 : WRITE (UNIT=output_unit, FMT="(/,T2,78('='))")
454 :
455 : WRITE (UNIT=output_unit, FMT="(T2,A,I0,A)") &
456 5 : "MTLR| Iteration ", u_iter, &
457 10 : " completed for all atomic kinds."
458 :
459 5 : WRITE (UNIT=output_unit, FMT="(T2,78('-'))")
460 :
461 : WRITE (UNIT=output_unit, FMT="(T2,A,ES16.8,A)") &
462 5 : "MTLR| Maximum change in U: ", &
463 10 : max_delta_u*evolt, " eV"
464 :
465 : WRITE (UNIT=output_unit, FMT="(T2,A,ES16.8,A)") &
466 5 : "MTLR| Maximum change in J: ", &
467 10 : max_delta_j*evolt, " eV"
468 :
469 : WRITE (UNIT=output_unit, FMT="(T2,A,ES16.8,A)") &
470 5 : "MTLR| Convergence threshold: ", &
471 10 : eps_u_j_loop*evolt, " eV"
472 :
473 5 : WRITE (UNIT=output_unit, FMT="(T2,78('-'))")
474 :
475 5 : IF (converged) THEN
476 : WRITE (UNIT=output_unit, FMT="(T2,A)") &
477 2 : "MTLR| U and J have converged."
478 :
479 : WRITE (UNIT=output_unit, FMT="(T2,A)") &
480 2 : "MTLR| Both maximum changes are below the convergence threshold."
481 :
482 3 : ELSE IF (u_iter < max_mtlr_iter) THEN
483 : WRITE (UNIT=output_unit, FMT="(T2,A)") &
484 3 : "MTLR| U and J have not yet converged."
485 :
486 : WRITE (UNIT=output_unit, FMT="(T2,A)") &
487 3 : "MTLR| Proceeding to the next linear-response U/J iteration."
488 :
489 : ELSE
490 : WRITE (UNIT=output_unit, FMT="(T2,A)") &
491 0 : "MTLR| U and J have not converged."
492 :
493 : WRITE (UNIT=output_unit, FMT="(T2,A)") &
494 0 : "MTLR| The maximum number of MTLR iterations has been reached."
495 : END IF
496 :
497 5 : WRITE (UNIT=output_unit, FMT="(T2,78('='))")
498 : END IF
499 :
500 10 : IF (converged) THEN
501 4 : IF (output_unit > 0) THEN
502 2 : WRITE (UNIT=output_unit, FMT="(/,T2,78('*'))")
503 : WRITE (UNIT=output_unit, FMT="(T2,A,I0,A)") &
504 2 : "MTLR| U and J converged after ", &
505 4 : u_iter, " iterations."
506 2 : WRITE (UNIT=output_unit, FMT="(T2,78('*'),/)")
507 : END IF
508 6 : ELSE IF (u_iter == max_mtlr_iter) THEN
509 0 : IF (output_unit > 0) THEN
510 0 : WRITE (UNIT=output_unit, FMT="(/,T2,78('*'))")
511 : WRITE (UNIT=output_unit, FMT="(T2,A,I0,A)") &
512 0 : "MTLR| U and J did not converge within the maximum of ", &
513 0 : max_mtlr_iter, " iterations."
514 : WRITE (UNIT=output_unit, FMT="(T2,A)") &
515 0 : "MTLR| Results from the final iteration will be reported."
516 0 : WRITE (UNIT=output_unit, FMT="(T2,78('*'),/)")
517 : END IF
518 : END IF
519 :
520 10 : IF (converged .OR. u_iter == max_mtlr_iter) THEN
521 4 : IF (output_unit > 0) THEN
522 2 : WRITE (UNIT=output_unit, FMT="(/,T2,78('*'))")
523 : WRITE (UNIT=output_unit, FMT="(T2,A,ES16.8,A)") &
524 2 : "MTLR| Maximum change in U: ", &
525 4 : max_delta_u*evolt, " eV"
526 : WRITE (UNIT=output_unit, FMT="(T2,A,ES16.8,A)") &
527 2 : "MTLR| Maximum change in J: ", &
528 4 : max_delta_j*evolt, " eV"
529 : WRITE (UNIT=output_unit, FMT="(T2,A,ES16.8,A)") &
530 2 : "MTLR| Maximum change over U and J: ", &
531 4 : MAX(max_delta_u, max_delta_j)*evolt, " eV"
532 2 : WRITE (UNIT=output_unit, FMT="(T2,78('*'))")
533 4 : DO ikind = 1, nkind
534 2 : IF (.NOT. mtlr_kind(ikind)) CYCLE
535 : WRITE (UNIT=output_unit, FMT="(/,T2,A,I0,A)") &
536 2 : "MTLR| Final parameters for KIND ", ikind, ":"
537 : WRITE (UNIT=output_unit, FMT="(T2,A,F16.8,A)") &
538 2 : "MTLR| Calculated Hubbard U: ", &
539 4 : u_new(ikind)*evolt, " eV"
540 : WRITE (UNIT=output_unit, FMT="(T2,A,F16.8,A)") &
541 2 : "MTLR| Calculated Hund J: ", &
542 4 : j_new(ikind)*evolt, " eV"
543 : WRITE (UNIT=output_unit, FMT="(T2,A)") &
544 2 : "MTLR| Recommended CP2K input parameters:"
545 : WRITE (UNIT=output_unit, FMT="(T2,A,F16.8)") &
546 2 : "MTLR| U_MINUS_J [eV] ", &
547 4 : (u_new(ikind) - j_new(ikind))*evolt
548 : WRITE (UNIT=output_unit, FMT="(T2,A,F16.8)") &
549 2 : "MTLR| J [eV] ", &
550 6 : j_new(ikind)*evolt
551 : END DO
552 2 : WRITE (UNIT=output_unit, FMT="(/,T2,78('*'),/)")
553 : END IF
554 : END IF
555 :
556 10 : IF (converged) THEN
557 : EXIT
558 : END IF
559 :
560 12 : u_old(:) = u_new
561 16 : j_old(:) = j_new
562 :
563 : END DO
564 :
565 8 : DO ikind = 1, nkind
566 4 : IF (.NOT. mtlr_kind(ikind)) CYCLE
567 4 : qs_kind_set(ikind)%dft_plus_u%u_minus_j = u_new(ikind) - j_new(ikind)
568 8 : qs_kind_set(ikind)%dft_plus_u%hund_j = j_new(ikind)
569 : END DO
570 4 : dft_control%mtlr_dft_with_perturbation = .FALSE.
571 4 : dft_control%perturbation_strength = 0.0_dp
572 :
573 4 : CALL force_env_calc_energy_force(force_env, calc_force=.FALSE.)
574 :
575 4 : DEALLOCATE (u_new)
576 4 : DEALLOCATE (j_new)
577 4 : DEALLOCATE (u_old)
578 4 : DEALLOCATE (j_old)
579 4 : DEALLOCATE (mtlr_kind)
580 :
581 4 : CALL timestop(handle)
582 :
583 12 : END SUBROUTINE do_mtlr_u_j
584 :
585 : END MODULE mtlr_u_j_methods
|
