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 2 : PROGRAM qs_charge_mixing_unittest
9 2 : USE ieee_arithmetic, ONLY: ieee_is_finite
10 : USE kinds, ONLY: dp
11 : USE message_passing, ONLY: mp_para_env_type, &
12 : mp_world_finalize, &
13 : mp_world_init
14 : USE qs_charge_mixing, ONLY: charge_mixing
15 : USE qs_density_mixing_types, ONLY: broyden_mixing_nr, &
16 : mixing_storage_release, &
17 : mixing_storage_type
18 : #if defined(__TBLITE)
19 : USE mctc_env, ONLY: error_type
20 : USE tblite_scf_mixer_type, ONLY: mixer_type
21 : USE tblite_scc_mixer, ONLY: new_cp2k_tblite_mixer
22 : #endif
23 :
24 : IMPLICIT NONE
25 :
26 : REAL(KIND=dp), DIMENSION(1, 1) :: charges
27 : TYPE(mp_para_env_type), POINTER :: para_env
28 : TYPE(mixing_storage_type), POINTER :: mixing_store
29 :
30 2 : NULLIFY (mixing_store, para_env)
31 2 : ALLOCATE (para_env)
32 2 : CALL mp_world_init(para_env)
33 :
34 8 : ALLOCATE (mixing_store)
35 2 : mixing_store%ncall = 0
36 2 : mixing_store%nbuffer = 4
37 2 : mixing_store%n_simple_mix = 0
38 2 : mixing_store%nskip_mixing = 0
39 2 : mixing_store%alpha = 0.1_dp
40 2 : mixing_store%broy_w0 = 0.01_dp
41 2 : mixing_store%nat_local = 1
42 2 : mixing_store%max_shell = 1
43 : ALLOCATE (mixing_store%acharge(1, 1, 4), mixing_store%dacharge(1, 1, 4), &
44 : mixing_store%dfbroy(1, 1, 4), mixing_store%ubroy(1, 1, 4), &
45 2 : mixing_store%abroy(4, 4), mixing_store%wbroy(4), mixing_store%atlist(1))
46 26 : mixing_store%acharge = 0.0_dp
47 26 : mixing_store%dacharge = 0.0_dp
48 26 : mixing_store%dfbroy = 0.0_dp
49 26 : mixing_store%ubroy = 0.0_dp
50 42 : mixing_store%abroy = 0.0_dp
51 10 : mixing_store%wbroy = 0.0_dp
52 4 : mixing_store%atlist = 1
53 :
54 6 : charges = 1.0_dp
55 2 : CALL charge_mixing(broyden_mixing_nr, mixing_store, charges, para_env, 1)
56 6 : charges = 0.5_dp
57 2 : CALL charge_mixing(broyden_mixing_nr, mixing_store, charges, para_env, 2)
58 :
59 : ! Reproduce a degenerate Broyden history by supplying the previous residual again.
60 2 : charges(1, 1) = mixing_store%acharge(1, 1, 2) - mixing_store%dacharge(1, 1, 1)
61 2 : CALL charge_mixing(broyden_mixing_nr, mixing_store, charges, para_env, 3)
62 :
63 6 : IF (.NOT. ALL(ieee_is_finite(charges))) THEN
64 0 : ERROR STOP "Broyden mixing produced non-finite charges for identical residuals"
65 : END IF
66 6 : IF (.NOT. ALL(mixing_store%dfbroy(:, :, 1) == 0.0_dp)) THEN
67 0 : ERROR STOP "Degenerate Broyden direction was not discarded"
68 : END IF
69 :
70 2 : CALL mixing_storage_release(mixing_store)
71 2 : DEALLOCATE (mixing_store)
72 : #if defined(__TBLITE)
73 2 : CALL test_tblite_history()
74 : #endif
75 2 : DEALLOCATE (para_env)
76 2 : CALL mp_world_finalize()
77 :
78 : #if defined(__TBLITE)
79 : CONTAINS
80 :
81 : ! **************************************************************************************************
82 : !> \brief A wrapped Broyden history must match a fresh replay of the retained input/output pairs.
83 : ! **************************************************************************************************
84 2 : SUBROUTINE test_tblite_history()
85 : INTEGER, PARAMETER :: nstep = 16, nvar = 4
86 :
87 : INTEGER :: config, i, j, memory, step
88 : REAL(KIND=dp) :: damping, max_weight, min_weight, omega0, &
89 : weight_factor
90 : REAL(KIND=dp), DIMENSION(nvar) :: actual, expected
91 : REAL(KIND=dp), DIMENSION(nvar, nstep) :: qin, qout
92 :
93 6 : CLASS(mixer_type), ALLOCATABLE :: mixer, replay
94 2 : TYPE(error_type), ALLOCATABLE :: error
95 :
96 34 : DO step = 1, nstep
97 162 : DO i = 1, nvar
98 128 : qin(i, step) = SIN(0.37_dp*REAL(i*step, dp))
99 160 : qout(i, step) = qin(i, step) + 0.05_dp*COS(0.21_dp*REAL(i + 2*step, dp))
100 : END DO
101 : END DO
102 6 : DO config = 1, 2
103 4 : damping = 0.4_dp
104 4 : omega0 = 0.01_dp
105 4 : min_weight = 1.0_dp
106 4 : max_weight = 100000.0_dp
107 4 : weight_factor = 0.01_dp
108 4 : IF (config == 2) THEN
109 2 : damping = 0.25_dp
110 2 : omega0 = 0.1_dp
111 2 : min_weight = 0.5_dp
112 2 : max_weight = 3.0_dp
113 2 : weight_factor = 0.2_dp
114 : END IF
115 22 : DO memory = 1, 4
116 30 : CALL new_cp2k_tblite_mixer(mixer, memory, nvar, damping, omega0, min_weight, max_weight, weight_factor)
117 276 : DO step = 1, nstep
118 256 : CALL mixer%set(qin(:, step))
119 256 : CALL mixer%diff(qout(:, step))
120 256 : CALL mixer%next(error)
121 256 : IF (ALLOCATED(error)) ERROR STOP "Wrapped tblite mixer failed"
122 256 : CALL mixer%get(actual)
123 :
124 : ! Replay at most memory + 1 pairs: this builds exactly the retained
125 : ! secant directions without ever overwriting a history slot.
126 510 : CALL new_cp2k_tblite_mixer(replay, memory, nvar, damping, omega0, min_weight, max_weight, weight_factor)
127 1072 : DO j = MAX(1, step - memory), step
128 816 : CALL replay%set(qin(:, j))
129 816 : CALL replay%diff(qout(:, j))
130 816 : CALL replay%next(error)
131 1072 : IF (ALLOCATED(error)) ERROR STOP "Replayed tblite mixer failed"
132 : END DO
133 256 : CALL replay%get(expected)
134 1280 : IF (.NOT. ALL(ieee_is_finite(actual))) ERROR STOP "Non-finite wrapped tblite history"
135 1296 : IF (MAXVAL(ABS(actual - expected)) > 1.0E-10_dp) THEN
136 0 : ERROR STOP "Wrapped tblite history differs from an unwrapped replay"
137 : END IF
138 : END DO
139 : END DO
140 : END DO
141 :
142 4 : END SUBROUTINE test_tblite_history
143 : #endif
144 :
145 : END PROGRAM qs_charge_mixing_unittest
146 : ! vim: set ts=3 sw=3 tw=132 :
|
