topology_snapshot.F

Go to the documentation of this file.

!--------------------------------------------------------------------------------------------------!
!   CP2K: A general program to perform molecular dynamics simulations                              !
!   Copyright 2000-2026 CP2K developers group <https://cp2k.org>                                   !
!                                                                                                  !
!   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
!--------------------------------------------------------------------------------------------------!
 
! **************************************************************************************************
!> \brief Reader and physical moving-basis links for version-1 topology snapshots.
! **************************************************************************************************
MODULE topology_snapshot
   USE ai_moments,                      ONLY: cossin
   USE cp_files,                        ONLY: close_file,&
                                              open_file
   USE ieee_arithmetic,                 ONLY: ieee_is_finite
   USE kinds,                           ONLY: dp
   USE orbital_pointers,                ONLY: indco,&
                                              init_orbital_pointers
#include "./base/base_uses.f90"
 
   IMPLICIT NONE
   PRIVATE
   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'topology_snapshot'
   REAL(KIND=dp), PARAMETER :: pi = 3.1415926535897932384626433832795_dp, cell_tol = 1.e-10_dp
   TYPE :: snapshot_shell_type
      INTEGER :: first = 0, count = 0, lmin = 0, lmax = 0
      REAL(KIND=dp) :: center(3) = 0.0_dp, radius = 0.0_dp
      REAL(KIND=dp), ALLOCATABLE :: exponent(:), radii(:), contraction(:, :)
   END TYPE snapshot_shell_type
   TYPE :: snapshot_atom_type
      INTEGER :: kind = 0, count = 0
      REAL(KIND=dp) :: center(3) = 0.0_dp
      TYPE(snapshot_shell_type), ALLOCATABLE :: shells(:)
   END TYPE snapshot_atom_type
   TYPE :: snapshot_type
      INTEGER :: nao = 0, rank = 0, nspin = 0, channel = 0
      INTEGER :: periodic(3) = 0
      REAL(kind=dp) :: cell(3, 3) = 0.0_dp, inverse(3, 3) = 0.0_dp, k(3) = 0.0_dp
      TYPE(snapshot_atom_type), ALLOCATABLE :: atoms(:)
      INTEGER, ALLOCATABLE :: bands(:)
      REAL(kind=dp), ALLOCATABLE :: energies(:)
      COMPLEX(KIND=dp), ALLOCATABLE :: coefficients(:, :)
   END TYPE snapshot_type
   PUBLIC :: snapshot_type, read_snapshot, snapshot_overlap, check_snapshot, check_snapshot_seam
CONTAINS
 
