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