glbopt_history.F

Go to the documentation of this file.

!--------------------------------------------------------------------------------------------------!
!   CP2K: A general program to perform molecular dynamics simulations                              !
!   Copyright 2000-2024 CP2K developers group <https://cp2k.org>                                   !
!                                                                                                  !
!   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
!--------------------------------------------------------------------------------------------------!
 
! **************************************************************************************************
!> \brief History of minima, calculates, stores and compares fingerprints of minima.
!>        Used by Minima Hopping and Minima Crawling.
!> \author Ole Schuett
! **************************************************************************************************
MODULE glbopt_history
   USE input_section_types,             ONLY: section_vals_type,&
                                              section_vals_val_get
   USE kinds,                           ONLY: dp
#include "../base/base_uses.f90"
 
   IMPLICIT NONE
   PRIVATE
 
   TYPE history_fingerprint_type
      PRIVATE
      REAL(KIND=dp)                            :: epot = 0.0
      REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: goedecker
   END TYPE history_fingerprint_type
 
   TYPE history_entry_type
      TYPE(history_fingerprint_type), POINTER :: p => null()
      INTEGER                                 :: id = -1
   END TYPE history_entry_type
 
   TYPE history_type
      PRIVATE
      TYPE(history_entry_type), DIMENSION(:), POINTER :: entries => null()
      INTEGER                              :: length = 0
      INTEGER                              :: iw = -1
      REAL(KIND=dp)                        :: e_precision = 0.0
      REAL(KIND=dp)                        :: fp_precision = 0.0
   END TYPE history_type
 
   PUBLIC :: history_type, history_fingerprint_type
   PUBLIC :: history_init, history_finalize
   PUBLIC :: history_add, history_lookup
   PUBLIC :: history_fingerprint
   PUBLIC :: history_fingerprint_match
 
   LOGICAL, PARAMETER                     :: debug = .false.
   INTEGER, PARAMETER                     :: history_grow_unit = 1000
CONTAINS
 
! **************************************************************************************************
!> \brief Initializes a history.
!> \param history ...
!> \param history_section ...
!> \param iw ...
!> \author Ole Schuett
! **************************************************************************************************
   SUBROUTINE history_init(history, history_section, iw)
      TYPE(history_type), INTENT(INOUT)                  :: history
      TYPE(section_vals_type), POINTER                   :: history_section
      INTEGER                                            :: iw
 
      ALLOCATE (history%entries(history_grow_unit))
      history%iw = iw
      CALL section_vals_val_get(history_section, "ENERGY_PRECISION", &
                                r_val=history%E_precision)
      CALL section_vals_val_get(history_section, "FINGERPRINT_PRECISION", &
                                r_val=history%FP_precision)
 
      IF (iw > 0) THEN
         WRITE (iw, '(A,T66,E15.3)') &
            " GLBOPT| History energy precision", history%E_precision
         WRITE (iw, '(A,T66,E15.3)') &
            " GLBOPT| History fingerprint precision", history%FP_precision
      END IF
   END SUBROUTINE history_init
 
! **************************************************************************************************
!> \brief Calculates a fingerprint for a given configuration.
!> \param Epot ...
!> \param pos ...
!> \return ...
!> \author Ole Schuett
! **************************************************************************************************
   FUNCTION history_fingerprint(Epot, pos) RESULT(fp)
      REAL(kind=dp), INTENT(IN)                          :: epot
      REAL(kind=dp), DIMENSION(:), INTENT(IN)            :: pos
      TYPE(history_fingerprint_type)                     :: fp
 
      INTEGER                                            :: handle
      REAL(kind=dp), DIMENSION(:), POINTER               :: tmp
 
      CALL timeset("glbopt_history_fingerprint", handle)
 
      NULLIFY (tmp)
      fp%Epot = epot
      CALL goedecker_fingerprint(pos, tmp)
 
      !copy pointer to allocatable
      ALLOCATE (fp%goedecker(SIZE(tmp)))
      fp%goedecker(:) = tmp
      DEALLOCATE (tmp)
 
      CALL timestop(handle)
   END FUNCTION history_fingerprint
 
