View Javadoc
1   // ******************************************************************************
2   //
3   // Title:       Force Field X.
4   // Description: Force Field X - Software for Molecular Biophysics.
5   // Copyright:   Copyright (c) Michael J. Schnieders 2001-2023.
6   //
7   // This file is part of Force Field X.
8   //
9   // Force Field X is free software; you can redistribute it and/or modify it
10  // under the terms of the GNU General Public License version 3 as published by
11  // the Free Software Foundation.
12  //
13  // Force Field X is distributed in the hope that it will be useful, but WITHOUT
14  // ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
15  // FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
16  // details.
17  //
18  // You should have received a copy of the GNU General Public License along with
19  // Force Field X; if not, write to the Free Software Foundation, Inc., 59 Temple
20  // Place, Suite 330, Boston, MA 02111-1307 USA
21  //
22  // Linking this library statically or dynamically with other modules is making a
23  // combined work based on this library. Thus, the terms and conditions of the
24  // GNU General Public License cover the whole combination.
25  //
26  // As a special exception, the copyright holders of this library give you
27  // permission to link this library with independent modules to produce an
28  // executable, regardless of the license terms of these independent modules, and
29  // to copy and distribute the resulting executable under terms of your choice,
30  // provided that you also meet, for each linked independent module, the terms
31  // and conditions of the license of that module. An independent module is a
32  // module which is not derived from or based on this library. If you modify this
33  // library, you may extend this exception to your version of the library, but
34  // you are not obligated to do so. If you do not wish to do so, delete this
35  // exception statement from your version.
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   * The Crystal class encapsulates the lattice parameters and space group that describe the geometry
76   * and symmetry of a crystal. Methods are available to apply the minimum image convention and space
77   * group symmetry operators.
78   *
79   * @author Michael J. Schnieders
80   * @see ReplicatesCrystal
81   * @since 1.0
82   */
83  public class Crystal {
84  
85    private static final Logger logger = Logger.getLogger(Crystal.class.getName());
86  
87    /** The space group of the crystal. */
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    /** Length of the cell edge in the direction of the <b>a</b> basis vector. */
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    /** Length of the cell edge in the direction of the <b>b</b> basis vector. */
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   /** Length of the cell edge in the direction of the <b>c</b> basis vector. */
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   /** The interaxial lattice angle between <b>b</b> and <b>c</b>. */
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   /** The interaxial lattice angle between <b>a</b> and <b>c</b>. */
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   /** The interaxial lattice angle between <b>a</b> and <b>b</b>. */
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   /** A mask equal to 0 for X-coordinates. */
123   private static final int XX = 0;
124   /** A mask equal to 1 for Y-coordinates. */
125   private static final int YY = 1;
126   /** A mask equal to 2 for Z-coordinates. */
127   private static final int ZZ = 2;
128 
129   /**
130    * Matrix to convert from fractional to Cartesian coordinates.
131    * <br>a-axis vector is the first row of A^(-1).
132    * <br>b-axis vector is the second row of A^(-1).
133    * <br>c-axis vector is the third row of A^(-1).
134    */
135   public final double[][] Ai = new double[3][3];
136   /** The direct space metric matrix. */
137   public final double[][] G = new double[3][3];
138   /** Reference to the space group crystal system. */
139   private final CrystalSystem crystalSystem;
140   /** Reference to the space group lattice system. */
141   private final LatticeSystem latticeSystem;
142   /** The crystal unit cell volume. */
143   public double volume;
144   /** Matrix to convert from Cartesian to fractional coordinates. */
145   public double[][] A;
146   /** Entry in the A matrix. */
147   public double A00, A01, A02, A10, A11, A12, A20, A21, A22;
148   /** Interfacial radius in the direction of the A-axis. */
149   public double interfacialRadiusA;
150   /** Interfacial radius in the direction of the B-axis. */
151   public double interfacialRadiusB;
152   /** Interfacial radius in the direction of the C-axis. */
153   public double interfacialRadiusC;
154   /**
155    * Anisotropic bulk solvent B-factor scaling (0 or 1 for each component).
156    */
157   public int[] scaleB = new int[6];
158   /**
159    * Number of bulk solvent B-factor components.
160    */
161   public int scaleN;
162   /** Entry in the Ai matrix. */
163   public double Ai00, Ai01, Ai02, Ai10, Ai11, Ai12, Ai20, Ai21, Ai22;
164   /** Change in the volume with respect to a. */
165   public double dVdA;
166   /** Change in the volume with respect to b. */
167   public double dVdB;
168   /** Change in the volume with respect to c. */
169   public double dVdC;
170   /**
171    * Change in the volume with respect to alpha (in Radians). This is set to zero if alpha is fixed.
172    */
173   public double dVdAlpha;
174   /**
175    * Change in the volume with respect to beta (in Radians). This is set to zero if beta is fixed.
176    */
177   public double dVdBeta;
178   /**
179    * Change in the volume with respect to gamma (in Radians). This is set to zero if gamma is fixed.
180    */
181   public double dVdGamma;
182   /**
183    * For some finite-difference calculations, it's currently necessary to remove lattice system
184    * restrictions.
185    */
186   boolean checkRestrictions = true;
187   /**
188    * An atom and one of its symmetry copies within the specialPositionCutoff should be flagged to be
189    * at a special position.
190    */
191   private double specialPositionCutoff = 0.3;
192   /** Copy of symmetry operators in Cartesian coordinates. */
193   private List<SymOp> symOpsCartesian;
194   /** The reciprocal space metric matrix. */
195   private double[][] Gstar;
196   /**
197    * SpecialPositionCutoff squared.
198    */
199   private double specialPositionCutoff2 = specialPositionCutoff * specialPositionCutoff;
200   /**
201    * Flag to indicate an aperiodic system.
202    */
203   private boolean aperiodic;
204 
205   /**
206    * The Crystal class encapsulates the lattice parameters and space group. Methods are available to
207    * apply the minimum image convention and to apply space group operators.
208    *
209    * @param a The a-axis length.
210    * @param b The b-axis length.
211    * @param c The c-axis length.
212    * @param alpha The alpha angle.
213    * @param beta The beta angle.
214    * @param gamma The gamma angle.
215    * @param sgNumber The space group number.
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    * The Crystal class encapsulates the lattice parameters and space group. Methods are available to
223    * apply the minimum image convention and to apply space group operators.
224    *
225    * @param a The a-axis length.
226    * @param b The b-axis length.
227    * @param c The c-axis length.
228    * @param alpha The alpha angle.
229    * @param beta The beta angle.
230    * @param gamma The gamma angle.
231    * @param sg The space group symbol.
232    */
233   public Crystal(double a, double b, double c, double alpha, double beta, double gamma, String sg) {
234     // Crystal SpaceGroup and LatticeSystem are final variables. Temp variable to delay assigning.
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       // Invalid parameters... Start error/warning log and try to fix.
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         // Try to convert between hexagonal and rhombohedral lattices to fix crystal.
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           // Successfully converted space group between hexagonal and rhombohedral.
286           logger.info(sb.toString());
287         }
288       } else {
289         // Invalid lattice parameters. Update Crystal as much as possible, then print error message.
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       // Valid parameters, update crystal and continue.
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         // determine unique axis
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           // rhombohedral
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    * checkProperties
405    *
406    * @param properties a {@link org.apache.commons.configuration2.CompositeConfiguration}
407    *     object.
408    * @return a {@link ffx.crystal.Crystal} object.
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     // check the space group name is valid
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    * invressq
440    *
441    * @param hkl a {@link ffx.crystal.HKL} object.
442    * @return a double.
443    */
444   public double invressq(HKL hkl) {
445     return hkl.quadForm(Gstar);
446   }
447 
448   /**
449    * res
450    *
451    * @param hkl a {@link ffx.crystal.HKL} object.
452    * @return a double.
453    */
454   public double res(HKL hkl) {
455     return 1.0 / sqrt(hkl.quadForm(Gstar));
456   }
457 
458   /**
459    * aperiodic
460    *
461    * @return a boolean.
462    */
463   public boolean aperiodic() {
464     return aperiodic;
465   }
466 
467   /**
468    * Apply a fractional symmetry operator to an array of Cartesian coordinates.
469    * <p>
470    * If the arrays x, y or z are null or not of length n, an exception is thrown.
471    * <p>
472    * If mateX, mateY or mateZ are null or not of length n, an exception is thrown.
473    *
474    * @param n Number of atoms.
475    * @param x Input fractional x-coordinates.
476    * @param y Input fractional y-coordinates.
477    * @param z Input fractional z-coordinates.
478    * @param mateX Output fractional x-coordinates.
479    * @param mateY Output fractional y-coordinates.
480    * @param mateZ Output fractional z-coordinates.
481    * @param symOp The fractional symmetry operator.
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       // Convert to fractional coordinates.
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       // Apply Symmetry Operator.
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       // Convert back to Cartesian coordinates.
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    * Apply a fractional symmetry operator to one set of cartesian coordinates.
533    *
534    * @param xyz Input cartesian coordinates.
535    * @param mate Symmetry mate cartesian coordinates.
536    * @param symOp The fractional symmetry operator.
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       // Convert to fractional coordinates.
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       // Apply Symmetry Operator.
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       // Convert back to Cartesian coordinates.
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    * Apply a fractional symmetry operator to one set of cartesian coordinates.
581    *
582    * @param xyz Input cartesian coordinates.
583    * @param mate Symmetry mate cartesian coordinates.
584    * @param symOp The fractional symmetry operator.
585    */
586   public void applySymRot(double[] xyz, double[] mate, SymOp symOp) {
587     double[][] rot = symOp.rot;
588     // Convert to fractional coordinates.
589     double xc = xyz[0];
590     double yc = xyz[1];
591     double zc = xyz[2];
592     // Convert to fractional coordinates.
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     // Apply Symmetry Operator.
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     // Convert back to Cartesian coordinates.
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    * Apply the transpose of a symmetry rotation to an array of Cartesian coordinates.
610    * <p>
611    * If the arrays x, y or z are null or not of length n, an exception is thrown.
612    * <p>
613    * If mateX, mateY or mateZ are null or not of length n, an exception is thrown.
614    *
615    * @param n Number of atoms.
616    * @param x Input x-coordinates.
617    * @param y Input y-coordinates.
618    * @param z Input z-coordinates.
619    * @param mateX Output x-coordinates.
620    * @param mateY Output y-coordinates.
621    * @param mateZ Output z-coordinates.
622    * @param symOp The symmetry operator.
623    * @param rotmat an array of double.
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     // The transformation operator R = ToCart * Rot * ToFrac
643     getTransformationOperator(symOp, rotmat);
644 
645     for (int i = 0; i < n; i++) {
646       // Apply R^T (its transpose).
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    * averageTensor
658    *
659    * @param m an array of double.
660    * @param r an array of double.
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    * averageTensor
684    *
685    * @param v an array of double.
686    * @param r an array of double.
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    * This method should be called to update the unit cell parameters of a crystal. The proposed
710    * parameters will only be accepted if symmetry restrictions are satisfied. If so, all Crystal
711    * variables that depend on the unit cell parameters will be updated.
712    *
713    * @param a length of the a-axis.
714    * @param b length of the b-axis.
715    * @param c length of the c-axis.
716    * @param alpha Angle between b-axis and c-axis.
717    * @param beta Angle between a-axis and c-axis.
718    * @param gamma Angle between a-axis and b-axis.
719    * @return The method return true if the parameters are accepted, false otherwise.
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    * This method should be called to update the unit cell parameters of a crystal. The proposed
756    * parameters will only be accepted if symmetry restrictions are satisfied. If so, all Crystal
757    * variables that depend on the unit cell parameters will be updated.
758    * <p>
759    * If the new parameters are accepted, the target asymmetric unit volume is achieved by uniformly
760    * scaling all lattice lengths.
761    *
762    * @param a length of the a-axis.
763    * @param b length of the b-axis.
764    * @param c length of the c-axis.
765    * @param alpha Angle between b-axis and c-axis.
766    * @param beta Angle between a-axis and c-axis.
767    * @param gamma Angle between a-axis and b-axis.
768    * @param targetAUVolume Target asymmetric unit volume.
769    * @return The method return true if the parameters are accepted, false otherwise.
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    * Two crystals are equal only if all unit cell parameters are exactly the same.
783    *
784    * @param o the Crystal to compare to.
785    * @return true if all unit cell parameters are exactly the same.
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    * Compute the density of the system.
815    *
816    * @param mass The total mass of the asymmetric unit.
817    * @return The density (g/cc)
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    * Return the number of symmetry operators for this crystal.
826    *
827    * @return The number of symmetry operators.
828    */
829   public int getNumSymOps() {
830     return spaceGroup.getNumberOfSymOps();
831   }
832 
833   /**
834    * Create a random Cartesian translation vector.
835    *
836    * <p>First, a random fractional translation is created. Second, the random fractional operator is
837    * converted to Cartesian coordinates.
838    *
839    * @return A random Cartesian translation vector.
840    */
841   public double[] getRandomCartTranslation() {
842     double[] coords = {random(), random(), random()};
843     toCartesianCoordinates(coords, coords);
844     return coords;
845   }
846 
847   /**
848    * Getter for the field <code>specialPositionCutoff</code>.
849    *
850    * @return a double.
851    */
852   public double getSpecialPositionCutoff() {
853     return specialPositionCutoff;
854   }
855 
856   /**
857    * Getter for the field <code>specialPositionCutoff2</code>.
858    *
859    * @return a double.
860    */
861   public double getSpecialPositionCutoff2() {
862     return specialPositionCutoff2;
863   }
864 
865   /**
866    * Setter for the field <code>specialPositionCutoff</code>.
867    *
868    * @param cutoff a double.
869    */
870   public void setSpecialPositionCutoff(double cutoff) {
871     specialPositionCutoff = cutoff;
872     specialPositionCutoff2 = cutoff * cutoff;
873   }
874 
875   /**
876    * Compute the total transformation operator R = ToCart * Rot * ToFrac.
877    *
878    * @param symOp Symmetry operator to apply.
879    * @param rotmat Resulting transformation operator R.
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    * The ReplicatesCrystal over-rides this method to return the unit cell rather than the
897    * ReplicateCell.
898    *
899    * @return The unit cell Crystal instance.
900    */
901   public Crystal getUnitCell() {
902     return this;
903   }
904 
905   /** {@inheritDoc} */
906   @Override
907   public int hashCode() {
908     return Objects.hash(a, b, c, alpha, beta, gamma, spaceGroup.number);
909   }
910 
911   /**
912    * Apply the minimum image convention.
913    *
914    * @param xyz the coordinates of the first atom.
915    * @param xyz2 the coordinates of the second atom.
916    * @return the output distance squared.
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    * Apply the minimum image convention.
927    *
928    * @param xyz input distances that are over-written.
929    * @return the output distance squared.
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    * Apply the minimum image convention.
955    *
956    * @param dx x-distance
957    * @param dy y-distance
958    * @param dz z-distance
959    * @return the output distance squared.
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    * Minimum distance between two coordinates over all symmetry operators.
979    *
980    * @param xyzA Coordinate A
981    * @param xyzB Coordinate B
982    * @return Minimum distance in crystal
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    * Strain the unit cell vectors.
1004    *
1005    * @param dStrain a 3x3 matrix of unitless Strain percentages.
1006    * @return True if the perturbation of cell vectors succeeds.
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           // Aij' = Sum over K ( Kron_jk + dStrain_jk ) * Aij
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    * Is this a finite system without periodic boundary conditions.
1046    *
1047    * @param aperiodic a boolean.
1048    */
1049   public void setAperiodic(boolean aperiodic) {
1050     this.aperiodic = aperiodic;
1051   }
1052 
1053   public double[] getCellParametersFromVectors(double[][] cellVectors) {
1054     // Update a-, b-, and c-axis lengths.
1055     double aa = length(cellVectors[0]);
1056     double bb = length(cellVectors[1]);
1057     double cc = length(cellVectors[2]);
1058 
1059     // Update alpha, beta and gamma angles.
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     // Load and return new cell parameters.
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    * Set the unit cell vectors.
1077    *
1078    * @param cellVectors 3x3 matrix of cell vectors.
1079    * @return True if the perturbation of cell vectors succeeds.
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    * Set the unit cell vectors. Scale lattice lengths if necessary to hit the target volume.
1088    *
1089    * @param cellVectors 3x3 matrix of cell vectors.
1090    * @param targetAUVolume the target volume for the new cell Vectors.
1091    * @return True if the perturbation of cell vectors succeeds.
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    * Return a CRYST1 record useful for writing a PDB file.
1115    *
1116    * @return The CRYST1 record.
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    * toCartesianCoordinates
1126    *
1127    * @param n an int.
1128    * @param xf an array of double for fractional x-coordinates.
1129    * @param yf an array of double for fractional y-coordinates.
1130    * @param zf an array of double for fractional z-coordinates.
1131    * @param x an array of double for cartesian x-coordinates.
1132    * @param y an array of double for cartesian y-coordinates.
1133    * @param z an array of double for cartesian z-coordinates.
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    * toCartesianCoordinates
1149    *
1150    * @param n an int.
1151    * @param frac an array of double for fractional coordinates.
1152    * @param cart an array of double for cartesian coordinates.
1153    */
1154   public void toCartesianCoordinates(int n, double[] frac, double[] cart) {
1155     int i3 = 0;
1156     for (int i = 0; i < n; i++) {
1157       // Convert to cartesian coordinates.
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    * toCartesianCoordinates
1173    *
1174    * @param xf an array of double for fractional coordinate.
1175    * @param x an array of double for cartesian coordinate.
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    * toFractionalCoordinates
1188    *
1189    * @param n an int.
1190    * @param x an array of double for cartesian x-coordinates.
1191    * @param y an array of double for cartesian y-coordinates.
1192    * @param z an array of double for cartesian z-coordinates.
1193    * @param xf an array of double for fractional x-coordinates.
1194    * @param yf an array of double for fractional y-coordinates.
1195    * @param zf an array of double for fractional z-coordinates.
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    * toFractionalCoordinates
1211    *
1212    * @param n an int.
1213    * @param cart an array of double for cartesian coordinates.
1214    * @param frac an array of double for fractional coordinates.
1215    */
1216   public void toFractionalCoordinates(int n, double[] cart, double[] frac) {
1217     int i3 = 0;
1218     for (int i = 0; i < n; i++) {
1219       // Convert to fractional coordinates.
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    * toFractionalCoordinates
1235    *
1236    * @param x an array of double for cartesian coordinate.
1237    * @param xf an array of double for fractional coordinate.
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    * toPrimaryCell
1250    *
1251    * @param in an array of double.
1252    * @param out an array of double.
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    * A String containing the unit cell parameters.
1264    *
1265    * @return A string with the unit cell parameters.
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   /** {@inheritDoc} */
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   /** Update all Crystal variables that are a function of unit cell parameters. */
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         // CUBIC, ORTHORHOMBIC, TETRAGONAL
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     // invert G to yield Gstar
1377     RealMatrix m = new Array2DRowRealMatrix(G, true);
1378     m = new LUDecomposition(m).getSolver().getInverse();
1379     Gstar = m.getData();
1380 
1381     // a is the first row of A^(-1).
1382     Ai[0][0] = a;
1383     Ai[0][1] = 0.0;
1384     Ai[0][2] = 0.0;
1385     // b is the second row of A^(-1).
1386     Ai[1][0] = b * cos_gamma;
1387     Ai[1][1] = b * sin_gamma;
1388     Ai[1][2] = 0.0;
1389     // c is the third row of A^(-1).
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     // Invert A^-1 to get A
1405     m = new Array2DRowRealMatrix(Ai, true);
1406     m = new LUDecomposition(m).getSolver().getInverse();
1407     A = m.getData();
1408 
1409     // The columns of A are the reciprocal basis vectors
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     // Reciprocal basis vector lengths
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     // Interfacial diameters from the dot product of the real and reciprocal vectors
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     // Divide by 2 to get radii.
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       // rot_c = A^(-1).rot_f.A
1459       RealMatrix rotMat = m.preMultiply(toCart.transpose()).multiply(toFrac.transpose());
1460       // tr_c = tr_f.A^(-1)
1461       double[] tr = toCart.preMultiply(symOp.tr);
1462       symOpsCartesian.add(new SymOp(rotMat.getData(), tr));
1463     }
1464   }
1465 }