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 2 : PROGRAM topology_symmetry_unittest
8 2 : USE ieee_arithmetic, ONLY: ieee_quiet_nan,&
9 : ieee_value
10 : USE iso_fortran_env, ONLY: output_unit
11 : USE kinds, ONLY: dp
12 : USE topology_symmetry, ONLY: character_multiplicities,&
13 : inversion_representation
14 : USE topology_tqc, ONLY: inversion_ebr_signature,&
15 : inversion_indicators
16 :
17 : IMPLICIT NONE
18 : COMPLEX(KIND=dp) :: table(6, 3), group_chars(6), metric(2, 2), coeff(4, 4), phase(2), tmp(4)
19 : REAL(KIND=dp) :: error, energy_error, energies(4)
20 : INTEGER :: multiplicities(3), status, counts(2), odd(8), ebr(2, 8), a, k, sign, strong, weak(3), z4, io_unit
21 : LOGICAL :: signed_atomic, nonnegative_atomic
22 :
23 2 : io_unit = output_unit
24 : ! Complete C3v/S3 character table per element, including a two-dimensional irrep.
25 14 : table(:, 1) = CMPLX([1, 1, 1, 1, 1, 1], 0, dp)
26 14 : table(:, 2) = CMPLX([1, 1, 1, -1, -1, -1], 0, dp)
27 14 : table(:, 3) = CMPLX([2, -1, -1, 0, 0, 0], 0, dp)
28 14 : group_chars = table(:, 1) + 2*table(:, 3)
29 2 : CALL character_multiplicities(group_chars, table, multiplicities, 1.e-10_dp, status)
30 8 : IF (status /= 0 .OR. ANY(multiplicities /= [1, 0, 2])) ERROR STOP 'C3v decomposition failed'
31 2 : group_chars(2) = group_chars(2) + 0.1_dp
32 2 : CALL character_multiplicities(group_chars, table, multiplicities, 1.e-10_dp, status)
33 2 : IF (status == 0) ERROR STOP 'Invalid character accepted'
34 2 : metric = CMPLX(0, 0, dp)
35 2 : metric(1, 1) = 2.0_dp
36 2 : metric(2, 2) = 0.5_dp
37 2 : coeff = CMPLX(0, 0, dp)
38 2 : coeff(1, 1) = SQRT(0.5_dp)
39 2 : coeff(2, 2) = SQRT(2.0_dp)
40 2 : coeff(3, 3) = SQRT(0.5_dp)
41 2 : coeff(4, 4) = SQRT(2.0_dp)
42 2 : phase = CMPLX([1, -1], 0, dp)
43 2 : energies = 0.0_dp
44 : CALL inversion_representation(metric, coeff, [1, 2], phase, energies, .TRUE., &
45 2 : 1.e-9_dp, 1.e-9_dp, counts, error, energy_error, status)
46 6 : IF (status /= 0 .OR. ANY(counts /= [2, 2])) ERROR STOP 'Metric/spinor parity failed'
47 10 : tmp = coeff(:, 1)
48 10 : coeff(:, 1) = (tmp + CMPLX(0, 1, dp)*coeff(:, 2))/SQRT(2.0_dp)
49 10 : coeff(:, 2) = (CMPLX(0, 1, dp)*tmp + coeff(:, 2))/SQRT(2.0_dp)
50 : CALL inversion_representation(metric, coeff, [1, 2], phase, energies, .TRUE., &
51 2 : 1.e-9_dp, 1.e-9_dp, counts, error, energy_error, status)
52 6 : IF (status /= 0 .OR. ANY(counts /= [2, 2])) ERROR STOP 'Degenerate-gauge parity failed'
53 : CALL inversion_representation(metric, coeff(:, :1), [1, 2], phase, energies(:1), .FALSE., &
54 2 : 1.e-9_dp, 1.e-9_dp, counts, error, energy_error, status)
55 2 : IF (status == 0) ERROR STOP 'Nonclosed parity subspace accepted'
56 : CALL inversion_representation(metric, coeff(:, :2), [1, 2], phase, energies(:2), .TRUE., &
57 2 : 1.e-9_dp, 1.e-9_dp, counts, error, energy_error, status)
58 2 : IF (status == 0) ERROR STOP 'Broken physical time reversal accepted'
59 2 : energies(2) = 1.0_dp
60 : CALL inversion_representation(metric, coeff, [1, 2], phase, energies, .TRUE., &
61 2 : 1.e-9_dp, 1.e-9_dp, counts, error, energy_error, status)
62 2 : IF (status == 0) ERROR STOP 'Noncommuting eigenspace accepted'
63 : ! Scalar two-site inversion with a genuinely off-diagonal AO metric.
64 14 : metric = CMPLX(0.2_dp, 0, dp)
65 2 : metric(1, 1) = 1.0_dp
66 2 : metric(2, 2) = 1.0_dp
67 6 : coeff(:2, 1) = CMPLX([1.0_dp, 1.0_dp]/SQRT(2.4_dp), 0, dp)
68 6 : coeff(:2, 2) = CMPLX([1.0_dp, -1.0_dp]/SQRT(1.6_dp), 0, dp)
69 6 : phase = CMPLX(1, 0, dp)
70 : CALL inversion_representation(metric, coeff(:2, :2), [2, 1], phase, energies(:2), .FALSE., &
71 2 : 1.e-9_dp, 1.e-9_dp, counts, error, energy_error, status)
72 6 : IF (status /= 0 .OR. ANY(counts /= [1, 1])) ERROR STOP 'Scalar two-site inversion failed'
73 : CALL inversion_representation(metric, coeff(:2, :1), [2, 1], -phase, energies(:1), .FALSE., &
74 6 : 1.e-9_dp, 1.e-9_dp, counts, error, energy_error, status)
75 6 : IF (status /= 0 .OR. ANY(counts /= [0, 1])) ERROR STOP 'Zone-boundary inversion phase failed'
76 2 : metric(1, 2) = CMPLX(0.2_dp, 0.01_dp, dp)
77 : CALL inversion_representation(metric, coeff(:2, :2), [2, 1], phase, energies(:2), .FALSE., &
78 2 : 1.e-9_dp, 1.e-9_dp, counts, error, energy_error, status)
79 2 : IF (status == 0) ERROR STOP 'Invalid metric accepted'
80 2 : metric(1, 2) = CMPLX(ieee_value(0.0_dp, ieee_quiet_nan), 0, dp)
81 : CALL inversion_representation(metric, coeff(:2, :2), [2, 1], phase, energies(:2), .FALSE., &
82 2 : 1.e-9_dp, 1.e-9_dp, counts, error, energy_error, status)
83 2 : IF (status == 0) ERROR STOP 'Nonfinite metric accepted'
84 18 : DO a = 0, 7
85 34 : DO sign = -1, 1, 2
86 288 : DO k = 0, 7
87 288 : odd(k + 1) = (1 - sign*(-1)**POPCNT(IAND(a, k)))/2
88 : END DO
89 32 : CALL inversion_indicators(odd, 1, 3, strong, weak, z4, status)
90 128 : IF (status /= 0 .OR. strong /= 0 .OR. ANY(weak /= 0) .OR. z4 /= 0) ERROR STOP 'Atomic indicator failed'
91 32 : CALL inversion_ebr_signature(odd, 1, signed_atomic, nonnegative_atomic, ebr, status)
92 32 : IF (status /= 0 .OR. .NOT. nonnegative_atomic .OR. .NOT. signed_atomic) ERROR STOP 'Atomic EBR failed'
93 816 : IF (SUM(ebr) /= 1) ERROR STOP 'EBR rank failed'
94 : END DO
95 : END DO
96 2 : odd = 0
97 2 : odd(1) = 1
98 2 : CALL inversion_indicators(odd, 1, 3, strong, weak, z4, status)
99 2 : IF (status /= 0 .OR. strong /= 1 .OR. z4 /= 1) ERROR STOP 'Strong index failed'
100 2 : CALL inversion_ebr_signature(odd, 1, signed_atomic, nonnegative_atomic, ebr, status)
101 2 : IF (status /= 0 .OR. signed_atomic .OR. nonnegative_atomic) ERROR STOP 'Stable obstruction missed'
102 2 : odd(1) = 2
103 2 : CALL inversion_indicators(odd, 2, 3, strong, weak, z4, status)
104 2 : IF (status /= 0 .OR. strong /= 0 .OR. z4 /= 2) ERROR STOP 'Inversion Z4=2 missed'
105 2 : odd(1) = 4
106 2 : CALL inversion_ebr_signature(odd, 4, signed_atomic, nonnegative_atomic, ebr, status)
107 2 : IF (status /= 0 .OR. .NOT. signed_atomic .OR. nonnegative_atomic) ERROR STOP 'Signed-only EBR missed'
108 18 : odd = -1
109 2 : CALL inversion_ebr_signature(odd, 4, signed_atomic, nonnegative_atomic, ebr, status)
110 2 : IF (status == 0) ERROR STOP 'Invalid multiplicity accepted'
111 2 : CALL exhaustive_signatures(2)
112 2 : CALL exhaustive_signatures(3)
113 2 : WRITE (io_unit, *) 'Native characters, parity, inversion indicators and EBR tests passed.'
114 : CONTAINS
115 :
116 : ! **************************************************************************************************
117 : !> \brief Compare the exact solver against independent enumeration of rank-two atomic sums.
118 : !> \param DIMENSION Dimension of the inversion subgroup
119 : ! **************************************************************************************************
120 4 : SUBROUTINE exhaustive_signatures(DIMENSION)
121 : INTEGER, INTENT(IN) :: dimension
122 :
123 : INTEGER :: code, digit, i, ierr, indicator, indices(3), j, n, nu, x, y, &
124 4 : values(2**DIMENSION), signatures(2**DIMENSION, 2**(DIMENSION + 1)), &
125 8 : reconstructed(2**DIMENSION), coefficients(2, 2**DIMENSION)
126 4 : LOGICAL :: indicators_zero, nonnegative_ok, possible(0:3**(2**DIMENSION) - 1), signed_ok
127 :
128 4 : n = 2**DIMENSION
129 28 : DO x = 0, n - 1
130 188 : DO i = 0, n - 1
131 160 : signatures(i + 1, 2*x + 1) = (1 - (-1)**POPCNT(IAND(x, i)))/2
132 184 : signatures(i + 1, 2*x + 2) = 1 - signatures(i + 1, 2*x + 1)
133 : END DO
134 : END DO
135 13288 : possible = .FALSE.
136 52 : DO x = 1, 2*n
137 692 : DO y = 1, 2*n
138 9856 : code = SUM((signatures(:, x) + signatures(:, y))*[(3**i, i=0, n - 1)])
139 688 : possible(code) = .TRUE.
140 : END DO
141 : END DO
142 13288 : DO code = 0, 3**n - 1
143 : digit = code
144 118908 : DO i = 1, n
145 105624 : values(i) = MOD(digit, 3)
146 118908 : digit = digit/3
147 : END DO
148 13284 : CALL inversion_ebr_signature(values, 2, signed_ok, nonnegative_ok, coefficients, ierr)
149 13284 : IF (ierr /= 0 .OR. (nonnegative_ok .NEQV. possible(code))) ERROR STOP 'Enumerated EBR mismatch'
150 13284 : CALL inversion_indicators(values, 2, DIMENSION, nu, indices, indicator, ierr)
151 13284 : IF (ierr /= 0) ERROR STOP 'Enumerated indicator failure'
152 13284 : indicators_zero = nu == 0
153 25110 : IF (DIMENSION == 3) indicators_zero = indicator == 0 .AND. ALL(indices == 0)
154 13284 : IF (signed_ok .NEQV. indicators_zero) ERROR STOP 'Indicator/atomic-lattice mismatch'
155 13284 : IF (.NOT. nonnegative_ok) CYCLE
156 2652 : reconstructed = 0
157 2652 : DO j = 1, n
158 : reconstructed = reconstructed + coefficients(1, j)*signatures(:, 2*j - 1) + &
159 20220 : coefficients(2, j)*signatures(:, 2*j)
160 : END DO
161 9964 : IF (ANY(reconstructed /= values) .OR. SUM(coefficients) /= 2) ERROR STOP 'EBR reconstruction failed'
162 : END DO
163 2 : END SUBROUTINE exhaustive_signatures
164 : END PROGRAM topology_symmetry_unittest
|
