Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 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 102 : 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 51 : REAL(KIND=dp), DIMENSION(:), POINTER :: tmp
93 :
94 51 : CALL timeset("glbopt_history_fingerprint", handle)
95 :
96 51 : NULLIFY (tmp)
97 51 : fp%Epot = Epot
98 51 : CALL goedecker_fingerprint(pos, tmp)
99 :
100 : !copy pointer to allocatable
101 153 : ALLOCATE (fp%goedecker(SIZE(tmp)))
102 561 : fp%goedecker(:) = tmp
103 51 : DEALLOCATE (tmp)
104 :
105 51 : CALL timestop(handle)
106 51 : 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 51 : 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 51 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: matrix, work
122 : REAL(KIND=dp), DIMENSION(3) :: d
123 :
124 51 : IF (ASSOCIATED(res)) CPABORT("goedecker_fingerprint: res already allocated")
125 51 : N = SIZE(pos)/3 ! number of atoms
126 :
127 306 : ALLOCATE (matrix(N, N), work(N, N))
128 561 : DO i = 1, N
129 510 : matrix(i, i) = 1.0
130 2856 : DO j = i + 1, N
131 9180 : d = pos(3*i - 2:3*i) - pos(3*j - 2:3*j)
132 9180 : d2 = SUM(d**2)
133 2295 : t = EXP(-0.5*d2)
134 2295 : matrix(i, j) = t
135 2805 : matrix(j, i) = t
136 : END DO
137 : END DO
138 153 : ALLOCATE (res(N))
139 : ! matrix values are garbage on exit because of jobz='N'
140 51 : CALL dsyev('N', 'U', N, matrix, N, res, work, N**2, info)
141 51 : IF (info /= 0) CPABORT("goedecker_fingerprint: DSYEV failed")
142 51 : 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 34 : PURE FUNCTION fingerprint_distance(fp1, fp2) RESULT(res)
171 : TYPE(history_fingerprint_type), INTENT(IN) :: fp1, fp2
172 : REAL(KIND=dp) :: res
173 :
174 374 : res = SQRT(SUM((fp1%goedecker - fp2%goedecker)**2)/SIZE(fp1%goedecker))
175 34 : 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 18 : 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 18 : TYPE(history_entry_type), DIMENSION(:), POINTER :: tmp
192 :
193 18 : CALL timeset("glbopt_history_add", handle)
194 :
195 18 : n = SIZE(history%entries)
196 18 : 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 18 : 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 17975 : DO i = n, k + 1, -1
211 17975 : history%entries(i) = history%entries(i - 1)
212 : END DO
213 :
214 18 : ALLOCATE (history%entries(k)%p)
215 18 : history%entries(k)%p = fingerprint
216 18 : IF (PRESENT(id)) &
217 11 : history%entries(k)%id = id
218 18 : history%length = history%length + 1
219 :
220 : IF (debug) THEN
221 : ! check history for correct order
222 : DO k = 1, history%length
223 : !WRITE(*,*) "history: ", k, "Epot",history%entries(k)%p%Epot
224 : IF (k > 1) THEN
225 : IF (history%entries(k - 1)%p%Epot > history%entries(k)%p%Epot) &
226 : CPABORT("history_add: history in wrong order")
227 : END IF
228 : END DO
229 : END IF
230 :
231 18 : CALL timestop(handle)
232 18 : END SUBROUTINE history_add
233 :
234 : ! **************************************************************************************************
235 : !> \brief Checks if a given fingerprints is contained in the history.
236 : !> \param history ...
237 : !> \param fingerprint ...
238 : !> \param found ...
239 : !> \param id ...
240 : !> \author Ole Schuett
241 : ! **************************************************************************************************
242 48 : SUBROUTINE history_lookup(history, fingerprint, found, id)
243 : TYPE(history_type), INTENT(IN) :: history
244 : TYPE(history_fingerprint_type), INTENT(IN) :: fingerprint
245 : LOGICAL, INTENT(OUT) :: found
246 : INTEGER, INTENT(OUT), OPTIONAL :: id
247 :
248 : INTEGER :: found_i, handle, i, k, k_max, k_min
249 : REAL(KIND=dp) :: best_match, dist, Epot
250 :
251 48 : CALL timeset("glbopt_history_lookup", handle)
252 :
253 48 : found = .FALSE.
254 48 : IF (PRESENT(id)) id = -1
255 48 : best_match = HUGE(1.0_dp)
256 :
257 48 : IF (history%length > 0) THEN
258 45 : Epot = fingerprint%Epot
259 45 : k = interpolation_search(history, fingerprint%Epot)
260 :
261 71 : DO k_min = k - 1, 1, -1
262 71 : IF (history%entries(k_min)%p%Epot < Epot - history%E_precision) EXIT
263 : END DO
264 :
265 50 : DO k_max = k, history%length
266 50 : IF (history%entries(k_max)%p%Epot > Epot + history%E_precision) EXIT
267 : END DO
268 :
269 45 : k_min = MAX(k_min + 1, 1)
270 45 : k_max = MIN(k_max - 1, history%length)
271 :
272 : IF (debug) found_i = -1
273 :
274 76 : DO i = k_min, k_max
275 31 : dist = fingerprint_distance(fingerprint, history%entries(i)%p)
276 : !WRITE(*,*) "entry ", i, " dist: ",dist
277 76 : IF (dist < history%fp_precision .AND. dist < best_match) THEN
278 30 : best_match = dist
279 30 : found = .TRUE.
280 30 : IF (PRESENT(id)) id = history%entries(i)%id
281 : IF (debug) found_i = i
282 : END IF
283 : END DO
284 :
285 : IF (debug) CALL verify_history_lookup(history, fingerprint, found_i)
286 : END IF
287 :
288 48 : CALL timestop(handle)
289 :
290 48 : END SUBROUTINE history_lookup
291 :
292 : ! **************************************************************************************************
293 : !> \brief Helper routine for history_lookup
294 : !> \param history ...
295 : !> \param Efind ...
296 : !> \return ...
297 : !> \author Ole Schuett
298 : ! **************************************************************************************************
299 63 : FUNCTION interpolation_search(history, Efind) RESULT(res)
300 : TYPE(history_type), INTENT(IN) :: history
301 : REAL(KIND=dp), INTENT(IN) :: Efind
302 : INTEGER :: res
303 :
304 : INTEGER :: high, low, mid
305 : REAL(KIND=dp) :: slope
306 :
307 63 : low = 1
308 63 : high = history%length
309 :
310 157 : DO WHILE (low < high)
311 : !linear interpolation
312 94 : slope = REAL(high - low, KIND=dp)/(history%entries(high)%p%Epot - history%entries(low)%p%Epot)
313 94 : mid = low + INT(slope*(Efind - history%entries(low)%p%Epot))
314 94 : mid = MIN(MAX(mid, low), high)
315 :
316 157 : IF (history%entries(mid)%p%Epot < Efind) THEN
317 46 : low = mid + 1
318 : ELSE
319 48 : high = mid - 1
320 : END IF
321 : END DO
322 :
323 63 : IF (0 < low .AND. low <= history%length) THEN
324 60 : IF (Efind > history%entries(low)%p%Epot) low = low + 1
325 : END IF
326 :
327 63 : res = low
328 63 : END FUNCTION interpolation_search
329 :
330 : ! **************************************************************************************************
331 : !> \brief Debugging routine, performs a slow (but robust) linear search.
332 : !> \param history ...
333 : !> \param fingerprint ...
334 : !> \param found_i_ref ...
335 : !> \author Ole Schuett
336 : ! **************************************************************************************************
337 0 : SUBROUTINE verify_history_lookup(history, fingerprint, found_i_ref)
338 : TYPE(history_type), INTENT(IN) :: history
339 : TYPE(history_fingerprint_type), INTENT(IN) :: fingerprint
340 : INTEGER, INTENT(IN) :: found_i_ref
341 :
342 : INTEGER :: found_i, i
343 : REAL(KIND=dp) :: best_fp_match, Epot_dist, fp_dist
344 :
345 0 : found_i = -1
346 0 : best_fp_match = HUGE(1.0_dp)
347 :
348 0 : DO i = 1, history%length
349 0 : Epot_dist = ABS(fingerprint%Epot - history%entries(i)%p%Epot)
350 0 : IF (Epot_dist > history%E_precision) CYCLE
351 0 : fp_dist = fingerprint_distance(fingerprint, history%entries(i)%p)
352 : !WRITE(*,*) "entry ", i, " dist: ",dist
353 0 : IF (fp_dist < history%fp_precision .AND. fp_dist < best_fp_match) THEN
354 0 : best_fp_match = fp_dist
355 0 : found_i = i
356 : END IF
357 : END DO
358 :
359 0 : IF (found_i /= found_i_ref) THEN
360 0 : WRITE (*, *) found_i, found_i_ref
361 0 : CPABORT("verify_history_lookup failed")
362 : END IF
363 :
364 0 : END SUBROUTINE verify_history_lookup
365 :
366 : ! **************************************************************************************************
367 : !> \brief Finalizes a history.
368 : !> \param history ...
369 : !> \author Ole Schuett
370 : ! **************************************************************************************************
371 3 : SUBROUTINE history_finalize(history)
372 : TYPE(history_type) :: history
373 :
374 : INTEGER :: i
375 :
376 21 : DO i = 1, history%length
377 18 : IF (ASSOCIATED(history%entries(i)%p)) &
378 21 : DEALLOCATE (history%entries(i)%p)
379 : END DO
380 :
381 3 : DEALLOCATE (history%entries)
382 :
383 3 : END SUBROUTINE history_finalize
384 :
385 0 : END MODULE glbopt_history
|
