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 CP2K-side tblite-compatible SCC Broyden mixer.
10 : ! **************************************************************************************************
11 : MODULE tblite_scc_mixer
12 : #if defined(__TBLITE)
13 : USE mctc_env, ONLY: error_type, &
14 : fatal_error, &
15 : wp
16 : USE tblite_lapack_getrf, ONLY: getrf => wrap_getrf
17 : USE tblite_lapack_getrs, ONLY: getrs => wrap_getrs
18 : USE tblite_scf_mixer_type, ONLY: mixer_type
19 : #endif
20 :
21 : IMPLICIT NONE
22 :
23 : PRIVATE
24 :
25 : #if defined(__TBLITE)
26 : PUBLIC :: new_cp2k_tblite_mixer
27 :
28 : TYPE, EXTENDS(mixer_type) :: cp2k_tblite_broyden_mixer_type
29 : INTEGER :: idif = 0
30 : INTEGER :: iget = 0
31 : INTEGER :: iset = 0
32 : INTEGER :: iter = 0
33 : INTEGER :: memory = 0
34 : INTEGER :: ndim = 0
35 : REAL(wp) :: damp = 0.0_wp
36 : REAL(wp) :: max_weight = 0.0_wp
37 : REAL(wp) :: min_weight = 0.0_wp
38 : REAL(wp) :: omega0 = 0.0_wp
39 : REAL(wp) :: weight_factor = 0.0_wp
40 : REAL(wp), ALLOCATABLE, DIMENSION(:, :) :: a, df, u
41 : REAL(wp), ALLOCATABLE, DIMENSION(:) :: dq, dqlast, omega, q_in, qlast_in
42 : CONTAINS
43 : PROCEDURE :: diff_1d => cp2k_tblite_mixer_diff_1d
44 : PROCEDURE :: get_1d => cp2k_tblite_mixer_get_1d
45 : PROCEDURE :: get_error => cp2k_tblite_mixer_get_error
46 : PROCEDURE :: next => cp2k_tblite_mixer_next
47 : PROCEDURE :: set_1d => cp2k_tblite_mixer_set_1d
48 : END TYPE cp2k_tblite_broyden_mixer_type
49 : #endif
50 :
51 : CONTAINS
52 :
53 : #if defined(__TBLITE)
54 : ! **************************************************************************************************
55 : !> \brief Create a CP2K-side tblite-compatible Broyden mixer.
56 : !> \param mixer mixer instance
57 : !> \param memory Broyden history length
58 : !> \param ndim number of SCC variables
59 : !> \param damping first-step damping
60 : !> \param omega0 Broyden regularization weight
61 : !> \param min_weight minimum dynamic Broyden weight
62 : !> \param max_weight maximum dynamic Broyden weight
63 : !> \param weight_factor residual-to-weight scaling factor
64 : ! **************************************************************************************************
65 2930 : SUBROUTINE new_cp2k_tblite_mixer(mixer, memory, ndim, damping, omega0, min_weight, max_weight, &
66 : weight_factor)
67 : CLASS(mixer_type), ALLOCATABLE, INTENT(OUT) :: mixer
68 : INTEGER, INTENT(IN) :: memory, ndim
69 : REAL(wp), INTENT(IN) :: damping, omega0, min_weight, max_weight, &
70 : weight_factor
71 :
72 2930 : TYPE(cp2k_tblite_broyden_mixer_type), ALLOCATABLE :: broyden
73 :
74 2930 : ALLOCATE (broyden)
75 2930 : broyden%ndim = ndim
76 2930 : broyden%memory = memory
77 2930 : broyden%damp = damping
78 2930 : broyden%omega0 = omega0
79 2930 : broyden%min_weight = min_weight
80 2930 : broyden%max_weight = max_weight
81 2930 : broyden%weight_factor = weight_factor
82 11720 : ALLOCATE (broyden%a(memory, memory))
83 11720 : ALLOCATE (broyden%df(ndim, memory))
84 8790 : ALLOCATE (broyden%u(ndim, memory))
85 8790 : ALLOCATE (broyden%dq(ndim))
86 5860 : ALLOCATE (broyden%dqlast(ndim))
87 8790 : ALLOCATE (broyden%omega(memory))
88 5860 : ALLOCATE (broyden%q_in(ndim))
89 5860 : ALLOCATE (broyden%qlast_in(ndim))
90 2930 : CALL MOVE_ALLOC(broyden, mixer)
91 :
92 2930 : END SUBROUTINE new_cp2k_tblite_mixer
93 :
94 : ! **************************************************************************************************
95 : !> \brief Set input SCC variables.
96 : !> \param self mixer
97 : !> \param qvec SCC-variable vector
98 : ! **************************************************************************************************
99 54442 : SUBROUTINE cp2k_tblite_mixer_set_1d(self, qvec)
100 : CLASS(cp2k_tblite_broyden_mixer_type), INTENT(INOUT) :: self
101 : REAL(wp), DIMENSION(:), INTENT(IN) :: qvec
102 :
103 1085464 : self%q_in(self%iset + 1:self%iset + SIZE(qvec)) = qvec
104 54442 : self%iset = self%iset + SIZE(qvec)
105 :
106 54442 : END SUBROUTINE cp2k_tblite_mixer_set_1d
107 :
108 : ! **************************************************************************************************
109 : !> \brief Set SCC residual.
110 : !> \param self mixer
111 : !> \param qvec output SCC-variable vector before mixing
112 : ! **************************************************************************************************
113 54454 : SUBROUTINE cp2k_tblite_mixer_diff_1d(self, qvec)
114 : CLASS(cp2k_tblite_broyden_mixer_type), INTENT(INOUT) :: self
115 : REAL(wp), DIMENSION(:), INTENT(IN) :: qvec
116 :
117 : self%dq(self%idif + 1:self%idif + SIZE(qvec)) = &
118 1085668 : qvec - self%q_in(self%idif + 1:self%idif + SIZE(qvec))
119 54454 : self%idif = self%idif + SIZE(qvec)
120 :
121 54454 : END SUBROUTINE cp2k_tblite_mixer_diff_1d
122 :
123 : ! **************************************************************************************************
124 : !> \brief Apply the next Broyden update.
125 : !> \param self mixer
126 : !> \param error error handling
127 : ! **************************************************************************************************
128 23472 : SUBROUTINE cp2k_tblite_mixer_next(self, error)
129 : CLASS(cp2k_tblite_broyden_mixer_type), INTENT(INOUT) :: self
130 : TYPE(error_type), ALLOCATABLE, INTENT(OUT) :: error
131 :
132 : INTEGER :: info
133 :
134 23472 : self%iset = 0
135 23472 : self%idif = 0
136 23472 : self%iget = 0
137 23472 : self%iter = self%iter + 1
138 : CALL cp2k_tblite_broyden(self%ndim, self%q_in, self%qlast_in, self%dq, self%dqlast, &
139 : self%iter, self%memory, self%damp, self%omega0, &
140 : self%min_weight, self%max_weight, self%weight_factor, &
141 23472 : self%omega, self%df, self%u, self%a, info)
142 23472 : IF (info /= 0) CALL fatal_error(error, "Broyden mixing failed to obtain next iteration")
143 :
144 23472 : END SUBROUTINE cp2k_tblite_mixer_next
145 :
146 : ! **************************************************************************************************
147 : !> \brief Get mixed SCC variables.
148 : !> \param self mixer
149 : !> \param qvec mixed SCC-variable vector
150 : ! **************************************************************************************************
151 53480 : SUBROUTINE cp2k_tblite_mixer_get_1d(self, qvec)
152 : CLASS(cp2k_tblite_broyden_mixer_type), INTENT(INOUT) :: self
153 : REAL(wp), DIMENSION(:), INTENT(OUT) :: qvec
154 :
155 1074682 : qvec(:) = self%q_in(self%iget + 1:self%iget + SIZE(qvec))
156 53480 : self%iget = self%iget + SIZE(qvec)
157 :
158 53480 : END SUBROUTINE cp2k_tblite_mixer_get_1d
159 :
160 : ! **************************************************************************************************
161 : !> \brief Get RMS SCC residual.
162 : !> \param self mixer
163 : !> \return residual
164 : ! **************************************************************************************************
165 22622 : PURE FUNCTION cp2k_tblite_mixer_get_error(self) RESULT(error)
166 : CLASS(cp2k_tblite_broyden_mixer_type), INTENT(IN) :: self
167 : REAL(wp) :: error
168 :
169 1049548 : error = SQRT(SUM(self%dq**2)/SIZE(self%dq))
170 :
171 22622 : END FUNCTION cp2k_tblite_mixer_get_error
172 :
173 : ! **************************************************************************************************
174 : !> \brief Modified Broyden update matching tblite's default algorithm, with explicit constants.
175 : !> \param n ...
176 : !> \param q ...
177 : !> \param qlast ...
178 : !> \param dq ...
179 : !> \param dqlast ...
180 : !> \param iter ...
181 : !> \param memory ...
182 : !> \param alpha ...
183 : !> \param omega0 ...
184 : !> \param minw ...
185 : !> \param maxw ...
186 : !> \param wfac ...
187 : !> \param omega ...
188 : !> \param df ...
189 : !> \param u ...
190 : !> \param a ...
191 : !> \param info ...
192 : ! **************************************************************************************************
193 23472 : SUBROUTINE cp2k_tblite_broyden(n, q, qlast, dq, dqlast, iter, memory, alpha, omega0, minw, &
194 23472 : maxw, wfac, omega, df, u, a, info)
195 : INTEGER, INTENT(IN) :: n
196 : REAL(wp), DIMENSION(n), INTENT(INOUT) :: q, qlast
197 : REAL(wp), DIMENSION(n), INTENT(IN) :: dq
198 : REAL(wp), DIMENSION(n), INTENT(INOUT) :: dqlast
199 : INTEGER, INTENT(IN) :: iter, memory
200 : REAL(wp), INTENT(IN) :: alpha, omega0, minw, maxw, wfac
201 : REAL(wp), DIMENSION(memory), INTENT(INOUT) :: omega
202 : REAL(wp), DIMENSION(n, memory), INTENT(INOUT) :: df, u
203 : REAL(wp), DIMENSION(memory, memory), INTENT(INOUT) :: a
204 : INTEGER, INTENT(OUT) :: info
205 :
206 : INTEGER :: first, i, ih, it1, itn, j, jh, nhist
207 23472 : INTEGER, ALLOCATABLE, DIMENSION(:) :: hist
208 : REAL(wp) :: inv
209 23472 : REAL(wp), ALLOCATABLE, DIMENSION(:, :) :: beta, c
210 :
211 23472 : info = 0
212 23472 : itn = iter - 1
213 :
214 23472 : IF (iter == 1) THEN
215 80970 : dqlast(:) = dq
216 80970 : qlast(:) = q
217 80970 : q(:) = q + alpha*dq
218 : RETURN
219 : END IF
220 :
221 20790 : it1 = MOD(itn - 1, memory) + 1
222 20790 : nhist = MIN(memory, itn)
223 20790 : first = MAX(1, itn - memory + 1)
224 166320 : ALLOCATE (hist(nhist), beta(nhist, nhist), c(nhist, 1))
225 : ! Keep the linear system in chronological order when the ring buffer wraps.
226 183960 : DO ih = 1, nhist
227 183960 : hist(ih) = MOD(first + ih - 2, memory) + 1
228 : END DO
229 :
230 965944 : omega(it1) = NORM2(dq)
231 20790 : IF (omega(it1) > (wfac/maxw)) THEN
232 11836 : omega(it1) = wfac/omega(it1)
233 : ELSE
234 8954 : omega(it1) = maxw
235 : END IF
236 20790 : IF (omega(it1) < minw) omega(it1) = minw
237 :
238 965944 : df(:, it1) = dq - dqlast
239 965944 : inv = MAX(NORM2(df(:, it1)), EPSILON(1.0_wp))
240 20790 : inv = 1.0_wp/inv
241 965944 : df(:, it1) = inv*df(:, it1)
242 :
243 183960 : DO ih = 1, nhist
244 163170 : i = hist(ih)
245 8669862 : a(i, it1) = DOT_PRODUCT(df(:, i), df(:, it1))
246 163170 : a(it1, i) = a(i, it1)
247 8690652 : c(ih, 1) = omega(i)*DOT_PRODUCT(df(:, i), dq)
248 : END DO
249 :
250 183960 : DO ih = 1, nhist
251 163170 : i = hist(ih)
252 3388352 : DO jh = 1, nhist
253 3225182 : j = hist(jh)
254 3388352 : beta(ih, jh) = omega(i)*omega(j)*a(i, j)
255 : END DO
256 183960 : beta(ih, ih) = beta(ih, ih) + omega0*omega0
257 : END DO
258 :
259 20790 : CALL cp2k_tblite_lineq(beta, c, info)
260 20790 : IF (info /= 0) RETURN
261 :
262 965944 : u(:, it1) = alpha*df(:, it1) + inv*(q - qlast)
263 965944 : dqlast(:) = dq
264 965944 : qlast(:) = q
265 965944 : q(:) = q + alpha*dq
266 :
267 183960 : DO ih = 1, nhist
268 163170 : i = hist(ih)
269 8690652 : q(:) = q - omega(i)*c(ih, 1)*u(:, i)
270 : END DO
271 :
272 44262 : END SUBROUTINE cp2k_tblite_broyden
273 :
274 : ! **************************************************************************************************
275 : !> \brief Solve the small Broyden linear system.
276 : !> \param a ...
277 : !> \param c ...
278 : !> \param info ...
279 : ! **************************************************************************************************
280 20790 : SUBROUTINE cp2k_tblite_lineq(a, c, info)
281 : REAL(wp), DIMENSION(:, :), INTENT(INOUT) :: a, c
282 : INTEGER, INTENT(OUT) :: info
283 :
284 20790 : INTEGER, ALLOCATABLE, DIMENSION(:) :: ipiv
285 :
286 62370 : ALLOCATE (ipiv(SIZE(a, 1)))
287 20790 : CALL getrf(a, ipiv, info)
288 20790 : IF (info == 0) CALL getrs(a, c, ipiv, info, trans="t")
289 :
290 20790 : END SUBROUTINE cp2k_tblite_lineq
291 : #endif
292 :
293 8790 : END MODULE tblite_scc_mixer
|
