graphcon.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 uses a combination of graphs and hashing to determine if two molecules
!>      are topologically equivalent, and if so, finds the one by one mapping
!> \note
!>      the graph isomorphism being solved is a computationally hard one
!>      and can not be solved in polynomial time in the general case
!>      http://mathworld.wolfram.com/IsomorphicGraphs.html
!>      the problem arises if many atoms are topologically equivalent
!>      the current algorithm is able to solve the problem for benzene (C6H6)
!>      but not for a fullerene (C60). Large systems are not really a problem (JAC).
!>      as almost all atoms are topologically unique.
!> \par History
!>      09.2006 [Joost VandeVondele]
!> \author Joost VandeVondele
! **************************************************************************************************
MODULE graphcon
 
   USE util,                            ONLY: sort
#include "./base/base_uses.f90"
 
   IMPLICIT NONE
 
   PRIVATE
   PUBLIC :: vertex, graph_type, reorder_graph, hash_molecule
 
   ! a molecule is an array of vertices, each vertex has a kind
   ! and a list of edges (bonds).
   ! (the number is the index of the other vertex in the array that builds the molecule)
! **************************************************************************************************
   TYPE graph_type
      TYPE(vertex), POINTER, DIMENSION(:) :: graph
   END TYPE graph_type
 
! **************************************************************************************************
   TYPE vertex
      INTEGER :: kind
      INTEGER, POINTER, DIMENSION(:) :: bonds
   END TYPE vertex
 
! **************************************************************************************************
   TYPE class
      INTEGER, DIMENSION(:), POINTER :: reference
      INTEGER, DIMENSION(:), POINTER :: unordered
      INTEGER :: kind
      INTEGER :: nele
      LOGICAL :: first
      INTEGER, DIMENSION(:), POINTER :: order
      INTEGER, DIMENSION(:), POINTER :: q
   END TYPE class
 
! **************************************************************************************************
   TYPE superclass
      INTEGER :: nele
      INTEGER, DIMENSION(:), POINTER :: classes
   END TYPE
 
CONTAINS
 
! **************************************************************************************************
!> \brief hashes a molecule to a number. Molecules that are the (topologically) the same
!>      have the same hash. However, there is a small chance that molecules with the same hash
!>      are different
!> \param reference IN  : molecule with atomic kinds and bonds
!> \param kind_ref OUT : an atomic hash which is the same for topologically equivalent atoms
!> \param hash OUT : a hash which is the same for topologically equivalent molecules
!> \par History
!>      09.2006 created [Joost VandeVondele]
!> \note
!>      Although relatively fast in general, might be quadratic with molecule size for
!>      some systems (e.g. linear alkanes)
! **************************************************************************************************
   SUBROUTINE hash_molecule(reference, kind_ref, hash)
      TYPE(vertex), DIMENSION(:), INTENT(IN)             :: reference
      INTEGER, DIMENSION(:), INTENT(OUT)                 :: kind_ref
      INTEGER, INTENT(OUT)                               :: hash
 
      INTEGER                                            :: i, ihash, n, nclasses, nclasses_old, &
                                                            old_class
      INTEGER, ALLOCATABLE, DIMENSION(:)                 :: index, kind_new
 
      n = SIZE(kind_ref)
      ALLOCATE (kind_new(n), index(n))
      kind_ref = reference%kind
      nclasses_old = 0
      DO ihash = 1, n
         ! generate a hash based on the the kind of each atom and the kind of its bonded atoms
         DO i = 1, n
            kind_new(i) = hash_kind(kind_ref(i), kind_ref(reference(i)%bonds))
         END DO
         kind_ref = kind_new
         ! find the number of equivalent atoms
         CALL sort(kind_new, n, index)
         nclasses = 1
         old_class = kind_new(1)
         DO i = 2, n
            IF (kind_new(i) .NE. old_class) THEN
               nclasses = nclasses + 1
               old_class = kind_new(i)
            END IF
         END DO
         ! if we have not generated new classes, we have presumably found all equivalence classes
         IF (nclasses == nclasses_old) EXIT
         nclasses_old = nclasses
         ! write(*,*) "Classes",Ihash, Nclasses
      END DO
      ! hash (sorted) kinds to a molecular hash
      hash = joaat_hash_i(kind_new)
      DEALLOCATE (kind_new, index)
   END SUBROUTINE hash_molecule
 
