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, ONLY: getrf, &
17 : 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 2334 : 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 2334 : TYPE(cp2k_tblite_broyden_mixer_type), ALLOCATABLE :: broyden
73 :
74 2334 : ALLOCATE (broyden)
75 2334 : broyden%ndim = ndim
76 2334 : broyden%memory = memory
77 2334 : broyden%damp = damping
78 2334 : broyden%omega0 = omega0
79 2334 : broyden%min_weight = min_weight
80 2334 : broyden%max_weight = max_weight
81 2334 : broyden%weight_factor = weight_factor
82 9336 : ALLOCATE (broyden%a(memory, memory))
83 9336 : ALLOCATE (broyden%df(ndim, memory))
84 7002 : ALLOCATE (broyden%u(ndim, memory))
85 7002 : ALLOCATE (broyden%dq(ndim))
86 4668 : ALLOCATE (broyden%dqlast(ndim))
87 7002 : ALLOCATE (broyden%omega(memory))
88 4668 : ALLOCATE (broyden%q_in(ndim))
89 4668 : ALLOCATE (broyden%qlast_in(ndim))
90 2334 : CALL MOVE_ALLOC(broyden, mixer)
91 :
92 2334 : 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 53622 : 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 1202592 : self%q_in(self%iset + 1:self%iset + SIZE(qvec)) = qvec
104 53622 : self%iset = self%iset + SIZE(qvec)
105 :
106 53622 : 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 53740 : 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 1204514 : qvec - self%q_in(self%idif + 1:self%idif + SIZE(qvec))
119 53740 : self%idif = self%idif + SIZE(qvec)
120 :
121 53740 : 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 21560 : 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 21560 : self%iset = 0
135 21560 : self%idif = 0
136 21560 : self%iget = 0
137 21560 : 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 21560 : self%omega, self%df, self%u, self%a, info)
142 21560 : IF (info /= 0) CALL fatal_error(error, "Broyden mixing failed to obtain next iteration")
143 :
144 21560 : 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 53412 : 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 1198980 : qvec(:) = self%q_in(self%iget + 1:self%iget + SIZE(qvec))
156 53412 : self%iget = self%iget + SIZE(qvec)
157 :
158 53412 : END SUBROUTINE cp2k_tblite_mixer_get_1d
159 :
160 : ! **************************************************************************************************
161 : !> \brief Get RMS SCC residual.
162 : !> \param self mixer
163 : !> \return residual
164 : ! **************************************************************************************************
165 21872 : 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 1177706 : error = SQRT(SUM(self%dq**2)/SIZE(self%dq))
170 :
171 21872 : 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 21560 : SUBROUTINE cp2k_tblite_broyden(n, q, qlast, dq, dqlast, iter, memory, alpha, omega0, minw, &
194 21560 : 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 :: i, it1, itn, j
207 : REAL(wp) :: inv
208 21560 : REAL(wp), ALLOCATABLE, DIMENSION(:, :) :: beta, c
209 :
210 21560 : info = 0
211 21560 : itn = iter - 1
212 21560 : it1 = MOD(itn - 1, memory) + 1
213 :
214 21560 : IF (iter == 1) THEN
215 74444 : dqlast(:) = dq
216 74444 : qlast(:) = q
217 74444 : q(:) = q + alpha*dq
218 4300 : RETURN
219 : END IF
220 :
221 116460 : ALLOCATE (beta(MIN(memory, itn), MIN(memory, itn)), c(MIN(memory, itn), 1))
222 :
223 1092684 : omega(it1) = SQRT(DOT_PRODUCT(dq, dq))
224 19410 : IF (omega(it1) > (wfac/maxw)) THEN
225 10859 : omega(it1) = wfac/omega(it1)
226 : ELSE
227 8551 : omega(it1) = maxw
228 : END IF
229 19410 : IF (omega(it1) < minw) omega(it1) = minw
230 :
231 1092684 : df(:, it1) = dq - dqlast
232 1092684 : inv = MAX(SQRT(DOT_PRODUCT(df(:, it1), df(:, it1))), EPSILON(1.0_wp))
233 19410 : inv = 1.0_wp/inv
234 1092684 : df(:, it1) = inv*df(:, it1)
235 :
236 204158 : DO j = MAX(1, itn - memory + 1), itn
237 184748 : i = MOD(j - 1, memory) + 1
238 13521520 : a(i, it1) = DOT_PRODUCT(df(:, i), df(:, it1))
239 184748 : a(it1, i) = a(i, it1)
240 13540930 : c(i, 1) = omega(i)*DOT_PRODUCT(df(:, i), dq)
241 : END DO
242 :
243 204158 : DO j = MAX(1, itn - memory + 1), itn
244 184748 : i = MOD(j - 1, memory) + 1
245 3540876 : beta(:it1, i) = omega(:it1)*omega(i)*a(:it1, i)
246 204158 : beta(i, i) = beta(i, i) + omega0*omega0
247 : END DO
248 :
249 19410 : CALL cp2k_tblite_lineq(beta, c, info)
250 19410 : IF (info /= 0) RETURN
251 :
252 1092684 : u(:, it1) = alpha*df(:, it1) + inv*(q - qlast)
253 1092684 : dqlast(:) = dq
254 1092684 : qlast(:) = q
255 1092684 : q(:) = q + alpha*dq
256 :
257 204158 : DO j = MAX(1, itn - memory + 1), itn
258 184748 : i = MOD(j - 1, memory) + 1
259 13540930 : q(:) = q - omega(i)*c(i, 1)*u(:, i)
260 : END DO
261 :
262 43120 : END SUBROUTINE cp2k_tblite_broyden
263 :
264 : ! **************************************************************************************************
265 : !> \brief Solve the small Broyden linear system.
266 : !> \param a ...
267 : !> \param c ...
268 : !> \param info ...
269 : ! **************************************************************************************************
270 19410 : SUBROUTINE cp2k_tblite_lineq(a, c, info)
271 : REAL(wp), DIMENSION(:, :), INTENT(INOUT) :: a, c
272 : INTEGER, INTENT(OUT) :: info
273 :
274 19410 : INTEGER, ALLOCATABLE, DIMENSION(:) :: ipiv
275 :
276 58230 : ALLOCATE (ipiv(SIZE(a, 1)))
277 19410 : CALL getrf(a, ipiv, info)
278 19410 : IF (info == 0) CALL getrs(a, c, ipiv, info, trans="t")
279 :
280 19410 : END SUBROUTINE cp2k_tblite_lineq
281 : #endif
282 :
283 7002 : END MODULE tblite_scc_mixer
|
