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