(git:9754b87)
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-2025 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("")
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) .GT. scf_control%outer_scf%bisect_trust_count) cycle
221 IF (scf_env%outer_scf%count(j) .GT. 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)) .LT. 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) .LT. 3) &
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 nvar = SIZE(scf_env%outer_scf%gradient, 1)
326 ALLOCATE (f(nvar, 1), x(nvar, 1))
327 DO i = 1, nvar
328 f(i, 1) = scf_env%outer_scf%gradient(i, ihistory) - scf_env%outer_scf%gradient(i, ihistory - 1)
329 x(i, 1) = scf_env%outer_scf%variables(i, ihistory) - scf_env%outer_scf%variables(i, ihistory - 1)
330 END DO
331 SELECT CASE (scf_control%outer_scf%cdft_opt_control%broyden_type)
332 CASE (broyden_type_1, broyden_type_1_explicit, broyden_type_1_ls, broyden_type_1_explicit_ls)
333 ! Broyden's 1st method
334 ! Denote: dx_n = \delta x_n; df_n = \delta f_n
335 ! 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)
336 scale = sum(matmul(transpose(x), matmul(inv_jacobian, f)))
337 scale = 1.0_dp/scale
338 IF (scale < 1.0e-12_dp) scale = 1.0e-12_dp
339 inv_jacobian = inv_jacobian + scale*matmul((x - matmul(inv_jacobian, f)), &
340 matmul(transpose(x), inv_jacobian))
341 CASE (broyden_type_2, broyden_type_2_explicit, broyden_type_2_ls, broyden_type_2_explicit_ls)
342 ! Broyden's 2nd method
343 ! J_(n+1)^(-1) = J_n^(-1) + (dx_n - J_n^(-1)*df_n)*(df_n^T)/(||df_n||^2)
344 scale = sum(matmul(transpose(f), f))
345 scale = 1.0_dp/scale
346 IF (scale < 1.0e-12_dp) scale = 1.0e-12_dp
347 inv_jacobian = inv_jacobian + scale*matmul((x - matmul(inv_jacobian, f)), transpose(inv_jacobian))
348 CASE DEFAULT
349 CALL cp_abort(__location__, &
350 "Unknown Broyden type: "// &
351 cp_to_string(scf_control%outer_scf%cdft_opt_control%broyden_type))
352 END SELECT
353 ! Clean up
354 DEALLOCATE (f, x)
355 END IF
356 ! Update variables x_(n+1) = x_n - J^(-1)*f(x_n)
357 scf_env%outer_scf%variables(:, ihistory + 1) = scf_env%outer_scf%variables(:, ihistory) - &
358 scf_control%outer_scf%cdft_opt_control%newton_step* &
359 matmul(inv_jacobian, scf_env%outer_scf%gradient(:, ihistory))
360 scf_control%outer_scf%cdft_opt_control%broyden_update = .true.
361 CASE (outer_scf_optimizer_newton, outer_scf_optimizer_newton_ls)
362 cpassert(ASSOCIATED(scf_env%outer_scf%inv_jacobian))
363 inv_jacobian => scf_env%outer_scf%inv_jacobian
364 scf_env%outer_scf%variables(:, ihistory + 1) = scf_env%outer_scf%variables(:, ihistory) - &
365 scf_control%outer_scf%cdft_opt_control%newton_step* &
366 matmul(inv_jacobian, scf_env%outer_scf%gradient(:, ihistory))
367 CASE DEFAULT
368 cpabort("")
369 END SELECT
370 EXIT
371 END DO
372 END IF
373
374 CALL timestop(handle)
375
376 END SUBROUTINE outer_loop_optimize
377
378! **************************************************************************************************
379!> \brief propagates the updated variables to wherever they need to be set in
380!> qs_env
381!> \param qs_env ...
382!> \param scf_env ...
383!> \par History
384!> 03.2006 created [Joost VandeVondele]
385! **************************************************************************************************
386 SUBROUTINE outer_loop_update_qs_env(qs_env, scf_env)
387 TYPE(qs_environment_type), POINTER :: qs_env
388 TYPE(qs_scf_env_type), POINTER :: scf_env
389
390 CHARACTER(LEN=*), PARAMETER :: routinen = 'outer_loop_update_qs_env'
391
392 INTEGER :: handle, ihistory, n
393 LOGICAL :: is_constraint
394 TYPE(cdft_control_type), POINTER :: cdft_control
395 TYPE(ddapc_restraint_type), POINTER :: ddapc_restraint_control
396 TYPE(dft_control_type), POINTER :: dft_control
397 TYPE(s2_restraint_type), POINTER :: s2_restraint_control
398 TYPE(scf_control_type), POINTER :: scf_control
399
400 CALL timeset(routinen, handle)
401
402 CALL get_qs_env(qs_env=qs_env, scf_control=scf_control, dft_control=dft_control)
403 ihistory = scf_env%outer_scf%iter_count
404
405 SELECT CASE (scf_control%outer_scf%type)
406 CASE (outer_scf_none)
407 ! do nothing
408 CASE (outer_scf_ddapc_constraint)
409 DO n = 1, SIZE(dft_control%qs_control%ddapc_restraint_control)
410 NULLIFY (ddapc_restraint_control)
411 ddapc_restraint_control => dft_control%qs_control%ddapc_restraint_control(n)
412 is_constraint = (ddapc_restraint_control%functional_form == do_ddapc_constraint)
413 IF (is_constraint) EXIT
414 END DO
415 ddapc_restraint_control%strength = scf_env%outer_scf%variables(1, ihistory + 1)
416 CASE (outer_scf_s2_constraint)
417 s2_restraint_control => dft_control%qs_control%s2_restraint_control
418 s2_restraint_control%strength = scf_env%outer_scf%variables(1, ihistory + 1)
419 CASE (outer_scf_cdft_constraint)
420 cdft_control => dft_control%qs_control%cdft_control
421 cdft_control%strength(:) = scf_env%outer_scf%variables(:, ihistory + 1)
422 CASE (outer_scf_basis_center_opt)
423 CALL qs_update_basis_center_pos(qs_env)
424 CASE DEFAULT
425 cpabort("")
426 END SELECT
427
428 CALL timestop(handle)
429
430 END SUBROUTINE outer_loop_update_qs_env
431
432! **************************************************************************************************
433!> \brief uses the outer_scf_history to extrapolate new values for the variables
434!> and updates their value in qs_env accordingly
435!> \param qs_env the qs_environment_type where to update the variables
436!> \par History
437!> 03.2006 created [Joost VandeVondele]
438!> \note
439!> it assumes that the current value of qs_env still needs to be added to the history
440!> simple multilinear extrapolation is employed
441! **************************************************************************************************
442 SUBROUTINE outer_loop_extrapolate(qs_env)
443 TYPE(qs_environment_type), POINTER :: qs_env
444
445 CHARACTER(LEN=*), PARAMETER :: routinen = 'outer_loop_extrapolate'
446
447 INTEGER :: handle, ihis, ivec, n, nhistory, &
448 nvariables, nvec, outer_scf_ihistory
449 LOGICAL :: is_constraint
450 REAL(kind=dp) :: alpha
451 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: extrapolation
452 REAL(kind=dp), DIMENSION(:, :), POINTER :: outer_scf_history
453 TYPE(ddapc_restraint_type), POINTER :: ddapc_restraint_control
454 TYPE(dft_control_type), POINTER :: dft_control
455 TYPE(scf_control_type), POINTER :: scf_control
456
457 CALL timeset(routinen, handle)
458
459 CALL get_qs_env(qs_env, outer_scf_history=outer_scf_history, &
460 outer_scf_ihistory=outer_scf_ihistory, &
461 scf_control=scf_control, dft_control=dft_control)
462
463 nvariables = SIZE(outer_scf_history, 1)
464 nhistory = SIZE(outer_scf_history, 2)
465 ALLOCATE (extrapolation(nvariables))
466 cpassert(nhistory > 0)
467
468 ! add the current version of qs_env to the history
469 outer_scf_ihistory = outer_scf_ihistory + 1
470 ivec = 1 + modulo(outer_scf_ihistory - 1, nhistory)
471 SELECT CASE (scf_control%outer_scf%type)
472 CASE (outer_scf_none)
473 outer_scf_history(1, ivec) = 0.0_dp
474 CASE (outer_scf_ddapc_constraint)
475 DO n = 1, SIZE(dft_control%qs_control%ddapc_restraint_control)
476 NULLIFY (ddapc_restraint_control)
477 ddapc_restraint_control => dft_control%qs_control%ddapc_restraint_control(n)
478 is_constraint = (ddapc_restraint_control%functional_form == do_ddapc_constraint)
479 IF (is_constraint) EXIT
480 END DO
481 outer_scf_history(1, ivec) = &
482 ddapc_restraint_control%strength
483 CASE (outer_scf_s2_constraint)
484 outer_scf_history(1, ivec) = &
485 dft_control%qs_control%s2_restraint_control%strength
486 CASE (outer_scf_cdft_constraint)
487 outer_scf_history(:, ivec) = &
488 dft_control%qs_control%cdft_control%strength(:)
489 CASE (outer_scf_basis_center_opt)
490 outer_scf_history(1, ivec) = 0.0_dp
491 CASE DEFAULT
492 cpabort("")
493 END SELECT
494 CALL set_qs_env(qs_env, outer_scf_ihistory=outer_scf_ihistory)
495 ! multilinear extrapolation
496 nvec = min(nhistory, outer_scf_ihistory)
497 alpha = nvec
498 ivec = 1 + modulo(outer_scf_ihistory - 1, nhistory)
499 extrapolation(:) = alpha*outer_scf_history(:, ivec)
500 DO ihis = 2, nvec
501 alpha = -1.0_dp*alpha*real(nvec - ihis + 1, dp)/real(ihis, dp)
502 ivec = 1 + modulo(outer_scf_ihistory - ihis, nhistory)
503 extrapolation(:) = extrapolation + alpha*outer_scf_history(:, ivec)
504 END DO
505
506 ! update qs_env to use this extrapolation
507 SELECT CASE (scf_control%outer_scf%type)
508 CASE (outer_scf_none)
509 ! nothing
510 CASE (outer_scf_ddapc_constraint)
511 ddapc_restraint_control%strength = extrapolation(1)
512 CASE (outer_scf_s2_constraint)
513 dft_control%qs_control%s2_restraint_control%strength = extrapolation(1)
514 CASE (outer_scf_cdft_constraint)
515 dft_control%qs_control%cdft_control%strength(:) = extrapolation(:)
516 CASE (outer_scf_basis_center_opt)
517 ! nothing to do
518 CASE DEFAULT
519 cpabort("")
520 END SELECT
521
522 DEALLOCATE (extrapolation)
523
524 CALL timestop(handle)
525
526 END SUBROUTINE outer_loop_extrapolate
527
528! **************************************************************************************************
529!> \brief switch between two outer_scf envs stored in cdft_control
530!> \param scf_env the scf_env where values need to be updated using cdft_control
531!> \param scf_control the scf_control where values need to be updated using cdft_control
532!> \param cdft_control container for the second outer_scf env
533!> \param dir determines what switching operation to perform
534!> \par History
535!> 12.2015 created [Nico Holmberg]
536! **************************************************************************************************
537
538 SUBROUTINE outer_loop_switch(scf_env, scf_control, cdft_control, dir)
539 TYPE(qs_scf_env_type), POINTER :: scf_env
540 TYPE(scf_control_type), POINTER :: scf_control
541 TYPE(cdft_control_type), POINTER :: cdft_control
542 INTEGER, INTENT(IN) :: dir
543
544 INTEGER :: nvariables
545
546 SELECT CASE (dir)
547 CASE (cdft2ot)
548 ! Constraint -> OT
549 ! Switch data in scf_control: first save values that might have changed
550 IF (ASSOCIATED(scf_control%outer_scf%cdft_opt_control)) THEN
551 cpassert(ASSOCIATED(cdft_control%constraint_control%cdft_opt_control))
552 CALL cdft_opt_type_copy(cdft_control%constraint_control%cdft_opt_control, &
553 scf_control%outer_scf%cdft_opt_control)
554 ! OT SCF does not need cdft_opt_control
555 CALL cdft_opt_type_release(scf_control%outer_scf%cdft_opt_control)
556 END IF
557 ! Now switch
558 scf_control%outer_scf%have_scf = cdft_control%ot_control%have_scf
559 scf_control%outer_scf%max_scf = cdft_control%ot_control%max_scf
560 scf_control%outer_scf%eps_scf = cdft_control%ot_control%eps_scf
561 scf_control%outer_scf%step_size = cdft_control%ot_control%step_size
562 scf_control%outer_scf%type = cdft_control%ot_control%type
563 scf_control%outer_scf%optimizer = cdft_control%ot_control%optimizer
564 scf_control%outer_scf%diis_buffer_length = cdft_control%ot_control%diis_buffer_length
565 scf_control%outer_scf%bisect_trust_count = cdft_control%ot_control%bisect_trust_count
566 ! Switch data in scf_env: first save current values for constraint
567 cdft_control%constraint%iter_count = scf_env%outer_scf%iter_count
568 cdft_control%constraint%energy = scf_env%outer_scf%energy
569 cdft_control%constraint%variables = scf_env%outer_scf%variables
570 cdft_control%constraint%gradient = scf_env%outer_scf%gradient
571 cdft_control%constraint%count = scf_env%outer_scf%count
572 cdft_control%constraint%deallocate_jacobian = scf_env%outer_scf%deallocate_jacobian
573 IF (ASSOCIATED(scf_env%outer_scf%inv_jacobian)) THEN
574 nvariables = SIZE(scf_env%outer_scf%inv_jacobian, 1)
575 IF (.NOT. ASSOCIATED(cdft_control%constraint%inv_jacobian)) &
576 ALLOCATE (cdft_control%constraint%inv_jacobian(nvariables, nvariables))
577 cdft_control%constraint%inv_jacobian = scf_env%outer_scf%inv_jacobian
578 END IF
579 ! Now switch
580 IF (ASSOCIATED(scf_env%outer_scf%energy)) &
581 DEALLOCATE (scf_env%outer_scf%energy)
582 ALLOCATE (scf_env%outer_scf%energy(scf_control%outer_scf%max_scf + 1))
583 scf_env%outer_scf%energy = 0.0_dp
584 IF (ASSOCIATED(scf_env%outer_scf%variables)) &
585 DEALLOCATE (scf_env%outer_scf%variables)
586 ALLOCATE (scf_env%outer_scf%variables(1, scf_control%outer_scf%max_scf + 1))
587 scf_env%outer_scf%variables = 0.0_dp
588 IF (ASSOCIATED(scf_env%outer_scf%gradient)) &
589 DEALLOCATE (scf_env%outer_scf%gradient)
590 ALLOCATE (scf_env%outer_scf%gradient(1, scf_control%outer_scf%max_scf + 1))
591 scf_env%outer_scf%gradient = 0.0_dp
592 IF (ASSOCIATED(scf_env%outer_scf%count)) &
593 DEALLOCATE (scf_env%outer_scf%count)
594 ALLOCATE (scf_env%outer_scf%count(scf_control%outer_scf%max_scf + 1))
595 scf_env%outer_scf%count = 0
596 ! OT SCF does not need Jacobian
597 scf_env%outer_scf%deallocate_jacobian = .true.
598 IF (ASSOCIATED(scf_env%outer_scf%inv_jacobian)) &
599 DEALLOCATE (scf_env%outer_scf%inv_jacobian)
600 CASE (ot2cdft)
601 ! OT -> constraint
602 scf_control%outer_scf%have_scf = cdft_control%constraint_control%have_scf
603 scf_control%outer_scf%max_scf = cdft_control%constraint_control%max_scf
604 scf_control%outer_scf%eps_scf = cdft_control%constraint_control%eps_scf
605 scf_control%outer_scf%step_size = cdft_control%constraint_control%step_size
606 scf_control%outer_scf%type = cdft_control%constraint_control%type
607 scf_control%outer_scf%optimizer = cdft_control%constraint_control%optimizer
608 scf_control%outer_scf%diis_buffer_length = cdft_control%constraint_control%diis_buffer_length
609 scf_control%outer_scf%bisect_trust_count = cdft_control%constraint_control%bisect_trust_count
610 CALL cdft_opt_type_copy(scf_control%outer_scf%cdft_opt_control, &
611 cdft_control%constraint_control%cdft_opt_control)
612 nvariables = SIZE(cdft_control%constraint%variables, 1)
613 IF (ASSOCIATED(scf_env%outer_scf%energy)) &
614 DEALLOCATE (scf_env%outer_scf%energy)
615 ALLOCATE (scf_env%outer_scf%energy(scf_control%outer_scf%max_scf + 1))
616 scf_env%outer_scf%energy = cdft_control%constraint%energy
617 IF (ASSOCIATED(scf_env%outer_scf%variables)) &
618 DEALLOCATE (scf_env%outer_scf%variables)
619 ALLOCATE (scf_env%outer_scf%variables(nvariables, scf_control%outer_scf%max_scf + 1))
620 scf_env%outer_scf%variables = cdft_control%constraint%variables
621 IF (ASSOCIATED(scf_env%outer_scf%gradient)) &
622 DEALLOCATE (scf_env%outer_scf%gradient)
623 ALLOCATE (scf_env%outer_scf%gradient(nvariables, scf_control%outer_scf%max_scf + 1))
624 scf_env%outer_scf%gradient = cdft_control%constraint%gradient
625 IF (ASSOCIATED(scf_env%outer_scf%count)) &
626 DEALLOCATE (scf_env%outer_scf%count)
627 ALLOCATE (scf_env%outer_scf%count(scf_control%outer_scf%max_scf + 1))
628 scf_env%outer_scf%count = cdft_control%constraint%count
629 scf_env%outer_scf%iter_count = cdft_control%constraint%iter_count
630 scf_env%outer_scf%deallocate_jacobian = cdft_control%constraint%deallocate_jacobian
631 IF (ASSOCIATED(cdft_control%constraint%inv_jacobian)) THEN
632 IF (ASSOCIATED(scf_env%outer_scf%inv_jacobian)) &
633 DEALLOCATE (scf_env%outer_scf%inv_jacobian)
634 ALLOCATE (scf_env%outer_scf%inv_jacobian(nvariables, nvariables))
635 scf_env%outer_scf%inv_jacobian = cdft_control%constraint%inv_jacobian
636 END IF
637 CASE DEFAULT
638 cpabort("")
639 END SELECT
640
641 END SUBROUTINE outer_loop_switch
642
643! **************************************************************************************************
644!> \brief purges outer_scf_history zeroing everything except
645!> the latest value of the outer_scf variable stored in qs_control
646!> \param qs_env the qs_environment_type where to purge
647!> \par History
648!> 05.2016 created [Nico Holmberg]
649! **************************************************************************************************
650 SUBROUTINE outer_loop_purge_history(qs_env)
651 TYPE(qs_environment_type), POINTER :: qs_env
652
653 CHARACTER(LEN=*), PARAMETER :: routinen = 'outer_loop_purge_history'
654
655 INTEGER :: handle, outer_scf_ihistory
656 REAL(kind=dp), DIMENSION(:, :), POINTER :: gradient_history, outer_scf_history, &
657 variable_history
658
659 CALL timeset(routinen, handle)
660
661 CALL get_qs_env(qs_env, outer_scf_history=outer_scf_history, &
662 outer_scf_ihistory=outer_scf_ihistory, &
663 gradient_history=gradient_history, &
664 variable_history=variable_history)
665 cpassert(SIZE(outer_scf_history, 2) > 0)
666 outer_scf_ihistory = 0
667 outer_scf_history = 0.0_dp
668 gradient_history = 0.0_dp
669 variable_history = 0.0_dp
670 CALL set_qs_env(qs_env, outer_scf_ihistory=outer_scf_ihistory)
671
672 CALL timestop(handle)
673
674 END SUBROUTINE outer_loop_purge_history
675
676END 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:120
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:17
input_constants::outer_scf_optimizer_sd
integer, parameter, public outer_scf_optimizer_sd
Definition input_constants.F:733
input_constants::do_s2_constraint
integer, parameter, public do_s2_constraint
Definition input_constants.F:760
input_constants::broyden_type_2_explicit_ls
integer, parameter, public broyden_type_2_explicit_ls
Definition input_constants.F:743
input_constants::outer_scf_optimizer_bisect
integer, parameter, public outer_scf_optimizer_bisect
Definition input_constants.F:733
input_constants::outer_scf_optimizer_secant
integer, parameter, public outer_scf_optimizer_secant
Definition input_constants.F:733
input_constants::outer_scf_cdft_constraint
integer, parameter, public outer_scf_cdft_constraint
Definition input_constants.F:724
input_constants::do_ddapc_constraint
integer, parameter, public do_ddapc_constraint
Definition input_constants.F:720
input_constants::broyden_type_1_explicit
integer, parameter, public broyden_type_1_explicit
Definition input_constants.F:743
input_constants::broyden_type_2_ls
integer, parameter, public broyden_type_2_ls
Definition input_constants.F:743
input_constants::broyden_type_1
integer, parameter, public broyden_type_1
Definition input_constants.F:743
input_constants::outer_scf_optimizer_broyden
integer, parameter, public outer_scf_optimizer_broyden
Definition input_constants.F:733
input_constants::broyden_type_1_explicit_ls
integer, parameter, public broyden_type_1_explicit_ls
Definition input_constants.F:743
input_constants::outer_scf_basis_center_opt
integer, parameter, public outer_scf_basis_center_opt
Definition input_constants.F:724
input_constants::broyden_type_2_explicit
integer, parameter, public broyden_type_2_explicit
Definition input_constants.F:743
input_constants::outer_scf_s2_constraint
integer, parameter, public outer_scf_s2_constraint
Definition input_constants.F:724
input_constants::broyden_type_2
integer, parameter, public broyden_type_2
Definition input_constants.F:743
input_constants::cdft2ot
integer, parameter, public cdft2ot
Definition input_constants.F:768
input_constants::outer_scf_ddapc_constraint
integer, parameter, public outer_scf_ddapc_constraint
Definition input_constants.F:724
input_constants::outer_scf_optimizer_newton_ls
integer, parameter, public outer_scf_optimizer_newton_ls
Definition input_constants.F:733
input_constants::outer_scf_optimizer_none
integer, parameter, public outer_scf_optimizer_none
Definition input_constants.F:733
input_constants::outer_scf_optimizer_newton
integer, parameter, public outer_scf_optimizer_newton
Definition input_constants.F:733
input_constants::outer_scf_none
integer, parameter, public outer_scf_none
Definition input_constants.F:724
input_constants::broyden_type_1_ls
integer, parameter, public broyden_type_1_ls
Definition input_constants.F:743
input_constants::outer_scf_optimizer_diis
integer, parameter, public outer_scf_optimizer_diis
Definition input_constants.F:733
input_constants::ot2cdft
integer, parameter, public ot2cdft
Definition input_constants.F:768
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:373
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:300
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, 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, 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, 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, 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, ecoul_1c, rho0_s_rs, rho0_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)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:514
qs_environment_types::set_qs_env
subroutine, public set_qs_env(qs_env, super_cell, mos, qmmm, qmmm_periodic, 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, 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)
Set the QUICKSTEP environment.
Definition qs_environment_types.F:1055
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:651
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:539
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:443
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:387
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:291
cp_control_types::dft_control_type
Definition cp_control_types.F:570
cp_control_types::s2_restraint_type
Definition cp_control_types.F:306
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:217
qs_scf_types::qs_scf_env_type
Definition qs_scf_types.F:97
scf_control_types::scf_control_type
Definition scf_control_types.F:124