! **************************************************************************************************
!> \brief Helper routine for history_fingerprint.
!>        Calculates a fingerprint based on inter-atomic distances.
!> \param pos ...
!> \param res ...
!> \author Stefan Goedecker
! **************************************************************************************************
   SUBROUTINE goedecker_fingerprint(pos, res)
      REAL(kind=dp), DIMENSION(:), INTENT(IN)            :: pos
      REAL(kind=dp), DIMENSION(:), POINTER               :: res
 
      INTEGER                                            :: i, info, j, n
      REAL(kind=dp)                                      :: d2, t
      REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :)        :: matrix, work
      REAL(kind=dp), DIMENSION(3)                        :: d
 
      IF (ASSOCIATED(res)) cpabort("goedecker_fingerprint: res already allocated")
      n = SIZE(pos)/3 ! number of atoms
 
      ALLOCATE (matrix(n, n), work(n, n))
      DO i = 1, n
         matrix(i, i) = 1.0
         DO j = i + 1, n
            d = pos(3*i - 2:3*i) - pos(3*j - 2:3*j)
            d2 = sum(d**2)
            t = exp(-0.5*d2)
            matrix(i, j) = t
            matrix(j, i) = t
         END DO
      END DO
      !TODO: call dsyv through cp2k wrappers
      !TODO: do we need to store lower triangle of matrix?
      ALLOCATE (res(n))
      CALL dsyev('N', 'U', n, matrix, n, res, work, n**2, info)
      IF (info /= 0) cpabort("goedecker_fingerprint: DSYEV failed")
   END SUBROUTINE goedecker_fingerprint
 
! **************************************************************************************************
!> \brief Checks if two given fingerprints match.
!> \param history ...
!> \param fp1 ...
!> \param fp2 ...
!> \return ...
!> \author Ole Schuett
! **************************************************************************************************
   FUNCTION history_fingerprint_match(history, fp1, fp2) RESULT(res)
      TYPE(history_type), INTENT(IN)                     :: history
      TYPE(history_fingerprint_type), INTENT(IN)         :: fp1, fp2
      LOGICAL                                            :: res
 
      res = (abs(fp1%Epot - fp2%Epot) < history%E_precision) .AND. &
            (fingerprint_distance(fp1, fp2) < history%fp_precision)
 
   END FUNCTION history_fingerprint_match
 
! **************************************************************************************************
!> \brief Helper routine for history_fingerprint_match
!>        Calculates the distance between two given fingerprints.
!> \param fp1 ...
!> \param fp2 ...
!> \return ...
!> \author Stefan Goedecker
! **************************************************************************************************
   PURE FUNCTION fingerprint_distance(fp1, fp2) RESULT(res)
      TYPE(history_fingerprint_type), INTENT(IN)         :: fp1, fp2
      REAL(kind=dp)                                      :: res
 
      res = sqrt(sum((fp1%goedecker - fp2%goedecker)**2)/SIZE(fp1%goedecker))
   END FUNCTION fingerprint_distance
 
! **************************************************************************************************
!> \brief Addes a new fingerprints to the history.
!>        Optionally, an abitrary id can be stored alongside the fingerprint.
!> \param history ...
!> \param fingerprint ...
!> \param id ...
!> \author Ole Schuett
! **************************************************************************************************
   SUBROUTINE history_add(history, fingerprint, id)
      TYPE(history_type), INTENT(INOUT)                  :: history
      TYPE(history_fingerprint_type), INTENT(IN)         :: fingerprint
      INTEGER, INTENT(IN), OPTIONAL                      :: id
 
      INTEGER                                            :: handle, i, k, n
      TYPE(history_entry_type), DIMENSION(:), POINTER    :: tmp
 
      CALL timeset("glbopt_history_add", handle)
 
      n = SIZE(history%entries)
      IF (n == history%length) THEN
         ! grow history%entries array
         tmp => history%entries
         ALLOCATE (history%entries(n + history_grow_unit))
         history%entries(1:n) = tmp(:)
         DEALLOCATE (tmp)
         n = n + history_grow_unit
      END IF
 
      k = interpolation_search(history, fingerprint%Epot)
 
      !history%entries(k+1:) = history%entries(k:n-1)
      !Workaround for an XLF bug - pointer array copy does
      !not work correctly
      DO i = n, k + 1, -1
         history%entries(i) = history%entries(i - 1)
      END DO
 
      ALLOCATE (history%entries(k)%p)
      history%entries(k)%p = fingerprint
      IF (PRESENT(id)) &
         history%entries(k)%id = id
      history%length = history%length + 1
 
      IF (debug) THEN
         ! check history for correct order
         DO k = 1, history%length
            !WRITE(*,*) "history: ", k, "Epot",history%entries(k)%p%Epot
            IF (k > 1) THEN
               IF (history%entries(k - 1)%p%Epot > history%entries(k)%p%Epot) &
                  cpabort("history_add: history in wrong order")
            END IF
         END DO
      END IF
 
      CALL timestop(handle)
   END SUBROUTINE history_add
 