! **************************************************************************************************
!> \brief Read a single point, retaining only its states; reject malformed/truncated exports.
!> \param filename Version-1 STATE_EXPORT file
!> \param point One-based point index
!> \param state Selected frame
! **************************************************************************************************
   SUBROUTINE read_snapshot(filename, point, state)
      CHARACTER(LEN=*), INTENT(IN)                       :: filename
      INTEGER, INTENT(IN)                                :: point
      TYPE(snapshot_type), INTENT(OUT)                   :: state
 
      CHARACTER(LEN=256)                                 :: header
      INTEGER                                            :: first, i, ia, ic, ios, ip, iset, ix, j, &
                                                            nat, nc, nk, np, ns, nt, offset, p, &
                                                            powers(3), refpowers(3), unit
      LOGICAL, ALLOCATABLE                               :: covered(:)
      REAL(kind=dp)                                      :: det, k(3), pair(2)
      REAL(kind=dp), ALLOCATABLE                         :: energies(:), row(:)
 
      CALL open_file(filename, unit_number=unit, file_status='OLD', file_action='READ')
      READ (unit, '(A)', iostat=ios) header
      cpassert(ios == 0)
      IF (trim(header) /= 'CP2K_TOPOLOGY_STATE 1') THEN
         cpabort('Unsupported topology snapshot version')
      END IF
      READ (unit, *, iostat=ios) nat, state%nao, state%rank, nk, state%nspin, nt, state%channel
      cpassert(ios == 0)
      IF (min(nat, state%nao, state%rank, nk, nt, state%channel) < 1) THEN
         cpabort('Invalid snapshot dimensions')
      END IF
      IF (state%rank > nt .OR. point < 1 .OR. point > nk) THEN
         cpabort('Invalid snapshot point or rank')
      END IF
      IF (state%nspin /= 1 .AND. state%nspin /= 2) THEN
         cpabort('Invalid snapshot spin components')
      END IF
      ALLOCATE (state%bands(state%rank), state%energies(nt), energies(nt), &
                state%coefficients(state%nao*state%nspin, state%rank), state%atoms(nat))
      READ (unit, *, iostat=ios) state%bands
      cpassert(ios == 0)
      IF (any(state%bands < 1) .OR. any(state%bands > nt)) THEN
         cpabort('Invalid snapshot band indices')
      END IF
      DO i = 2, state%rank
         IF (state%bands(i) <= state%bands(i - 1)) THEN
            cpabort('Snapshot bands must be strictly ordered')
         END IF
      END DO
      DO j = 1, 3
         READ (unit, *, iostat=ios) state%cell(:, j)
         cpassert(ios == 0)
      END DO
      IF (.NOT. all(ieee_is_finite(state%cell))) THEN
         cpabort('Nonfinite snapshot cell')
      END IF
      associate(a => state%cell, b => state%inverse)
         b(1, :) = [a(2, 2)*a(3, 3) - a(2, 3)*a(3, 2), a(1, 3)*a(3, 2) - a(1, 2)*a(3, 3), a(1, 2)*a(2, 3) - a(1, 3)*a(2, 2)]
         b(2, :) = [a(2, 3)*a(3, 1) - a(2, 1)*a(3, 3), a(1, 1)*a(3, 3) - a(1, 3)*a(3, 1), a(1, 3)*a(2, 1) - a(1, 1)*a(2, 3)]
         b(3, :) = [a(2, 1)*a(3, 2) - a(2, 2)*a(3, 1), a(1, 2)*a(3, 1) - a(1, 1)*a(3, 2), a(1, 1)*a(2, 2) - a(1, 2)*a(2, 1)]
         det = dot_product(a(:, 1), b(1, :))
         IF (abs(det) < cell_tol) THEN
            cpabort('Singular snapshot cell')
         END IF
         b(:, :) = b/det
      END associate
      READ (unit, *, iostat=ios) state%periodic
      cpassert(ios == 0)
      IF (any(state%periodic < 0) .OR. any(state%periodic > 1)) THEN
         cpabort('Invalid snapshot periodicity')
      END IF
      offset = 0
      DO ia = 1, nat
         associate(atom => state%atoms(ia))
            READ (unit, *, iostat=ios) ix, atom%kind, ns, atom%count
            cpassert(ios == 0)
            IF (ix /= ia .OR. min(atom%kind, ns, atom%count) < 1) THEN
               cpabort('Invalid snapshot atom')
            END IF
            READ (unit, *, iostat=ios) atom%center
            cpassert(ios == 0)
            IF (.NOT. all(ieee_is_finite(atom%center))) THEN
               cpabort('Nonfinite atom center')
            END IF
            ALLOCATE (atom%shells(ns), covered(atom%count))
            covered(:) = .false.
            DO iset = 1, ns
               associate(s => atom%shells(iset))
                  READ (unit, *, iostat=ios) first, s%count, np, nc, s%lmin, s%radius
                  cpassert(ios == 0)
                  IF (min(first, s%count, np, nc) < 1 .OR. s%lmin < 0) THEN
                     cpabort('Invalid Gaussian shell')
                  END IF
                  IF (.NOT. ieee_is_finite(s%radius) .OR. s%radius <= 0) THEN
                     cpabort('Invalid shell screening radius')
                  END IF
                  IF (first + s%count - 1 > atom%count) THEN
                     cpabort('Shell AO range exceeds atom')
                  END IF
                  IF (any(covered(first:first + s%count - 1))) THEN
                     cpabort('Overlapping shell AO ranges')
                  END IF
                  covered(first:first + s%count - 1) = .true.
                  s%lmax = 0
                  DO WHILE ((s%lmax + 1)*(s%lmax + 2)*(s%lmax + 3)/6 < nc)
                     s%lmax = s%lmax + 1
                  END DO
                  IF ((s%lmax + 1)*(s%lmax + 2)*(s%lmax + 3)/6 /= nc .OR. s%lmin > s%lmax) THEN
                     cpabort('Invalid Cartesian Gaussian count')
                  END IF
                  CALL init_orbital_pointers(s%lmax + 1)
                  s%first = offset + first
                  s%center(:) = atom%center
                  ALLOCATE (s%exponent(np), s%radii(np), s%contraction(nc*np, s%count), row(s%count))
                  DO ip = 1, np
                     READ (unit, *, iostat=ios) s%exponent(ip), s%radii(ip)
                     cpassert(ios == 0)
                     IF (.NOT. ieee_is_finite(s%exponent(ip)) .OR. .NOT. ieee_is_finite(s%radii(ip))) THEN
                        cpabort('Nonfinite primitive Gaussian')
                     END IF
                     IF (s%exponent(ip) <= 0 .OR. s%radii(ip) <= 0) THEN
                        cpabort('Invalid primitive Gaussian')
                     END IF
                     DO ic = 1, nc
                        READ (unit, *, iostat=ios) powers, row
                        cpassert(ios == 0)
                        refpowers(:) = indco(:, ic)
                        IF (any(powers /= refpowers)) THEN
                           cpabort('Unexpected Cartesian Gaussian ordering')
                        END IF
                        IF (.NOT. all(ieee_is_finite(row))) THEN
                           cpabort('Nonfinite contraction')
                        END IF
                        s%contraction((ip - 1)*nc + ic, :) = row
                     END DO
                  END DO
                  DEALLOCATE (row)
               END associate
            END DO
            IF (.NOT. all(covered)) THEN
               cpabort('Incomplete atom AO coverage')
            END IF
            DEALLOCATE (covered)
            offset = offset + atom%count
         END associate
      END DO
      IF (offset /= state%nao) THEN
         cpabort('Snapshot AO count mismatch')
      END IF
      DO p = 1, nk
         READ (unit, *, iostat=ios) ix, k
         cpassert(ios == 0)
         IF (ix /= p .OR. .NOT. all(ieee_is_finite(k))) THEN
            cpabort('Invalid snapshot k-point')
         END IF
         READ (unit, *, iostat=ios) energies
         cpassert(ios == 0)
         IF (.NOT. all(ieee_is_finite(energies))) THEN
            cpabort('Nonfinite snapshot spectrum')
         END IF
         IF (any(energies(2:) < energies(:nt - 1))) THEN
            cpabort('Unordered snapshot spectrum')
         END IF
         IF (p == point) THEN
            state%k(:) = k
            state%energies(:) = energies
         END IF
         DO j = 1, state%rank
            DO i = 1, state%nao*state%nspin
               READ (unit, *, iostat=ios) pair
               cpassert(ios == 0)
               IF (.NOT. all(ieee_is_finite(pair))) THEN
                  cpabort('Nonfinite state coefficient')
               END IF
               IF (p == point) state%coefficients(i, j) = cmplx(pair(1), pair(2), dp)
            END DO
         END DO
      END DO
      DO
         READ (unit, '(A)', iostat=ios) header
         IF (ios < 0) EXIT
         IF (ios /= 0 .OR. len_trim(header) /= 0) THEN
            cpabort('Unexpected trailing snapshot content')
         END IF
      END DO
      CALL close_file(unit)
   END SUBROUTINE read_snapshot
 
