Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 2 : PROGRAM topology_curvature_unittest
9 2 : USE kinds, ONLY: dp
10 : USE topology_curvature, ONLY: curvature_density,&
11 : link_plaquette,&
12 : polar_link
13 :
14 : IMPLICIT NONE
15 : REAL(KIND=dp), PARAMETER :: pi = 3.1415926535897932384626433832795_dp, tolerance = 1.e-11_dp
16 : COMPLEX(KIND=dp) :: p(2, 2, 6), q(2, 2, 6), g(2, 2), link(2, 2)
17 : REAL(KIND=dp) :: value, maximum, minimum, dirac_value
18 : INTEGER :: i, status
19 2 : p(:, :, :) = 0.0_dp
20 14 : DO i = 1, 6
21 12 : p(1, 1, i) = 1.0_dp
22 14 : p(2, 2, i) = 1.0_dp
23 : END DO
24 2 : p(1, 1, 1) = EXP(CMPLX(0.0_dp, 0.2_dp, dp))
25 2 : p(2, 2, 1) = EXP(CMPLX(0.0_dp, -0.1_dp, dp))
26 2 : p(1, 1, 6) = EXP(CMPLX(0.0_dp, 0.3_dp, dp))
27 2 : p(2, 2, 6) = EXP(CMPLX(0.0_dp, 0.4_dp, dp))
28 2 : CALL curvature_density(p, 4, pi/2.0_dp, value, maximum, status)
29 2 : IF (status /= 0 .OR. ABS(value - 0.02_dp/(4.0_dp*pi*pi)) > tolerance) ERROR STOP 'Incorrect C2 product trace'
30 6 : g(1, :) = [CMPLX(1.0_dp, 0.0_dp, dp), CMPLX(0.0_dp, 1.0_dp, dp)]/SQRT(2.0_dp)
31 6 : g(2, :) = [CMPLX(0.0_dp, 1.0_dp, dp), CMPLX(1.0_dp, 0.0_dp, dp)]/SQRT(2.0_dp)
32 14 : DO i = 1, 6
33 434 : q(:, :, i) = MATMUL(CONJG(TRANSPOSE(g)), MATMUL(p(:, :, i), g))
34 : END DO
35 2 : CALL curvature_density(q, 4, pi/2.0_dp, value, maximum, status)
36 2 : IF (status /= 0 .OR. ABS(value - 0.02_dp/(4.0_dp*pi*pi)) > tolerance) ERROR STOP 'C2 gauge covariance failed'
37 2 : CALL curvature_density(p(:, :, 1:1), 2, pi/2.0_dp, value, maximum, status)
38 2 : IF (status /= 0 .OR. ABS(value + 0.1_dp/(2.0_dp*pi)) > tolerance) ERROR STOP 'Incorrect determinant C1 phase'
39 2 : CALL curvature_density(p, 4, 0.1_dp, value, maximum, status)
40 2 : IF (status /= -3) ERROR STOP 'Unresolved C2 phase accepted'
41 14 : p(:, :, 1) = 0.0_dp
42 2 : CALL polar_link(p(:, :, 1), link, minimum, status, tolerance)
43 2 : IF (status == 0) ERROR STOP 'Singular link accepted'
44 2 : CALL curvature_density(p, 4, pi/2.0_dp, value, maximum, status)
45 2 : IF (status /= -2) ERROR STOP 'Nonunitary plaquette accepted'
46 2 : CALL dirac_pairing(8, dirac_value)
47 2 : IF (ABS(dirac_value + 0.742412846261164_dp) > 1.e-6_dp) ERROR STOP 'Incorrect non-Abelian Dirac C2'
48 : CONTAINS
49 :
50 : ! **************************************************************************************************
51 : !> \brief Non-Abelian occupied doublet of the 4D Dirac model, exact C2=-1 at mass -3.
52 : !> \param n ...
53 : !> \param RESULT ...
54 : ! **************************************************************************************************
55 2 : SUBROUTINE dirac_pairing(n, RESULT)
56 : INTEGER, INTENT(IN) :: n
57 : REAL(KIND=dp), INTENT(OUT) :: result
58 :
59 : COMPLEX(KIND=dp) :: gamma(4, 4, 5), h(4, 4), id(2, 2), &
60 : overlap(2, 2), plaquettes(2, 2, 6), &
61 : sx(2, 2), sy(2, 2), sz(2, 2), work(32)
62 2 : COMPLEX(KIND=dp), ALLOCATABLE :: links(:, :, :, :), states(:, :, :)
63 : INTEGER :: at_mu, at_nu, d, INDEX(4), mu, nu, &
64 : other, p, rem, status, v
65 : REAL(KIND=dp) :: evals(4), k(4), local_value, maximum, &
66 : minimum, rwork(12)
67 :
68 2 : sx = 0.0_dp
69 2 : sx(1, 2) = 1.0_dp
70 2 : sx(2, 1) = 1.0_dp
71 14 : sy = CMPLX(0.0_dp, 1.0_dp, dp)*sx
72 2 : sy(1, 2) = -sy(1, 2)
73 2 : sz = 0.0_dp
74 2 : sz(1, 1) = 1.0_dp
75 2 : sz(2, 2) = -1.0_dp
76 2 : id = 0.0_dp
77 2 : id(1, 1) = 1.0_dp
78 2 : id(2, 2) = 1.0_dp
79 2 : CALL tensor(sx, sx, gamma(:, :, 1))
80 2 : CALL tensor(sx, sy, gamma(:, :, 2))
81 2 : CALL tensor(sx, sz, gamma(:, :, 3))
82 2 : CALL tensor(sy, id, gamma(:, :, 4))
83 2 : CALL tensor(sz, id, gamma(:, :, 5))
84 10 : ALLOCATE (states(4, 2, n**4), links(2, 2, 4, n**4))
85 8194 : DO v = 1, n**4
86 8192 : rem = v - 1
87 40960 : DO d = 1, 4
88 32768 : INDEX(d) = MOD(rem, n)
89 40960 : rem = rem/n
90 : END DO
91 40960 : k = 2.0_dp*pi*REAL(index, dp)/REAL(n, dp)
92 8192 : h = CMPLX(0.0_dp, 0.0_dp, dp)
93 40960 : DO d = 1, 4
94 696320 : h = h + SIN(k(d))*gamma(:, :, d)
95 : END DO
96 204800 : h = h + (-3.0_dp + SUM(COS(k)))*gamma(:, :, 5)
97 8192 : CALL zheev('V', 'U', 4, h, 4, evals, work, SIZE(work), rwork, status)
98 8192 : IF (status /= 0) ERROR STOP 'Dirac eigensolver failed'
99 90114 : states(:, :, v) = h(:, 1:2)
100 : END DO
101 8194 : DO v = 1, n**4
102 40962 : DO d = 1, 4
103 32768 : other = neighbour(v, d, n)
104 753664 : overlap = MATMUL(CONJG(TRANSPOSE(states(:, :, v))), states(:, :, other))
105 32768 : CALL polar_link(overlap, links(:, :, d, v), minimum, status, tolerance)
106 40960 : IF (status /= 0) ERROR STOP 'Singular Dirac link'
107 : END DO
108 : END DO
109 2 : RESULT = 0.0_dp
110 8194 : DO v = 1, n**4
111 8192 : p = 0
112 40960 : DO mu = 1, 4
113 90112 : DO nu = mu + 1, 4
114 49152 : p = p + 1
115 49152 : at_mu = neighbour(v, mu, n)
116 49152 : at_nu = neighbour(v, nu, n)
117 : CALL link_plaquette(links(:, :, mu, v), links(:, :, nu, at_mu), links(:, :, mu, at_nu), &
118 81920 : links(:, :, nu, v), plaquettes(:, :, p))
119 : END DO
120 : END DO
121 8192 : CALL curvature_density(plaquettes, 4, pi/2.0_dp, local_value, maximum, status)
122 8192 : IF (status /= 0) ERROR STOP 'Invalid Dirac plaquette'
123 8194 : RESULT = RESULT + local_value
124 : END DO
125 4 : END SUBROUTINE dirac_pairing
126 :
127 : ! **************************************************************************************************
128 : !> \brief Periodic forward neighbour in a first-axis-fastest mesh.
129 : !> \param v ...
130 : !> \param d ...
131 : !> \param n ...
132 : !> \return ...
133 : ! **************************************************************************************************
134 131072 : INTEGER FUNCTION neighbour(v, d, n) RESULT(other)
135 : INTEGER, INTENT(IN) :: v, d, n
136 :
137 : INTEGER :: coordinate, stride
138 :
139 131072 : stride = n**(d - 1)
140 131072 : coordinate = MOD((v - 1)/stride, n)
141 131072 : other = v + stride
142 131072 : IF (coordinate == n - 1) other = other - n*stride
143 131072 : END FUNCTION neighbour
144 :
145 : ! **************************************************************************************************
146 : !> \brief Kronecker product of two 2x2 matrices.
147 : !> \param a ...
148 : !> \param b ...
149 : !> \param c ...
150 : ! **************************************************************************************************
151 10 : SUBROUTINE tensor(a, b, c)
152 : COMPLEX(KIND=dp), INTENT(IN) :: a(2, 2), b(2, 2)
153 : COMPLEX(KIND=dp), INTENT(OUT) :: c(4, 4)
154 :
155 : INTEGER :: i, j
156 :
157 30 : DO j = 1, 2
158 70 : DO i = 1, 2
159 300 : c(2*i - 1:2*i, 2*j - 1:2*j) = a(i, j)*b
160 : END DO
161 : END DO
162 10 : END SUBROUTINE tensor
163 : END PROGRAM topology_curvature_unittest
|
