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 streaming mesh analysis with physical endpoint sewing.
10 : ! **************************************************************************************************
11 : MODULE topology_phason_analysis
12 : USE cp_files, ONLY: close_file,&
13 : open_file
14 : USE ieee_arithmetic, ONLY: ieee_is_finite
15 : USE iso_fortran_env, ONLY: int64
16 : USE kinds, ONLY: default_path_length,&
17 : dp
18 : USE topology_curvature, ONLY: curvature_density,&
19 : link_plaquette,&
20 : polar_link
21 : USE topology_snapshot, ONLY: check_snapshot,&
22 : check_snapshot_seam,&
23 : read_snapshot,&
24 : snapshot_overlap,&
25 : snapshot_type
26 : #include "./base/base_uses.f90"
27 :
28 : IMPLICIT NONE
29 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'topology_phason_analysis'
30 : PRIVATE
31 : TYPE :: phason_result_type
32 : INTEGER :: dimension = 0, rank = 0, SHAPE(4) = 0, nat = 0, point = 0
33 : REAL(KIND=dp) :: value = 0.0_dp, origin(3) = 0.0_dp, directions(3, 4) = 0.0_dp
34 : CHARACTER(LEN=default_path_length) :: first_file = ''
35 : END TYPE phason_result_type
36 : PUBLIC :: analyze_mesh, phason_result_type, compare_refinements
37 :
38 : CONTAINS
39 :
40 : ! **************************************************************************************************
41 : !> \brief Flatten a zero-based index with the first coordinate varying fastest.
42 : !> \param index Zero-based coordinate along each axis
43 : !> \param shape Number of entries along each axis
44 : !> \return One-based linear index
45 : ! **************************************************************************************************
46 0 : INTEGER FUNCTION flat(index, shape) RESULT(value)
47 : INTEGER, INTENT(IN) :: index(:), SHAPE(:)
48 :
49 : INTEGER :: d, stride
50 :
51 0 : value = 1
52 0 : stride = 1
53 0 : DO d = 1, SIZE(shape)
54 0 : value = value + INDEX(d)*stride
55 0 : stride = stride*SHAPE(d)
56 : END DO
57 0 : END FUNCTION flat
58 :
59 : ! **************************************************************************************************
60 : !> \brief Inverse of flat; first coordinate varies fastest.
61 : !> \param value One-based linear index
62 : !> \param shape Number of entries along each axis
63 : !> \param index Zero-based coordinate along each axis
64 : ! **************************************************************************************************
65 0 : SUBROUTINE unflat(value, shape, index)
66 : INTEGER, INTENT(IN) :: value, SHAPE(:)
67 : INTEGER, INTENT(OUT) :: index(:)
68 :
69 : INTEGER :: d, remainder
70 :
71 0 : remainder = value - 1
72 0 : DO d = 1, SIZE(shape)
73 0 : INDEX(d) = MOD(remainder, SHAPE(d))
74 0 : remainder = remainder/SHAPE(d)
75 : END DO
76 0 : END SUBROUTINE unflat
77 :
78 : ! **************************************************************************************************
79 : !> \brief Build physical links or read reference links; stream plaquettes through an anonymous file.
80 : !> \param path Manifest or reference-link filename
81 : !> \param physical True for physical STATE_EXPORT frames, false for reference links
82 : !> \param report Mesh metadata and unrounded Chern estimate for refinement comparison
83 : !> \param io_unit Unit receiving mesh diagnostics
84 : ! **************************************************************************************************
85 0 : SUBROUTINE analyze_mesh(path, physical, report, io_unit)
86 : CHARACTER(LEN=*), INTENT(IN) :: path
87 : LOGICAL, INTENT(IN) :: physical
88 : TYPE(phason_result_type), INTENT(OUT) :: report
89 : INTEGER, INTENT(IN) :: io_unit
90 : CHARACTER(LEN=256) :: header
91 0 : CHARACTER(LEN=default_path_length), ALLOCATABLE :: files(:)
92 : INTEGER :: input, store, stat, dim, rank, nvertex, nframe, nat, vertex, axis, mu, nu, p, i, j, na, nb, nbase, idx, record_size
93 0 : INTEGER, ALLOCATABLE :: shape(:), extended(:), INDEX(:), next(:), zero(:), points(:), permutation(:, :)
94 0 : LOGICAL, ALLOCATABLE :: checked(:)
95 0 : REAL(KIND=dp), ALLOCATABLE :: offsets(:), directions(:, :)
96 : REAL(KIND=dp) :: origin(3), ka(3), kb(3), tol_metric, tol_gap, tol_sv, phase_limit, minimum, local_min, &
97 : max_metric, min_gap, err, gap, value, contribution, max_phase, local_phase, seam_error, &
98 : valence_max, conduction_min, energy_shift, compensation, y, t
99 0 : COMPLEX(KIND=dp), ALLOCATABLE :: raw(:, :), link(:, :), sewing(:, :), u(:, :, :), pqs(:, :, :)
100 0 : TYPE(snapshot_type), ALLOCATABLE :: a, b, base
101 : LOGICAL :: common_reference, prefix
102 : INTEGER(KIND=int64) :: total
103 :
104 0 : CALL open_file(path, unit_number=input, file_status='OLD', file_action='READ')
105 0 : READ (input, '(A)', IOSTAT=stat) header
106 0 : CPASSERT(stat == 0)
107 0 : IF (physical) THEN
108 0 : IF (TRIM(header) /= 'CP2K_PHASON_MESH 1') THEN
109 0 : CPABORT('Expected CP2K_PHASON_MESH 1')
110 : END IF
111 : ELSE
112 0 : IF (TRIM(header) /= 'CP2K_TOPOLOGY_LINKS 1') THEN
113 0 : CPABORT('Expected CP2K_TOPOLOGY_LINKS 1')
114 : END IF
115 : END IF
116 0 : READ (input, *, IOSTAT=stat) dim, rank
117 0 : CPASSERT(stat == 0)
118 0 : IF ((dim /= 2 .AND. dim /= 4) .OR. rank < 1) THEN
119 0 : CPABORT('Invalid mesh dimension or rank')
120 : END IF
121 0 : ALLOCATE (SHAPE(dim), extended(dim), INDEX(dim), next(dim), zero(dim))
122 0 : READ (input, *, IOSTAT=stat) shape
123 0 : CPASSERT(stat == 0)
124 0 : IF (ANY(shape < 3)) THEN
125 0 : CPABORT('At least three samples per mesh axis are required')
126 : END IF
127 0 : total = 1_int64
128 0 : DO axis = 1, dim
129 0 : IF (total > INT(HUGE(1), int64)/INT(SHAPE(axis) + 1, int64)) THEN
130 0 : CPABORT('Mesh index overflow')
131 : END IF
132 0 : total = total*INT(SHAPE(axis) + 1, int64)
133 : END DO
134 0 : extended(:) = shape + 1
135 0 : nvertex = PRODUCT(shape)
136 0 : nframe = INT(total)
137 0 : report%dimension = dim
138 0 : report%rank = rank
139 0 : report%shape(1:dim) = shape
140 0 : READ (input, *, IOSTAT=stat) tol_metric, tol_gap, tol_sv, phase_limit
141 0 : CPASSERT(stat == 0)
142 0 : IF (.NOT. ALL(ieee_is_finite([tol_metric, tol_gap, tol_sv, phase_limit]))) THEN
143 0 : CPABORT('Nonfinite tolerances')
144 : END IF
145 0 : IF (MIN(tol_metric, tol_gap, tol_sv, phase_limit) <= 0.0_dp .OR. MAX(tol_metric, tol_sv) >= 1.0_dp .OR. &
146 : phase_limit >= ACOS(-1.0_dp)) THEN
147 0 : CPABORT('Invalid mesh tolerances')
148 : END IF
149 0 : ALLOCATE (raw(rank, rank), link(rank, rank), sewing(rank, rank), u(rank, rank, 4), pqs(rank, rank, dim*(dim - 1)/2))
150 : ! Anonymous scratch storage avoids persistent, stale link caches.
151 : CALL open_file('', unit_number=store, file_status='SCRATCH', file_access='STREAM', &
152 0 : file_form='UNFORMATTED', file_action='READWRITE')
153 0 : INQUIRE (IOLENGTH=record_size) link
154 0 : minimum = HUGE(1.0_dp)
155 0 : max_metric = 0.0_dp
156 0 : min_gap = HUGE(1.0_dp)
157 0 : seam_error = 0.0_dp
158 0 : valence_max = -HUGE(1.0_dp)
159 0 : conduction_min = HUGE(1.0_dp)
160 0 : common_reference = .FALSE.
161 0 : prefix = .TRUE.
162 0 : IF (physical) THEN
163 0 : ALLOCATE (a, b, base)
164 0 : ALLOCATE (files(nframe), points(nframe), offsets(nframe), checked(nframe), directions(3, dim))
165 0 : checked(:) = .FALSE.
166 0 : READ (input, *, IOSTAT=stat) origin
167 0 : CPASSERT(stat == 0)
168 0 : IF (.NOT. ALL(ieee_is_finite(origin))) THEN
169 0 : CPABORT('Nonfinite mesh origin')
170 : END IF
171 0 : DO axis = 1, dim
172 0 : READ (input, *, IOSTAT=stat) directions(:, axis)
173 0 : CPASSERT(stat == 0)
174 : END DO
175 0 : IF (.NOT. ALL(ieee_is_finite(directions))) THEN
176 0 : CPABORT('Nonfinite reciprocal periods')
177 : END IF
178 0 : IF (MAXVAL(ABS(directions - ANINT(directions))) > 1.e-10_dp) THEN
179 0 : CPABORT('Noninteger reciprocal periods')
180 : END IF
181 0 : READ (input, *, IOSTAT=stat) nat, common_reference
182 0 : CPASSERT(stat == 0)
183 0 : IF (nat < 1) THEN
184 0 : CPABORT('Invalid mesh atom count')
185 : END IF
186 0 : ALLOCATE (permutation(nat, dim))
187 0 : DO axis = 1, dim
188 0 : READ (input, *, IOSTAT=stat) permutation(:, axis)
189 0 : CPASSERT(stat == 0)
190 0 : IF (ANY(permutation(:, axis) < 1) .OR. ANY(permutation(:, axis) > nat)) THEN
191 0 : CPABORT('Invalid seam permutation')
192 : END IF
193 0 : DO idx = 1, nat
194 0 : IF (COUNT(permutation(:, axis) == idx) /= 1) THEN
195 0 : CPABORT('Non-bijective seam map')
196 : END IF
197 : END DO
198 : END DO
199 0 : DO mu = 1, dim
200 0 : DO nu = mu + 1, dim
201 0 : IF (ANY(permutation(permutation(:, mu), nu) /= permutation(permutation(:, nu), mu))) THEN
202 0 : CPABORT('The declared torus atom permutations do not commute')
203 : END IF
204 : END DO
205 : END DO
206 0 : DO idx = 1, nframe
207 0 : READ (input, *, IOSTAT=stat) files(idx), points(idx), offsets(idx)
208 0 : CPASSERT(stat == 0)
209 0 : IF (points(idx) < 1 .OR. .NOT. ieee_is_finite(offsets(idx))) THEN
210 0 : CPABORT('Invalid frame record')
211 : END IF
212 : END DO
213 0 : report%nat = nat
214 0 : report%origin(:) = origin
215 0 : report%directions(:, 1:dim) = directions
216 0 : report%first_file = files(1)
217 0 : report%point = points(1)
218 : END IF
219 0 : DO vertex = 1, nvertex
220 0 : CALL unflat(vertex, shape, index)
221 0 : IF (physical) THEN
222 0 : na = flat(index, extended)
223 0 : CALL frame(na, a)
224 0 : ka(:) = origin + MATMUL(directions, REAL(index, dp)/REAL(shape, dp))
225 : END IF
226 0 : DO axis = 1, dim
227 0 : IF (physical) THEN
228 0 : next(:) = index
229 0 : next(axis) = next(axis) + 1
230 0 : nb = flat(next, extended)
231 0 : CALL frame(nb, b)
232 0 : kb(:) = origin + MATMUL(directions, REAL(next, dp)/REAL(shape, dp))
233 0 : CALL snapshot_overlap(a, b, ka, kb, raw)
234 0 : IF (next(axis) == SHAPE(axis)) THEN
235 0 : zero(:) = next
236 0 : zero(axis) = 0
237 0 : nbase = flat(zero, extended)
238 0 : CALL frame(nbase, base)
239 0 : CALL check_snapshot_seam(base, b, permutation(:, axis), tol_metric)
240 0 : CALL snapshot_overlap(b, base, kb, kb, sewing)
241 : ! Both endpoint frames must span the same physical subspace, not merely have full rank.
242 0 : link(:, :) = MATMUL(CONJG(TRANSPOSE(sewing)), sewing)
243 0 : DO idx = 1, rank
244 0 : link(idx, idx) = link(idx, idx) - 1.0_dp
245 : END DO
246 0 : err = MAXVAL(ABS(link))
247 0 : seam_error = MAX(seam_error, err)
248 0 : IF (err > tol_metric) THEN
249 0 : CPABORT('Endpoint selected subspaces do not close')
250 : END IF
251 0 : CALL polar_link(sewing, link, local_min, stat, tol_sv)
252 0 : IF (stat /= 0) THEN
253 0 : CPABORT('Invalid endpoint sewing link')
254 : END IF
255 0 : raw(:, :) = MATMUL(raw, link)
256 : END IF
257 : ELSE
258 0 : DO j = 1, rank
259 0 : DO i = 1, rank
260 0 : READ (input, *, IOSTAT=stat) ka(1:2)
261 0 : CPASSERT(stat == 0)
262 0 : raw(i, j) = CMPLX(ka(1), ka(2), dp)
263 : END DO
264 : END DO
265 : END IF
266 0 : CALL polar_link(raw, link, local_min, stat, tol_sv)
267 0 : IF (stat /= 0) THEN
268 0 : CPABORT('Singular/nonfinite mesh link; refine or check the subspace')
269 : END IF
270 0 : minimum = MIN(minimum, local_min)
271 0 : WRITE (store, IOSTAT=stat) link
272 0 : CPASSERT(stat == 0)
273 : END DO
274 : END DO
275 : DO
276 0 : READ (input, '(A)', IOSTAT=stat) header
277 0 : IF (stat < 0) EXIT
278 0 : IF (stat /= 0 .OR. LEN_TRIM(header) /= 0) THEN
279 0 : CPABORT('Unexpected trailing mesh data')
280 : END IF
281 : END DO
282 0 : CALL close_file(input)
283 0 : value = 0.0_dp
284 0 : compensation = 0.0_dp
285 0 : max_phase = 0.0_dp
286 0 : DO vertex = 1, nvertex
287 0 : CALL unflat(vertex, shape, index)
288 0 : p = 0
289 0 : DO mu = 1, dim
290 0 : DO nu = mu + 1, dim
291 0 : p = p + 1
292 0 : CALL read_link(vertex, mu, u(:, :, 1))
293 0 : next(:) = index
294 0 : next(mu) = MOD(next(mu) + 1, SHAPE(mu))
295 0 : CALL read_link(flat(next, shape), nu, u(:, :, 2))
296 0 : next(:) = index
297 0 : next(nu) = MOD(next(nu) + 1, SHAPE(nu))
298 0 : CALL read_link(flat(next, shape), mu, u(:, :, 3))
299 0 : CALL read_link(vertex, nu, u(:, :, 4))
300 0 : CALL link_plaquette(u(:, :, 1), u(:, :, 2), u(:, :, 3), u(:, :, 4), pqs(:, :, p))
301 : END DO
302 : END DO
303 0 : CALL curvature_density(pqs, dim, phase_limit, contribution, local_phase, stat)
304 0 : IF (stat /= 0) THEN
305 0 : CPABORT('Unresolved/nonunitary plaquette; refine the mesh')
306 : END IF
307 0 : max_phase = MAX(max_phase, local_phase)
308 0 : y = contribution - compensation
309 0 : t = value + y
310 0 : compensation = (t - value) - y
311 0 : value = t
312 : END DO
313 0 : CALL close_file(store)
314 0 : report%value = value
315 0 : WRITE (io_unit, '(A,I0)') 'PARAMETER_DIMENSION ', dim
316 0 : WRITE (io_unit, '(A,*(I0,1X))') 'MESH ', shape
317 0 : WRITE (io_unit, '(A,ES26.17)') 'CHERN_RAW ', value
318 0 : WRITE (io_unit, '(A,ES26.17)') 'MINIMUM_LINK_SINGULAR_VALUE ', minimum
319 0 : WRITE (io_unit, '(A,ES26.17)') 'MAXIMUM_PLAQUETTE_PHASE ', max_phase
320 0 : IF (physical) THEN
321 0 : WRITE (io_unit, '(A,ES26.17)') 'MAXIMUM_METRIC_ERROR ', max_metric
322 0 : WRITE (io_unit, '(A,ES26.17)') 'MAXIMUM_SEAM_ERROR ', seam_error
323 0 : WRITE (io_unit, '(A,ES26.17)') 'MINIMUM_DIRECT_GAP_HA ', min_gap
324 0 : IF (common_reference .AND. prefix) THEN
325 0 : WRITE (io_unit, '(A,ES26.17)') 'SAMPLED_INDIRECT_GAP_HA ', conduction_min - valence_max
326 : ELSE
327 0 : WRITE (io_unit, '(A)') 'INDIRECT_GAP_NOT_EVALUATED: no asserted common energy reference or not a lowest-band prefix'
328 : END IF
329 : END IF
330 0 : WRITE (io_unit, '(A)') 'MESH_CONVERGENCE_NOT_ESTABLISHED: repeat on finer meshes; retain unrounded C2 values.'
331 : CONTAINS
332 : ! **************************************************************************************************
333 : !> \brief Read one frame and update validation diagnostics on its first use.
334 : !> \param number One-based frame index in the extended mesh
335 : !> \param state Physical frame with selected coefficients
336 : ! **************************************************************************************************
337 0 : SUBROUTINE frame(number, state)
338 : INTEGER, INTENT(IN) :: number
339 : TYPE(snapshot_type), INTENT(OUT) :: state
340 :
341 : INTEGER :: iband
342 :
343 0 : CALL read_snapshot(TRIM(files(number)), points(number), state)
344 0 : IF (state%rank /= rank .OR. SIZE(state%atoms) /= nat) THEN
345 0 : CPABORT('Frame rank/atom count differs from manifest')
346 : END IF
347 0 : IF (.NOT. checked(number)) THEN
348 : CALL check_snapshot(state, tol_metric, tol_gap, err, gap)
349 0 : max_metric = MAX(max_metric, err)
350 0 : min_gap = MIN(min_gap, gap)
351 0 : checked(number) = .TRUE.
352 0 : DO iband = 1, rank
353 0 : IF (state%bands(iband) /= iband) prefix = .FALSE.
354 : END DO
355 0 : IF (rank >= SIZE(state%energies)) prefix = .FALSE.
356 0 : IF (prefix) THEN
357 0 : energy_shift = offsets(number)
358 0 : valence_max = MAX(valence_max, state%energies(rank) + energy_shift)
359 0 : conduction_min = MIN(conduction_min, state%energies(rank + 1) + energy_shift)
360 : END IF
361 : END IF
362 0 : END SUBROUTINE frame
363 : ! **************************************************************************************************
364 : !> \brief Read a previously validated link from the scratch stream.
365 : !> \param vertex_index One-based mesh vertex
366 : !> \param direction One-based parameter axis
367 : !> \param matrix Unitary selected-state link
368 : ! **************************************************************************************************
369 0 : SUBROUTINE read_link(vertex_index, direction, matrix)
370 : INTEGER, INTENT(IN) :: vertex_index, direction
371 : COMPLEX(KIND=dp), INTENT(OUT) :: matrix(:, :)
372 :
373 : INTEGER(KIND=int64) :: position
374 :
375 0 : position = (INT(vertex_index - 1, int64)*INT(dim, int64) + INT(direction - 1, int64))*INT(record_size, int64) + 1_int64
376 0 : READ (store, POS=position, IOSTAT=stat) matrix
377 0 : CPASSERT(stat == 0)
378 0 : END SUBROUTINE read_link
379 : END SUBROUTINE analyze_mesh
380 :
381 : ! **************************************************************************************************
382 : !> \brief Compare two jointly refined meshes of the same declared family; never round the estimate.
383 : !> \param coarse Coarser mesh report
384 : !> \param fine Finer mesh report
385 : !> \param tolerance Maximum inter-mesh difference and distance to an integer
386 : !> \param io_unit Unit receiving convergence diagnostics
387 : ! **************************************************************************************************
388 0 : SUBROUTINE compare_refinements(coarse, fine, tolerance, io_unit)
389 : TYPE(phason_result_type), INTENT(IN) :: coarse, fine
390 : REAL(KIND=dp), INTENT(IN) :: tolerance
391 : INTEGER, INTENT(IN) :: io_unit
392 0 : TYPE(snapshot_type), ALLOCATABLE :: a, b
393 0 : COMPLEX(KIND=dp), ALLOCATABLE :: m(:, :), metric(:, :)
394 : INTEGER :: dim, i
395 : REAL(KIND=dp) :: difference, residual
396 0 : IF (.NOT. ieee_is_finite(tolerance)) THEN
397 0 : CPABORT('Nonfinite convergence tolerance')
398 : END IF
399 0 : IF (tolerance <= 0.0_dp .OR. tolerance >= 0.5_dp) THEN
400 0 : CPABORT('Invalid convergence tolerance')
401 : END IF
402 0 : dim = coarse%dimension
403 0 : IF (fine%dimension /= dim .OR. fine%rank /= coarse%rank .OR. fine%nat /= coarse%nat) THEN
404 0 : CPABORT('Incompatible refinement reports')
405 : END IF
406 0 : DO i = 1, dim
407 0 : IF (fine%shape(i) < 2*coarse%shape(i) .OR. MOD(fine%shape(i), coarse%shape(i)) /= 0) THEN
408 0 : CPABORT('Refine every mesh axis by an integer factor of at least two')
409 : END IF
410 : END DO
411 0 : IF (MAXVAL(ABS(coarse%origin - fine%origin)) > 1.e-10_dp .OR. &
412 : MAXVAL(ABS(coarse%directions - fine%directions)) > 1.e-10_dp) THEN
413 0 : CPABORT('Refinement changes the parameter origin or reciprocal periods')
414 : END IF
415 0 : IF (coarse%nat > 0) THEN
416 0 : ALLOCATE (a, b)
417 0 : CALL read_snapshot(TRIM(coarse%first_file), coarse%point, a)
418 0 : CALL read_snapshot(TRIM(fine%first_file), fine%point, b)
419 0 : ALLOCATE (m(a%rank, a%rank), metric(a%rank, a%rank))
420 0 : CALL snapshot_overlap(a, b, a%k, a%k, m)
421 0 : metric(:, :) = MATMUL(CONJG(TRANSPOSE(m)), m)
422 0 : DO i = 1, a%rank
423 0 : metric(i, i) = metric(i, i) - 1.0_dp
424 : END DO
425 0 : IF (MAXVAL(ABS(metric)) > 1.e-7_dp) THEN
426 0 : CPABORT('Refinement changes the initial physical subspace')
427 : END IF
428 : END IF
429 0 : difference = ABS(fine%value - coarse%value)
430 0 : residual = ABS(fine%value - ANINT(fine%value))
431 0 : WRITE (io_unit, '(A,ES26.17)') 'REFINEMENT_CHANGE ', difference
432 0 : WRITE (io_unit, '(A,ES26.17)') 'INTEGER_RESIDUAL ', residual
433 0 : IF (difference > tolerance .OR. residual > tolerance) THEN
434 0 : CPABORT('Chern mesh refinement has not converged to the requested tolerance')
435 : END IF
436 0 : WRITE (io_unit, '(A)') 'REFINEMENT_CHECK_PASSED: numerical check, not a certificate of a bulk gap.'
437 0 : END SUBROUTINE compare_refinements
438 0 : END MODULE topology_phason_analysis
|
