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_dbcsr_operations, ONLY: copy_fm_to_dbcsr
19 : USE cp_log_handling, ONLY: cp_get_default_logger,&
20 : cp_logger_get_default_io_unit,&
21 : cp_logger_type
22 : USE force_env_methods, ONLY: force_env_calc_energy_force
23 : USE force_env_types, ONLY: force_env_type
24 : USE input_constants, ONLY: mtlr_atomic_perturbations,&
25 : mtlr_reference_from_atomic,&
26 : mtlr_reference_from_restart
27 : USE kinds, ONLY: dp
28 : USE physcon, ONLY: evolt
29 : USE qs_environment_types, ONLY: get_qs_env
30 : USE qs_kind_types, ONLY: qs_kind_type
31 : USE qs_mo_types, ONLY: deallocate_mo_set,&
32 : duplicate_mo_set,&
33 : mo_set_type,&
34 : reassign_allocated_mos
35 : #include "./base/base_uses.f90"
36 :
37 : IMPLICIT NONE
38 :
39 : PRIVATE
40 : PUBLIC :: do_mtlr_u_j
41 :
42 : CONTAINS
43 : ! **************************************************************************************************
44 : !> \brief Driver for self-consistent MTLR U/J iteration.
45 : !> Each outer iteration performs:
46 : !> 1) one standard ENERGY SCF
47 : !> 2) one MTLR evaluation of U and J
48 : !> 3) one update of the Hubbard parameters
49 : !> until U and J are converged.
50 : !> using a method based on Lowdin charges
51 : !> \f[Q = S^{1/2} P S^{1/2}\f]
52 : !> where \b P and \b S are the density and the
53 : !> overlap matrix, respectively.
54 : !> \param[in,out] force_env ...
55 : !> \date 29.07.2026
56 : !> \author Ziwei Chai
57 : !> \version 1.0
58 : ! **************************************************************************************************
59 16 : SUBROUTINE do_mtlr_u_j(force_env)
60 :
61 : TYPE(force_env_type), INTENT(INOUT), POINTER :: force_env
62 :
63 : CHARACTER(LEN=*), PARAMETER :: routineN = 'do_mtlr_u_j'
64 :
65 : INTEGER :: handle, ikind, iset, k, max_mtlr_iter, &
66 : n, nkind, output_unit, p_iter, u_iter
67 16 : INTEGER, DIMENSION(:), POINTER :: atom_list
68 : LOGICAL :: any_dft_plus_u, any_mtlr_kind, &
69 : converged, do_reference_scf
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 16 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: j_new, j_old, perturbation_strength, &
77 16 : trq_minus, trq_plus, u_new, u_old, &
78 16 : vhxc_minus, vhxc_plus
79 16 : 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 16 : TYPE(mo_set_type), ALLOCATABLE, DIMENSION(:) :: reference_mos
83 16 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
84 16 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
85 :
86 16 : CALL timeset(routineN, handle)
87 :
88 16 : NULLIFY (atom_list, qs_kind_set, dft_control, logger, mos)
89 :
90 16 : logger => cp_get_default_logger()
91 16 : output_unit = cp_logger_get_default_io_unit(logger)
92 :
93 16 : CPASSERT(ASSOCIATED(force_env))
94 16 : 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 : mos=mos, &
100 16 : atomic_kind_set=atomic_kind_set)
101 :
102 16 : CPASSERT(ASSOCIATED(atomic_kind_set))
103 16 : CPASSERT(ASSOCIATED(dft_control))
104 16 : CPASSERT(ASSOCIATED(mos))
105 16 : CPASSERT(ASSOCIATED(qs_kind_set))
106 :
107 16 : nkind = SIZE(atomic_kind_set)
108 16 : IF (SIZE(qs_kind_set) /= nkind) THEN
109 0 : CPABORT("The atomic-kind and Quickstep-kind arrays have inconsistent sizes.")
110 : END IF
111 48 : ALLOCATE (mtlr_kind(nkind))
112 16 : mtlr_kind(:) = .FALSE.
113 16 : any_dft_plus_u = .FALSE.
114 16 : any_mtlr_kind = .FALSE.
115 36 : DO ikind = 1, nkind
116 20 : IF (.NOT. ASSOCIATED(qs_kind_set(ikind)%dft_plus_u)) CYCLE
117 16 : any_dft_plus_u = .TRUE.
118 32 : IF (qs_kind_set(ikind)%dft_plus_u%do_mtlr) THEN
119 16 : mtlr_kind(ikind) = .TRUE.
120 16 : any_mtlr_kind = .TRUE.
121 : END IF
122 : END DO
123 16 : 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 16 : 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 16 : SELECT CASE (dft_control%mtlr_initialization_mode)
133 : CASE (mtlr_atomic_perturbations)
134 12 : do_reference_scf = .FALSE.
135 : CASE (mtlr_reference_from_atomic, mtlr_reference_from_restart)
136 12 : do_reference_scf = .TRUE.
137 : CASE DEFAULT
138 16 : CPABORT("The MTLR SCF initialization mode was not resolved.")
139 : END SELECT
140 :
141 16 : IF (output_unit > 0) THEN
142 8 : WRITE (UNIT=output_unit, FMT="(/,T2,78('='))")
143 : WRITE (UNIT=output_unit, FMT="(T2,A)") &
144 8 : "MTLR| SCF initialization"
145 8 : WRITE (UNIT=output_unit, FMT="(T2,78('-'))")
146 10 : SELECT CASE (dft_control%mtlr_initialization_mode)
147 : CASE (mtlr_atomic_perturbations)
148 : WRITE (UNIT=output_unit, FMT="(T2,A,T66,A13)") &
149 2 : "MTLR| Reference SCF:", "OFF"
150 : WRITE (UNIT=output_unit, FMT="(T2,A,T66,A13)") &
151 2 : "MTLR| Perturbation initial guess:", "ATOMIC"
152 : CASE (mtlr_reference_from_atomic)
153 : WRITE (UNIT=output_unit, FMT="(T2,A,T66,A13)") &
154 4 : "MTLR| Reference SCF:", "ON"
155 : WRITE (UNIT=output_unit, FMT="(T2,A,T66,A13)") &
156 4 : "MTLR| Reference SCF initial guess:", "ATOMIC"
157 : WRITE (UNIT=output_unit, FMT="(T2,A,T66,A13)") &
158 4 : "MTLR| Perturbation initial guess:", "REFERENCE MOs"
159 : CASE (mtlr_reference_from_restart)
160 : WRITE (UNIT=output_unit, FMT="(T2,A,T66,A13)") &
161 2 : "MTLR| Reference SCF:", "ON"
162 : WRITE (UNIT=output_unit, FMT="(T2,A,T66,A13)") &
163 2 : "MTLR| Reference SCF initial guess:", "RESTART"
164 : WRITE (UNIT=output_unit, FMT="(T2,A,T66,A13)") &
165 10 : "MTLR| Perturbation initial guess:", "REFERENCE MOs"
166 : END SELECT
167 8 : WRITE (UNIT=output_unit, FMT="(T2,78('='))")
168 : END IF
169 :
170 48 : ALLOCATE (u_new(nkind))
171 32 : ALLOCATE (j_new(nkind))
172 32 : ALLOCATE (u_old(nkind))
173 32 : ALLOCATE (j_old(nkind))
174 16 : u_new(:) = 0.0_dp
175 16 : j_new(:) = 0.0_dp
176 16 : u_old(:) = 0.0_dp
177 16 : j_old(:) = 0.0_dp
178 16 : converged = .FALSE.
179 16 : eps_u_j_loop = dft_control%eps_u_j_loop
180 16 : max_mtlr_iter = dft_control%max_mtlr_iter
181 16 : IF (dft_control%nspins /= 2) THEN
182 0 : CPABORT("Unrestricted KS has to be used (the number of spin channels should be 2).")
183 : END IF
184 16 : IF (max_mtlr_iter < 1) THEN
185 0 : CPABORT("MAX_MTLR_LOOP must be at least one.")
186 : END IF
187 16 : IF (eps_u_j_loop <= 0.0_dp) THEN
188 0 : CPABORT("EPS_U_J_LOOP must be positive.")
189 : END IF
190 :
191 : ! Ensure that the DFT+U+J machinery remains active during the
192 : ! initial MTLR calculation, even when a parameter starts from zero.
193 36 : DO ikind = 1, nkind
194 20 : IF (.NOT. mtlr_kind(ikind)) CYCLE
195 16 : IF (qs_kind_set(ikind)%dft_plus_u%u_minus_j == 0.0_dp) THEN
196 0 : qs_kind_set(ikind)%dft_plus_u%u_minus_j = EPSILON(1.0_dp)
197 : END IF
198 16 : IF (qs_kind_set(ikind)%dft_plus_u%hund_j == 0.0_dp) THEN
199 0 : qs_kind_set(ikind)%dft_plus_u%hund_j = EPSILON(1.0_dp)
200 : END IF
201 16 : j_old(ikind) = qs_kind_set(ikind)%dft_plus_u%hund_j
202 : u_old(ikind) = qs_kind_set(ikind)%dft_plus_u%u_minus_j + &
203 36 : qs_kind_set(ikind)%dft_plus_u%hund_j
204 : END DO
205 :
206 44 : DO u_iter = 1, max_mtlr_iter
207 :
208 44 : dft_control%mtlr_dft_with_perturbation = .FALSE.
209 :
210 44 : IF (do_reference_scf) THEN
211 32 : IF (output_unit > 0) THEN
212 16 : WRITE (UNIT=output_unit, FMT="(/,T2,78('='))")
213 : WRITE (UNIT=output_unit, FMT="(T2,A)") &
214 16 : "MTLR| Starting the unperturbed reference SCF."
215 : WRITE (UNIT=output_unit, FMT="(T2,A,T72,I8)") &
216 16 : "MTLR| U/J iteration:", u_iter
217 16 : WRITE (UNIT=output_unit, FMT="(T2,78('='))")
218 : END IF
219 32 : CALL force_env_calc_energy_force(force_env, calc_force=.FALSE.)
220 :
221 32 : IF (.NOT. ALLOCATED(reference_mos)) THEN
222 60 : ALLOCATE (reference_mos(SIZE(mos)))
223 : ELSE
224 60 : DO iset = 1, SIZE(reference_mos)
225 60 : CALL deallocate_mo_set(reference_mos(iset))
226 : END DO
227 : END IF
228 96 : DO iset = 1, SIZE(mos)
229 96 : CALL duplicate_mo_set(reference_mos(iset), mos(iset))
230 : END DO
231 : END IF
232 :
233 100 : DO ikind = 1, nkind
234 :
235 56 : IF (.NOT. mtlr_kind(ikind)) CYCLE
236 44 : CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom_list)
237 44 : IF (.NOT. ANY(atom_list == qs_kind_set(ikind)%dft_plus_u%lr_atom)) THEN
238 0 : CPABORT("INDEX_PERTURBED_ATOM does not belong to the KIND containing the MTLR section.")
239 : END IF
240 44 : IF (.NOT. ALLOCATED( &
241 : qs_kind_set(ikind)%dft_plus_u%perturbation_strength)) THEN
242 0 : CPABORT("MTLR target does not contain perturbation strengths.")
243 : END IF
244 :
245 44 : dft_control%mtlr_ikind = ikind
246 :
247 44 : n = SIZE(qs_kind_set(ikind)%dft_plus_u%perturbation_strength)
248 44 : IF (n < 3) THEN
249 0 : CPABORT("MTLR linear regression requires at least three perturbation strengths.")
250 : END IF
251 44 : l = REAL(n, dp)
252 132 : ALLOCATE (perturbation_strength(n))
253 88 : ALLOCATE (trq_plus(n))
254 88 : ALLOCATE (vhxc_plus(n))
255 88 : ALLOCATE (trq_minus(n))
256 88 : ALLOCATE (vhxc_minus(n))
257 44 : trq_plus = 0.0_dp
258 44 : trq_minus = 0.0_dp
259 44 : vhxc_plus = 0.0_dp
260 44 : vhxc_minus = 0.0_dp
261 44 : sum_trq_plus = 0.0_dp
262 44 : sum_vhxc_plus = 0.0_dp
263 44 : sum_trq_plus_x_trq_plus = 0.0_dp
264 44 : sum_trq_plus_x_vhxc_plus = 0.0_dp
265 44 : sum_trq_minus = 0.0_dp
266 44 : sum_vhxc_minus = 0.0_dp
267 44 : sum_trq_minus_x_trq_minus = 0.0_dp
268 44 : sum_trq_minus_x_vhxc_minus = 0.0_dp
269 44 : std_plus = 0.0_dp
270 44 : std_minus = 0.0_dp
271 264 : perturbation_strength(:) = qs_kind_set(ikind)%dft_plus_u%perturbation_strength(:)
272 :
273 264 : DO p_iter = 1, n
274 :
275 220 : IF (output_unit > 0) THEN
276 110 : WRITE (UNIT=output_unit, FMT="(/,T2,78('='))")
277 : WRITE (UNIT=output_unit, FMT="(T2,A,T72,I8)") &
278 110 : "MTLR| U/J iteration:", u_iter
279 : WRITE (UNIT=output_unit, FMT="(T2,A,T74,I1,A4,I1)") &
280 110 : "MTLR| Perturbation SCF:", p_iter, " of ", n
281 : WRITE (UNIT=output_unit, FMT="(T2,A,T72,I8)") &
282 110 : "MTLR| Target KIND index:", ikind
283 : WRITE (UNIT=output_unit, FMT="(T2,A,T72,I8)") &
284 110 : "MTLR| Target atom index:", &
285 220 : qs_kind_set(ikind)%dft_plus_u%lr_atom
286 : WRITE (UNIT=output_unit, FMT="(T2,A,T63,F14.8,A3)") &
287 110 : "MTLR| Perturbation strength:", &
288 220 : perturbation_strength(p_iter)*evolt, " eV"
289 110 : WRITE (UNIT=output_unit, FMT="(T2,78('='))")
290 : END IF
291 :
292 220 : dft_control%perturbation_strength = perturbation_strength(p_iter)
293 220 : dft_control%mtlr_dft_with_perturbation = .TRUE.
294 :
295 220 : IF (do_reference_scf) CALL restore_reference_mos(mos, reference_mos)
296 220 : CALL force_env_calc_energy_force(force_env, calc_force=.FALSE.)
297 :
298 220 : trq_plus(p_iter) = dft_control%trq(1) + dft_control%trq(2)
299 220 : trq_minus(p_iter) = dft_control%trq(1) - dft_control%trq(2)
300 220 : vhxc_plus(p_iter) = dft_control%vhxc(1) + dft_control%vhxc(2)
301 264 : vhxc_minus(p_iter) = dft_control%vhxc(1) - dft_control%vhxc(2)
302 :
303 : END DO
304 :
305 : !Calculate all the number about vhxc and trq.
306 264 : DO k = 1, n
307 220 : sum_trq_plus_x_vhxc_plus = sum_trq_plus_x_vhxc_plus + trq_plus(k)*vhxc_plus(k)
308 220 : sum_trq_plus = sum_trq_plus + trq_plus(k)
309 220 : sum_vhxc_plus = sum_vhxc_plus + vhxc_plus(k)
310 264 : sum_trq_plus_x_trq_plus = sum_trq_plus_x_trq_plus + trq_plus(k)*trq_plus(k)
311 : END DO
312 :
313 264 : DO k = 1, n
314 220 : sum_trq_minus_x_vhxc_minus = sum_trq_minus_x_vhxc_minus + trq_minus(k)*vhxc_minus(k)
315 220 : sum_trq_minus = sum_trq_minus + trq_minus(k)
316 220 : sum_vhxc_minus = sum_vhxc_minus + vhxc_minus(k)
317 264 : sum_trq_minus_x_trq_minus = sum_trq_minus_x_trq_minus + trq_minus(k)*trq_minus(k)
318 : END DO
319 :
320 44 : denominator_plus = l*sum_trq_plus_x_trq_plus - sum_trq_plus**2.0_dp
321 44 : denominator_minus = l*sum_trq_minus_x_trq_minus - sum_trq_minus**2.0_dp
322 44 : IF (ABS(denominator_plus) < 100.0_dp*EPSILON(1.0_dp)) THEN
323 0 : CPABORT("MTLR regression is singular.")
324 : END IF
325 44 : IF (ABS(denominator_minus) < 100.0_dp*EPSILON(1.0_dp)) THEN
326 0 : CPABORT("MTLR regression is singular.")
327 : END IF
328 : fhxc_minus = (l*sum_trq_minus_x_vhxc_minus - sum_trq_minus*sum_vhxc_minus) &
329 44 : /denominator_minus
330 : fhxc_plus = (l*sum_trq_plus_x_vhxc_plus - sum_trq_plus*sum_vhxc_plus) &
331 44 : /denominator_plus
332 44 : intercept_minus = (sum_vhxc_minus - fhxc_minus*sum_trq_minus)/l
333 44 : intercept_plus = (sum_vhxc_plus - fhxc_plus*sum_trq_plus)/l
334 264 : DO k = 1, n
335 264 : std_minus = std_minus + (vhxc_minus(k) - (trq_minus(k)*fhxc_minus + intercept_minus))**2/(l - 2)
336 : END DO
337 264 : DO k = 1, n
338 264 : std_plus = std_plus + (vhxc_plus(k) - (trq_plus(k)*fhxc_plus + intercept_plus))**2/(l - 2)
339 : END DO
340 44 : std_minus = SQRT(std_minus)
341 44 : std_plus = SQRT(std_plus)
342 :
343 44 : u_new(ikind) = 0.5_dp*fhxc_plus
344 44 : j_new(ikind) = -0.5_dp*fhxc_minus
345 44 : delta_u = u_new(ikind) - u_old(ikind)
346 44 : delta_j = j_new(ikind) - j_old(ikind)
347 44 : std_u = 0.5_dp*std_plus
348 44 : std_j = 0.5_dp*std_minus
349 :
350 44 : IF (output_unit > 0) THEN
351 22 : WRITE (UNIT=output_unit, FMT="(/,T2,78('='))")
352 : WRITE (UNIT=output_unit, FMT="(T2,A,I0,A,I0,A)") &
353 22 : "MTLR| U/J iteration ", u_iter, &
354 44 : " results for KIND ", ikind, "."
355 22 : WRITE (UNIT=output_unit, FMT="(T2,78('-'))")
356 : WRITE (UNIT=output_unit, FMT="(T2,A)") &
357 22 : "MTLR| Linear-regression data for Hubbard U"
358 : WRITE (UNIT=output_unit, &
359 : FMT="(T2,A3,T8,A10,T22,A12,T37,A13,T52,A13,T67,A13)") &
360 22 : "Pt", &
361 22 : "Pert.[eV]", &
362 22 : "trq(+)", &
363 22 : "vhxc(+)[eV]", &
364 22 : "vhxc(+)-fit", &
365 44 : "vhxc(+)-resid"
366 132 : DO k = 1, n
367 : WRITE (UNIT=output_unit, &
368 : FMT="(T2,I3,T8,F10.4,T22,ES12.4,T37,ES13.5,"// &
369 : " T52,ES13.5,T67,ES13.4)") &
370 110 : k, &
371 110 : perturbation_strength(k)*evolt, &
372 110 : trq_plus(k), &
373 110 : vhxc_plus(k)*evolt, &
374 110 : (fhxc_plus*trq_plus(k) + intercept_plus)*evolt, &
375 : (vhxc_plus(k) - &
376 242 : (fhxc_plus*trq_plus(k) + intercept_plus))*evolt
377 : END DO
378 22 : WRITE (UNIT=output_unit, FMT="(T2,78('-'))")
379 : WRITE (UNIT=output_unit, FMT="(T2,A)") &
380 22 : "MTLR| Linear-regression data for Hund J"
381 : WRITE (UNIT=output_unit, &
382 : FMT="(T2,A3,T8,A10,T22,A12,T37,A13,T52,A13,T67,A13)") &
383 22 : "Pt", &
384 22 : "Pert.[eV]", &
385 22 : "trq(-)", &
386 22 : "vhxc(-)[eV]", &
387 22 : "vhxc(-)-fit", &
388 44 : "vhxc(-)-resid"
389 132 : DO k = 1, n
390 : WRITE (UNIT=output_unit, &
391 : FMT="(T2,I3,T8,F10.4,T22,ES12.4,T37,ES13.5,"// &
392 : " T52,ES13.5,T67,ES13.4)") &
393 110 : k, &
394 110 : perturbation_strength(k)*evolt, &
395 110 : trq_minus(k), &
396 110 : vhxc_minus(k)*evolt, &
397 110 : (fhxc_minus*trq_minus(k) + intercept_minus)*evolt, &
398 : (vhxc_minus(k) - &
399 242 : (fhxc_minus*trq_minus(k) + intercept_minus))*evolt
400 : END DO
401 22 : WRITE (UNIT=output_unit, FMT="(T2,78('-'))")
402 : WRITE (UNIT=output_unit, &
403 : FMT="(T2,A,T22,A16,T43,A16,T64,A16)") &
404 22 : "MTLR| Parameter", &
405 22 : "Old [eV]", &
406 22 : "New [eV]", &
407 44 : "Change [eV]"
408 22 : WRITE (UNIT=output_unit, FMT="(T2,78('-'))")
409 : WRITE (UNIT=output_unit, &
410 : FMT="(T2,A,T22,ES16.8,T43,ES16.8,T64,ES16.8)") &
411 22 : "MTLR| Hubbard U", &
412 22 : u_old(ikind)*evolt, &
413 22 : u_new(ikind)*evolt, &
414 44 : (u_new(ikind) - u_old(ikind))*evolt
415 : WRITE (UNIT=output_unit, &
416 : FMT="(T2,A,T22,ES16.8,T43,ES16.8,T64,ES16.8)") &
417 22 : "MTLR| Hund J", &
418 22 : j_old(ikind)*evolt, &
419 22 : j_new(ikind)*evolt, &
420 44 : (j_new(ikind) - j_old(ikind))*evolt
421 22 : WRITE (UNIT=output_unit, FMT="(T2,78('-'))")
422 : WRITE (UNIT=output_unit, FMT="(T2,A,T61,ES16.8,A3)") &
423 22 : "MTLR| Absolute change in U:", &
424 44 : ABS(u_new(ikind) - u_old(ikind))*evolt, " eV"
425 : WRITE (UNIT=output_unit, FMT="(T2,A,T61,ES16.8,A3)") &
426 22 : "MTLR| Absolute change in J:", &
427 44 : ABS(j_new(ikind) - j_old(ikind))*evolt, " eV"
428 : WRITE (UNIT=output_unit, FMT="(T2,A,T61,ES16.8,A3)") &
429 22 : "MTLR| U fit residual:", &
430 44 : std_u*evolt, " eV"
431 : WRITE (UNIT=output_unit, FMT="(T2,A,T61,ES16.8,A3)") &
432 22 : "MTLR| J fit residual:", &
433 44 : std_j*evolt, " eV"
434 22 : WRITE (UNIT=output_unit, FMT="(T2,78('='))")
435 : END IF
436 :
437 44 : DEALLOCATE (perturbation_strength)
438 44 : DEALLOCATE (trq_plus)
439 44 : DEALLOCATE (vhxc_plus)
440 44 : DEALLOCATE (trq_minus)
441 100 : DEALLOCATE (vhxc_minus)
442 :
443 : END DO
444 :
445 44 : IF (do_reference_scf) CALL restore_reference_mos(mos, reference_mos)
446 :
447 100 : DO ikind = 1, nkind
448 56 : IF (.NOT. mtlr_kind(ikind)) CYCLE
449 44 : qs_kind_set(ikind)%dft_plus_u%u_minus_j = u_new(ikind) - j_new(ikind)
450 100 : qs_kind_set(ikind)%dft_plus_u%hund_j = j_new(ikind)
451 : END DO
452 :
453 144 : max_delta_u = MAXVAL(ABS(PACK(u_new(:) - u_old(:), mtlr_kind)))
454 144 : max_delta_j = MAXVAL(ABS(PACK(j_new(:) - j_old(:), mtlr_kind)))
455 44 : IF (max_delta_u < eps_u_j_loop .AND. max_delta_j < eps_u_j_loop) THEN
456 16 : converged = .TRUE.
457 : END IF
458 :
459 44 : IF (output_unit > 0) THEN
460 22 : WRITE (UNIT=output_unit, FMT="(/,T2,78('='))")
461 :
462 : WRITE (UNIT=output_unit, FMT="(T2,A,I0,A)") &
463 22 : "MTLR| Iteration ", u_iter, &
464 44 : " completed for all atomic kinds."
465 :
466 22 : WRITE (UNIT=output_unit, FMT="(T2,78('-'))")
467 :
468 : WRITE (UNIT=output_unit, FMT="(T2,A,ES16.8,A)") &
469 22 : "MTLR| Maximum change in U: ", &
470 44 : max_delta_u*evolt, " eV"
471 :
472 : WRITE (UNIT=output_unit, FMT="(T2,A,ES16.8,A)") &
473 22 : "MTLR| Maximum change in J: ", &
474 44 : max_delta_j*evolt, " eV"
475 :
476 : WRITE (UNIT=output_unit, FMT="(T2,A,ES16.8,A)") &
477 22 : "MTLR| Convergence threshold: ", &
478 44 : eps_u_j_loop*evolt, " eV"
479 :
480 22 : WRITE (UNIT=output_unit, FMT="(T2,78('-'))")
481 :
482 22 : IF (converged) THEN
483 : WRITE (UNIT=output_unit, FMT="(T2,A)") &
484 8 : "MTLR| U and J have converged."
485 :
486 : WRITE (UNIT=output_unit, FMT="(T2,A)") &
487 8 : "MTLR| Both maximum changes are below the convergence threshold."
488 :
489 14 : ELSE IF (u_iter < max_mtlr_iter) THEN
490 : WRITE (UNIT=output_unit, FMT="(T2,A)") &
491 14 : "MTLR| U and J have not yet converged."
492 :
493 : WRITE (UNIT=output_unit, FMT="(T2,A)") &
494 14 : "MTLR| Proceeding to the next linear response U/J iteration."
495 :
496 : ELSE
497 : WRITE (UNIT=output_unit, FMT="(T2,A)") &
498 0 : "MTLR| U and J have not converged."
499 :
500 : WRITE (UNIT=output_unit, FMT="(T2,A)") &
501 0 : "MTLR| The maximum number of MTLR iterations has been reached."
502 : END IF
503 :
504 22 : WRITE (UNIT=output_unit, FMT="(T2,78('='))")
505 : END IF
506 :
507 44 : IF (converged) THEN
508 16 : IF (output_unit > 0) THEN
509 8 : WRITE (UNIT=output_unit, FMT="(/,T2,78('*'))")
510 : WRITE (UNIT=output_unit, FMT="(T2,A,I0,A)") &
511 8 : "MTLR| U and J converged after ", &
512 16 : u_iter, " iterations."
513 8 : WRITE (UNIT=output_unit, FMT="(T2,78('*'),/)")
514 : END IF
515 28 : ELSE IF (u_iter == max_mtlr_iter) THEN
516 0 : IF (output_unit > 0) THEN
517 0 : WRITE (UNIT=output_unit, FMT="(/,T2,78('*'))")
518 : WRITE (UNIT=output_unit, FMT="(T2,A,I0,A)") &
519 0 : "MTLR| U and J did not converge within the maximum of ", &
520 0 : max_mtlr_iter, " iterations."
521 : WRITE (UNIT=output_unit, FMT="(T2,A)") &
522 0 : "MTLR| Results from the final iteration will be reported."
523 0 : WRITE (UNIT=output_unit, FMT="(T2,78('*'),/)")
524 : END IF
525 : END IF
526 :
527 44 : IF (converged .OR. u_iter == max_mtlr_iter) THEN
528 16 : IF (output_unit > 0) THEN
529 8 : WRITE (UNIT=output_unit, FMT="(/,T2,78('*'))")
530 : WRITE (UNIT=output_unit, FMT="(T2,A,ES16.8,A)") &
531 8 : "MTLR| Maximum change in U: ", &
532 16 : max_delta_u*evolt, " eV"
533 : WRITE (UNIT=output_unit, FMT="(T2,A,ES16.8,A)") &
534 8 : "MTLR| Maximum change in J: ", &
535 16 : max_delta_j*evolt, " eV"
536 : WRITE (UNIT=output_unit, FMT="(T2,A,ES16.8,A)") &
537 8 : "MTLR| Maximum change over U and J: ", &
538 16 : MAX(max_delta_u, max_delta_j)*evolt, " eV"
539 8 : WRITE (UNIT=output_unit, FMT="(T2,78('*'))")
540 18 : DO ikind = 1, nkind
541 10 : IF (.NOT. mtlr_kind(ikind)) CYCLE
542 : WRITE (UNIT=output_unit, FMT="(/,T2,A,I0,A)") &
543 8 : "MTLR| Final parameters for KIND ", ikind, ":"
544 : WRITE (UNIT=output_unit, FMT="(T2,A,F16.8,A)") &
545 8 : "MTLR| Calculated Hubbard U: ", &
546 16 : u_new(ikind)*evolt, " eV"
547 : WRITE (UNIT=output_unit, FMT="(T2,A,F16.8,A)") &
548 8 : "MTLR| Calculated Hund J: ", &
549 16 : j_new(ikind)*evolt, " eV"
550 : WRITE (UNIT=output_unit, FMT="(T2,A)") &
551 8 : "MTLR| Recommended CP2K input parameters:"
552 : WRITE (UNIT=output_unit, FMT="(T2,A,F16.8)") &
553 8 : "MTLR| U_MINUS_J [eV] ", &
554 16 : (u_new(ikind) - j_new(ikind))*evolt
555 : WRITE (UNIT=output_unit, FMT="(T2,A,F16.8)") &
556 8 : "MTLR| J [eV] ", &
557 26 : j_new(ikind)*evolt
558 : END DO
559 8 : WRITE (UNIT=output_unit, FMT="(/,T2,78('*'),/)")
560 : END IF
561 : END IF
562 :
563 44 : IF (converged) THEN
564 : EXIT
565 : END IF
566 :
567 64 : u_old(:) = u_new
568 80 : j_old(:) = j_new
569 :
570 : END DO
571 :
572 36 : DO ikind = 1, nkind
573 20 : IF (.NOT. mtlr_kind(ikind)) CYCLE
574 16 : qs_kind_set(ikind)%dft_plus_u%u_minus_j = u_new(ikind) - j_new(ikind)
575 36 : qs_kind_set(ikind)%dft_plus_u%hund_j = j_new(ikind)
576 : END DO
577 16 : dft_control%mtlr_dft_with_perturbation = .FALSE.
578 16 : dft_control%perturbation_strength = 0.0_dp
579 :
580 16 : CALL force_env_calc_energy_force(force_env, calc_force=.FALSE.)
581 :
582 16 : IF (ALLOCATED(reference_mos)) THEN
583 36 : DO iset = 1, SIZE(reference_mos)
584 36 : CALL deallocate_mo_set(reference_mos(iset))
585 : END DO
586 12 : DEALLOCATE (reference_mos)
587 : END IF
588 16 : DEALLOCATE (u_new)
589 16 : DEALLOCATE (j_new)
590 16 : DEALLOCATE (u_old)
591 16 : DEALLOCATE (j_old)
592 16 : DEALLOCATE (mtlr_kind)
593 :
594 16 : CALL timestop(handle)
595 :
596 32 : END SUBROUTINE do_mtlr_u_j
597 :
598 : ! **************************************************************************************************
599 : !> \brief Restore the fixed reference MOs before an MTLR perturbation.
600 : !> \param[in,out] mos ...
601 : !> \param[in,out] reference_mos ...
602 : ! **************************************************************************************************
603 192 : SUBROUTINE restore_reference_mos(mos, reference_mos)
604 :
605 : TYPE(mo_set_type), DIMENSION(:), INTENT(INOUT) :: mos, reference_mos
606 :
607 : INTEGER :: iset
608 :
609 192 : CPASSERT(SIZE(mos) == SIZE(reference_mos))
610 :
611 576 : DO iset = 1, SIZE(mos)
612 384 : CPASSERT(mos(iset)%use_mo_coeff_b .EQV. reference_mos(iset)%use_mo_coeff_b)
613 384 : CALL reassign_allocated_mos(mos(iset), reference_mos(iset))
614 576 : IF (mos(iset)%use_mo_coeff_b) THEN
615 312 : CPASSERT(ASSOCIATED(mos(iset)%mo_coeff_b))
616 312 : CALL copy_fm_to_dbcsr(mos(iset)%mo_coeff, mos(iset)%mo_coeff_b)
617 : END IF
618 : END DO
619 :
620 192 : END SUBROUTINE restore_reference_mos
621 :
622 : END MODULE mtlr_u_j_methods
|
