topology_phason_analysis.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 Native streaming mesh analysis with physical endpoint sewing.
! **************************************************************************************************
MODULE topology_phason_analysis
   USE cp_files,                        ONLY: close_file,&
                                              open_file
   USE ieee_arithmetic,                 ONLY: ieee_is_finite
   USE iso_fortran_env,                 ONLY: int64,&
                                              output_unit
   USE kinds,                           ONLY: default_path_length,&
                                              dp
   USE topology_curvature,              ONLY: curvature_density,&
                                              link_plaquette,&
                                              polar_link
   USE topology_snapshot,               ONLY: check_snapshot,&
                                              check_snapshot_seam,&
                                              read_snapshot,&
                                              snapshot_overlap,&
                                              snapshot_type
#include "./base/base_uses.f90"
 
   IMPLICIT NONE
   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'topology_phason_analysis'
   PRIVATE
   TYPE :: phason_result_type
      INTEGER :: dimension = 0, rank = 0, shape(4) = 0, nat = 0, point = 0
      REAL(kind=dp) :: value = 0.0_dp, origin(3) = 0.0_dp, directions(3, 4) = 0.0_dp
      CHARACTER(LEN=default_path_length) :: first_file = ''
   END TYPE phason_result_type
   PUBLIC :: analyze_mesh, phason_result_type, compare_refinements
 
CONTAINS
 
! **************************************************************************************************
!> \brief Flatten a zero-based index with the first coordinate varying fastest.
!> \param index Zero-based coordinate along each axis
!> \param shape Number of entries along each axis
!> \return One-based linear index
! **************************************************************************************************
   INTEGER FUNCTION flat(index, shape) RESULT(value)
      INTEGER, INTENT(IN)                                :: index(:), shape(:)
 
      INTEGER                                            :: d, stride
 
      value = 1
      stride = 1
      DO d = 1, SIZE(shape)
         value = value + index(d)*stride
         stride = stride*shape(d)
      END DO
   END FUNCTION flat
 
! **************************************************************************************************
!> \brief Inverse of flat; first coordinate varies fastest.
!> \param value One-based linear index
!> \param shape Number of entries along each axis
!> \param index Zero-based coordinate along each axis
! **************************************************************************************************
   SUBROUTINE unflat(value, shape, index)
      INTEGER, INTENT(IN)                                :: value, shape(:)
      INTEGER, INTENT(OUT)                               :: index(:)
 
      INTEGER                                            :: d, remainder
 
      remainder = value - 1
      DO d = 1, SIZE(shape)
         index(d) = mod(remainder, shape(d))
         remainder = remainder/shape(d)
      END DO
   END SUBROUTINE unflat
 