! **************************************************************************************************
!> \brief If two molecules are topologically the same, finds the ordering that maps
!>      the unordered one on the ordered one.
!> \param reference molecular description (see type definition)
!> \param unordered molecular description (see type definition)
!> \param order the mapping reference=order(unordred) if matches=.TRUE.
!>                             undefined if matches=.FALSE.
!> \param matches .TRUE. = the ordering was found
!> \par History
!>      09.2006 created [Joost VandeVondele]
!> \note
!>      See not at the top of the file about why this algorithm might consider
!>      molecules with a large number of equivalent atoms as different
!>      despite the fact that an ordering could exist for which they are the same
! **************************************************************************************************
   SUBROUTINE reorder_graph(reference, unordered, order, matches)
      TYPE(vertex), DIMENSION(:), INTENT(IN)             :: reference, unordered
      INTEGER, DIMENSION(:), INTENT(OUT)                 :: order
      LOGICAL, INTENT(OUT)                               :: matches
 
      INTEGER, PARAMETER                                 :: max_tries = 1000000
 
      INTEGER                                            :: hash_re, hash_un, i, iclass, iele, &
                                                            isuperclass, itries, j, n, nclasses, &
                                                            nele, old_class
      INTEGER, ALLOCATABLE, DIMENSION(:)                 :: class_of_atom, index_ref, index_un, &
                                                            kind_ref, kind_ref_ordered, kind_un, &
                                                            kind_un_ordered, superclass_of_atom
      TYPE(class), ALLOCATABLE, DIMENSION(:)             :: classes
      TYPE(superclass), ALLOCATABLE, DIMENSION(:)        :: superclasses
 
