topology_symmetry_unittest.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                                                      !
!--------------------------------------------------------------------------------------------------!
PROGRAM topology_symmetry_unittest
   USE ieee_arithmetic,                 ONLY: ieee_quiet_nan,&
                                              ieee_value
   USE iso_fortran_env,                 ONLY: output_unit
   USE kinds,                           ONLY: dp
   USE topology_symmetry,               ONLY: character_multiplicities,&
                                              inversion_representation
   USE topology_tqc,                    ONLY: inversion_ebr_signature,&
                                              inversion_indicators
 
   IMPLICIT NONE
   COMPLEX(KIND=dp) :: table(6, 3), group_chars(6), metric(2, 2), coeff(4, 4), phase(2), tmp(4)
   REAL(kind=dp) :: error, energy_error, energies(4)
   INTEGER :: multiplicities(3), status, counts(2), odd(8), ebr(2, 8), a, k, sign, strong, weak(3), z4, io_unit
   LOGICAL :: signed_atomic, nonnegative_atomic
 
   io_unit = output_unit
   ! Complete C3v/S3 character table per element, including a two-dimensional irrep.
   table(:, 1) = cmplx([1, 1, 1, 1, 1, 1], 0, dp)
   table(:, 2) = cmplx([1, 1, 1, -1, -1, -1], 0, dp)
   table(:, 3) = cmplx([2, -1, -1, 0, 0, 0], 0, dp)
   group_chars = table(:, 1) + 2*table(:, 3)
   CALL character_multiplicities(group_chars, table, multiplicities, 1.e-10_dp, status)
   IF (status /= 0 .OR. any(multiplicities /= [1, 0, 2])) error stop 'C3v decomposition failed'
   group_chars(2) = group_chars(2) + 0.1_dp
   CALL character_multiplicities(group_chars, table, multiplicities, 1.e-10_dp, status)
   IF (status == 0) error stop 'Invalid character accepted'
   metric = cmplx(0, 0, dp)
   metric(1, 1) = 2.0_dp
   metric(2, 2) = 0.5_dp
   coeff = cmplx(0, 0, dp)
   coeff(1, 1) = sqrt(0.5_dp)
   coeff(2, 2) = sqrt(2.0_dp)
   coeff(3, 3) = sqrt(0.5_dp)
   coeff(4, 4) = sqrt(2.0_dp)
   phase = cmplx([1, -1], 0, dp)
   energies = 0.0_dp
   CALL inversion_representation(metric, coeff, [1, 2], phase, energies, .true., &
                                 1.e-9_dp, 1.e-9_dp, counts, error, energy_error, status)
   IF (status /= 0 .OR. any(counts /= [2, 2])) error stop 'Metric/spinor parity failed'
   tmp = coeff(:, 1)
   coeff(:, 1) = (tmp + cmplx(0, 1, dp)*coeff(:, 2))/sqrt(2.0_dp)
   coeff(:, 2) = (cmplx(0, 1, dp)*tmp + coeff(:, 2))/sqrt(2.0_dp)
   CALL inversion_representation(metric, coeff, [1, 2], phase, energies, .true., &
                                 1.e-9_dp, 1.e-9_dp, counts, error, energy_error, status)
   IF (status /= 0 .OR. any(counts /= [2, 2])) error stop 'Degenerate-gauge parity failed'
   CALL inversion_representation(metric, coeff(:, :1), [1, 2], phase, energies(:1), .false., &
                                 1.e-9_dp, 1.e-9_dp, counts, error, energy_error, status)
   IF (status == 0) error stop 'Nonclosed parity subspace accepted'
   CALL inversion_representation(metric, coeff(:, :2), [1, 2], phase, energies(:2), .true., &
                                 1.e-9_dp, 1.e-9_dp, counts, error, energy_error, status)
   IF (status == 0) error stop 'Broken physical time reversal accepted'
   energies(2) = 1.0_dp
   CALL inversion_representation(metric, coeff, [1, 2], phase, energies, .true., &
                                 1.e-9_dp, 1.e-9_dp, counts, error, energy_error, status)
   IF (status == 0) error stop 'Noncommuting eigenspace accepted'
   ! Scalar two-site inversion with a genuinely off-diagonal AO metric.
   metric = cmplx(0.2_dp, 0, dp)
   metric(1, 1) = 1.0_dp
   metric(2, 2) = 1.0_dp
   coeff(:2, 1) = cmplx([1.0_dp, 1.0_dp]/sqrt(2.4_dp), 0, dp)
   coeff(:2, 2) = cmplx([1.0_dp, -1.0_dp]/sqrt(1.6_dp), 0, dp)
   phase = cmplx(1, 0, dp)
   CALL inversion_representation(metric, coeff(:2, :2), [2, 1], phase, energies(:2), .false., &
                                 1.e-9_dp, 1.e-9_dp, counts, error, energy_error, status)
   IF (status /= 0 .OR. any(counts /= [1, 1])) error stop 'Scalar two-site inversion failed'
   CALL inversion_representation(metric, coeff(:2, :1), [2, 1], -phase, energies(:1), .false., &
                                 1.e-9_dp, 1.e-9_dp, counts, error, energy_error, status)
   IF (status /= 0 .OR. any(counts /= [0, 1])) error stop 'Zone-boundary inversion phase failed'
   metric(1, 2) = cmplx(0.2_dp, 0.01_dp, dp)
   CALL inversion_representation(metric, coeff(:2, :2), [2, 1], phase, energies(:2), .false., &
                                 1.e-9_dp, 1.e-9_dp, counts, error, energy_error, status)
   IF (status == 0) error stop 'Invalid metric accepted'
   metric(1, 2) = cmplx(ieee_value(0.0_dp, ieee_quiet_nan), 0, dp)
   CALL inversion_representation(metric, coeff(:2, :2), [2, 1], phase, energies(:2), .false., &
                                 1.e-9_dp, 1.e-9_dp, counts, error, energy_error, status)
   IF (status == 0) error stop 'Nonfinite metric accepted'
   DO a = 0, 7
      DO sign = -1, 1, 2
         DO k = 0, 7
            odd(k + 1) = (1 - sign*(-1)**popcnt(iand(a, k)))/2
         END DO
         CALL inversion_indicators(odd, 1, 3, strong, weak, z4, status)
         IF (status /= 0 .OR. strong /= 0 .OR. any(weak /= 0) .OR. z4 /= 0) error stop 'Atomic indicator failed'
         CALL inversion_ebr_signature(odd, 1, signed_atomic, nonnegative_atomic, ebr, status)
         IF (status /= 0 .OR. .NOT. nonnegative_atomic .OR. .NOT. signed_atomic) error stop 'Atomic EBR failed'
         IF (sum(ebr) /= 1) error stop 'EBR rank failed'
      END DO
   END DO
   odd = 0
   odd(1) = 1
   CALL inversion_indicators(odd, 1, 3, strong, weak, z4, status)
   IF (status /= 0 .OR. strong /= 1 .OR. z4 /= 1) error stop 'Strong index failed'
   CALL inversion_ebr_signature(odd, 1, signed_atomic, nonnegative_atomic, ebr, status)
   IF (status /= 0 .OR. signed_atomic .OR. nonnegative_atomic) error stop 'Stable obstruction missed'
   odd(1) = 2
   CALL inversion_indicators(odd, 2, 3, strong, weak, z4, status)
   IF (status /= 0 .OR. strong /= 0 .OR. z4 /= 2) error stop 'Inversion Z4=2 missed'
   odd(1) = 4
   CALL inversion_ebr_signature(odd, 4, signed_atomic, nonnegative_atomic, ebr, status)
   IF (status /= 0 .OR. .NOT. signed_atomic .OR. nonnegative_atomic) error stop 'Signed-only EBR missed'
   odd = -1
   CALL inversion_ebr_signature(odd, 4, signed_atomic, nonnegative_atomic, ebr, status)
   IF (status == 0) error stop 'Invalid multiplicity accepted'
   CALL exhaustive_signatures(2)
   CALL exhaustive_signatures(3)
   WRITE (io_unit, *) 'Native characters, parity, inversion indicators and EBR tests passed.'
