Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 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 4856 : 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 4856 : 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 4856 : res = 0
91 : END SELECT
92 :
93 4856 : 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 5654 : 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 5654 : CALL timeset(routineN, handle)
118 :
119 : CALL get_qs_env(qs_env=qs_env, scf_control=scf_control, &
120 5654 : dft_control=dft_control, energy=energy)
121 5654 : CPASSERT(scf_control%outer_scf%have_scf)
122 :
123 5654 : ihistory = scf_env%outer_scf%iter_count
124 5654 : CPASSERT(ihistory <= SIZE(scf_env%outer_scf%energy, 1))
125 :
126 5654 : scf_env%outer_scf%energy(ihistory) = energy%total
127 :
128 10606 : 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 4952 : scf_env%outer_scf%variables(1, ihistory) = scf_env%iter_delta
132 4952 : 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 626 : CPASSERT(dft_control%qs_control%cdft)
157 626 : cdft_control => dft_control%qs_control%cdft_control
158 1386 : DO ivar = 1, SIZE(scf_env%outer_scf%gradient, 1)
159 760 : scf_env%outer_scf%variables(ivar, ihistory) = cdft_control%strength(ivar)
160 : scf_env%outer_scf%gradient(ivar, ihistory) = cdft_control%value(ivar) - &
161 1386 : 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 5654 : CPABORT("")
169 :
170 : END SELECT
171 :
172 5654 : CALL timestop(handle)
173 :
174 5654 : 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 1089 : 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 1089 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: ev
196 1089 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: a, b, f, x
197 1089 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: inv_jacobian
198 :
199 1089 : CALL timeset(routineN, handle)
200 :
201 1089 : ihistory = scf_env%outer_scf%iter_count
202 1089 : optimizer_type = scf_control%outer_scf%optimizer
203 1089 : NULLIFY (inv_jacobian)
204 :
205 1089 : 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) .GT. scf_control%outer_scf%bisect_trust_count) CYCLE
221 12 : 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 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 1706 : 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 176 : scf_control%outer_scf%step_size*scf_env%outer_scf%gradient(:, ihistory)
251 : CASE (outer_scf_optimizer_diis)
252 86 : CPASSERT(scf_control%outer_scf%diis_buffer_length > 0)
253 : ! set up DIIS matrix
254 86 : nb = MIN(ihistory, scf_control%outer_scf%diis_buffer_length)
255 86 : IF (nb < 2) THEN
256 : optimizer_type = outer_scf_optimizer_sd
257 : CYCLE
258 : ELSE
259 256 : ALLOCATE (b(nb + 1, nb + 1), a(nb + 1, nb + 1), ev(nb + 1))
260 98 : DO I = 1, nb
261 200 : DO J = I, nb
262 102 : ibuf = ihistory - nb + i
263 102 : jbuf = ihistory - nb + j
264 : b(I, J) = DOT_PRODUCT(scf_env%outer_scf%gradient(:, ibuf), &
265 204 : scf_env%outer_scf%gradient(:, jbuf))
266 168 : b(J, I) = b(I, J)
267 : END DO
268 : END DO
269 130 : b(nb + 1, :) = -1.0_dp
270 130 : b(:, nb + 1) = -1.0_dp
271 32 : b(nb + 1, nb + 1) = 0.0_dp
272 :
273 32 : CALL diamat_all(b, ev)
274 432 : a(:, :) = b
275 130 : DO I = 1, nb + 1
276 130 : IF (ABS(ev(I)) .LT. 1.0E-12_dp) THEN
277 10 : a(:, I) = 0.0_dp
278 : ELSE
279 390 : a(:, I) = a(:, I)/ev(I)
280 : END IF
281 : END DO
282 724 : ev(:) = -MATMUL(a, b(nb + 1, :))
283 :
284 64 : scf_env%outer_scf%variables(:, ihistory + 1) = 0.0_dp
285 98 : DO i = 1, nb
286 66 : ibuf = ihistory - nb + i
287 : scf_env%outer_scf%variables(:, ihistory + 1) = scf_env%outer_scf%variables(:, ihistory + 1) + &
288 164 : ev(i)*scf_env%outer_scf%variables(:, ibuf)
289 : END DO
290 32 : 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 8 : 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) .LT. 3) &
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 6 : nvar = SIZE(scf_env%outer_scf%gradient, 1)
326 24 : ALLOCATE (f(nvar, 1), x(nvar, 1))
327 12 : DO i = 1, nvar
328 6 : f(i, 1) = scf_env%outer_scf%gradient(i, ihistory) - scf_env%outer_scf%gradient(i, ihistory - 1)
329 12 : x(i, 1) = scf_env%outer_scf%variables(i, ihistory) - scf_env%outer_scf%variables(i, ihistory - 1)
330 : END DO
331 10 : 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 92 : scale = SUM(MATMUL(TRANSPOSE(x), MATMUL(inv_jacobian, f)))
337 4 : scale = 1.0_dp/scale
338 4 : IF (scale < 1.0E-12_dp) scale = 1.0E-12_dp
339 28 : inv_jacobian = inv_jacobian + scale*MATMUL((x - MATMUL(inv_jacobian, f)), &
340 120 : 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 22 : scale = SUM(MATMUL(TRANSPOSE(f), f))
345 2 : scale = 1.0_dp/scale
346 2 : IF (scale < 1.0E-12_dp) scale = 1.0E-12_dp
347 52 : 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 6 : cp_to_string(scf_control%outer_scf%cdft_opt_control%broyden_type))
352 : END SELECT
353 : ! Clean up
354 6 : 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 160 : MATMUL(inv_jacobian, scf_env%outer_scf%gradient(:, ihistory))
360 16 : scf_control%outer_scf%cdft_opt_control%broyden_update = .TRUE.
361 : CASE (outer_scf_optimizer_newton, outer_scf_optimizer_newton_ls)
362 94 : CPASSERT(ASSOCIATED(scf_env%outer_scf%inv_jacobian))
363 94 : 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 1612 : MATMUL(inv_jacobian, scf_env%outer_scf%gradient(:, ihistory))
367 : CASE DEFAULT
368 1129 : CPABORT("")
369 : END SELECT
370 : EXIT
371 : END DO
372 : END IF
373 :
374 1089 : CALL timestop(handle)
375 :
376 2178 : 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 1203 : 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 1203 : CALL timeset(routineN, handle)
401 :
402 1203 : CALL get_qs_env(qs_env=qs_env, scf_control=scf_control, dft_control=dft_control)
403 1203 : ihistory = scf_env%outer_scf%iter_count
404 :
405 1253 : SELECT CASE (scf_control%outer_scf%type)
406 : CASE (outer_scf_none)
407 : ! do nothing
408 : CASE (outer_scf_ddapc_constraint)
409 50 : DO n = 1, SIZE(dft_control%qs_control%ddapc_restraint_control)
410 50 : NULLIFY (ddapc_restraint_control)
411 50 : ddapc_restraint_control => dft_control%qs_control%ddapc_restraint_control(n)
412 50 : is_constraint = (ddapc_restraint_control%functional_form == do_ddapc_constraint)
413 50 : IF (is_constraint) EXIT
414 : END DO
415 50 : ddapc_restraint_control%strength = scf_env%outer_scf%variables(1, ihistory + 1)
416 : CASE (outer_scf_s2_constraint)
417 0 : s2_restraint_control => dft_control%qs_control%s2_restraint_control
418 0 : s2_restraint_control%strength = scf_env%outer_scf%variables(1, ihistory + 1)
419 : CASE (outer_scf_cdft_constraint)
420 300 : cdft_control => dft_control%qs_control%cdft_control
421 1408 : cdft_control%strength(:) = scf_env%outer_scf%variables(:, ihistory + 1)
422 : CASE (outer_scf_basis_center_opt)
423 0 : CALL qs_update_basis_center_pos(qs_env)
424 : CASE DEFAULT
425 1203 : CPABORT("")
426 : END SELECT
427 :
428 1203 : CALL timestop(handle)
429 :
430 1203 : 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 3831 : 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 3831 : REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: extrapolation
452 3831 : 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 3831 : 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 3831 : scf_control=scf_control, dft_control=dft_control)
462 :
463 3831 : nvariables = SIZE(outer_scf_history, 1)
464 3831 : nhistory = SIZE(outer_scf_history, 2)
465 11493 : ALLOCATE (extrapolation(nvariables))
466 3831 : CPASSERT(nhistory > 0)
467 :
468 : ! add the current version of qs_env to the history
469 3831 : outer_scf_ihistory = outer_scf_ihistory + 1
470 3831 : ivec = 1 + MODULO(outer_scf_ihistory - 1, nhistory)
471 7310 : SELECT CASE (scf_control%outer_scf%type)
472 : CASE (outer_scf_none)
473 3479 : outer_scf_history(1, ivec) = 0.0_dp
474 : CASE (outer_scf_ddapc_constraint)
475 26 : DO n = 1, SIZE(dft_control%qs_control%ddapc_restraint_control)
476 26 : NULLIFY (ddapc_restraint_control)
477 26 : ddapc_restraint_control => dft_control%qs_control%ddapc_restraint_control(n)
478 26 : is_constraint = (ddapc_restraint_control%functional_form == do_ddapc_constraint)
479 26 : IF (is_constraint) EXIT
480 : END DO
481 : outer_scf_history(1, ivec) = &
482 26 : ddapc_restraint_control%strength
483 : CASE (outer_scf_s2_constraint)
484 : outer_scf_history(1, ivec) = &
485 0 : dft_control%qs_control%s2_restraint_control%strength
486 : CASE (outer_scf_cdft_constraint)
487 : outer_scf_history(:, ivec) = &
488 1364 : dft_control%qs_control%cdft_control%strength(:)
489 : CASE (outer_scf_basis_center_opt)
490 0 : outer_scf_history(1, ivec) = 0.0_dp
491 : CASE DEFAULT
492 3831 : CPABORT("")
493 : END SELECT
494 3831 : CALL set_qs_env(qs_env, outer_scf_ihistory=outer_scf_ihistory)
495 : ! multilinear extrapolation
496 3831 : nvec = MIN(nhistory, outer_scf_ihistory)
497 3831 : alpha = nvec
498 3831 : ivec = 1 + MODULO(outer_scf_ihistory - 1, nhistory)
499 7692 : extrapolation(:) = alpha*outer_scf_history(:, ivec)
500 8461 : DO ihis = 2, nvec
501 4630 : alpha = -1.0_dp*alpha*REAL(nvec - ihis + 1, dp)/REAL(ihis, dp)
502 4630 : ivec = 1 + MODULO(outer_scf_ihistory - ihis, nhistory)
503 13099 : extrapolation(:) = extrapolation + alpha*outer_scf_history(:, ivec)
504 : END DO
505 :
506 : ! update qs_env to use this extrapolation
507 3857 : SELECT CASE (scf_control%outer_scf%type)
508 : CASE (outer_scf_none)
509 : ! nothing
510 : CASE (outer_scf_ddapc_constraint)
511 26 : ddapc_restraint_control%strength = extrapolation(1)
512 : CASE (outer_scf_s2_constraint)
513 0 : dft_control%qs_control%s2_restraint_control%strength = extrapolation(1)
514 : CASE (outer_scf_cdft_constraint)
515 682 : dft_control%qs_control%cdft_control%strength(:) = extrapolation(:)
516 : CASE (outer_scf_basis_center_opt)
517 : ! nothing to do
518 : CASE DEFAULT
519 3831 : CPABORT("")
520 : END SELECT
521 :
522 3831 : DEALLOCATE (extrapolation)
523 :
524 3831 : CALL timestop(handle)
525 :
526 3831 : 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 1516 : 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 2142 : SELECT CASE (dir)
547 : CASE (cdft2ot)
548 : ! Constraint -> OT
549 : ! Switch data in scf_control: first save values that might have changed
550 626 : IF (ASSOCIATED(scf_control%outer_scf%cdft_opt_control)) THEN
551 344 : CPASSERT(ASSOCIATED(cdft_control%constraint_control%cdft_opt_control))
552 : CALL cdft_opt_type_copy(cdft_control%constraint_control%cdft_opt_control, &
553 344 : scf_control%outer_scf%cdft_opt_control)
554 : ! OT SCF does not need cdft_opt_control
555 344 : CALL cdft_opt_type_release(scf_control%outer_scf%cdft_opt_control)
556 : END IF
557 : ! Now switch
558 626 : scf_control%outer_scf%have_scf = cdft_control%ot_control%have_scf
559 626 : scf_control%outer_scf%max_scf = cdft_control%ot_control%max_scf
560 626 : scf_control%outer_scf%eps_scf = cdft_control%ot_control%eps_scf
561 626 : scf_control%outer_scf%step_size = cdft_control%ot_control%step_size
562 626 : scf_control%outer_scf%type = cdft_control%ot_control%type
563 626 : scf_control%outer_scf%optimizer = cdft_control%ot_control%optimizer
564 626 : scf_control%outer_scf%diis_buffer_length = cdft_control%ot_control%diis_buffer_length
565 626 : 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 626 : cdft_control%constraint%iter_count = scf_env%outer_scf%iter_count
568 9008 : cdft_control%constraint%energy = scf_env%outer_scf%energy
569 17784 : cdft_control%constraint%variables = scf_env%outer_scf%variables
570 17784 : cdft_control%constraint%gradient = scf_env%outer_scf%gradient
571 9008 : cdft_control%constraint%count = scf_env%outer_scf%count
572 626 : cdft_control%constraint%deallocate_jacobian = scf_env%outer_scf%deallocate_jacobian
573 626 : IF (ASSOCIATED(scf_env%outer_scf%inv_jacobian)) THEN
574 152 : nvariables = SIZE(scf_env%outer_scf%inv_jacobian, 1)
575 152 : IF (.NOT. ASSOCIATED(cdft_control%constraint%inv_jacobian)) &
576 192 : ALLOCATE (cdft_control%constraint%inv_jacobian(nvariables, nvariables))
577 1504 : cdft_control%constraint%inv_jacobian = scf_env%outer_scf%inv_jacobian
578 : END IF
579 : ! Now switch
580 626 : IF (ASSOCIATED(scf_env%outer_scf%energy)) &
581 626 : DEALLOCATE (scf_env%outer_scf%energy)
582 1878 : ALLOCATE (scf_env%outer_scf%energy(scf_control%outer_scf%max_scf + 1))
583 2180 : scf_env%outer_scf%energy = 0.0_dp
584 626 : IF (ASSOCIATED(scf_env%outer_scf%variables)) &
585 626 : DEALLOCATE (scf_env%outer_scf%variables)
586 1878 : ALLOCATE (scf_env%outer_scf%variables(1, scf_control%outer_scf%max_scf + 1))
587 3734 : scf_env%outer_scf%variables = 0.0_dp
588 626 : IF (ASSOCIATED(scf_env%outer_scf%gradient)) &
589 626 : DEALLOCATE (scf_env%outer_scf%gradient)
590 1878 : ALLOCATE (scf_env%outer_scf%gradient(1, scf_control%outer_scf%max_scf + 1))
591 3734 : scf_env%outer_scf%gradient = 0.0_dp
592 626 : IF (ASSOCIATED(scf_env%outer_scf%count)) &
593 626 : DEALLOCATE (scf_env%outer_scf%count)
594 1878 : ALLOCATE (scf_env%outer_scf%count(scf_control%outer_scf%max_scf + 1))
595 2180 : scf_env%outer_scf%count = 0
596 : ! OT SCF does not need Jacobian
597 626 : scf_env%outer_scf%deallocate_jacobian = .TRUE.
598 626 : IF (ASSOCIATED(scf_env%outer_scf%inv_jacobian)) &
599 152 : DEALLOCATE (scf_env%outer_scf%inv_jacobian)
600 : CASE (ot2cdft)
601 : ! OT -> constraint
602 890 : scf_control%outer_scf%have_scf = cdft_control%constraint_control%have_scf
603 890 : scf_control%outer_scf%max_scf = cdft_control%constraint_control%max_scf
604 890 : scf_control%outer_scf%eps_scf = cdft_control%constraint_control%eps_scf
605 890 : scf_control%outer_scf%step_size = cdft_control%constraint_control%step_size
606 890 : scf_control%outer_scf%type = cdft_control%constraint_control%type
607 890 : scf_control%outer_scf%optimizer = cdft_control%constraint_control%optimizer
608 890 : scf_control%outer_scf%diis_buffer_length = cdft_control%constraint_control%diis_buffer_length
609 890 : 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 890 : cdft_control%constraint_control%cdft_opt_control)
612 890 : nvariables = SIZE(cdft_control%constraint%variables, 1)
613 890 : IF (ASSOCIATED(scf_env%outer_scf%energy)) &
614 890 : DEALLOCATE (scf_env%outer_scf%energy)
615 2670 : ALLOCATE (scf_env%outer_scf%energy(scf_control%outer_scf%max_scf + 1))
616 12424 : scf_env%outer_scf%energy = cdft_control%constraint%energy
617 890 : IF (ASSOCIATED(scf_env%outer_scf%variables)) &
618 890 : DEALLOCATE (scf_env%outer_scf%variables)
619 3560 : ALLOCATE (scf_env%outer_scf%variables(nvariables, scf_control%outer_scf%max_scf + 1))
620 24236 : scf_env%outer_scf%variables = cdft_control%constraint%variables
621 890 : IF (ASSOCIATED(scf_env%outer_scf%gradient)) &
622 890 : DEALLOCATE (scf_env%outer_scf%gradient)
623 3560 : ALLOCATE (scf_env%outer_scf%gradient(nvariables, scf_control%outer_scf%max_scf + 1))
624 24236 : scf_env%outer_scf%gradient = cdft_control%constraint%gradient
625 890 : IF (ASSOCIATED(scf_env%outer_scf%count)) &
626 890 : DEALLOCATE (scf_env%outer_scf%count)
627 2670 : ALLOCATE (scf_env%outer_scf%count(scf_control%outer_scf%max_scf + 1))
628 12424 : scf_env%outer_scf%count = cdft_control%constraint%count
629 890 : scf_env%outer_scf%iter_count = cdft_control%constraint%iter_count
630 890 : scf_env%outer_scf%deallocate_jacobian = cdft_control%constraint%deallocate_jacobian
631 890 : IF (ASSOCIATED(cdft_control%constraint%inv_jacobian)) THEN
632 188 : IF (ASSOCIATED(scf_env%outer_scf%inv_jacobian)) &
633 0 : DEALLOCATE (scf_env%outer_scf%inv_jacobian)
634 752 : ALLOCATE (scf_env%outer_scf%inv_jacobian(nvariables, nvariables))
635 1864 : scf_env%outer_scf%inv_jacobian = cdft_control%constraint%inv_jacobian
636 : END IF
637 : CASE DEFAULT
638 1516 : CPABORT("")
639 : END SELECT
640 :
641 1516 : 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 0 : 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 0 : REAL(kind=dp), DIMENSION(:, :), POINTER :: gradient_history, outer_scf_history, &
657 0 : variable_history
658 :
659 0 : 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 0 : variable_history=variable_history)
665 0 : CPASSERT(SIZE(outer_scf_history, 2) > 0)
666 0 : outer_scf_ihistory = 0
667 0 : outer_scf_history = 0.0_dp
668 0 : gradient_history = 0.0_dp
669 0 : variable_history = 0.0_dp
670 0 : CALL set_qs_env(qs_env, outer_scf_ihistory=outer_scf_ihistory)
671 :
672 0 : CALL timestop(handle)
673 :
674 0 : END SUBROUTINE outer_loop_purge_history
675 :
676 154 : END MODULE qs_outer_scf
|