! allows for worst case matching of two benzenes ... (6!)*(6!)/6=86400
! with some margin for other molecules
! molecules with no symmetry e.g. JAC need less than 500 tries
! catch the cases where the molecules are trivially different
 
      IF (SIZE(reference) .NE. SIZE(unordered)) THEN
         matches = .false.
         RETURN
      END IF
 
      ! catch the case where the molecules are already in the right order
      n = SIZE(order)
      order = (/(i, i=1, n)/)
      IF (matrix_equal(reference, unordered, order)) THEN
         matches = .true.
         RETURN
      END IF
 
      ! determine the kind of each atom, and the hash of the whole molecule
      ALLOCATE (kind_ref(n), kind_un(n), index_ref(n), index_un(n), &
                kind_ref_ordered(n), kind_un_ordered(n), &
                class_of_atom(n), superclass_of_atom(n))
      CALL hash_molecule(reference, kind_ref, hash_re)
      CALL hash_molecule(unordered, kind_un, hash_un)
      IF (hash_re .NE. hash_un) THEN
         matches = .false.
         RETURN
      END IF
 
      ! generate the classes of equivalent atoms, i.e. the groups of atoms of the same topological kind
      kind_ref_ordered(:) = kind_ref
      CALL sort(kind_ref_ordered, n, index_ref)
      kind_un_ordered(:) = kind_un
      CALL sort(kind_un_ordered, n, index_un)
      IF (any(kind_ref_ordered .NE. kind_un_ordered)) THEN
         matches = .false.
         RETURN
      END IF
 
      ! count different classes, assign their kinds, and the number of elements
      nclasses = 1
      old_class = kind_ref_ordered(1)
      DO i = 2, n
         IF (kind_ref_ordered(i) .NE. old_class) THEN
            nclasses = nclasses + 1
            old_class = kind_ref_ordered(i)
         END IF
      END DO
      ALLOCATE (classes(nclasses))
      classes(1)%kind = kind_ref_ordered(1)
      nclasses = 1
      classes(1)%Nele = 1
      DO i = 2, n
         IF (kind_ref_ordered(i) .NE. classes(nclasses)%kind) THEN
            nclasses = nclasses + 1
            classes(nclasses)%kind = kind_ref_ordered(i)
            classes(nclasses)%Nele = 1
         ELSE
            classes(nclasses)%Nele = classes(nclasses)%Nele + 1
         END IF
      END DO
 
      ! assign the atoms to their classes
      iele = 0
      DO i = 1, nclasses
         nele = classes(i)%Nele
         ALLOCATE (classes(i)%reference(nele))
         ALLOCATE (classes(i)%unordered(nele))
         DO j = 1, nele
            iele = iele + 1
            classes(i)%reference(j) = index_ref(iele)
            classes(i)%unordered(j) = index_un(iele)
         END DO
         class_of_atom(classes(i)%reference) = i
         ALLOCATE (classes(i)%order(nele))
         ALLOCATE (classes(i)%q(nele))
         classes(i)%order = (/(j, j=1, nele)/)
         classes(i)%first = .true.
      END DO
 
      ! find which groups of classes (superclasses) that can be solved independently.
      ! only classes with more than one element that are connected need to be reordered simultaniously
 
      ! find these connected components in a recursive way
      superclass_of_atom = -1
      isuperclass = 0
      DO i = 1, n
         ! this atom belongs to a class with several equivalent atoms, and has not yet been found
         IF (superclass_of_atom(i) .EQ. -1 .AND. classes(class_of_atom(i))%Nele > 1) THEN
            isuperclass = isuperclass + 1
            CALL spread_superclass(i, isuperclass, superclass_of_atom, class_of_atom, classes, reference)
         END IF
      END DO
 
      ! put classes into superclasses
      ALLOCATE (superclasses(isuperclass))
      superclasses%Nele = 0
      DO i = 1, nclasses
         j = superclass_of_atom(classes(i)%reference(1))
         IF (j > 0) superclasses(j)%Nele = superclasses(j)%Nele + 1
      END DO
      DO i = 1, isuperclass
         ALLOCATE (superclasses(i)%classes(superclasses(i)%Nele))
         superclasses(i)%Nele = 0
      END DO
      DO i = 1, nclasses
         j = superclass_of_atom(classes(i)%reference(1))
         IF (j > 0) THEN
            superclasses(j)%Nele = superclasses(j)%Nele + 1
            superclasses(j)%classes(superclasses(j)%Nele) = i
         END IF
      END DO
 
      ! write(*,*) "Class generation time",t2-t1
      ! WRITE(*,*) "Nclasses, max size, total-non-1 ",Nclasses,MAXVAL(classes%Nele),COUNT(classes%Nele>1)
      ! write(*,*) "isuperclass ",isuperclass
 
      ! assign the order array to their initial value
      DO iclass = 1, nclasses
         order(classes(iclass)%unordered) = classes(iclass)%reference(classes(iclass)%order)
      END DO
 
      ! reorder the atoms superclass after superclass
      itries = 0
      DO i = 1, isuperclass
         DO
            itries = itries + 1
 
            ! assign the current order
            DO iclass = 1, superclasses(i)%Nele
               j = superclasses(i)%classes(iclass)
               order(classes(j)%unordered) = classes(j)%reference(classes(j)%order)
            END DO
 
            ! check for matches within this superclass only, be happy if we have a match
            matches = matrix_superclass_equal(reference, unordered, order, superclasses(i), classes)
            IF (itries > max_tries) THEN
               WRITE (*, *) "Could not find the 1-to-1 mapping to prove graph isomorphism"
               WRITE (*, *) "Reordering failed, assuming these molecules are different"
               EXIT
            END IF
            IF (matches) EXIT
 
            ! generate next permutation within this superclass
            DO iclass = 1, superclasses(i)%Nele
               j = superclasses(i)%classes(iclass)
               CALL all_permutations(classes(j)%order, classes(j)%Nele, &
                                     classes(j)%q, classes(j)%first)
               IF (.NOT. classes(j)%first) EXIT
            END DO
 
            ! we are back at the original permutation so we're unable to match this superclass.
            IF (iclass .EQ. superclasses(i)%Nele .AND. &
                classes(superclasses(i)%classes(superclasses(i)%Nele))%first) EXIT
         END DO
         ! failed in this superblock, can exit now
         IF (.NOT. matches) EXIT
      END DO
 
      ! the final check, just to be sure
      matches = matrix_equal(reference, unordered, order)
 
      DO iclass = 1, nclasses
         DEALLOCATE (classes(iclass)%reference)
         DEALLOCATE (classes(iclass)%unordered)
         DEALLOCATE (classes(iclass)%order)
         DEALLOCATE (classes(iclass)%q)
      END DO
      DEALLOCATE (classes)
      DO i = 1, isuperclass
         DEALLOCATE (superclasses(i)%classes)
      END DO
      DEALLOCATE (superclasses)
   END SUBROUTINE reorder_graph
 
! **************************************************************************************************
!> \brief spreads the superclass over all atoms of this class and all their bonded atoms
!>      provided that the latter belong to a class which contains more than one element
!> \param I ...
!> \param isuperclass ...
!> \param superclass_of_atom ...
!> \param class_of_atom ...
!> \param classes ...
!> \param reference ...
!> \par History
!>      09.2006 created [Joost VandeVondele]
! **************************************************************************************************
   RECURSIVE SUBROUTINE spread_superclass(I, isuperclass, superclass_of_atom, class_of_atom, &
                                          classes, reference)
      INTEGER, INTENT(IN)                                :: i, isuperclass
      INTEGER, DIMENSION(:), INTENT(INOUT)               :: superclass_of_atom
      INTEGER, DIMENSION(:), INTENT(IN)                  :: class_of_atom
      TYPE(class), DIMENSION(:), INTENT(IN)              :: classes
      TYPE(vertex), DIMENSION(:), INTENT(IN)             :: reference
 
      INTEGER                                            :: j
 
      IF (superclass_of_atom(i) .EQ. -1 .AND. classes(class_of_atom(i))%Nele > 1) THEN
         superclass_of_atom(i) = isuperclass
         DO j = 1, classes(class_of_atom(i))%Nele
            CALL spread_superclass(classes(class_of_atom(i))%reference(j), isuperclass, &
                                   superclass_of_atom, class_of_atom, classes, reference)
         END DO
         DO j = 1, SIZE(reference(i)%bonds)
            CALL spread_superclass(reference(i)%bonds(j), isuperclass, &
                                   superclass_of_atom, class_of_atom, classes, reference)
         END DO
      END IF
   END SUBROUTINE spread_superclass
 
