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, coarse, fine
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(6, coarse)
47 2 : CALL dirac_pairing(8, fine)
48 2 : IF (ABS(fine + 1.0_dp) >= ABS(coarse + 1.0_dp)) ERROR STOP 'C2 Dirac refinement failed'
49 2 : IF (ABS(fine + 0.742412846261164_dp) > 1.e-6_dp) ERROR STOP 'Incorrect non-Abelian Dirac C2'
50 : CONTAINS
51 :
52 : ! **************************************************************************************************
53 : !> \brief Non-Abelian occupied doublet of the 4D Dirac model, exact C2=-1 at mass -3.
54 : !> \param n ...
55 : !> \param RESULT ...
56 : ! **************************************************************************************************
57 4 : SUBROUTINE dirac_pairing(n, RESULT)
58 : INTEGER, INTENT(IN) :: n
59 : REAL(KIND=dp), INTENT(OUT) :: result
60 :
61 : COMPLEX(KIND=dp) :: gamma(4, 4, 5), h(4, 4), id(2, 2), &
62 : overlap(2, 2), plaquettes(2, 2, 6), &
63 : sx(2, 2), sy(2, 2), sz(2, 2), work(32)
64 4 : COMPLEX(KIND=dp), ALLOCATABLE :: links(:, :, :, :), states(:, :, :)
65 : INTEGER :: at_mu, at_nu, d, INDEX(4), mu, nu, &
66 : other, p, rem, status, v
67 : REAL(KIND=dp) :: evals(4), k(4), local_value, maximum, &
68 : minimum, rwork(12)
69 :
70 4 : sx = 0.0_dp
71 4 : sx(1, 2) = 1.0_dp
72 4 : sx(2, 1) = 1.0_dp
73 28 : sy = CMPLX(0.0_dp, 1.0_dp, dp)*sx
74 4 : sy(1, 2) = -sy(1, 2)
75 4 : sz = 0.0_dp
76 4 : sz(1, 1) = 1.0_dp
77 4 : sz(2, 2) = -1.0_dp
78 4 : id = 0.0_dp
79 4 : id(1, 1) = 1.0_dp
80 4 : id(2, 2) = 1.0_dp
81 4 : CALL tensor(sx, sx, gamma(:, :, 1))
82 4 : CALL tensor(sx, sy, gamma(:, :, 2))
83 4 : CALL tensor(sx, sz, gamma(:, :, 3))
84 4 : CALL tensor(sy, id, gamma(:, :, 4))
85 4 : CALL tensor(sz, id, gamma(:, :, 5))
86 20 : ALLOCATE (states(4, 2, n**4), links(2, 2, 4, n**4))
87 10788 : DO v = 1, n**4
88 10784 : rem = v - 1
89 53920 : DO d = 1, 4
90 43136 : INDEX(d) = MOD(rem, n)
91 53920 : rem = rem/n
92 : END DO
93 53920 : k = 2.0_dp*pi*REAL(index, dp)/REAL(n, dp)
94 10784 : h = CMPLX(0.0_dp, 0.0_dp, dp)
95 53920 : DO d = 1, 4
96 916640 : h = h + SIN(k(d))*gamma(:, :, d)
97 : END DO
98 269600 : h = h + (-3.0_dp + SUM(COS(k)))*gamma(:, :, 5)
99 10784 : CALL zheev('V', 'U', 4, h, 4, evals, work, SIZE(work), rwork, status)
100 10784 : IF (status /= 0) ERROR STOP 'Dirac eigensolver failed'
101 118628 : states(:, :, v) = h(:, 1:2)
102 : END DO
103 10788 : DO v = 1, n**4
104 53924 : DO d = 1, 4
105 43136 : other = neighbour(v, d, n)
106 992128 : overlap = MATMUL(CONJG(TRANSPOSE(states(:, :, v))), states(:, :, other))
107 43136 : CALL polar_link(overlap, links(:, :, d, v), minimum, status, tolerance)
108 53920 : IF (status /= 0) ERROR STOP 'Singular Dirac link'
109 : END DO
110 : END DO
111 4 : RESULT = 0.0_dp
112 10788 : DO v = 1, n**4
113 10784 : p = 0
114 53920 : DO mu = 1, 4
115 118624 : DO nu = mu + 1, 4
116 64704 : p = p + 1
117 64704 : at_mu = neighbour(v, mu, n)
118 64704 : at_nu = neighbour(v, nu, n)
119 : CALL link_plaquette(links(:, :, mu, v), links(:, :, nu, at_mu), links(:, :, mu, at_nu), &
120 107840 : links(:, :, nu, v), plaquettes(:, :, p))
121 : END DO
122 : END DO
123 10784 : CALL curvature_density(plaquettes, 4, pi/2.0_dp, local_value, maximum, status)
124 10784 : IF (status /= 0) ERROR STOP 'Invalid Dirac plaquette'
125 10788 : RESULT = RESULT + local_value
126 : END DO
127 6 : END SUBROUTINE dirac_pairing
128 :
129 : ! **************************************************************************************************
130 : !> \brief Periodic forward neighbour in a first-axis-fastest mesh.
131 : !> \param v ...
132 : !> \param d ...
133 : !> \param n ...
134 : !> \return ...
135 : ! **************************************************************************************************
136 172544 : INTEGER FUNCTION neighbour(v, d, n) RESULT(other)
137 : INTEGER, INTENT(IN) :: v, d, n
138 :
139 : INTEGER :: coordinate, stride
140 :
141 172544 : stride = n**(d - 1)
142 172544 : coordinate = MOD((v - 1)/stride, n)
143 172544 : other = v + stride
144 172544 : IF (coordinate == n - 1) other = other - n*stride
145 172544 : END FUNCTION neighbour
146 :
147 : ! **************************************************************************************************
148 : !> \brief Kronecker product of two 2x2 matrices.
149 : !> \param a ...
150 : !> \param b ...
151 : !> \param c ...
152 : ! **************************************************************************************************
153 20 : SUBROUTINE tensor(a, b, c)
154 : COMPLEX(KIND=dp), INTENT(IN) :: a(2, 2), b(2, 2)
155 : COMPLEX(KIND=dp), INTENT(OUT) :: c(4, 4)
156 :
157 : INTEGER :: i, j
158 :
159 60 : DO j = 1, 2
160 140 : DO i = 1, 2
161 600 : c(2*i - 1:2*i, 2*j - 1:2*j) = a(i, j)*b
162 : END DO
163 : END DO
164 20 : END SUBROUTINE tensor
165 : END PROGRAM topology_curvature_unittest
|