! **************************************************************************************************
!> \brief Apply screened Gaussian cross-geometry operator by atom blocks, without a dense AO matrix.
!> \param a Left physical state
!> \param b Right physical state
!> \param ka Unwrapped left fractional k-point
!> \param kb Unwrapped right fractional k-point
!> \param overlap Selected-state overlap including all spinor components
! **************************************************************************************************
   SUBROUTINE snapshot_overlap(a, b, ka, kb, overlap)
      TYPE(snapshot_type), INTENT(IN)                    :: a, b
      REAL(kind=dp), INTENT(IN)                          :: ka(3), kb(3)
      COMPLEX(KIND=dp), INTENT(OUT)                      :: overlap(:, :)
 
      COMPLEX(KIND=dp) :: part(size(overlap, 1), size(overlap, 2))
      INTEGER                                            :: ia
 
      IF (any(shape(overlap) /= [a%rank, b%rank])) THEN
         cpabort('Incorrect overlap output dimensions')
      END IF
      IF (a%nspin /= b%nspin .OR. a%channel /= b%channel .OR. a%rank /= b%rank) THEN
         cpabort('Changed spin or rank')
      END IF
      IF (any(a%bands /= b%bands)) THEN
         cpabort('Changed selected band indices')
      END IF
      IF (SIZE(a%atoms) /= SIZE(b%atoms) .OR. a%nao /= b%nao) THEN
         cpabort('Changed atom or AO count')
      END IF
      IF (maxval(abs(a%cell - b%cell)) > cell_tol .OR. any(a%periodic /= b%periodic)) THEN
         cpabort('Changed cell')
      END IF
      IF (.NOT. all(ieee_is_finite(ka)) .OR. .NOT. all(ieee_is_finite(kb))) THEN
         cpabort('Nonfinite k-point')
      END IF
      IF (maxval(abs(ka - a%k - anint(ka - a%k))) > cell_tol .OR. &
          maxval(abs(kb - b%k - anint(kb - b%k))) > cell_tol) THEN
         cpabort('Frame k-point does not match request')
      END IF
      overlap(:, :) = 0.0_dp