! **************************************************************************************************
!> \brief determines of the vertices of this superclass have the same edges
!> \param reference ...
!> \param unordered ...
!> \param order ...
!> \param super ...
!> \param classes ...
!> \return ...
!> \par History
!>      09.2006 created [Joost VandeVondele]
! **************************************************************************************************
   FUNCTION matrix_superclass_equal(reference, unordered, order, super, classes) RESULT(res)
      TYPE(vertex), DIMENSION(:), INTENT(IN)             :: reference, unordered
      INTEGER, DIMENSION(:), INTENT(IN)                  :: order
      TYPE(superclass), INTENT(IN)                       :: super
      TYPE(class), DIMENSION(:), INTENT(IN)              :: classes
      LOGICAL                                            :: res
 
      INTEGER                                            :: i, iclass, iele, j
 
! I is the atom in the unordered set
 
      loop: DO iclass = 1, super%Nele
         DO iele = 1, classes(super%classes(iclass))%Nele
            i = classes(super%classes(iclass))%unordered(iele)
            res = (reference(order(i))%kind == unordered(i)%kind .AND. &
                   SIZE(reference(order(i))%bonds) == SIZE(unordered(i)%bonds))
            IF (res) THEN
               DO j = 1, SIZE(reference(order(i))%bonds)
                  IF (all(reference(order(i))%bonds(:) .NE. order(unordered(i)%bonds(j)))) THEN
                     res = .false.
                     EXIT loop
                  END IF
               END DO
            ELSE
               EXIT loop
            END IF
         END DO
      END DO loop
   END FUNCTION matrix_superclass_equal
 
! **************************************************************************************************
!> \brief determines of the vertices of the full set is equal, i.e.
!>      we have the same connectivity graph
!> \param reference ...
!> \param unordered ...
!> \param order ...
!> \return ...
!> \par History
!>      09.2006 created [Joost VandeVondele]
! **************************************************************************************************
   FUNCTION matrix_equal(reference, unordered, order) RESULT(res)
      TYPE(vertex), DIMENSION(:), INTENT(IN)             :: reference, unordered
      INTEGER, DIMENSION(:), INTENT(IN)                  :: order
      LOGICAL                                            :: res
 
      INTEGER                                            :: i, j
 
      loop: DO i = 1, SIZE(reference)
         res = (reference(order(i))%kind == unordered(i)%kind .AND. &
                SIZE(reference(order(i))%bonds) == SIZE(unordered(i)%bonds))
         IF (res) THEN
            DO j = 1, SIZE(reference(order(i))%bonds)
               IF (all(reference(order(i))%bonds(:) .NE. order(unordered(i)%bonds(j)))) THEN
                  res = .false.
                  EXIT loop
               END IF
            END DO
         ELSE
            EXIT loop
         END IF
      END DO loop
   END FUNCTION matrix_equal
 
! **************************************************************************************************
!> \brief creates a hash for an atom based on its own kind and on the kinds
!>       of its bonded neighbors
!> \param me ...
!> \param bonds ...
!> \return ...
!> \par History
!>      09.2006 created [Joost VandeVondele]
!> \note
!>       bonds are sorted so that the order of neighbors appearing in the bonded list
!>       is not important
! **************************************************************************************************
   FUNCTION hash_kind(me, bonds) RESULT(res)
      INTEGER, INTENT(IN)                                :: me
      INTEGER, DIMENSION(:), INTENT(IN)                  :: bonds
      INTEGER                                            :: res
 
      INTEGER                                            :: i, n
      INTEGER, ALLOCATABLE, DIMENSION(:)                 :: index, ordered_bonds
 
      n = SIZE(bonds)
      ALLOCATE (ordered_bonds(n + 1), index(n))
      DO i = 1, n
         ordered_bonds(i) = bonds(i)
      END DO
      ordered_bonds(n + 1) = me
      ! N: only sort the bonds, not me
      CALL sort(ordered_bonds, n, index)
      res = joaat_hash_i(ordered_bonds)
   END FUNCTION hash_kind
 
