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 History of minima, calculates, stores and compares fingerprints of minima.
10 : !> Used by Minima Hopping and Minima Crawling.
11 : !> \author Ole Schuett
12 : ! **************************************************************************************************
13 : MODULE glbopt_history
14 : USE input_section_types, ONLY: section_vals_type,&
15 : section_vals_val_get
16 : USE kinds, ONLY: dp
17 : #include "../base/base_uses.f90"
18 :
19 : IMPLICIT NONE
20 : PRIVATE
21 :
22 : TYPE history_fingerprint_type
23 : PRIVATE
24 : REAL(KIND=dp) :: Epot = 0.0
25 : REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: goedecker
26 : END TYPE history_fingerprint_type
27 :
28 : TYPE history_entry_type
29 : TYPE(history_fingerprint_type), POINTER :: p => Null()
30 : INTEGER :: id = -1
31 : END TYPE history_entry_type
32 :
33 : TYPE history_type
34 : PRIVATE
35 : TYPE(history_entry_type), DIMENSION(:), POINTER :: entries => Null()
36 : INTEGER :: length = 0
37 : INTEGER :: iw = -1
38 : REAL(KIND=dp) :: E_precision = 0.0
39 : REAL(KIND=dp) :: FP_precision = 0.0
40 : END TYPE history_type
41 :
42 : PUBLIC :: history_type, history_fingerprint_type
43 : PUBLIC :: history_init, history_finalize
44 : PUBLIC :: history_add, history_lookup
45 : PUBLIC :: history_fingerprint
46 : PUBLIC :: history_fingerprint_match
47 :
48 : LOGICAL, PARAMETER :: debug = .FALSE.
49 : INTEGER, PARAMETER :: history_grow_unit = 1000
50 : CONTAINS
51 :
52 : ! **************************************************************************************************
53 : !> \brief Initializes a history.
54 : !> \param history ...
55 : !> \param history_section ...
56 : !> \param iw ...
57 : !> \author Ole Schuett
58 : ! **************************************************************************************************
59 3 : SUBROUTINE history_init(history, history_section, iw)
60 : TYPE(history_type), INTENT(INOUT) :: history
61 : TYPE(section_vals_type), POINTER :: history_section
62 : INTEGER :: iw
63 :
64 3003 : ALLOCATE (history%entries(history_grow_unit))
65 3 : history%iw = iw
66 : CALL section_vals_val_get(history_section, "ENERGY_PRECISION", &
67 3 : r_val=history%E_precision)
68 : CALL section_vals_val_get(history_section, "FINGERPRINT_PRECISION", &
69 3 : r_val=history%FP_precision)
70 :
71 3 : IF (iw > 0) THEN
72 : WRITE (iw, '(A,T66,E15.3)') &
73 3 : " GLBOPT| History energy precision", history%E_precision
74 : WRITE (iw, '(A,T66,E15.3)') &
75 3 : " GLBOPT| History fingerprint precision", history%FP_precision
76 : END IF
77 3 : END SUBROUTINE history_init
78 :
79 : ! **************************************************************************************************
80 : !> \brief Calculates a fingerprint for a given configuration.
81 : !> \param Epot ...
82 : !> \param pos ...
83 : !> \return ...
84 : !> \author Ole Schuett
85 : ! **************************************************************************************************
86 60 : FUNCTION history_fingerprint(Epot, pos) RESULT(fp)
87 : REAL(KIND=dp), INTENT(IN) :: Epot
88 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: pos
89 : TYPE(history_fingerprint_type) :: fp
90 :
91 : INTEGER :: handle
92 30 : REAL(KIND=dp), DIMENSION(:), POINTER :: tmp
93 :
94 30 : CALL timeset("glbopt_history_fingerprint", handle)
95 :
96 30 : NULLIFY (tmp)
97 30 : fp%Epot = Epot
98 30 : CALL goedecker_fingerprint(pos, tmp)
99 :
100 : !copy pointer to allocatable
101 90 : ALLOCATE (fp%goedecker(SIZE(tmp)))
102 330 : fp%goedecker(:) = tmp
103 30 : DEALLOCATE (tmp)
104 :
105 30 : CALL timestop(handle)
106 30 : END FUNCTION history_fingerprint
107 :
108 : ! **************************************************************************************************
109 : !> \brief Helper routine for history_fingerprint.
110 : !> Calculates a fingerprint based on inter-atomic distances.
111 : !> \param pos ...
112 : !> \param res ...
113 : !> \author Stefan Goedecker
114 : ! **************************************************************************************************
115 30 : SUBROUTINE goedecker_fingerprint(pos, res)
116 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: pos
117 : REAL(KIND=dp), DIMENSION(:), POINTER :: res
118 :
119 : INTEGER :: i, info, j, N
120 : REAL(KIND=dp) :: d2, t
121 30 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: matrix, work
122 : REAL(KIND=dp), DIMENSION(3) :: d
123 :
124 30 : IF (ASSOCIATED(res)) CPABORT("goedecker_fingerprint: res already allocated")
125 30 : N = SIZE(pos)/3 ! number of atoms
126 :
127 180 : ALLOCATE (matrix(N, N), work(N, N))
128 330 : DO i = 1, N
129 300 : matrix(i, i) = 1.0
130 1680 : DO j = i + 1, N
131 5400 : d = pos(3*i - 2:3*i) - pos(3*j - 2:3*j)
132 5400 : d2 = SUM(d**2)
133 1350 : t = EXP(-0.5*d2)
134 1350 : matrix(i, j) = t
135 1650 : matrix(j, i) = t
136 : END DO
137 : END DO
138 90 : ALLOCATE (res(N))
139 : ! matrix values are garbage on exit because of jobz='N'
140 30 : CALL dsyev('N', 'U', N, matrix, N, res, work, N**2, info)
141 30 : IF (info /= 0) CPABORT("goedecker_fingerprint: DSYEV failed")
142 30 : END SUBROUTINE goedecker_fingerprint
143 :
144 : ! **************************************************************************************************
145 : !> \brief Checks if two given fingerprints match.
146 : !> \param history ...
147 : !> \param fp1 ...
148 : !> \param fp2 ...
149 : !> \return ...
150 : !> \author Ole Schuett
151 : ! **************************************************************************************************
152 11 : FUNCTION history_fingerprint_match(history, fp1, fp2) RESULT(res)
153 : TYPE(history_type), INTENT(IN) :: history
154 : TYPE(history_fingerprint_type), INTENT(IN) :: fp1, fp2
155 : LOGICAL :: res
156 :
157 : res = (ABS(fp1%Epot - fp2%Epot) < history%E_precision) .AND. &
158 11 : (fingerprint_distance(fp1, fp2) < history%fp_precision)
159 :
160 11 : END FUNCTION history_fingerprint_match
161 :
162 : ! **************************************************************************************************
163 : !> \brief Helper routine for history_fingerprint_match
164 : !> Calculates the distance between two given fingerprints.
165 : !> \param fp1 ...
166 : !> \param fp2 ...
167 : !> \return ...
168 : !> \author Stefan Goedecker
169 : ! **************************************************************************************************
170 12 : PURE FUNCTION fingerprint_distance(fp1, fp2) RESULT(res)
171 : TYPE(history_fingerprint_type), INTENT(IN) :: fp1, fp2
172 : REAL(KIND=dp) :: res
173 :
174 132 : res = SQRT(SUM((fp1%goedecker - fp2%goedecker)**2)/SIZE(fp1%goedecker))
175 12 : END FUNCTION fingerprint_distance
176 :
177 : ! **************************************************************************************************
178 : !> \brief Addes a new fingerprints to the history.
179 : !> Optionally, an abitrary id can be stored alongside the fingerprint.
180 : !> \param history ...
181 : !> \param fingerprint ...
182 : !> \param id ...
183 : !> \author Ole Schuett
184 : ! **************************************************************************************************
185 19 : SUBROUTINE history_add(history, fingerprint, id)
186 : TYPE(history_type), INTENT(INOUT) :: history
187 : TYPE(history_fingerprint_type), INTENT(IN) :: fingerprint
188 : INTEGER, INTENT(IN), OPTIONAL :: id
189 :
190 : INTEGER :: handle, i, k, n
191 19 : TYPE(history_entry_type), DIMENSION(:), POINTER :: tmp
192 :
193 19 : CALL timeset("glbopt_history_add", handle)
194 :
195 19 : n = SIZE(history%entries)
196 19 : IF (n == history%length) THEN
197 : ! grow history%entries array
198 0 : tmp => history%entries
199 0 : ALLOCATE (history%entries(n + history_grow_unit))
200 0 : history%entries(1:n) = tmp(:)
201 0 : DEALLOCATE (tmp)
202 0 : n = n + history_grow_unit
203 : END IF
204 :
205 19 : k = interpolation_search(history, fingerprint%Epot)
206 :
207 : !history%entries(k+1:) = history%entries(k:n-1)
208 : !Workaround for an XLF bug - pointer array copy does
209 : !not work correctly
210 18975 : DO i = n, k + 1, -1
211 18975 : history%entries(i) = history%entries(i - 1)
212 : END DO
213 :
214 19 : ALLOCATE (history%entries(k)%p)
215 19 : history%entries(k)%p = fingerprint
216 19 : IF (PRESENT(id)) THEN
217 12 : history%entries(k)%id = id
218 : END IF
219 19 : history%length = history%length + 1
220 :
221 : IF (debug) THEN
222 : ! check history for correct order
223 : DO k = 1, history%length
224 : !WRITE(*,*) "history: ", k, "Epot",history%entries(k)%p%Epot
225 : IF (k > 1) THEN
226 : IF (history%entries(k - 1)%p%Epot > history%entries(k)%p%Epot) THEN
227 : CPABORT("history_add: history in wrong order")
228 : END IF
229 : END IF
230 : END DO
231 : END IF
232 :
233 19 : CALL timestop(handle)
234 19 : END SUBROUTINE history_add
235 :
236 : ! **************************************************************************************************
237 : !> \brief Checks if a given fingerprints is contained in the history.
238 : !> \param history ...
239 : !> \param fingerprint ...
240 : !> \param found ...
241 : !> \param id ...
242 : !> \author Ole Schuett
243 : ! **************************************************************************************************
244 27 : SUBROUTINE history_lookup(history, fingerprint, found, id)
245 : TYPE(history_type), INTENT(IN) :: history
246 : TYPE(history_fingerprint_type), INTENT(IN) :: fingerprint
247 : LOGICAL, INTENT(OUT) :: found
248 : INTEGER, INTENT(OUT), OPTIONAL :: id
249 :
250 : INTEGER :: found_i, handle, i, k, k_max, k_min
251 : REAL(KIND=dp) :: best_match, dist, Epot
252 :
253 27 : CALL timeset("glbopt_history_lookup", handle)
254 :
255 27 : found = .FALSE.
256 27 : IF (PRESENT(id)) id = -1
257 27 : best_match = HUGE(1.0_dp)
258 :
259 27 : IF (history%length > 0) THEN
260 24 : Epot = fingerprint%Epot
261 24 : k = interpolation_search(history, fingerprint%Epot)
262 :
263 31 : DO k_min = k - 1, 1, -1
264 31 : IF (history%entries(k_min)%p%Epot < Epot - history%E_precision) EXIT
265 : END DO
266 :
267 26 : DO k_max = k, history%length
268 26 : IF (history%entries(k_max)%p%Epot > Epot + history%E_precision) EXIT
269 : END DO
270 :
271 24 : k_min = MAX(k_min + 1, 1)
272 24 : k_max = MIN(k_max - 1, history%length)
273 :
274 : IF (debug) found_i = -1
275 :
276 33 : DO i = k_min, k_max
277 9 : dist = fingerprint_distance(fingerprint, history%entries(i)%p)
278 : !WRITE(*,*) "entry ", i, " dist: ",dist
279 33 : IF (dist < history%fp_precision .AND. dist < best_match) THEN
280 8 : best_match = dist
281 8 : found = .TRUE.
282 8 : IF (PRESENT(id)) id = history%entries(i)%id
283 : IF (debug) found_i = i
284 : END IF
285 : END DO
286 :
287 : IF (debug) CALL verify_history_lookup(history, fingerprint, found_i)
288 : END IF
289 :
290 27 : CALL timestop(handle)
291 :
292 27 : END SUBROUTINE history_lookup
293 :
294 : ! **************************************************************************************************
295 : !> \brief Helper routine for history_lookup
296 : !> \param history ...
297 : !> \param Efind ...
298 : !> \return ...
299 : !> \author Ole Schuett
300 : ! **************************************************************************************************
301 43 : FUNCTION interpolation_search(history, Efind) RESULT(res)
302 : TYPE(history_type), INTENT(IN) :: history
303 : REAL(KIND=dp), INTENT(IN) :: Efind
304 : INTEGER :: res
305 :
306 : INTEGER :: high, low, mid
307 : REAL(KIND=dp) :: slope
308 :
309 43 : low = 1
310 43 : high = history%length
311 :
312 97 : DO WHILE (low < high)
313 : !linear interpolation
314 54 : slope = REAL(high - low, KIND=dp)/(history%entries(high)%p%Epot - history%entries(low)%p%Epot)
315 54 : mid = low + INT(slope*(Efind - history%entries(low)%p%Epot))
316 54 : mid = MIN(MAX(mid, low), high)
317 :
318 97 : IF (history%entries(mid)%p%Epot < Efind) THEN
319 26 : low = mid + 1
320 : ELSE
321 28 : high = mid - 1
322 : END IF
323 : END DO
324 :
325 43 : IF (0 < low .AND. low <= history%length) THEN
326 40 : IF (Efind > history%entries(low)%p%Epot) low = low + 1
327 : END IF
328 :
329 43 : res = low
330 43 : END FUNCTION interpolation_search
331 :
332 : ! **************************************************************************************************
333 : !> \brief Debugging routine, performs a slow (but robust) linear search.
334 : !> \param history ...
335 : !> \param fingerprint ...
336 : !> \param found_i_ref ...
337 : !> \author Ole Schuett
338 : ! **************************************************************************************************
339 0 : SUBROUTINE verify_history_lookup(history, fingerprint, found_i_ref)
340 : TYPE(history_type), INTENT(IN) :: history
341 : TYPE(history_fingerprint_type), INTENT(IN) :: fingerprint
342 : INTEGER, INTENT(IN) :: found_i_ref
343 :
344 : INTEGER :: found_i, i
345 : REAL(KIND=dp) :: best_fp_match, Epot_dist, fp_dist
346 :
347 0 : found_i = -1
348 0 : best_fp_match = HUGE(1.0_dp)
349 :
350 0 : DO i = 1, history%length
351 0 : Epot_dist = ABS(fingerprint%Epot - history%entries(i)%p%Epot)
352 0 : IF (Epot_dist > history%E_precision) CYCLE
353 0 : fp_dist = fingerprint_distance(fingerprint, history%entries(i)%p)
354 : !WRITE(*,*) "entry ", i, " dist: ",dist
355 0 : IF (fp_dist < history%fp_precision .AND. fp_dist < best_fp_match) THEN
356 0 : best_fp_match = fp_dist
357 0 : found_i = i
358 : END IF
359 : END DO
360 :
361 0 : IF (found_i /= found_i_ref) THEN
362 0 : WRITE (*, *) found_i, found_i_ref
363 0 : CPABORT("verify_history_lookup failed")
364 : END IF
365 :
366 0 : END SUBROUTINE verify_history_lookup
367 :
368 : ! **************************************************************************************************
369 : !> \brief Finalizes a history.
370 : !> \param history ...
371 : !> \author Ole Schuett
372 : ! **************************************************************************************************
373 3 : SUBROUTINE history_finalize(history)
374 : TYPE(history_type) :: history
375 :
376 : INTEGER :: i
377 :
378 22 : DO i = 1, history%length
379 22 : IF (ASSOCIATED(history%entries(i)%p)) THEN
380 19 : DEALLOCATE (history%entries(i)%p)
381 : END IF
382 : END DO
383 :
384 3 : DEALLOCATE (history%entries)
385 :
386 3 : END SUBROUTINE history_finalize
387 :
388 0 : END MODULE glbopt_history
|
