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 : ! **************************************************************************************************
9 : !> \brief Gauge-covariant C1/C2 plaquette kernels, independent of the state representation.
10 : !> \note C2 is the second Chern character pairing, not generally the second Chern class.
11 : ! **************************************************************************************************
12 : MODULE topology_curvature
13 : USE ieee_arithmetic, ONLY: ieee_is_finite
14 : USE kinds, ONLY: dp
15 : USE topology_wilson, ONLY: wilson_step
16 :
17 : IMPLICIT NONE
18 : PRIVATE
19 : REAL(KIND=dp), PARAMETER :: pi = 3.1415926535897932384626433832795_dp, unitary_tol = 1.e-8_dp
20 : PUBLIC :: polar_link, link_plaquette, curvature_density
21 : CONTAINS
22 :
23 : ! **************************************************************************************************
24 : !> \brief Polar-unitarize a physical link, retaining its minimum singular value.
25 : !> \param overlap Raw selected-state overlap matrix
26 : !> \param link Unitary polar factor
27 : !> \param minimum Smallest singular value of the raw overlap
28 : !> \param status Zero on success; nonzero for an invalid or singular link
29 : !> \param tolerance Minimum allowed singular value
30 : ! **************************************************************************************************
31 43138 : SUBROUTINE polar_link(overlap, link, minimum, status, tolerance)
32 : COMPLEX(KIND=dp), INTENT(IN) :: overlap(:, :)
33 : COMPLEX(KIND=dp), INTENT(OUT) :: link(:, :)
34 : REAL(KIND=dp), INTENT(OUT) :: minimum
35 : INTEGER, INTENT(OUT) :: status
36 : REAL(KIND=dp), INTENT(IN) :: tolerance
37 :
38 : INTEGER :: i
39 :
40 301966 : link(:, :) = 0.0_dp
41 129414 : DO i = 1, MIN(SIZE(link, 1), SIZE(link, 2))
42 129414 : link(i, i) = 1.0_dp
43 : END DO
44 43138 : CALL wilson_step(link, overlap, minimum, status, tolerance)
45 43138 : END SUBROUTINE polar_link
46 :
47 : ! **************************************************************************************************
48 : !> \brief Oriented plaquette U_mu(x) U_nu(x+mu) U_mu(x+nu)^dagger U_nu(x)^dagger.
49 : !> \param u_mu Link in direction mu based at x
50 : !> \param u_nu_at_mu Link in direction nu based at x+mu
51 : !> \param u_mu_at_nu Link in direction mu based at x+nu
52 : !> \param u_nu Link in direction nu based at x
53 : !> \param p Oriented unitary plaquette based at x
54 : ! **************************************************************************************************
55 64704 : SUBROUTINE link_plaquette(u_mu, u_nu_at_mu, u_mu_at_nu, u_nu, p)
56 : COMPLEX(KIND=dp), INTENT(IN) :: u_mu(:, :), u_nu_at_mu(:, :), &
57 : u_mu_at_nu(:, :), u_nu(:, :)
58 : COMPLEX(KIND=dp), INTENT(OUT) :: p(:, :)
59 :
60 3882240 : p(:, :) = MATMUL(MATMUL(u_mu, u_nu_at_mu), MATMUL(CONJG(TRANSPOSE(u_mu_at_nu)), CONJG(TRANSPOSE(u_nu))))
61 64704 : END SUBROUTINE link_plaquette
62 :
63 : ! **************************************************************************************************
64 : !> \brief Local contribution to C1 or C2. Report unrounded values and reject unresolved phases.
65 : !> \param plaquettes Order 12 for C1; 12,13,14,23,24,34 for C2; all based at the same vertex
66 : !> \param ndim Parameter dimension, 2 or 4
67 : !> \param phase_limit Maximum allowed principal plaquette phase, strictly between zero and pi
68 : !> \param value Unrounded local contribution
69 : !> \param maximum Maximum absolute plaquette phase
70 : !> \param status Zero on success, negative invalid/unresolved input, positive LAPACK failure
71 : ! **************************************************************************************************
72 10794 : SUBROUTINE curvature_density(plaquettes, ndim, phase_limit, value, maximum, status)
73 : COMPLEX(KIND=dp), INTENT(IN) :: plaquettes(:, :, :)
74 : INTEGER, INTENT(IN) :: ndim
75 : REAL(KIND=dp), INTENT(IN) :: phase_limit
76 : REAL(KIND=dp), INTENT(OUT) :: value, maximum
77 : INTEGER, INTENT(OUT) :: status
78 :
79 10794 : COMPLEX(KIND=dp), ALLOCATABLE :: a(:, :), e(:), f(:, :, :), term(:, :), &
80 10794 : v(:, :), work(:)
81 : INTEGER :: i, j, n, np, sdim
82 10794 : LOGICAL, ALLOCATABLE :: bwork(:)
83 : REAL(KIND=dp) :: determinant_phase
84 10794 : REAL(KIND=dp), ALLOCATABLE :: phase(:), rwork(:)
85 :
86 10794 : status = -1
87 10794 : value = 0.0_dp
88 10794 : maximum = 0.0_dp
89 10794 : IF (ndim /= 2 .AND. ndim /= 4) RETURN
90 10794 : IF (.NOT. ieee_is_finite(phase_limit)) RETURN
91 10794 : IF (phase_limit <= 0.0_dp .OR. phase_limit >= pi) RETURN
92 10794 : n = SIZE(plaquettes, 1)
93 10794 : np = ndim*(ndim - 1)/2
94 10794 : IF (n < 1 .OR. SIZE(plaquettes, 2) /= n .OR. SIZE(plaquettes, 3) /= np) RETURN
95 464072 : IF (.NOT. ALL(ieee_is_finite(REAL(plaquettes, dp)))) RETURN
96 464072 : IF (.NOT. ALL(ieee_is_finite(AIMAG(plaquettes)))) RETURN
97 226674 : ALLOCATE (a(n, n), v(n, n), e(n), work(4*n), rwork(n), phase(n), bwork(n), f(n, n, np), term(n, n))
98 75524 : DO j = 1, np
99 453138 : a(:, :) = plaquettes(:, :, j)
100 971010 : term(:, :) = MATMUL(CONJG(TRANSPOSE(a)), a)
101 194202 : DO i = 1, n
102 194202 : term(i, i) = term(i, i) - 1.0_dp
103 : END DO
104 453138 : IF (MAXVAL(ABS(term)) > unitary_tol) THEN
105 2 : status = -2
106 2 : RETURN
107 : END IF
108 64732 : CALL zgees('V', 'N', select_none, n, a, n, sdim, e, v, n, work, SIZE(work), rwork, bwork, status)
109 64732 : IF (status /= 0) RETURN
110 194196 : phase(:) = ATAN2(AIMAG(e), REAL(e, dp))
111 64732 : IF (ndim == 2) THEN
112 10 : determinant_phase = ATAN2(SIN(SUM(phase)), COS(SUM(phase)))
113 2 : maximum = ABS(determinant_phase)
114 2 : value = -determinant_phase/(2.0_dp*pi)
115 : ELSE
116 194190 : maximum = MAX(maximum, MAXVAL(ABS(phase)))
117 453110 : a(:, :) = v
118 194190 : DO i = 1, n
119 453110 : a(:, i) = CMPLX(0.0_dp, phase(i), dp)*v(:, i)
120 : END DO
121 1359330 : f(:, :, j) = MATMUL(a, CONJG(TRANSPOSE(v)))
122 : END IF
123 75522 : IF (maximum >= phase_limit) THEN
124 2 : status = -3
125 2 : RETURN
126 : END IF
127 : END DO
128 10790 : IF (ndim == 4) THEN
129 614916 : term(:, :) = MATMUL(f(:, :, 1), f(:, :, 6)) - MATMUL(f(:, :, 2), f(:, :, 5)) + MATMUL(f(:, :, 3), f(:, :, 4))
130 32364 : DO i = 1, n
131 32364 : value = value - REAL(term(i, i), dp)/(4.0_dp*pi*pi)
132 : END DO
133 : END IF
134 10790 : status = 0
135 10794 : END SUBROUTINE curvature_density
136 :
137 : ! **************************************************************************************************
138 : !> \brief Unused eigenvalue-selection callback for an unsorted complex Schur decomposition.
139 : !> \param z Eigenvalue supplied by LAPACK
140 : !> \return Always false; no eigenvalue reordering is requested
141 : ! **************************************************************************************************
142 0 : LOGICAL FUNCTION select_none(z)
143 : COMPLEX(KIND=dp), INTENT(IN) :: z
144 :
145 : select_none = .FALSE.
146 : IF (.NOT. ieee_is_finite(REAL(z, dp))) select_none = .FALSE.
147 0 : END FUNCTION select_none
148 75492 : END MODULE topology_curvature
|