! **************************************************************************************************
!> \brief generates the hash of an array of integers and the index in the table
!> \param key an integer array of any length
!> \return ...
!> \par History
!>       09.2006 created [Joost VandeVondele]
!> \note
!>       http://en.wikipedia.org/wiki/Hash_table
!>       http://www.burtleburtle.net/bob/hash/doobs.html
!>       However, since fortran doesn't have an unsigned 4 byte int
!>       we compute it using an integer with the appropriate range
!>       we return already the index in the table as a final result
! **************************************************************************************************
   FUNCTION joaat_hash_i(key) RESULT(hash_index)
      INTEGER, DIMENSION(:), INTENT(IN)                  :: key
      INTEGER                                            :: hash_index
 
      INTEGER, PARAMETER                                 :: k64 = selected_int_kind(10)
      INTEGER(KIND=k64), PARAMETER                       :: b32 = 2_k64**32 - 1_k64
 
      INTEGER                                            :: i
      INTEGER(KIND=k64)                                  :: hash
 
      hash = 0_k64
      DO i = 1, SIZE(key)
         hash = iand(hash + ibits(key(i), 0, 8), b32)
         hash = iand(hash + iand(ishft(hash, 10), b32), b32)
         hash = iand(ieor(hash, iand(ishft(hash, -6), b32)), b32)
         hash = iand(hash + ibits(key(i), 8, 8), b32)
         hash = iand(hash + iand(ishft(hash, 10), b32), b32)
         hash = iand(ieor(hash, iand(ishft(hash, -6), b32)), b32)
         hash = iand(hash + ibits(key(i), 16, 8), b32)
         hash = iand(hash + iand(ishft(hash, 10), b32), b32)
         hash = iand(ieor(hash, iand(ishft(hash, -6), b32)), b32)
         hash = iand(hash + ibits(key(i), 24, 8), b32)
         hash = iand(hash + iand(ishft(hash, 10), b32), b32)
         hash = iand(ieor(hash, iand(ishft(hash, -6), b32)), b32)
      END DO
      hash = iand(hash + iand(ishft(hash, 3), b32), b32)
      hash = iand(ieor(hash, iand(ishft(hash, -11), b32)), b32)
      hash = iand(hash + iand(ishft(hash, 15), b32), b32)
      ! hash is the real 32bit hash value of the string,
      ! hash_index is an index in the hash_table
      hash_index = int(mod(hash, int(huge(hash_index), kind=k64)), kind=kind(hash_index))
   END FUNCTION joaat_hash_i
 
!===ACM Algorithm 323, Generation of Permutations in Lexicographic
!   Order (G6) by R. J. Ord-Smith, CACM 11 (Feb. 1968):117
!   Original Algorithm modified via Certification by I.M. Leitch,
!   17 March 1969.
! Algol to Fortran 77 by H.D.Knoble <hdkLESS at SPAM psu dot edu>,
!                                          May 1995.
!   x = initial values (/1...n/), first=.TRUE.
!   q = scratch
!   first = .TRUE. if you're back at the original order
! **************************************************************************************************
!> \brief ...
!> \param x ...
!> \param n ...
!> \param q ...
!> \param first ...
! **************************************************************************************************
   SUBROUTINE all_permutations(x, n, q, first)
      INTEGER                                            :: n, x(n), q(n)
      LOGICAL                                            :: first
 
      INTEGER                                            :: k, m, t
 
      IF (n == 1) RETURN
      IF (first) THEN
         first = .false.
         DO m = 1, n - 1
            q(m) = n
         END DO
      END IF
      IF (q(n - 1) .EQ. n) THEN
         q(n - 1) = n - 1
         t = x(n)
         x(n) = x(n - 1)
         x(n - 1) = t
         RETURN
      END IF
      DO k = n - 1, 1, -1
         IF (q(k) .EQ. k) THEN
            q(k) = n
         ELSE
            go to 1
         END IF
      END DO
      first = .true.
      k = 1
      GOTO 2
1     m = q(k)
      t = x(m)
      x(m) = x(k)
      x(k) = t
      q(k) = m - 1
      k = k + 1
2     m = n
      t = x(m)
      x(m) = x(k)
      x(k) = t
      m = m - 1
      k = k + 1
      DO WHILE (k .LT. m)
         t = x(m)
         x(m) = x(k)
         x(k) = t
         m = m - 1
         k = k + 1
      END DO
   END SUBROUTINE
END MODULE graphcon