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 Unit test for the CELL_OPT constraint parser and guard predicates
10 : !> \par Testing principle
11 : !> Every check computes its expectation with a method that differs from
12 : !> the implementation: the checks enumerate bits one by one, while the
13 : !> implementation uses IAND set intersections. A bug on either side
14 : !> breaks the comparison. The projection checks take their expectation
15 : !> from the guard predicates and their behavior from
16 : !> apply_cell_constraints, so the hand-mirrored coupling tables and the
17 : !> executable projections verify each other.
18 : !> \par Coverage
19 : !> All 64 constraint masks, every token of the grammar, every crystal
20 : !> system, and the parse(label(m)) round trip. About 1800 combinations
21 : !> execute in total, plus the guard decisions checked against the real
22 : !> projections: 64 masks on 14 cell shapes with four gradient draws
23 : !> each, every draw passed through apply_cell_constraints twice, and
24 : !> the pressure-criterion predicate over 64 masks, both cell shapes
25 : !> and both KEEP_ANGLES settings.
26 : !> \par Failure behavior
27 : !> The first mismatch prints the offending mask or token and stops with
28 : !> exit code 1. do_regtest runs this program from the UNIT_TESTS
29 : !> registry and reads only the exit code.
30 : ! **************************************************************************************************
31 2 : PROGRAM cell_opt_utils_TEST
32 2 : USE cell_methods, ONLY: cell_create
33 : USE cell_opt_utils, ONLY: &
34 : apply_cell_constraints, cell_opt_constraint_label, comp_ax, comp_bx, comp_by, comp_cx, &
35 : comp_cy, comp_cz, fixes_partially, keep_angles_conflict, keep_symmetry_conflict, &
36 : mask_ab_plane, mask_all, mask_b, mask_c, mask_diag, mask_diag_ab, mask_diag_ac, &
37 : mask_diag_bc, mask_x, mask_xy, mask_xz, mask_y, mask_yz, mask_z, parse_cell_constraint, &
38 : parse_token, pres_crit_shear_pinned, slot_is_fixed
39 : USE cell_types, ONLY: &
40 : cell_release, cell_sym_cubic, cell_sym_hexagonal_gamma_120, cell_sym_hexagonal_gamma_60, &
41 : cell_sym_monoclinic, cell_sym_monoclinic_gamma_ab, cell_sym_none, cell_sym_orthorhombic, &
42 : cell_sym_rhombohedral, cell_sym_tetragonal_ab, cell_sym_tetragonal_ac, &
43 : cell_sym_tetragonal_bc, cell_sym_triclinic, cell_type
44 : USE kinds, ONLY: default_string_length,&
45 : dp
46 :
47 : IMPLICIT NONE
48 :
49 2 : CALL check_parse_tokens()
50 2 : CALL check_parse_2d_completion()
51 2 : CALL check_fixes_partially()
52 2 : CALL check_slot_is_fixed()
53 2 : CALL check_label_parse_roundtrip()
54 2 : CALL check_label_legacy_strings()
55 2 : CALL check_keep_angles_conflict()
56 2 : CALL check_keep_symmetry_conflict()
57 2 : CALL check_parse_token_rejections()
58 2 : CALL check_overlap_classification()
59 2 : CALL check_projection_guard_consistency()
60 2 : CALL check_pres_crit_shear_pinned()
61 :
62 2 : PRINT *, "cell_opt_utils_unittest: all checks passed"
63 : CONTAINS
64 :
65 : ! **************************************************************************************************
66 : !> \brief Parse every token of the grammar and compare against the mask table.
67 : !> Covers unions, free case and spacing, and the redundancy warning path.
68 : ! **************************************************************************************************
69 2 : SUBROUTINE check_parse_tokens()
70 : TYPE(cell_type), POINTER :: cell
71 :
72 2 : CALL make_cell3d(cell)
73 : ! Every call parses with explicit=.TRUE., so the 2D auto-completion
74 : ! stays out of the way and the check sees the raw token semantics.
75 2 : CALL parse_and_check("NONE", 0, cell)
76 2 : CALL parse_and_check("none", 0, cell)
77 2 : CALL parse_and_check("A", comp_ax, cell)
78 2 : CALL parse_and_check("AX", comp_ax, cell)
79 2 : CALL parse_and_check("B", mask_b, cell)
80 2 : CALL parse_and_check("BX", comp_bx, cell)
81 2 : CALL parse_and_check("BY", comp_by, cell)
82 2 : CALL parse_and_check("C", mask_c, cell)
83 2 : CALL parse_and_check("Z", mask_z, cell)
84 2 : CALL parse_and_check("CX", comp_cx, cell)
85 2 : CALL parse_and_check("CY", comp_cy, cell)
86 2 : CALL parse_and_check("CZ", comp_cz, cell)
87 2 : CALL parse_and_check("X", mask_x, cell)
88 2 : CALL parse_and_check("Y", mask_y, cell)
89 2 : CALL parse_and_check("XY", mask_xy, cell)
90 2 : CALL parse_and_check("XZ", mask_xz, cell)
91 2 : CALL parse_and_check("YZ", mask_yz, cell)
92 : ! Union semantics: several tokens accumulate. Spelling, case and
93 : ! spacing are free. The last line repeats tokens on purpose and only
94 : ! expects a warning, not a failure of the parse.
95 2 : CALL parse_and_check("AX BY CZ", IOR(IOR(comp_ax, comp_by), comp_cz), cell)
96 2 : CALL parse_and_check("x y", IOR(mask_x, mask_y), cell)
97 2 : CALL parse_and_check(" z ", mask_z, cell)
98 2 : CALL parse_and_check("a b ax", IOR(comp_ax, mask_b), cell)
99 2 : CALL cell_release(cell)
100 2 : PRINT *, "check_parse_tokens: OK"
101 2 : END SUBROUTINE check_parse_tokens
102 :
103 : ! **************************************************************************************************
104 : !> \brief Parse one token string and compare against the expected mask
105 : !> \param token_string ...
106 : !> \param expected_mask ...
107 : !> \param cell ...
108 : ! **************************************************************************************************
109 42 : SUBROUTINE parse_and_check(token_string, expected_mask, cell)
110 : CHARACTER(LEN=*), INTENT(IN) :: token_string
111 : INTEGER, INTENT(IN) :: expected_mask
112 : TYPE(cell_type), POINTER :: cell
113 :
114 : CHARACTER(LEN=default_string_length) :: tokens(16)
115 : INTEGER :: ntok, parsed_mask
116 :
117 42 : CALL split_tokens(token_string, tokens, ntok)
118 42 : CALL parse_cell_constraint(tokens(:ntok), cell, .TRUE., parsed_mask)
119 42 : IF (parsed_mask /= expected_mask) THEN
120 0 : PRINT *, "check_parse_tokens: unexpected mask for token "//TRIM(token_string)
121 0 : ERROR STOP 1
122 : END IF
123 42 : END SUBROUTINE parse_and_check
124 :
125 : ! **************************************************************************************************
126 : !> \brief Splits a blank-separated string into tokens, mimicking how the input
127 : !> tokenizer delivers the value of a char_t list keyword
128 : !> \param string ...
129 : !> \param tokens ...
130 : !> \param ntok ...
131 : ! **************************************************************************************************
132 340 : SUBROUTINE split_tokens(string, tokens, ntok)
133 : CHARACTER(LEN=*), INTENT(IN) :: string
134 : CHARACTER(LEN=default_string_length), INTENT(OUT) :: tokens(:)
135 : INTEGER, INTENT(OUT) :: ntok
136 :
137 170 : CHARACTER(LEN=LEN(string)) :: rest
138 : INTEGER :: isep
139 :
140 170 : rest = ADJUSTL(string)
141 170 : ntok = 0
142 538 : DO WHILE (LEN_TRIM(rest) > 0)
143 402 : isep = INDEX(rest, " ")
144 402 : ntok = ntok + 1
145 402 : IF (isep == 0) THEN
146 34 : tokens(ntok) = TRIM(rest)
147 : EXIT
148 : END IF
149 368 : tokens(ntok) = rest(:isep - 1)
150 504 : rest = ADJUSTL(rest(isep + 1:))
151 : END DO
152 170 : END SUBROUTINE split_tokens
153 :
154 : ! **************************************************************************************************
155 : !> \brief For 2D periodicity, an unset constraint auto-fixes the non-periodic
156 : !> direction. Explicit constraints that cover the direction stay unchanged.
157 : ! **************************************************************************************************
158 2 : SUBROUTINE check_parse_2d_completion()
159 : INTEGER :: parsed_mask
160 : TYPE(cell_type), POINTER :: cell
161 :
162 2 : CALL make_cell2d(cell)
163 : ! The keyword is unset. The parser must complete the mask so that the
164 : ! non-periodic z direction is fixed, and warn.
165 4 : CALL parse_cell_constraint(["NONE"], cell, .FALSE., parsed_mask)
166 2 : IF (parsed_mask /= mask_z) THEN
167 0 : PRINT *, "check_parse_2d_completion: auto-set did not fix the non-periodic z"
168 0 : ERROR STOP 1
169 : END IF
170 : ! The keyword is set and covers z. The parser must keep the mask as
171 : ! given; YZ covers z through the C vector.
172 4 : CALL parse_cell_constraint(["Z"], cell, .TRUE., parsed_mask)
173 2 : IF (parsed_mask /= mask_z) THEN
174 0 : PRINT *, "check_parse_2d_completion: explicit Z changed"
175 0 : ERROR STOP 1
176 : END IF
177 4 : CALL parse_cell_constraint(["YZ"], cell, .TRUE., parsed_mask)
178 2 : IF (parsed_mask /= mask_yz) THEN
179 0 : PRINT *, "check_parse_2d_completion: explicit YZ changed"
180 0 : ERROR STOP 1
181 : END IF
182 2 : CALL cell_release(cell)
183 2 : PRINT *, "check_parse_2d_completion: OK"
184 2 : END SUBROUTINE check_parse_2d_completion
185 :
186 : ! **************************************************************************************************
187 : !> \brief Compare fixes_partially against a per-bit enumeration for all 64 masks
188 : !> and the eight coupled slot sets (512 combinations).
189 : !> \par Method
190 : !> For each mask and slot set, walk the six bit positions of the set.
191 : !> Track two flags: does the mask fix at least one bit of the set, and
192 : !> does it leave at least one bit of the set free. Partial fixation
193 : !> means both flags are true. This enumeration shares no code with the
194 : !> IAND form inside fixes_partially.
195 : ! **************************************************************************************************
196 2 : SUBROUTINE check_fixes_partially()
197 : INTEGER, PARAMETER :: slot_sets(8) = [mask_b, mask_c, mask_diag, mask_ab_plane, mask_diag_ab,&
198 : mask_diag_ac, mask_diag_bc, mask_all]
199 :
200 : INTEGER :: bit_index, coupled_set, mask_value, &
201 : set_index
202 : LOGICAL :: expected_partial, set_has_fixed_bit, &
203 : set_has_free_bit
204 :
205 18 : DO set_index = 1, SIZE(slot_sets)
206 16 : coupled_set = slot_sets(set_index)
207 1042 : DO mask_value = 0, 63
208 : ! Independent expectation: enumerate the bits of the set.
209 : set_has_fixed_bit = .FALSE.
210 : set_has_free_bit = .FALSE.
211 7168 : DO bit_index = 0, 5
212 7168 : IF (BTEST(coupled_set, bit_index)) THEN
213 2944 : IF (BTEST(mask_value, bit_index)) THEN
214 : set_has_fixed_bit = .TRUE.
215 : ELSE
216 1472 : set_has_free_bit = .TRUE.
217 : END IF
218 : END IF
219 : END DO
220 1024 : expected_partial = set_has_fixed_bit .AND. set_has_free_bit
221 1040 : IF (fixes_partially(mask_value, coupled_set) .NEQV. expected_partial) THEN
222 0 : PRINT *, "check_fixes_partially: predicate disagrees with the bit enumeration"
223 0 : ERROR STOP 1
224 : END IF
225 : END DO
226 : END DO
227 2 : PRINT *, "check_fixes_partially: OK (64 masks x 8 slot sets)"
228 2 : END SUBROUTINE check_fixes_partially
229 :
230 : ! **************************************************************************************************
231 : !> \brief Compare slot_is_fixed against a modulo-power decomposition for all 64
232 : !> masks and the six slots (384 combinations).
233 : !> \par Method
234 : !> Bit (slot - 1) of the mask states whether slot is fixed. Extract that
235 : !> bit arithmetically: the remainder of the mask modulo 2**slot drops
236 : !> all bits at or above slot, so comparing it against 2**(slot - 1)
237 : !> reads exactly bit (slot - 1). BTEST is never used on this side.
238 : ! **************************************************************************************************
239 2 : SUBROUTINE check_slot_is_fixed()
240 : INTEGER :: mask_value, slot_index
241 :
242 130 : DO mask_value = 0, 63
243 898 : DO slot_index = 1, 6
244 768 : IF (slot_is_fixed(mask_value, slot_index) .NEQV. &
245 128 : MOD(mask_value, 2**slot_index) >= 2**(slot_index - 1)) THEN
246 0 : PRINT *, "check_slot_is_fixed: disagrees with the arithmetic decomposition"
247 0 : ERROR STOP 1
248 : END IF
249 : END DO
250 : END DO
251 2 : PRINT *, "check_slot_is_fixed: OK (64 masks x 6 slots)"
252 2 : END SUBROUTINE check_slot_is_fixed
253 :
254 : ! **************************************************************************************************
255 : !> \brief For every mask m, parsing the label of m must reproduce m. This ties the
256 : !> spelling of the grammar to its semantics in both directions.
257 : ! **************************************************************************************************
258 2 : SUBROUTINE check_label_parse_roundtrip()
259 : CHARACTER(LEN=default_string_length) :: tokens(16)
260 : INTEGER :: mask_value, ntok, parsed_mask
261 : TYPE(cell_type), POINTER :: cell
262 :
263 2 : CALL make_cell3d(cell)
264 : ! The label function maps a mask to a token string, the parser maps a
265 : ! token string back to a mask. Requiring the composition to be the
266 : ! identity for all 64 masks proves both directions at once.
267 130 : DO mask_value = 0, 63
268 128 : CALL split_tokens(cell_opt_constraint_label(mask_value), tokens, ntok)
269 128 : CALL parse_cell_constraint(tokens(:ntok), cell, .TRUE., parsed_mask)
270 130 : IF (parsed_mask /= mask_value) THEN
271 0 : PRINT *, "check_label_parse_roundtrip: parse(label(m)) /= m"
272 0 : ERROR STOP 1
273 : END IF
274 : END DO
275 2 : CALL cell_release(cell)
276 2 : PRINT *, "check_label_parse_roundtrip: OK (64 masks)"
277 2 : END SUBROUTINE check_label_parse_roundtrip
278 :
279 : ! **************************************************************************************************
280 : !> \brief The seven legacy constraint labels and NONE must not change. Downstream
281 : !> tools read these strings.
282 : ! **************************************************************************************************
283 2 : SUBROUTINE check_label_legacy_strings()
284 : ! The labels carry leading blanks for column alignment in the banner,
285 : ! so the comparison strips them with ADJUSTL before TRIM.
286 2 : IF (TRIM(ADJUSTL(cell_opt_constraint_label(0))) /= "NONE") THEN
287 0 : PRINT *, "check_label_legacy_strings: NONE label changed"
288 2 : ERROR STOP 1
289 : END IF
290 2 : IF (TRIM(ADJUSTL(cell_opt_constraint_label(mask_x))) /= "X") THEN
291 0 : PRINT *, "check_label_legacy_strings: X label changed"
292 2 : ERROR STOP 1
293 : END IF
294 2 : IF (TRIM(ADJUSTL(cell_opt_constraint_label(mask_y))) /= "Y") THEN
295 0 : PRINT *, "check_label_legacy_strings: Y label changed"
296 2 : ERROR STOP 1
297 : END IF
298 2 : IF (TRIM(ADJUSTL(cell_opt_constraint_label(mask_z))) /= "Z") THEN
299 0 : PRINT *, "check_label_legacy_strings: Z label changed"
300 2 : ERROR STOP 1
301 : END IF
302 2 : IF (TRIM(ADJUSTL(cell_opt_constraint_label(mask_xy))) /= "XY") THEN
303 0 : PRINT *, "check_label_legacy_strings: XY label changed"
304 2 : ERROR STOP 1
305 : END IF
306 2 : IF (TRIM(ADJUSTL(cell_opt_constraint_label(mask_xz))) /= "XZ") THEN
307 0 : PRINT *, "check_label_legacy_strings: XZ label changed"
308 2 : ERROR STOP 1
309 : END IF
310 2 : IF (TRIM(ADJUSTL(cell_opt_constraint_label(mask_yz))) /= "YZ") THEN
311 0 : PRINT *, "check_label_legacy_strings: YZ label changed"
312 2 : ERROR STOP 1
313 : END IF
314 2 : PRINT *, "check_label_legacy_strings: OK"
315 2 : END SUBROUTINE check_label_legacy_strings
316 :
317 : ! **************************************************************************************************
318 : !> \brief Compare keep_angles_conflict against a per-bit enumeration for all 64
319 : !> masks on orthorhombic and non-orthorhombic cells (128 combinations).
320 : !> \par Method
321 : !> An orthorhombic cell can never conflict. On a non-orthorhombic cell,
322 : !> the conflict is a partial fixation of the b vector group {BX,BY} or
323 : !> of the c vector group {CX,CY,CZ}. The expectation enumerates the
324 : !> bits of each group independently of the IAND form.
325 : ! **************************************************************************************************
326 2 : SUBROUTINE check_keep_angles_conflict()
327 : INTEGER :: bit_index, coupled_set, group_index, &
328 : mask_value, shape_index
329 : LOGICAL :: expected_conflict, is_orthorhombic, &
330 : set_has_fixed_bit, set_has_free_bit
331 :
332 6 : DO shape_index = 0, 1
333 4 : is_orthorhombic = (shape_index == 1)
334 262 : DO mask_value = 0, 63
335 256 : expected_conflict = .FALSE.
336 256 : IF (.NOT. is_orthorhombic) THEN
337 : ! Conflict = partial fixation of the b group or the c group.
338 384 : DO group_index = 1, 2
339 256 : IF (group_index == 1) THEN
340 : coupled_set = mask_b
341 : ELSE
342 128 : coupled_set = mask_c
343 : END IF
344 256 : set_has_fixed_bit = .FALSE.
345 256 : set_has_free_bit = .FALSE.
346 1792 : DO bit_index = 0, 5
347 1792 : IF (BTEST(coupled_set, bit_index)) THEN
348 640 : IF (BTEST(mask_value, bit_index)) THEN
349 : set_has_fixed_bit = .TRUE.
350 : ELSE
351 320 : set_has_free_bit = .TRUE.
352 : END IF
353 : END IF
354 : END DO
355 384 : IF (set_has_fixed_bit .AND. set_has_free_bit) expected_conflict = .TRUE.
356 : END DO
357 : END IF
358 260 : IF (keep_angles_conflict(mask_value, is_orthorhombic) .NEQV. expected_conflict) THEN
359 0 : PRINT *, "check_keep_angles_conflict: mismatch, mask", mask_value, &
360 0 : "orthorhombic", is_orthorhombic
361 0 : ERROR STOP 1
362 : END IF
363 : END DO
364 : END DO
365 2 : PRINT *, "check_keep_angles_conflict: OK (64 masks x 2 cell shapes)"
366 2 : END SUBROUTINE check_keep_angles_conflict
367 :
368 : ! **************************************************************************************************
369 : !> \brief Compare keep_symmetry_conflict against a per-bit enumeration for all 64
370 : !> masks and all twelve crystal systems (768 combinations).
371 : !> \par Method
372 : !> Each crystal system couples a set of cell components; the coupled set
373 : !> per system mirrors the table in the KEEP_SYMMETRY guard comment. The
374 : !> expectation is a partial fixation of that set, computed by bit
375 : !> enumeration. Systems without coupling map to an empty set and can
376 : !> never conflict.
377 : ! **************************************************************************************************
378 2 : SUBROUTINE check_keep_symmetry_conflict()
379 : INTEGER, PARAMETER :: crystal_system_ids(12) = [cell_sym_none, cell_sym_triclinic, &
380 : cell_sym_monoclinic, cell_sym_monoclinic_gamma_ab, cell_sym_orthorhombic, &
381 : cell_sym_tetragonal_ab, cell_sym_tetragonal_ac, cell_sym_tetragonal_bc, &
382 : cell_sym_rhombohedral, cell_sym_hexagonal_gamma_60, cell_sym_hexagonal_gamma_120, &
383 : cell_sym_cubic]
384 :
385 : INTEGER :: bit_index, coupled_set, mask_value, &
386 : system_index
387 : LOGICAL :: expected_conflict, set_has_fixed_bit, &
388 : set_has_free_bit
389 :
390 26 : DO system_index = 1, SIZE(crystal_system_ids)
391 : ! Mirror of the coupled-set table in the KEEP_SYMMETRY guard.
392 24 : SELECT CASE (crystal_system_ids(system_index))
393 : CASE (cell_sym_cubic)
394 2 : coupled_set = mask_diag
395 : CASE (cell_sym_rhombohedral)
396 2 : coupled_set = mask_all
397 : CASE (cell_sym_hexagonal_gamma_60, cell_sym_hexagonal_gamma_120, &
398 : cell_sym_monoclinic_gamma_ab)
399 6 : coupled_set = mask_ab_plane
400 : CASE (cell_sym_tetragonal_ab)
401 2 : coupled_set = mask_diag_ab
402 : CASE (cell_sym_tetragonal_ac)
403 2 : coupled_set = mask_diag_ac
404 : CASE (cell_sym_tetragonal_bc)
405 2 : coupled_set = mask_diag_bc
406 : CASE DEFAULT
407 : ! No coupling: triclinic, monoclinic, orthorhombic, none.
408 24 : coupled_set = 0
409 : END SELECT
410 1562 : DO mask_value = 0, 63
411 : set_has_fixed_bit = .FALSE.
412 : set_has_free_bit = .FALSE.
413 10752 : DO bit_index = 0, 5
414 10752 : IF (BTEST(coupled_set, bit_index)) THEN
415 3072 : IF (BTEST(mask_value, bit_index)) THEN
416 : set_has_fixed_bit = .TRUE.
417 : ELSE
418 1536 : set_has_free_bit = .TRUE.
419 : END IF
420 : END IF
421 : END DO
422 1536 : expected_conflict = set_has_fixed_bit .AND. set_has_free_bit
423 1536 : IF (keep_symmetry_conflict(mask_value, crystal_system_ids(system_index)) .NEQV. &
424 24 : expected_conflict) THEN
425 0 : PRINT *, "check_keep_symmetry_conflict: mismatch, mask", mask_value, &
426 0 : "system", crystal_system_ids(system_index)
427 0 : ERROR STOP 1
428 : END IF
429 : END DO
430 : END DO
431 2 : PRINT *, "check_keep_symmetry_conflict: OK (64 masks x 12 systems)"
432 2 : END SUBROUTINE check_keep_symmetry_conflict
433 :
434 : ! **************************************************************************************************
435 : !> \brief Classify the rejecting tokens without triggering an abort
436 : !> \par Method
437 : !> The token classification lives in the pure parse_token, so the
438 : !> rejection cases are testable in-process: AY, AZ and BZ must be
439 : !> rejected with the always-zero reason, the angle tokens with the
440 : !> not-supported reason, and a near-miss with the unknown-token reason
441 : !> naming the accepted vocabulary.
442 : ! **************************************************************************************************
443 10 : SUBROUTINE check_parse_token_rejections()
444 : CHARACTER(LEN=default_string_length) :: reason
445 : INTEGER :: mask_value
446 : LOGICAL :: accepted
447 :
448 : ! AY, AZ, BZ: the cell convention places a along x and b in the xy
449 : ! plane, so these components are identically zero. The reason says so.
450 2 : CALL parse_token("AY", accepted, mask_value, reason)
451 2 : IF (accepted .OR. INDEX(reason, "always zero") == 0) THEN
452 0 : PRINT *, "check_parse_token_rejections: AY not rejected with the always-zero reason"
453 0 : ERROR STOP 1
454 : END IF
455 2 : CALL parse_token("AZ", accepted, mask_value, reason)
456 2 : IF (accepted .OR. INDEX(reason, "always zero") == 0) THEN
457 0 : PRINT *, "check_parse_token_rejections: AZ not rejected with the always-zero reason"
458 0 : ERROR STOP 1
459 : END IF
460 2 : CALL parse_token("BZ", accepted, mask_value, reason)
461 2 : IF (accepted .OR. INDEX(reason, "always zero") == 0) THEN
462 0 : PRINT *, "check_parse_token_rejections: BZ not rejected with the always-zero reason"
463 0 : ERROR STOP 1
464 : END IF
465 : ! The angle tokens name what is missing rather than pretending support.
466 2 : CALL parse_token("ALPHA", accepted, mask_value, reason)
467 2 : IF (accepted .OR. INDEX(reason, "not supported") == 0) THEN
468 0 : PRINT *, "check_parse_token_rejections: ALPHA not rejected as unsupported"
469 0 : ERROR STOP 1
470 : END IF
471 : ! A near-miss gets the unknown-token reason with the accepted vocabulary.
472 2 : CALL parse_token("YX", accepted, mask_value, reason)
473 2 : IF (accepted .OR. INDEX(reason, "Unknown CELL_OPT CONSTRAINT token") == 0 .OR. &
474 : INDEX(reason, "Accepted tokens") == 0) THEN
475 0 : PRINT *, "check_parse_token_rejections: YX not rejected with the vocabulary"
476 0 : ERROR STOP 1
477 : END IF
478 2 : PRINT *, "check_parse_token_rejections: OK (AY AZ BZ ALPHA YX)"
479 2 : END SUBROUTINE check_parse_token_rejections
480 :
481 : ! **************************************************************************************************
482 : !> \brief The overlap warning branches on full versus partial coverage
483 : !> \par Method
484 : !> The warning classification is a pure decision: overlap equal to the
485 : !> token is a full redundancy; a non-empty proper subset is a partial
486 : !> overlap, and its additions are the complement of the old mask inside
487 : !> the token, which must be non-empty. Checked for all 64 x 63 pairs.
488 : ! **************************************************************************************************
489 2 : SUBROUTINE check_overlap_classification()
490 : INTEGER :: old_mask_value, overlap_value, &
491 : token_mask_value
492 :
493 130 : DO old_mask_value = 0, 63
494 8194 : DO token_mask_value = 1, 63
495 8064 : overlap_value = IAND(token_mask_value, old_mask_value)
496 8192 : IF (overlap_value /= token_mask_value .AND. overlap_value /= 0) THEN
497 : ! Partial overlap: the additions the warning names must be
498 : ! non-empty, else the token would be fully redundant.
499 5404 : IF (IAND(NOT(old_mask_value), token_mask_value) == 0) THEN
500 0 : PRINT *, "check_overlap_classification: partial overlap with empty additions"
501 0 : ERROR STOP 1
502 : END IF
503 : END IF
504 : END DO
505 : END DO
506 2 : PRINT *, "check_overlap_classification: OK (64 x 63 mask pairs)"
507 2 : END SUBROUTINE check_overlap_classification
508 :
509 : ! **************************************************************************************************
510 : !> \brief The guard decisions must agree with the projections they protect
511 : !> \par Method
512 : !> A guard rejects a mask because the projection and the hard zeroing
513 : !> would fight over the coupled slots. The signature of that fight is
514 : !> non-idempotence: with a consistent mask, applying
515 : !> apply_cell_constraints to its own output reproduces it, and the
516 : !> output lies in the projected subspace with every fixed slot exactly
517 : !> zero. With a conflicting mask, the projection pushes components back
518 : !> into slots the zeroing has emptied, so a second application moves
519 : !> the result for a generic gradient. The expectation comes from the
520 : !> guard predicate, the behavior from the projection code: the coupled
521 : !> set table and the SELECT CASE in apply_cell_constraints verify each
522 : !> other instead of a third hand-written copy.
523 : !> \note The monoclinic_gamma_ab projection is idempotent only on a cell
524 : !> with equal a and b lengths, and the KEEP_ANGLES projections divide
525 : !> by the norms of b and c. The shaped cells below satisfy both.
526 : ! **************************************************************************************************
527 2 : SUBROUTINE check_projection_guard_consistency()
528 : INTEGER, PARAMETER :: crystal_system_ids(12) = [cell_sym_none, cell_sym_triclinic, &
529 : cell_sym_monoclinic, cell_sym_monoclinic_gamma_ab, cell_sym_orthorhombic, &
530 : cell_sym_tetragonal_ab, cell_sym_tetragonal_ac, cell_sym_tetragonal_bc, &
531 : cell_sym_rhombohedral, cell_sym_hexagonal_gamma_60, cell_sym_hexagonal_gamma_120, &
532 : cell_sym_cubic]
533 :
534 : INTEGER :: mask_value, shape_index, system_index
535 : LOGICAL :: expected_conflict, is_orthorhombic
536 : TYPE(cell_type), POINTER :: cell
537 :
538 : ! Part one: KEEP_ANGLES against keep_angles_conflict, on an
539 : ! orthorhombic and a triclinic cell. The triclinic cell has nonzero
540 : ! in-plane components in b and c, so a partial fixation of either
541 : ! vector is visible in the projections.
542 6 : DO shape_index = 1, 2
543 4 : is_orthorhombic = (shape_index == 2)
544 4 : IF (is_orthorhombic) THEN
545 2 : CALL make_projection_cell(cell, .TRUE., cell_sym_orthorhombic)
546 : ELSE
547 2 : CALL make_projection_cell(cell, .FALSE., cell_sym_triclinic)
548 : END IF
549 260 : DO mask_value = 0, 63
550 256 : expected_conflict = keep_angles_conflict(mask_value, is_orthorhombic)
551 : CALL projection_case(cell, .TRUE., .FALSE., mask_value, expected_conflict, &
552 260 : "check_projection_guard_consistency: KEEP_ANGLES")
553 : END DO
554 6 : CALL cell_release(cell)
555 : END DO
556 :
557 : ! Part two: KEEP_SYMMETRY against keep_symmetry_conflict, on one
558 : ! shaped cell per crystal system.
559 26 : DO system_index = 1, SIZE(crystal_system_ids)
560 24 : CALL make_projection_cell(cell, .FALSE., crystal_system_ids(system_index))
561 1560 : DO mask_value = 0, 63
562 1536 : expected_conflict = keep_symmetry_conflict(mask_value, crystal_system_ids(system_index))
563 : CALL projection_case(cell, .FALSE., .TRUE., mask_value, expected_conflict, &
564 1560 : "check_projection_guard_consistency: KEEP_SYMMETRY")
565 : END DO
566 26 : CALL cell_release(cell)
567 : END DO
568 2 : PRINT *, "check_projection_guard_consistency: OK (64 masks x 14 cell shapes x 4 draws)"
569 2 : END SUBROUTINE check_projection_guard_consistency
570 :
571 : ! **************************************************************************************************
572 : !> \brief One cell and mask: a stationary projection iff the guard sees no conflict
573 : !> \param cell shaped cell carrying orthorhombic flag and symmetry id
574 : !> \param use_keep_angles run the KEEP_ANGLES projection
575 : !> \param use_keep_symmetry run the KEEP_SYMMETRY projection
576 : !> \param mask_value constraint mask under test
577 : !> \param expect_conflict guard decision for this mask and cell
578 : !> \param label check name prefix for failure messages
579 : ! **************************************************************************************************
580 1792 : SUBROUTINE projection_case(cell, use_keep_angles, use_keep_symmetry, mask_value, &
581 : expect_conflict, label)
582 : TYPE(cell_type), POINTER :: cell
583 : LOGICAL, INTENT(IN) :: use_keep_angles, use_keep_symmetry
584 : INTEGER, INTENT(IN) :: mask_value
585 : LOGICAL, INTENT(IN) :: expect_conflict
586 : CHARACTER(LEN=*), INTENT(IN) :: label
587 :
588 : INTEGER :: draw, slot
589 : LOGICAL :: any_draw_moved
590 : REAL(KIND=dp) :: scale
591 1792 : REAL(KIND=dp), DIMENSION(:), POINTER :: g1, g2
592 :
593 1792 : ALLOCATE (g1(6), g2(6))
594 1792 : any_draw_moved = .FALSE.
595 8960 : DO draw = 1, 4
596 50176 : DO slot = 1, 6
597 : ! Deterministic entries of order one, varied per draw so an
598 : ! accidental symmetry of one gradient cannot mask a moved
599 : ! result. The fractional part keeps entries away from zero.
600 50176 : g1(slot) = MOD(draw*7 + slot*13, 19) - 9 + 0.37_dp*slot + 0.11_dp*draw
601 : END DO
602 7168 : CALL apply_cell_constraints(g1, cell, use_keep_angles, use_keep_symmetry, mask_value)
603 93184 : g2 = g1
604 7168 : CALL apply_cell_constraints(g2, cell, use_keep_angles, use_keep_symmetry, mask_value)
605 50176 : scale = MAX(1.0_dp, MAXVAL(ABS(g1)))
606 : ! The hard zeroing runs last, so fixed slots are exactly zero in
607 : ! the output of the first application, whatever the projections
608 : ! did before it.
609 50176 : DO slot = 1, 6
610 50176 : IF (slot_is_fixed(mask_value, slot) .AND. g1(slot) /= 0.0_dp) THEN
611 0 : PRINT *, TRIM(label)//": fixed slot nonzero, mask", mask_value, "slot", slot
612 0 : ERROR STOP 1
613 : END IF
614 : END DO
615 8960 : IF (expect_conflict) THEN
616 22736 : IF (MAXVAL(ABS(g2 - g1)) > 1.0E-6_dp*scale) any_draw_moved = .TRUE.
617 : ELSE
618 : ! Consistent mask: the second application reproduces the first
619 : ! up to rounding, and the output lies in the projected
620 : ! subspace the projections define.
621 27440 : IF (MAXVAL(ABS(g2 - g1)) > 1.0E-12_dp*scale) THEN
622 0 : PRINT *, TRIM(label)//": projection not stationary, mask", mask_value
623 0 : ERROR STOP 1
624 : END IF
625 3920 : IF (use_keep_angles) CALL check_angles_output(g1, cell, scale, label, mask_value)
626 3920 : IF (use_keep_symmetry) CALL check_symmetry_output(g1, cell, scale, label, mask_value)
627 : END IF
628 : END DO
629 1792 : IF (expect_conflict .AND. .NOT. any_draw_moved) THEN
630 0 : PRINT *, TRIM(label)//": conflict mask never moved the projection, mask", mask_value
631 0 : ERROR STOP 1
632 : END IF
633 1792 : DEALLOCATE (g1, g2)
634 1792 : END SUBROUTINE projection_case
635 :
636 : ! **************************************************************************************************
637 : !> \brief A KEEP_ANGLES output must follow the cell vectors it projects onto
638 : !> \param g projected gradient (six cell components)
639 : !> \param cell shaped cell (hmat and orthorhombic flag are read)
640 : !> \param scale magnitude scale for the comparison tolerances
641 : !> \param label check name prefix for failure messages
642 : !> \param mask_value constraint mask under test
643 : ! **************************************************************************************************
644 576 : SUBROUTINE check_angles_output(g, cell, scale, label, mask_value)
645 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: g
646 : TYPE(cell_type), POINTER :: cell
647 : REAL(KIND=dp), INTENT(IN) :: scale
648 : CHARACTER(LEN=*), INTENT(IN) :: label
649 : INTEGER, INTENT(IN) :: mask_value
650 :
651 576 : IF (cell%orthorhombic) THEN
652 : ! An orthorhombic cell keeps its off-axis slots exactly zero.
653 512 : IF (g(2) /= 0.0_dp .OR. g(4) /= 0.0_dp .OR. g(5) /= 0.0_dp) THEN
654 0 : PRINT *, TRIM(label)//": orthorhombic off-axis slot nonzero, mask", mask_value
655 0 : ERROR STOP 1
656 : END IF
657 : ELSE
658 : ! b lies in the xy plane: the projected pair is parallel to b when
659 : ! the two-dimensional cross product vanishes.
660 64 : IF (ABS(g(2)*cell%hmat(2, 2) - g(3)*cell%hmat(1, 2)) > 1.0E-12_dp*scale) THEN
661 0 : PRINT *, TRIM(label)//": b pair not parallel to b, mask", mask_value
662 0 : ERROR STOP 1
663 : END IF
664 : ! c spans all three directions: every cross-product component
665 : ! with c must vanish.
666 : IF (ABS(g(5)*cell%hmat(3, 3) - g(6)*cell%hmat(2, 3)) > 1.0E-12_dp*scale .OR. &
667 64 : ABS(g(6)*cell%hmat(1, 3) - g(4)*cell%hmat(3, 3)) > 1.0E-12_dp*scale .OR. &
668 : ABS(g(4)*cell%hmat(2, 3) - g(5)*cell%hmat(1, 3)) > 1.0E-12_dp*scale) THEN
669 0 : PRINT *, TRIM(label)//": c triple not parallel to c, mask", mask_value
670 0 : ERROR STOP 1
671 : END IF
672 : END IF
673 576 : END SUBROUTINE check_angles_output
674 :
675 : ! **************************************************************************************************
676 : !> \brief A KEEP_SYMMETRY output must satisfy the coupling of its crystal system
677 : !> \param g projected gradient (six cell components)
678 : !> \param cell shaped cell (symmetry id is read)
679 : !> \param scale magnitude scale for the comparison tolerances
680 : !> \param label check name prefix for failure messages
681 : !> \param mask_value constraint mask under test
682 : !> \note Equalities are asserted where the projection enforces equal slots
683 : !> (cubic, tetragonal) or fixed ratios (hexagonal). Rhombohedral
684 : !> couples all six slots with angle-dependent ratios and is carried
685 : !> by the stationarity check alone. Systems without coupling
686 : !> (none, triclinic) impose nothing.
687 : ! **************************************************************************************************
688 3344 : SUBROUTINE check_symmetry_output(g, cell, scale, label, mask_value)
689 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: g
690 : TYPE(cell_type), POINTER :: cell
691 : REAL(KIND=dp), INTENT(IN) :: scale
692 : CHARACTER(LEN=*), INTENT(IN) :: label
693 : INTEGER, INTENT(IN) :: mask_value
694 :
695 3472 : SELECT CASE (cell%symmetry_id)
696 : CASE (cell_sym_cubic)
697 : IF (ABS(g(1) - g(3)) > 1.0E-12_dp*scale .OR. ABS(g(1) - g(6)) > 1.0E-12_dp*scale .OR. &
698 128 : g(2) /= 0.0_dp .OR. g(4) /= 0.0_dp .OR. g(5) /= 0.0_dp) THEN
699 0 : PRINT *, TRIM(label)//": cubic coupling violated, mask", mask_value
700 0 : ERROR STOP 1
701 : END IF
702 : CASE (cell_sym_tetragonal_ab)
703 : IF (ABS(g(1) - g(3)) > 1.0E-12_dp*scale .OR. g(2) /= 0.0_dp .OR. &
704 256 : g(4) /= 0.0_dp .OR. g(5) /= 0.0_dp) THEN
705 0 : PRINT *, TRIM(label)//": tetragonal_ab coupling violated, mask", mask_value
706 0 : ERROR STOP 1
707 : END IF
708 : CASE (cell_sym_tetragonal_ac)
709 : IF (ABS(g(1) - g(6)) > 1.0E-12_dp*scale .OR. g(2) /= 0.0_dp .OR. &
710 256 : g(4) /= 0.0_dp .OR. g(5) /= 0.0_dp) THEN
711 0 : PRINT *, TRIM(label)//": tetragonal_ac coupling violated, mask", mask_value
712 0 : ERROR STOP 1
713 : END IF
714 : CASE (cell_sym_tetragonal_bc)
715 : IF (ABS(g(3) - g(6)) > 1.0E-12_dp*scale .OR. g(2) /= 0.0_dp .OR. &
716 256 : g(4) /= 0.0_dp .OR. g(5) /= 0.0_dp) THEN
717 0 : PRINT *, TRIM(label)//": tetragonal_bc coupling violated, mask", mask_value
718 0 : ERROR STOP 1
719 : END IF
720 : CASE (cell_sym_orthorhombic)
721 512 : IF (g(2) /= 0.0_dp .OR. g(4) /= 0.0_dp .OR. g(5) /= 0.0_dp) THEN
722 0 : PRINT *, TRIM(label)//": orthorhombic off-axis slots nonzero, mask", mask_value
723 0 : ERROR STOP 1
724 : END IF
725 : CASE (cell_sym_monoclinic)
726 512 : IF (g(2) /= 0.0_dp .OR. g(5) /= 0.0_dp) THEN
727 0 : PRINT *, TRIM(label)//": monoclinic off-axis slots nonzero, mask", mask_value
728 0 : ERROR STOP 1
729 : END IF
730 : CASE (cell_sym_hexagonal_gamma_60)
731 : IF (ABS(g(2) - 0.5_dp*g(1)) > 1.0E-12_dp*scale .OR. &
732 : ABS(g(3) - SQRT(3.0_dp)*g(2)) > 1.0E-12_dp*scale .OR. &
733 128 : g(4) /= 0.0_dp .OR. g(5) /= 0.0_dp) THEN
734 0 : PRINT *, TRIM(label)//": hexagonal_gamma_60 coupling violated, mask", mask_value
735 0 : ERROR STOP 1
736 : END IF
737 : CASE (cell_sym_hexagonal_gamma_120)
738 : IF (ABS(g(2) + 0.5_dp*g(1)) > 1.0E-12_dp*scale .OR. &
739 : ABS(g(3) + SQRT(3.0_dp)*g(2)) > 1.0E-12_dp*scale .OR. &
740 128 : g(4) /= 0.0_dp .OR. g(5) /= 0.0_dp) THEN
741 0 : PRINT *, TRIM(label)//": hexagonal_gamma_120 coupling violated, mask", mask_value
742 0 : ERROR STOP 1
743 : END IF
744 : CASE (cell_sym_monoclinic_gamma_ab)
745 3344 : IF (g(4) /= 0.0_dp .OR. g(5) /= 0.0_dp) THEN
746 0 : PRINT *, TRIM(label)//": monoclinic_gamma_ab off-axis slots nonzero, mask", mask_value
747 0 : ERROR STOP 1
748 : END IF
749 : CASE DEFAULT
750 : ! none, triclinic, rhombohedral: carried by the stationarity
751 : ! check.
752 : END SELECT
753 3344 : END SUBROUTINE check_symmetry_output
754 :
755 : ! **************************************************************************************************
756 : !> \brief The pressure-criterion predicate must fire exactly on the pinned-shear domain
757 : !> \par Method
758 : !> The oracle enumerates slot bits by position: the diagonal slots
759 : !> AX, BY, CZ are bits 0, 2 and 5, the shear slots BX, CX, CY are
760 : !> bits 1, 3 and 4. The criterion is expected to fire on a
761 : !> non-orthorhombic cell when some but not all diagonal slots are
762 : !> fixed and either the mask pins a shear slot or KEEP_ANGLES pins
763 : !> the angles together with a fixed diagonal slot. Zero or all
764 : !> diagonal slots fixed never fire (the hydrostatic trace, resp. the
765 : !> external-pressure bypass). The firing counts per cell shape and
766 : !> KEEP_ANGLES setting are asserted as well, and mask NONE never
767 : !> fires, so existing KEEP_ANGLE regtests without a constraint stay
768 : !> warning-free.
769 : ! **************************************************************************************************
770 2 : SUBROUTINE check_pres_crit_shear_pinned()
771 : INTEGER :: count_fired, ka_index, mask_value, &
772 : shape_index
773 : LOGICAL :: all_diag_fixed, diag_hit, expected, &
774 : is_orthorhombic, shear_hit, &
775 : use_keep_angles
776 :
777 6 : DO shape_index = 0, 1
778 4 : is_orthorhombic = (shape_index == 1)
779 14 : DO ka_index = 0, 1
780 8 : use_keep_angles = (ka_index == 1)
781 8 : count_fired = 0
782 520 : DO mask_value = 0, 63
783 : ! Independent expectation from the bit positions.
784 : all_diag_fixed = BTEST(mask_value, 0) .AND. BTEST(mask_value, 2) .AND. &
785 512 : BTEST(mask_value, 5)
786 : shear_hit = BTEST(mask_value, 1) .OR. BTEST(mask_value, 3) .OR. &
787 512 : BTEST(mask_value, 4)
788 : diag_hit = BTEST(mask_value, 0) .OR. BTEST(mask_value, 2) .OR. &
789 512 : BTEST(mask_value, 5)
790 : expected = .NOT. is_orthorhombic .AND. .NOT. all_diag_fixed .AND. &
791 512 : (shear_hit .OR. (use_keep_angles .AND. diag_hit))
792 : IF (pres_crit_shear_pinned(mask_value, is_orthorhombic, use_keep_angles) &
793 512 : .NEQV. expected) THEN
794 0 : PRINT *, "check_pres_crit_shear_pinned: mismatch, mask", mask_value, &
795 0 : "orthorhombic", is_orthorhombic, "keep_angles", use_keep_angles
796 0 : ERROR STOP 1
797 : END IF
798 520 : IF (expected) count_fired = count_fired + 1
799 : END DO
800 : ! Firing counts of the pinned-shear domain: an orthorhombic
801 : ! cell never fires; on an oblique cell the mask domain has 56
802 : ! members of which the all-diagonals-fixed ones are excluded,
803 : ! and KEEP_ANGLES adds every mask with a fixed diagonal slot.
804 8 : IF (is_orthorhombic) THEN
805 4 : IF (count_fired /= 0) THEN
806 0 : PRINT *, "check_pres_crit_shear_pinned: fired on an orthorhombic cell"
807 0 : ERROR STOP 1
808 : END IF
809 : CYCLE
810 : END IF
811 8 : IF (use_keep_angles) THEN
812 2 : IF (count_fired /= 55) THEN
813 0 : PRINT *, "check_pres_crit_shear_pinned: expected 55 masks, got", count_fired
814 0 : ERROR STOP 1
815 : END IF
816 : ELSE
817 2 : IF (count_fired /= 49) THEN
818 0 : PRINT *, "check_pres_crit_shear_pinned: expected 49 masks, got", count_fired
819 0 : ERROR STOP 1
820 : END IF
821 : END IF
822 : END DO
823 : END DO
824 : ! Effective warned set under KEEP_ANGLES: the merged abort guard
825 : ! (keep_angles_conflict) rejects every mask that fixes part of b or c,
826 : ! so only whole-vector and a-only masks reach the warning; of those,
827 : ! NONE stays outside the domain and ABC takes the bypass.
828 2 : count_fired = 0
829 130 : DO mask_value = 0, 63
830 128 : IF (pres_crit_shear_pinned(mask_value, .FALSE., .TRUE.) .AND. &
831 2 : .NOT. keep_angles_conflict(mask_value, .FALSE.)) THEN
832 12 : count_fired = count_fired + 1
833 42 : IF (.NOT. ANY(mask_value == [1, 6, 7, 56, 57, 62])) THEN
834 0 : PRINT *, "check_pres_crit_shear_pinned: unexpected survivor", mask_value
835 0 : ERROR STOP 1
836 : END IF
837 : END IF
838 : END DO
839 2 : IF (count_fired /= 6) THEN
840 0 : PRINT *, "check_pres_crit_shear_pinned: expected 6 effective masks, got", count_fired
841 0 : ERROR STOP 1
842 : END IF
843 2 : PRINT *, "check_pres_crit_shear_pinned: OK (64 masks x 2 cell shapes x 2 flags, 6 effective)"
844 2 : END SUBROUTINE check_pres_crit_shear_pinned
845 :
846 : ! **************************************************************************************************
847 : !> \brief Create a 3D-periodic cell. parse_cell_constraint reads only the perd field.
848 : !> \param cell ...
849 : ! **************************************************************************************************
850 4 : SUBROUTINE make_cell3d(cell)
851 : TYPE(cell_type), POINTER :: cell
852 :
853 4 : NULLIFY (cell)
854 4 : CALL cell_create(cell)
855 16 : cell%perd(:) = 1
856 4 : END SUBROUTINE make_cell3d
857 :
858 : ! **************************************************************************************************
859 : !> \brief Create a cell periodic in x and y. The non-periodic z direction must be
860 : !> constrained.
861 : !> \param cell ...
862 : ! **************************************************************************************************
863 2 : SUBROUTINE make_cell2d(cell)
864 : TYPE(cell_type), POINTER :: cell
865 :
866 2 : NULLIFY (cell)
867 2 : CALL cell_create(cell)
868 8 : cell%perd(:) = [1, 1, 0]
869 2 : END SUBROUTINE make_cell2d
870 :
871 : ! **************************************************************************************************
872 : !> \brief Create a 3D-periodic cell with a lattice consistent with symmetry_id
873 : !> \param cell ...
874 : !> \param is_orthorhombic sets the orthorhombic flag read by KEEP_ANGLES
875 : !> \param symmetry_id crystal system read by KEEP_SYMMETRY
876 : !> \note apply_cell_constraints reads hmat, the orthorhombic flag and
877 : !> symmetry_id. Each system whose projection reads hmat gets a
878 : !> consistent lattice: monoclinic_gamma_ab needs equal a and b
879 : !> lengths, hexagonal needs the named gamma, rhombohedral needs equal
880 : !> pairwise angles. The other projections are lattice independent and
881 : !> share the generic triclinic lattice.
882 : ! **************************************************************************************************
883 28 : SUBROUTINE make_projection_cell(cell, is_orthorhombic, symmetry_id)
884 : TYPE(cell_type), POINTER :: cell
885 : LOGICAL, INTENT(IN) :: is_orthorhombic
886 : INTEGER, INTENT(IN) :: symmetry_id
887 :
888 28 : NULLIFY (cell)
889 28 : CALL cell_create(cell)
890 112 : cell%perd(:) = 1
891 28 : cell%orthorhombic = is_orthorhombic
892 28 : cell%symmetry_id = symmetry_id
893 4 : SELECT CASE (symmetry_id)
894 : CASE (cell_sym_orthorhombic)
895 : cell%hmat = RESHAPE([10.0_dp, 0.0_dp, 0.0_dp, 0.0_dp, 9.0_dp, 0.0_dp, 0.0_dp, 0.0_dp, &
896 52 : 8.0_dp], [3, 3])
897 : CASE (cell_sym_monoclinic_gamma_ab, cell_sym_hexagonal_gamma_60)
898 : ! a and b of length 2 with gamma = 60 degree, c orthogonal.
899 : cell%hmat = RESHAPE([2.0_dp, 0.0_dp, 0.0_dp, 1.0_dp, 1.7320508075688772_dp, 0.0_dp, &
900 52 : 0.0_dp, 0.0_dp, 8.0_dp], [3, 3])
901 : CASE (cell_sym_hexagonal_gamma_120)
902 : cell%hmat = RESHAPE([2.0_dp, 0.0_dp, 0.0_dp, -1.0_dp, 1.7320508075688772_dp, 0.0_dp, &
903 26 : 0.0_dp, 0.0_dp, 8.0_dp], [3, 3])
904 : CASE (cell_sym_rhombohedral)
905 : ! Three vectors of length 1 with pairwise angles of 60 degree.
906 : cell%hmat = RESHAPE([1.0_dp, 0.0_dp, 0.0_dp, 0.5_dp, 0.8660254037844386_dp, 0.0_dp, &
907 26 : 0.5_dp, 0.2886751345948129_dp, 0.8164965809277260_dp], [3, 3])
908 : CASE DEFAULT
909 : ! Generic triclinic lattice: b and c carry nonzero in-plane
910 : ! components, so partial fixations show up in the projections.
911 : cell%hmat = RESHAPE([10.0_dp, 0.0_dp, 0.0_dp, 1.0_dp, 9.0_dp, 0.0_dp, 0.5_dp, 1.5_dp, &
912 220 : 8.0_dp], [3, 3])
913 : END SELECT
914 28 : END SUBROUTINE make_projection_cell
915 :
916 : END PROGRAM cell_opt_utils_TEST
|
