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.MatrixMath.mat3Inverse;
42  import static ffx.numerics.math.MatrixMath.mat4Mat4;
43  import static ffx.numerics.math.ScalarMath.mod;
44  import static java.lang.Double.parseDouble;
45  import static java.lang.String.format;
46  import static org.apache.commons.math3.util.FastMath.PI;
47  import static org.apache.commons.math3.util.FastMath.cos;
48  import static org.apache.commons.math3.util.FastMath.random;
49  import static org.apache.commons.math3.util.FastMath.rint;
50  import static org.apache.commons.math3.util.FastMath.sin;
51  import static org.apache.commons.math3.util.FastMath.sqrt;
52  
53  /**
54   * The SymOp class defines the rotation and translation of a single symmetry operator.
55   *
56   * @author Michael J. Schnieders
57   * @see SpaceGroup
58   * @since 1.0
59   */
60  public class SymOp {
61  
62    /** Constant <code>ZERO = 0.0</code> */
63    private static final double ZERO = 0.0;
64    /** Constant <code>ZERO_ROTATION = {{1.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0.0, 0.0, 1.0}}</code> */
65    public static final double[][] ZERO_ROTATION = {
66        {1.0, ZERO, ZERO},
67        {ZERO, 1.0, ZERO},
68        {ZERO, ZERO, 1.0}};
69    /** Constant <code>Tr_0_0_0={ZERO, ZERO, ZERO}</code> */
70    public static final double[] Tr_0_0_0 = {ZERO, ZERO, ZERO};
71    /** Constant <code>f12 = 1.0 / 2.0</code> */
72    private static final double f12 = 1.0 / 2.0;
73    /** Constant <code>Tr_12_0_12={f12, ZERO, f12}</code> */
74    static final double[] Tr_12_0_12 = {f12, ZERO, f12};
75    /** Constant <code>Tr_0_12_12={ZERO, f12, f12}</code> */
76    static final double[] Tr_0_12_12 = {ZERO, f12, f12};
77    /** Constant <code>Tr_12_12_0={f12, f12, ZERO}</code> */
78    static final double[] Tr_12_12_0 = {f12, f12, ZERO};
79    /** Constant <code>Tr_12_12_12={f12, f12, f12}</code> */
80    static final double[] Tr_12_12_12 = {f12, f12, f12};
81    /** Constant <code>Tr_0_12_0={ZERO, f12, ZERO}</code> */
82    static final double[] Tr_0_12_0 = {ZERO, f12, ZERO};
83    /** Constant <code>Tr_12_0_0={f12, ZERO, ZERO}</code> */
84    static final double[] Tr_12_0_0 = {f12, ZERO, ZERO};
85    /** Constant <code>Tr_0_0_12={ZERO, ZERO, f12}</code> */
86    static final double[] Tr_0_0_12 = {ZERO, ZERO, f12};
87    /** Constant <code>f13 = 1.0 / 3.0</code> */
88    private static final double f13 = 1.0 / 3.0;
89    /** Constant <code>Tr_0_0_13={ZERO, ZERO, f13}</code> */
90    static final double[] Tr_0_0_13 = {ZERO, ZERO, f13};
91    /** Constant <code>f23 = 2.0 / 3.0</code> */
92    private static final double f23 = 2.0 / 3.0;
93    /** Constant <code>Tr_23_13_13={f23, f13, f13}</code> */
94    static final double[] Tr_23_13_13 = {f23, f13, f13};
95    /** Constant <code>Tr_13_23_23={f13, f23, f23}</code> */
96    static final double[] Tr_13_23_23 = {f13, f23, f23};
97    /** Constant <code>Tr_0_0_23={ZERO, ZERO, f23}</code> */
98    static final double[] Tr_0_0_23 = {ZERO, ZERO, f23};
99    /** Constant <code>f14 = 1.0 / 4.0</code> */
100   private static final double f14 = 1.0 / 4.0;
101   /** Constant <code>Tr_12_0_14={f12, ZERO, f14}</code> */
102   static final double[] Tr_12_0_14 = {f12, ZERO, f14};
103   /** Constant <code>Tr_0_12_14={ZERO, f12, f14}</code> */
104   static final double[] Tr_0_12_14 = {ZERO, f12, f14};
105   /** Constant <code>Tr_14_14_14={f14, f14, f14}</code> */
106   static final double[] Tr_14_14_14 = {f14, f14, f14};
107   /** Constant <code>Tr_12_12_14={f12, f12, f14}</code> */
108   static final double[] Tr_12_12_14 = {f12, f12, f14};
109   /** Constant <code>Tr_0_0_14={ZERO, ZERO, f14}</code> */
110   static final double[] Tr_0_0_14 = {ZERO, ZERO, f14};
111   /** Constant <code>f34 = 3.0 / 4.0</code> */
112   private static final double f34 = 3.0 / 4.0;
113   /** Constant <code>Tr_0_0_34={ZERO, ZERO, f34}</code> */
114   static final double[] Tr_0_0_34 = {ZERO, ZERO, f34};
115   /** Constant <code>Tr_12_0_34={f12, ZERO, f34}</code> */
116   static final double[] Tr_12_0_34 = {f12, ZERO, f34};
117   /** Constant <code>Tr_0_12_34={ZERO, f12, f34}</code> */
118   static final double[] Tr_0_12_34 = {ZERO, f12, f34};
119   /** Constant <code>Tr_34_14_14={f34, f14, f14}</code> */
120   static final double[] Tr_34_14_14 = {f34, f14, f14};
121   /** Constant <code>Tr_14_14_34={f14, f14, f34}</code> */
122   static final double[] Tr_14_14_34 = {f14, f14, f34};
123   /** Constant <code>Tr_14_34_14={f14, f34, f14}</code> */
124   static final double[] Tr_14_34_14 = {f14, f34, f14};
125   /** Constant <code>Tr_12_12_34={f12, f12, f34}</code> */
126   static final double[] Tr_12_12_34 = {f12, f12, f34};
127   /** Constant <code>Tr_14_34_34={f14, f34, f34}</code> */
128   static final double[] Tr_14_34_34 = {f14, f34, f34};
129   /** Constant <code>Tr_34_34_14={f34, f34, f14}</code> */
130   static final double[] Tr_34_34_14 = {f34, f34, f14};
131   /** Constant <code>Tr_34_34_34={f34, f34, f34}</code> */
132   static final double[] Tr_34_34_34 = {f34, f34, f34};
133   /** Constant <code>Tr_34_14_34={f34, f14, f34}</code> */
134   static final double[] Tr_34_14_34 = {f34, f14, f34};
135   /** Constant <code>f16 = 1.0 / 6.0</code> */
136   private static final double f16 = 1.0 / 6.0;
137   /** Constant <code>Tr_13_23_16={f13, f23, f16}</code> */
138   static final double[] Tr_13_23_16 = {f13, f23, f16};
139   /** Constant <code>Tr_0_0_16={ZERO, ZERO, f16}</code> */
140   static final double[] Tr_0_0_16 = {ZERO, ZERO, f16};
141   /** Constant <code>f56 = 5.0 / 6.0</code> */
142   private static final double f56 = 5.0 / 6.0;
143   /** Constant <code>Tr_0_0_56={ZERO, ZERO, f56}</code> */
144   static final double[] Tr_0_0_56 = {ZERO, ZERO, f56};
145   /** Constant <code>Tr_23_13_56={f23, f13, f56}</code> */
146   static final double[] Tr_23_13_56 = {f23, f13, f56};
147   /** Constant <code>X = {1.0, ZERO, ZERO}</code> */
148   private static final double[] X = {1.0, ZERO, ZERO};
149   /** Constant <code>Y = {ZERO, 1.0, ZERO}</code> */
150   private static final double[] Y = {ZERO, 1.0, ZERO};
151   /** Constant <code>Z = {ZERO, ZERO, 1.0}</code> */
152   private static final double[] Z = {ZERO, ZERO, 1.0};
153   /** Constant <code>Rot_Y_Z_X={Y, Z, X}</code> */
154   static final double[][] Rot_Y_Z_X = {Y, Z, X};
155   /** Constant <code>Rot_X_Y_Z={X, Y, Z}</code> */
156   static final double[][] Rot_X_Y_Z = {X, Y, Z};
157   /** Constant <code>Rot_Z_X_Y={Z, X, Y}</code> */
158   static final double[][] Rot_Z_X_Y = {Z, X, Y};
159   /** Constant <code>Rot_X_Z_Y={X, Z, Y}</code> */
160   static final double[][] Rot_X_Z_Y = {X, Z, Y};
161   /** Constant <code>Rot_Z_Y_X={Z, Y, X}</code> */
162   static final double[][] Rot_Z_Y_X = {Z, Y, X};
163   /** Constant <code>Rot_Y_X_Z={Y, X, Z}</code> */
164   static final double[][] Rot_Y_X_Z = {Y, X, Z};
165   /** Constant <code>mX = {-1.0, ZERO, ZERO}</code> */
166   private static final double[] mX = {-1.0, ZERO, ZERO};
167   /** Constant <code>Rot_Y_mX_Z={Y, mX, Z}</code> */
168   static final double[][] Rot_Y_mX_Z = {Y, mX, Z};
169   /** Constant <code>Rot_mX_Z_Y={mX, Z, Y}</code> */
170   static final double[][] Rot_mX_Z_Y = {mX, Z, Y};
171   /** Constant <code>Rot_Y_Z_mX={Y, Z, mX}</code> */
172   static final double[][] Rot_Y_Z_mX = {Y, Z, mX};
173   /** Constant <code>Rot_mX_Y_Z={mX, Y, Z}</code> */
174   static final double[][] Rot_mX_Y_Z = {mX, Y, Z};
175   /** Constant <code>Rot_Z_Y_mX={Z, Y, mX}</code> */
176   static final double[][] Rot_Z_Y_mX = {Z, Y, mX};
177   /** Constant <code>Rot_Z_mX_Y={Z, mX, Y}</code> */
178   static final double[][] Rot_Z_mX_Y = {Z, mX, Y};
179   /** Constant <code>mY = {ZERO, -1.0, ZERO}</code> */
180   private static final double[] mY = {ZERO, -1.0, ZERO};
181   /** Constant <code>Rot_Z_mY_X={Z, mY, X}</code> */
182   static final double[][] Rot_Z_mY_X = {Z, mY, X};
183   /** Constant <code>Rot_X_Z_mY={X, Z, mY}</code> */
184   static final double[][] Rot_X_Z_mY = {X, Z, mY};
185   /** Constant <code>Rot_mY_X_Z={mY, X, Z}</code> */
186   static final double[][] Rot_mY_X_Z = {mY, X, Z};
187   /** Constant <code>Rot_mY_Z_mX={mY, Z, mX}</code> */
188   static final double[][] Rot_mY_Z_mX = {mY, Z, mX};
189   /** Constant <code>Rot_mY_Z_X={mY, Z, X}</code> */
190   static final double[][] Rot_mY_Z_X = {mY, Z, X};
191   /** Constant <code>Rot_Z_X_mY={Z, X, mY}</code> */
192   static final double[][] Rot_Z_X_mY = {Z, X, mY};
193   /** Constant <code>Rot_Z_mX_mY={Z, mX, mY}</code> */
194   static final double[][] Rot_Z_mX_mY = {Z, mX, mY};
195   /** Constant <code>Rot_mX_Z_mY={mX, Z, mY}</code> */
196   static final double[][] Rot_mX_Z_mY = {mX, Z, mY};
197   /** Constant <code>Rot_X_mY_Z={X, mY, Z}</code> */
198   static final double[][] Rot_X_mY_Z = {X, mY, Z};
199   /** Constant <code>Rot_mY_mX_Z={mY, mX, Z}</code> */
200   static final double[][] Rot_mY_mX_Z = {mY, mX, Z};
201   /** Constant <code>Rot_Z_mY_mX={Z, mY, mX}</code> */
202   static final double[][] Rot_Z_mY_mX = {Z, mY, mX};
203   /** Constant <code>Rot_mX_mY_Z={mX, mY, Z}</code> */
204   static final double[][] Rot_mX_mY_Z = {mX, mY, Z};
205   /** Constant <code>mZ = {ZERO, ZERO, -1.0}</code> */
206   private static final double[] mZ = {ZERO, ZERO, -1.0};
207   /** Constant <code>Rot_Y_mX_mZ={Y, mX, mZ}</code> */
208   static final double[][] Rot_Y_mX_mZ = {Y, mX, mZ};
209   /** Constant <code>Rot_mX_Y_mZ={mX, Y, mZ}</code> */
210   static final double[][] Rot_mX_Y_mZ = {mX, Y, mZ};
211   /** Constant <code>Rot_X_mZ_Y={X, mZ, Y}</code> */
212   static final double[][] Rot_X_mZ_Y = {X, mZ, Y};
213   /** Constant <code>Rot_mY_mZ_X={mY, mZ, X}</code> */
214   static final double[][] Rot_mY_mZ_X = {mY, mZ, X};
215   /** Constant <code>Rot_Y_X_mZ={Y, X, mZ}</code> */
216   static final double[][] Rot_Y_X_mZ = {Y, X, mZ};
217   /** Constant <code>Rot_Y_mZ_X={Y, mZ, X}</code> */
218   static final double[][] Rot_Y_mZ_X = {Y, mZ, X};
219   /** Constant <code>Rot_mX_mY_mZ={mX, mY, mZ}</code> */
220   static final double[][] Rot_mX_mY_mZ = {mX, mY, mZ};
221   /** Constant <code>Rot_X_Y_mZ={X, Y, mZ}</code> */
222   static final double[][] Rot_X_Y_mZ = {X, Y, mZ};
223   /** Constant <code>Rot_mZ_mY_mX={mZ, mY, mX}</code> */
224   static final double[][] Rot_mZ_mY_mX = {mZ, mY, mX};
225   /** Constant <code>Rot_X_mZ_mY={X, mZ, mY}</code> */
226   static final double[][] Rot_X_mZ_mY = {X, mZ, mY};
227   /** Constant <code>Rot_mY_mX_mZ={mY, mX, mZ}</code> */
228   static final double[][] Rot_mY_mX_mZ = {mY, mX, mZ};
229   /** Constant <code>Rot_mY_X_mZ={mY, X, mZ}</code> */
230   static final double[][] Rot_mY_X_mZ = {mY, X, mZ};
231   /** Constant <code>Rot_mX_mZ_mY={mX, mZ, mY}</code> */
232   static final double[][] Rot_mX_mZ_mY = {mX, mZ, mY};
233   /** Constant <code>Rot_mZ_mX_mY={mZ, mX, mY}</code> */
234   static final double[][] Rot_mZ_mX_mY = {mZ, mX, mY};
235   /** Constant <code>Rot_mZ_mY_X={mZ, mY, X}</code> */
236   static final double[][] Rot_mZ_mY_X = {mZ, mY, X};
237   /** Constant <code>Rot_mZ_Y_mX={mZ, Y, mX}</code> */
238   static final double[][] Rot_mZ_Y_mX = {mZ, Y, mX};
239   /** Constant <code>Rot_mZ_mX_Y={mZ, mX, Y}</code> */
240   static final double[][] Rot_mZ_mX_Y = {mZ, mX, Y};
241   /** Constant <code>Rot_mX_mZ_Y={mX, mZ, Y}</code> */
242   static final double[][] Rot_mX_mZ_Y = {mX, mZ, Y};
243   /** Constant <code>Rot_X_mY_mZ={X, mY, mZ}</code> */
244   static final double[][] Rot_X_mY_mZ = {X, mY, mZ};
245   /** Constant <code>Rot_mZ_X_Y={mZ, X, Y}</code> */
246   static final double[][] Rot_mZ_X_Y = {mZ, X, Y};
247   /** Constant <code>Rot_Y_mZ_mX={Y, mZ, mX}</code> */
248   static final double[][] Rot_Y_mZ_mX = {Y, mZ, mX};
249   /** Constant <code>Rot_mY_mZ_mX={mY, mZ, mX}</code> */
250   static final double[][] Rot_mY_mZ_mX = {mY, mZ, mX};
251   /** Constant <code>Rot_mZ_Y_X={mZ, Y, X}</code> */
252   static final double[][] Rot_mZ_Y_X = {mZ, Y, X};
253   /** Constant <code>Rot_mZ_X_mY={mZ, X, mY}</code> */
254   static final double[][] Rot_mZ_X_mY = {mZ, X, mY};
255   /** Constant <code>XmY = {1.0, -1.0, ZERO}</code> */
256   private static final double[] XmY = {1.0, -1.0, ZERO};
257   /** Constant <code>Rot_XmY_X_mZ={XmY, X, mZ}</code> */
258   static final double[][] Rot_XmY_X_mZ = {XmY, X, mZ};
259   /** Constant <code>Rot_XmY_X_Z={XmY, X, Z}</code> */
260   static final double[][] Rot_XmY_X_Z = {XmY, X, Z};
261   /** Constant <code>Rot_XmY_mY_Z={XmY, mY, Z}</code> */
262   static final double[][] Rot_XmY_mY_Z = {XmY, mY, Z};
263   /** Constant <code>Rot_X_XmY_Z={X, XmY, Z}</code> */
264   static final double[][] Rot_X_XmY_Z = {X, XmY, Z};
265   /** Constant <code>Rot_X_XmY_mZ={X, XmY, mZ}</code> */
266   static final double[][] Rot_X_XmY_mZ = {X, XmY, mZ};
267   /** Constant <code>Rot_mY_XmY_mZ={mY, XmY, mZ}</code> */
268   static final double[][] Rot_mY_XmY_mZ = {mY, XmY, mZ};
269   /** Constant <code>Rot_mY_XmY_Z={mY, XmY, Z}</code> */
270   static final double[][] Rot_mY_XmY_Z = {mY, XmY, Z};
271   /** Constant <code>Rot_XmY_mY_mZ={XmY, mY, mZ}</code> */
272   static final double[][] Rot_XmY_mY_mZ = {XmY, mY, mZ};
273   /** Constant <code>mXY = {-1.0, 1.0, ZERO}</code> */
274   private static final double[] mXY = {-1.0, 1.0, ZERO};
275   /** Constant <code>Rot_Y_mXY_Z={Y, mXY, Z}</code> */
276   static final double[][] Rot_Y_mXY_Z = {Y, mXY, Z};
277   /** Constant <code>Rot_mX_mXY_mZ={mX, mXY, mZ}</code> */
278   static final double[][] Rot_mX_mXY_mZ = {mX, mXY, mZ};
279   /** Constant <code>Rot_mXY_Y_Z={mXY, Y, Z}</code> */
280   static final double[][] Rot_mXY_Y_Z = {mXY, Y, Z};
281   /** Constant <code>Rot_mXY_mX_Z={mXY, mX, Z}</code> */
282   static final double[][] Rot_mXY_mX_Z = {mXY, mX, Z};
283   /** Constant <code>Rot_mXY_Y_mZ={mXY, Y, mZ}</code> */
284   static final double[][] Rot_mXY_Y_mZ = {mXY, Y, mZ};
285   /** Constant <code>Rot_mXY_mX_mZ={mXY, mX, mZ}</code> */
286   static final double[][] Rot_mXY_mX_mZ = {mXY, mX, mZ};
287   /** Constant <code>Rot_mX_mXY_Z={mX, mXY, Z}</code> */
288   static final double[][] Rot_mX_mXY_Z = {mX, mXY, Z};
289   /** Constant <code>Rot_Y_mXY_mZ={Y, mXY, mZ}</code> */
290   static final double[][] Rot_Y_mXY_mZ = {Y, mXY, mZ};
291   /** The rotation matrix in fractional coordinates. */
292   public final double[][] rot;
293   /** The translation vector in fractional coordinates. */
294   public final double[] tr;
295   /** A mask equal to 0 for X-coordinates. */
296   private static final int XX = 0;
297   /** A mask equal to 1 for Y-coordinates. */
298   private static final int YY = 1;
299   /** A mask equal to 2 for Z-coordinates. */
300   private static final int ZZ = 2;
301   /** Replicates position for the given symmetry operator (LxMxN). */
302   public final int[] replicatesVector;
303 
304   /**
305    * The SymOp constructor using a rotation matrix and translation vector.
306    *
307    * @param rot The rotation matrix.
308    * @param tr The translation vector.
309    */
310   public SymOp(double[][] rot, double[] tr) {
311     this.rot = rot;
312     this.tr = tr;
313     replicatesVector = new int[]{0, 0, 0};
314   }
315 
316   /**
317    * The SymOp constructor using a rotation matrix and translation vector.
318    *
319    * @param rot The rotation matrix.
320    * @param tr The translation vector.
321    * @param replicatesVector Describes symmetry operators location within replicates crystal.
322    */
323   public SymOp(double[][] rot, double[] tr, int[] replicatesVector) {
324     this.rot = rot;
325     this.tr = tr;
326     this.replicatesVector = replicatesVector;
327   }
328 
329 
330   /**
331    * The SymOp constructor using a 4x4 matrix.
332    *
333    * @param m The rotation matrix and translation vector as a 4x4 matrix.
334    */
335   public SymOp(double[][] m) {
336     this.rot = new double[3][3];
337     rot[0][0] = m[0][0];
338     rot[0][1] = m[0][1];
339     rot[0][2] = m[0][2];
340     rot[1][0] = m[1][0];
341     rot[1][1] = m[1][1];
342     rot[1][2] = m[1][2];
343     rot[2][0] = m[2][0];
344     rot[2][1] = m[2][1];
345     rot[2][2] = m[2][2];
346 
347     this.tr = new double[3];
348     tr[0] = m[0][3] / m[3][3];
349     tr[1] = m[1][3] / m[3][3];
350     tr[2] = m[2][3] / m[3][3];
351 
352     replicatesVector = new int[]{0, 0, 0};
353   }
354 
355   /**
356    * symPhaseShift
357    *
358    * @param hkl an array of double.
359    * @return a double.
360    */
361   public double symPhaseShift(double[] hkl) {
362     // Apply translation
363     return -2.0 * PI * (hkl[0] * tr[0] + hkl[1] * tr[1] + hkl[2] * tr[2]);
364   }
365 
366   /**
367    * symPhaseShift
368    *
369    * @param hkl a {@link HKL} object.
370    * @return a double.
371    */
372   public double symPhaseShift(HKL hkl) {
373     // Apply translation
374     return -2.0 * PI * (hkl.getH() * tr[0] + hkl.getK() * tr[1] + hkl.getL() * tr[2]);
375   }
376 
377   /**
378    * Return the SymOp as a 4x4 matrix.
379    *
380    * @return A 4x4 matrix representation of the SymOp.
381    */
382   public double[][] asMatrix() {
383     return new double[][] {
384         {rot[0][0], rot[0][1], rot[0][2], tr[0]},
385         {rot[1][0], rot[1][1], rot[1][2], tr[1]},
386         {rot[2][0], rot[2][1], rot[2][2], tr[2]},
387         {0.0, 0.0, 0.0, 1.0}};
388   }
389 
390   /**
391    * Return the combined SymOp that is equivalent to first applying <code>this</code> SymOp and then
392    * the argument. Note: Applied as rotation then translation.
393    * <code>X' = S_arg(S_this(X))</code>
394    * <code>X' = S_combined(X)</code>
395    *
396    * @param symOp The SymOp to append to <code>this</code> SymOp.
397    * @return The combined SymOp.
398    */
399   public SymOp append(SymOp symOp) {
400     return new SymOp(mat4Mat4(symOp.asMatrix(), asMatrix()));
401   }
402 
403   /**
404    * Return the combined SymOp that is equivalent to first applying the argument and then
405    * <code>this</code> SymOp. Note: Applied as rotation then translation.
406    * <code>X' = S_this(S_arg(X))</code>
407    * <code>X' = S_combined(X)</code>
408    *
409    * @param symOp The SymOp to prepend to <code>this</code> SymOp.
410    * @return The combined SymOp.
411    */
412   public SymOp prepend(SymOp symOp) {
413     return new SymOp(mat4Mat4(asMatrix(), symOp.asMatrix()));
414   }
415 
416   /**
417    * Return the combined SymOp that is equivalent to first applying symOp1 and then SymOp2. Note:
418    * Applied as rotation then translation.
419    * <p>
420    * <code>X' = S_2(S_1(X))</code>
421    * <code>X' = S_combined(X)</code>
422    *
423    * @param symOp1 The fist SymOp.
424    * @param symOp2 The second SymOp.
425    * @return The combined SymOp.
426    */
427   public static SymOp combineSymOps(SymOp symOp1, SymOp symOp2) {
428     return new SymOp(mat4Mat4(symOp2.asMatrix(), symOp1.asMatrix()));
429   }
430 
431   /**
432    * Apply a Cartesian symmetry operator to an array of Cartesian coordinates. If the arrays x, y or
433    * z are null or not of length n, the method returns immediately. If mateX, mateY or mateZ are null
434    * or not of length n, the method returns immediately.
435    *
436    * @param n Number of atoms.
437    * @param x Input cartesian x-coordinates.
438    * @param y Input cartesian y-coordinates.
439    * @param z Input cartesian z-coordinates.
440    * @param mateX Output cartesian x-coordinates.
441    * @param mateY Output cartesian y-coordinates.
442    * @param mateZ Output cartesian z-coordinates.
443    * @param symOp The cartesian symmetry operator.
444    */
445   public static void applyCartSymOp(int n, double[] x, double[] y, double[] z,
446       double[] mateX, double[] mateY, double[] mateZ, SymOp symOp) {
447     if (x == null || y == null || z == null) {
448       throw new IllegalArgumentException("The input arrays x, y and z must not be null.");
449     }
450     if (x.length < n || y.length < n || z.length < n) {
451       throw new IllegalArgumentException("The input arrays x, y and z must be of length n: " + n);
452     }
453     if (mateX == null || mateY == null || mateZ == null) {
454       throw new IllegalArgumentException("The output arrays mateX, mateY and mateZ must not be null.");
455     }
456     if (mateX.length < n || mateY.length < n || mateZ.length < n) {
457       throw new IllegalArgumentException("The output arrays mateX, mateY and mateZ must be of length n: " + n);
458     }
459 
460     final double[][] rot = symOp.rot;
461     final double[] trans = symOp.tr;
462 
463     final double rot00 = rot[0][0];
464     final double rot10 = rot[1][0];
465     final double rot20 = rot[2][0];
466     final double rot01 = rot[0][1];
467     final double rot11 = rot[1][1];
468     final double rot21 = rot[2][1];
469     final double rot02 = rot[0][2];
470     final double rot12 = rot[1][2];
471     final double rot22 = rot[2][2];
472     final double t0 = trans[0];
473     final double t1 = trans[1];
474     final double t2 = trans[2];
475     for (int i = 0; i < n; i++) {
476       double xc = x[i];
477       double yc = y[i];
478       double zc = z[i];
479       // Apply Symmetry Operator.
480       mateX[i] = rot00 * xc + rot01 * yc + rot02 * zc + t0;
481       mateY[i] = rot10 * xc + rot11 * yc + rot12 * zc + t1;
482       mateZ[i] = rot20 * xc + rot21 * yc + rot22 * zc + t2;
483     }
484   }
485 
486   /**
487    * Apply a cartesian symmetry operator to an array of coordinates.
488    *
489    * @param xyz Input  cartesian coordinates.
490    * @param mate Symmetry mate  cartesian coordinates.
491    * @param symOp The cartesian symmetry operator.
492    */
493   public static void applyCartesianSymOp(double[] xyz, double[] mate, SymOp symOp) {
494     applyCartesianSymOp(xyz, mate, symOp, null);
495   }
496 
497   /**
498    * Apply a cartesian symmetry operator to an array of coordinates.
499    *
500    * @param xyz Input  cartesian coordinates.
501    * @param mate Symmetry mate  cartesian coordinates.
502    * @param symOp The cartesian symmetry operator.
503    * @param mask Only apply the SymOp if the per atom mask is true.
504    */
505   public static void applyCartesianSymOp(double[] xyz, double[] mate, SymOp symOp, boolean[] mask) {
506     var rot = symOp.rot;
507     var trans = symOp.tr;
508 
509     assert (xyz.length % 3 == 0);
510     assert (xyz.length == mate.length);
511     int len = xyz.length / 3;
512 
513     // Load the SymOp into local variables.
514     var rot00 = rot[0][0];
515     var rot10 = rot[1][0];
516     var rot20 = rot[2][0];
517     var rot01 = rot[0][1];
518     var rot11 = rot[1][1];
519     var rot21 = rot[2][1];
520     var rot02 = rot[0][2];
521     var rot12 = rot[1][2];
522     var rot22 = rot[2][2];
523     var tx = trans[0];
524     var ty = trans[1];
525     var tz = trans[2];
526 
527     if (mask == null) {
528       for (int i = 0; i < len; i++) {
529         int index = i * 3;
530         var xc = xyz[index + XX];
531         var yc = xyz[index + YY];
532         var zc = xyz[index + ZZ];
533         // Apply Symmetry Operator.
534         mate[index + XX] = rot00 * xc + rot01 * yc + rot02 * zc + tx;
535         mate[index + YY] = rot10 * xc + rot11 * yc + rot12 * zc + ty;
536         mate[index + ZZ] = rot20 * xc + rot21 * yc + rot22 * zc + tz;
537       }
538     } else {
539       for (int i = 0; i < len; i++) {
540         int index = i * 3;
541         var xc = xyz[index + XX];
542         var yc = xyz[index + YY];
543         var zc = xyz[index + ZZ];
544         if (mask[i]) {
545           // Apply Symmetry Operator.
546           mate[index + XX] = rot00 * xc + rot01 * yc + rot02 * zc + tx;
547           mate[index + YY] = rot10 * xc + rot11 * yc + rot12 * zc + ty;
548           mate[index + ZZ] = rot20 * xc + rot21 * yc + rot22 * zc + tz;
549         } else {
550           mate[index + XX] = xc;
551           mate[index + YY] = yc;
552           mate[index + ZZ] = zc;
553         }
554       }
555     }
556   }
557 
558   /**
559    * Apply a fractional symmetry operator to one set of coordinates.
560    *
561    * @param xyz Input fractional coordinates.
562    * @param mate Symmetry mate fractional coordinates.
563    * @param symOp The fractional symmetry operator.
564    */
565   public static void applyFracSymOp(double[] xyz, double[] mate, SymOp symOp) {
566     var rot = symOp.rot;
567     var trans = symOp.tr;
568     var xf = xyz[0];
569     var yf = xyz[1];
570     var zf = xyz[2];
571     // Apply Symmetry Operator.
572     mate[0] = rot[0][0] * xf + rot[0][1] * yf + rot[0][2] * zf + trans[0];
573     mate[1] = rot[1][0] * xf + rot[1][1] * yf + rot[1][2] * zf + trans[1];
574     mate[2] = rot[2][0] * xf + rot[2][1] * yf + rot[2][2] * zf + trans[2];
575   }
576 
577   /**
578    * Apply a symmetry operator to one set of coordinates.
579    *
580    * @param h Input coordinates.
581    * @param k Input coordinates.
582    * @param l Input coordinates.
583    * @param mate Symmetry mate coordinates.
584    * @param symOp The symmetry operator.
585    * @param nx number of unit cell translations
586    * @param ny number of unit cell translations
587    * @param nz number of unit cell translations
588    */
589   public static void applySymOp(int h, int k, int l, int[] mate, SymOp symOp, int nx, int ny,
590       int nz) {
591     var rot = symOp.rot;
592     var trans = symOp.tr;
593     // Apply Symmetry Operator.
594     mate[0] =
595         (int) rot[0][0] * h + (int) rot[0][1] * k + (int) rot[0][2] * l + (int) rint(nx * trans[0]);
596     mate[1] =
597         (int) rot[1][0] * h + (int) rot[1][1] * k + (int) rot[1][2] * l + (int) rint(ny * trans[1]);
598     mate[2] =
599         (int) rot[2][0] * h + (int) rot[2][1] * k + (int) rot[2][2] * l + (int) rint(nz * trans[2]);
600     mate[0] = mod(mate[0], nx);
601     mate[1] = mod(mate[1], ny);
602     mate[2] = mod(mate[2], nz);
603   }
604 
605   /**
606    * Apply a symmetry operator to one HKL.
607    *
608    * @param hkl Input HKL.
609    * @param mate Symmetry mate HKL.
610    * @param symOp The symmetry operator.
611    */
612   public static void applySymRot(HKL hkl, HKL mate, SymOp symOp) {
613     var rot = symOp.rot;
614     double h = hkl.getH();
615     double k = hkl.getK();
616     double l = hkl.getL();
617     double hs = rot[0][0] * h + rot[0][1] * k + rot[0][2] * l;
618     double ks = rot[1][0] * h + rot[1][1] * k + rot[1][2] * l;
619     double ls = rot[2][0] * h + rot[2][1] * k + rot[2][2] * l;
620     // Convert back to HKL
621     mate.setH((int) rint(hs));
622     mate.setK((int) rint(ks));
623     mate.setL((int) rint(ls));
624   }
625 
626   /**
627    * Apply a Cartesian symmetry rotation to an array of Cartesian coordinates. The length of xyz must
628    * be divisible by 3 and mate must have the same length.
629    *
630    * @param xyz Input cartesian x, y, z-coordinates.
631    * @param mate Output cartesian x, y, z-coordinates.
632    * @param symOp The fractional symmetry operator.
633    */
634   public static void applyCartesianSymRot(double[] xyz, double[] mate, SymOp symOp) {
635     applyCartesianSymRot(xyz, mate, symOp, null);
636   }
637 
638   /**
639    * Apply a Cartesian symmetry rotation to an array of Cartesian coordinates. The length of xyz must
640    * be divisible by 3 and mate must have the same length.
641    *
642    * @param xyz Input cartesian x, y, z-coordinates.
643    * @param mate Output cartesian x, y, z-coordinates.
644    * @param symOp The fractional symmetry operator.
645    * @param mask Only apply the SymOp if the per atom mask is true.
646    */
647   public static void applyCartesianSymRot(double[] xyz, double[] mate, SymOp symOp, boolean[] mask) {
648     int l = xyz.length;
649     int n = l / 3;
650     assert (l % 3 == 0);
651     assert (mate.length == l);
652 
653     // Load the rotation matrix
654     var rot = symOp.rot;
655     var rot00 = rot[0][0];
656     var rot10 = rot[1][0];
657     var rot20 = rot[2][0];
658     var rot01 = rot[0][1];
659     var rot11 = rot[1][1];
660     var rot21 = rot[2][1];
661     var rot02 = rot[0][2];
662     var rot12 = rot[1][2];
663     var rot22 = rot[2][2];
664 
665     if (mask == null) {
666       for (int i = 0; i < n; i++) {
667         int index = i * 3;
668         var xi = xyz[index + XX];
669         var yi = xyz[index + YY];
670         var zi = xyz[index + ZZ];
671         // Apply Symmetry Operator.
672         mate[index + XX] = rot00 * xi + rot01 * yi + rot02 * zi;
673         mate[index + YY] = rot10 * xi + rot11 * yi + rot12 * zi;
674         mate[index + ZZ] = rot20 * xi + rot21 * yi + rot22 * zi;
675       }
676     } else {
677       for (int i = 0; i < n; i++) {
678         int index = i * 3;
679         var xi = xyz[index + XX];
680         var yi = xyz[index + YY];
681         var zi = xyz[index + ZZ];
682         if (mask[i]) {
683           // Apply Symmetry Operator.
684           mate[index + XX] = rot00 * xi + rot01 * yi + rot02 * zi;
685           mate[index + YY] = rot10 * xi + rot11 * yi + rot12 * zi;
686           mate[index + ZZ] = rot20 * xi + rot21 * yi + rot22 * zi;
687         } else {
688           mate[index + XX] = xi;
689           mate[index + YY] = yi;
690           mate[index + ZZ] = zi;
691         }
692       }
693     }
694   }
695 
696   /**
697    * Apply a transpose rotation symmetry operator to one HKL.
698    *
699    * @param hkl Input HKL.
700    * @param mate Symmetry mate HKL.
701    * @param symOp The symmetry operator.
702    */
703   public static void applyTransSymRot(HKL hkl, HKL mate, SymOp symOp) {
704     double[][] rot = symOp.rot;
705     double h = hkl.getH();
706     double k = hkl.getK();
707     double l = hkl.getL();
708     // Apply transpose Symmetry Operator.
709     double hs = rot[0][0] * h + rot[1][0] * k + rot[2][0] * l;
710     double ks = rot[0][1] * h + rot[1][1] * k + rot[2][1] * l;
711     double ls = rot[0][2] * h + rot[1][2] * k + rot[2][2] * l;
712     // Convert back to HKL
713     mate.setH((int) rint(hs));
714     mate.setK((int) rint(ks));
715     mate.setL((int) rint(ls));
716   }
717 
718   /**
719    * Print a Sym Op matrix as a continued line string.
720    *
721    * @param symOp Symmetry operation to print.
722    * @return Continued line string.
723    */
724   public static String asMatrixString(SymOp symOp) {
725     double[][] values = symOp.asMatrix();
726     return format("""
727                      %14.8f %14.8f %14.8f \\
728                                %14.8f %14.8f %14.8f \\
729                                %14.8f %14.8f %14.8f \\
730                                %14.8f %14.8f %14.8f \
731                     """, values[0][0],
732             values[0][1], values[0][2], values[1][0], values[1][1], values[1][2], values[2][0],
733             values[2][1], values[2][2], values[0][3], values[1][3], values[2][3]);
734   }
735 
736   /**
737    * Invert a symmetry operator.
738    *
739    * @param symOp Original symmetry operator of which the inverse is desired.
740    * @return SymOp The inverse symmetry operator of the one supplied.
741    */
742   public static SymOp invertSymOp(SymOp symOp) {
743     var tr = symOp.tr;
744     var rot = symOp.rot;
745     var inv = mat3Inverse(rot);
746     return new SymOp(inv,
747         new double[] {-dot(inv[0], tr), -dot(inv[1], tr), -dot(inv[2], tr)});
748   }
749 
750   /**
751    * Create a SymOp from an input String.
752    *
753    * @param s Input <code>String</code> containing 12 white-space delimited double values.
754    * @return The SymOp.
755    */
756   public static SymOp parse(String s) {
757     String[] tokens = s.split(" +");
758     if (tokens.length < 12) {
759       return null;
760     }
761     return new SymOp(new double[][] {
762         {parseDouble(tokens[0]), parseDouble(tokens[1]), parseDouble(tokens[2])},
763         {parseDouble(tokens[3]), parseDouble(tokens[4]), parseDouble(tokens[5])},
764         {parseDouble(tokens[6]), parseDouble(tokens[7]), parseDouble(tokens[8])}},
765         new double[] {parseDouble(tokens[9]), parseDouble(tokens[10]), parseDouble(tokens[11])});
766   }
767 
768   /**
769    * Generate a random Cartesian Symmetry Operator.
770    *
771    * @param scalar The range of translations will be from -scalar/2 .. scalar/2.
772    * @return A Cartesian SymOp with a random rotation and translation.
773    */
774   public static SymOp randomSymOpFactory(double scalar) {
775     double[] tr = {scalar * (random() - 0.5), scalar * (random() - 0.5), scalar * (random() - 0.5)};
776     return randomSymOpFactory(tr);
777   }
778 
779   /**
780    * Generate a random Cartesian Symmetry Operator.
781    *
782    * <p>The random rotation matrix is derived from: Arvo, James (1992), "Fast random rotation
783    * matrices", in David Kirk, Graphics Gems III, San Diego: Academic Press Professional, pp.
784    * 117–120, ISBN 978-0-12-409671-4
785    *
786    * @param tr The translations to apply.
787    * @return A Cartesian SymOp with a random rotation and translation.
788    */
789   public static SymOp randomSymOpFactory(double[] tr) {
790     double[][] rot = new double[3][3];
791     double PI2 = 2.0 * PI;
792     double[] x = new double[3];
793     x[0] = random();
794     x[1] = random();
795     x[2] = random();
796     /* Rotation about the pole (Z).      */
797     double theta = x[0] * PI2;
798     /* For direction of pole deflection. */
799     double phi = x[1] * PI2;
800     /* For magnitude of pole deflection. */
801     double z = x[2] * 2.0;
802     /*
803      Compute a vector V used for distributing points over the sphere via
804      the reflection I - V Transpose(V). This formulation of V will
805      guarantee that if x[1] and x[2] are uniformly distributed, the
806      reflected points will be uniform on the sphere. Note that V has
807      length sqrt(2) to eliminate the 2 in the Householder matrix.
808     */
809     double r = sqrt(z);
810     double Vx = sin(phi) * r;
811     double Vy = cos(phi) * r;
812     double Vz = sqrt(2.0 - z);
813     /*
814      Compute the row vector S = Transpose(V) * R, where R is a simple
815      rotation by theta about the z-axis. No need to compute Sz since it's
816      just Vz.
817     */
818     double st = sin(theta);
819     double ct = cos(theta);
820     double Sx = Vx * ct - Vy * st;
821     double Sy = Vx * st + Vy * ct;
822 
823     // Construct the rotation matrix ( V Transpose(V) - I ) R, which is equivalent to V S - R.
824     rot[0][0] = Vx * Sx - ct;
825     rot[0][1] = Vx * Sy - st;
826     rot[0][2] = Vx * Vz;
827     rot[1][0] = Vy * Sx + st;
828     rot[1][1] = Vy * Sy - ct;
829     rot[1][2] = Vy * Vz;
830     rot[2][0] = Vz * Sx;
831     rot[2][1] = Vz * Sy;
832     // This equals Vz * Vz - 1.0
833     rot[2][2] = 1.0 - z;
834 
835     return new SymOp(rot, tr);
836   }
837 
838   /**
839    * trtoString
840    *
841    * @param tr a double.
842    * @return a {@link java.lang.String} object.
843    */
844   private static String trtoString(double tr) {
845     if (tr == f12) {
846       return "+1/2";
847     }
848     if (tr == f13) {
849       return "+1/3";
850     }
851     if (tr == f23) {
852       return "+2/3";
853     }
854     if (tr == f14) {
855       return "+1/4";
856     }
857     if (tr == f34) {
858       return "+3/4";
859     }
860     if (tr == f16) {
861       return "+1/6";
862     }
863     if (tr == f56) {
864       return "+5/6";
865     }
866 
867     return "";
868   }
869 
870   /** {@inheritDoc} */
871   @Override
872   public String toString() {
873     StringBuilder sb = new StringBuilder(" Rotation operator:\n");
874     sb.append(format(" [[%4.1f,%4.1f,%4.1f]\n  [%4.1f,%4.1f,%4.1f]\n  [%4.1f,%4.1f,%4.1f]]\n",
875         rot[0][0], rot[0][1], rot[0][2],
876         rot[1][0], rot[1][1], rot[1][2],
877         rot[2][0], rot[2][1], rot[2][2]));
878     sb.append(" Translation:\n");
879     sb.append(format(" [%4.2f,%4.2f,%4.2f]", tr[0], tr[1], tr[2]));
880     return sb.toString();
881   }
882 
883   /**
884    * Print the symmetry operator with double precision.
885    *
886    * @return String of rotation/translation with double precision.
887    */
888   public String toStringPrecise() {
889     StringBuilder sb = new StringBuilder(" Rotation operator:\n");
890     sb.append(format(
891         " [[%18.16e,%18.16e,%18.16e]\n  [%18.16e,%18.16e,%18.16e]\n  [%18.16e,%18.16e,%18.16e]]\n",
892         rot[0][0], rot[0][1], rot[0][2],
893         rot[1][0], rot[1][1], rot[1][2],
894         rot[2][0], rot[2][1], rot[2][2]));
895     sb.append(" Translation:\n");
896     sb.append(format(" [%18.16e,%18.16e,%18.16e]", tr[0], tr[1], tr[2]));
897     return sb.toString();
898   }
899 
900   /**
901    * toXYZString
902    *
903    * @return a {@link java.lang.String} object.
904    */
905   public String toXYZString() {
906     StringBuilder sb = new StringBuilder();
907     for (int i = 0; i < 3; i++) {
908       boolean s = false;
909       if (rot[i][0] < 0.0) {
910         sb.append("-X");
911         s = true;
912       } else if (rot[i][0] > 0.0) {
913         sb.append("X");
914         s = true;
915       }
916       if (rot[i][1] < 0.0) {
917         sb.append("-Y");
918         s = true;
919       } else if (rot[i][1] > 0.0) {
920         sb.append(s ? "+Y" : "Y");
921         s = true;
922       }
923       if (rot[i][2] < 0.0) {
924         sb.append("-Z");
925       } else if (rot[i][2] > 0.0) {
926         sb.append(s ? "+Z" : "Z");
927       }
928 
929       if (tr[i] > 0.0) {
930         sb.append(trtoString(tr[i]));
931       }
932 
933       if (i < 2) {
934         sb.append(", ");
935       }
936     }
937     return sb.toString();
938   }
939 }