(git:f2099e5)
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_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
42CONTAINS
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 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 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 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(mo_set_type), ALLOCATABLE, DIMENSION(:) :: reference_mos
83 TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
84 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
85
86 CALL timeset(routinen, handle)
87
88 NULLIFY (atom_list, qs_kind_set, dft_control, logger, mos)
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 mos=mos, &
100 atomic_kind_set=atomic_kind_set)
101
102 cpassert(ASSOCIATED(atomic_kind_set))
103 cpassert(ASSOCIATED(dft_control))
104 cpassert(ASSOCIATED(mos))
105 cpassert(ASSOCIATED(qs_kind_set))
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 SELECT CASE (dft_control%mtlr_initialization_mode)
133 CASE (mtlr_atomic_perturbations)
134 do_reference_scf = .false.
135 CASE (mtlr_reference_from_atomic, mtlr_reference_from_restart)
136 do_reference_scf = .true.
137 CASE DEFAULT
138 cpabort("The MTLR SCF initialization mode was not resolved.")
139 END SELECT
140
141 IF (output_unit > 0) THEN
142 WRITE (unit=output_unit, fmt="(/,T2,78('='))")
143 WRITE (unit=output_unit, fmt="(T2,A)") &
144 "MTLR| SCF initialization"
145 WRITE (unit=output_unit, fmt="(T2,78('-'))")
146 SELECT CASE (dft_control%mtlr_initialization_mode)
147 CASE (mtlr_atomic_perturbations)
148 WRITE (unit=output_unit, fmt="(T2,A,T66,A13)") &
149 "MTLR| Reference SCF:", "OFF"
150 WRITE (unit=output_unit, fmt="(T2,A,T66,A13)") &
151 "MTLR| Perturbation initial guess:", "ATOMIC"
152 CASE (mtlr_reference_from_atomic)
153 WRITE (unit=output_unit, fmt="(T2,A,T66,A13)") &
154 "MTLR| Reference SCF:", "ON"
155 WRITE (unit=output_unit, fmt="(T2,A,T66,A13)") &
156 "MTLR| Reference SCF initial guess:", "ATOMIC"
157 WRITE (unit=output_unit, fmt="(T2,A,T66,A13)") &
158 "MTLR| Perturbation initial guess:", "REFERENCE MOs"
159 CASE (mtlr_reference_from_restart)
160 WRITE (unit=output_unit, fmt="(T2,A,T66,A13)") &
161 "MTLR| Reference SCF:", "ON"
162 WRITE (unit=output_unit, fmt="(T2,A,T66,A13)") &
163 "MTLR| Reference SCF initial guess:", "RESTART"
164 WRITE (unit=output_unit, fmt="(T2,A,T66,A13)") &
165 "MTLR| Perturbation initial guess:", "REFERENCE MOs"
166 END SELECT
167 WRITE (unit=output_unit, fmt="(T2,78('='))")
168 END IF
169
170 ALLOCATE (u_new(nkind))
171 ALLOCATE (j_new(nkind))
172 ALLOCATE (u_old(nkind))
173 ALLOCATE (j_old(nkind))
174 u_new(:) = 0.0_dp
175 j_new(:) = 0.0_dp
176 u_old(:) = 0.0_dp
177 j_old(:) = 0.0_dp
178 converged = .false.
179 eps_u_j_loop = dft_control%eps_u_j_loop
180 max_mtlr_iter = dft_control%max_mtlr_iter
181 IF (dft_control%nspins /= 2) THEN
182 cpabort("Unrestricted KS has to be used (the number of spin channels should be 2).")
183 END IF
184 IF (max_mtlr_iter < 1) THEN
185 cpabort("MAX_MTLR_LOOP must be at least one.")
186 END IF
187 IF (eps_u_j_loop <= 0.0_dp) THEN
188 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 DO ikind = 1, nkind
194 IF (.NOT. mtlr_kind(ikind)) cycle
195 IF (qs_kind_set(ikind)%dft_plus_u%u_minus_j == 0.0_dp) THEN
196 qs_kind_set(ikind)%dft_plus_u%u_minus_j = epsilon(1.0_dp)
197 END IF
198 IF (qs_kind_set(ikind)%dft_plus_u%hund_j == 0.0_dp) THEN
199 qs_kind_set(ikind)%dft_plus_u%hund_j = epsilon(1.0_dp)
200 END IF
201 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 qs_kind_set(ikind)%dft_plus_u%hund_j
204 END DO
205
206 DO u_iter = 1, max_mtlr_iter
207
208 dft_control%mtlr_dft_with_perturbation = .false.
209
210 IF (do_reference_scf) THEN
211 IF (output_unit > 0) THEN
212 WRITE (unit=output_unit, fmt="(/,T2,78('='))")
213 WRITE (unit=output_unit, fmt="(T2,A)") &
214 "MTLR| Starting the unperturbed reference SCF."
215 WRITE (unit=output_unit, fmt="(T2,A,T72,I8)") &
216 "MTLR| U/J iteration:", u_iter
217 WRITE (unit=output_unit, fmt="(T2,78('='))")
218 END IF
219 CALL force_env_calc_energy_force(force_env, calc_force=.false.)
220
221 IF (.NOT. ALLOCATED(reference_mos)) THEN
222 ALLOCATE (reference_mos(SIZE(mos)))
223 ELSE
224 DO iset = 1, SIZE(reference_mos)
225 CALL deallocate_mo_set(reference_mos(iset))
226 END DO
227 END IF
228 DO iset = 1, SIZE(mos)
229 CALL duplicate_mo_set(reference_mos(iset), mos(iset))
230 END DO
231 END IF
232
233 DO ikind = 1, nkind
234
235 IF (.NOT. mtlr_kind(ikind)) cycle
236 CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom_list)
237 IF (.NOT. any(atom_list == qs_kind_set(ikind)%dft_plus_u%lr_atom)) THEN
238 cpabort("INDEX_PERTURBED_ATOM does not belong to the KIND containing the MTLR section.")
239 END IF
240 IF (.NOT. ALLOCATED( &
241 qs_kind_set(ikind)%dft_plus_u%perturbation_strength)) THEN
242 cpabort("MTLR target does not contain perturbation strengths.")
243 END IF
244
245 dft_control%mtlr_ikind = ikind
246
247 n = SIZE(qs_kind_set(ikind)%dft_plus_u%perturbation_strength)
248 IF (n < 3) THEN
249 cpabort("MTLR linear regression requires at least three perturbation strengths.")
250 END IF
251 l = real(n, dp)
252 ALLOCATE (perturbation_strength(n))
253 ALLOCATE (trq_plus(n))
254 ALLOCATE (vhxc_plus(n))
255 ALLOCATE (trq_minus(n))
256 ALLOCATE (vhxc_minus(n))
257 trq_plus = 0.0_dp
258 trq_minus = 0.0_dp
259 vhxc_plus = 0.0_dp
260 vhxc_minus = 0.0_dp
261 sum_trq_plus = 0.0_dp
262 sum_vhxc_plus = 0.0_dp
263 sum_trq_plus_x_trq_plus = 0.0_dp
264 sum_trq_plus_x_vhxc_plus = 0.0_dp
265 sum_trq_minus = 0.0_dp
266 sum_vhxc_minus = 0.0_dp
267 sum_trq_minus_x_trq_minus = 0.0_dp
268 sum_trq_minus_x_vhxc_minus = 0.0_dp
269 std_plus = 0.0_dp
270 std_minus = 0.0_dp
271 perturbation_strength(:) = qs_kind_set(ikind)%dft_plus_u%perturbation_strength(:)
272
273 DO p_iter = 1, n
274
275 IF (output_unit > 0) THEN
276 WRITE (unit=output_unit, fmt="(/,T2,78('='))")
277 WRITE (unit=output_unit, fmt="(T2,A,T72,I8)") &
278 "MTLR| U/J iteration:", u_iter
279 WRITE (unit=output_unit, fmt="(T2,A,T74,I1,A4,I1)") &
280 "MTLR| Perturbation SCF:", p_iter, " of ", n
281 WRITE (unit=output_unit, fmt="(T2,A,T72,I8)") &
282 "MTLR| Target KIND index:", ikind
283 WRITE (unit=output_unit, fmt="(T2,A,T72,I8)") &
284 "MTLR| Target atom index:", &
285 qs_kind_set(ikind)%dft_plus_u%lr_atom
286 WRITE (unit=output_unit, fmt="(T2,A,T63,F14.8,A3)") &
287 "MTLR| Perturbation strength:", &
288 perturbation_strength(p_iter)*evolt, " eV"
289 WRITE (unit=output_unit, fmt="(T2,78('='))")
290 END IF
291
292 dft_control%perturbation_strength = perturbation_strength(p_iter)
293 dft_control%mtlr_dft_with_perturbation = .true.
294
295 IF (do_reference_scf) CALL restore_reference_mos(mos, reference_mos)
296 CALL force_env_calc_energy_force(force_env, calc_force=.false.)
297
298 trq_plus(p_iter) = dft_control%trq(1) + dft_control%trq(2)
299 trq_minus(p_iter) = dft_control%trq(1) - dft_control%trq(2)
300 vhxc_plus(p_iter) = dft_control%vhxc(1) + dft_control%vhxc(2)
301 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 DO k = 1, n
307 sum_trq_plus_x_vhxc_plus = sum_trq_plus_x_vhxc_plus + trq_plus(k)*vhxc_plus(k)
308 sum_trq_plus = sum_trq_plus + trq_plus(k)
309 sum_vhxc_plus = sum_vhxc_plus + vhxc_plus(k)
310 sum_trq_plus_x_trq_plus = sum_trq_plus_x_trq_plus + trq_plus(k)*trq_plus(k)
311 END DO
312
313 DO k = 1, n
314 sum_trq_minus_x_vhxc_minus = sum_trq_minus_x_vhxc_minus + trq_minus(k)*vhxc_minus(k)
315 sum_trq_minus = sum_trq_minus + trq_minus(k)
316 sum_vhxc_minus = sum_vhxc_minus + vhxc_minus(k)
317 sum_trq_minus_x_trq_minus = sum_trq_minus_x_trq_minus + trq_minus(k)*trq_minus(k)
318 END DO
319
320 denominator_plus = l*sum_trq_plus_x_trq_plus - sum_trq_plus**2.0_dp
321 denominator_minus = l*sum_trq_minus_x_trq_minus - sum_trq_minus**2.0_dp
322 IF (abs(denominator_plus) < 100.0_dp*epsilon(1.0_dp)) THEN
323 cpabort("MTLR regression is singular.")
324 END IF
325 IF (abs(denominator_minus) < 100.0_dp*epsilon(1.0_dp)) THEN
326 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 /denominator_minus
330 fhxc_plus = (l*sum_trq_plus_x_vhxc_plus - sum_trq_plus*sum_vhxc_plus) &
331 /denominator_plus
332 intercept_minus = (sum_vhxc_minus - fhxc_minus*sum_trq_minus)/l
333 intercept_plus = (sum_vhxc_plus - fhxc_plus*sum_trq_plus)/l
334 DO k = 1, n
335 std_minus = std_minus + (vhxc_minus(k) - (trq_minus(k)*fhxc_minus + intercept_minus))**2/(l - 2)
336 END DO
337 DO k = 1, n
338 std_plus = std_plus + (vhxc_plus(k) - (trq_plus(k)*fhxc_plus + intercept_plus))**2/(l - 2)
339 END DO
340 std_minus = sqrt(std_minus)
341 std_plus = sqrt(std_plus)
342
343 u_new(ikind) = 0.5_dp*fhxc_plus
344 j_new(ikind) = -0.5_dp*fhxc_minus
345 delta_u = u_new(ikind) - u_old(ikind)
346 delta_j = j_new(ikind) - j_old(ikind)
347 std_u = 0.5_dp*std_plus
348 std_j = 0.5_dp*std_minus
349
350 IF (output_unit > 0) THEN
351 WRITE (unit=output_unit, fmt="(/,T2,78('='))")
352 WRITE (unit=output_unit, fmt="(T2,A,I0,A,I0,A)") &
353 "MTLR| U/J iteration ", u_iter, &
354 " results for KIND ", ikind, "."
355 WRITE (unit=output_unit, fmt="(T2,78('-'))")
356 WRITE (unit=output_unit, fmt="(T2,A)") &
357 "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 "Pt", &
361 "Pert.[eV]", &
362 "trq(+)", &
363 "vhxc(+)[eV]", &
364 "vhxc(+)-fit", &
365 "vhxc(+)-resid"
366 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 k, &
371 perturbation_strength(k)*evolt, &
372 trq_plus(k), &
373 vhxc_plus(k)*evolt, &
374 (fhxc_plus*trq_plus(k) + intercept_plus)*evolt, &
375 (vhxc_plus(k) - &
376 (fhxc_plus*trq_plus(k) + intercept_plus))*evolt
377 END DO
378 WRITE (unit=output_unit, fmt="(T2,78('-'))")
379 WRITE (unit=output_unit, fmt="(T2,A)") &
380 "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 "Pt", &
384 "Pert.[eV]", &
385 "trq(-)", &
386 "vhxc(-)[eV]", &
387 "vhxc(-)-fit", &
388 "vhxc(-)-resid"
389 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 k, &
394 perturbation_strength(k)*evolt, &
395 trq_minus(k), &
396 vhxc_minus(k)*evolt, &
397 (fhxc_minus*trq_minus(k) + intercept_minus)*evolt, &
398 (vhxc_minus(k) - &
399 (fhxc_minus*trq_minus(k) + intercept_minus))*evolt
400 END DO
401 WRITE (unit=output_unit, fmt="(T2,78('-'))")
402 WRITE (unit=output_unit, &
403 fmt="(T2,A,T22,A16,T43,A16,T64,A16)") &
404 "MTLR| Parameter", &
405 "Old [eV]", &
406 "New [eV]", &
407 "Change [eV]"
408 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 "MTLR| Hubbard U", &
412 u_old(ikind)*evolt, &
413 u_new(ikind)*evolt, &
414 (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 "MTLR| Hund J", &
418 j_old(ikind)*evolt, &
419 j_new(ikind)*evolt, &
420 (j_new(ikind) - j_old(ikind))*evolt
421 WRITE (unit=output_unit, fmt="(T2,78('-'))")
422 WRITE (unit=output_unit, fmt="(T2,A,T61,ES16.8,A3)") &
423 "MTLR| Absolute change in U:", &
424 abs(u_new(ikind) - u_old(ikind))*evolt, " eV"
425 WRITE (unit=output_unit, fmt="(T2,A,T61,ES16.8,A3)") &
426 "MTLR| Absolute change in J:", &
427 abs(j_new(ikind) - j_old(ikind))*evolt, " eV"
428 WRITE (unit=output_unit, fmt="(T2,A,T61,ES16.8,A3)") &
429 "MTLR| U fit residual:", &
430 std_u*evolt, " eV"
431 WRITE (unit=output_unit, fmt="(T2,A,T61,ES16.8,A3)") &
432 "MTLR| J fit residual:", &
433 std_j*evolt, " eV"
434 WRITE (unit=output_unit, fmt="(T2,78('='))")
435 END IF
436
437 DEALLOCATE (perturbation_strength)
438 DEALLOCATE (trq_plus)
439 DEALLOCATE (vhxc_plus)
440 DEALLOCATE (trq_minus)
441 DEALLOCATE (vhxc_minus)
442
443 END DO
444
445 IF (do_reference_scf) CALL restore_reference_mos(mos, reference_mos)
446
447 DO ikind = 1, nkind
448 IF (.NOT. mtlr_kind(ikind)) cycle
449 qs_kind_set(ikind)%dft_plus_u%u_minus_j = u_new(ikind) - j_new(ikind)
450 qs_kind_set(ikind)%dft_plus_u%hund_j = j_new(ikind)
451 END DO
452
453 max_delta_u = maxval(abs(pack(u_new(:) - u_old(:), mtlr_kind)))
454 max_delta_j = maxval(abs(pack(j_new(:) - j_old(:), mtlr_kind)))
455 IF (max_delta_u < eps_u_j_loop .AND. max_delta_j < eps_u_j_loop) THEN
456 converged = .true.
457 END IF
458
459 IF (output_unit > 0) THEN
460 WRITE (unit=output_unit, fmt="(/,T2,78('='))")
461
462 WRITE (unit=output_unit, fmt="(T2,A,I0,A)") &
463 "MTLR| Iteration ", u_iter, &
464 " completed for all atomic kinds."
465
466 WRITE (unit=output_unit, fmt="(T2,78('-'))")
467
468 WRITE (unit=output_unit, fmt="(T2,A,ES16.8,A)") &
469 "MTLR| Maximum change in U: ", &
470 max_delta_u*evolt, " eV"
471
472 WRITE (unit=output_unit, fmt="(T2,A,ES16.8,A)") &
473 "MTLR| Maximum change in J: ", &
474 max_delta_j*evolt, " eV"
475
476 WRITE (unit=output_unit, fmt="(T2,A,ES16.8,A)") &
477 "MTLR| Convergence threshold: ", &
478 eps_u_j_loop*evolt, " eV"
479
480 WRITE (unit=output_unit, fmt="(T2,78('-'))")
481
482 IF (converged) THEN
483 WRITE (unit=output_unit, fmt="(T2,A)") &
484 "MTLR| U and J have converged."
485
486 WRITE (unit=output_unit, fmt="(T2,A)") &
487 "MTLR| Both maximum changes are below the convergence threshold."
488
489 ELSE IF (u_iter < max_mtlr_iter) THEN
490 WRITE (unit=output_unit, fmt="(T2,A)") &
491 "MTLR| U and J have not yet converged."
492
493 WRITE (unit=output_unit, fmt="(T2,A)") &
494 "MTLR| Proceeding to the next linear response U/J iteration."
495
496 ELSE
497 WRITE (unit=output_unit, fmt="(T2,A)") &
498 "MTLR| U and J have not converged."
499
500 WRITE (unit=output_unit, fmt="(T2,A)") &
501 "MTLR| The maximum number of MTLR iterations has been reached."
502 END IF
503
504 WRITE (unit=output_unit, fmt="(T2,78('='))")
505 END IF
506
507 IF (converged) THEN
508 IF (output_unit > 0) THEN
509 WRITE (unit=output_unit, fmt="(/,T2,78('*'))")
510 WRITE (unit=output_unit, fmt="(T2,A,I0,A)") &
511 "MTLR| U and J converged after ", &
512 u_iter, " iterations."
513 WRITE (unit=output_unit, fmt="(T2,78('*'),/)")
514 END IF
515 ELSE IF (u_iter == max_mtlr_iter) THEN
516 IF (output_unit > 0) THEN
517 WRITE (unit=output_unit, fmt="(/,T2,78('*'))")
518 WRITE (unit=output_unit, fmt="(T2,A,I0,A)") &
519 "MTLR| U and J did not converge within the maximum of ", &
520 max_mtlr_iter, " iterations."
521 WRITE (unit=output_unit, fmt="(T2,A)") &
522 "MTLR| Results from the final iteration will be reported."
523 WRITE (unit=output_unit, fmt="(T2,78('*'),/)")
524 END IF
525 END IF
526
527 IF (converged .OR. u_iter == max_mtlr_iter) THEN
528 IF (output_unit > 0) THEN
529 WRITE (unit=output_unit, fmt="(/,T2,78('*'))")
530 WRITE (unit=output_unit, fmt="(T2,A,ES16.8,A)") &
531 "MTLR| Maximum change in U: ", &
532 max_delta_u*evolt, " eV"
533 WRITE (unit=output_unit, fmt="(T2,A,ES16.8,A)") &
534 "MTLR| Maximum change in J: ", &
535 max_delta_j*evolt, " eV"
536 WRITE (unit=output_unit, fmt="(T2,A,ES16.8,A)") &
537 "MTLR| Maximum change over U and J: ", &
538 max(max_delta_u, max_delta_j)*evolt, " eV"
539 WRITE (unit=output_unit, fmt="(T2,78('*'))")
540 DO ikind = 1, nkind
541 IF (.NOT. mtlr_kind(ikind)) cycle
542 WRITE (unit=output_unit, fmt="(/,T2,A,I0,A)") &
543 "MTLR| Final parameters for KIND ", ikind, ":"
544 WRITE (unit=output_unit, fmt="(T2,A,F16.8,A)") &
545 "MTLR| Calculated Hubbard U: ", &
546 u_new(ikind)*evolt, " eV"
547 WRITE (unit=output_unit, fmt="(T2,A,F16.8,A)") &
548 "MTLR| Calculated Hund J: ", &
549 j_new(ikind)*evolt, " eV"
550 WRITE (unit=output_unit, fmt="(T2,A)") &
551 "MTLR| Recommended CP2K input parameters:"
552 WRITE (unit=output_unit, fmt="(T2,A,F16.8)") &
553 "MTLR| U_MINUS_J [eV] ", &
554 (u_new(ikind) - j_new(ikind))*evolt
555 WRITE (unit=output_unit, fmt="(T2,A,F16.8)") &
556 "MTLR| J [eV] ", &
557 j_new(ikind)*evolt
558 END DO
559 WRITE (unit=output_unit, fmt="(/,T2,78('*'),/)")
560 END IF
561 END IF
562
563 IF (converged) THEN
564 EXIT
565 END IF
566
567 u_old(:) = u_new
568 j_old(:) = j_new
569
570 END DO
571
572 DO ikind = 1, nkind
573 IF (.NOT. mtlr_kind(ikind)) cycle
574 qs_kind_set(ikind)%dft_plus_u%u_minus_j = u_new(ikind) - j_new(ikind)
575 qs_kind_set(ikind)%dft_plus_u%hund_j = j_new(ikind)
576 END DO
577 dft_control%mtlr_dft_with_perturbation = .false.
578 dft_control%perturbation_strength = 0.0_dp
579
580 CALL force_env_calc_energy_force(force_env, calc_force=.false.)
581
582 IF (ALLOCATED(reference_mos)) THEN
583 DO iset = 1, SIZE(reference_mos)
584 CALL deallocate_mo_set(reference_mos(iset))
585 END DO
586 DEALLOCATE (reference_mos)
587 END IF
588 DEALLOCATE (u_new)
589 DEALLOCATE (j_new)
590 DEALLOCATE (u_old)
591 DEALLOCATE (j_old)
592 DEALLOCATE (mtlr_kind)
593
594 CALL timestop(handle)
595
596 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 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 cpassert(SIZE(mos) == SIZE(reference_mos))
610
611 DO iset = 1, SIZE(mos)
612 cpassert(mos(iset)%use_mo_coeff_b .EQV. reference_mos(iset)%use_mo_coeff_b)
613 CALL reassign_allocated_mos(mos(iset), reference_mos(iset))
614 IF (mos(iset)%use_mo_coeff_b) THEN
615 cpassert(ASSOCIATED(mos(iset)%mo_coeff_b))
616 CALL copy_fm_to_dbcsr(mos(iset)%mo_coeff, mos(iset)%mo_coeff_b)
617 END IF
618 END DO
619
620 END SUBROUTINE restore_reference_mos
621
622END 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_dbcsr_operations
DBCSR operations in CP2K.
Definition cp_dbcsr_operations.F:18
cp_dbcsr_operations::copy_fm_to_dbcsr
subroutine, public copy_fm_to_dbcsr(fm, matrix, keep_sparsity)
Copy a BLACS matrix to a dbcsr matrix.
Definition cp_dbcsr_operations.F:129
cp_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
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:18
input_constants::mtlr_reference_from_restart
integer, parameter, public mtlr_reference_from_restart
Definition input_constants.F:226
input_constants::mtlr_atomic_perturbations
integer, parameter, public mtlr_atomic_perturbations
Definition input_constants.F:226
input_constants::mtlr_reference_from_atomic
integer, parameter, public mtlr_reference_from_atomic
Definition input_constants.F:226
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:60
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:545
qs_kind_types
Define the quickstep kind type and their sub types.
Definition qs_kind_types.F:23
qs_mo_types
Definition and initialisation of the mo data type.
Definition qs_mo_types.F:22
qs_mo_types::duplicate_mo_set
subroutine, public duplicate_mo_set(mo_set_new, mo_set_old)
allocate a new mo_set, and copy the old data
Definition qs_mo_types.F:162
qs_mo_types::deallocate_mo_set
subroutine, public deallocate_mo_set(mo_set)
Deallocate a wavefunction data structure.
Definition qs_mo_types.F:394
qs_mo_types::reassign_allocated_mos
subroutine, public reassign_allocated_mos(mo_set_new, mo_set_old)
reassign an already allocated mo_set
Definition qs_mo_types.F:117
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:780
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
qs_kind_types::qs_kind_type
Provides all information about a quickstep kind.
Definition qs_kind_types.F:190
qs_mo_types::mo_set_type
Definition qs_mo_types.F:53