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.numerics.math.DoubleMath.dot;
41 import static ffx.numerics.math.DoubleMath.length;
42 import static ffx.numerics.math.MatrixMath.mat3Mat3;
43 import static ffx.numerics.math.MatrixMath.mat3SymVec6;
44 import static ffx.numerics.math.MatrixMath.transpose3;
45 import static ffx.numerics.math.ScalarMath.mod;
46 import static ffx.utilities.Constants.AVOGADRO;
47 import static ffx.utilities.KeywordGroup.UnitCellAndSpaceGroup;
48 import static ffx.utilities.StringUtils.padRight;
49 import static java.lang.String.format;
50 import static org.apache.commons.math3.util.FastMath.abs;
51 import static org.apache.commons.math3.util.FastMath.acos;
52 import static org.apache.commons.math3.util.FastMath.cbrt;
53 import static org.apache.commons.math3.util.FastMath.cos;
54 import static org.apache.commons.math3.util.FastMath.floor;
55 import static org.apache.commons.math3.util.FastMath.min;
56 import static org.apache.commons.math3.util.FastMath.random;
57 import static org.apache.commons.math3.util.FastMath.signum;
58 import static org.apache.commons.math3.util.FastMath.sin;
59 import static org.apache.commons.math3.util.FastMath.sqrt;
60 import static org.apache.commons.math3.util.FastMath.toDegrees;
61 import static org.apache.commons.math3.util.FastMath.toRadians;
62
63 import ffx.utilities.FFXKeyword;
64 import java.util.ArrayList;
65 import java.util.List;
66 import java.util.Objects;
67 import java.util.logging.Level;
68 import java.util.logging.Logger;
69 import org.apache.commons.configuration2.CompositeConfiguration;
70 import org.apache.commons.math3.linear.Array2DRowRealMatrix;
71 import org.apache.commons.math3.linear.LUDecomposition;
72 import org.apache.commons.math3.linear.RealMatrix;
73
74
75
76
77
78
79
80
81
82
83 public class Crystal {
84
85 private static final Logger logger = Logger.getLogger(Crystal.class.getName());
86
87
88 @FFXKeyword(name = "SpaceGroup", keywordGroup = UnitCellAndSpaceGroup, clazz = SpaceGroup.class, defaultValue = "P1",
89 description = "This keyword selects the space group to be used in manipulation of crystal unit cells and asymmetric units.")
90 public final SpaceGroup spaceGroup;
91
92
93 @FFXKeyword(name = "a-axis", keywordGroup = UnitCellAndSpaceGroup, defaultValue = "None",
94 description = "Sets the value of the a-axis length for a crystal unit cell, or, equivalently, the X-axis length for a periodic box (Angstroms).")
95 public double a;
96
97
98 @FFXKeyword(name = "b-axis", keywordGroup = UnitCellAndSpaceGroup, defaultValue = "A-Axis",
99 description = "Sets the value of the b-axis length for a crystal unit cell, or, equivalently, the Y-axis length for a periodic box (Angstroms).")
100 public double b;
101
102
103 @FFXKeyword(name = "c-axis", keywordGroup = UnitCellAndSpaceGroup, defaultValue = "A-Axis",
104 description = "Sets the value of the c-axis length for a crystal unit cell, or, equivalently, the Z-axis length for a periodic box (Angstroms).")
105 public double c;
106
107
108 @FFXKeyword(name = "alpha", keywordGroup = UnitCellAndSpaceGroup, defaultValue = "90.0",
109 description = "Sets the value of the α-angle of a crystal unit cell, i.e., the angle between the b-axis and c-axis of a unit cell, or, equivalently, the angle between the Y-axis and Z-axis of a periodic box.")
110 public double alpha;
111
112
113 @FFXKeyword(name = "beta", keywordGroup = UnitCellAndSpaceGroup, defaultValue = "Alpha",
114 description = "Sets the value of the β-angle of a crystal unit cell, i.e., the angle between the a-axis and c-axis of a unit cell, or, equivalently, the angle between the X-axis and Z-axis of a periodic box.")
115 public double beta;
116
117
118 @FFXKeyword(name = "gamma", keywordGroup = UnitCellAndSpaceGroup, defaultValue = "Alpha",
119 description = "Sets the value of the γ-angle of a crystal unit cell, i.e., the angle between the a-axis and b-axis of a unit cell, or, equivalently, the angle between the X-axis and Y-axis of a periodic box.")
120 public double gamma;
121
122
123 private static final int XX = 0;
124
125 private static final int YY = 1;
126
127 private static final int ZZ = 2;
128
129
130
131
132
133
134
135 public final double[][] Ai = new double[3][3];
136
137 public final double[][] G = new double[3][3];
138
139 private final CrystalSystem crystalSystem;
140
141 private final LatticeSystem latticeSystem;
142
143 public double volume;
144
145 public double[][] A;
146
147 public double A00, A01, A02, A10, A11, A12, A20, A21, A22;
148
149 public double interfacialRadiusA;
150
151 public double interfacialRadiusB;
152
153 public double interfacialRadiusC;
154
155
156
157 public int[] scaleB = new int[6];
158
159
160
161 public int scaleN;
162
163 public double Ai00, Ai01, Ai02, Ai10, Ai11, Ai12, Ai20, Ai21, Ai22;
164
165 public double dVdA;
166
167 public double dVdB;
168
169 public double dVdC;
170
171
172
173 public double dVdAlpha;
174
175
176
177 public double dVdBeta;
178
179
180
181 public double dVdGamma;
182
183
184
185
186 boolean checkRestrictions = true;
187
188
189
190
191 private double specialPositionCutoff = 0.3;
192
193 private List<SymOp> symOpsCartesian;
194
195 private double[][] Gstar;
196
197
198
199 private double specialPositionCutoff2 = specialPositionCutoff * specialPositionCutoff;
200
201
202
203 private boolean aperiodic;
204
205
206
207
208
209
210
211
212
213
214
215
216
217 public Crystal(double a, double b, double c, double alpha, double beta, double gamma, int sgNumber) {
218 this(a, b, c, alpha, beta, gamma, SpaceGroupDefinitions.spaceGroupFactory(sgNumber).pdbName);
219 }
220
221
222
223
224
225
226
227
228
229
230
231
232
233 public Crystal(double a, double b, double c, double alpha, double beta, double gamma, String sg) {
234
235 SpaceGroup tempSG = SpaceGroupDefinitions.spaceGroupFactory(sg);
236 LatticeSystem tempLS = tempSG.latticeSystem;
237 this.a = a;
238 this.b = b;
239 this.c = c;
240 this.alpha = alpha;
241 this.beta = beta;
242 this.gamma = gamma;
243 aperiodic = false;
244 if (!tempLS.validParameters(a, b, c, alpha, beta, gamma)) {
245
246 StringBuilder sb = new StringBuilder(format(
247 " The %s lattice parameters do not satisfy the %s lattice system restrictions.\n",
248 tempSG.pdbName, tempLS));
249 sb.append(format(" A-axis: %18.15e\n", a));
250 sb.append(format(" B-axis: %18.15e\n", b));
251 sb.append(format(" C-axis: %18.15e\n", c));
252 sb.append(format(" Alpha: %18.15e\n", alpha));
253 sb.append(format(" Beta: %18.15e\n", beta));
254 sb.append(format(" Gamma: %18.15e\n", gamma));
255 logger.info(sb.toString());
256 sb = new StringBuilder();
257 if (tempLS == LatticeSystem.HEXAGONAL_LATTICE
258 || tempLS == LatticeSystem.RHOMBOHEDRAL_LATTICE) {
259
260 Crystal convertedCrystal = SpaceGroupConversions.hrConversion(a, b, c, alpha, beta, gamma,
261 tempSG);
262 this.a = convertedCrystal.a;
263 this.b = convertedCrystal.b;
264 this.c = convertedCrystal.c;
265 this.alpha = convertedCrystal.alpha;
266 this.beta = convertedCrystal.beta;
267 this.gamma = convertedCrystal.gamma;
268 spaceGroup = convertedCrystal.spaceGroup;
269 crystalSystem = spaceGroup.crystalSystem;
270 latticeSystem = spaceGroup.latticeSystem;
271 sb.append(" Converted ").append(tempSG.pdbName).append(" to ").append(spaceGroup.pdbName);
272 if (!latticeSystem.validParameters(this.a, this.b, this.c, this.alpha, this.beta,
273 this.gamma)) {
274 sb.append(format(
275 " The %s lattice parameters do not satisfy the %s lattice system restrictions.\n",
276 spaceGroup.pdbName, latticeSystem));
277 sb.append(format(" A-axis: %18.15e\n", this.a));
278 sb.append(format(" B-axis: %18.15e\n", this.b));
279 sb.append(format(" C-axis: %18.15e\n", this.c));
280 sb.append(format(" Alpha: %18.15e\n", this.alpha));
281 sb.append(format(" Beta: %18.15e\n", this.beta));
282 sb.append(format(" Gamma: %18.15e\n", this.gamma));
283 logger.severe(sb.toString());
284 } else {
285
286 logger.info(sb.toString());
287 }
288 } else {
289
290 this.a = a;
291 this.b = b;
292 this.c = c;
293 this.alpha = alpha;
294 this.beta = beta;
295 this.gamma = gamma;
296 spaceGroup = tempSG;
297 crystalSystem = spaceGroup.crystalSystem;
298 latticeSystem = spaceGroup.latticeSystem;
299 logger.severe(sb.toString());
300 }
301 } else {
302
303 this.a = a;
304 this.b = b;
305 this.c = c;
306 this.alpha = alpha;
307 this.beta = beta;
308 this.gamma = gamma;
309 spaceGroup = tempSG;
310 crystalSystem = spaceGroup.crystalSystem;
311 latticeSystem = spaceGroup.latticeSystem;
312 }
313
314 for (int i = 0; i < 6; i++) {
315 scaleB[i] = -1;
316 }
317
318 SymOp symop;
319 double[][] rot;
320 int index = 0;
321 switch (crystalSystem) {
322 case TRICLINIC:
323 for (int i = 0; i < 6; i++) {
324 scaleB[i] = index++;
325 }
326 break;
327 case MONOCLINIC:
328 scaleB[0] = index++;
329 scaleB[1] = index++;
330 scaleB[2] = index++;
331
332 symop = spaceGroup.symOps.get(1);
333 rot = symop.rot;
334 if (rot[0][0] > 0) {
335 scaleB[5] = index++;
336 } else if (rot[1][1] > 0) {
337 scaleB[4] = index++;
338 } else {
339 scaleB[3] = index++;
340 }
341 break;
342 case ORTHORHOMBIC:
343 scaleB[0] = index++;
344 scaleB[1] = index++;
345 scaleB[2] = index++;
346 break;
347 case TETRAGONAL:
348 scaleB[0] = index++;
349 scaleB[1] = scaleB[0];
350 scaleB[2] = index++;
351 break;
352 case TRIGONAL:
353 case HEXAGONAL:
354 boolean hexagonal = false;
355 for (int i = 1; i < spaceGroup.symOps.size(); i++) {
356 symop = spaceGroup.symOps.get(i);
357 rot = symop.rot;
358 index = 0;
359 if ((rot[1][1] * rot[1][2] == -1)
360 || (rot[2][1] * rot[2][2] == -1)
361 || (rot[1][1] * rot[1][2] == 1)
362 || (rot[2][1] * rot[2][2] == 1)) {
363 scaleB[0] = index++;
364 scaleB[1] = index++;
365 scaleB[2] = scaleB[1];
366 hexagonal = true;
367 } else if ((rot[0][0] * rot[0][2] == -1)
368 || (rot[2][0] * rot[2][2] == -1)
369 || (rot[0][0] * rot[0][2] == 1)
370 || (rot[2][0] * rot[2][2] == 1)) {
371 scaleB[0] = index++;
372 scaleB[1] = index++;
373 scaleB[2] = scaleB[0];
374 hexagonal = true;
375 } else if ((rot[0][0] * rot[0][1] == -1)
376 || (rot[1][0] * rot[1][1] == -1)
377 || (rot[0][0] * rot[0][1] == 1)
378 || (rot[1][0] * rot[1][1] == 1)) {
379 scaleB[0] = index++;
380 scaleB[1] = scaleB[0];
381 scaleB[2] = index++;
382 hexagonal = true;
383 }
384 if (hexagonal) {
385 break;
386 }
387 }
388 if (!hexagonal) {
389
390 scaleB[3] = index++;
391 scaleB[4] = scaleB[3];
392 scaleB[5] = scaleB[3];
393 }
394 break;
395 case CUBIC:
396 break;
397 }
398 scaleN = index;
399
400 updateCrystal();
401 }
402
403
404
405
406
407
408
409
410 public static Crystal checkProperties(CompositeConfiguration properties) {
411 double a = properties.getDouble("a-axis", -1.0);
412 double b = properties.getDouble("b-axis", -1.0);
413 double c = properties.getDouble("c-axis", -1.0);
414 double alpha = properties.getDouble("alpha", -1.0);
415 double beta = properties.getDouble("beta", -1.0);
416 double gamma = properties.getDouble("gamma", -1.0);
417 String sg = properties.getString("spacegroup", null);
418
419 sg = SpaceGroupInfo.pdb2ShortName(sg);
420
421 if (a < 0.0 || b < 0.0 || c < 0.0 || alpha < 0.0 || beta < 0.0 || gamma < 0.0 || sg == null) {
422 return null;
423 }
424
425
426 SpaceGroup spaceGroup = SpaceGroupDefinitions.spaceGroupFactory(sg);
427 if (spaceGroup == null) {
428 sg = sg.replaceAll(" ", "");
429 spaceGroup = SpaceGroupDefinitions.spaceGroupFactory(sg);
430 if (spaceGroup == null) {
431 return null;
432 }
433 }
434
435 return new Crystal(a, b, c, alpha, beta, gamma, sg);
436 }
437
438
439
440
441
442
443
444 public double invressq(HKL hkl) {
445 return hkl.quadForm(Gstar);
446 }
447
448
449
450
451
452
453
454 public double res(HKL hkl) {
455 return 1.0 / sqrt(hkl.quadForm(Gstar));
456 }
457
458
459
460
461
462
463 public boolean aperiodic() {
464 return aperiodic;
465 }
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483 public void applySymOp(int n, double[] x, double[] y, double[] z,
484 double[] mateX, double[] mateY, double[] mateZ, SymOp symOp) {
485 if (x == null || y == null || z == null) {
486 throw new IllegalArgumentException("The input arrays x, y and z must not be null.");
487 }
488 if (x.length < n || y.length < n || z.length < n) {
489 throw new IllegalArgumentException("The input arrays x, y and z must be of length n: " + n);
490 }
491 if (mateX == null || mateY == null || mateZ == null) {
492 throw new IllegalArgumentException("The output arrays mateX, mateY and mateZ must not be null.");
493 }
494 if (mateX.length < n || mateY.length < n || mateZ.length < n) {
495 throw new IllegalArgumentException("The output arrays mateX, mateY and mateZ must be of length n: " + n);
496 }
497 final double[][] rot = symOp.rot;
498 final double[] trans = symOp.tr;
499
500 final double rot00 = rot[0][0];
501 final double rot10 = rot[1][0];
502 final double rot20 = rot[2][0];
503 final double rot01 = rot[0][1];
504 final double rot11 = rot[1][1];
505 final double rot21 = rot[2][1];
506 final double rot02 = rot[0][2];
507 final double rot12 = rot[1][2];
508 final double rot22 = rot[2][2];
509 final double t0 = trans[0];
510 final double t1 = trans[1];
511 final double t2 = trans[2];
512 for (int i = 0; i < n; i++) {
513 double xc = x[i];
514 double yc = y[i];
515 double zc = z[i];
516
517 double xi = xc * A00 + yc * A10 + zc * A20;
518 double yi = xc * A01 + yc * A11 + zc * A21;
519 double zi = xc * A02 + yc * A12 + zc * A22;
520
521 double fx = rot00 * xi + rot01 * yi + rot02 * zi + t0;
522 double fy = rot10 * xi + rot11 * yi + rot12 * zi + t1;
523 double fz = rot20 * xi + rot21 * yi + rot22 * zi + t2;
524
525 mateX[i] = fx * Ai00 + fy * Ai10 + fz * Ai20;
526 mateY[i] = fx * Ai01 + fy * Ai11 + fz * Ai21;
527 mateZ[i] = fx * Ai02 + fy * Ai12 + fz * Ai22;
528 }
529 }
530
531
532
533
534
535
536
537
538 public void applySymOp(double[] xyz, double[] mate, SymOp symOp) {
539 assert (xyz.length % 3 == 0);
540 assert (xyz.length == mate.length);
541
542 var rot = symOp.rot;
543 var r00 = rot[0][0];
544 var r01 = rot[0][1];
545 var r02 = rot[0][2];
546 var r10 = rot[1][0];
547 var r11 = rot[1][1];
548 var r12 = rot[1][2];
549 var r20 = rot[2][0];
550 var r21 = rot[2][1];
551 var r22 = rot[2][2];
552
553 var trans = symOp.tr;
554 var t0 = trans[0];
555 var t1 = trans[1];
556 var t2 = trans[2];
557
558 int len = xyz.length / 3;
559 for (int i = 0; i < len; i++) {
560 int index = i * 3;
561 var xc = xyz[index + XX];
562 var yc = xyz[index + YY];
563 var zc = xyz[index + ZZ];
564
565 var xi = xc * A00 + yc * A10 + zc * A20;
566 var yi = xc * A01 + yc * A11 + zc * A21;
567 var zi = xc * A02 + yc * A12 + zc * A22;
568
569 var fx = r00 * xi + r01 * yi + r02 * zi + t0;
570 var fy = r10 * xi + r11 * yi + r12 * zi + t1;
571 var fz = r20 * xi + r21 * yi + r22 * zi + t2;
572
573 mate[index + XX] = fx * Ai00 + fy * Ai10 + fz * Ai20;
574 mate[index + YY] = fx * Ai01 + fy * Ai11 + fz * Ai21;
575 mate[index + ZZ] = fx * Ai02 + fy * Ai12 + fz * Ai22;
576 }
577 }
578
579
580
581
582
583
584
585
586 public void applySymRot(double[] xyz, double[] mate, SymOp symOp) {
587 double[][] rot = symOp.rot;
588
589 double xc = xyz[0];
590 double yc = xyz[1];
591 double zc = xyz[2];
592
593 double xi = xc * A00 + yc * A10 + zc * A20;
594 double yi = xc * A01 + yc * A11 + zc * A21;
595 double zi = xc * A02 + yc * A12 + zc * A22;
596
597
598 double fx = rot[0][0] * xi + rot[0][1] * yi + rot[0][2] * zi;
599 double fy = rot[1][0] * xi + rot[1][1] * yi + rot[1][2] * zi;
600 double fz = rot[2][0] * xi + rot[2][1] * yi + rot[2][2] * zi;
601
602
603 mate[0] = fx * Ai00 + fy * Ai10 + fz * Ai20;
604 mate[1] = fx * Ai01 + fy * Ai11 + fz * Ai21;
605 mate[2] = fx * Ai02 + fy * Ai12 + fz * Ai22;
606 }
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625 public void applyTransSymRot(
626 int n, double[] x, double[] y, double[] z,
627 double[] mateX, double[] mateY, double[] mateZ, SymOp symOp, double[][] rotmat) {
628
629 if (x == null || y == null || z == null) {
630 throw new IllegalArgumentException("The input arrays x, y and z must not be null.");
631 }
632 if (x.length < n || y.length < n || z.length < n) {
633 throw new IllegalArgumentException("The input arrays x, y and z must be of length n: " + n);
634 }
635 if (mateX == null || mateY == null || mateZ == null) {
636 throw new IllegalArgumentException("The output arrays mateX, mateY and mateZ must not be null.");
637 }
638 if (mateX.length < n || mateY.length < n || mateZ.length < n) {
639 throw new IllegalArgumentException("The output arrays mateX, mateY and mateZ must be of length n: " + n);
640 }
641
642
643 getTransformationOperator(symOp, rotmat);
644
645 for (int i = 0; i < n; i++) {
646
647 double xc = x[i];
648 double yc = y[i];
649 double zc = z[i];
650 mateX[i] = xc * rotmat[0][0] + yc * rotmat[1][0] + zc * rotmat[2][0];
651 mateY[i] = xc * rotmat[0][1] + yc * rotmat[1][1] + zc * rotmat[2][1];
652 mateZ[i] = xc * rotmat[0][2] + yc * rotmat[1][2] + zc * rotmat[2][2];
653 }
654 }
655
656
657
658
659
660
661
662 public void averageTensor(double[][] m, double[][] r) {
663 int n = spaceGroup.symOps.size();
664 for (int i = 0; i < n; i++) {
665 SymOp symop = spaceGroup.symOps.get(i);
666 double[][] rot = symop.rot;
667 double[][] rt = transpose3(rot);
668 double[][] rmrt = mat3Mat3(mat3Mat3(rot, m), rt);
669 for (int j = 0; j < 3; j++) {
670 for (int k = 0; k < 3; k++) {
671 r[j][k] += rmrt[j][k];
672 }
673 }
674 }
675 for (int j = 0; j < 3; j++) {
676 for (int k = 0; k < 3; k++) {
677 r[j][k] /= 6.0;
678 }
679 }
680 }
681
682
683
684
685
686
687
688 public void averageTensor(double[] v, double[][] r) {
689 int n = spaceGroup.symOps.size();
690 for (int i = 0; i < n; i++) {
691 SymOp symop = spaceGroup.symOps.get(i);
692 double[][] rot = symop.rot;
693 double[][] rt = transpose3(rot);
694 double[][] rmrt = mat3Mat3(mat3SymVec6(rot, v), rt);
695 for (int j = 0; j < 3; j++) {
696 for (int k = 0; k < 3; k++) {
697 r[j][k] += rmrt[j][k];
698 }
699 }
700 }
701 for (int j = 0; j < 3; j++) {
702 for (int k = 0; k < 3; k++) {
703 r[j][k] /= 6.0;
704 }
705 }
706 }
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721 public boolean changeUnitCellParameters(
722 double a, double b, double c, double alpha, double beta, double gamma) {
723 if (checkRestrictions) {
724 if (!latticeSystem.validParameters(a, b, c, alpha, beta, gamma)) {
725 if (logger.isLoggable(Level.FINE)) {
726 StringBuilder sb = new StringBuilder(
727 " The proposed lattice parameters for " + spaceGroup.pdbName
728 + " do not satisfy the " + latticeSystem +
729 " lattice system restrictions and were ignored.\n");
730 sb.append(format(" A-axis: %18.15e\n", a));
731 sb.append(format(" B-axis: %18.15e\n", b));
732 sb.append(format(" C-axis: %18.15e\n", c));
733 sb.append(format(" Alpha: %18.15e\n", alpha));
734 sb.append(format(" Beta: %18.15e\n", beta));
735 sb.append(format(" Gamma: %18.15e\n", gamma));
736 logger.fine(sb.toString());
737 }
738 return false;
739 }
740 }
741
742 this.a = a;
743 this.b = b;
744 this.c = c;
745 this.alpha = alpha;
746 this.beta = beta;
747 this.gamma = gamma;
748
749 updateCrystal();
750
751 return true;
752 }
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771 public boolean changeUnitCellParametersAndVolume(
772 double a, double b, double c, double alpha, double beta, double gamma, double targetAUVolume) {
773 if (changeUnitCellParameters(a, b, c, alpha, beta, gamma)) {
774 double currentAUVolume = volume / getNumSymOps();
775 double scale = cbrt(targetAUVolume / currentAUVolume);
776 return changeUnitCellParameters(scale * a, scale * b, scale * c, alpha, beta, gamma);
777 }
778 return false;
779 }
780
781
782
783
784
785
786
787 @Override
788 public boolean equals(Object o) {
789 if (this == o) {
790 return true;
791 }
792 if (o == null || getClass() != o.getClass()) {
793 return false;
794 }
795 Crystal crystal = (Crystal) o;
796 return (a == crystal.a
797 && b == crystal.b
798 && c == crystal.c
799 && alpha == crystal.alpha
800 && beta == crystal.beta
801 && gamma == crystal.gamma
802 && spaceGroup.number == crystal.spaceGroup.number);
803 }
804
805 public boolean getCheckRestrictions() {
806 return checkRestrictions;
807 }
808
809 public void setCheckRestrictions(boolean checkRestrictions) {
810 this.checkRestrictions = checkRestrictions;
811 }
812
813
814
815
816
817
818
819 public double getDensity(double mass) {
820 int nSymm = spaceGroup.symOps.size();
821 return (mass * nSymm / AVOGADRO) * (1.0e24 / volume);
822 }
823
824
825
826
827
828
829 public int getNumSymOps() {
830 return spaceGroup.getNumberOfSymOps();
831 }
832
833
834
835
836
837
838
839
840
841 public double[] getRandomCartTranslation() {
842 double[] coords = {random(), random(), random()};
843 toCartesianCoordinates(coords, coords);
844 return coords;
845 }
846
847
848
849
850
851
852 public double getSpecialPositionCutoff() {
853 return specialPositionCutoff;
854 }
855
856
857
858
859
860
861 public double getSpecialPositionCutoff2() {
862 return specialPositionCutoff2;
863 }
864
865
866
867
868
869
870 public void setSpecialPositionCutoff(double cutoff) {
871 specialPositionCutoff = cutoff;
872 specialPositionCutoff2 = cutoff * cutoff;
873 }
874
875
876
877
878
879
880
881 public void getTransformationOperator(SymOp symOp, double[][] rotmat) {
882 double[][] rot = symOp.rot;
883 for (int i = 0; i < 3; i++) {
884 for (int j = 0; j < 3; j++) {
885 rotmat[i][j] = 0.0;
886 for (int k = 0; k < 3; k++) {
887 for (int l = 0; l < 3; l++) {
888 rotmat[i][j] += Ai[k][i] * rot[k][l] * A[j][l];
889 }
890 }
891 }
892 }
893 }
894
895
896
897
898
899
900
901 public Crystal getUnitCell() {
902 return this;
903 }
904
905
906 @Override
907 public int hashCode() {
908 return Objects.hash(a, b, c, alpha, beta, gamma, spaceGroup.number);
909 }
910
911
912
913
914
915
916
917
918 public double image(double[] xyz, double[] xyz2) {
919 double dx = xyz[0] - xyz2[0];
920 double dy = xyz[1] - xyz2[1];
921 double dz = xyz[2] - xyz2[2];
922 return image(dx, dy, dz);
923 }
924
925
926
927
928
929
930
931 public double image(final double[] xyz) {
932 double x = xyz[0];
933 double y = xyz[1];
934 double z = xyz[2];
935 if (aperiodic) {
936 return x * x + y * y + z * z;
937 }
938 double xf = x * A00 + y * A10 + z * A20;
939 double yf = x * A01 + y * A11 + z * A21;
940 double zf = x * A02 + y * A12 + z * A22;
941 xf = floor(abs(xf) + 0.5) * signum(-xf) + xf;
942 yf = floor(abs(yf) + 0.5) * signum(-yf) + yf;
943 zf = floor(abs(zf) + 0.5) * signum(-zf) + zf;
944 x = xf * Ai00 + yf * Ai10 + zf * Ai20;
945 y = xf * Ai01 + yf * Ai11 + zf * Ai21;
946 z = xf * Ai02 + yf * Ai12 + zf * Ai22;
947 xyz[0] = x;
948 xyz[1] = y;
949 xyz[2] = z;
950 return x * x + y * y + z * z;
951 }
952
953
954
955
956
957
958
959
960
961 public double image(double dx, double dy, double dz) {
962 if (aperiodic) {
963 return dx * dx + dy * dy + dz * dz;
964 }
965 double xf = dx * A00 + dy * A10 + dz * A20;
966 double yf = dx * A01 + dy * A11 + dz * A21;
967 double zf = dx * A02 + dy * A12 + dz * A22;
968 xf = floor(abs(xf) + 0.5) * signum(-xf) + xf;
969 yf = floor(abs(yf) + 0.5) * signum(-yf) + yf;
970 zf = floor(abs(zf) + 0.5) * signum(-zf) + zf;
971 dx = xf * Ai00 + yf * Ai10 + zf * Ai20;
972 dy = xf * Ai01 + yf * Ai11 + zf * Ai21;
973 dz = xf * Ai02 + yf * Ai12 + zf * Ai22;
974 return dx * dx + dy * dy + dz * dz;
975 }
976
977
978
979
980
981
982
983
984 public double minDistOverSymOps(double[] xyzA, double[] xyzB) {
985 double dist = 0;
986 for (int i = 0; i < 3; i++) {
987 double dx = xyzA[i] - xyzB[i];
988 dist += (dx * dx);
989 }
990 var symB = new double[3];
991 for (SymOp symOp : spaceGroup.symOps) {
992 applySymOp(xyzB, symB, symOp);
993 for (int i = 0; i < 3; i++) {
994 symB[i] -= xyzA[i];
995 }
996 double d = image(symB);
997 dist = min(d, dist);
998 }
999 return sqrt(dist);
1000 }
1001
1002
1003
1004
1005
1006
1007
1008 public boolean perturbCellVectors(double[][] dStrain) {
1009
1010 double[][] AA = new double[3][3];
1011 double[][] newAi = new double[3][3];
1012
1013 for (int i = 0; i < 3; i++) {
1014 for (int j = 0; j < 3; j++) {
1015 AA[i][j] = getUnitCell().Ai[i][j];
1016 }
1017 }
1018
1019 for (int i = 0; i < 3; i++) {
1020 for (int j = 0; j < 3; j++) {
1021 for (int k = 0; k < 3; k++) {
1022 double Kronecker = 0.0;
1023 if (j == k) {
1024 Kronecker = 1.0;
1025 }
1026
1027 newAi[i][j] += (Kronecker + dStrain[j][k]) * AA[i][k];
1028 }
1029 }
1030 }
1031
1032 return setCellVectors(newAi);
1033 }
1034
1035 public boolean randomParameters(double dens, double mass) {
1036 double[] params = latticeSystem.resetUnitCellParams();
1037 boolean succeed = changeUnitCellParameters(params[0], params[1], params[2], params[3], params[4], params[5]);
1038 if (succeed) {
1039 setDensity(dens, mass);
1040 }
1041 return succeed;
1042 }
1043
1044
1045
1046
1047
1048
1049 public void setAperiodic(boolean aperiodic) {
1050 this.aperiodic = aperiodic;
1051 }
1052
1053 public double[] getCellParametersFromVectors(double[][] cellVectors) {
1054
1055 double aa = length(cellVectors[0]);
1056 double bb = length(cellVectors[1]);
1057 double cc = length(cellVectors[2]);
1058
1059
1060 double aalpha = toDegrees(acos(dot(cellVectors[1], cellVectors[2]) / (bb * cc)));
1061 double bbeta = toDegrees(acos(dot(cellVectors[0], cellVectors[2]) / (aa * cc)));
1062 double ggamma = toDegrees(acos(dot(cellVectors[0], cellVectors[1]) / (aa * bb)));
1063
1064
1065 double[] params = new double[6];
1066 params[0] = aa;
1067 params[1] = bb;
1068 params[2] = cc;
1069 params[3] = aalpha;
1070 params[4] = bbeta;
1071 params[5] = ggamma;
1072 return params;
1073 }
1074
1075
1076
1077
1078
1079
1080
1081 public boolean setCellVectors(double[][] cellVectors) {
1082 double[] p = getCellParametersFromVectors(cellVectors);
1083 return changeUnitCellParameters(p[0], p[1], p[2], p[3], p[4], p[5]);
1084 }
1085
1086
1087
1088
1089
1090
1091
1092
1093 public boolean setCellVectorsAndVolume(double[][] cellVectors, double targetAUVolume) {
1094 if (setCellVectors(cellVectors)) {
1095 double currentAUVolume = volume / getNumSymOps();
1096 double scale = cbrt(targetAUVolume / currentAUVolume);
1097 return changeUnitCellParameters(scale * a, scale * b, scale * c, alpha, beta, gamma);
1098 } else {
1099 return false;
1100 }
1101 }
1102
1103 public void setDensity(double dens, double mass) {
1104 double currentDensity = getDensity(mass);
1105 double scale = cbrt(currentDensity / dens);
1106 changeUnitCellParameters(a * scale, b * scale, c * scale, alpha, beta, gamma);
1107 currentDensity = getDensity(mass);
1108 if (logger.isLoggable(Level.FINE)) {
1109 logger.fine(format(" Updated density %6.3f (g/cc) with unit cell %s.", currentDensity, toShortString()));
1110 }
1111 }
1112
1113
1114
1115
1116
1117
1118 public String toCRYST1() {
1119 return format(
1120 "CRYST1%9.3f%9.3f%9.3f%7.2f%7.2f%7.2f %10s\n",
1121 a, b, c, alpha, beta, gamma, padRight(spaceGroup.pdbName, 10));
1122 }
1123
1124
1125
1126
1127
1128
1129
1130
1131
1132
1133
1134
1135 public void toCartesianCoordinates(
1136 int n, double[] xf, double[] yf, double[] zf, double[] x, double[] y, double[] z) {
1137 for (int i = 0; i < n; i++) {
1138 double xi = xf[i];
1139 double yi = yf[i];
1140 double zi = zf[i];
1141 x[i] = xi * Ai00 + yi * Ai10 + zi * Ai20;
1142 y[i] = xi * Ai01 + yi * Ai11 + zi * Ai21;
1143 z[i] = xi * Ai02 + yi * Ai12 + zi * Ai22;
1144 }
1145 }
1146
1147
1148
1149
1150
1151
1152
1153
1154 public void toCartesianCoordinates(int n, double[] frac, double[] cart) {
1155 int i3 = 0;
1156 for (int i = 0; i < n; i++) {
1157
1158 int iX = i3 + XX;
1159 int iY = i3 + YY;
1160 int iZ = i3 + ZZ;
1161 i3 += 3;
1162 double xf = frac[iX];
1163 double yf = frac[iY];
1164 double zf = frac[iZ];
1165 cart[iX] = xf * Ai00 + yf * Ai10 + zf * Ai20;
1166 cart[iY] = xf * Ai01 + yf * Ai11 + zf * Ai21;
1167 cart[iZ] = xf * Ai02 + yf * Ai12 + zf * Ai22;
1168 }
1169 }
1170
1171
1172
1173
1174
1175
1176
1177 public void toCartesianCoordinates(double[] xf, double[] x) {
1178 double fx = xf[0];
1179 double fy = xf[1];
1180 double fz = xf[2];
1181 x[0] = fx * Ai00 + fy * Ai10 + fz * Ai20;
1182 x[1] = fx * Ai01 + fy * Ai11 + fz * Ai21;
1183 x[2] = fx * Ai02 + fy * Ai12 + fz * Ai22;
1184 }
1185
1186
1187
1188
1189
1190
1191
1192
1193
1194
1195
1196
1197 public void toFractionalCoordinates(
1198 int n, double[] x, double[] y, double[] z, double[] xf, double[] yf, double[] zf) {
1199 for (int i = 0; i < n; i++) {
1200 double xc = x[i];
1201 double yc = y[i];
1202 double zc = z[i];
1203 xf[i] = xc * A00 + yc * A10 + zc * A20;
1204 yf[i] = xc * A01 + yc * A11 + zc * A21;
1205 zf[i] = xc * A02 + yc * A12 + zc * A22;
1206 }
1207 }
1208
1209
1210
1211
1212
1213
1214
1215
1216 public void toFractionalCoordinates(int n, double[] cart, double[] frac) {
1217 int i3 = 0;
1218 for (int i = 0; i < n; i++) {
1219
1220 int iX = i3 + XX;
1221 int iY = i3 + YY;
1222 int iZ = i3 + ZZ;
1223 i3 += 3;
1224 double xc = cart[iX];
1225 double yc = cart[iY];
1226 double zc = cart[iZ];
1227 frac[iX] = xc * A00 + yc * A10 + zc * A20;
1228 frac[iY] = xc * A01 + yc * A11 + zc * A21;
1229 frac[iZ] = xc * A02 + yc * A12 + zc * A22;
1230 }
1231 }
1232
1233
1234
1235
1236
1237
1238
1239 public void toFractionalCoordinates(double[] x, double[] xf) {
1240 double xc = x[0];
1241 double yc = x[1];
1242 double zc = x[2];
1243 xf[0] = xc * A00 + yc * A10 + zc * A20;
1244 xf[1] = xc * A01 + yc * A11 + zc * A21;
1245 xf[2] = xc * A02 + yc * A12 + zc * A22;
1246 }
1247
1248
1249
1250
1251
1252
1253
1254 public void toPrimaryCell(double[] in, double[] out) {
1255 toFractionalCoordinates(in, out);
1256 out[0] = mod(out[0], 1.0);
1257 out[1] = mod(out[1], 1.0);
1258 out[2] = mod(out[2], 1.0);
1259 toCartesianCoordinates(out, out);
1260 }
1261
1262
1263
1264
1265
1266
1267 public String toShortString() {
1268 return format("%6.2f %6.2f %6.2f %6.2f %6.2f %6.2f", a, b, c, alpha, beta, gamma);
1269 }
1270
1271
1272 @Override
1273 public String toString() {
1274 StringBuilder sb = new StringBuilder("\n Unit Cell\n");
1275 sb.append(
1276 format(" A-axis: %8.3f (%8.3f, %8.3f, %8.3f)\n", a, Ai00, Ai01,
1277 Ai02));
1278 sb.append(
1279 format(" B-axis: %8.3f (%8.3f, %8.3f, %8.3f)\n", b, Ai10, Ai11,
1280 Ai12));
1281 sb.append(
1282 format(" C-axis: %8.3f (%8.3f, %8.3f, %8.3f)\n", c, Ai20, Ai21,
1283 Ai22));
1284 sb.append(format(" Alpha: %8.3f\n", alpha));
1285 sb.append(format(" Beta: %8.3f\n", beta));
1286 sb.append(format(" Gamma: %8.3f\n", gamma));
1287 sb.append(" Space group\n");
1288 sb.append(format(" Number: %3d\n", spaceGroup.number));
1289 sb.append(format(" Symbol: %8s\n", spaceGroup.shortName));
1290 sb.append(
1291 format(" Number of Symmetry Operators: %3d", spaceGroup.getNumberOfSymOps()));
1292 return sb.toString();
1293 }
1294
1295
1296 public void updateCrystal() {
1297
1298 double cos_alpha;
1299 double sin_beta;
1300 double cos_beta;
1301 double sin_gamma;
1302 double cos_gamma;
1303 double beta_term;
1304 double gamma_term;
1305
1306 switch (crystalSystem) {
1307 default -> {
1308
1309 cos_alpha = 0.0;
1310 cos_beta = 0.0;
1311 sin_gamma = 1.0;
1312 cos_gamma = 0.0;
1313 beta_term = 0.0;
1314 gamma_term = 1.0;
1315 volume = a * b * c;
1316 dVdA = b * c;
1317 dVdB = a * c;
1318 dVdC = a * b;
1319 dVdAlpha = 0.0;
1320 dVdBeta = 0.0;
1321 dVdGamma = 0.0;
1322 }
1323 case MONOCLINIC -> {
1324 cos_alpha = 0.0;
1325 sin_beta = sin(toRadians(beta));
1326 cos_beta = cos(toRadians(beta));
1327 sin_gamma = 1.0;
1328 cos_gamma = 0.0;
1329 beta_term = 0.0;
1330 gamma_term = sin_beta;
1331 volume = sin_beta * a * b * c;
1332 dVdA = sin_beta * b * c;
1333 dVdB = sin_beta * a * c;
1334 dVdC = sin_beta * a * b;
1335 dVdAlpha = 0.0;
1336 dVdBeta = cos_beta * a * b * c;
1337 dVdGamma = 0.0;
1338 }
1339 case HEXAGONAL, TRICLINIC, TRIGONAL -> {
1340 double sin_alpha = sin(toRadians(alpha));
1341 cos_alpha = cos(toRadians(alpha));
1342 sin_beta = sin(toRadians(beta));
1343 cos_beta = cos(toRadians(beta));
1344 sin_gamma = sin(toRadians(gamma));
1345 cos_gamma = cos(toRadians(gamma));
1346 beta_term = (cos_alpha - cos_beta * cos_gamma) / sin_gamma;
1347 gamma_term = sqrt(sin_beta * sin_beta - beta_term * beta_term);
1348 volume = sin_gamma * gamma_term * a * b * c;
1349 dVdA = sin_gamma * gamma_term * b * c;
1350 dVdB = sin_gamma * gamma_term * a * c;
1351 dVdC = sin_gamma * gamma_term * a * b;
1352 double dbeta = 2.0 * sin_beta * cos_beta
1353 - (2.0 * (cos_alpha - cos_beta * cos_gamma) * sin_beta * cos_gamma)
1354 / (sin_gamma * sin_gamma);
1355 double dgamma1 = -2.0 * (cos_alpha - cos_beta * cos_gamma) * cos_beta / sin_gamma;
1356 double dgamma2 = cos_alpha - cos_beta * cos_gamma;
1357 dgamma2 *= dgamma2 * 2.0 * cos_gamma / (sin_gamma * sin_gamma * sin_gamma);
1358 dVdAlpha = (cos_alpha - cos_beta * cos_gamma) * sin_alpha
1359 / (sin_gamma * gamma_term) * a * b * c;
1360 dVdBeta = 0.5 * sin_gamma * dbeta / gamma_term * a * b * c;
1361 dVdGamma = (cos_gamma * gamma_term + 0.5 * sin_gamma * (dgamma1 + dgamma2) / gamma_term)
1362 * a * b * c;
1363 }
1364 }
1365
1366 G[0][0] = a * a;
1367 G[0][1] = a * b * cos_gamma;
1368 G[0][2] = a * c * cos_beta;
1369 G[1][0] = G[0][1];
1370 G[1][1] = b * b;
1371 G[1][2] = b * c * cos_alpha;
1372 G[2][0] = G[0][2];
1373 G[2][1] = G[1][2];
1374 G[2][2] = c * c;
1375
1376
1377 RealMatrix m = new Array2DRowRealMatrix(G, true);
1378 m = new LUDecomposition(m).getSolver().getInverse();
1379 Gstar = m.getData();
1380
1381
1382 Ai[0][0] = a;
1383 Ai[0][1] = 0.0;
1384 Ai[0][2] = 0.0;
1385
1386 Ai[1][0] = b * cos_gamma;
1387 Ai[1][1] = b * sin_gamma;
1388 Ai[1][2] = 0.0;
1389
1390 Ai[2][0] = c * cos_beta;
1391 Ai[2][1] = c * beta_term;
1392 Ai[2][2] = c * gamma_term;
1393
1394 Ai00 = Ai[0][0];
1395 Ai01 = Ai[0][1];
1396 Ai02 = Ai[0][2];
1397 Ai10 = Ai[1][0];
1398 Ai11 = Ai[1][1];
1399 Ai12 = Ai[1][2];
1400 Ai20 = Ai[2][0];
1401 Ai21 = Ai[2][1];
1402 Ai22 = Ai[2][2];
1403
1404
1405 m = new Array2DRowRealMatrix(Ai, true);
1406 m = new LUDecomposition(m).getSolver().getInverse();
1407 A = m.getData();
1408
1409
1410 A00 = A[0][0];
1411 A10 = A[1][0];
1412 A20 = A[2][0];
1413 A01 = A[0][1];
1414 A11 = A[1][1];
1415 A21 = A[2][1];
1416 A02 = A[0][2];
1417 A12 = A[1][2];
1418 A22 = A[2][2];
1419
1420
1421 double aStar = 1.0 / sqrt(A00 * A00 + A10 * A10 + A20 * A20);
1422 double bStar = 1.0 / sqrt(A01 * A01 + A11 * A11 + A21 * A21);
1423 double cStar = 1.0 / sqrt(A02 * A02 + A12 * A12 + A22 * A22);
1424 if (logger.isLoggable(Level.FINEST)) {
1425 logger.finest(
1426 format(" Reciprocal Lattice Lengths: (%8.3f, %8.3f, %8.3f)", aStar, bStar, cStar));
1427 }
1428
1429
1430 interfacialRadiusA = (Ai00 * A00 + Ai01 * A10 + Ai02 * A20) * aStar;
1431 interfacialRadiusB = (Ai10 * A01 + Ai11 * A11 + Ai12 * A21) * bStar;
1432 interfacialRadiusC = (Ai20 * A02 + Ai21 * A12 + Ai22 * A22) * cStar;
1433
1434
1435 interfacialRadiusA /= 2.0;
1436 interfacialRadiusB /= 2.0;
1437 interfacialRadiusC /= 2.0;
1438
1439 if (logger.isLoggable(Level.FINEST)) {
1440 logger.finest(
1441 format(
1442 " Interfacial radii: (%8.3f, %8.3f, %8.3f)",
1443 interfacialRadiusA, interfacialRadiusB, interfacialRadiusC));
1444 }
1445
1446 List<SymOp> symOps = spaceGroup.symOps;
1447 int nSymm = symOps.size();
1448 if (symOpsCartesian == null) {
1449 symOpsCartesian = new ArrayList<>(nSymm);
1450 } else {
1451 symOpsCartesian.clear();
1452 }
1453
1454 RealMatrix toFrac = new Array2DRowRealMatrix(A);
1455 RealMatrix toCart = new Array2DRowRealMatrix(Ai);
1456 for (SymOp symOp : symOps) {
1457 m = new Array2DRowRealMatrix(symOp.rot);
1458
1459 RealMatrix rotMat = m.preMultiply(toCart.transpose()).multiply(toFrac.transpose());
1460
1461 double[] tr = toCart.preMultiply(symOp.tr);
1462 symOpsCartesian.add(new SymOp(rotMat.getData(), tr));
1463 }
1464 }
1465 }