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 Reader and physical moving-basis links for version-1 topology snapshots.
10 : ! **************************************************************************************************
11 : MODULE topology_snapshot
12 : USE ai_moments, ONLY: cossin
13 : USE cp_files, ONLY: close_file,&
14 : open_file
15 : USE ieee_arithmetic, ONLY: ieee_is_finite
16 : USE kinds, ONLY: dp
17 : USE orbital_pointers, ONLY: indco,&
18 : init_orbital_pointers
19 : #include "./base/base_uses.f90"
20 :
21 : IMPLICIT NONE
22 : PRIVATE
23 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'topology_snapshot'
24 : REAL(KIND=dp), PARAMETER :: pi = 3.1415926535897932384626433832795_dp, cell_tol = 1.e-10_dp
25 : TYPE :: snapshot_shell_type
26 : INTEGER :: first = 0, count = 0, lmin = 0, lmax = 0
27 : REAL(KIND=dp) :: center(3) = 0.0_dp, radius = 0.0_dp
28 : REAL(KIND=dp), ALLOCATABLE :: exponent(:), radii(:), contraction(:, :)
29 : END TYPE snapshot_shell_type
30 : TYPE :: snapshot_atom_type
31 : INTEGER :: kind = 0, count = 0
32 : REAL(KIND=dp) :: center(3) = 0.0_dp
33 : TYPE(snapshot_shell_type), ALLOCATABLE :: shells(:)
34 : END TYPE snapshot_atom_type
35 : TYPE :: snapshot_type
36 : INTEGER :: nao = 0, rank = 0, nspin = 0, channel = 0
37 : INTEGER :: periodic(3) = 0
38 : REAL(KIND=dp) :: cell(3, 3) = 0.0_dp, inverse(3, 3) = 0.0_dp, k(3) = 0.0_dp
39 : TYPE(snapshot_atom_type), ALLOCATABLE :: atoms(:)
40 : INTEGER, ALLOCATABLE :: bands(:)
41 : REAL(KIND=dp), ALLOCATABLE :: energies(:)
42 : COMPLEX(KIND=dp), ALLOCATABLE :: coefficients(:, :)
43 : END TYPE snapshot_type
44 : PUBLIC :: snapshot_type, read_snapshot, snapshot_overlap, check_snapshot, check_snapshot_seam
45 : CONTAINS
46 :
47 : ! **************************************************************************************************
48 : !> \brief Read a single point, retaining only its states; reject malformed/truncated exports.
49 : !> \param filename Version-1 STATE_EXPORT file
50 : !> \param point One-based point index
51 : !> \param state Selected frame
52 : ! **************************************************************************************************
53 0 : SUBROUTINE read_snapshot(filename, point, state)
54 : CHARACTER(LEN=*), INTENT(IN) :: filename
55 : INTEGER, INTENT(IN) :: point
56 : TYPE(snapshot_type), INTENT(OUT) :: state
57 :
58 : CHARACTER(LEN=256) :: header
59 : INTEGER :: first, i, ia, ic, ios, ip, iset, ix, j, &
60 : nat, nc, nk, np, ns, nt, offset, p, &
61 : powers(3), refpowers(3), unit
62 0 : LOGICAL, ALLOCATABLE :: covered(:)
63 : REAL(KIND=dp) :: det, k(3), pair(2)
64 0 : REAL(KIND=dp), ALLOCATABLE :: energies(:), row(:)
65 :
66 0 : CALL open_file(filename, unit_number=unit, file_status='OLD', file_action='READ')
67 0 : READ (unit, '(A)', IOSTAT=ios) header
68 0 : CPASSERT(ios == 0)
69 0 : IF (TRIM(header) /= 'CP2K_TOPOLOGY_STATE 1') THEN
70 0 : CPABORT('Unsupported topology snapshot version')
71 : END IF
72 0 : READ (unit, *, IOSTAT=ios) nat, state%nao, state%rank, nk, state%nspin, nt, state%channel
73 0 : CPASSERT(ios == 0)
74 0 : IF (MIN(nat, state%nao, state%rank, nk, nt, state%channel) < 1) THEN
75 0 : CPABORT('Invalid snapshot dimensions')
76 : END IF
77 0 : IF (state%rank > nt .OR. point < 1 .OR. point > nk) THEN
78 0 : CPABORT('Invalid snapshot point or rank')
79 : END IF
80 0 : IF (state%nspin /= 1 .AND. state%nspin /= 2) THEN
81 0 : CPABORT('Invalid snapshot spin components')
82 : END IF
83 0 : ALLOCATE (state%bands(state%rank), state%energies(nt), energies(nt), &
84 0 : state%coefficients(state%nao*state%nspin, state%rank), state%atoms(nat))
85 0 : READ (unit, *, IOSTAT=ios) state%bands
86 0 : CPASSERT(ios == 0)
87 0 : IF (ANY(state%bands < 1) .OR. ANY(state%bands > nt)) THEN
88 0 : CPABORT('Invalid snapshot band indices')
89 : END IF
90 0 : DO i = 2, state%rank
91 0 : IF (state%bands(i) <= state%bands(i - 1)) THEN
92 0 : CPABORT('Snapshot bands must be strictly ordered')
93 : END IF
94 : END DO
95 0 : DO j = 1, 3
96 0 : READ (unit, *, IOSTAT=ios) state%cell(:, j)
97 0 : CPASSERT(ios == 0)
98 : END DO
99 0 : IF (.NOT. ALL(ieee_is_finite(state%cell))) THEN
100 0 : CPABORT('Nonfinite snapshot cell')
101 : END IF
102 : ASSOCIATE (a => state%cell, b => state%inverse)
103 0 : b(1, :) = [a(2, 2)*a(3, 3) - a(2, 3)*a(3, 2), a(1, 3)*a(3, 2) - a(1, 2)*a(3, 3), a(1, 2)*a(2, 3) - a(1, 3)*a(2, 2)]
104 0 : b(2, :) = [a(2, 3)*a(3, 1) - a(2, 1)*a(3, 3), a(1, 1)*a(3, 3) - a(1, 3)*a(3, 1), a(1, 3)*a(2, 1) - a(1, 1)*a(2, 3)]
105 0 : b(3, :) = [a(2, 1)*a(3, 2) - a(2, 2)*a(3, 1), a(1, 2)*a(3, 1) - a(1, 1)*a(3, 2), a(1, 1)*a(2, 2) - a(1, 2)*a(2, 1)]
106 0 : det = DOT_PRODUCT(a(:, 1), b(1, :))
107 0 : IF (ABS(det) < cell_tol) THEN
108 0 : CPABORT('Singular snapshot cell')
109 : END IF
110 0 : b(:, :) = b/det
111 : END ASSOCIATE
112 0 : READ (unit, *, IOSTAT=ios) state%periodic
113 0 : CPASSERT(ios == 0)
114 0 : IF (ANY(state%periodic < 0) .OR. ANY(state%periodic > 1)) THEN
115 0 : CPABORT('Invalid snapshot periodicity')
116 : END IF
117 : offset = 0
118 0 : DO ia = 1, nat
119 0 : ASSOCIATE (atom => state%atoms(ia))
120 0 : READ (unit, *, IOSTAT=ios) ix, atom%kind, ns, atom%count
121 0 : CPASSERT(ios == 0)
122 0 : IF (ix /= ia .OR. MIN(atom%kind, ns, atom%count) < 1) THEN
123 0 : CPABORT('Invalid snapshot atom')
124 : END IF
125 0 : READ (unit, *, IOSTAT=ios) atom%center
126 0 : CPASSERT(ios == 0)
127 0 : IF (.NOT. ALL(ieee_is_finite(atom%center))) THEN
128 0 : CPABORT('Nonfinite atom center')
129 : END IF
130 0 : ALLOCATE (atom%shells(ns), covered(atom%count))
131 0 : covered(:) = .FALSE.
132 0 : DO iset = 1, ns
133 0 : ASSOCIATE (s => atom%shells(iset))
134 0 : READ (unit, *, IOSTAT=ios) first, s%count, np, nc, s%lmin, s%radius
135 0 : CPASSERT(ios == 0)
136 0 : IF (MIN(first, s%count, np, nc) < 1 .OR. s%lmin < 0) THEN
137 0 : CPABORT('Invalid Gaussian shell')
138 : END IF
139 0 : IF (.NOT. ieee_is_finite(s%radius) .OR. s%radius <= 0) THEN
140 0 : CPABORT('Invalid shell screening radius')
141 : END IF
142 0 : IF (first + s%count - 1 > atom%count) THEN
143 0 : CPABORT('Shell AO range exceeds atom')
144 : END IF
145 0 : IF (ANY(covered(first:first + s%count - 1))) THEN
146 0 : CPABORT('Overlapping shell AO ranges')
147 : END IF
148 0 : covered(first:first + s%count - 1) = .TRUE.
149 0 : s%lmax = 0
150 0 : DO WHILE ((s%lmax + 1)*(s%lmax + 2)*(s%lmax + 3)/6 < nc)
151 0 : s%lmax = s%lmax + 1
152 : END DO
153 0 : IF ((s%lmax + 1)*(s%lmax + 2)*(s%lmax + 3)/6 /= nc .OR. s%lmin > s%lmax) THEN
154 0 : CPABORT('Invalid Cartesian Gaussian count')
155 : END IF
156 0 : CALL init_orbital_pointers(s%lmax + 1)
157 0 : s%first = offset + first
158 0 : s%center(:) = atom%center
159 0 : ALLOCATE (s%exponent(np), s%radii(np), s%contraction(nc*np, s%count), row(s%count))
160 0 : DO ip = 1, np
161 0 : READ (unit, *, IOSTAT=ios) s%exponent(ip), s%radii(ip)
162 0 : CPASSERT(ios == 0)
163 0 : IF (.NOT. ieee_is_finite(s%exponent(ip)) .OR. .NOT. ieee_is_finite(s%radii(ip))) THEN
164 0 : CPABORT('Nonfinite primitive Gaussian')
165 : END IF
166 0 : IF (s%exponent(ip) <= 0 .OR. s%radii(ip) <= 0) THEN
167 0 : CPABORT('Invalid primitive Gaussian')
168 : END IF
169 0 : DO ic = 1, nc
170 0 : READ (unit, *, IOSTAT=ios) powers, row
171 0 : CPASSERT(ios == 0)
172 0 : refpowers(:) = indco(:, ic)
173 0 : IF (ANY(powers /= refpowers)) THEN
174 0 : CPABORT('Unexpected Cartesian Gaussian ordering')
175 : END IF
176 0 : IF (.NOT. ALL(ieee_is_finite(row))) THEN
177 0 : CPABORT('Nonfinite contraction')
178 : END IF
179 0 : s%contraction((ip - 1)*nc + ic, :) = row
180 : END DO
181 : END DO
182 0 : DEALLOCATE (row)
183 : END ASSOCIATE
184 : END DO
185 0 : IF (.NOT. ALL(covered)) THEN
186 0 : CPABORT('Incomplete atom AO coverage')
187 : END IF
188 0 : DEALLOCATE (covered)
189 0 : offset = offset + atom%count
190 : END ASSOCIATE
191 : END DO
192 0 : IF (offset /= state%nao) THEN
193 0 : CPABORT('Snapshot AO count mismatch')
194 : END IF
195 0 : DO p = 1, nk
196 0 : READ (unit, *, IOSTAT=ios) ix, k
197 0 : CPASSERT(ios == 0)
198 0 : IF (ix /= p .OR. .NOT. ALL(ieee_is_finite(k))) THEN
199 0 : CPABORT('Invalid snapshot k-point')
200 : END IF
201 0 : READ (unit, *, IOSTAT=ios) energies
202 0 : CPASSERT(ios == 0)
203 0 : IF (.NOT. ALL(ieee_is_finite(energies))) THEN
204 0 : CPABORT('Nonfinite snapshot spectrum')
205 : END IF
206 0 : IF (ANY(energies(2:) < energies(:nt - 1))) THEN
207 0 : CPABORT('Unordered snapshot spectrum')
208 : END IF
209 0 : IF (p == point) THEN
210 0 : state%k(:) = k
211 0 : state%energies(:) = energies
212 : END IF
213 0 : DO j = 1, state%rank
214 0 : DO i = 1, state%nao*state%nspin
215 0 : READ (unit, *, IOSTAT=ios) pair
216 0 : CPASSERT(ios == 0)
217 0 : IF (.NOT. ALL(ieee_is_finite(pair))) THEN
218 0 : CPABORT('Nonfinite state coefficient')
219 : END IF
220 0 : IF (p == point) state%coefficients(i, j) = CMPLX(pair(1), pair(2), dp)
221 : END DO
222 : END DO
223 : END DO
224 : DO
225 0 : READ (unit, '(A)', IOSTAT=ios) header
226 0 : IF (ios < 0) EXIT
227 0 : IF (ios /= 0 .OR. LEN_TRIM(header) /= 0) THEN
228 0 : CPABORT('Unexpected trailing snapshot content')
229 : END IF
230 : END DO
231 0 : CALL close_file(unit)
232 0 : END SUBROUTINE read_snapshot
233 :
234 : ! **************************************************************************************************
235 : !> \brief Apply screened Gaussian cross-geometry operator by atom blocks, without a dense AO matrix.
236 : !> \param a Left physical state
237 : !> \param b Right physical state
238 : !> \param ka Unwrapped left fractional k-point
239 : !> \param kb Unwrapped right fractional k-point
240 : !> \param overlap Selected-state overlap including all spinor components
241 : ! **************************************************************************************************
242 8 : SUBROUTINE snapshot_overlap(a, b, ka, kb, overlap)
243 : TYPE(snapshot_type), INTENT(IN) :: a, b
244 : REAL(KIND=dp), INTENT(IN) :: ka(3), kb(3)
245 : COMPLEX(KIND=dp), INTENT(OUT) :: overlap(:, :)
246 :
247 8 : COMPLEX(KIND=dp) :: part(SIZE(overlap, 1), SIZE(overlap, 2))
248 : INTEGER :: ia
249 :
250 24 : IF (ANY(SHAPE(overlap) /= [a%rank, b%rank])) THEN
251 0 : CPABORT('Incorrect overlap output dimensions')
252 : END IF
253 8 : IF (a%nspin /= b%nspin .OR. a%channel /= b%channel .OR. a%rank /= b%rank) THEN
254 0 : CPABORT('Changed spin or rank')
255 : END IF
256 16 : IF (ANY(a%bands /= b%bands)) THEN
257 0 : CPABORT('Changed selected band indices')
258 : END IF
259 8 : IF (SIZE(a%atoms) /= SIZE(b%atoms) .OR. a%nao /= b%nao) THEN
260 0 : CPABORT('Changed atom or AO count')
261 : END IF
262 136 : IF (MAXVAL(ABS(a%cell - b%cell)) > cell_tol .OR. ANY(a%periodic /= b%periodic)) THEN
263 0 : CPABORT('Changed cell')
264 : END IF
265 64 : IF (.NOT. ALL(ieee_is_finite(ka)) .OR. .NOT. ALL(ieee_is_finite(kb))) THEN
266 0 : CPABORT('Nonfinite k-point')
267 : END IF
268 56 : IF (MAXVAL(ABS(ka - a%k - ANINT(ka - a%k))) > cell_tol .OR. &
269 : MAXVAL(ABS(kb - b%k - ANINT(kb - b%k))) > cell_tol) THEN
270 0 : CPABORT('Frame k-point does not match request')
271 : END IF
272 24 : overlap(:, :) = 0.0_dp
273 8 : !$OMP PARALLEL DO DEFAULT(NONE) SHARED(A,b,ka,kb) PRIVATE(ia,part) REDUCTION(+:overlap) SCHEDULE(DYNAMIC)
274 : DO ia = 1, SIZE(a%atoms)
275 : CALL atom_overlap(a, b, ia, ka, kb, part)
276 : overlap(:, :) = overlap + part
277 : END DO
278 : !$OMP END PARALLEL DO
279 8 : END SUBROUTINE snapshot_overlap
280 :
281 : ! **************************************************************************************************
282 : !> \brief Screened atom-row contribution using CP2K's cossin primitive integral recurrence.
283 : !> \param a Left physical state
284 : !> \param b Right physical state
285 : !> \param ia Left atom whose AO rows are contracted
286 : !> \param ka Unwrapped left fractional k-point
287 : !> \param kb Unwrapped right fractional k-point
288 : !> \param atom_link This atom's contribution to the selected-state overlap
289 : ! **************************************************************************************************
290 8 : SUBROUTINE atom_overlap(a, b, ia, ka, kb, atom_link)
291 : TYPE(snapshot_type), INTENT(IN) :: a, b
292 : INTEGER, INTENT(IN) :: ia
293 : REAL(KIND=dp), INTENT(IN) :: ka(3), kb(3)
294 : COMPLEX(KIND=dp), INTENT(OUT) :: atom_link(:, :)
295 :
296 : COMPLEX(KIND=dp) :: phase
297 8 : COMPLEX(KIND=dp), ALLOCATABLE :: block(:, :), oc(:, :, :)
298 : INTEGER :: af, bf, hi(3), i, ib, image(3), j, k, &
299 : lo(3), nc_a, nc_b, s, sa, sb
300 : REAL(KIND=dp) :: bounds(3), cutoff, disp(3), q(3), rb(3)
301 8 : REAL(KIND=dp), ALLOCATABLE :: cosine(:, :), sine(:, :)
302 :
303 64 : q(:) = 2.0_dp*pi*MATMUL(TRANSPOSE(a%inverse), kb - ka)
304 24 : atom_link(:, :) = 0.0_dp
305 16 : DO sa = 1, SIZE(a%atoms(ia)%shells)
306 8 : ASSOCIATE (left => a%atoms(ia)%shells(sa))
307 40 : ALLOCATE (oc(left%count, b%rank, a%nspin))
308 8 : oc(:, :, :) = 0.0_dp
309 8 : nc_a = SIZE(left%contraction, 1)
310 16 : DO ib = 1, SIZE(b%atoms)
311 24 : DO sb = 1, SIZE(b%atoms(ib)%shells)
312 8 : ASSOCIATE (right => b%atoms(ib)%shells(sb))
313 8 : cutoff = left%radius + right%radius
314 128 : disp(:) = MATMUL(a%inverse, right%center - left%center)
315 32 : DO i = 1, 3
316 104 : bounds(i) = cutoff*NORM2(a%inverse(i, :))
317 : END DO
318 32 : lo(:) = CEILING(-disp - bounds)
319 32 : hi(:) = FLOOR(-disp + bounds)
320 80 : WHERE (a%periodic == 0)
321 : lo = 0
322 : hi = 0
323 : END WHERE
324 8 : nc_b = SIZE(right%contraction, 1)
325 72 : ALLOCATE (cosine(nc_a, nc_b), sine(nc_a, nc_b), block(left%count, right%count))
326 8 : block(:, :) = 0.0_dp
327 16 : DO k = lo(3), hi(3)
328 24 : DO j = lo(2), hi(2)
329 24 : DO i = lo(1), hi(1)
330 32 : image(:) = [i, j, k]
331 152 : rb(:) = right%center + MATMUL(a%cell, REAL(image, dp))
332 32 : IF (NORM2(rb - left%center) > cutoff) CYCLE
333 : CALL cossin(left%lmax, SIZE(left%exponent), left%exponent, left%radii, left%lmin, &
334 : right%lmax, SIZE(right%exponent), right%exponent, right%radii, right%lmin, &
335 8 : left%center, rb, q, cosine, sine)
336 32 : phase = EXP(CMPLX(0.0_dp, 2.0_dp*pi*DOT_PRODUCT(kb, REAL(image, dp)), dp))
337 32 : block(:, :) = block + phase*MATMUL(TRANSPOSE(left%contraction), &
338 200 : MATMUL(CMPLX(cosine, -sine, dp), right%contraction))
339 : END DO
340 : END DO
341 : END DO
342 16 : DO s = 1, a%nspin
343 8 : bf = right%first + (s - 1)*b%nao
344 80 : oc(:, :, s) = oc(:, :, s) + MATMUL(block, b%coefficients(bf:bf + right%count - 1, :))
345 : END DO
346 16 : DEALLOCATE (cosine, sine, block)
347 : END ASSOCIATE
348 : END DO
349 : END DO
350 16 : DO s = 1, a%nspin
351 8 : af = left%first + (s - 1)*a%nao
352 80 : atom_link(:, :) = atom_link + MATMUL(CONJG(TRANSPOSE(a%coefficients(af:af + left%count - 1, :))), oc(:, :, s))
353 : END DO
354 16 : DEALLOCATE (oc)
355 : END ASSOCIATE
356 : END DO
357 8 : END SUBROUTINE atom_overlap
358 :
359 : ! **************************************************************************************************
360 : !> \brief Check AO-metric normalization and separation at every selected/excluded boundary.
361 : !> \param a Physical frame to validate
362 : !> \param metric_tol Maximum entrywise error in C^dagger S C minus identity
363 : !> \param gap_tol Minimum sampled selected/excluded-band separation in hartree
364 : !> \param metric_error Observed normalization error
365 : !> \param gap Smallest sampled selected/excluded-band separation in hartree
366 : ! **************************************************************************************************
367 2 : SUBROUTINE check_snapshot(a, metric_tol, gap_tol, metric_error, gap)
368 : TYPE(snapshot_type), INTENT(IN) :: a
369 : REAL(KIND=dp), INTENT(IN) :: metric_tol, gap_tol
370 : REAL(KIND=dp), INTENT(OUT) :: metric_error, gap
371 :
372 : INTEGER :: i
373 4 : LOGICAL :: selected(SIZE(a%energies))
374 4 : COMPLEX(KIND=dp) :: metric(a%rank, a%rank)
375 :
376 2 : CALL snapshot_overlap(a, a, a%k, a%k, metric)
377 4 : DO i = 1, a%rank
378 4 : metric(i, i) = metric(i, i) - 1.0_dp
379 : END DO
380 6 : metric_error = MAXVAL(ABS(metric))
381 2 : IF (.NOT. ieee_is_finite(metric_error) .OR. metric_error > metric_tol) THEN
382 0 : CPABORT('Frame is not AO-metric normalized')
383 : END IF
384 6 : selected(:) = .FALSE.
385 6 : selected(a%bands) = .TRUE.
386 2 : gap = HUGE(1.0_dp)
387 4 : DO i = 1, SIZE(a%energies) - 1
388 4 : IF (selected(i) .NEQV. selected(i + 1)) gap = MIN(gap, a%energies(i + 1) - a%energies(i))
389 : END DO
390 2 : IF (gap == HUGE(1.0_dp) .OR. gap <= gap_tol) THEN
391 0 : CPABORT('No verified sampled subspace isolation')
392 : END IF
393 2 : END SUBROUTINE check_snapshot
394 :
395 : ! **************************************************************************************************
396 : !> \brief Check explicitly prescribed atom permutation, periodic translations and basis at a seam.
397 : !> \param a Start frame
398 : !> \param b Endpoint frame
399 : !> \param permutation One-based endpoint-atom to start-atom map
400 : !> \param tolerance Absolute tolerance for geometry and basis agreement
401 : ! **************************************************************************************************
402 2 : SUBROUTINE check_snapshot_seam(a, b, permutation, tolerance)
403 : TYPE(snapshot_type), INTENT(IN) :: a, b
404 : INTEGER, INTENT(IN) :: permutation(:)
405 : REAL(KIND=dp), INTENT(IN) :: tolerance
406 :
407 : INTEGER :: i, j, n, s
408 : REAL(KIND=dp) :: delta(3), shift(3)
409 :
410 2 : n = SIZE(a%atoms)
411 2 : IF (SIZE(b%atoms) /= n .OR. SIZE(permutation) /= n) THEN
412 0 : CPABORT('Seam atom count mismatch')
413 : END IF
414 8 : IF (ANY(permutation < 1) .OR. ANY(permutation > n)) THEN
415 0 : CPABORT('Invalid seam permutation')
416 : END IF
417 4 : DO i = 1, n
418 6 : IF (COUNT(permutation == i) /= 1) THEN
419 0 : CPABORT('Seam atom map is not a permutation')
420 : END IF
421 : END DO
422 34 : IF (MAXVAL(ABS(a%cell - b%cell)) > tolerance .OR. ANY(a%periodic /= b%periodic)) THEN
423 0 : CPABORT('Seam cell mismatch')
424 : END IF
425 4 : DO i = 1, n
426 2 : j = permutation(i)
427 2 : ASSOCIATE (left => a%atoms(j), right => b%atoms(i))
428 2 : IF (left%kind /= right%kind .OR. left%count /= right%count) THEN
429 0 : CPABORT('Seam changes atom kind or AO count')
430 : END IF
431 2 : IF (SIZE(left%shells) /= SIZE(right%shells)) THEN
432 0 : CPABORT('Seam changes basis')
433 : END IF
434 8 : delta(:) = right%center - left%center
435 32 : shift(:) = ANINT(MATMUL(a%inverse, delta))
436 8 : WHERE (a%periodic == 0) shift = 0.0_dp
437 32 : IF (MAXVAL(ABS(delta - MATMUL(a%cell, shift))) > tolerance) THEN
438 0 : CPABORT('Geometry does not close at seam')
439 : END IF
440 6 : DO s = 1, SIZE(left%shells)
441 2 : ASSOCIATE (x => left%shells(s), y => right%shells(s))
442 2 : IF (x%count /= y%count .OR. x%lmin /= y%lmin .OR. x%lmax /= y%lmax) THEN
443 0 : CPABORT('Seam changes shell')
444 : END IF
445 6 : IF (ANY(SHAPE(x%contraction) /= SHAPE(y%contraction))) THEN
446 0 : CPABORT('Seam changes contractions')
447 : END IF
448 2 : IF (SIZE(x%exponent) /= SIZE(y%exponent)) THEN
449 0 : CPABORT('Seam changes primitives')
450 : END IF
451 4 : IF (MAXVAL(ABS(x%radii - y%radii)) > tolerance .OR. ABS(x%radius - y%radius) > tolerance) THEN
452 0 : CPABORT('Seam changes Gaussian screening radii')
453 : END IF
454 8 : IF (MAXVAL(ABS(x%exponent - y%exponent)) > tolerance .OR. &
455 2 : MAXVAL(ABS(x%contraction - y%contraction)) > tolerance) THEN
456 0 : CPABORT('Seam changes physical basis')
457 : END IF
458 : END ASSOCIATE
459 : END DO
460 : END ASSOCIATE
461 : END DO
462 2 : END SUBROUTINE check_snapshot_seam
463 24 : END MODULE topology_snapshot
|