CONTAINS
 
! **************************************************************************************************
!> \brief Compare the exact solver against independent enumeration of rank-two atomic sums.
!> \param DIMENSION Dimension of the inversion subgroup
! **************************************************************************************************
   SUBROUTINE exhaustive_signatures(DIMENSION)
      INTEGER, INTENT(IN)                                :: dimension
 
      INTEGER :: code, digit, i, ierr, indicator, indices(3), j, n, nu, x, y, &
         values(2**DIMENSION), signatures(2**DIMENSION, 2**(DIMENSION + 1)), &
         reconstructed(2**DIMENSION), coefficients(2, 2**DIMENSION)
      LOGICAL :: indicators_zero, nonnegative_ok, possible(0:3**(2**DIMENSION) - 1), signed_ok
 
      n = 2**dimension
      DO x = 0, n - 1
         DO i = 0, n - 1
            signatures(i + 1, 2*x + 1) = (1 - (-1)**popcnt(iand(x, i)))/2
            signatures(i + 1, 2*x + 2) = 1 - signatures(i + 1, 2*x + 1)
         END DO
      END DO
      possible = .false.
      DO x = 1, 2*n
         DO y = 1, 2*n
            code = sum((signatures(:, x) + signatures(:, y))*[(3**i, i=0, n - 1)])
            possible(code) = .true.
         END DO
      END DO
      DO code = 0, 3**n - 1
         digit = code
         DO i = 1, n
            values(i) = mod(digit, 3)
            digit = digit/3
         END DO
         CALL inversion_ebr_signature(values, 2, signed_ok, nonnegative_ok, coefficients, ierr)
         IF (ierr /= 0 .OR. (nonnegative_ok .NEQV. possible(code))) error stop 'Enumerated EBR mismatch'
         CALL inversion_indicators(values, 2, dimension, nu, indices, indicator, ierr)
         IF (ierr /= 0) error stop 'Enumerated indicator failure'
         indicators_zero = nu == 0
         IF (dimension == 3) indicators_zero = indicator == 0 .AND. all(indices == 0)
         IF (signed_ok .NEQV. indicators_zero) error stop 'Indicator/atomic-lattice mismatch'
         IF (.NOT. nonnegative_ok) cycle
         reconstructed = 0
         DO j = 1, n
            reconstructed = reconstructed + coefficients(1, j)*signatures(:, 2*j - 1) + &
                            coefficients(2, j)*signatures(:, 2*j)
         END DO
         IF (any(reconstructed /= values) .OR. sum(coefficients) /= 2) error stop 'EBR reconstruction failed'
      END DO
   END SUBROUTINE exhaustive_signatures
END PROGRAM topology_symmetry_unittest