Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Routines for performing an outer scf loop
10 : !> \par History
11 : !> Created [2006.03]
12 : !> \author Joost VandeVondele
13 : ! **************************************************************************************************
14 : MODULE 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 :
57 : CONTAINS
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 5552 : 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 5552 : SELECT CASE (scf_control%outer_scf%type)
75 : CASE (outer_scf_ddapc_constraint)
76 : res = 1
77 : CASE (outer_scf_s2_constraint)
78 62 : res = 1
79 : CASE (outer_scf_cdft_constraint)
80 62 : IF (PRESENT(cdft_control)) THEN
81 62 : 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 0 : res = 1
89 : CASE DEFAULT
90 5552 : res = 0
91 : END SELECT
92 :
93 5552 : 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 6241 : 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 6241 : CALL timeset(routineN, handle)
118 :
119 : CALL get_qs_env(qs_env=qs_env, scf_control=scf_control, &
120 6241 : dft_control=dft_control, energy=energy)
121 6241 : CPASSERT(scf_control%outer_scf%have_scf)
122 :
123 6241 : ihistory = scf_env%outer_scf%iter_count
124 6241 : CPASSERT(ihistory <= SIZE(scf_env%outer_scf%energy, 1))
125 :
126 6241 : scf_env%outer_scf%energy(ihistory) = energy%total
127 :
128 11724 : 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 5483 : scf_env%outer_scf%variables(1, ihistory) = scf_env%iter_delta
132 5483 : scf_env%outer_scf%gradient(1, ihistory) = scf_env%iter_delta
133 : CASE (outer_scf_ddapc_constraint)
134 76 : CPASSERT(dft_control%qs_control%ddapc_restraint)
135 76 : DO n = 1, SIZE(dft_control%qs_control%ddapc_restraint_control)
136 76 : NULLIFY (ddapc_restraint_control)
137 76 : ddapc_restraint_control => dft_control%qs_control%ddapc_restraint_control(n)
138 76 : is_constraint = (ddapc_restraint_control%functional_form == do_ddapc_constraint)
139 76 : IF (is_constraint) EXIT
140 : END DO
141 76 : CPASSERT(is_constraint)
142 :
143 152 : scf_env%outer_scf%variables(:, ihistory) = ddapc_restraint_control%strength
144 : scf_env%outer_scf%gradient(:, ihistory) = ddapc_restraint_control%ddapc_order_p - &
145 152 : ddapc_restraint_control%target
146 : CASE (outer_scf_s2_constraint)
147 0 : CPASSERT(dft_control%qs_control%s2_restraint)
148 0 : s2_restraint_control => dft_control%qs_control%s2_restraint_control
149 0 : is_constraint = (s2_restraint_control%functional_form == do_s2_constraint)
150 0 : CPASSERT(is_constraint)
151 :
152 0 : scf_env%outer_scf%variables(:, ihistory) = s2_restraint_control%strength
153 : scf_env%outer_scf%gradient(:, ihistory) = s2_restraint_control%s2_order_p - &
154 0 : s2_restraint_control%target
155 : CASE (outer_scf_cdft_constraint)
156 682 : CPASSERT(dft_control%qs_control%cdft)
157 682 : cdft_control => dft_control%qs_control%cdft_control
158 1498 : DO ivar = 1, SIZE(scf_env%outer_scf%gradient, 1)
159 816 : scf_env%outer_scf%variables(ivar, ihistory) = cdft_control%strength(ivar)
160 : scf_env%outer_scf%gradient(ivar, ihistory) = cdft_control%value(ivar) - &
161 1498 : cdft_control%target(ivar)
162 : END DO
163 : CASE (outer_scf_basis_center_opt)
164 0 : CALL qs_basis_center_gradient(qs_env)
165 0 : scf_env%outer_scf%gradient(:, ihistory) = return_basis_center_gradient_norm(qs_env)
166 :
167 : CASE DEFAULT
168 6241 : CPABORT("Unknown outer SCF type")
169 :
170 : END SELECT
171 :
172 6241 : CALL timestop(handle)
173 :
174 6241 : 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 1150 : 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 1150 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: ev
196 1150 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: a, b, f, x
197 1150 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: inv_jacobian
198 :
199 1150 : CALL timeset(routineN, handle)
200 :
201 1150 : ihistory = scf_env%outer_scf%iter_count
202 1150 : optimizer_type = scf_control%outer_scf%optimizer
203 1150 : NULLIFY (inv_jacobian)
204 :
205 1150 : IF (scf_control%outer_scf%type == outer_scf_basis_center_opt) THEN
206 0 : 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 44 : SELECT CASE (optimizer_type)
211 : CASE (outer_scf_optimizer_bisect) ! bisection on the gradient, needs to be 1D
212 44 : 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 44 : ilow = -1
215 44 : ihigh = -1
216 44 : interval = HUGE(interval)
217 100 : DO i = 1, ihistory
218 112 : DO j = i + 1, ihistory
219 : ! distrust often used points
220 12 : IF (scf_env%outer_scf%count(i) > scf_control%outer_scf%bisect_trust_count) CYCLE
221 12 : 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 12 : IF (scf_env%outer_scf%gradient(1, i)* &
225 56 : scf_env%outer_scf%gradient(1, j) < 0.0_dp) THEN
226 4 : tmp = ABS(scf_env%outer_scf%variables(1, i) - scf_env%outer_scf%variables(1, j))
227 4 : IF (tmp < interval) THEN
228 4 : ilow = i
229 4 : ihigh = j
230 4 : interval = tmp
231 : END IF
232 : END IF
233 : END DO
234 : END DO
235 44 : 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 4 : scf_env%outer_scf%count(ilow) = scf_env%outer_scf%count(ilow) + 1
240 4 : 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 8 : scf_env%outer_scf%variables(:, ihigh))
243 : CASE (outer_scf_optimizer_none)
244 1752 : 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 212 : scf_control%outer_scf%step_size*scf_env%outer_scf%gradient(:, ihistory)
251 : CASE (outer_scf_optimizer_diis)
252 124 : CPASSERT(scf_control%outer_scf%diis_buffer_length > 0)
253 : ! set up DIIS matrix
254 124 : nb = MIN(ihistory, scf_control%outer_scf%diis_buffer_length)
255 124 : IF (nb < 2) THEN
256 : optimizer_type = outer_scf_optimizer_sd
257 : CYCLE
258 : ELSE
259 416 : ALLOCATE (b(nb + 1, nb + 1), a(nb + 1, nb + 1), ev(nb + 1))
260 168 : DO I = 1, nb
261 360 : DO J = I, nb
262 192 : ibuf = ihistory - nb + i
263 192 : jbuf = ihistory - nb + j
264 : b(I, J) = DOT_PRODUCT(scf_env%outer_scf%gradient(:, ibuf), &
265 384 : scf_env%outer_scf%gradient(:, jbuf))
266 308 : b(J, I) = b(I, J)
267 : END DO
268 : END DO
269 220 : b(nb + 1, :) = -1.0_dp
270 220 : b(:, nb + 1) = -1.0_dp
271 52 : b(nb + 1, nb + 1) = 0.0_dp
272 :
273 52 : CALL diamat_all(b, ev)
274 772 : a(:, :) = b
275 220 : DO I = 1, nb + 1
276 220 : IF (ABS(ev(I)) < 1.0E-12_dp) THEN
277 60 : a(:, I) = 0.0_dp
278 : ELSE
279 660 : a(:, I) = a(:, I)/ev(I)
280 : END IF
281 : END DO
282 1044 : ev(:) = -MATMUL(a, b(nb + 1, :))
283 :
284 104 : scf_env%outer_scf%variables(:, ihistory + 1) = 0.0_dp
285 168 : DO i = 1, nb
286 116 : ibuf = ihistory - nb + i
287 : scf_env%outer_scf%variables(:, ihistory + 1) = scf_env%outer_scf%variables(:, ihistory + 1) + &
288 284 : ev(i)*scf_env%outer_scf%variables(:, ibuf)
289 : END DO
290 52 : DEALLOCATE (a, b, ev)
291 : END IF
292 : CASE (outer_scf_optimizer_secant)
293 4 : CPASSERT(SIZE(scf_env%outer_scf%gradient, 2) >= 3)
294 4 : CPASSERT(SIZE(scf_env%outer_scf%gradient, 1) == 1)
295 4 : nvar = SIZE(scf_env%outer_scf%gradient, 1)
296 4 : 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 2 : scf_env%outer_scf%gradient(1, ihistory)
308 : CASE (outer_scf_optimizer_broyden)
309 24 : IF (.NOT. ASSOCIATED(scf_env%outer_scf%inv_jacobian)) THEN
310 : ! Inverse Jacobian not yet built, switch to sd
311 106 : optimizer_type = outer_scf_optimizer_sd
312 : CYCLE
313 : END IF
314 16 : inv_jacobian => scf_env%outer_scf%inv_jacobian
315 16 : IF (ihistory < 2) THEN
316 : ! Cannot perform a Broyden update without enough SCF history on this energy evaluation
317 2 : scf_control%outer_scf%cdft_opt_control%broyden_update = .FALSE.
318 : END IF
319 16 : IF (scf_control%outer_scf%cdft_opt_control%broyden_update) THEN
320 : ! Perform a Broyden update of the inverse Jacobian J^(-1)
321 6 : IF (SIZE(scf_env%outer_scf%gradient, 2) < 3) THEN
322 : CALL cp_abort(__LOCATION__, &
323 : "Keyword EXTRAPOLATION_ORDER in section OUTER_SCF "// &
324 0 : "must be greater than or equal to 3 for Broyden optimizers.")
325 : END IF
326 6 : nvar = SIZE(scf_env%outer_scf%gradient, 1)
327 24 : ALLOCATE (f(nvar, 1), x(nvar, 1))
328 12 : DO i = 1, nvar
329 6 : f(i, 1) = scf_env%outer_scf%gradient(i, ihistory) - scf_env%outer_scf%gradient(i, ihistory - 1)
330 12 : x(i, 1) = scf_env%outer_scf%variables(i, ihistory) - scf_env%outer_scf%variables(i, ihistory - 1)
331 : END DO
332 10 : 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 64 : scale = SUM(MATMUL(TRANSPOSE(x), MATMUL(inv_jacobian, f)))
338 4 : scale = 1.0_dp/scale
339 : IF (scale < 1.0E-12_dp) scale = 1.0E-12_dp
340 28 : inv_jacobian = inv_jacobian + scale*MATMUL((x - MATMUL(inv_jacobian, f)), &
341 84 : 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 14 : scale = SUM(MATMUL(TRANSPOSE(f), f))
346 2 : scale = 1.0_dp/scale
347 : IF (scale < 1.0E-12_dp) scale = 1.0E-12_dp
348 36 : 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 6 : cp_to_string(scf_control%outer_scf%cdft_opt_control%broyden_type))
353 : END SELECT
354 : ! Clean up
355 6 : 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 128 : MATMUL(inv_jacobian, scf_env%outer_scf%gradient(:, ihistory))
361 16 : scf_control%outer_scf%cdft_opt_control%broyden_update = .TRUE.
362 : CASE (outer_scf_optimizer_newton, outer_scf_optimizer_newton_ls)
363 94 : CPASSERT(ASSOCIATED(scf_env%outer_scf%inv_jacobian))
364 94 : 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 1368 : MATMUL(inv_jacobian, scf_env%outer_scf%gradient(:, ihistory))
368 : CASE DEFAULT
369 1190 : CPABORT("Unknown outer SCF optimizer")
370 : END SELECT
371 : EXIT
372 : END DO
373 : END IF
374 :
375 1150 : CALL timestop(handle)
376 :
377 2300 : 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 1264 : 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 1264 : CALL timeset(routineN, handle)
402 :
403 1264 : CALL get_qs_env(qs_env=qs_env, scf_control=scf_control, dft_control=dft_control)
404 1264 : ihistory = scf_env%outer_scf%iter_count
405 :
406 1314 : SELECT CASE (scf_control%outer_scf%type)
407 : CASE (outer_scf_none)
408 : ! do nothing
409 : CASE (outer_scf_ddapc_constraint)
410 50 : DO n = 1, SIZE(dft_control%qs_control%ddapc_restraint_control)
411 50 : NULLIFY (ddapc_restraint_control)
412 50 : ddapc_restraint_control => dft_control%qs_control%ddapc_restraint_control(n)
413 50 : is_constraint = (ddapc_restraint_control%functional_form == do_ddapc_constraint)
414 50 : IF (is_constraint) EXIT
415 : END DO
416 50 : ddapc_restraint_control%strength = scf_env%outer_scf%variables(1, ihistory + 1)
417 : CASE (outer_scf_s2_constraint)
418 0 : s2_restraint_control => dft_control%qs_control%s2_restraint_control
419 0 : s2_restraint_control%strength = scf_env%outer_scf%variables(1, ihistory + 1)
420 : CASE (outer_scf_cdft_constraint)
421 338 : cdft_control => dft_control%qs_control%cdft_control
422 1560 : cdft_control%strength(:) = scf_env%outer_scf%variables(:, ihistory + 1)
423 : CASE (outer_scf_basis_center_opt)
424 0 : CALL qs_update_basis_center_pos(qs_env)
425 : CASE DEFAULT
426 1264 : CPABORT("Unknown outer SCF type")
427 : END SELECT
428 :
429 1264 : CALL timestop(handle)
430 :
431 1264 : 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 4301 : 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 4301 : REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: extrapolation
453 4301 : 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 4301 : 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 4301 : scf_control=scf_control, dft_control=dft_control)
463 :
464 4301 : nvariables = SIZE(outer_scf_history, 1)
465 4301 : nhistory = SIZE(outer_scf_history, 2)
466 12903 : ALLOCATE (extrapolation(nvariables))
467 4301 : CPASSERT(nhistory > 0)
468 :
469 : ! add the current version of qs_env to the history
470 4301 : outer_scf_ihistory = outer_scf_ihistory + 1
471 4301 : ivec = 1 + MODULO(outer_scf_ihistory - 1, nhistory)
472 8232 : SELECT CASE (scf_control%outer_scf%type)
473 : CASE (outer_scf_none)
474 3931 : outer_scf_history(1, ivec) = 0.0_dp
475 : CASE (outer_scf_ddapc_constraint)
476 26 : DO n = 1, SIZE(dft_control%qs_control%ddapc_restraint_control)
477 26 : NULLIFY (ddapc_restraint_control)
478 26 : ddapc_restraint_control => dft_control%qs_control%ddapc_restraint_control(n)
479 26 : is_constraint = (ddapc_restraint_control%functional_form == do_ddapc_constraint)
480 26 : IF (is_constraint) EXIT
481 : END DO
482 : outer_scf_history(1, ivec) = &
483 26 : ddapc_restraint_control%strength
484 : CASE (outer_scf_s2_constraint)
485 : outer_scf_history(1, ivec) = &
486 0 : dft_control%qs_control%s2_restraint_control%strength
487 : CASE (outer_scf_cdft_constraint)
488 : outer_scf_history(:, ivec) = &
489 1436 : dft_control%qs_control%cdft_control%strength(:)
490 : CASE (outer_scf_basis_center_opt)
491 0 : outer_scf_history(1, ivec) = 0.0_dp
492 : CASE DEFAULT
493 4301 : CPABORT("Unknown outer SCF type")
494 : END SELECT
495 4301 : CALL set_qs_env(qs_env, outer_scf_ihistory=outer_scf_ihistory)
496 : ! multilinear extrapolation
497 4301 : nvec = MIN(nhistory, outer_scf_ihistory)
498 4301 : alpha = nvec
499 4301 : ivec = 1 + MODULO(outer_scf_ihistory - 1, nhistory)
500 8632 : extrapolation(:) = alpha*outer_scf_history(:, ivec)
501 9323 : DO ihis = 2, nvec
502 5022 : alpha = -1.0_dp*alpha*REAL(nvec - ihis + 1, dp)/REAL(ihis, dp)
503 5022 : ivec = 1 + MODULO(outer_scf_ihistory - ihis, nhistory)
504 14353 : extrapolation(:) = extrapolation + alpha*outer_scf_history(:, ivec)
505 : END DO
506 :
507 : ! update qs_env to use this extrapolation
508 4327 : SELECT CASE (scf_control%outer_scf%type)
509 : CASE (outer_scf_none)
510 : ! nothing
511 : CASE (outer_scf_ddapc_constraint)
512 26 : ddapc_restraint_control%strength = extrapolation(1)
513 : CASE (outer_scf_s2_constraint)
514 0 : dft_control%qs_control%s2_restraint_control%strength = extrapolation(1)
515 : CASE (outer_scf_cdft_constraint)
516 718 : dft_control%qs_control%cdft_control%strength(:) = extrapolation(:)
517 : CASE (outer_scf_basis_center_opt)
518 : ! nothing to do
519 : CASE DEFAULT
520 4301 : CPABORT("Unknown outer SCF type")
521 : END SELECT
522 :
523 4301 : DEALLOCATE (extrapolation)
524 :
525 4301 : CALL timestop(handle)
526 :
527 4301 : 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 1646 : 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 2328 : SELECT CASE (dir)
548 : CASE (cdft2ot)
549 : ! Constraint -> OT
550 : ! Switch data in scf_control: first save values that might have changed
551 682 : IF (ASSOCIATED(scf_control%outer_scf%cdft_opt_control)) THEN
552 344 : CPASSERT(ASSOCIATED(cdft_control%constraint_control%cdft_opt_control))
553 : CALL cdft_opt_type_copy(cdft_control%constraint_control%cdft_opt_control, &
554 344 : scf_control%outer_scf%cdft_opt_control)
555 : ! OT SCF does not need cdft_opt_control
556 344 : CALL cdft_opt_type_release(scf_control%outer_scf%cdft_opt_control)
557 : END IF
558 : ! Now switch
559 682 : scf_control%outer_scf%have_scf = cdft_control%ot_control%have_scf
560 682 : scf_control%outer_scf%max_scf = cdft_control%ot_control%max_scf
561 682 : scf_control%outer_scf%eps_scf = cdft_control%ot_control%eps_scf
562 682 : scf_control%outer_scf%step_size = cdft_control%ot_control%step_size
563 682 : scf_control%outer_scf%type = cdft_control%ot_control%type
564 682 : scf_control%outer_scf%optimizer = cdft_control%ot_control%optimizer
565 682 : scf_control%outer_scf%diis_buffer_length = cdft_control%ot_control%diis_buffer_length
566 682 : 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 682 : cdft_control%constraint%iter_count = scf_env%outer_scf%iter_count
569 10064 : cdft_control%constraint%energy = scf_env%outer_scf%energy
570 19784 : cdft_control%constraint%variables = scf_env%outer_scf%variables
571 19784 : cdft_control%constraint%gradient = scf_env%outer_scf%gradient
572 10064 : cdft_control%constraint%count = scf_env%outer_scf%count
573 682 : cdft_control%constraint%deallocate_jacobian = scf_env%outer_scf%deallocate_jacobian
574 682 : IF (ASSOCIATED(scf_env%outer_scf%inv_jacobian)) THEN
575 152 : nvariables = SIZE(scf_env%outer_scf%inv_jacobian, 1)
576 152 : IF (.NOT. ASSOCIATED(cdft_control%constraint%inv_jacobian)) THEN
577 192 : ALLOCATE (cdft_control%constraint%inv_jacobian(nvariables, nvariables))
578 : END IF
579 1504 : cdft_control%constraint%inv_jacobian = scf_env%outer_scf%inv_jacobian
580 : END IF
581 : ! Now switch
582 682 : IF (ASSOCIATED(scf_env%outer_scf%energy)) THEN
583 682 : DEALLOCATE (scf_env%outer_scf%energy)
584 : END IF
585 2046 : ALLOCATE (scf_env%outer_scf%energy(scf_control%outer_scf%max_scf + 1))
586 2388 : scf_env%outer_scf%energy = 0.0_dp
587 682 : IF (ASSOCIATED(scf_env%outer_scf%variables)) THEN
588 682 : DEALLOCATE (scf_env%outer_scf%variables)
589 : END IF
590 2046 : ALLOCATE (scf_env%outer_scf%variables(1, scf_control%outer_scf%max_scf + 1))
591 4094 : scf_env%outer_scf%variables = 0.0_dp
592 682 : IF (ASSOCIATED(scf_env%outer_scf%gradient)) THEN
593 682 : DEALLOCATE (scf_env%outer_scf%gradient)
594 : END IF
595 2046 : ALLOCATE (scf_env%outer_scf%gradient(1, scf_control%outer_scf%max_scf + 1))
596 4094 : scf_env%outer_scf%gradient = 0.0_dp
597 682 : IF (ASSOCIATED(scf_env%outer_scf%count)) THEN
598 682 : DEALLOCATE (scf_env%outer_scf%count)
599 : END IF
600 2046 : ALLOCATE (scf_env%outer_scf%count(scf_control%outer_scf%max_scf + 1))
601 2388 : scf_env%outer_scf%count = 0
602 : ! OT SCF does not need Jacobian
603 682 : scf_env%outer_scf%deallocate_jacobian = .TRUE.
604 682 : IF (ASSOCIATED(scf_env%outer_scf%inv_jacobian)) THEN
605 152 : DEALLOCATE (scf_env%outer_scf%inv_jacobian)
606 : END IF
607 : CASE (ot2cdft)
608 : ! OT -> constraint
609 964 : scf_control%outer_scf%have_scf = cdft_control%constraint_control%have_scf
610 964 : scf_control%outer_scf%max_scf = cdft_control%constraint_control%max_scf
611 964 : scf_control%outer_scf%eps_scf = cdft_control%constraint_control%eps_scf
612 964 : scf_control%outer_scf%step_size = cdft_control%constraint_control%step_size
613 964 : scf_control%outer_scf%type = cdft_control%constraint_control%type
614 964 : scf_control%outer_scf%optimizer = cdft_control%constraint_control%optimizer
615 964 : scf_control%outer_scf%diis_buffer_length = cdft_control%constraint_control%diis_buffer_length
616 964 : 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 964 : cdft_control%constraint_control%cdft_opt_control)
619 964 : nvariables = SIZE(cdft_control%constraint%variables, 1)
620 964 : IF (ASSOCIATED(scf_env%outer_scf%energy)) THEN
621 964 : DEALLOCATE (scf_env%outer_scf%energy)
622 : END IF
623 2892 : ALLOCATE (scf_env%outer_scf%energy(scf_control%outer_scf%max_scf + 1))
624 13768 : scf_env%outer_scf%energy = cdft_control%constraint%energy
625 964 : IF (ASSOCIATED(scf_env%outer_scf%variables)) THEN
626 964 : DEALLOCATE (scf_env%outer_scf%variables)
627 : END IF
628 3856 : ALLOCATE (scf_env%outer_scf%variables(nvariables, scf_control%outer_scf%max_scf + 1))
629 26776 : scf_env%outer_scf%variables = cdft_control%constraint%variables
630 964 : IF (ASSOCIATED(scf_env%outer_scf%gradient)) THEN
631 964 : DEALLOCATE (scf_env%outer_scf%gradient)
632 : END IF
633 3856 : ALLOCATE (scf_env%outer_scf%gradient(nvariables, scf_control%outer_scf%max_scf + 1))
634 26776 : scf_env%outer_scf%gradient = cdft_control%constraint%gradient
635 964 : IF (ASSOCIATED(scf_env%outer_scf%count)) THEN
636 964 : DEALLOCATE (scf_env%outer_scf%count)
637 : END IF
638 2892 : ALLOCATE (scf_env%outer_scf%count(scf_control%outer_scf%max_scf + 1))
639 13768 : scf_env%outer_scf%count = cdft_control%constraint%count
640 964 : scf_env%outer_scf%iter_count = cdft_control%constraint%iter_count
641 964 : scf_env%outer_scf%deallocate_jacobian = cdft_control%constraint%deallocate_jacobian
642 964 : IF (ASSOCIATED(cdft_control%constraint%inv_jacobian)) THEN
643 188 : IF (ASSOCIATED(scf_env%outer_scf%inv_jacobian)) THEN
644 0 : DEALLOCATE (scf_env%outer_scf%inv_jacobian)
645 : END IF
646 752 : ALLOCATE (scf_env%outer_scf%inv_jacobian(nvariables, nvariables))
647 1864 : scf_env%outer_scf%inv_jacobian = cdft_control%constraint%inv_jacobian
648 : END IF
649 : CASE DEFAULT
650 1646 : CPABORT("Switching direction should be either ot2cdft or cdft2ot")
651 : END SELECT
652 :
653 1646 : 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 0 : 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 0 : REAL(kind=dp), DIMENSION(:, :), POINTER :: gradient_history, outer_scf_history, &
669 0 : variable_history
670 :
671 0 : 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 0 : variable_history=variable_history)
677 0 : CPASSERT(SIZE(outer_scf_history, 2) > 0)
678 0 : outer_scf_ihistory = 0
679 0 : outer_scf_history = 0.0_dp
680 0 : gradient_history = 0.0_dp
681 0 : variable_history = 0.0_dp
682 0 : CALL set_qs_env(qs_env, outer_scf_ihistory=outer_scf_ihistory)
683 :
684 0 : CALL timestop(handle)
685 :
686 0 : END SUBROUTINE outer_loop_purge_history
687 :
688 174 : END MODULE qs_outer_scf
|
