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 : CONTAINS
19 :
20 : ! **************************************************************************************************
21 : !> \brief Decompose a character against a complete unitary character table.
22 : !> \param CHARACTER One character per group element (not per conjugacy class)
23 : !> \param table Characters, group element first and irrep second; identity first
24 : !> \param multiplicities Nonnegative integer multiplicities, invalid unless status=0
25 : !> \param tolerance Absolute numerical tolerance
26 : !> \param status Zero on success, negative for invalid/incomplete/nonintegral data
27 : !> \note All projective irreps must use the same factor system. This kernel does not
28 : !> infer a space group or certify the supplied character table's provenance.
29 : ! **************************************************************************************************
30 52 : SUBROUTINE character_multiplicities(CHARACTER, table, multiplicities, tolerance, status)
31 : COMPLEX(KIND=dp), INTENT(IN) :: character(:), table(:, :)
32 : INTEGER, INTENT(OUT) :: multiplicities(:)
33 : REAL(KIND=dp), INTENT(IN) :: tolerance
34 : INTEGER, INTENT(OUT) :: status
35 :
36 : COMPLEX(KIND=dp) :: value
37 : INTEGER :: dimension, i, j, nirrep, order, total
38 :
39 52 : status = -1
40 160 : multiplicities = -1
41 52 : order = SIZE(CHARACTER)
42 52 : nirrep = SIZE(table, 2)
43 52 : IF (order < 1 .OR. order > 1024 .OR. nirrep < 1 .OR. SIZE(table, 1) /= order) RETURN
44 52 : IF (SIZE(multiplicities) /= nirrep .OR. tolerance <= 0.0_dp) RETURN
45 52 : IF (.NOT. ieee_is_finite(tolerance)) RETURN
46 172 : IF (.NOT. ALL(ieee_is_finite(REAL(CHARACTER, dp)))) RETURN
47 172 : IF (.NOT. ALL(ieee_is_finite(AIMAG(CHARACTER)))) RETURN
48 424 : IF (.NOT. ALL(ieee_is_finite(REAL(table, dp)))) RETURN
49 424 : IF (.NOT. ALL(ieee_is_finite(AIMAG(table)))) RETURN
50 424 : IF (MAXVAL(ABS(table)) > REAL(order, dp) + tolerance) RETURN
51 172 : IF (MAXVAL(ABS(CHARACTER)) > 1.e6_dp) RETURN
52 : total = 0
53 160 : DO i = 1, nirrep
54 108 : DIMENSION = NINT(REAL(table(1, i), dp))
55 108 : IF (DIMENSION < 1 .OR. DIMENSION > order) RETURN
56 108 : IF (ABS(table(1, i) - DIMENSION) > tolerance) RETURN
57 108 : total = total + DIMENSION**2
58 108 : IF (total > order) RETURN
59 388 : DO j = 1, nirrep
60 828 : value = SUM(CONJG(table(:, i))*table(:, j))/REAL(order, dp)
61 228 : IF (i == j) value = value - 1.0_dp
62 336 : IF (ABS(value) > tolerance) RETURN
63 : END DO
64 : END DO
65 52 : IF (total /= order) RETURN
66 52 : status = -2
67 146 : DO i = 1, nirrep
68 332 : value = SUM(CONJG(table(:, i))*CHARACTER)/REAL(order, dp)
69 100 : j = NINT(REAL(value, dp))
70 100 : IF (j < 0 .OR. ABS(value - j) > tolerance) RETURN
71 140 : multiplicities(i) = j
72 : END DO
73 638 : IF (MAXVAL(ABS(MATMUL(table, CMPLX(multiplicities, 0, dp)) - CHARACTER)) > tolerance) RETURN
74 46 : status = 0
75 : END SUBROUTINE character_multiplicities
76 :
77 : ! **************************************************************************************************
78 : !> \brief Inversion irreps of an isolated scalar or spinor eigenspace in an AO metric.
79 : !> \param metric Scalar AO overlap at a TRIM
80 : !> \param coeff Selected coefficients, spinor components stacked by AO
81 : !> \param mapping Target AO index under inversion
82 : !> \param phase Source-AO multiplier, including orbital parity and lattice phase
83 : !> \param energies Selected eigenvalues in Hartree
84 : !> \param check_tr Check physical spinful time reversal as well as inversion
85 : !> \param tolerance Dimensionless metric/subspace/character tolerance
86 : !> \param energy_tolerance Eigenvalue-commutator tolerance in Hartree
87 : !> \param counts Multiplicities of even/odd inversion irreps, counting states not pairs
88 : !> \param error Largest dimensionless residual
89 : !> \param energy_error Largest projected symmetry/eigenvalue commutator
90 : !> \param status Zero on success; invalid output for nonzero status
91 : !> \param diagnostics Optional metric, normalization, inversion and time-reversal residuals
92 : ! **************************************************************************************************
93 50 : SUBROUTINE inversion_representation(metric, coeff, mapping, phase, energies, check_tr, &
94 : tolerance, energy_tolerance, counts, error, energy_error, status, diagnostics)
95 : COMPLEX(KIND=dp), INTENT(IN) :: metric(:, :), coeff(:, :)
96 : INTEGER, INTENT(IN) :: mapping(:)
97 : COMPLEX(KIND=dp), INTENT(IN) :: phase(:)
98 : REAL(KIND=dp), INTENT(IN) :: energies(:)
99 : LOGICAL, INTENT(IN) :: check_tr
100 : REAL(KIND=dp), INTENT(IN) :: tolerance, energy_tolerance
101 : INTEGER, INTENT(OUT) :: counts(2)
102 : REAL(KIND=dp), INTENT(OUT) :: error, energy_error
103 : INTEGER, INTENT(OUT) :: status
104 : REAL(KIND=dp), INTENT(OUT), OPTIONAL :: diagnostics(4)
105 :
106 : COMPLEX(KIND=dp) :: characters(2), table(2, 2)
107 50 : COMPLEX(KIND=dp), ALLOCATABLE :: pc(:, :), projected(:, :), sc(:, :), &
108 50 : trc(:, :), work(:, :)
109 : INTEGER :: first, i, j, last, nao, nspin, rank, s
110 : REAL(KIND=dp) :: metric_scale
111 :
112 50 : status = -1
113 150 : counts = -1
114 50 : error = HUGE(1.0_dp)
115 50 : energy_error = HUGE(1.0_dp)
116 178 : IF (PRESENT(diagnostics)) diagnostics = HUGE(1.0_dp)
117 50 : nao = SIZE(metric, 1)
118 50 : rank = SIZE(coeff, 2)
119 50 : IF (nao < 1 .OR. rank < 1 .OR. SIZE(metric, 2) /= nao) RETURN
120 50 : IF (MOD(SIZE(coeff, 1), nao) /= 0) RETURN
121 50 : nspin = SIZE(coeff, 1)/nao
122 50 : IF (nspin < 1 .OR. nspin > 2) RETURN
123 50 : IF (check_tr .AND. (nspin /= 2 .OR. MOD(rank, 2) /= 0)) RETURN
124 50 : IF (SIZE(mapping) /= nao .OR. SIZE(phase) /= nao .OR. SIZE(energies) /= rank) RETURN
125 50 : IF (tolerance <= 0.0_dp .OR. energy_tolerance <= 0.0_dp) RETURN
126 50 : IF (.NOT. ieee_is_finite(tolerance) .OR. .NOT. ieee_is_finite(energy_tolerance)) RETURN
127 20278 : IF (.NOT. ALL(ieee_is_finite(REAL(metric, dp))) .OR. .NOT. ALL(ieee_is_finite(AIMAG(metric)))) RETURN
128 17608 : IF (.NOT. ALL(ieee_is_finite(REAL(coeff, dp))) .OR. .NOT. ALL(ieee_is_finite(AIMAG(coeff)))) RETURN
129 1200 : IF (.NOT. ALL(ieee_is_finite(REAL(phase, dp))) .OR. .NOT. ALL(ieee_is_finite(AIMAG(phase)))) RETURN
130 344 : IF (.NOT. ALL(ieee_is_finite(energies))) RETURN
131 1200 : IF (ANY(mapping < 1) .OR. ANY(mapping > nao)) RETURN
132 600 : IF (ANY(ABS(ABS(phase) - 1.0_dp) > tolerance)) RETURN
133 600 : DO i = 1, nao
134 552 : IF (mapping(mapping(i)) /= i) RETURN
135 600 : IF (ABS(phase(i)*phase(mapping(i)) - 1.0_dp) > tolerance) RETURN
136 : END DO
137 10128 : metric_scale = MAX(1.0_dp, MAXVAL(ABS(metric)))
138 10128 : error = MAXVAL(ABS(metric - CONJG(TRANSPOSE(metric))))/metric_scale
139 600 : DO j = 1, nao
140 10128 : DO i = 1, nao
141 : error = MAX(error, ABS(CONJG(phase(i))*metric(mapping(i), mapping(j))*phase(j) - metric(i, j))/ &
142 10080 : metric_scale)
143 : END DO
144 : END DO
145 48 : IF (PRESENT(diagnostics)) diagnostics(1) = error
146 528 : ALLOCATE (sc(nao*nspin, rank), pc(nao*nspin, rank), projected(rank, rank), work(rank, rank))
147 138 : DO s = 1, nspin
148 90 : first = (s - 1)*nao
149 90 : last = s*nao
150 169878 : sc(first + 1:last, :) = MATMUL(metric, coeff(first + 1:last, :))
151 1230 : DO i = 1, nao
152 9642 : pc(first + mapping(i), :) = phase(i)*coeff(first + i, :)
153 : END DO
154 : END DO
155 69536 : work(:, :) = MATMUL(CONJG(TRANSPOSE(coeff)), sc)
156 344 : DO i = 1, rank
157 344 : work(i, i) = work(i, i) - 1.0_dp
158 : END DO
159 2516 : error = MAX(error, MAXVAL(ABS(work)))
160 2352 : IF (PRESENT(diagnostics)) diagnostics(2) = MAXVAL(ABS(work))
161 69536 : projected(:, :) = MATMUL(CONJG(TRANSPOSE(sc)), pc)
162 2516 : IF (.NOT. ALL(ieee_is_finite(REAL(projected, dp)))) RETURN
163 2516 : IF (.NOT. ALL(ieee_is_finite(AIMAG(projected)))) RETURN
164 19336 : work(:, :) = MATMUL(CONJG(TRANSPOSE(projected)), projected)
165 344 : DO i = 1, rank
166 344 : work(i, i) = work(i, i) - 1.0_dp
167 : END DO
168 5032 : error = MAX(error, MAXVAL(ABS(work)), MAXVAL(ABS(projected - CONJG(TRANSPOSE(projected)))))
169 48 : IF (PRESENT(diagnostics)) diagnostics(3) = MAX(MAXVAL(ABS(work)), &
170 4672 : MAXVAL(ABS(projected - CONJG(TRANSPOSE(projected)))))
171 : IF (PRESENT(diagnostics)) THEN
172 32 : IF (.NOT. check_tr) diagnostics(4) = 0.0_dp
173 : END IF
174 48 : energy_error = 0.0_dp
175 344 : DO j = 1, rank
176 2516 : DO i = 1, rank
177 2468 : energy_error = MAX(energy_error, ABS(projected(i, j)*(energies(i) - energies(j))))
178 : END DO
179 : END DO
180 48 : characters(1) = CMPLX(rank, 0, dp)
181 48 : characters(2) = CMPLX(0, 0, dp)
182 344 : DO i = 1, rank
183 344 : characters(2) = characters(2) + projected(i, i)
184 : END DO
185 144 : table(:, 1) = CMPLX([1, 1], 0, dp)
186 144 : table(:, 2) = CMPLX([1, -1], 0, dp)
187 48 : CALL character_multiplicities(characters, table, counts, tolerance*rank, status)
188 48 : IF (status /= 0) RETURN
189 44 : status = -3
190 44 : IF (check_tr) THEN
191 : ! In a real Gaussian Bloch basis at a TRIM, Theta=(i sigma_y) K.
192 20144 : error = MAX(error, MAXVAL(ABS(AIMAG(metric)))/MAX(1.0_dp, MAXVAL(ABS(metric))))
193 160 : ALLOCATE (trc(2*nao, rank))
194 4540 : trc(:nao, :) = CONJG(coeff(nao + 1:, :))
195 4540 : trc(nao + 1:, :) = -CONJG(coeff(:nao, :))
196 69452 : projected(:, :) = MATMUL(CONJG(TRANSPOSE(sc)), trc)
197 19260 : work(:, :) = MATMUL(CONJG(TRANSPOSE(projected)), projected)
198 324 : DO i = 1, rank
199 324 : work(i, i) = work(i, i) - 1.0_dp
200 : END DO
201 4952 : error = MAX(error, MAXVAL(ABS(work)), MAXVAL(ABS(projected + TRANSPOSE(projected))))
202 40 : IF (PRESENT(diagnostics)) diagnostics(4) = MAX(MAXVAL(ABS(work)), &
203 4672 : MAXVAL(ABS(projected + TRANSPOSE(projected))))
204 324 : DO j = 1, rank
205 2476 : DO i = 1, rank
206 2436 : energy_error = MAX(energy_error, ABS(projected(i, j)*(energies(i) - energies(j))))
207 : END DO
208 : END DO
209 116 : IF (ANY(MOD(counts, 2) /= 0)) RETURN
210 : END IF
211 42 : IF (.NOT. ieee_is_finite(error) .OR. .NOT. ieee_is_finite(energy_error)) RETURN
212 42 : IF (error > tolerance .OR. energy_error > energy_tolerance) RETURN
213 40 : status = 0
214 50 : END SUBROUTINE inversion_representation
215 46 : END MODULE topology_symmetry
|
