topology_curvature_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_curvature_unittest
   USE kinds,                           ONLY: dp
   USE topology_curvature,              ONLY: curvature_density,&
                                              link_plaquette,&
                                              polar_link
 
   IMPLICIT NONE
   REAL(kind=dp), PARAMETER :: pi = 3.1415926535897932384626433832795_dp, tolerance = 1.e-11_dp
   COMPLEX(KIND=dp) :: p(2, 2, 6), q(2, 2, 6), g(2, 2), link(2, 2)
   REAL(kind=dp) :: value, maximum, minimum, coarse, fine
   INTEGER :: i, status
   p(:, :, :) = 0.0_dp
   DO i = 1, 6
      p(1, 1, i) = 1.0_dp
      p(2, 2, i) = 1.0_dp
   END DO
   p(1, 1, 1) = exp(cmplx(0.0_dp, 0.2_dp, dp))
   p(2, 2, 1) = exp(cmplx(0.0_dp, -0.1_dp, dp))
   p(1, 1, 6) = exp(cmplx(0.0_dp, 0.3_dp, dp))
   p(2, 2, 6) = exp(cmplx(0.0_dp, 0.4_dp, dp))
   CALL curvature_density(p, 4, pi/2.0_dp, value, maximum, status)
   IF (status /= 0 .OR. abs(value - 0.02_dp/(4.0_dp*pi*pi)) > tolerance) error stop 'Incorrect C2 product trace'
   g(1, :) = [cmplx(1.0_dp, 0.0_dp, dp), cmplx(0.0_dp, 1.0_dp, dp)]/sqrt(2.0_dp)
   g(2, :) = [cmplx(0.0_dp, 1.0_dp, dp), cmplx(1.0_dp, 0.0_dp, dp)]/sqrt(2.0_dp)
   DO i = 1, 6
      q(:, :, i) = matmul(conjg(transpose(g)), matmul(p(:, :, i), g))
   END DO
   CALL curvature_density(q, 4, pi/2.0_dp, value, maximum, status)
   IF (status /= 0 .OR. abs(value - 0.02_dp/(4.0_dp*pi*pi)) > tolerance) error stop 'C2 gauge covariance failed'
   CALL curvature_density(p(:, :, 1:1), 2, pi/2.0_dp, value, maximum, status)
   IF (status /= 0 .OR. abs(value + 0.1_dp/(2.0_dp*pi)) > tolerance) error stop 'Incorrect determinant C1 phase'
   CALL curvature_density(p, 4, 0.1_dp, value, maximum, status)
   IF (status /= -3) error stop 'Unresolved C2 phase accepted'
   p(:, :, 1) = 0.0_dp
   CALL polar_link(p(:, :, 1), link, minimum, status, tolerance)
   IF (status == 0) error stop 'Singular link accepted'
   CALL curvature_density(p, 4, pi/2.0_dp, value, maximum, status)
   IF (status /= -2) error stop 'Nonunitary plaquette accepted'
   CALL dirac_pairing(6, coarse)
   CALL dirac_pairing(8, fine)
   IF (abs(fine + 1.0_dp) >= abs(coarse + 1.0_dp)) error stop 'C2 Dirac refinement failed'
   IF (abs(fine + 0.742412846261164_dp) > 1.e-6_dp) error stop 'Incorrect non-Abelian Dirac C2'
CONTAINS
 