! **************************************************************************************************
!> \brief Build physical links or read reference links; stream plaquettes through an anonymous file.
!> \param path Manifest or reference-link filename
!> \param physical True for physical STATE_EXPORT frames, false for reference links
!> \param report Mesh metadata and unrounded Chern estimate for refinement comparison
! **************************************************************************************************
   SUBROUTINE analyze_mesh(path, physical, report)
      CHARACTER(LEN=*), INTENT(IN) :: path
      LOGICAL, INTENT(IN) :: physical
      TYPE(phason_result_type), INTENT(OUT) :: report
      CHARACTER(LEN=256) :: header
      CHARACTER(LEN=default_path_length), ALLOCATABLE :: files(:)
      INTEGER :: input, store, stat, dim, rank, nvertex, nframe, nat, vertex, axis, mu, nu, p, i, j, na, nb, nbase, idx, record_size
      INTEGER, ALLOCATABLE :: shape(:), extended(:), index(:), next(:), zero(:), points(:), permutation(:, :)
      LOGICAL, ALLOCATABLE :: checked(:)
      REAL(kind=dp), ALLOCATABLE :: offsets(:), directions(:, :)
      REAL(kind=dp) :: origin(3), ka(3), kb(3), tol_metric, tol_gap, tol_sv, phase_limit, minimum, local_min, &
                       max_metric, min_gap, err, gap, value, contribution, max_phase, local_phase, seam_error, &
                       valence_max, conduction_min, energy_shift, compensation, y, t
      COMPLEX(KIND=dp), ALLOCATABLE :: raw(:, :), link(:, :), sewing(:, :), u(:, :, :), pqs(:, :, :)
      TYPE(snapshot_type), ALLOCATABLE :: a, b, base
      LOGICAL :: common_reference, prefix
      INTEGER(KIND=int64) :: total
 
      CALL open_file(path, unit_number=input, file_status='OLD', file_action='READ')
      READ (input, '(A)', iostat=stat) header
      cpassert(stat == 0)
      IF (physical) THEN
         IF (trim(header) /= 'CP2K_PHASON_MESH 1') THEN
            cpabort('Expected CP2K_PHASON_MESH 1')
         END IF
      ELSE
         IF (trim(header) /= 'CP2K_TOPOLOGY_LINKS 1') THEN
            cpabort('Expected CP2K_TOPOLOGY_LINKS 1')
         END IF
      END IF
      READ (input, *, iostat=stat) dim, rank
      cpassert(stat == 0)
      IF ((dim /= 2 .AND. dim /= 4) .OR. rank < 1) THEN
         cpabort('Invalid mesh dimension or rank')
      END IF
      ALLOCATE (shape(dim), extended(dim), index(dim), next(dim), zero(dim))
      READ (input, *, iostat=stat) shape
      cpassert(stat == 0)
      IF (any(shape < 3)) THEN
         cpabort('At least three samples per mesh axis are required')
      END IF
      total = 1_int64
      DO axis = 1, dim
         IF (total > int(huge(1), int64)/int(shape(axis) + 1, int64)) THEN
            cpabort('Mesh index overflow')
         END IF
         total = total*int(shape(axis) + 1, int64)
      END DO
      extended(:) = shape + 1
      nvertex = product(shape)
      nframe = int(total)
      report%dimension = dim
      report%rank = rank
      report%shape(1:dim) = shape
      READ (input, *, iostat=stat) tol_metric, tol_gap, tol_sv, phase_limit
      cpassert(stat == 0)
      IF (.NOT. all(ieee_is_finite([tol_metric, tol_gap, tol_sv, phase_limit]))) THEN
         cpabort('Nonfinite tolerances')
      END IF
      IF (min(tol_metric, tol_gap, tol_sv, phase_limit) <= 0.0_dp .OR. max(tol_metric, tol_sv) >= 1.0_dp .OR. &
          phase_limit >= acos(-1.0_dp)) THEN
         cpabort('Invalid mesh tolerances')
      END IF
      ALLOCATE (raw(rank, rank), link(rank, rank), sewing(rank, rank), u(rank, rank, 4), pqs(rank, rank, dim*(dim - 1)/2))
      ! cp_files requires a filename; an anonymous scratch stream avoids persistent, stale link caches.
      OPEN (newunit=store, status='SCRATCH', access='STREAM', form='UNFORMATTED', iostat=stat)
      cpassert(stat == 0)
      INQUIRE (iolength=record_size) link
      minimum = huge(1.0_dp)
      max_metric = 0.0_dp
      min_gap = huge(1.0_dp)
      seam_error = 0.0_dp
      valence_max = -huge(1.0_dp)
      conduction_min = huge(1.0_dp)
      common_reference = .false.
      prefix = .true.
      IF (physical) THEN
         ALLOCATE (a, b, base)
         ALLOCATE (files(nframe), points(nframe), offsets(nframe), checked(nframe), directions(3, dim))
         checked(:) = .false.
         READ (input, *, iostat=stat) origin
         cpassert(stat == 0)
         IF (.NOT. all(ieee_is_finite(origin))) THEN
            cpabort('Nonfinite mesh origin')
         END IF
         DO axis = 1, dim
            READ (input, *, iostat=stat) directions(:, axis)
            cpassert(stat == 0)
         END DO
         IF (.NOT. all(ieee_is_finite(directions))) THEN
            cpabort('Nonfinite reciprocal periods')
         END IF
         IF (maxval(abs(directions - anint(directions))) > 1.e-10_dp) THEN
            cpabort('Noninteger reciprocal periods')
         END IF
         READ (input, *, iostat=stat) nat, common_reference
         cpassert(stat == 0)
         IF (nat < 1) THEN
            cpabort('Invalid mesh atom count')
         END IF
         ALLOCATE (permutation(nat, dim))
         DO axis = 1, dim
            READ (input, *, iostat=stat) permutation(:, axis)
            cpassert(stat == 0)
            IF (any(permutation(:, axis) < 1) .OR. any(permutation(:, axis) > nat)) THEN
               cpabort('Invalid seam permutation')
            END IF
            DO idx = 1, nat
               IF (count(permutation(:, axis) == idx) /= 1) THEN
                  cpabort('Non-bijective seam map')
               END IF
            END DO
         END DO
         DO mu = 1, dim
            DO nu = mu + 1, dim
               IF (any(permutation(permutation(:, mu), nu) /= permutation(permutation(:, nu), mu))) THEN
                  cpabort('The declared torus atom permutations do not commute')
               END IF
            END DO
         END DO
         DO idx = 1, nframe
            READ (input, *, iostat=stat) files(idx), points(idx), offsets(idx)
            cpassert(stat == 0)
            IF (points(idx) < 1 .OR. .NOT. ieee_is_finite(offsets(idx))) THEN
               cpabort('Invalid frame record')
            END IF
         END DO
         report%nat = nat
         report%origin(:) = origin
         report%directions(:, 1:dim) = directions
         report%first_file = files(1)
         report%point = points(1)
      END IF
      DO vertex = 1, nvertex
         CALL unflat(vertex, shape, index)
         IF (physical) THEN
            na = flat(index, extended)
            CALL frame(na, a)
            ka(:) = origin + matmul(directions, real(index, dp)/real(shape, dp))
         END IF
         DO axis = 1, dim
            IF (physical) THEN
               next(:) = index
               next(axis) = next(axis) + 1
               nb = flat(next, extended)
               CALL frame(nb, b)
               kb(:) = origin + matmul(directions, real(next, dp)/real(shape, dp))
               CALL snapshot_overlap(a, b, ka, kb, raw)
               IF (next(axis) == shape(axis)) THEN
                  zero(:) = next
                  zero(axis) = 0
                  nbase = flat(zero, extended)
                  CALL frame(nbase, base)
                  CALL check_snapshot_seam(base, b, permutation(:, axis), tol_metric)
                  CALL snapshot_overlap(b, base, kb, kb, sewing)
                  ! Both endpoint frames must span the same physical subspace, not merely have full rank.
                  link(:, :) = matmul(conjg(transpose(sewing)), sewing)
                  DO idx = 1, rank
                     link(idx, idx) = link(idx, idx) - 1.0_dp
                  END DO
                  err = maxval(abs(link))
                  seam_error = max(seam_error, err)
                  IF (err > tol_metric) THEN
                     cpabort('Endpoint selected subspaces do not close')
                  END IF
                  CALL polar_link(sewing, link, local_min, stat, tol_sv)
                  IF (stat /= 0) THEN
                     cpabort('Invalid endpoint sewing link')
                  END IF
                  raw(:, :) = matmul(raw, link)
               END IF
            ELSE
               DO j = 1, rank
                  DO i = 1, rank
                     READ (input, *, iostat=stat) ka(1:2)
                     cpassert(stat == 0)
                     raw(i, j) = cmplx(ka(1), ka(2), dp)
                  END DO
               END DO
            END IF
            CALL polar_link(raw, link, local_min, stat, tol_sv)
            IF (stat /= 0) THEN
               cpabort('Singular/nonfinite mesh link; refine or check the subspace')
            END IF
            minimum = min(minimum, local_min)
            WRITE (store, iostat=stat) link
            cpassert(stat == 0)
         END DO
      END DO
      DO
         READ (input, '(A)', iostat=stat) header
         IF (stat < 0) EXIT
         IF (stat /= 0 .OR. len_trim(header) /= 0) THEN
            cpabort('Unexpected trailing mesh data')
         END IF
      END DO
      CALL close_file(input)
      value = 0.0_dp
      compensation = 0.0_dp
      max_phase = 0.0_dp
      DO vertex = 1, nvertex
         CALL unflat(vertex, shape, index)
         p = 0
         DO mu = 1, dim
            DO nu = mu + 1, dim
               p = p + 1
               CALL read_link(vertex, mu, u(:, :, 1))
               next(:) = index
               next(mu) = mod(next(mu) + 1, shape(mu))
               CALL read_link(flat(next, shape), nu, u(:, :, 2))
               next(:) = index
               next(nu) = mod(next(nu) + 1, shape(nu))
               CALL read_link(flat(next, shape), mu, u(:, :, 3))
               CALL read_link(vertex, nu, u(:, :, 4))
               CALL link_plaquette(u(:, :, 1), u(:, :, 2), u(:, :, 3), u(:, :, 4), pqs(:, :, p))
            END DO
         END DO
         CALL curvature_density(pqs, dim, phase_limit, contribution, local_phase, stat)
         IF (stat /= 0) THEN
            cpabort('Unresolved/nonunitary plaquette; refine the mesh')
         END IF
         max_phase = max(max_phase, local_phase)
         y = contribution - compensation
         t = value + y
         compensation = (t - value) - y
         value = t
      END DO
      CLOSE (store)
      report%value = value
      WRITE (output_unit, '(A,I0)') 'PARAMETER_DIMENSION ', dim
      WRITE (output_unit, '(A,*(I0,1X))') 'MESH ', shape
      WRITE (output_unit, '(A,ES26.17)') 'CHERN_RAW ', value
      WRITE (output_unit, '(A,ES26.17)') 'MINIMUM_LINK_SINGULAR_VALUE ', minimum
      WRITE (output_unit, '(A,ES26.17)') 'MAXIMUM_PLAQUETTE_PHASE ', max_phase
      IF (physical) THEN
         WRITE (output_unit, '(A,ES26.17)') 'MAXIMUM_METRIC_ERROR ', max_metric
         WRITE (output_unit, '(A,ES26.17)') 'MAXIMUM_SEAM_ERROR ', seam_error
         WRITE (output_unit, '(A,ES26.17)') 'MINIMUM_DIRECT_GAP_HA ', min_gap
         IF (common_reference .AND. prefix) THEN
            WRITE (output_unit, '(A,ES26.17)') 'SAMPLED_INDIRECT_GAP_HA ', conduction_min - valence_max
         ELSE
            WRITE (output_unit, '(A)') 'INDIRECT_GAP_NOT_EVALUATED: no asserted common energy reference or not a lowest-band prefix'
         END IF
      END IF
      WRITE (output_unit, '(A)') 'MESH_CONVERGENCE_NOT_ESTABLISHED: repeat on finer meshes; retain unrounded C2 values.'
   CONTAINS
