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