!$OMP PARALLEL DO DEFAULT(NONE) SHARED(A,b,ka,kb) PRIVATE(ia,part) REDUCTION(+:overlap) SCHEDULE(DYNAMIC)
      DO ia = 1, SIZE(a%atoms)
         CALL atom_overlap(a, b, ia, ka, kb, part)
         overlap(:, :) = overlap + part
      END DO
!$OMP END PARALLEL DO
   END SUBROUTINE snapshot_overlap
 
! **************************************************************************************************
!> \brief Screened atom-row contribution using CP2K's cossin primitive integral recurrence.
!> \param a Left physical state
!> \param b Right physical state
!> \param ia Left atom whose AO rows are contracted
!> \param ka Unwrapped left fractional k-point
!> \param kb Unwrapped right fractional k-point
!> \param atom_link This atom's contribution to the selected-state overlap
! **************************************************************************************************
   SUBROUTINE atom_overlap(a, b, ia, ka, kb, atom_link)
      TYPE(snapshot_type), INTENT(IN)                    :: a, b
      INTEGER, INTENT(IN)                                :: ia
      REAL(kind=dp), INTENT(IN)                          :: ka(3), kb(3)
      COMPLEX(KIND=dp), INTENT(OUT)                      :: atom_link(:, :)
 
      COMPLEX(KIND=dp)                                   :: phase
      COMPLEX(KIND=dp), ALLOCATABLE                      :: block(:, :), oc(:, :, :)
      INTEGER                                            :: af, bf, hi(3), i, ib, image(3), j, k, &
                                                            lo(3), nc_a, nc_b, s, sa, sb
      REAL(kind=dp)                                      :: bounds(3), cutoff, disp(3), q(3), rb(3)
      REAL(kind=dp), ALLOCATABLE                         :: cosine(:, :), sine(:, :)
 
      q(:) = 2.0_dp*pi*matmul(transpose(a%inverse), kb - ka)
      atom_link(:, :) = 0.0_dp
      DO sa = 1, SIZE(a%atoms(ia)%shells)
         associate(left => a%atoms(ia)%shells(sa))
            ALLOCATE (oc(left%count, b%rank, a%nspin))
            oc(:, :, :) = 0.0_dp
            nc_a = SIZE(left%contraction, 1)
            DO ib = 1, SIZE(b%atoms)
               DO sb = 1, SIZE(b%atoms(ib)%shells)
                  associate(right => b%atoms(ib)%shells(sb))
                     cutoff = left%radius + right%radius
                     disp(:) = matmul(a%inverse, right%center - left%center)
                     DO i = 1, 3
                        bounds(i) = cutoff*norm2(a%inverse(i, :))
                     END DO
                     lo(:) = ceiling(-disp - bounds)
                     hi(:) = floor(-disp + bounds)
                     WHERE (a%periodic == 0)
                        lo = 0
                        hi = 0
                     END WHERE
                     nc_b = SIZE(right%contraction, 1)
                     ALLOCATE (cosine(nc_a, nc_b), sine(nc_a, nc_b), block(left%count, right%count))
                     block(:, :) = 0.0_dp
                     DO k = lo(3), hi(3)
                        DO j = lo(2), hi(2)
                           DO i = lo(1), hi(1)
                              image(:) = [i, j, k]
                              rb(:) = right%center + matmul(a%cell, real(image, dp))
                              IF (norm2(rb - left%center) > cutoff) cycle
                              CALL cossin(left%lmax, SIZE(left%exponent), left%exponent, left%radii, left%lmin, &
                                          right%lmax, SIZE(right%exponent), right%exponent, right%radii, right%lmin, &
                                          left%center, rb, q, cosine, sine)
                              phase = exp(cmplx(0.0_dp, 2.0_dp*pi*dot_product(kb, real(image, dp)), dp))
                              block(:, :) = block + phase*matmul(transpose(left%contraction), &
                                                                 matmul(cmplx(cosine, -sine, dp), right%contraction))
                           END DO
                        END DO
                     END DO
                     DO s = 1, a%nspin
                        bf = right%first + (s - 1)*b%nao
                        oc(:, :, s) = oc(:, :, s) + matmul(block, b%coefficients(bf:bf + right%count - 1, :))
                     END DO
                     DEALLOCATE (cosine, sine, block)
                  END associate
               END DO
            END DO
            DO s = 1, a%nspin
               af = left%first + (s - 1)*a%nao
               atom_link(:, :) = atom_link + matmul(conjg(transpose(a%coefficients(af:af + left%count - 1, :))), oc(:, :, s))
            END DO
            DEALLOCATE (oc)
         END associate
      END DO
   END SUBROUTINE atom_overlap
 
