(git:98357aa)
Loading...
Searching...
No Matches
qs_outer_scf.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! **************************************************************************************************
9!> \brief Routines for performing an outer scf loop
10!> \par History
11!> Created [2006.03]
12!> \author Joost VandeVondele
13! **************************************************************************************************
14MODULE qs_outer_scf
15 USE cp_control_types, ONLY: ddapc_restraint_type,&
16 dft_control_type,&
17 s2_restraint_type
18 USE cp_log_handling, ONLY: cp_to_string
19 USE input_constants, ONLY: &
20 broyden_type_1, broyden_type_1_explicit, broyden_type_1_explicit_ls, broyden_type_1_ls, &
21 broyden_type_2, broyden_type_2_explicit, broyden_type_2_explicit_ls, broyden_type_2_ls, &
22 cdft2ot, do_ddapc_constraint, do_s2_constraint, ot2cdft, outer_scf_basis_center_opt, &
23 outer_scf_cdft_constraint, outer_scf_ddapc_constraint, outer_scf_none, &
24 outer_scf_optimizer_bisect, outer_scf_optimizer_broyden, outer_scf_optimizer_diis, &
25 outer_scf_optimizer_newton, outer_scf_optimizer_newton_ls, outer_scf_optimizer_none, &
26 outer_scf_optimizer_sd, outer_scf_optimizer_secant, outer_scf_s2_constraint
27 USE kinds, ONLY: dp
28 USE mathlib, ONLY: diamat_all
29 USE qs_basis_gradient, ONLY: qs_basis_center_gradient,&
30 qs_update_basis_center_pos,&
31 return_basis_center_gradient_norm
32 USE qs_cdft_opt_types, ONLY: cdft_opt_type_copy,&
33 cdft_opt_type_release
34 USE qs_cdft_types, ONLY: cdft_control_type
35 USE qs_energy_types, ONLY: qs_energy_type
36 USE qs_environment_types, ONLY: get_qs_env,&
37 qs_environment_type,&
38 set_qs_env
39 USE qs_scf_types, ONLY: qs_scf_env_type
40 USE scf_control_types, ONLY: scf_control_type
41#include "./base/base_uses.f90"
42
43 IMPLICIT NONE
44
45 PRIVATE
46
47! *** Global parameters ***
48
49 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_outer_scf'
50
51! *** Public subroutines ***
52
53 PUBLIC :: outer_loop_gradient, outer_loop_optimize, outer_loop_update_qs_env, &
54 outer_loop_variables_count, outer_loop_extrapolate, &
55 outer_loop_switch, outer_loop_purge_history
56
57CONTAINS
58
59! **************************************************************************************************
60!> \brief returns the number of variables that is employed in the outer loop. with a CDFT constraint
61!> this value is returned by the cdft_control type
62!> \param scf_control the outer loop control type
63!> \param cdft_control the cdft loop control type
64!> \return the number of variables
65!> \par History
66!> 03.2006 created [Joost VandeVondele]
67! **************************************************************************************************
68 FUNCTION outer_loop_variables_count(scf_control, cdft_control) RESULT(res)
69 TYPE(scf_control_type), POINTER :: scf_control
70 TYPE(cdft_control_type), INTENT(IN), OPTIONAL, &
71 POINTER :: cdft_control
72 INTEGER :: res
73
74 SELECT CASE (scf_control%outer_scf%type)
75 CASE (outer_scf_ddapc_constraint)
76 res = 1
77 CASE (outer_scf_s2_constraint)
78 res = 1
79 CASE (outer_scf_cdft_constraint)
80 IF (PRESENT(cdft_control)) THEN
81 res = SIZE(cdft_control%target)
82 ELSE
83 res = 1
84 END IF
85 CASE (outer_scf_basis_center_opt)
86 res = 1
87 CASE (outer_scf_none) ! just needed to communicate the gradient criterion
88 res = 1
89 CASE DEFAULT
90 res = 0
91 END SELECT
92
93 END FUNCTION outer_loop_variables_count
94
95! **************************************************************************************************
96!> \brief computes the gradient wrt to the outer loop variables
97!> \param qs_env ...
98!> \param scf_env ...
99!> \par History
100!> 03.2006 created [Joost VandeVondele]
101! **************************************************************************************************
102 SUBROUTINE outer_loop_gradient(qs_env, scf_env)
103 TYPE(qs_environment_type), POINTER :: qs_env
104 TYPE(qs_scf_env_type), POINTER :: scf_env
105
106 CHARACTER(LEN=*), PARAMETER :: routinen = 'outer_loop_gradient'
107
108 INTEGER :: handle, ihistory, ivar, n
109 LOGICAL :: is_constraint
110 TYPE(cdft_control_type), POINTER :: cdft_control
111 TYPE(ddapc_restraint_type), POINTER :: ddapc_restraint_control
112 TYPE(dft_control_type), POINTER :: dft_control
113 TYPE(qs_energy_type), POINTER :: energy
114 TYPE(s2_restraint_type), POINTER :: s2_restraint_control
115 TYPE(scf_control_type), POINTER :: scf_control
116
117 CALL timeset(routinen, handle)
118
119 CALL get_qs_env(qs_env=qs_env, scf_control=scf_control, &
120 dft_control=dft_control, energy=energy)
121 cpassert(scf_control%outer_scf%have_scf)
122
123 ihistory = scf_env%outer_scf%iter_count
124 cpassert(ihistory <= SIZE(scf_env%outer_scf%energy, 1))
125
126 scf_env%outer_scf%energy(ihistory) = energy%total
127
128 SELECT CASE (scf_control%outer_scf%type)
129 CASE (outer_scf_none)
130 ! just pass the inner loop scf criterion to the outer loop one
131 scf_env%outer_scf%variables(1, ihistory) = scf_env%iter_delta
132 scf_env%outer_scf%gradient(1, ihistory) = scf_env%iter_delta
133 CASE (outer_scf_ddapc_constraint)
134 cpassert(dft_control%qs_control%ddapc_restraint)
135 DO n = 1, SIZE(dft_control%qs_control%ddapc_restraint_control)
136 NULLIFY (ddapc_restraint_control)
137 ddapc_restraint_control => dft_control%qs_control%ddapc_restraint_control(n)
138 is_constraint = (ddapc_restraint_control%functional_form == do_ddapc_constraint)
139 IF (is_constraint) EXIT
140 END DO
141 cpassert(is_constraint)
142
143 scf_env%outer_scf%variables(:, ihistory) = ddapc_restraint_control%strength
144 scf_env%outer_scf%gradient(:, ihistory) = ddapc_restraint_control%ddapc_order_p - &
145 ddapc_restraint_control%target
146 CASE (outer_scf_s2_constraint)
147 cpassert(dft_control%qs_control%s2_restraint)
148 s2_restraint_control => dft_control%qs_control%s2_restraint_control
149 is_constraint = (s2_restraint_control%functional_form == do_s2_constraint)
150 cpassert(is_constraint)
151
152 scf_env%outer_scf%variables(:, ihistory) = s2_restraint_control%strength
153 scf_env%outer_scf%gradient(:, ihistory) = s2_restraint_control%s2_order_p - &
154 s2_restraint_control%target
155 CASE (outer_scf_cdft_constraint)
156 cpassert(dft_control%qs_control%cdft)
157 cdft_control => dft_control%qs_control%cdft_control
158 DO ivar = 1, SIZE(scf_env%outer_scf%gradient, 1)
159 scf_env%outer_scf%variables(ivar, ihistory) = cdft_control%strength(ivar)
160 scf_env%outer_scf%gradient(ivar, ihistory) = cdft_control%value(ivar) - &
161 cdft_control%target(ivar)
162 END DO
163 CASE (outer_scf_basis_center_opt)
164 CALL qs_basis_center_gradient(qs_env)
165 scf_env%outer_scf%gradient(:, ihistory) = return_basis_center_gradient_norm(qs_env)
166
167 CASE DEFAULT
168 cpabort("Unknown outer SCF type")
169
170 END SELECT
171
172 CALL timestop(handle)
173
174 END SUBROUTINE outer_loop_gradient
175
176! **************************************************************************************************
177!> \brief optimizes the parameters of the outer_scf
178!> \param scf_env the scf_env where to optimize the parameters
179!> \param scf_control control parameters for the optimization
180!> \par History
181!> 03.2006 created [Joost VandeVondele]
182!> 01.2017 added Broyden and Newton optimizers [Nico Holmberg]
183!> \note
184!> ought to be general, and independent of the actual kind of variables
185! **************************************************************************************************
186 SUBROUTINE outer_loop_optimize(scf_env, scf_control)
187 TYPE(qs_scf_env_type), POINTER :: scf_env
188 TYPE(scf_control_type), POINTER :: scf_control
189
190 CHARACTER(LEN=*), PARAMETER :: routinen = 'outer_loop_optimize'
191
192 INTEGER :: handle, i, ibuf, ihigh, ihistory, ilow, &
193 j, jbuf, nb, nvar, optimizer_type
194 REAL(kind=dp) :: interval, scale, tmp
195 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: ev
196 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: a, b, f, x
197 REAL(kind=dp), DIMENSION(:, :), POINTER :: inv_jacobian
198
199 CALL timeset(routinen, handle)
200
201 ihistory = scf_env%outer_scf%iter_count
202 optimizer_type = scf_control%outer_scf%optimizer
203 NULLIFY (inv_jacobian)
204
205 IF (scf_control%outer_scf%type == outer_scf_basis_center_opt) THEN
206 scf_env%outer_scf%variables(:, ihistory + 1) = scf_env%outer_scf%variables(:, ihistory)
207 ELSE
208 DO WHILE (.true.) ! if we need a different run type we'll restart here
209
210 SELECT CASE (optimizer_type)
211 CASE (outer_scf_optimizer_bisect) ! bisection on the gradient, needs to be 1D
212 cpassert(SIZE(scf_env%outer_scf%gradient(:, 1)) == 1)
213 ! find the pair of points that bracket a zero of the gradient, with the smallest interval possible
214 ilow = -1
215 ihigh = -1
216 interval = huge(interval)
217 DO i = 1, ihistory
218 DO j = i + 1, ihistory
219 ! distrust often used points
220 IF (scf_env%outer_scf%count(i) > scf_control%outer_scf%bisect_trust_count) cycle
221 IF (scf_env%outer_scf%count(j) > scf_control%outer_scf%bisect_trust_count) cycle
222
223 ! if they bracket a zero use them
224 IF (scf_env%outer_scf%gradient(1, i)* &
225 scf_env%outer_scf%gradient(1, j) < 0.0_dp) THEN
226 tmp = abs(scf_env%outer_scf%variables(1, i) - scf_env%outer_scf%variables(1, j))
227 IF (tmp < interval) THEN
228 ilow = i
229 ihigh = j
230 interval = tmp
231 END IF
232 END IF
233 END DO
234 END DO
235 IF (ilow == -1) THEN ! we didn't bracket a minimum yet, try something else
236 optimizer_type = outer_scf_optimizer_diis
237 cycle
238 END IF
239 scf_env%outer_scf%count(ilow) = scf_env%outer_scf%count(ilow) + 1
240 scf_env%outer_scf%count(ihigh) = scf_env%outer_scf%count(ihigh) + 1
241 scf_env%outer_scf%variables(:, ihistory + 1) = 0.5_dp*(scf_env%outer_scf%variables(:, ilow) + &
242 scf_env%outer_scf%variables(:, ihigh))
243 CASE (outer_scf_optimizer_none)
244 scf_env%outer_scf%variables(:, ihistory + 1) = scf_env%outer_scf%variables(:, ihistory)
245 CASE (outer_scf_optimizer_sd)
246 ! Notice that we are just trying to find a stationary point
247 ! e.g. the ddpac_constraint, one maximizes the function, so the stepsize might have
248 ! to be negative
249 scf_env%outer_scf%variables(:, ihistory + 1) = scf_env%outer_scf%variables(:, ihistory) - &
250 scf_control%outer_scf%step_size*scf_env%outer_scf%gradient(:, ihistory)
251 CASE (outer_scf_optimizer_diis)
252 cpassert(scf_control%outer_scf%diis_buffer_length > 0)
253 ! set up DIIS matrix
254 nb = min(ihistory, scf_control%outer_scf%diis_buffer_length)
255 IF (nb < 2) THEN
256 optimizer_type = outer_scf_optimizer_sd
257 cycle
258 ELSE
259 ALLOCATE (b(nb + 1, nb + 1), a(nb + 1, nb + 1), ev(nb + 1))
260 DO i = 1, nb
261 DO j = i, nb
262 ibuf = ihistory - nb + i
263 jbuf = ihistory - nb + j
264 b(i, j) = dot_product(scf_env%outer_scf%gradient(:, ibuf), &
265 scf_env%outer_scf%gradient(:, jbuf))
266 b(j, i) = b(i, j)
267 END DO
268 END DO
269 b(nb + 1, :) = -1.0_dp
270 b(:, nb + 1) = -1.0_dp
271 b(nb + 1, nb + 1) = 0.0_dp
272
273 CALL diamat_all(b, ev)
274 a(:, :) = b
275 DO i = 1, nb + 1
276 IF (abs(ev(i)) < 1.0e-12_dp) THEN
277 a(:, i) = 0.0_dp
278 ELSE
279 a(:, i) = a(:, i)/ev(i)
280 END IF
281 END DO
282 ev(:) = -matmul(a, b(nb + 1, :))
283
284 scf_env%outer_scf%variables(:, ihistory + 1) = 0.0_dp
285 DO i = 1, nb
286 ibuf = ihistory - nb + i
287 scf_env%outer_scf%variables(:, ihistory + 1) = scf_env%outer_scf%variables(:, ihistory + 1) + &
288 ev(i)*scf_env%outer_scf%variables(:, ibuf)
289 END DO
290 DEALLOCATE (a, b, ev)
291 END IF
292 CASE (outer_scf_optimizer_secant)
293 cpassert(SIZE(scf_env%outer_scf%gradient, 2) >= 3)
294 cpassert(SIZE(scf_env%outer_scf%gradient, 1) == 1)
295 nvar = SIZE(scf_env%outer_scf%gradient, 1)
296 IF (ihistory < 2) THEN
297 ! Need two history values to use secant, switch to sd
298 optimizer_type = outer_scf_optimizer_sd
299 cycle
300 END IF
301 ! secant update
302 scf_env%outer_scf%variables(1, ihistory + 1) = scf_env%outer_scf%variables(1, ihistory) - &
303 (scf_env%outer_scf%variables(1, ihistory) - &
304 scf_env%outer_scf%variables(1, ihistory - 1))/ &
305 (scf_env%outer_scf%gradient(1, ihistory) - &
306 scf_env%outer_scf%gradient(1, ihistory - 1))* &
307 scf_env%outer_scf%gradient(1, ihistory)
308 CASE (outer_scf_optimizer_broyden)
309 IF (.NOT. ASSOCIATED(scf_env%outer_scf%inv_jacobian)) THEN
310 ! Inverse Jacobian not yet built, switch to sd
311 optimizer_type = outer_scf_optimizer_sd
312 cycle
313 END IF
314 inv_jacobian => scf_env%outer_scf%inv_jacobian
315 IF (ihistory < 2) THEN
316 ! Cannot perform a Broyden update without enough SCF history on this energy evaluation
317 scf_control%outer_scf%cdft_opt_control%broyden_update = .false.
318 END IF
319 IF (scf_control%outer_scf%cdft_opt_control%broyden_update) THEN
320 ! Perform a Broyden update of the inverse Jacobian J^(-1)
321 IF (SIZE(scf_env%outer_scf%gradient, 2) < 3) THEN
322 CALL cp_abort(__location__, &
323 "Keyword EXTRAPOLATION_ORDER in section OUTER_SCF "// &
324 "must be greater than or equal to 3 for Broyden optimizers.")
325 END IF
326 nvar = SIZE(scf_env%outer_scf%gradient, 1)
327 ALLOCATE (f(nvar, 1), x(nvar, 1))
328 DO i = 1, nvar
329 f(i, 1) = scf_env%outer_scf%gradient(i, ihistory) - scf_env%outer_scf%gradient(i, ihistory - 1)
330 x(i, 1) = scf_env%outer_scf%variables(i, ihistory) - scf_env%outer_scf%variables(i, ihistory - 1)
331 END DO
332 SELECT CASE (scf_control%outer_scf%cdft_opt_control%broyden_type)
333 CASE (broyden_type_1, broyden_type_1_explicit, broyden_type_1_ls, broyden_type_1_explicit_ls)
334 ! Broyden's 1st method
335 ! Denote: dx_n = \delta x_n; df_n = \delta f_n
336 ! J_(n+1)^(-1) = J_n^(-1) + (dx_n - J_n^(-1)*df_n)*(dx_n^T * J_n^(-1))/(dx_n^T * J_n^(-1) * df_n)
337 scale = sum(matmul(transpose(x), matmul(inv_jacobian, f)))
338 scale = 1.0_dp/scale
339 IF (scale < 1.0e-12_dp) scale = 1.0e-12_dp
340 inv_jacobian = inv_jacobian + scale*matmul((x - matmul(inv_jacobian, f)), &
341 matmul(transpose(x), inv_jacobian))
342 CASE (broyden_type_2, broyden_type_2_explicit, broyden_type_2_ls, broyden_type_2_explicit_ls)
343 ! Broyden's 2nd method
344 ! J_(n+1)^(-1) = J_n^(-1) + (dx_n - J_n^(-1)*df_n)*(df_n^T)/(||df_n||^2)
345 scale = sum(matmul(transpose(f), f))
346 scale = 1.0_dp/scale
347 IF (scale < 1.0e-12_dp) scale = 1.0e-12_dp
348 inv_jacobian = inv_jacobian + scale*matmul((x - matmul(inv_jacobian, f)), transpose(inv_jacobian))
349 CASE DEFAULT
350 CALL cp_abort(__location__, &
351 "Unknown Broyden type: "// &
352 cp_to_string(scf_control%outer_scf%cdft_opt_control%broyden_type))
353 END SELECT
354 ! Clean up
355 DEALLOCATE (f, x)
356 END IF
357 ! Update variables x_(n+1) = x_n - J^(-1)*f(x_n)
358 scf_env%outer_scf%variables(:, ihistory + 1) = scf_env%outer_scf%variables(:, ihistory) - &
359 scf_control%outer_scf%cdft_opt_control%newton_step* &
360 matmul(inv_jacobian, scf_env%outer_scf%gradient(:, ihistory))
361 scf_control%outer_scf%cdft_opt_control%broyden_update = .true.
362 CASE (outer_scf_optimizer_newton, outer_scf_optimizer_newton_ls)
363 cpassert(ASSOCIATED(scf_env%outer_scf%inv_jacobian))
364 inv_jacobian => scf_env%outer_scf%inv_jacobian
365 scf_env%outer_scf%variables(:, ihistory + 1) = scf_env%outer_scf%variables(:, ihistory) - &
366 scf_control%outer_scf%cdft_opt_control%newton_step* &
367 matmul(inv_jacobian, scf_env%outer_scf%gradient(:, ihistory))
368 CASE DEFAULT
369 cpabort("Unknown outer SCF optimizer")
370 END SELECT
371 EXIT
372 END DO
373 END IF
374
375 CALL timestop(handle)
376
377 END SUBROUTINE outer_loop_optimize
378
379! **************************************************************************************************
380!> \brief propagates the updated variables to wherever they need to be set in
381!> qs_env
382!> \param qs_env ...
383!> \param scf_env ...
384!> \par History
385!> 03.2006 created [Joost VandeVondele]
386! **************************************************************************************************
387 SUBROUTINE outer_loop_update_qs_env(qs_env, scf_env)
388 TYPE(qs_environment_type), POINTER :: qs_env
389 TYPE(qs_scf_env_type), POINTER :: scf_env
390
391 CHARACTER(LEN=*), PARAMETER :: routinen = 'outer_loop_update_qs_env'
392
393 INTEGER :: handle, ihistory, n
394 LOGICAL :: is_constraint
395 TYPE(cdft_control_type), POINTER :: cdft_control
396 TYPE(ddapc_restraint_type), POINTER :: ddapc_restraint_control
397 TYPE(dft_control_type), POINTER :: dft_control
398 TYPE(s2_restraint_type), POINTER :: s2_restraint_control
399 TYPE(scf_control_type), POINTER :: scf_control
400
401 CALL timeset(routinen, handle)
402
403 CALL get_qs_env(qs_env=qs_env, scf_control=scf_control, dft_control=dft_control)
404 ihistory = scf_env%outer_scf%iter_count
405
406 SELECT CASE (scf_control%outer_scf%type)
407 CASE (outer_scf_none)
408 ! do nothing
409 CASE (outer_scf_ddapc_constraint)
410 DO n = 1, SIZE(dft_control%qs_control%ddapc_restraint_control)
411 NULLIFY (ddapc_restraint_control)
412 ddapc_restraint_control => dft_control%qs_control%ddapc_restraint_control(n)
413 is_constraint = (ddapc_restraint_control%functional_form == do_ddapc_constraint)
414 IF (is_constraint) EXIT
415 END DO
416 ddapc_restraint_control%strength = scf_env%outer_scf%variables(1, ihistory + 1)
417 CASE (outer_scf_s2_constraint)
418 s2_restraint_control => dft_control%qs_control%s2_restraint_control
419 s2_restraint_control%strength = scf_env%outer_scf%variables(1, ihistory + 1)
420 CASE (outer_scf_cdft_constraint)
421 cdft_control => dft_control%qs_control%cdft_control
422 cdft_control%strength(:) = scf_env%outer_scf%variables(:, ihistory + 1)
423 CASE (outer_scf_basis_center_opt)
424 CALL qs_update_basis_center_pos(qs_env)
425 CASE DEFAULT
426 cpabort("Unknown outer SCF type")
427 END SELECT
428
429 CALL timestop(handle)
430
431 END SUBROUTINE outer_loop_update_qs_env
432
433! **************************************************************************************************
434!> \brief uses the outer_scf_history to extrapolate new values for the variables
435!> and updates their value in qs_env accordingly
436!> \param qs_env the qs_environment_type where to update the variables
437!> \par History
438!> 03.2006 created [Joost VandeVondele]
439!> \note
440!> it assumes that the current value of qs_env still needs to be added to the history
441!> simple multilinear extrapolation is employed
442! **************************************************************************************************
443 SUBROUTINE outer_loop_extrapolate(qs_env)
444 TYPE(qs_environment_type), POINTER :: qs_env
445
446 CHARACTER(LEN=*), PARAMETER :: routinen = 'outer_loop_extrapolate'
447
448 INTEGER :: handle, ihis, ivec, n, nhistory, &
449 nvariables, nvec, outer_scf_ihistory
450 LOGICAL :: is_constraint
451 REAL(kind=dp) :: alpha
452 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: extrapolation
453 REAL(kind=dp), DIMENSION(:, :), POINTER :: outer_scf_history
454 TYPE(ddapc_restraint_type), POINTER :: ddapc_restraint_control
455 TYPE(dft_control_type), POINTER :: dft_control
456 TYPE(scf_control_type), POINTER :: scf_control
457
458 CALL timeset(routinen, handle)
459
460 CALL get_qs_env(qs_env, outer_scf_history=outer_scf_history, &
461 outer_scf_ihistory=outer_scf_ihistory, &
462 scf_control=scf_control, dft_control=dft_control)
463
464 nvariables = SIZE(outer_scf_history, 1)
465 nhistory = SIZE(outer_scf_history, 2)
466 ALLOCATE (extrapolation(nvariables))
467 cpassert(nhistory > 0)
468
469 ! add the current version of qs_env to the history
470 outer_scf_ihistory = outer_scf_ihistory + 1
471 ivec = 1 + modulo(outer_scf_ihistory - 1, nhistory)
472 SELECT CASE (scf_control%outer_scf%type)
473 CASE (outer_scf_none)
474 outer_scf_history(1, ivec) = 0.0_dp
475 CASE (outer_scf_ddapc_constraint)
476 DO n = 1, SIZE(dft_control%qs_control%ddapc_restraint_control)
477 NULLIFY (ddapc_restraint_control)
478 ddapc_restraint_control => dft_control%qs_control%ddapc_restraint_control(n)
479 is_constraint = (ddapc_restraint_control%functional_form == do_ddapc_constraint)
480 IF (is_constraint) EXIT
481 END DO
482 outer_scf_history(1, ivec) = &
483 ddapc_restraint_control%strength
484 CASE (outer_scf_s2_constraint)
485 outer_scf_history(1, ivec) = &
486 dft_control%qs_control%s2_restraint_control%strength
487 CASE (outer_scf_cdft_constraint)
488 outer_scf_history(:, ivec) = &
489 dft_control%qs_control%cdft_control%strength(:)
490 CASE (outer_scf_basis_center_opt)
491 outer_scf_history(1, ivec) = 0.0_dp
492 CASE DEFAULT
493 cpabort("Unknown outer SCF type")
494 END SELECT
495 CALL set_qs_env(qs_env, outer_scf_ihistory=outer_scf_ihistory)
496 ! multilinear extrapolation
497 nvec = min(nhistory, outer_scf_ihistory)
498 alpha = nvec
499 ivec = 1 + modulo(outer_scf_ihistory - 1, nhistory)
500 extrapolation(:) = alpha*outer_scf_history(:, ivec)
501 DO ihis = 2, nvec
502 alpha = -1.0_dp*alpha*real(nvec - ihis + 1, dp)/real(ihis, dp)
503 ivec = 1 + modulo(outer_scf_ihistory - ihis, nhistory)
504 extrapolation(:) = extrapolation + alpha*outer_scf_history(:, ivec)
505 END DO
506
507 ! update qs_env to use this extrapolation
508 SELECT CASE (scf_control%outer_scf%type)
509 CASE (outer_scf_none)
510 ! nothing
511 CASE (outer_scf_ddapc_constraint)
512 ddapc_restraint_control%strength = extrapolation(1)
513 CASE (outer_scf_s2_constraint)
514 dft_control%qs_control%s2_restraint_control%strength = extrapolation(1)
515 CASE (outer_scf_cdft_constraint)
516 dft_control%qs_control%cdft_control%strength(:) = extrapolation(:)
517 CASE (outer_scf_basis_center_opt)
518 ! nothing to do
519 CASE DEFAULT
520 cpabort("Unknown outer SCF type")
521 END SELECT
522
523 DEALLOCATE (extrapolation)
524
525 CALL timestop(handle)
526
527 END SUBROUTINE outer_loop_extrapolate
528
529! **************************************************************************************************
530!> \brief switch between two outer_scf envs stored in cdft_control
531!> \param scf_env the scf_env where values need to be updated using cdft_control
532!> \param scf_control the scf_control where values need to be updated using cdft_control
533!> \param cdft_control container for the second outer_scf env
534!> \param dir determines what switching operation to perform
535!> \par History
536!> 12.2015 created [Nico Holmberg]
537! **************************************************************************************************
538
539 SUBROUTINE outer_loop_switch(scf_env, scf_control, cdft_control, dir)
540 TYPE(qs_scf_env_type), POINTER :: scf_env
541 TYPE(scf_control_type), POINTER :: scf_control
542 TYPE(cdft_control_type), POINTER :: cdft_control
543 INTEGER, INTENT(IN) :: dir
544
545 INTEGER :: nvariables
546
547 SELECT CASE (dir)
548 CASE (cdft2ot)
549 ! Constraint -> OT
550 ! Switch data in scf_control: first save values that might have changed
551 IF (ASSOCIATED(scf_control%outer_scf%cdft_opt_control)) THEN
552 cpassert(ASSOCIATED(cdft_control%constraint_control%cdft_opt_control))
553 CALL cdft_opt_type_copy(cdft_control%constraint_control%cdft_opt_control, &
554 scf_control%outer_scf%cdft_opt_control)
555 ! OT SCF does not need cdft_opt_control
556 CALL cdft_opt_type_release(scf_control%outer_scf%cdft_opt_control)
557 END IF
558 ! Now switch
559 scf_control%outer_scf%have_scf = cdft_control%ot_control%have_scf
560 scf_control%outer_scf%max_scf = cdft_control%ot_control%max_scf
561 scf_control%outer_scf%eps_scf = cdft_control%ot_control%eps_scf
562 scf_control%outer_scf%step_size = cdft_control%ot_control%step_size
563 scf_control%outer_scf%type = cdft_control%ot_control%type
564 scf_control%outer_scf%optimizer = cdft_control%ot_control%optimizer
565 scf_control%outer_scf%diis_buffer_length = cdft_control%ot_control%diis_buffer_length
566 scf_control%outer_scf%bisect_trust_count = cdft_control%ot_control%bisect_trust_count
567 ! Switch data in scf_env: first save current values for constraint
568 cdft_control%constraint%iter_count = scf_env%outer_scf%iter_count
569 cdft_control%constraint%energy = scf_env%outer_scf%energy
570 cdft_control%constraint%variables = scf_env%outer_scf%variables
571 cdft_control%constraint%gradient = scf_env%outer_scf%gradient
572 cdft_control%constraint%count = scf_env%outer_scf%count
573 cdft_control%constraint%deallocate_jacobian = scf_env%outer_scf%deallocate_jacobian
574 IF (ASSOCIATED(scf_env%outer_scf%inv_jacobian)) THEN
575 nvariables = SIZE(scf_env%outer_scf%inv_jacobian, 1)
576 IF (.NOT. ASSOCIATED(cdft_control%constraint%inv_jacobian)) THEN
577 ALLOCATE (cdft_control%constraint%inv_jacobian(nvariables, nvariables))
578 END IF
579 cdft_control%constraint%inv_jacobian = scf_env%outer_scf%inv_jacobian
580 END IF
581 ! Now switch
582 IF (ASSOCIATED(scf_env%outer_scf%energy)) THEN
583 DEALLOCATE (scf_env%outer_scf%energy)
584 END IF
585 ALLOCATE (scf_env%outer_scf%energy(scf_control%outer_scf%max_scf + 1))
586 scf_env%outer_scf%energy = 0.0_dp
587 IF (ASSOCIATED(scf_env%outer_scf%variables)) THEN
588 DEALLOCATE (scf_env%outer_scf%variables)
589 END IF
590 ALLOCATE (scf_env%outer_scf%variables(1, scf_control%outer_scf%max_scf + 1))
591 scf_env%outer_scf%variables = 0.0_dp
592 IF (ASSOCIATED(scf_env%outer_scf%gradient)) THEN
593 DEALLOCATE (scf_env%outer_scf%gradient)
594 END IF
595 ALLOCATE (scf_env%outer_scf%gradient(1, scf_control%outer_scf%max_scf + 1))
596 scf_env%outer_scf%gradient = 0.0_dp
597 IF (ASSOCIATED(scf_env%outer_scf%count)) THEN
598 DEALLOCATE (scf_env%outer_scf%count)
599 END IF
600 ALLOCATE (scf_env%outer_scf%count(scf_control%outer_scf%max_scf + 1))
601 scf_env%outer_scf%count = 0
602 ! OT SCF does not need Jacobian
603 scf_env%outer_scf%deallocate_jacobian = .true.
604 IF (ASSOCIATED(scf_env%outer_scf%inv_jacobian)) THEN
605 DEALLOCATE (scf_env%outer_scf%inv_jacobian)
606 END IF
607 CASE (ot2cdft)
608 ! OT -> constraint
609 scf_control%outer_scf%have_scf = cdft_control%constraint_control%have_scf
610 scf_control%outer_scf%max_scf = cdft_control%constraint_control%max_scf
611 scf_control%outer_scf%eps_scf = cdft_control%constraint_control%eps_scf
612 scf_control%outer_scf%step_size = cdft_control%constraint_control%step_size
613 scf_control%outer_scf%type = cdft_control%constraint_control%type
614 scf_control%outer_scf%optimizer = cdft_control%constraint_control%optimizer
615 scf_control%outer_scf%diis_buffer_length = cdft_control%constraint_control%diis_buffer_length
616 scf_control%outer_scf%bisect_trust_count = cdft_control%constraint_control%bisect_trust_count
617 CALL cdft_opt_type_copy(scf_control%outer_scf%cdft_opt_control, &
618 cdft_control%constraint_control%cdft_opt_control)
619 nvariables = SIZE(cdft_control%constraint%variables, 1)
620 IF (ASSOCIATED(scf_env%outer_scf%energy)) THEN
621 DEALLOCATE (scf_env%outer_scf%energy)
622 END IF
623 ALLOCATE (scf_env%outer_scf%energy(scf_control%outer_scf%max_scf + 1))
624 scf_env%outer_scf%energy = cdft_control%constraint%energy
625 IF (ASSOCIATED(scf_env%outer_scf%variables)) THEN
626 DEALLOCATE (scf_env%outer_scf%variables)
627 END IF
628 ALLOCATE (scf_env%outer_scf%variables(nvariables, scf_control%outer_scf%max_scf + 1))
629 scf_env%outer_scf%variables = cdft_control%constraint%variables
630 IF (ASSOCIATED(scf_env%outer_scf%gradient)) THEN
631 DEALLOCATE (scf_env%outer_scf%gradient)
632 END IF
633 ALLOCATE (scf_env%outer_scf%gradient(nvariables, scf_control%outer_scf%max_scf + 1))
634 scf_env%outer_scf%gradient = cdft_control%constraint%gradient
635 IF (ASSOCIATED(scf_env%outer_scf%count)) THEN
636 DEALLOCATE (scf_env%outer_scf%count)
637 END IF
638 ALLOCATE (scf_env%outer_scf%count(scf_control%outer_scf%max_scf + 1))
639 scf_env%outer_scf%count = cdft_control%constraint%count
640 scf_env%outer_scf%iter_count = cdft_control%constraint%iter_count
641 scf_env%outer_scf%deallocate_jacobian = cdft_control%constraint%deallocate_jacobian
642 IF (ASSOCIATED(cdft_control%constraint%inv_jacobian)) THEN
643 IF (ASSOCIATED(scf_env%outer_scf%inv_jacobian)) THEN
644 DEALLOCATE (scf_env%outer_scf%inv_jacobian)
645 END IF
646 ALLOCATE (scf_env%outer_scf%inv_jacobian(nvariables, nvariables))
647 scf_env%outer_scf%inv_jacobian = cdft_control%constraint%inv_jacobian
648 END IF
649 CASE DEFAULT
650 cpabort("Switching direction should be either ot2cdft or cdft2ot")
651 END SELECT
652
653 END SUBROUTINE outer_loop_switch
654
655! **************************************************************************************************
656!> \brief purges outer_scf_history zeroing everything except
657!> the latest value of the outer_scf variable stored in qs_control
658!> \param qs_env the qs_environment_type where to purge
659!> \par History
660!> 05.2016 created [Nico Holmberg]
661! **************************************************************************************************
662 SUBROUTINE outer_loop_purge_history(qs_env)
663 TYPE(qs_environment_type), POINTER :: qs_env
664
665 CHARACTER(LEN=*), PARAMETER :: routinen = 'outer_loop_purge_history'
666
667 INTEGER :: handle, outer_scf_ihistory
668 REAL(kind=dp), DIMENSION(:, :), POINTER :: gradient_history, outer_scf_history, &
669 variable_history
670
671 CALL timeset(routinen, handle)
672
673 CALL get_qs_env(qs_env, outer_scf_history=outer_scf_history, &
674 outer_scf_ihistory=outer_scf_ihistory, &
675 gradient_history=gradient_history, &
676 variable_history=variable_history)
677 cpassert(SIZE(outer_scf_history, 2) > 0)
678 outer_scf_ihistory = 0
679 outer_scf_history = 0.0_dp
680 gradient_history = 0.0_dp
681 variable_history = 0.0_dp
682 CALL set_qs_env(qs_env, outer_scf_ihistory=outer_scf_ihistory)
683
684 CALL timestop(handle)
685
686 END SUBROUTINE outer_loop_purge_history
687
688END MODULE qs_outer_scf
modulo
static GRID_HOST_DEVICE int modulo(int a, int m)
Equivalent of Fortran's MODULO, which always return a positive number. https://gcc....
Definition grid_common.h:125
cp_log_handling::cp_to_string
Definition cp_log_handling.F:90
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
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:18
input_constants::outer_scf_optimizer_sd
integer, parameter, public outer_scf_optimizer_sd
Definition input_constants.F:763
input_constants::do_s2_constraint
integer, parameter, public do_s2_constraint
Definition input_constants.F:790
input_constants::broyden_type_2_explicit_ls
integer, parameter, public broyden_type_2_explicit_ls
Definition input_constants.F:773
input_constants::outer_scf_optimizer_bisect
integer, parameter, public outer_scf_optimizer_bisect
Definition input_constants.F:763
input_constants::outer_scf_optimizer_secant
integer, parameter, public outer_scf_optimizer_secant
Definition input_constants.F:763
input_constants::outer_scf_cdft_constraint
integer, parameter, public outer_scf_cdft_constraint
Definition input_constants.F:754
input_constants::do_ddapc_constraint
integer, parameter, public do_ddapc_constraint
Definition input_constants.F:750
input_constants::broyden_type_1_explicit
integer, parameter, public broyden_type_1_explicit
Definition input_constants.F:773
input_constants::broyden_type_2_ls
integer, parameter, public broyden_type_2_ls
Definition input_constants.F:773
input_constants::broyden_type_1
integer, parameter, public broyden_type_1
Definition input_constants.F:773
input_constants::outer_scf_optimizer_broyden
integer, parameter, public outer_scf_optimizer_broyden
Definition input_constants.F:763
input_constants::broyden_type_1_explicit_ls
integer, parameter, public broyden_type_1_explicit_ls
Definition input_constants.F:773
input_constants::outer_scf_basis_center_opt
integer, parameter, public outer_scf_basis_center_opt
Definition input_constants.F:754
input_constants::broyden_type_2_explicit
integer, parameter, public broyden_type_2_explicit
Definition input_constants.F:773
input_constants::outer_scf_s2_constraint
integer, parameter, public outer_scf_s2_constraint
Definition input_constants.F:754
input_constants::broyden_type_2
integer, parameter, public broyden_type_2
Definition input_constants.F:773
input_constants::cdft2ot
integer, parameter, public cdft2ot
Definition input_constants.F:798
input_constants::outer_scf_ddapc_constraint
integer, parameter, public outer_scf_ddapc_constraint
Definition input_constants.F:754
input_constants::outer_scf_optimizer_newton_ls
integer, parameter, public outer_scf_optimizer_newton_ls
Definition input_constants.F:763
input_constants::outer_scf_optimizer_none
integer, parameter, public outer_scf_optimizer_none
Definition input_constants.F:763
input_constants::outer_scf_optimizer_newton
integer, parameter, public outer_scf_optimizer_newton
Definition input_constants.F:763
input_constants::outer_scf_none
integer, parameter, public outer_scf_none
Definition input_constants.F:754
input_constants::broyden_type_1_ls
integer, parameter, public broyden_type_1_ls
Definition input_constants.F:773
input_constants::outer_scf_optimizer_diis
integer, parameter, public outer_scf_optimizer_diis
Definition input_constants.F:763
input_constants::ot2cdft
integer, parameter, public ot2cdft
Definition input_constants.F:798
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
mathlib
Collection of simple mathematical functions and subroutines.
Definition mathlib.F:15
mathlib::diamat_all
subroutine, public diamat_all(a, eigval, dac)
Diagonalize the symmetric n by n matrix a using the LAPACK library. Only the upper triangle of matrix...
Definition mathlib.F:381
qs_basis_gradient
...
Definition qs_basis_gradient.F:12
qs_basis_gradient::qs_basis_center_gradient
subroutine, public qs_basis_center_gradient(qs_env)
...
Definition qs_basis_gradient.F:84
qs_basis_gradient::return_basis_center_gradient_norm
real(kind=dp) function, public return_basis_center_gradient_norm(qs_env)
... returns the norm of the gradient vector, taking only floating components into account
Definition qs_basis_gradient.F:198
qs_basis_gradient::qs_update_basis_center_pos
subroutine, public qs_update_basis_center_pos(qs_env)
move atoms with kind float according to gradient
Definition qs_basis_gradient.F:233
qs_cdft_opt_types
Control parameters for optimizers that work with CDFT constraints.
Definition qs_cdft_opt_types.F:14
qs_cdft_opt_types::cdft_opt_type_release
subroutine, public cdft_opt_type_release(cdft_opt_control)
releases the CDFT optimizer control object
Definition qs_cdft_opt_types.F:134
qs_cdft_opt_types::cdft_opt_type_copy
subroutine, public cdft_opt_type_copy(new, old)
copies settings between two CDFT optimizer control objects retaining both
Definition qs_cdft_opt_types.F:307
qs_cdft_types
Defines CDFT control structures.
Definition qs_cdft_types.F:14
qs_energy_types
Definition qs_energy_types.F:14
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_environment_types::set_qs_env
subroutine, public set_qs_env(qs_env, super_cell, mos, qmmm, qmmm_periodic, mimic, ewald_env, ewald_pw, mpools, rho_external, external_vxc, mask, scf_control, rel_control, qs_charges, ks_env, ks_qmmm_env, wf_history, scf_env, active_space, input, oce, rho_atom_set, rho0_atom_set, rho0_mpole, run_rtp, rtp, rhoz_set, rhoz_tot, ecoul_1c, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, efield, rhoz_cneo_set, linres_control, xas_env, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, ls_scf_env, do_transport, transport_env, lri_env, lri_density, exstate_env, ec_env, dispersion_env, harris_env, gcp_env, mp2_env, bs_env, kg_env, force, kpoints, wanniercentres, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs, do_rixs, tb_tblite)
Set the QUICKSTEP environment.
Definition qs_environment_types.F:1120
qs_outer_scf
Routines for performing an outer scf loop.
Definition qs_outer_scf.F:14
qs_outer_scf::outer_loop_purge_history
subroutine, public outer_loop_purge_history(qs_env)
purges outer_scf_history zeroing everything except the latest value of the outer_scf variable stored ...
Definition qs_outer_scf.F:663
qs_outer_scf::outer_loop_switch
subroutine, public outer_loop_switch(scf_env, scf_control, cdft_control, dir)
switch between two outer_scf envs stored in cdft_control
Definition qs_outer_scf.F:540
qs_outer_scf::outer_loop_optimize
subroutine, public outer_loop_optimize(scf_env, scf_control)
optimizes the parameters of the outer_scf
Definition qs_outer_scf.F:187
qs_outer_scf::outer_loop_extrapolate
subroutine, public outer_loop_extrapolate(qs_env)
uses the outer_scf_history to extrapolate new values for the variables and updates their value in qs_...
Definition qs_outer_scf.F:444
qs_outer_scf::outer_loop_update_qs_env
subroutine, public outer_loop_update_qs_env(qs_env, scf_env)
propagates the updated variables to wherever they need to be set in qs_env
Definition qs_outer_scf.F:388
qs_outer_scf::outer_loop_variables_count
integer function, public outer_loop_variables_count(scf_control, cdft_control)
returns the number of variables that is employed in the outer loop. with a CDFT constraint this value...
Definition qs_outer_scf.F:69
qs_outer_scf::outer_loop_gradient
subroutine, public outer_loop_gradient(qs_env, scf_env)
computes the gradient wrt to the outer loop variables
Definition qs_outer_scf.F:103
qs_scf_types
module that contains the definitions of the scf types
Definition qs_scf_types.F:14
scf_control_types
parameters that control an scf iteration
Definition scf_control_types.F:17
cp_control_types::ddapc_restraint_type
Definition cp_control_types.F:432
cp_control_types::dft_control_type
Definition cp_control_types.F:757
cp_control_types::s2_restraint_type
Definition cp_control_types.F:447
qs_cdft_types::cdft_control_type
Definition qs_cdft_types.F:219
qs_energy_types::qs_energy_type
Definition qs_energy_types.F:25
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:222
qs_scf_types::qs_scf_env_type
Definition qs_scf_types.F:97
scf_control_types::scf_control_type
Definition scf_control_types.F:126