! **************************************************************************************************
!> \brief Read one frame and update validation diagnostics on its first use.
!> \param number One-based frame index in the extended mesh
!> \param state Physical frame with selected coefficients
! **************************************************************************************************
      SUBROUTINE frame(number, state)
      INTEGER, INTENT(IN)                                :: number
      TYPE(snapshot_type), INTENT(OUT)                   :: state
 
      INTEGER                                            :: iband
 
         CALL read_snapshot(trim(files(number)), points(number), state)
         IF (state%rank /= rank .OR. SIZE(state%atoms) /= nat) THEN
            cpabort('Frame rank/atom count differs from manifest')
         END IF
         IF (.NOT. checked(number)) THEN
            CALL check_snapshot(state, tol_metric, tol_gap, err, gap)
            max_metric = max(max_metric, err)
            min_gap = min(min_gap, gap)
            checked(number) = .true.
            DO iband = 1, rank
               IF (state%bands(iband) /= iband) prefix = .false.
            END DO
            IF (rank >= SIZE(state%energies)) prefix = .false.
            IF (prefix) THEN
               energy_shift = offsets(number)
               valence_max = max(valence_max, state%energies(rank) + energy_shift)
               conduction_min = min(conduction_min, state%energies(rank + 1) + energy_shift)
            END IF
         END IF
      END SUBROUTINE frame