! **************************************************************************************************
!> \brief Check AO-metric normalization and separation at every selected/excluded boundary.
!> \param a Physical frame to validate
!> \param metric_tol Maximum entrywise error in C^dagger S C minus identity
!> \param gap_tol Minimum sampled selected/excluded-band separation in hartree
!> \param metric_error Observed normalization error
!> \param gap Smallest sampled selected/excluded-band separation in hartree
! **************************************************************************************************
   SUBROUTINE check_snapshot(a, metric_tol, gap_tol, metric_error, gap)
      TYPE(snapshot_type), INTENT(IN)                    :: a
      REAL(kind=dp), INTENT(IN)                          :: metric_tol, gap_tol
      REAL(kind=dp), INTENT(OUT)                         :: metric_error, gap
 
      INTEGER                                            :: i
      LOGICAL                                            :: selected(size(a%energies))
      COMPLEX(KIND=dp)                                   :: metric(a%rank, a%rank)
 
      CALL snapshot_overlap(a, a, a%k, a%k, metric)
      DO i = 1, a%rank
         metric(i, i) = metric(i, i) - 1.0_dp
      END DO
      metric_error = maxval(abs(metric))
      IF (.NOT. ieee_is_finite(metric_error) .OR. metric_error > metric_tol) THEN
         cpabort('Frame is not AO-metric normalized')
      END IF
      selected(:) = .false.
      selected(a%bands) = .true.
      gap = huge(1.0_dp)
      DO i = 1, SIZE(a%energies) - 1
         IF (selected(i) .NEQV. selected(i + 1)) gap = min(gap, a%energies(i + 1) - a%energies(i))
      END DO
      IF (gap == huge(1.0_dp) .OR. gap <= gap_tol) THEN
         cpabort('No verified sampled subspace isolation')
      END IF
   END SUBROUTINE check_snapshot
 
