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 Native character decomposition and metric-aware inversion representations.
10 : ! **************************************************************************************************
11 : MODULE topology_symmetry
12 : USE ieee_arithmetic, ONLY: ieee_is_finite
13 : USE kinds, ONLY: dp
14 :
15 : IMPLICIT NONE
16 : PRIVATE
17 : PUBLIC :: character_multiplicities, inversion_representation
18 : PUBLIC :: spin_rotation, transform_spinor_frame
19 : CONTAINS
20 :
21 : ! **************************************************************************************************
22 : !> \brief Lift a Cartesian orthogonal operation to the physical spin-1/2 representation.
23 : !> \param rotation Cartesian polar-vector operation, including improper operations
24 : !> \param spin SU(2) lift of det(rotation)*rotation; inversion acts trivially on spin
25 : !> \param tolerance Orthogonality tolerance
26 : !> \param status Zero on success, negative for invalid geometry
27 : !> \note The quaternion sign fixes one representative, not a single-valued SO(3) representation.
28 : !> Products must retain the resulting double-group factor system. Orbital parity is separate.
29 : ! **************************************************************************************************
30 12 : SUBROUTINE spin_rotation(rotation, spin, tolerance, status)
31 : REAL(dp), INTENT(IN) :: rotation(3, 3)
32 : COMPLEX(dp), INTENT(OUT) :: spin(2, 2)
33 : REAL(dp), INTENT(IN) :: tolerance
34 : INTEGER, INTENT(OUT) :: status
35 :
36 : INTEGER :: branch, i, j, k
37 : REAL(dp) :: candidates(0:3), det, q(0:3), r(3, 3), &
38 : scale, work(3, 3)
39 :
40 12 : spin = 0.0_dp
41 12 : status = -1
42 156 : IF (.NOT. ALL(ieee_is_finite(rotation)) .OR. .NOT. ieee_is_finite(tolerance)) RETURN
43 12 : IF (tolerance <= 0.0_dp .OR. tolerance >= 1.0_dp) RETURN
44 480 : work = MATMUL(TRANSPOSE(rotation), rotation)
45 48 : DO i = 1, 3
46 48 : work(i, i) = work(i, i) - 1.0_dp
47 : END DO
48 156 : IF (MAXVAL(ABS(work)) > tolerance) RETURN
49 : det = rotation(1, 1)*(rotation(2, 2)*rotation(3, 3) - rotation(2, 3)*rotation(3, 2)) &
50 : - rotation(1, 2)*(rotation(2, 1)*rotation(3, 3) - rotation(2, 3)*rotation(3, 1)) &
51 12 : + rotation(1, 3)*(rotation(2, 1)*rotation(3, 2) - rotation(2, 2)*rotation(3, 1))
52 12 : IF (ABS(ABS(det) - 1.0_dp) > 3*tolerance) RETURN
53 156 : r = SIGN(1.0_dp, det)*rotation
54 12 : candidates(0) = 1.0_dp + r(1, 1) + r(2, 2) + r(3, 3)
55 48 : DO i = 1, 3
56 48 : candidates(i) = 2.0_dp*(1.0_dp + r(i, i)) - candidates(0)
57 : END DO
58 60 : branch = MAXLOC(candidates, DIM=1) - 1
59 12 : scale = 2.0_dp*SQRT(MAX(0.0_dp, candidates(branch)))
60 12 : q(branch) = 0.25_dp*scale
61 12 : IF (branch == 0) THEN
62 4 : q(1) = (r(3, 2) - r(2, 3))/scale
63 4 : q(2) = (r(1, 3) - r(3, 1))/scale
64 4 : q(3) = (r(2, 1) - r(1, 2))/scale
65 : ELSE
66 8 : i = branch
67 8 : j = MODULO(i, 3) + 1
68 8 : k = MODULO(i + 1, 3) + 1
69 8 : q(0) = (r(k, j) - r(j, k))/scale
70 8 : q(j) = (r(j, i) + r(i, j))/scale
71 8 : q(k) = (r(k, i) + r(i, k))/scale
72 : END IF
73 108 : q = q/SQRT(SUM(q**2))
74 12 : spin(1, 1) = CMPLX(q(0), -q(3), dp)
75 12 : spin(2, 2) = CONJG(spin(1, 1))
76 12 : spin(1, 2) = CMPLX(-q(2), -q(1), dp)
77 12 : spin(2, 1) = -CONJG(spin(1, 2))
78 12 : status = 0
79 : END SUBROUTINE spin_rotation
80 :
81 : ! **************************************************************************************************
82 : !> \brief Compose orbital, Bloch, physical-spin and optional time-reversal transformations.
83 : !> \param orbital Scalar AO action from the source to the target k, including all lattice phases
84 : !> \param spin Unitary spin matrix from spin_rotation
85 : !> \param source Source spin-major frame
86 : !> \param antiunitary Apply physical Theta=i*sigma_y*K before the spatial operation
87 : !> \param image Transformed spin-major frame in the target AO basis
88 : !> \param tolerance Spin unitarity tolerance
89 : !> \param status Zero on success, negative for malformed data
90 : !> \note Orbital metric/Hamiltonian covariance and k mapping must be checked by the caller.
91 : ! **************************************************************************************************
92 12 : SUBROUTINE transform_spinor_frame(orbital, spin, source, antiunitary, image, tolerance, status)
93 : COMPLEX(dp), INTENT(IN) :: orbital(:, :), spin(2, 2), source(:, :)
94 : LOGICAL, INTENT(IN) :: antiunitary
95 : COMPLEX(dp), INTENT(OUT) :: image(:, :)
96 : REAL(dp), INTENT(IN) :: tolerance
97 : INTEGER, INTENT(OUT) :: status
98 :
99 : COMPLEX(dp) :: gram(2, 2)
100 12 : COMPLEX(dp), ALLOCATABLE :: down(:, :), up(:, :)
101 : INTEGER :: m, n
102 :
103 12 : status = -1
104 28976 : image = 0.0_dp
105 12 : n = SIZE(orbital, 2)
106 12 : m = SIZE(orbital, 1)
107 12 : IF (n < 1 .OR. m < 1 .OR. SIZE(source, 1) /= 2*n) RETURN
108 36 : IF (ANY(SHAPE(image) /= [2*m, SIZE(source, 2)])) RETURN
109 12 : IF (.NOT. ieee_is_finite(tolerance) .OR. tolerance <= 0.0_dp .OR. tolerance >= 1.0_dp) RETURN
110 14792 : IF (.NOT. ALL(ieee_is_finite(REAL(orbital, dp))) .OR. .NOT. ALL(ieee_is_finite(AIMAG(orbital)))) RETURN
111 57952 : IF (.NOT. ALL(ieee_is_finite(REAL(source, dp))) .OR. .NOT. ALL(ieee_is_finite(AIMAG(source)))) RETURN
112 168 : IF (.NOT. ALL(ieee_is_finite(REAL(spin, dp))) .OR. .NOT. ALL(ieee_is_finite(AIMAG(spin)))) RETURN
113 180 : gram = MATMUL(CONJG(TRANSPOSE(spin)), spin)
114 12 : gram(1, 1) = gram(1, 1) - 1.0_dp
115 12 : gram(2, 2) = gram(2, 2) - 1.0_dp
116 84 : IF (MAXVAL(ABS(gram)) > tolerance) RETURN
117 72 : ALLOCATE (up(m, SIZE(source, 2)), down(m, SIZE(source, 2)))
118 12 : IF (antiunitary) THEN
119 9518 : up(:, :) = MATMUL(orbital, CONJG(source(n + 1:, :)))
120 10254 : down(:, :) = -MATMUL(orbital, CONJG(source(:n, :)))
121 : ELSE
122 10 : up(:, :) = MATMUL(orbital, source(:n, :))
123 10 : down(:, :) = MATMUL(orbital, source(n + 1:, :))
124 : END IF
125 14780 : image(:m, :) = spin(1, 1)*up + spin(1, 2)*down
126 14780 : image(m + 1:, :) = spin(2, 1)*up + spin(2, 2)*down
127 57952 : IF (.NOT. ALL(ieee_is_finite(REAL(image, dp))) .OR. .NOT. ALL(ieee_is_finite(AIMAG(image)))) RETURN
128 12 : status = 0
129 12 : END SUBROUTINE transform_spinor_frame
130 :
131 : ! **************************************************************************************************
132 : !> \brief Decompose a character against a complete unitary character table.
133 : !> \param CHARACTER One character per group element (not per conjugacy class)
134 : !> \param table Characters, group element first and irrep second; identity first
135 : !> \param multiplicities Nonnegative integer multiplicities, invalid unless status=0
136 : !> \param tolerance Absolute numerical tolerance
137 : !> \param status Zero on success, negative for invalid/incomplete/nonintegral data
138 : !> \note All projective irreps must use the same factor system. This kernel does not
139 : !> infer a space group or certify the supplied character table's provenance.
140 : ! **************************************************************************************************
141 56 : SUBROUTINE character_multiplicities(CHARACTER, table, multiplicities, tolerance, status)
142 : COMPLEX(KIND=dp), INTENT(IN) :: character(:), table(:, :)
143 : INTEGER, INTENT(OUT) :: multiplicities(:)
144 : REAL(KIND=dp), INTENT(IN) :: tolerance
145 : INTEGER, INTENT(OUT) :: status
146 :
147 : COMPLEX(KIND=dp) :: value
148 : INTEGER :: dimension, i, j, nirrep, order, total
149 :
150 56 : status = -1
151 172 : multiplicities = -1
152 56 : order = SIZE(CHARACTER)
153 56 : nirrep = SIZE(table, 2)
154 56 : IF (order < 1 .OR. order > 1024 .OR. nirrep < 1 .OR. SIZE(table, 1) /= order) RETURN
155 56 : IF (SIZE(multiplicities) /= nirrep .OR. tolerance <= 0.0_dp) RETURN
156 56 : IF (.NOT. ieee_is_finite(tolerance)) RETURN
157 184 : IF (.NOT. ALL(ieee_is_finite(REAL(CHARACTER, dp)))) RETURN
158 184 : IF (.NOT. ALL(ieee_is_finite(AIMAG(CHARACTER)))) RETURN
159 452 : IF (.NOT. ALL(ieee_is_finite(REAL(table, dp)))) RETURN
160 452 : IF (.NOT. ALL(ieee_is_finite(AIMAG(table)))) RETURN
161 452 : IF (MAXVAL(ABS(table)) > REAL(order, dp) + tolerance) RETURN
162 184 : IF (MAXVAL(ABS(CHARACTER)) > 1.e6_dp) RETURN
163 : total = 0
164 172 : DO i = 1, nirrep
165 116 : DIMENSION = NINT(REAL(table(1, i), dp))
166 116 : IF (DIMENSION < 1 .OR. DIMENSION > order) RETURN
167 116 : IF (ABS(table(1, i) - DIMENSION) > tolerance) RETURN
168 116 : total = total + DIMENSION**2
169 116 : IF (total > order) RETURN
170 416 : DO j = 1, nirrep
171 876 : value = SUM(CONJG(table(:, i))*table(:, j))/REAL(order, dp)
172 244 : IF (i == j) value = value - 1.0_dp
173 360 : IF (ABS(value) > tolerance) RETURN
174 : END DO
175 : END DO
176 56 : IF (total /= order) RETURN
177 56 : status = -2
178 158 : DO i = 1, nirrep
179 356 : value = SUM(CONJG(table(:, i))*CHARACTER)/REAL(order, dp)
180 108 : j = NINT(REAL(value, dp))
181 108 : IF (j < 0 .OR. ABS(value - j) > tolerance) RETURN
182 152 : multiplicities(i) = j
183 : END DO
184 690 : IF (MAXVAL(ABS(MATMUL(table, CMPLX(multiplicities, 0, dp)) - CHARACTER)) > tolerance) RETURN
185 50 : status = 0
186 : END SUBROUTINE character_multiplicities
187 :
188 : ! **************************************************************************************************
189 : !> \brief Inversion irreps of an isolated scalar or spinor eigenspace in an AO metric.
190 : !> \param metric Scalar AO overlap at a TRIM
191 : !> \param coeff Selected coefficients, spinor components stacked by AO
192 : !> \param mapping Target AO index under inversion
193 : !> \param phase Source-AO multiplier, including orbital parity and lattice phase
194 : !> \param energies Selected eigenvalues in Hartree
195 : !> \param check_tr Check physical spinful time reversal as well as inversion
196 : !> \param tolerance Dimensionless metric/subspace/character tolerance
197 : !> \param energy_tolerance Eigenvalue-commutator tolerance in Hartree
198 : !> \param counts Multiplicities of even/odd inversion irreps, counting states not pairs
199 : !> \param error Largest dimensionless residual
200 : !> \param energy_error Largest projected symmetry/eigenvalue commutator
201 : !> \param status Zero on success; invalid output for nonzero status
202 : !> \param diagnostics Optional metric, normalization, inversion and time-reversal residuals
203 : ! **************************************************************************************************
204 54 : SUBROUTINE inversion_representation(metric, coeff, mapping, phase, energies, check_tr, &
205 : tolerance, energy_tolerance, counts, error, energy_error, status, diagnostics)
206 : COMPLEX(KIND=dp), INTENT(IN) :: metric(:, :), coeff(:, :)
207 : INTEGER, INTENT(IN) :: mapping(:)
208 : COMPLEX(KIND=dp), INTENT(IN) :: phase(:)
209 : REAL(KIND=dp), INTENT(IN) :: energies(:)
210 : LOGICAL, INTENT(IN) :: check_tr
211 : REAL(KIND=dp), INTENT(IN) :: tolerance, energy_tolerance
212 : INTEGER, INTENT(OUT) :: counts(2)
213 : REAL(KIND=dp), INTENT(OUT) :: error, energy_error
214 : INTEGER, INTENT(OUT) :: status
215 : REAL(KIND=dp), INTENT(OUT), OPTIONAL :: diagnostics(4)
216 :
217 : COMPLEX(KIND=dp) :: characters(2), table(2, 2)
218 54 : COMPLEX(KIND=dp), ALLOCATABLE :: pc(:, :), projected(:, :), sc(:, :), &
219 54 : trc(:, :), work(:, :)
220 : INTEGER :: first, i, j, last, nao, nspin, rank, s
221 : REAL(KIND=dp) :: metric_scale
222 :
223 54 : status = -1
224 162 : counts = -1
225 54 : error = HUGE(1.0_dp)
226 54 : energy_error = HUGE(1.0_dp)
227 198 : IF (PRESENT(diagnostics)) diagnostics = HUGE(1.0_dp)
228 54 : nao = SIZE(metric, 1)
229 54 : rank = SIZE(coeff, 2)
230 54 : IF (nao < 1 .OR. rank < 1 .OR. SIZE(metric, 2) /= nao) RETURN
231 54 : IF (MOD(SIZE(coeff, 1), nao) /= 0) RETURN
232 54 : nspin = SIZE(coeff, 1)/nao
233 54 : IF (nspin < 1 .OR. nspin > 2) RETURN
234 54 : IF (check_tr .AND. (nspin /= 2 .OR. MOD(rank, 2) /= 0)) RETURN
235 54 : IF (SIZE(mapping) /= nao .OR. SIZE(phase) /= nao .OR. SIZE(energies) /= rank) RETURN
236 54 : IF (tolerance <= 0.0_dp .OR. energy_tolerance <= 0.0_dp) RETURN
237 54 : IF (.NOT. ieee_is_finite(tolerance) .OR. .NOT. ieee_is_finite(energy_tolerance)) RETURN
238 20526 : IF (.NOT. ALL(ieee_is_finite(REAL(metric, dp))) .OR. .NOT. ALL(ieee_is_finite(AIMAG(metric)))) RETURN
239 17664 : IF (.NOT. ALL(ieee_is_finite(REAL(coeff, dp))) .OR. .NOT. ALL(ieee_is_finite(AIMAG(coeff)))) RETURN
240 1248 : IF (.NOT. ALL(ieee_is_finite(REAL(phase, dp))) .OR. .NOT. ALL(ieee_is_finite(AIMAG(phase)))) RETURN
241 352 : IF (.NOT. ALL(ieee_is_finite(energies))) RETURN
242 1248 : IF (ANY(mapping < 1) .OR. ANY(mapping > nao)) RETURN
243 624 : IF (ANY(ABS(ABS(phase) - 1.0_dp) > tolerance)) RETURN
244 624 : DO i = 1, nao
245 572 : IF (mapping(mapping(i)) /= i) RETURN
246 624 : IF (ABS(phase(i)*phase(mapping(i)) - 1.0_dp) > tolerance) RETURN
247 : END DO
248 10252 : metric_scale = MAX(1.0_dp, MAXVAL(ABS(metric)))
249 10252 : error = MAXVAL(ABS(metric - CONJG(TRANSPOSE(metric))))/metric_scale
250 624 : DO j = 1, nao
251 10252 : DO i = 1, nao
252 : error = MAX(error, ABS(CONJG(phase(i))*metric(mapping(i), mapping(j))*phase(j) - metric(i, j))/ &
253 10200 : metric_scale)
254 : END DO
255 : END DO
256 52 : IF (PRESENT(diagnostics)) diagnostics(1) = error
257 572 : ALLOCATE (sc(nao*nspin, rank), pc(nao*nspin, rank), projected(rank, rank), work(rank, rank))
258 146 : DO s = 1, nspin
259 94 : first = (s - 1)*nao
260 94 : last = s*nao
261 170030 : sc(first + 1:last, :) = MATMUL(metric, coeff(first + 1:last, :))
262 1258 : DO i = 1, nao
263 9686 : pc(first + mapping(i), :) = phase(i)*coeff(first + i, :)
264 : END DO
265 : END DO
266 69568 : work(:, :) = MATMUL(CONJG(TRANSPOSE(coeff)), sc)
267 352 : DO i = 1, rank
268 352 : work(i, i) = work(i, i) - 1.0_dp
269 : END DO
270 2528 : error = MAX(error, MAXVAL(ABS(work)))
271 2364 : IF (PRESENT(diagnostics)) diagnostics(2) = MAXVAL(ABS(work))
272 69568 : projected(:, :) = MATMUL(CONJG(TRANSPOSE(sc)), pc)
273 2528 : IF (.NOT. ALL(ieee_is_finite(REAL(projected, dp)))) RETURN
274 2528 : IF (.NOT. ALL(ieee_is_finite(AIMAG(projected)))) RETURN
275 19352 : work(:, :) = MATMUL(CONJG(TRANSPOSE(projected)), projected)
276 352 : DO i = 1, rank
277 352 : work(i, i) = work(i, i) - 1.0_dp
278 : END DO
279 5056 : error = MAX(error, MAXVAL(ABS(work)), MAXVAL(ABS(projected - CONJG(TRANSPOSE(projected)))))
280 52 : IF (PRESENT(diagnostics)) diagnostics(3) = MAX(MAXVAL(ABS(work)), &
281 4696 : MAXVAL(ABS(projected - CONJG(TRANSPOSE(projected)))))
282 : IF (PRESENT(diagnostics)) THEN
283 36 : IF (.NOT. check_tr) diagnostics(4) = 0.0_dp
284 : END IF
285 52 : energy_error = 0.0_dp
286 352 : DO j = 1, rank
287 2528 : DO i = 1, rank
288 2476 : energy_error = MAX(energy_error, ABS(projected(i, j)*(energies(i) - energies(j))))
289 : END DO
290 : END DO
291 52 : characters(1) = CMPLX(rank, 0, dp)
292 52 : characters(2) = CMPLX(0, 0, dp)
293 352 : DO i = 1, rank
294 352 : characters(2) = characters(2) + projected(i, i)
295 : END DO
296 156 : table(:, 1) = CMPLX([1, 1], 0, dp)
297 156 : table(:, 2) = CMPLX([1, -1], 0, dp)
298 52 : CALL character_multiplicities(characters, table, counts, tolerance*rank, status)
299 52 : IF (status /= 0) RETURN
300 48 : status = -3
301 48 : IF (check_tr) THEN
302 : ! In a real Gaussian Bloch basis at a TRIM, Theta=(i sigma_y) K.
303 20144 : error = MAX(error, MAXVAL(ABS(AIMAG(metric)))/MAX(1.0_dp, MAXVAL(ABS(metric))))
304 160 : ALLOCATE (trc(2*nao, rank))
305 4540 : trc(:nao, :) = CONJG(coeff(nao + 1:, :))
306 4540 : trc(nao + 1:, :) = -CONJG(coeff(:nao, :))
307 69452 : projected(:, :) = MATMUL(CONJG(TRANSPOSE(sc)), trc)
308 19260 : work(:, :) = MATMUL(CONJG(TRANSPOSE(projected)), projected)
309 324 : DO i = 1, rank
310 324 : work(i, i) = work(i, i) - 1.0_dp
311 : END DO
312 4952 : error = MAX(error, MAXVAL(ABS(work)), MAXVAL(ABS(projected + TRANSPOSE(projected))))
313 40 : IF (PRESENT(diagnostics)) diagnostics(4) = MAX(MAXVAL(ABS(work)), &
314 4672 : MAXVAL(ABS(projected + TRANSPOSE(projected))))
315 324 : DO j = 1, rank
316 2476 : DO i = 1, rank
317 2436 : energy_error = MAX(energy_error, ABS(projected(i, j)*(energies(i) - energies(j))))
318 : END DO
319 : END DO
320 116 : IF (ANY(MOD(counts, 2) /= 0)) RETURN
321 : END IF
322 46 : IF (.NOT. ieee_is_finite(error) .OR. .NOT. ieee_is_finite(energy_error)) RETURN
323 46 : IF (error > tolerance .OR. energy_error > energy_tolerance) RETURN
324 44 : status = 0
325 54 : END SUBROUTINE inversion_representation
326 52 : END MODULE topology_symmetry
|