! **************************************************************************************************
!> \brief Read a previously validated link from the scratch stream.
!> \param vertex_index One-based mesh vertex
!> \param direction One-based parameter axis
!> \param matrix Unitary selected-state link
! **************************************************************************************************
      SUBROUTINE read_link(vertex_index, direction, matrix)
      INTEGER, INTENT(IN)                                :: vertex_index, direction
      COMPLEX(KIND=dp), INTENT(OUT)                      :: matrix(:, :)
 
      INTEGER(KIND=int64)                                :: position
 
         position = (int(vertex_index - 1, int64)*int(dim, int64) + int(direction - 1, int64))*int(record_size, int64) + 1_int64
         READ (store, pos=position, iostat=stat) matrix
         cpassert(stat == 0)
      END SUBROUTINE read_link
   END SUBROUTINE analyze_mesh
 
! **************************************************************************************************
!> \brief Compare two jointly refined meshes of the same declared family; never round the estimate.
!> \param coarse Coarser mesh report
!> \param fine Finer mesh report
!> \param tolerance Maximum inter-mesh difference and distance to an integer
! **************************************************************************************************
   SUBROUTINE compare_refinements(coarse, fine, tolerance)
      TYPE(phason_result_type), INTENT(IN) :: coarse, fine
      REAL(kind=dp), INTENT(IN) :: tolerance
      TYPE(snapshot_type), ALLOCATABLE :: a, b
      COMPLEX(KIND=dp), ALLOCATABLE :: m(:, :), metric(:, :)
      INTEGER :: dim, i
      REAL(kind=dp) :: difference, residual
      IF (.NOT. ieee_is_finite(tolerance)) THEN
         cpabort('Nonfinite convergence tolerance')
      END IF
      IF (tolerance <= 0.0_dp .OR. tolerance >= 0.5_dp) THEN
         cpabort('Invalid convergence tolerance')
      END IF
      dim = coarse%dimension
      IF (fine%dimension /= dim .OR. fine%rank /= coarse%rank .OR. fine%nat /= coarse%nat) THEN
         cpabort('Incompatible refinement reports')
      END IF
      DO i = 1, dim
         IF (fine%shape(i) < 2*coarse%shape(i) .OR. mod(fine%shape(i), coarse%shape(i)) /= 0) THEN
            cpabort('Refine every mesh axis by an integer factor of at least two')
         END IF
      END DO
      IF (maxval(abs(coarse%origin - fine%origin)) > 1.e-10_dp .OR. &
          maxval(abs(coarse%directions - fine%directions)) > 1.e-10_dp) THEN
         cpabort('Refinement changes the parameter origin or reciprocal periods')
      END IF
      IF (coarse%nat > 0) THEN
         ALLOCATE (a, b)
         CALL read_snapshot(trim(coarse%first_file), coarse%point, a)
         CALL read_snapshot(trim(fine%first_file), fine%point, b)
         ALLOCATE (m(a%rank, a%rank), metric(a%rank, a%rank))
         CALL snapshot_overlap(a, b, a%k, a%k, m)
         metric(:, :) = matmul(conjg(transpose(m)), m)
         DO i = 1, a%rank
            metric(i, i) = metric(i, i) - 1.0_dp
         END DO
         IF (maxval(abs(metric)) > 1.e-7_dp) THEN
            cpabort('Refinement changes the initial physical subspace')
         END IF
      END IF
      difference = abs(fine%value - coarse%value)
      residual = abs(fine%value - anint(fine%value))
      WRITE (output_unit, '(A,ES26.17)') 'REFINEMENT_CHANGE ', difference
      WRITE (output_unit, '(A,ES26.17)') 'INTEGER_RESIDUAL ', residual
      IF (difference > tolerance .OR. residual > tolerance) THEN
         cpabort('Chern mesh refinement has not converged to the requested tolerance')
      END IF
      WRITE (output_unit, '(A)') 'REFINEMENT_CHECK_PASSED: numerical check, not a certificate of a bulk gap.'
   END SUBROUTINE compare_refinements
END MODULE topology_phason_analysis