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 2 : PROGRAM cell_opt_utils_TEST
9 2 : USE cell_methods, ONLY: cell_create
10 : USE cell_opt_utils, ONLY: &
11 : apply_cell_constraints, cell_opt_constraint_label, comp_ax, comp_bx, comp_by, comp_cz, &
12 : keep_angles_conflict, keep_symmetry_conflict, mask_all, mask_b, mask_c, &
13 : parse_cell_constraint, parse_token, pres_crit_shear_pinned
14 : USE cell_types, ONLY: &
15 : cell_release, cell_sym_cubic, cell_sym_hexagonal_gamma_120, cell_sym_hexagonal_gamma_60, &
16 : cell_sym_monoclinic, cell_sym_monoclinic_gamma_ab, cell_sym_none, cell_sym_orthorhombic, &
17 : cell_sym_rhombohedral, cell_sym_tetragonal_ab, cell_sym_tetragonal_ac, &
18 : cell_sym_tetragonal_bc, cell_sym_triclinic, cell_type
19 : USE kinds, ONLY: default_string_length,&
20 : dp
21 : USE string_utilities, ONLY: remove_word
22 :
23 : IMPLICIT NONE
24 :
25 : REAL(KIND=dp), PARAMETER :: initial(6) = [1.0_dp, 2.0_dp, 3.0_dp, 4.0_dp, 5.0_dp, 6.0_dp]
26 : REAL(KIND=dp), DIMENSION(:), POINTER :: gradient => NULL()
27 : TYPE(cell_type), POINTER :: cell => NULL()
28 :
29 2 : CALL cell_create(cell)
30 2 : ALLOCATE (gradient(6))
31 2 : CALL check_parser()
32 2 : CALL check_projections()
33 2 : CALL check_guards()
34 2 : DEALLOCATE (gradient)
35 2 : CALL cell_release(cell)
36 :
37 : CONTAINS
38 :
39 : ! **************************************************************************************************
40 : !> \brief Check token masks, label round trips, rejections and 2D defaults.
41 : ! **************************************************************************************************
42 2 : SUBROUTINE check_parser()
43 : CHARACTER(LEN=4), PARAMETER :: &
44 : legacy(7) = ["NONE", " X", " Y", " Z", " XY", " XZ", " YZ"], names(18) = [ &
45 : CHARACTER(LEN=4) :: "NONE", "X", "Y", "Z", "XY", "XZ", "YZ", "a", "AX", "B", "BX", "BY", &
46 : "C", "CX", "CY", "CZ", "none", " z "]
47 : CHARACTER(LEN=5), PARAMETER :: &
48 : invalid(5) = [CHARACTER(LEN=5) :: "AY", "AZ", "BZ", "ALPHA", "YX"]
49 : INTEGER, PARAMETER :: &
50 : masks(18) = [0, 11, 22, 56, 31, 59, 62, 1, 1, 6, 2, 4, 56, 8, 16, 32, 0, 56]
51 :
52 : CHARACTER(LEN=default_string_length) :: label, reason, tokens(6)
53 : INTEGER :: i, mask, n, wanted
54 : LOGICAL :: accepted
55 :
56 : ! Bits 0..5 represent AX, BX, BY, CX, CY, CZ, respectively.
57 :
58 8 : cell%perd = 1
59 38 : DO i = 1, SIZE(names)
60 72 : CALL parse_cell_constraint([names(i)], cell, .TRUE., mask)
61 36 : IF (mask /= masks(i)) ERROR STOP "Incorrect constraint token mask"
62 38 : IF (i <= SIZE(legacy)) THEN
63 14 : IF (cell_opt_constraint_label(mask) /= legacy(i)) ERROR STOP "Legacy constraint label changed"
64 : END IF
65 : END DO
66 8 : CALL parse_cell_constraint(["AX", "BY", "CZ"], cell, .TRUE., mask)
67 2 : IF (mask /= 37) ERROR STOP "Constraint components were not combined"
68 8 : CALL parse_cell_constraint(["a ", "b ", "ax"], cell, .TRUE., mask)
69 2 : IF (mask /= 7) ERROR STOP "Redundant constraint changed the mask"
70 6 : CALL parse_cell_constraint(["X", "B"], cell, .TRUE., mask)
71 2 : IF (mask /= 15) ERROR STOP "Partially overlapping constraint lost components"
72 :
73 130 : DO wanted = 0, 63
74 128 : label = ADJUSTL(cell_opt_constraint_label(wanted))
75 128 : n = 0
76 478 : DO WHILE (LEN_TRIM(label) > 0)
77 350 : n = n + 1
78 350 : READ (label, *) tokens(n)
79 350 : CALL remove_word(label)
80 478 : label = ADJUSTL(label)
81 : END DO
82 128 : CALL parse_cell_constraint(tokens(:n), cell, .TRUE., mask)
83 130 : IF (mask /= wanted) ERROR STOP "Constraint label round trip failed"
84 : END DO
85 12 : DO i = 1, SIZE(invalid)
86 10 : CALL parse_token(invalid(i), accepted, mask, reason)
87 12 : IF (accepted .OR. LEN_TRIM(reason) == 0) ERROR STOP "Invalid constraint token accepted"
88 : END DO
89 8 : DO i = 1, 3
90 24 : cell%perd = 1
91 6 : cell%perd(i) = 0
92 12 : CALL parse_cell_constraint(["NONE"], cell, .FALSE., mask)
93 6 : IF (mask /= masks(i + 1)) ERROR STOP "Incorrect 2D default constraint"
94 12 : CALL parse_cell_constraint([names(i + 1)], cell, .TRUE., mask)
95 6 : IF (mask /= masks(i + 1)) ERROR STOP "Explicit 2D constraint changed"
96 24 : CALL parse_cell_constraint(["A", "B", "C"], cell, .TRUE., mask)
97 8 : IF (mask /= mask_all) ERROR STOP "Explicit 2D constraint superset changed"
98 : END DO
99 2 : END SUBROUTINE check_parser
100 :
101 : ! **************************************************************************************************
102 : !> \brief Check each projection against a known gradient and apply it twice.
103 : ! **************************************************************************************************
104 2 : SUBROUTINE check_projections()
105 : REAL(KIND=dp) :: a, normal(3), reference(6), root3
106 :
107 2 : root3 = SQRT(3.0_dp)
108 8 : cell%perd = 1
109 2 : cell%orthorhombic = .FALSE.
110 26 : cell%hmat = RESHAPE([2.0_dp, 0.0_dp, 0.0_dp, 1.0_dp, 2.0_dp, 0.0_dp, 1.0_dp, 2.0_dp, 2.0_dp], [3, 3])
111 2 : CALL check_projection(cell_sym_none, initial, .FALSE., 0)
112 2 : CALL check_projection(cell_sym_triclinic, initial, .FALSE., 0)
113 : CALL check_projection(cell_sym_triclinic, &
114 2 : [1.0_dp, 8.0_dp/5, 16.0_dp/5, 26.0_dp/9, 52.0_dp/9, 52.0_dp/9], .TRUE., 0)
115 : CALL check_projection(cell_sym_triclinic, [0.0_dp, 2.0_dp, 0.0_dp, 4.0_dp, 5.0_dp, 0.0_dp], &
116 2 : .FALSE., comp_ax + comp_by + comp_cz)
117 : CALL check_projection(cell_sym_triclinic, [0.0_dp, 8.0_dp/5, 16.0_dp/5, 0.0_dp, 0.0_dp, 0.0_dp], &
118 2 : .TRUE., comp_ax + mask_c)
119 26 : cell%hmat = RESHAPE([2.0_dp, 0.0_dp, 0.0_dp, 0.0_dp, 2.0_dp, 0.0_dp, 0.0_dp, 0.0_dp, 2.0_dp], [3, 3])
120 2 : cell%orthorhombic = .TRUE.
121 2 : CALL check_projection(cell_sym_cubic, [10.0_dp/3, 0.0_dp, 10.0_dp/3, 0.0_dp, 0.0_dp, 10.0_dp/3], .FALSE., 0)
122 2 : CALL check_projection(cell_sym_tetragonal_ab, [2.0_dp, 0.0_dp, 2.0_dp, 0.0_dp, 0.0_dp, 6.0_dp], .FALSE., 0)
123 2 : CALL check_projection(cell_sym_tetragonal_ac, [3.5_dp, 0.0_dp, 3.0_dp, 0.0_dp, 0.0_dp, 3.5_dp], .FALSE., 0)
124 2 : CALL check_projection(cell_sym_tetragonal_bc, [1.0_dp, 0.0_dp, 4.5_dp, 0.0_dp, 0.0_dp, 4.5_dp], .FALSE., 0)
125 2 : CALL check_projection(cell_sym_monoclinic, [1.0_dp, 0.0_dp, 3.0_dp, 4.0_dp, 0.0_dp, 6.0_dp], .FALSE., 0)
126 44 : cell%hmat = RESHAPE([1.0_dp, 0.0_dp, 0.0_dp, 0.5_dp, root3/2, 0.0_dp, 0.0_dp, 0.0_dp, 2.0_dp], [3, 3])
127 2 : cell%orthorhombic = .FALSE.
128 2 : a = 1.0_dp + 0.75_dp*root3
129 14 : CALL check_projection(cell_sym_hexagonal_gamma_60, [a, a/2, root3*a/2, 0.0_dp, 0.0_dp, 6.0_dp], .FALSE., 0)
130 2 : a = 0.75_dp*root3
131 2 : cell%hmat(1, 2) = -0.5_dp
132 14 : CALL check_projection(cell_sym_hexagonal_gamma_120, [a, -a/2, root3*a/2, 0.0_dp, 0.0_dp, 6.0_dp], .FALSE., 0)
133 :
134 : cell%hmat = RESHAPE([1.0_dp, 0.0_dp, 0.0_dp, 0.5_dp, root3/2, 0.0_dp, &
135 44 : 0.5_dp, root3/6, SQRT(2.0_dp/3)], [3, 3])
136 : ! Equal lengths and 60-degree angles leave one radial direction.
137 14 : reference = [1.0_dp, 0.5_dp, root3/2, 0.5_dp, root3/6, SQRT(2.0_dp/3)]
138 28 : reference = reference*DOT_PRODUCT(reference, initial)/3.0_dp
139 2 : CALL check_projection(cell_sym_rhombohedral, reference, .FALSE., 0)
140 8 : cell%hmat(:, 3) = [0.0_dp, 0.0_dp, 2.0_dp]
141 : ! Equal a/b lengths require -dAX + cos(gamma)*dBX + sin(gamma)*dBY = 0.
142 8 : normal = [-1.0_dp, 0.5_dp, root3/2]
143 2 : reference = initial
144 16 : reference(1:3) = initial(1:3) - normal*DOT_PRODUCT(normal, initial(1:3))/2.0_dp
145 6 : reference(4:5) = 0.0_dp
146 2 : CALL check_projection(cell_sym_monoclinic_gamma_ab, reference, .FALSE., 0)
147 :
148 26 : cell%hmat = 0.0_dp
149 2 : cell%hmat(1, 1) = 2.0_dp
150 2 : cell%hmat(2, 2) = 3.0_dp
151 2 : cell%hmat(3, 3) = 4.0_dp
152 2 : cell%orthorhombic = .TRUE.
153 2 : CALL check_projection(cell_sym_orthorhombic, [1.0_dp, 0.0_dp, 3.0_dp, 0.0_dp, 0.0_dp, 6.0_dp], .FALSE., 0)
154 2 : CALL check_projection(cell_sym_orthorhombic, [1.0_dp, 0.0_dp, 3.0_dp, 0.0_dp, 0.0_dp, 6.0_dp], .TRUE., 0)
155 2 : CALL check_projection(cell_sym_cubic, [0.0_dp, 0.0_dp, 0.0_dp, 0.0_dp, 0.0_dp, 0.0_dp], .TRUE., mask_all)
156 2 : END SUBROUTINE check_projections
157 :
158 : ! **************************************************************************************************
159 : !> \brief Compare projected and expected gradients, including idempotence.
160 : !> \param symmetry cell symmetry
161 : !> \param expected projected gradient
162 : !> \param keep_angles apply the angle projection instead of the symmetry projection
163 : !> \param mask fixed components
164 : ! **************************************************************************************************
165 34 : SUBROUTINE check_projection(symmetry, expected, keep_angles, mask)
166 : INTEGER, INTENT(IN) :: symmetry
167 : REAL(KIND=dp), INTENT(IN) :: expected(6)
168 : LOGICAL, INTENT(IN) :: keep_angles
169 : INTEGER, INTENT(IN) :: mask
170 :
171 : INTEGER :: repeat
172 :
173 34 : cell%symmetry_id = symmetry
174 238 : gradient = initial
175 102 : DO repeat = 1, 2
176 68 : CALL apply_cell_constraints(gradient, cell, keep_angles,.NOT. keep_angles, mask)
177 510 : IF (.NOT. ALL(ABS(gradient - expected) <= 1.0E-12_dp)) THEN
178 0 : PRINT *, "Incorrect cell projection: symmetry, mask, repetition:", symmetry, mask, repeat
179 0 : ERROR STOP 1
180 : END IF
181 : END DO
182 34 : END SUBROUTINE check_projection
183 :
184 : ! **************************************************************************************************
185 : !> \brief Check representative legal/conflicting constraints and pressure warnings.
186 : ! **************************************************************************************************
187 2 : SUBROUTINE check_guards()
188 : INTEGER, PARAMETER :: &
189 : masks(8) = [0, comp_ax, mask_b, mask_c, comp_bx, comp_by, comp_cz, mask_all], &
190 : partial(12) = [1, 1, 1, 4, 2, 2, 2, 32, 1, 1, 1, 1], systems(12) = [cell_sym_cubic, &
191 : cell_sym_tetragonal_ab, cell_sym_tetragonal_ac, cell_sym_tetragonal_bc, &
192 : cell_sym_hexagonal_gamma_60, cell_sym_hexagonal_gamma_120, cell_sym_monoclinic_gamma_ab, &
193 : cell_sym_rhombohedral, cell_sym_orthorhombic, cell_sym_monoclinic, cell_sym_triclinic, &
194 : cell_sym_none]
195 : LOGICAL, PARAMETER :: angle_conflict(8) = [.FALSE., .FALSE., .FALSE., .FALSE., .TRUE., .TRUE.&
196 : , .TRUE., .FALSE.], &
197 : pressure(8) = [.FALSE., .FALSE., .TRUE., .TRUE., .TRUE., .FALSE., .FALSE., .FALSE.]
198 :
199 : INTEGER :: i
200 :
201 26 : DO i = 1, SIZE(systems)
202 24 : IF (keep_symmetry_conflict(0, systems(i))) ERROR STOP "Unconstrained symmetry rejected"
203 24 : IF (keep_symmetry_conflict(mask_all, systems(i))) ERROR STOP "Fully fixed symmetry rejected"
204 26 : IF (keep_symmetry_conflict(partial(i), systems(i)) .NEQV. (i <= 8)) ERROR STOP "Incorrect symmetry conflict"
205 : END DO
206 2 : IF (keep_symmetry_conflict(comp_ax + comp_by, cell_sym_tetragonal_ab)) ERROR STOP "Complete ab constraint rejected"
207 18 : DO i = 1, SIZE(masks)
208 16 : IF (keep_angles_conflict(masks(i), .FALSE.) .NEQV. angle_conflict(i)) ERROR STOP "Incorrect angle conflict"
209 16 : IF (keep_angles_conflict(masks(i), .TRUE.)) ERROR STOP "Orthorhombic angle constraint rejected"
210 16 : IF (pres_crit_shear_pinned(masks(i), .FALSE., .FALSE.) .NEQV. pressure(i)) ERROR STOP "Incorrect pressure warning"
211 18 : IF (pres_crit_shear_pinned(masks(i), .TRUE., .TRUE.)) ERROR STOP "Orthorhombic pressure warning"
212 : END DO
213 2 : IF (.NOT. pres_crit_shear_pinned(comp_ax, .FALSE., .TRUE.)) ERROR STOP "Missing pinned-angle pressure warning"
214 2 : IF (pres_crit_shear_pinned(37 + comp_bx, .FALSE., .TRUE.)) ERROR STOP "Pressure warning with all diagonals fixed"
215 2 : END SUBROUTINE check_guards
216 :
217 : END PROGRAM cell_opt_utils_TEST
|