! **************************************************************************************************
!> \brief Non-Abelian occupied doublet of the 4D Dirac model, exact C2=-1 at mass -3.
!> \param n ...
!> \param RESULT ...
! **************************************************************************************************
   SUBROUTINE dirac_pairing(n, RESULT)
      INTEGER, INTENT(IN)                                :: n
      REAL(KIND=dp), INTENT(OUT)                         :: result
 
      COMPLEX(KIND=dp)                                   :: gamma(4, 4, 5), h(4, 4), id(2, 2), &
                                                            overlap(2, 2), plaquettes(2, 2, 6), &
                                                            sx(2, 2), sy(2, 2), sz(2, 2), work(32)
      COMPLEX(KIND=dp), ALLOCATABLE                      :: links(:, :, :, :), states(:, :, :)
      INTEGER                                            :: at_mu, at_nu, d, INDEX(4), mu, nu, &
                                                            other, p, rem, status, v
      REAL(KIND=dp)                                      :: evals(4), k(4), local_value, maximum, &
                                                            minimum, rwork(12)
 
      sx = 0.0_dp
      sx(1, 2) = 1.0_dp
      sx(2, 1) = 1.0_dp
      sy = cmplx(0.0_dp, 1.0_dp, dp)*sx
      sy(1, 2) = -sy(1, 2)
      sz = 0.0_dp
      sz(1, 1) = 1.0_dp
      sz(2, 2) = -1.0_dp
      id = 0.0_dp
      id(1, 1) = 1.0_dp
      id(2, 2) = 1.0_dp
      CALL tensor(sx, sx, gamma(:, :, 1))
      CALL tensor(sx, sy, gamma(:, :, 2))
      CALL tensor(sx, sz, gamma(:, :, 3))
      CALL tensor(sy, id, gamma(:, :, 4))
      CALL tensor(sz, id, gamma(:, :, 5))
      ALLOCATE (states(4, 2, n**4), links(2, 2, 4, n**4))
      DO v = 1, n**4
         rem = v - 1
         DO d = 1, 4
            index(d) = mod(rem, n)
            rem = rem/n
         END DO
         k = 2.0_dp*pi*real(index, dp)/real(n, dp)
         h = cmplx(0.0_dp, 0.0_dp, dp)
         DO d = 1, 4
            h = h + sin(k(d))*gamma(:, :, d)
         END DO
         h = h + (-3.0_dp + sum(cos(k)))*gamma(:, :, 5)
         CALL zheev('V', 'U', 4, h, 4, evals, work, SIZE(work), rwork, status)
         IF (status /= 0) error stop 'Dirac eigensolver failed'
         states(:, :, v) = h(:, 1:2)
      END DO
      DO v = 1, n**4
         DO d = 1, 4
            other = neighbour(v, d, n)
            overlap = matmul(conjg(transpose(states(:, :, v))), states(:, :, other))
            CALL polar_link(overlap, links(:, :, d, v), minimum, status, tolerance)
            IF (status /= 0) error stop 'Singular Dirac link'
         END DO
      END DO
      result = 0.0_dp
      DO v = 1, n**4
         p = 0
         DO mu = 1, 4
            DO nu = mu + 1, 4
               p = p + 1
               at_mu = neighbour(v, mu, n)
               at_nu = neighbour(v, nu, n)
               CALL link_plaquette(links(:, :, mu, v), links(:, :, nu, at_mu), links(:, :, mu, at_nu), &
                                   links(:, :, nu, v), plaquettes(:, :, p))
            END DO
         END DO
         CALL curvature_density(plaquettes, 4, pi/2.0_dp, local_value, maximum, status)
         IF (status /= 0) error stop 'Invalid Dirac plaquette'
         result = result + local_value
      END DO
   END SUBROUTINE dirac_pairing
 
! **************************************************************************************************
!> \brief Periodic forward neighbour in a first-axis-fastest mesh.
!> \param v ...
!> \param d ...
!> \param n ...
!> \return ...
! **************************************************************************************************
   INTEGER FUNCTION neighbour(v, d, n) RESULT(other)
      INTEGER, INTENT(IN)                                :: v, d, n
 
      INTEGER                                            :: coordinate, stride
 
      stride = n**(d - 1)
      coordinate = mod((v - 1)/stride, n)
      other = v + stride
      IF (coordinate == n - 1) other = other - n*stride
   END FUNCTION neighbour
 
! **************************************************************************************************
!> \brief Kronecker product of two 2x2 matrices.
!> \param a ...
!> \param b ...
!> \param c ...
! **************************************************************************************************
   SUBROUTINE tensor(a, b, c)
      COMPLEX(KIND=dp), INTENT(IN)                       :: a(2, 2), b(2, 2)
      COMPLEX(KIND=dp), INTENT(OUT)                      :: c(4, 4)
 
      INTEGER                                            :: i, j
 
      DO j = 1, 2
         DO i = 1, 2
            c(2*i - 1:2*i, 2*j - 1:2*j) = a(i, j)*b
         END DO
      END DO
   END SUBROUTINE tensor
END PROGRAM topology_curvature_unittest