! **************************************************************************************************
!> \brief Check explicitly prescribed atom permutation, periodic translations and basis at a seam.
!> \param a Start frame
!> \param b Endpoint frame
!> \param permutation One-based endpoint-atom to start-atom map
!> \param tolerance Absolute tolerance for geometry and basis agreement
! **************************************************************************************************
   SUBROUTINE check_snapshot_seam(a, b, permutation, tolerance)
      TYPE(snapshot_type), INTENT(IN)                    :: a, b
      INTEGER, INTENT(IN)                                :: permutation(:)
      REAL(kind=dp), INTENT(IN)                          :: tolerance
 
      INTEGER                                            :: i, j, n, s
      REAL(kind=dp)                                      :: delta(3), shift(3)
 
      n = SIZE(a%atoms)
      IF (SIZE(b%atoms) /= n .OR. SIZE(permutation) /= n) THEN
         cpabort('Seam atom count mismatch')
      END IF
      IF (any(permutation < 1) .OR. any(permutation > n)) THEN
         cpabort('Invalid seam permutation')
      END IF
      DO i = 1, n
         IF (count(permutation == i) /= 1) THEN
            cpabort('Seam atom map is not a permutation')
         END IF
      END DO
      IF (maxval(abs(a%cell - b%cell)) > tolerance .OR. any(a%periodic /= b%periodic)) THEN
         cpabort('Seam cell mismatch')
      END IF
      DO i = 1, n
         j = permutation(i)
         associate(left => a%atoms(j), right => b%atoms(i))
            IF (left%kind /= right%kind .OR. left%count /= right%count) THEN
               cpabort('Seam changes atom kind or AO count')
            END IF
            IF (SIZE(left%shells) /= SIZE(right%shells)) THEN
               cpabort('Seam changes basis')
            END IF
            delta(:) = right%center - left%center
            shift(:) = anint(matmul(a%inverse, delta))
            WHERE (a%periodic == 0) shift = 0.0_dp
            IF (maxval(abs(delta - matmul(a%cell, shift))) > tolerance) THEN
               cpabort('Geometry does not close at seam')
            END IF
            DO s = 1, SIZE(left%shells)
               associate(x => left%shells(s), y => right%shells(s))
                  IF (x%count /= y%count .OR. x%lmin /= y%lmin .OR. x%lmax /= y%lmax) THEN
                     cpabort('Seam changes shell')
                  END IF
                  IF (any(shape(x%contraction) /= shape(y%contraction))) THEN
                     cpabort('Seam changes contractions')
                  END IF
                  IF (SIZE(x%exponent) /= SIZE(y%exponent)) THEN
                     cpabort('Seam changes primitives')
                  END IF
                  IF (maxval(abs(x%radii - y%radii)) > tolerance .OR. abs(x%radius - y%radius) > tolerance) THEN
                     cpabort('Seam changes Gaussian screening radii')
                  END IF
                  IF (maxval(abs(x%exponent - y%exponent)) > tolerance .OR. &
                      maxval(abs(x%contraction - y%contraction)) > tolerance) THEN
                     cpabort('Seam changes physical basis')
                  END IF
               END associate
            END DO
         END associate
      END DO
   END SUBROUTINE check_snapshot_seam
END MODULE topology_snapshot