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.PropertyGroup.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.FFXProperty;
64
65 import java.util.ArrayList;
66 import java.util.List;
67 import java.util.Objects;
68 import java.util.logging.Level;
69 import java.util.logging.Logger;
70
71 import org.apache.commons.configuration2.CompositeConfiguration;
72 import org.apache.commons.math3.linear.Array2DRowRealMatrix;
73 import org.apache.commons.math3.linear.LUDecomposition;
74 import org.apache.commons.math3.linear.RealMatrix;
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, SpaceGroupDefinitions.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 = SpaceGroupDefinitions.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 = SpaceGroupDefinitions.spaceGroupFactory(sg);
502 if (spaceGroup == null) {
503 sg = sg.replaceAll(" ", "");
504 spaceGroup = SpaceGroupDefinitions.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 public void averageTensor(double[][] m, double[][] r) {
738 int n = spaceGroup.symOps.size();
739 for (int i = 0; i < n; i++) {
740 SymOp symop = spaceGroup.symOps.get(i);
741 double[][] rot = symop.rot;
742 double[][] rt = transpose3(rot);
743 double[][] rmrt = mat3Mat3(mat3Mat3(rot, m), rt);
744 for (int j = 0; j < 3; j++) {
745 for (int k = 0; k < 3; k++) {
746 r[j][k] += rmrt[j][k];
747 }
748 }
749 }
750 for (int j = 0; j < 3; j++) {
751 for (int k = 0; k < 3; k++) {
752 r[j][k] /= 6.0;
753 }
754 }
755 }
756
757
758
759
760
761
762
763 public void averageTensor(double[] v, double[][] r) {
764 int n = spaceGroup.symOps.size();
765 for (int i = 0; i < n; i++) {
766 SymOp symop = spaceGroup.symOps.get(i);
767 double[][] rot = symop.rot;
768 double[][] rt = transpose3(rot);
769 double[][] rmrt = mat3Mat3(mat3SymVec6(rot, v), rt);
770 for (int j = 0; j < 3; j++) {
771 for (int k = 0; k < 3; k++) {
772 r[j][k] += rmrt[j][k];
773 }
774 }
775 }
776 for (int j = 0; j < 3; j++) {
777 for (int k = 0; k < 3; k++) {
778 r[j][k] /= 6.0;
779 }
780 }
781 }
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796 public boolean changeUnitCellParameters(
797 double a, double b, double c, double alpha, double beta, double gamma) {
798 if (checkRestrictions) {
799 if (!latticeSystem.validParameters(a, b, c, alpha, beta, gamma)) {
800 if (logger.isLoggable(Level.FINE)) {
801 StringBuilder sb = new StringBuilder(
802 " The proposed lattice parameters for " + spaceGroup.pdbName
803 + " do not satisfy the " + latticeSystem +
804 " lattice system restrictions and were ignored.\n");
805 sb.append(format(" A-axis: %18.15e\n", a));
806 sb.append(format(" B-axis: %18.15e\n", b));
807 sb.append(format(" C-axis: %18.15e\n", c));
808 sb.append(format(" Alpha: %18.15e\n", alpha));
809 sb.append(format(" Beta: %18.15e\n", beta));
810 sb.append(format(" Gamma: %18.15e\n", gamma));
811 logger.fine(sb.toString());
812 }
813 return false;
814 }
815 }
816
817 this.a = a;
818 this.b = b;
819 this.c = c;
820 this.alpha = alpha;
821 this.beta = beta;
822 this.gamma = gamma;
823
824 updateCrystal();
825
826 return true;
827 }
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846 public boolean changeUnitCellParametersAndVolume(
847 double a, double b, double c, double alpha, double beta, double gamma, double targetAUVolume) {
848 if (changeUnitCellParameters(a, b, c, alpha, beta, gamma)) {
849 double currentAUVolume = volume / getNumSymOps();
850 double scale = cbrt(targetAUVolume / currentAUVolume);
851 return changeUnitCellParameters(scale * a, scale * b, scale * c, alpha, beta, gamma);
852 }
853 return false;
854 }
855
856
857
858
859
860
861
862 @Override
863 public boolean equals(Object o) {
864 if (this == o) {
865 return true;
866 }
867 if (o == null || getClass() != o.getClass()) {
868 return false;
869 }
870 Crystal crystal = (Crystal) o;
871 return (a == crystal.a
872 && b == crystal.b
873 && c == crystal.c
874 && alpha == crystal.alpha
875 && beta == crystal.beta
876 && gamma == crystal.gamma
877 && spaceGroup.number == crystal.spaceGroup.number);
878 }
879
880 public boolean getCheckRestrictions() {
881 return checkRestrictions;
882 }
883
884 public void setCheckRestrictions(boolean checkRestrictions) {
885 this.checkRestrictions = checkRestrictions;
886 }
887
888
889
890
891
892
893
894 public double getDensity(double mass) {
895 int nSymm = spaceGroup.symOps.size();
896 return (mass * nSymm / AVOGADRO) * (1.0e24 / volume);
897 }
898
899
900
901
902
903
904 public int getNumSymOps() {
905 return spaceGroup.getNumberOfSymOps();
906 }
907
908
909
910
911
912
913
914
915
916 public double[] getRandomCartTranslation() {
917 double[] coords = {random(), random(), random()};
918 toCartesianCoordinates(coords, coords);
919 return coords;
920 }
921
922
923
924
925
926
927 public double getSpecialPositionCutoff() {
928 return specialPositionCutoff;
929 }
930
931
932
933
934
935
936 public double getSpecialPositionCutoff2() {
937 return specialPositionCutoff2;
938 }
939
940
941
942
943
944
945 public void setSpecialPositionCutoff(double cutoff) {
946 specialPositionCutoff = cutoff;
947 specialPositionCutoff2 = cutoff * cutoff;
948 }
949
950
951
952
953
954
955
956 public void getTransformationOperator(SymOp symOp, double[][] rotmat) {
957 double[][] rot = symOp.rot;
958 for (int i = 0; i < 3; i++) {
959 for (int j = 0; j < 3; j++) {
960 rotmat[i][j] = 0.0;
961 for (int k = 0; k < 3; k++) {
962 for (int l = 0; l < 3; l++) {
963 rotmat[i][j] += Ai[k][i] * rot[k][l] * A[j][l];
964 }
965 }
966 }
967 }
968 }
969
970
971
972
973
974
975
976 public Crystal getUnitCell() {
977 return this;
978 }
979
980
981
982
983 @Override
984 public int hashCode() {
985 return Objects.hash(a, b, c, alpha, beta, gamma, spaceGroup.number);
986 }
987
988
989
990
991
992
993
994
995 public double image(double[] xyz, double[] xyz2) {
996 double dx = xyz[0] - xyz2[0];
997 double dy = xyz[1] - xyz2[1];
998 double dz = xyz[2] - xyz2[2];
999 return image(dx, dy, dz);
1000 }
1001
1002
1003
1004
1005
1006
1007
1008 public double image(final double[] xyz) {
1009 double x = xyz[0];
1010 double y = xyz[1];
1011 double z = xyz[2];
1012 if (aperiodic) {
1013 return x * x + y * y + z * z;
1014 }
1015 double xf = x * A00 + y * A10 + z * A20;
1016 double yf = x * A01 + y * A11 + z * A21;
1017 double zf = x * A02 + y * A12 + z * A22;
1018 xf = floor(abs(xf) + 0.5) * signum(-xf) + xf;
1019 yf = floor(abs(yf) + 0.5) * signum(-yf) + yf;
1020 zf = floor(abs(zf) + 0.5) * signum(-zf) + zf;
1021 x = xf * Ai00 + yf * Ai10 + zf * Ai20;
1022 y = xf * Ai01 + yf * Ai11 + zf * Ai21;
1023 z = xf * Ai02 + yf * Ai12 + zf * Ai22;
1024 xyz[0] = x;
1025 xyz[1] = y;
1026 xyz[2] = z;
1027 return x * x + y * y + z * z;
1028 }
1029
1030
1031
1032
1033
1034
1035
1036
1037
1038 public double image(double dx, double dy, double dz) {
1039 if (aperiodic) {
1040 return dx * dx + dy * dy + dz * dz;
1041 }
1042 double xf = dx * A00 + dy * A10 + dz * A20;
1043 double yf = dx * A01 + dy * A11 + dz * A21;
1044 double zf = dx * A02 + dy * A12 + dz * A22;
1045 xf = floor(abs(xf) + 0.5) * signum(-xf) + xf;
1046 yf = floor(abs(yf) + 0.5) * signum(-yf) + yf;
1047 zf = floor(abs(zf) + 0.5) * signum(-zf) + zf;
1048 dx = xf * Ai00 + yf * Ai10 + zf * Ai20;
1049 dy = xf * Ai01 + yf * Ai11 + zf * Ai21;
1050 dz = xf * Ai02 + yf * Ai12 + zf * Ai22;
1051 return dx * dx + dy * dy + dz * dz;
1052 }
1053
1054
1055
1056
1057
1058
1059
1060
1061 public double minDistOverSymOps(double[] xyzA, double[] xyzB) {
1062 double dist = 0;
1063 for (int i = 0; i < 3; i++) {
1064 double dx = xyzA[i] - xyzB[i];
1065 dist += (dx * dx);
1066 }
1067 var symB = new double[3];
1068 for (SymOp symOp : spaceGroup.symOps) {
1069 applySymOp(xyzB, symB, symOp);
1070 for (int i = 0; i < 3; i++) {
1071 symB[i] -= xyzA[i];
1072 }
1073 double d = image(symB);
1074 dist = min(d, dist);
1075 }
1076 return sqrt(dist);
1077 }
1078
1079
1080
1081
1082
1083
1084
1085 public boolean perturbCellVectors(double[][] dStrain) {
1086
1087 double[][] AA = new double[3][3];
1088 double[][] newAi = new double[3][3];
1089
1090 for (int i = 0; i < 3; i++) {
1091 for (int j = 0; j < 3; j++) {
1092 AA[i][j] = getUnitCell().Ai[i][j];
1093 }
1094 }
1095
1096 for (int i = 0; i < 3; i++) {
1097 for (int j = 0; j < 3; j++) {
1098 for (int k = 0; k < 3; k++) {
1099 double Kronecker = 0.0;
1100 if (j == k) {
1101 Kronecker = 1.0;
1102 }
1103
1104 newAi[i][j] += (Kronecker + dStrain[j][k]) * AA[i][k];
1105 }
1106 }
1107 }
1108
1109 return setCellVectors(newAi);
1110 }
1111
1112 public boolean randomParameters(double dens, double mass) {
1113 double[] params = latticeSystem.resetUnitCellParams();
1114 boolean succeed = changeUnitCellParameters(params[0], params[1], params[2], params[3], params[4], params[5]);
1115 if (succeed) {
1116 setDensity(dens, mass);
1117 }
1118 return succeed;
1119 }
1120
1121
1122
1123
1124
1125
1126 public void setAperiodic(boolean aperiodic) {
1127 this.aperiodic = aperiodic;
1128 }
1129
1130 public double[] getCellParametersFromVectors(double[][] cellVectors) {
1131
1132 double aa = length(cellVectors[0]);
1133 double bb = length(cellVectors[1]);
1134 double cc = length(cellVectors[2]);
1135
1136
1137 double aalpha = toDegrees(acos(dot(cellVectors[1], cellVectors[2]) / (bb * cc)));
1138 double bbeta = toDegrees(acos(dot(cellVectors[0], cellVectors[2]) / (aa * cc)));
1139 double ggamma = toDegrees(acos(dot(cellVectors[0], cellVectors[1]) / (aa * bb)));
1140
1141
1142 double[] params = new double[6];
1143 params[0] = aa;
1144 params[1] = bb;
1145 params[2] = cc;
1146 params[3] = aalpha;
1147 params[4] = bbeta;
1148 params[5] = ggamma;
1149 return params;
1150 }
1151
1152
1153
1154
1155
1156
1157
1158 public boolean setCellVectors(double[][] cellVectors) {
1159 double[] p = getCellParametersFromVectors(cellVectors);
1160 return changeUnitCellParameters(p[0], p[1], p[2], p[3], p[4], p[5]);
1161 }
1162
1163
1164
1165
1166
1167
1168
1169
1170 public boolean setCellVectorsAndVolume(double[][] cellVectors, double targetAUVolume) {
1171 if (setCellVectors(cellVectors)) {
1172 double currentAUVolume = volume / getNumSymOps();
1173 double scale = cbrt(targetAUVolume / currentAUVolume);
1174 return changeUnitCellParameters(scale * a, scale * b, scale * c, alpha, beta, gamma);
1175 } else {
1176 return false;
1177 }
1178 }
1179
1180 public void setDensity(double dens, double mass) {
1181 double currentDensity = getDensity(mass);
1182 double scale = cbrt(currentDensity / dens);
1183 changeUnitCellParameters(a * scale, b * scale, c * scale, alpha, beta, gamma);
1184 currentDensity = getDensity(mass);
1185 if (logger.isLoggable(Level.FINE)) {
1186 logger.fine(format(" Updated density %6.3f (g/cc) with unit cell %s.", currentDensity, toShortString()));
1187 }
1188 }
1189
1190
1191
1192
1193
1194
1195 public String toCRYST1() {
1196 return format(
1197 "CRYST1%9.3f%9.3f%9.3f%7.2f%7.2f%7.2f %10s\n",
1198 a, b, c, alpha, beta, gamma, padRight(spaceGroup.pdbName, 10));
1199 }
1200
1201
1202
1203
1204
1205
1206
1207
1208
1209
1210
1211
1212 public void toCartesianCoordinates(
1213 int n, double[] xf, double[] yf, double[] zf, double[] x, double[] y, double[] z) {
1214 for (int i = 0; i < n; i++) {
1215 double xi = xf[i];
1216 double yi = yf[i];
1217 double zi = zf[i];
1218 x[i] = xi * Ai00 + yi * Ai10 + zi * Ai20;
1219 y[i] = xi * Ai01 + yi * Ai11 + zi * Ai21;
1220 z[i] = xi * Ai02 + yi * Ai12 + zi * Ai22;
1221 }
1222 }
1223
1224
1225
1226
1227
1228
1229
1230
1231 public void toCartesianCoordinates(int n, double[] frac, double[] cart) {
1232 int i3 = 0;
1233 for (int i = 0; i < n; i++) {
1234
1235 int iX = i3 + XX;
1236 int iY = i3 + YY;
1237 int iZ = i3 + ZZ;
1238 i3 += 3;
1239 double xf = frac[iX];
1240 double yf = frac[iY];
1241 double zf = frac[iZ];
1242 cart[iX] = xf * Ai00 + yf * Ai10 + zf * Ai20;
1243 cart[iY] = xf * Ai01 + yf * Ai11 + zf * Ai21;
1244 cart[iZ] = xf * Ai02 + yf * Ai12 + zf * Ai22;
1245 }
1246 }
1247
1248
1249
1250
1251
1252
1253
1254 public void toCartesianCoordinates(double[] xf, double[] x) {
1255 double fx = xf[0];
1256 double fy = xf[1];
1257 double fz = xf[2];
1258 x[0] = fx * Ai00 + fy * Ai10 + fz * Ai20;
1259 x[1] = fx * Ai01 + fy * Ai11 + fz * Ai21;
1260 x[2] = fx * Ai02 + fy * Ai12 + fz * Ai22;
1261 }
1262
1263
1264
1265
1266
1267
1268
1269
1270
1271
1272
1273
1274 public void toFractionalCoordinates(
1275 int n, double[] x, double[] y, double[] z, double[] xf, double[] yf, double[] zf) {
1276 for (int i = 0; i < n; i++) {
1277 double xc = x[i];
1278 double yc = y[i];
1279 double zc = z[i];
1280 xf[i] = xc * A00 + yc * A10 + zc * A20;
1281 yf[i] = xc * A01 + yc * A11 + zc * A21;
1282 zf[i] = xc * A02 + yc * A12 + zc * A22;
1283 }
1284 }
1285
1286
1287
1288
1289
1290
1291
1292
1293 public void toFractionalCoordinates(int n, double[] cart, double[] frac) {
1294 int i3 = 0;
1295 for (int i = 0; i < n; i++) {
1296
1297 int iX = i3 + XX;
1298 int iY = i3 + YY;
1299 int iZ = i3 + ZZ;
1300 i3 += 3;
1301 double xc = cart[iX];
1302 double yc = cart[iY];
1303 double zc = cart[iZ];
1304 frac[iX] = xc * A00 + yc * A10 + zc * A20;
1305 frac[iY] = xc * A01 + yc * A11 + zc * A21;
1306 frac[iZ] = xc * A02 + yc * A12 + zc * A22;
1307 }
1308 }
1309
1310
1311
1312
1313
1314
1315
1316 public void toFractionalCoordinates(double[] x, double[] xf) {
1317 double xc = x[0];
1318 double yc = x[1];
1319 double zc = x[2];
1320 xf[0] = xc * A00 + yc * A10 + zc * A20;
1321 xf[1] = xc * A01 + yc * A11 + zc * A21;
1322 xf[2] = xc * A02 + yc * A12 + zc * A22;
1323 }
1324
1325
1326
1327
1328
1329
1330
1331 public void toPrimaryCell(double[] in, double[] out) {
1332 toFractionalCoordinates(in, out);
1333 out[0] = mod(out[0], 1.0);
1334 out[1] = mod(out[1], 1.0);
1335 out[2] = mod(out[2], 1.0);
1336 toCartesianCoordinates(out, out);
1337 }
1338
1339
1340
1341
1342
1343
1344 public String toShortString() {
1345 return format("%6.2f %6.2f %6.2f %6.2f %6.2f %6.2f", a, b, c, alpha, beta, gamma);
1346 }
1347
1348
1349
1350
1351 @Override
1352 public String toString() {
1353 StringBuilder sb = new StringBuilder("\n Unit Cell\n");
1354 sb.append(
1355 format(" A-axis: %8.3f (%8.3f, %8.3f, %8.3f)\n", a, Ai00, Ai01,
1356 Ai02));
1357 sb.append(
1358 format(" B-axis: %8.3f (%8.3f, %8.3f, %8.3f)\n", b, Ai10, Ai11,
1359 Ai12));
1360 sb.append(
1361 format(" C-axis: %8.3f (%8.3f, %8.3f, %8.3f)\n", c, Ai20, Ai21,
1362 Ai22));
1363 sb.append(format(" Alpha: %8.3f\n", alpha));
1364 sb.append(format(" Beta: %8.3f\n", beta));
1365 sb.append(format(" Gamma: %8.3f\n", gamma));
1366 sb.append(" Space group\n");
1367 sb.append(format(" Number: %3d\n", spaceGroup.number));
1368 sb.append(format(" Symbol: %8s\n", spaceGroup.shortName));
1369 sb.append(
1370 format(" Number of Symmetry Operators: %3d", spaceGroup.getNumberOfSymOps()));
1371 return sb.toString();
1372 }
1373
1374
1375
1376
1377 public final void updateCrystal() {
1378
1379 double cos_alpha;
1380 double sin_beta;
1381 double cos_beta;
1382 double sin_gamma;
1383 double cos_gamma;
1384 double beta_term;
1385 double gamma_term;
1386
1387 switch (crystalSystem) {
1388 default -> {
1389
1390 cos_alpha = 0.0;
1391 cos_beta = 0.0;
1392 sin_gamma = 1.0;
1393 cos_gamma = 0.0;
1394 beta_term = 0.0;
1395 gamma_term = 1.0;
1396 volume = a * b * c;
1397 dVdA = b * c;
1398 dVdB = a * c;
1399 dVdC = a * b;
1400 dVdAlpha = 0.0;
1401 dVdBeta = 0.0;
1402 dVdGamma = 0.0;
1403 }
1404 case MONOCLINIC -> {
1405 cos_alpha = 0.0;
1406 sin_beta = sin(toRadians(beta));
1407 cos_beta = cos(toRadians(beta));
1408 sin_gamma = 1.0;
1409 cos_gamma = 0.0;
1410 beta_term = 0.0;
1411 gamma_term = sin_beta;
1412 volume = sin_beta * a * b * c;
1413 dVdA = sin_beta * b * c;
1414 dVdB = sin_beta * a * c;
1415 dVdC = sin_beta * a * b;
1416 dVdAlpha = 0.0;
1417 dVdBeta = cos_beta * a * b * c;
1418 dVdGamma = 0.0;
1419 }
1420 case HEXAGONAL, TRICLINIC, TRIGONAL -> {
1421 double sin_alpha = sin(toRadians(alpha));
1422 cos_alpha = cos(toRadians(alpha));
1423 sin_beta = sin(toRadians(beta));
1424 cos_beta = cos(toRadians(beta));
1425 sin_gamma = sin(toRadians(gamma));
1426 cos_gamma = cos(toRadians(gamma));
1427 beta_term = (cos_alpha - cos_beta * cos_gamma) / sin_gamma;
1428 gamma_term = sqrt(sin_beta * sin_beta - beta_term * beta_term);
1429 volume = sin_gamma * gamma_term * a * b * c;
1430 dVdA = sin_gamma * gamma_term * b * c;
1431 dVdB = sin_gamma * gamma_term * a * c;
1432 dVdC = sin_gamma * gamma_term * a * b;
1433 double dbeta = 2.0 * sin_beta * cos_beta
1434 - (2.0 * (cos_alpha - cos_beta * cos_gamma) * sin_beta * cos_gamma)
1435 / (sin_gamma * sin_gamma);
1436 double dgamma1 = -2.0 * (cos_alpha - cos_beta * cos_gamma) * cos_beta / sin_gamma;
1437 double dgamma2 = cos_alpha - cos_beta * cos_gamma;
1438 dgamma2 *= dgamma2 * 2.0 * cos_gamma / (sin_gamma * sin_gamma * sin_gamma);
1439 dVdAlpha = (cos_alpha - cos_beta * cos_gamma) * sin_alpha
1440 / (sin_gamma * gamma_term) * a * b * c;
1441 dVdBeta = 0.5 * sin_gamma * dbeta / gamma_term * a * b * c;
1442 dVdGamma = (cos_gamma * gamma_term + 0.5 * sin_gamma * (dgamma1 + dgamma2) / gamma_term)
1443 * a * b * c;
1444 }
1445 }
1446
1447 G[0][0] = a * a;
1448 G[0][1] = a * b * cos_gamma;
1449 G[0][2] = a * c * cos_beta;
1450 G[1][0] = G[0][1];
1451 G[1][1] = b * b;
1452 G[1][2] = b * c * cos_alpha;
1453 G[2][0] = G[0][2];
1454 G[2][1] = G[1][2];
1455 G[2][2] = c * c;
1456
1457
1458 RealMatrix m = new Array2DRowRealMatrix(G, true);
1459 m = new LUDecomposition(m).getSolver().getInverse();
1460 Gstar = m.getData();
1461
1462
1463 Ai[0][0] = a;
1464 Ai[0][1] = 0.0;
1465 Ai[0][2] = 0.0;
1466
1467 Ai[1][0] = b * cos_gamma;
1468 Ai[1][1] = b * sin_gamma;
1469 Ai[1][2] = 0.0;
1470
1471 Ai[2][0] = c * cos_beta;
1472 Ai[2][1] = c * beta_term;
1473 Ai[2][2] = c * gamma_term;
1474
1475 Ai00 = Ai[0][0];
1476 Ai01 = Ai[0][1];
1477 Ai02 = Ai[0][2];
1478 Ai10 = Ai[1][0];
1479 Ai11 = Ai[1][1];
1480 Ai12 = Ai[1][2];
1481 Ai20 = Ai[2][0];
1482 Ai21 = Ai[2][1];
1483 Ai22 = Ai[2][2];
1484
1485
1486 m = new Array2DRowRealMatrix(Ai, true);
1487 m = new LUDecomposition(m).getSolver().getInverse();
1488 A = m.getData();
1489
1490
1491 A00 = A[0][0];
1492 A10 = A[1][0];
1493 A20 = A[2][0];
1494 A01 = A[0][1];
1495 A11 = A[1][1];
1496 A21 = A[2][1];
1497 A02 = A[0][2];
1498 A12 = A[1][2];
1499 A22 = A[2][2];
1500
1501
1502 double aStar = 1.0 / sqrt(A00 * A00 + A10 * A10 + A20 * A20);
1503 double bStar = 1.0 / sqrt(A01 * A01 + A11 * A11 + A21 * A21);
1504 double cStar = 1.0 / sqrt(A02 * A02 + A12 * A12 + A22 * A22);
1505 if (logger.isLoggable(Level.FINEST)) {
1506 logger.finest(
1507 format(" Reciprocal Lattice Lengths: (%8.3f, %8.3f, %8.3f)", aStar, bStar, cStar));
1508 }
1509
1510
1511 interfacialRadiusA = (Ai00 * A00 + Ai01 * A10 + Ai02 * A20) * aStar;
1512 interfacialRadiusB = (Ai10 * A01 + Ai11 * A11 + Ai12 * A21) * bStar;
1513 interfacialRadiusC = (Ai20 * A02 + Ai21 * A12 + Ai22 * A22) * cStar;
1514
1515
1516 interfacialRadiusA /= 2.0;
1517 interfacialRadiusB /= 2.0;
1518 interfacialRadiusC /= 2.0;
1519
1520 if (logger.isLoggable(Level.FINEST)) {
1521 logger.finest(
1522 format(
1523 " Interfacial radii: (%8.3f, %8.3f, %8.3f)",
1524 interfacialRadiusA, interfacialRadiusB, interfacialRadiusC));
1525 }
1526
1527 List<SymOp> symOps = spaceGroup.symOps;
1528 int nSymm = symOps.size();
1529 if (symOpsCartesian == null) {
1530 symOpsCartesian = new ArrayList<>(nSymm);
1531 } else {
1532 symOpsCartesian.clear();
1533 }
1534
1535 RealMatrix toFrac = new Array2DRowRealMatrix(A);
1536 RealMatrix toCart = new Array2DRowRealMatrix(Ai);
1537 for (SymOp symOp : symOps) {
1538 m = new Array2DRowRealMatrix(symOp.rot);
1539
1540 RealMatrix rotMat = m.preMultiply(toCart.transpose()).multiply(toFrac.transpose());
1541
1542 double[] tr = toCart.preMultiply(symOp.tr);
1543 symOpsCartesian.add(new SymOp(rotMat.getData(), tr));
1544 }
1545 }
1546 }