1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38 package ffx.crystal;
39
40 import static ffx.crystal.LatticeSystem.check;
41 import static org.junit.Assert.assertEquals;
42
43 import java.util.Arrays;
44 import java.util.Collection;
45
46 import ffx.utilities.FFXTest;
47 import org.junit.Test;
48 import org.junit.runner.RunWith;
49 import org.junit.runners.Parameterized;
50 import org.junit.runners.Parameterized.Parameters;
51
52
53
54
55
56
57 @RunWith(Parameterized.class)
58 public class CrystalVolumeDerivativeTest extends FFXTest {
59
60 private final String info;
61 private final Crystal crystal;
62
63 private double eps = 0.00001;
64
65 private double epsD = Math.toDegrees(eps);
66 private double eps2 = 2.0 * eps;
67 private double tolerance = eps * 10.0;
68
69 public CrystalVolumeDerivativeTest(
70 String info,
71 double a,
72 double b,
73 double c,
74 double alpha,
75 double beta,
76 double gamma,
77 String sg) {
78 this.info = info;
79 this.crystal = new Crystal(a, b, c, alpha, beta, gamma, sg);
80 }
81
82 @Parameters
83 public static Collection<Object[]> data() {
84 return Arrays.asList(
85 new Object[][] {
86 {"Triclinic (3TRW)", 28.38, 31.73, 36.75, 90.12, 99.61, 96.52, "P-1"},
87 {"Cubic (3RN4)", 92.69, 92.69, 92.69, 90.0, 90.0, 90.0, "P23"},
88 {"Orthorhombic (2WLD)", 79.09, 94.81, 100.85, 90.0, 90.0, 90.0, "P222"},
89 {"Monoclinic (3V0E)", 50.85, 38.60, 89.83, 90.0, 103.99, 90.0, "P2"},
90 {"Hexagonal (4DAC)", 63.67, 63.67, 40.40, 90.0, 90.0, 120.0, "P6"},
91 {"Tetragonal (3TI8)", 112.56, 112.56, 66.81, 90.0, 90.0, 90.0, "P4"},
92 {"Trigonal (3TNZ)", 108.99, 108.99, 49.40, 90.0, 90.0, 120.0, "P3"}
93 });
94 }
95
96 @Test
97 public void finiteDifferenceTest() {
98
99
100 double a = crystal.a;
101 double b = crystal.b;
102 double c = crystal.c;
103 double alpha = crystal.alpha;
104 double beta = crystal.beta;
105 double gamma = crystal.gamma;
106
107
108 double dVdA = crystal.dVdA;
109 double dVdB = crystal.dVdB;
110 double dVdC = crystal.dVdC;
111 double dVdAlpha = crystal.dVdAlpha;
112 double dVdBeta = crystal.dVdBeta;
113 double dVdGamma = crystal.dVdGamma;
114 double dV;
115
116 switch (crystal.spaceGroup.crystalSystem) {
117 case TRICLINIC:
118 testdVdA();
119 testdVdB();
120 testdVdC();
121 testdVdAlpha();
122 testdVdBeta();
123 testdVdGamma();
124 break;
125 case MONOCLINIC:
126
127 testdVdA();
128 testdVdB();
129 testdVdC();
130 testdVdBeta();
131 break;
132 case ORTHORHOMBIC:
133
134 testdVdA();
135 testdVdB();
136 testdVdC();
137 break;
138 case TETRAGONAL:
139
140 crystal.changeUnitCellParameters(a + eps, b + eps, c, alpha, beta, gamma);
141 dV = crystal.volume;
142 crystal.changeUnitCellParameters(a - eps, b - eps, c, alpha, beta, gamma);
143 dV -= crystal.volume;
144 assertEquals(
145 info + " analytic and FD volume derivatives with respect to A & B should agree.",
146 dVdA + dVdB,
147 dV / eps2,
148 tolerance);
149 assertEquals(
150 info + " analytic volume derivatives with respect to A and B should be equal.",
151 dVdA,
152 dVdB,
153 tolerance);
154 testdVdC();
155 break;
156 case TRIGONAL:
157 if (check(a, b) && check(b, c) && check(alpha, beta) && check(beta, gamma)) {
158
159 crystal.changeUnitCellParameters(a + eps, b + eps, c + eps, alpha, beta, gamma);
160 dV = crystal.volume;
161 crystal.changeUnitCellParameters(a - eps, b - eps, c - eps, alpha, beta, gamma);
162 dV -= crystal.volume;
163 double dTot = dVdA + dVdB + dVdC;
164 assertEquals(
165 info + " analytic and FD volume derivatives with respect to the sides should agree.",
166 dTot,
167 dV / eps2,
168 tolerance);
169 assertEquals(
170 info + " analytic volume derivatives with respect to A and B should be equal.",
171 dVdA,
172 dVdB,
173 tolerance);
174 crystal.changeUnitCellParameters(a, b, c, alpha + epsD, beta + epsD, gamma + epsD);
175 dV = crystal.volume;
176 crystal.changeUnitCellParameters(a, b, c, alpha - epsD, beta - epsD, gamma - epsD);
177 dV -= crystal.volume;
178 dTot = dVdAlpha + dVdBeta + dVdGamma;
179 assertEquals(
180 info + " analytic and FD volume derivatives with respect to the angles should agree.",
181 dTot,
182 dV / eps2,
183 tolerance);
184 assertEquals(
185 info + " analytic volume derivatives with respect to Alpha and Beta should be equal.",
186 dVdAlpha,
187 dVdBeta,
188 tolerance);
189 } else if (check(a, b) && check(alpha, 90.0) && check(beta, 90.0) && check(gamma, 120.0)) {
190
191 crystal.changeUnitCellParameters(a + eps, b + eps, c, alpha, beta, gamma);
192 dV = crystal.volume;
193 crystal.changeUnitCellParameters(a - eps, b - eps, c, alpha, beta, gamma);
194 dV -= crystal.volume;
195 assertEquals(
196 info + " analytic and FD volume derivatives with respect to A & B should agree.",
197 dVdA + dVdB,
198 dV / eps2,
199 tolerance);
200 assertEquals(
201 info + " analytic volume derivatives with respect to A and B should be equal.",
202 dVdA,
203 dVdB,
204 tolerance);
205 testdVdC();
206 }
207 break;
208 case HEXAGONAL:
209
210 crystal.changeUnitCellParameters(a + eps, b + eps, c, alpha, beta, gamma);
211 dV = crystal.volume;
212 crystal.changeUnitCellParameters(a - eps, b - eps, c, alpha, beta, gamma);
213 dV -= crystal.volume;
214 assertEquals(
215 info + " analytic and FD volume derivatives with respect to A & B should agree.",
216 dVdA + dVdB,
217 dV / eps2,
218 tolerance);
219 assertEquals(
220 info + " analytic volume derivatives with respect to A and B should be equal.",
221 dVdA,
222 dVdB,
223 tolerance);
224 testdVdC();
225 break;
226 case CUBIC:
227
228 crystal.changeUnitCellParameters(a + eps, b + eps, c + eps, alpha, beta, gamma);
229 dV = crystal.volume;
230 crystal.changeUnitCellParameters(a - eps, b - eps, c - eps, alpha, beta, gamma);
231 dV -= crystal.volume;
232 double dTot = dVdA + dVdB + dVdC;
233 assertEquals(
234 info + " analytic and FD volume derivatives with respect to A, B & C should agree.",
235 dTot,
236 dV / eps2,
237 tolerance);
238 assertEquals(
239 info + " analytic volume derivatives with respect to A and B should be equal.",
240 dVdA,
241 dVdB,
242 tolerance);
243 break;
244 }
245 }
246
247 private void testdVdA() {
248 double dVdA = crystal.dVdA;
249 double a = crystal.a;
250 double b = crystal.b;
251 double c = crystal.c;
252 double alpha = crystal.alpha;
253 double beta = crystal.beta;
254 double gamma = crystal.gamma;
255 crystal.changeUnitCellParameters(a + eps, b, c, alpha, beta, gamma);
256 double dV = crystal.volume;
257 crystal.changeUnitCellParameters(a - eps, b, c, alpha, beta, gamma);
258 dV -= crystal.volume;
259 assertEquals(
260 info + " analytic and FD volume derivatives with respect to A should agree.",
261 dVdA,
262 dV / eps2,
263 tolerance);
264 }
265
266 private void testdVdB() {
267 double dVdB = crystal.dVdB;
268 double a = crystal.a;
269 double b = crystal.b;
270 double c = crystal.c;
271 double alpha = crystal.alpha;
272 double beta = crystal.beta;
273 double gamma = crystal.gamma;
274 crystal.changeUnitCellParameters(a, b + eps, c, alpha, beta, gamma);
275 double dV = crystal.volume;
276 crystal.changeUnitCellParameters(a, b - eps, c, alpha, beta, gamma);
277 dV -= crystal.volume;
278 assertEquals(
279 info + " analytic and FD volume derivatives with respect to B should agree.",
280 dVdB,
281 dV / eps2,
282 tolerance);
283 }
284
285 private void testdVdC() {
286 double dVdC = crystal.dVdC;
287 double a = crystal.a;
288 double b = crystal.b;
289 double c = crystal.c;
290 double alpha = crystal.alpha;
291 double beta = crystal.beta;
292 double gamma = crystal.gamma;
293 crystal.changeUnitCellParameters(a, b, c + eps, alpha, beta, gamma);
294 double dV = crystal.volume;
295 crystal.changeUnitCellParameters(a, b, c - eps, alpha, beta, gamma);
296 dV -= crystal.volume;
297 assertEquals(
298 info + " analytic and FD volume derivatives with respect to C should agree.",
299 dVdC,
300 dV / eps2,
301 tolerance);
302 }
303
304 private void testdVdAlpha() {
305 double dVdAlpha = crystal.dVdAlpha;
306 double a = crystal.a;
307 double b = crystal.b;
308 double c = crystal.c;
309 double alpha = crystal.alpha;
310 double beta = crystal.beta;
311 double gamma = crystal.gamma;
312 crystal.changeUnitCellParameters(a, b, c, alpha + epsD, beta, gamma);
313 double dV = crystal.volume;
314 crystal.changeUnitCellParameters(a, b, c, alpha - epsD, beta, gamma);
315 dV -= crystal.volume;
316 assertEquals(
317 info + " analytic and FD volume derivatives with respect to Alpha should agree.",
318 dVdAlpha,
319 dV / eps2,
320 tolerance);
321 }
322
323 private void testdVdBeta() {
324 double dVdBeta = crystal.dVdBeta;
325 double a = crystal.a;
326 double b = crystal.b;
327 double c = crystal.c;
328 double alpha = crystal.alpha;
329 double beta = crystal.beta;
330 double gamma = crystal.gamma;
331 crystal.changeUnitCellParameters(a, b, c, alpha, beta + epsD, gamma);
332 double dV = crystal.volume;
333 crystal.changeUnitCellParameters(a, b, c, alpha, beta - epsD, gamma);
334 dV -= crystal.volume;
335 assertEquals(
336 info + " analytic and FD volume derivatives with respect to Beta should agree.",
337 dVdBeta,
338 dV / eps2,
339 tolerance);
340 }
341
342 private void testdVdGamma() {
343 double dVdGamma = crystal.dVdGamma;
344 double a = crystal.a;
345 double b = crystal.b;
346 double c = crystal.c;
347 double alpha = crystal.alpha;
348 double beta = crystal.beta;
349 double gamma = crystal.gamma;
350 crystal.changeUnitCellParameters(a, b, c, alpha, beta, gamma + epsD);
351 double dV = crystal.volume;
352 crystal.changeUnitCellParameters(a, b, c, alpha, beta, gamma - epsD);
353 dV -= crystal.volume;
354 assertEquals(
355 info + " analytic and FD volume derivatives with respect to Gamma should agree.",
356 dVdGamma,
357 dV / eps2,
358 tolerance);
359 }
360 }