Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 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 126 : 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 63 : REAL(KIND=dp), DIMENSION(:), POINTER :: tmp
93 :
94 63 : CALL timeset("glbopt_history_fingerprint", handle)
95 :
96 63 : NULLIFY (tmp)
97 63 : fp%Epot = Epot
98 63 : CALL goedecker_fingerprint(pos, tmp)
99 :
100 : !copy pointer to allocatable
101 189 : ALLOCATE (fp%goedecker(SIZE(tmp)))
102 693 : fp%goedecker(:) = tmp
103 63 : DEALLOCATE (tmp)
104 :
105 63 : CALL timestop(handle)
106 63 : 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 63 : 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 63 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: matrix, work
122 : REAL(KIND=dp), DIMENSION(3) :: d
123 :
124 63 : IF (ASSOCIATED(res)) CPABORT("goedecker_fingerprint: res already allocated")
125 63 : N = SIZE(pos)/3 ! number of atoms
126 :
127 441 : ALLOCATE (matrix(N, N), work(N, N))
128 693 : DO i = 1, N
129 630 : matrix(i, i) = 1.0
130 3528 : DO j = i + 1, N
131 11340 : d = pos(3*i - 2:3*i) - pos(3*j - 2:3*j)
132 11340 : d2 = SUM(d**2)
133 2835 : t = EXP(-0.5*d2)
134 2835 : matrix(i, j) = t
135 3465 : matrix(j, i) = t
136 : END DO
137 : END DO
138 : !TODO: call dsyv through cp2k wrappers
139 : !TODO: do we need to store lower triangle of matrix?
140 189 : ALLOCATE (res(N))
141 63 : CALL DSYEV('N', 'U', N, matrix, N, res, work, N**2, info)
142 63 : IF (info /= 0) CPABORT("goedecker_fingerprint: DSYEV failed")
143 63 : END SUBROUTINE goedecker_fingerprint
144 :
145 : ! **************************************************************************************************
146 : !> \brief Checks if two given fingerprints match.
147 : !> \param history ...
148 : !> \param fp1 ...
149 : !> \param fp2 ...
150 : !> \return ...
151 : !> \author Ole Schuett
152 : ! **************************************************************************************************
153 11 : FUNCTION history_fingerprint_match(history, fp1, fp2) RESULT(res)
154 : TYPE(history_type), INTENT(IN) :: history
155 : TYPE(history_fingerprint_type), INTENT(IN) :: fp1, fp2
156 : LOGICAL :: res
157 :
158 : res = (ABS(fp1%Epot - fp2%Epot) < history%E_precision) .AND. &
159 11 : (fingerprint_distance(fp1, fp2) < history%fp_precision)
160 :
161 11 : END FUNCTION history_fingerprint_match
162 :
163 : ! **************************************************************************************************
164 : !> \brief Helper routine for history_fingerprint_match
165 : !> Calculates the distance between two given fingerprints.
166 : !> \param fp1 ...
167 : !> \param fp2 ...
168 : !> \return ...
169 : !> \author Stefan Goedecker
170 : ! **************************************************************************************************
171 45 : PURE FUNCTION fingerprint_distance(fp1, fp2) RESULT(res)
172 : TYPE(history_fingerprint_type), INTENT(IN) :: fp1, fp2
173 : REAL(KIND=dp) :: res
174 :
175 495 : res = SQRT(SUM((fp1%goedecker - fp2%goedecker)**2)/SIZE(fp1%goedecker))
176 45 : END FUNCTION fingerprint_distance
177 :
178 : ! **************************************************************************************************
179 : !> \brief Addes a new fingerprints to the history.
180 : !> Optionally, an abitrary id can be stored alongside the fingerprint.
181 : !> \param history ...
182 : !> \param fingerprint ...
183 : !> \param id ...
184 : !> \author Ole Schuett
185 : ! **************************************************************************************************
186 19 : SUBROUTINE history_add(history, fingerprint, id)
187 : TYPE(history_type), INTENT(INOUT) :: history
188 : TYPE(history_fingerprint_type), INTENT(IN) :: fingerprint
189 : INTEGER, INTENT(IN), OPTIONAL :: id
190 :
191 : INTEGER :: handle, i, k, n
192 19 : TYPE(history_entry_type), DIMENSION(:), POINTER :: tmp
193 :
194 19 : CALL timeset("glbopt_history_add", handle)
195 :
196 19 : n = SIZE(history%entries)
197 19 : IF (n == history%length) THEN
198 : ! grow history%entries array
199 0 : tmp => history%entries
200 0 : ALLOCATE (history%entries(n + history_grow_unit))
201 0 : history%entries(1:n) = tmp(:)
202 0 : DEALLOCATE (tmp)
203 0 : n = n + history_grow_unit
204 : END IF
205 :
206 19 : k = interpolation_search(history, fingerprint%Epot)
207 :
208 : !history%entries(k+1:) = history%entries(k:n-1)
209 : !Workaround for an XLF bug - pointer array copy does
210 : !not work correctly
211 18968 : DO i = n, k + 1, -1
212 18968 : history%entries(i) = history%entries(i - 1)
213 : END DO
214 :
215 19 : ALLOCATE (history%entries(k)%p)
216 19 : history%entries(k)%p = fingerprint
217 19 : IF (PRESENT(id)) &
218 12 : history%entries(k)%id = id
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) &
227 : CPABORT("history_add: history in wrong order")
228 : END IF
229 : END DO
230 : END IF
231 :
232 19 : CALL timestop(handle)
233 19 : END SUBROUTINE history_add
234 :
235 : ! **************************************************************************************************
236 : !> \brief Checks if a given fingerprints is contained in the history.
237 : !> \param history ...
238 : !> \param fingerprint ...
239 : !> \param found ...
240 : !> \param id ...
241 : !> \author Ole Schuett
242 : ! **************************************************************************************************
243 60 : SUBROUTINE history_lookup(history, fingerprint, found, id)
244 : TYPE(history_type), INTENT(IN) :: history
245 : TYPE(history_fingerprint_type), INTENT(IN) :: fingerprint
246 : LOGICAL, INTENT(OUT) :: found
247 : INTEGER, INTENT(OUT), OPTIONAL :: id
248 :
249 : INTEGER :: found_i, handle, i, k, k_max, k_min
250 : REAL(KIND=dp) :: best_match, dist, Epot
251 :
252 60 : CALL timeset("glbopt_history_lookup", handle)
253 :
254 60 : found = .FALSE.
255 60 : IF (PRESENT(id)) id = -1
256 60 : best_match = HUGE(1.0_dp)
257 :
258 60 : IF (history%length > 0) THEN
259 57 : Epot = fingerprint%Epot
260 57 : k = interpolation_search(history, fingerprint%Epot)
261 :
262 77 : DO k_min = k - 1, 1, -1
263 77 : IF (history%entries(k_min)%p%Epot < Epot - history%E_precision) EXIT
264 : END DO
265 :
266 79 : DO k_max = k, history%length
267 79 : IF (history%entries(k_max)%p%Epot > Epot + history%E_precision) EXIT
268 : END DO
269 :
270 57 : k_min = MAX(k_min + 1, 1)
271 57 : k_max = MIN(k_max - 1, history%length)
272 :
273 : IF (debug) found_i = -1
274 :
275 99 : DO i = k_min, k_max
276 42 : dist = fingerprint_distance(fingerprint, history%entries(i)%p)
277 : !WRITE(*,*) "entry ", i, " dist: ",dist
278 99 : IF (dist < history%fp_precision .AND. dist < best_match) THEN
279 41 : best_match = dist
280 41 : found = .TRUE.
281 41 : IF (PRESENT(id)) id = history%entries(i)%id
282 : IF (debug) found_i = i
283 : END IF
284 : END DO
285 :
286 : IF (debug) CALL verify_history_lookup(history, fingerprint, found_i)
287 : END IF
288 :
289 60 : CALL timestop(handle)
290 :
291 60 : END SUBROUTINE history_lookup
292 :
293 : ! **************************************************************************************************
294 : !> \brief Helper routine for history_lookup
295 : !> \param history ...
296 : !> \param Efind ...
297 : !> \return ...
298 : !> \author Ole Schuett
299 : ! **************************************************************************************************
300 76 : FUNCTION interpolation_search(history, Efind) RESULT(res)
301 : TYPE(history_type), INTENT(IN) :: history
302 : REAL(KIND=dp), INTENT(IN) :: Efind
303 : INTEGER :: res
304 :
305 : INTEGER :: high, low, mid
306 : REAL(KIND=dp) :: slope
307 :
308 76 : low = 1
309 76 : high = history%length
310 :
311 193 : DO WHILE (low < high)
312 : !linear interpolation
313 117 : slope = REAL(high - low, KIND=dp)/(history%entries(high)%p%Epot - history%entries(low)%p%Epot)
314 117 : mid = low + INT(slope*(Efind - history%entries(low)%p%Epot))
315 117 : mid = MIN(MAX(mid, low), high)
316 :
317 193 : IF (history%entries(mid)%p%Epot < Efind) THEN
318 58 : low = mid + 1
319 : ELSE
320 59 : high = mid - 1
321 : END IF
322 : END DO
323 :
324 76 : IF (0 < low .AND. low <= history%length) THEN
325 71 : IF (Efind > history%entries(low)%p%Epot) low = low + 1
326 : END IF
327 :
328 76 : res = low
329 76 : END FUNCTION interpolation_search
330 :
331 : ! **************************************************************************************************
332 : !> \brief Debugging routine, performs a slow (but robust) linear search.
333 : !> \param history ...
334 : !> \param fingerprint ...
335 : !> \param found_i_ref ...
336 : !> \author Ole Schuett
337 : ! **************************************************************************************************
338 0 : SUBROUTINE verify_history_lookup(history, fingerprint, found_i_ref)
339 : TYPE(history_type), INTENT(IN) :: history
340 : TYPE(history_fingerprint_type), INTENT(IN) :: fingerprint
341 : INTEGER, INTENT(IN) :: found_i_ref
342 :
343 : INTEGER :: found_i, i
344 : REAL(KIND=dp) :: best_fp_match, Epot_dist, fp_dist
345 :
346 0 : found_i = -1
347 0 : best_fp_match = HUGE(1.0_dp)
348 :
349 0 : DO i = 1, history%length
350 0 : Epot_dist = ABS(fingerprint%Epot - history%entries(i)%p%Epot)
351 0 : IF (Epot_dist > history%E_precision) CYCLE
352 0 : fp_dist = fingerprint_distance(fingerprint, history%entries(i)%p)
353 : !WRITE(*,*) "entry ", i, " dist: ",dist
354 0 : IF (fp_dist < history%fp_precision .AND. fp_dist < best_fp_match) THEN
355 0 : best_fp_match = fp_dist
356 0 : found_i = i
357 : END IF
358 : END DO
359 :
360 0 : IF (found_i /= found_i_ref) THEN
361 0 : WRITE (*, *) found_i, found_i_ref
362 0 : CPABORT("verify_history_lookup failed")
363 : END IF
364 :
365 0 : END SUBROUTINE verify_history_lookup
366 :
367 : ! **************************************************************************************************
368 : !> \brief Finalizes a history.
369 : !> \param history ...
370 : !> \author Ole Schuett
371 : ! **************************************************************************************************
372 3 : SUBROUTINE history_finalize(history)
373 : TYPE(history_type) :: history
374 :
375 : INTEGER :: i
376 :
377 22 : DO i = 1, history%length
378 19 : IF (ASSOCIATED(history%entries(i)%p)) &
379 22 : DEALLOCATE (history%entries(i)%p)
380 : END DO
381 :
382 3 : DEALLOCATE (history%entries)
383 :
384 3 : END SUBROUTINE history_finalize
385 :
386 0 : END MODULE glbopt_history
|