! **************************************************************************************************
!> \brief Checks if a given fingerprints is contained in the history.
!> \param history ...
!> \param fingerprint ...
!> \param found ...
!> \param id ...
!> \author Ole Schuett
! **************************************************************************************************
   SUBROUTINE history_lookup(history, fingerprint, found, id)
      TYPE(history_type), INTENT(IN)                     :: history
      TYPE(history_fingerprint_type), INTENT(IN)         :: fingerprint
      LOGICAL, INTENT(OUT)                               :: found
      INTEGER, INTENT(OUT), OPTIONAL                     :: id
 
      INTEGER                                            :: found_i, handle, i, k, k_max, k_min
      REAL(kind=dp)                                      :: best_match, dist, epot
 
      CALL timeset("glbopt_history_lookup", handle)
 
      found = .false.
      IF (PRESENT(id)) id = -1
      best_match = huge(1.0_dp)
 
      IF (history%length > 0) THEN
         epot = fingerprint%Epot
         k = interpolation_search(history, fingerprint%Epot)
 
         DO k_min = k - 1, 1, -1
            IF (history%entries(k_min)%p%Epot < epot - history%E_precision) EXIT
         END DO
 
         DO k_max = k, history%length
            IF (history%entries(k_max)%p%Epot > epot + history%E_precision) EXIT
         END DO
 
         k_min = max(k_min + 1, 1)
         k_max = min(k_max - 1, history%length)
 
         IF (debug) found_i = -1
 
         DO i = k_min, k_max
            dist = fingerprint_distance(fingerprint, history%entries(i)%p)
            !WRITE(*,*) "entry ", i, " dist: ",dist
            IF (dist < history%fp_precision .AND. dist < best_match) THEN
               best_match = dist
               found = .true.
               IF (PRESENT(id)) id = history%entries(i)%id
               IF (debug) found_i = i
            END IF
         END DO
 
         IF (debug) CALL verify_history_lookup(history, fingerprint, found_i)
      END IF
 
      CALL timestop(handle)
 
   END SUBROUTINE history_lookup
 
! **************************************************************************************************
!> \brief Helper routine for history_lookup
!> \param history ...
!> \param Efind ...
!> \return ...
!> \author Ole Schuett
! **************************************************************************************************
   FUNCTION interpolation_search(history, Efind) RESULT(res)
      TYPE(history_type), INTENT(IN)                     :: history
      REAL(kind=dp), INTENT(IN)                          :: efind
      INTEGER                                            :: res
 
      INTEGER                                            :: high, low, mid
      REAL(kind=dp)                                      :: slope
 
      low = 1
      high = history%length
 
      DO WHILE (low < high)
         !linear interpolation
         slope = real(high - low, kind=dp)/(history%entries(high)%p%Epot - history%entries(low)%p%Epot)
         mid = low + int(slope*(efind - history%entries(low)%p%Epot))
         mid = min(max(mid, low), high)
 
         IF (history%entries(mid)%p%Epot < efind) THEN
            low = mid + 1
         ELSE
            high = mid - 1
         END IF
      END DO
 
      IF (0 < low .AND. low <= history%length) THEN
         IF (efind > history%entries(low)%p%Epot) low = low + 1
      END IF
 
      res = low
   END FUNCTION interpolation_search
 
! **************************************************************************************************
!> \brief Debugging routine, performs a slow (but robust) linear search.
!> \param history ...
!> \param fingerprint ...
!> \param found_i_ref ...
!> \author Ole Schuett
! **************************************************************************************************
   SUBROUTINE verify_history_lookup(history, fingerprint, found_i_ref)
      TYPE(history_type), INTENT(IN)                     :: history
      TYPE(history_fingerprint_type), INTENT(IN)         :: fingerprint
      INTEGER, INTENT(IN)                                :: found_i_ref
 
      INTEGER                                            :: found_i, i
      REAL(kind=dp)                                      :: best_fp_match, epot_dist, fp_dist
 
      found_i = -1
      best_fp_match = huge(1.0_dp)
 
      DO i = 1, history%length
         epot_dist = abs(fingerprint%Epot - history%entries(i)%p%Epot)
         IF (epot_dist > history%E_precision) cycle
         fp_dist = fingerprint_distance(fingerprint, history%entries(i)%p)
         !WRITE(*,*) "entry ", i, " dist: ",dist
         IF (fp_dist < history%fp_precision .AND. fp_dist < best_fp_match) THEN
            best_fp_match = fp_dist
            found_i = i
         END IF
      END DO
 
      IF (found_i /= found_i_ref) THEN
         WRITE (*, *) found_i, found_i_ref
         cpabort("verify_history_lookup failed")
      END IF
 
   END SUBROUTINE verify_history_lookup
 
! **************************************************************************************************
!> \brief Finalizes a history.
!> \param history ...
!> \author Ole Schuett
! **************************************************************************************************
   SUBROUTINE history_finalize(history)
      TYPE(history_type)                                 :: history
 
      INTEGER                                            :: i
 
      DO i = 1, history%length
         IF (ASSOCIATED(history%entries(i)%p)) &
            DEALLOCATE (history%entries(i)%p)
      END DO
 
      DEALLOCATE (history%entries)
 
   END SUBROUTINE history_finalize
 
END MODULE glbopt_history