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.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 
302   /**
303    * The SymOp constructor using a rotation matrix and translation vector.
304    *
305    * @param rot The rotation matrix.
306    * @param tr The translation vector.
307    */
308   public SymOp(double[][] rot, double[] tr) {
309     this.rot = rot;
310     this.tr = tr;
311   }
312 
313 
314   /**
315    * The SymOp constructor using a 4x4 matrix.
316    *
317    * @param m The rotation matrix and translation vector as a 4x4 matrix.
318    */
319   public SymOp(double[][] m) {
320     this.rot = new double[3][3];
321     rot[0][0] = m[0][0];
322     rot[0][1] = m[0][1];
323     rot[0][2] = m[0][2];
324     rot[1][0] = m[1][0];
325     rot[1][1] = m[1][1];
326     rot[1][2] = m[1][2];
327     rot[2][0] = m[2][0];
328     rot[2][1] = m[2][1];
329     rot[2][2] = m[2][2];
330 
331     this.tr = new double[3];
332     tr[0] = m[0][3] / m[3][3];
333     tr[1] = m[1][3] / m[3][3];
334     tr[2] = m[2][3] / m[3][3];
335   }
336 
337   /**
338    * symPhaseShift
339    *
340    * @param hkl an array of double.
341    * @return a double.
342    */
343   public double symPhaseShift(double[] hkl) {
344     // Apply translation
345     return -2.0 * PI * (hkl[0] * tr[0] + hkl[1] * tr[1] + hkl[2] * tr[2]);
346   }
347 
348   /**
349    * symPhaseShift
350    *
351    * @param hkl a {@link HKL} object.
352    * @return a double.
353    */
354   public double symPhaseShift(HKL hkl) {
355     // Apply translation
356     return -2.0 * PI * (hkl.getH() * tr[0] + hkl.getK() * tr[1] + hkl.getL() * tr[2]);
357   }
358 
359   /**
360    * Return the SymOp as a 4x4 matrix.
361    *
362    * @return A 4x4 matrix representation of the SymOp.
363    */
364   public double[][] asMatrix() {
365     return new double[][] {
366         {rot[0][0], rot[0][1], rot[0][2], tr[0]},
367         {rot[1][0], rot[1][1], rot[1][2], tr[1]},
368         {rot[2][0], rot[2][1], rot[2][2], tr[2]},
369         {0.0, 0.0, 0.0, 1.0}};
370   }
371 
372   /**
373    * Return the combined SymOp that is equivalent to first applying <code>this</code> SymOp and then
374    * the argument. Note: Applied as rotation then translation.
375    * <code>X' = S_arg(S_this(X))</code>
376    * <code>X' = S_combined(X)</code>
377    *
378    * @param symOp The SymOp to append to <code>this</code> SymOp.
379    * @return The combined SymOp.
380    */
381   public SymOp append(SymOp symOp) {
382     return new SymOp(mat4Mat4(symOp.asMatrix(), asMatrix()));
383   }
384 
385   /**
386    * Return the combined SymOp that is equivalent to first applying the argument and then
387    * <code>this</code> SymOp. Note: Applied as rotation then translation.
388    * <code>X' = S_this(S_arg(X))</code>
389    * <code>X' = S_combined(X)</code>
390    *
391    * @param symOp The SymOp to prepend to <code>this</code> Symop.
392    * @return The combined SymOp.
393    */
394   public SymOp prepend(SymOp symOp) {
395     return new SymOp(mat4Mat4(asMatrix(), symOp.asMatrix()));
396   }
397 
398   /**
399    * Return the combined SymOp that is equivalent to first applying symOp1 and then SymOp2. Note:
400    * Applied as rotation then translation.
401    *
402    * <code>X' = S_2(S_1(X))</code>
403    * <code>X' = S_combined(X)</code>
404    *
405    * @param symOp1 The fist SymOp.
406    * @param symOp2 The second SymOp.
407    * @return The combined SymOp.
408    */
409   public static SymOp combineSymOps(SymOp symOp1, SymOp symOp2) {
410     return new SymOp(mat4Mat4(symOp2.asMatrix(), symOp1.asMatrix()));
411   }
412 
413   /**
414    * Apply a Cartesian symmetry operator to an array of Cartesian coordinates. If the arrays x, y or
415    * z are null or not of length n, the method returns immediately. If mateX, mateY or mateZ are null
416    * or not of length n, new arrays are allocated.
417    *
418    * @param n Number of atoms.
419    * @param x Input cartesian x-coordinates.
420    * @param y Input cartesian y-coordinates.
421    * @param z Input cartesian z-coordinates.
422    * @param mateX Output cartesian x-coordinates.
423    * @param mateY Output cartesian y-coordinates.
424    * @param mateZ Output cartesian z-coordinates.
425    * @param symOp The cartesian symmetry operator.
426    */
427   public static void applyCartSymOp(
428       int n,
429       double[] x,
430       double[] y,
431       double[] z,
432       double[] mateX,
433       double[] mateY,
434       double[] mateZ,
435       SymOp symOp) {
436     if (x == null || y == null || z == null) {
437       return;
438     }
439     if (x.length < n || y.length < n || z.length < n) {
440       return;
441     }
442     if (mateX == null || mateX.length < n) {
443       mateX = new double[n];
444     }
445     if (mateY == null || mateY.length < n) {
446       mateY = new double[n];
447     }
448     if (mateZ == null || mateZ.length < n) {
449       mateZ = new double[n];
450     }
451 
452     final double[][] rot = symOp.rot;
453     final double[] trans = symOp.tr;
454 
455     final double rot00 = rot[0][0];
456     final double rot10 = rot[1][0];
457     final double rot20 = rot[2][0];
458     final double rot01 = rot[0][1];
459     final double rot11 = rot[1][1];
460     final double rot21 = rot[2][1];
461     final double rot02 = rot[0][2];
462     final double rot12 = rot[1][2];
463     final double rot22 = rot[2][2];
464     final double t0 = trans[0];
465     final double t1 = trans[1];
466     final double t2 = trans[2];
467     for (int i = 0; i < n; i++) {
468       double xc = x[i];
469       double yc = y[i];
470       double zc = z[i];
471       // Apply Symmetry Operator.
472       mateX[i] = rot00 * xc + rot01 * yc + rot02 * zc + t0;
473       mateY[i] = rot10 * xc + rot11 * yc + rot12 * zc + t1;
474       mateZ[i] = rot20 * xc + rot21 * yc + rot22 * zc + t2;
475     }
476   }
477 
478   /**
479    * Apply a cartesian symmetry operator to an array of coordinates.
480    *
481    * @param xyz Input  cartesian coordinates.
482    * @param mate Symmetry mate  cartesian coordinates.
483    * @param symOp The cartesian symmetry operator.
484    */
485   public static void applyCartesianSymOp(double[] xyz, double[] mate, SymOp symOp) {
486     applyCartesianSymOp(xyz, mate, symOp, null);
487   }
488 
489   /**
490    * Apply a cartesian symmetry operator to an array of coordinates.
491    *
492    * @param xyz Input  cartesian coordinates.
493    * @param mate Symmetry mate  cartesian coordinates.
494    * @param symOp The cartesian symmetry operator.
495    * @param mask Only apply the SymOp if the per atom mask is true.
496    */
497   public static void applyCartesianSymOp(double[] xyz, double[] mate, SymOp symOp, boolean[] mask) {
498     var rot = symOp.rot;
499     var trans = symOp.tr;
500 
501     assert (xyz.length % 3 == 0);
502     assert (xyz.length == mate.length);
503     int len = xyz.length / 3;
504 
505     // Load the SymOp into local variables.
506     var rot00 = rot[0][0];
507     var rot10 = rot[1][0];
508     var rot20 = rot[2][0];
509     var rot01 = rot[0][1];
510     var rot11 = rot[1][1];
511     var rot21 = rot[2][1];
512     var rot02 = rot[0][2];
513     var rot12 = rot[1][2];
514     var rot22 = rot[2][2];
515     var tx = trans[0];
516     var ty = trans[1];
517     var tz = trans[2];
518 
519     if (mask == null) {
520       for (int i = 0; i < len; i++) {
521         int index = i * 3;
522         var xc = xyz[index + XX];
523         var yc = xyz[index + YY];
524         var zc = xyz[index + ZZ];
525         // Apply Symmetry Operator.
526         mate[index + XX] = rot00 * xc + rot01 * yc + rot02 * zc + tx;
527         mate[index + YY] = rot10 * xc + rot11 * yc + rot12 * zc + ty;
528         mate[index + ZZ] = rot20 * xc + rot21 * yc + rot22 * zc + tz;
529       }
530     } else {
531       for (int i = 0; i < len; i++) {
532         int index = i * 3;
533         var xc = xyz[index + XX];
534         var yc = xyz[index + YY];
535         var zc = xyz[index + ZZ];
536         if (mask[i]) {
537           // Apply Symmetry Operator.
538           mate[index + XX] = rot00 * xc + rot01 * yc + rot02 * zc + tx;
539           mate[index + YY] = rot10 * xc + rot11 * yc + rot12 * zc + ty;
540           mate[index + ZZ] = rot20 * xc + rot21 * yc + rot22 * zc + tz;
541         } else {
542           mate[index + XX] = xc;
543           mate[index + YY] = yc;
544           mate[index + ZZ] = zc;
545         }
546       }
547     }
548   }
549 
550   /**
551    * Apply a fractional symmetry operator to one set of coordinates.
552    *
553    * @param xyz Input fractional coordinates.
554    * @param mate Symmetry mate fractional coordinates.
555    * @param symOp The fractional symmetry operator.
556    */
557   public static void applyFracSymOp(double[] xyz, double[] mate, SymOp symOp) {
558     var rot = symOp.rot;
559     var trans = symOp.tr;
560     var xf = xyz[0];
561     var yf = xyz[1];
562     var zf = xyz[2];
563     // Apply Symmetry Operator.
564     mate[0] = rot[0][0] * xf + rot[0][1] * yf + rot[0][2] * zf + trans[0];
565     mate[1] = rot[1][0] * xf + rot[1][1] * yf + rot[1][2] * zf + trans[1];
566     mate[2] = rot[2][0] * xf + rot[2][1] * yf + rot[2][2] * zf + trans[2];
567   }
568 
569   /**
570    * Apply a symmetry operator to one set of coordinates.
571    *
572    * @param h Input coordinates.
573    * @param k Input coordinates.
574    * @param l Input coordinates.
575    * @param mate Symmetry mate coordinates.
576    * @param symOp The symmetry operator.
577    * @param nx number of unit cell translations
578    * @param ny number of unit cell translations
579    * @param nz number of unit cell translations
580    */
581   public static void applySymOp(int h, int k, int l, int[] mate, SymOp symOp, int nx, int ny,
582       int nz) {
583     var rot = symOp.rot;
584     var trans = symOp.tr;
585     // Apply Symmetry Operator.
586     mate[0] =
587         (int) rot[0][0] * h + (int) rot[0][1] * k + (int) rot[0][2] * l + (int) rint(nx * trans[0]);
588     mate[1] =
589         (int) rot[1][0] * h + (int) rot[1][1] * k + (int) rot[1][2] * l + (int) rint(ny * trans[1]);
590     mate[2] =
591         (int) rot[2][0] * h + (int) rot[2][1] * k + (int) rot[2][2] * l + (int) rint(nz * trans[2]);
592     mate[0] = mod(mate[0], nx);
593     mate[1] = mod(mate[1], ny);
594     mate[2] = mod(mate[2], nz);
595   }
596 
597   /**
598    * Apply a symmetry operator to one HKL.
599    *
600    * @param hkl Input HKL.
601    * @param mate Symmetry mate HKL.
602    * @param symOp The symmetry operator.
603    */
604   public static void applySymRot(HKL hkl, HKL mate, SymOp symOp) {
605     var rot = symOp.rot;
606     double h = hkl.getH();
607     double k = hkl.getK();
608     double l = hkl.getL();
609     double hs = rot[0][0] * h + rot[0][1] * k + rot[0][2] * l;
610     double ks = rot[1][0] * h + rot[1][1] * k + rot[1][2] * l;
611     double ls = rot[2][0] * h + rot[2][1] * k + rot[2][2] * l;
612     // Convert back to HKL
613     mate.setH((int) rint(hs));
614     mate.setK((int) rint(ks));
615     mate.setL((int) rint(ls));
616   }
617 
618   /**
619    * Apply a Cartesian symmetry rotation to an array of Cartesian coordinates. The length of xyz must
620    * be divisible by 3 and mate must have the same length.
621    *
622    * @param xyz Input cartesian x, y, z-coordinates.
623    * @param mate Output cartesian x, y, z-coordinates.
624    * @param symOp The fractional symmetry operator.
625    */
626   public static void applyCartesianSymRot(double[] xyz, double[] mate, SymOp symOp) {
627     applyCartesianSymRot(xyz, mate, symOp, null);
628   }
629 
630   /**
631    * Apply a Cartesian symmetry rotation to an array of Cartesian coordinates. The length of xyz must
632    * be divisible by 3 and mate must have the same length.
633    *
634    * @param xyz Input cartesian x, y, z-coordinates.
635    * @param mate Output cartesian x, y, z-coordinates.
636    * @param symOp The fractional symmetry operator.
637    * @param mask Only apply the SymOp if the per atom mask is true.
638    */
639   public static void applyCartesianSymRot(double[] xyz, double[] mate, SymOp symOp, boolean[] mask) {
640     int l = xyz.length;
641     int n = l / 3;
642     assert (l % 3 == 0);
643     assert (mate.length == l);
644 
645     // Load the rotation matrix
646     var rot = symOp.rot;
647     var rot00 = rot[0][0];
648     var rot10 = rot[1][0];
649     var rot20 = rot[2][0];
650     var rot01 = rot[0][1];
651     var rot11 = rot[1][1];
652     var rot21 = rot[2][1];
653     var rot02 = rot[0][2];
654     var rot12 = rot[1][2];
655     var rot22 = rot[2][2];
656 
657     if (mask == null) {
658       for (int i = 0; i < n; i++) {
659         int index = i * 3;
660         var xi = xyz[index + XX];
661         var yi = xyz[index + YY];
662         var zi = xyz[index + ZZ];
663         // Apply Symmetry Operator.
664         mate[index + XX] = rot00 * xi + rot01 * yi + rot02 * zi;
665         mate[index + YY] = rot10 * xi + rot11 * yi + rot12 * zi;
666         mate[index + ZZ] = rot20 * xi + rot21 * yi + rot22 * zi;
667       }
668     } else {
669       for (int i = 0; i < n; i++) {
670         int index = i * 3;
671         var xi = xyz[index + XX];
672         var yi = xyz[index + YY];
673         var zi = xyz[index + ZZ];
674         if (mask[i]) {
675           // Apply Symmetry Operator.
676           mate[index + XX] = rot00 * xi + rot01 * yi + rot02 * zi;
677           mate[index + YY] = rot10 * xi + rot11 * yi + rot12 * zi;
678           mate[index + ZZ] = rot20 * xi + rot21 * yi + rot22 * zi;
679         } else {
680           mate[index + XX] = xi;
681           mate[index + YY] = yi;
682           mate[index + ZZ] = zi;
683         }
684       }
685     }
686   }
687 
688   /**
689    * Apply a transpose rotation symmetry operator to one HKL.
690    *
691    * @param hkl Input HKL.
692    * @param mate Symmetry mate HKL.
693    * @param symOp The symmetry operator.
694    */
695   public static void applyTransSymRot(HKL hkl, HKL mate, SymOp symOp) {
696     double[][] rot = symOp.rot;
697     double h = hkl.getH();
698     double k = hkl.getK();
699     double l = hkl.getL();
700     // Apply transpose Symmetry Operator.
701     double hs = rot[0][0] * h + rot[1][0] * k + rot[2][0] * l;
702     double ks = rot[0][1] * h + rot[1][1] * k + rot[2][1] * l;
703     double ls = rot[0][2] * h + rot[1][2] * k + rot[2][2] * l;
704     // Convert back to HKL
705     mate.setH((int) rint(hs));
706     mate.setK((int) rint(ks));
707     mate.setL((int) rint(ls));
708   }
709 
710   /**
711    * Invert a symmetry operator.
712    *
713    * @param symOp Original symmetry operator of which the inverse is desired.
714    * @return SymOp The inverse symmetry operator of the one supplied.
715    */
716   public static SymOp invertSymOp(SymOp symOp) {
717     var tr = symOp.tr;
718     var rot = symOp.rot;
719     var inv = mat3Inverse(rot);
720     return new SymOp(inv,
721         new double[] {-dot(inv[0], tr), -dot(inv[1], tr), -dot(inv[2], tr)});
722   }
723 
724   /**
725    * Create a SymOp from an input String.
726    *
727    * @param s Input <code>String</code> containing 12 white-space delimited double values.
728    * @return The SymOp.
729    */
730   public static SymOp parse(String s) {
731     String[] tokens = s.split(" +");
732     if (tokens.length < 12) {
733       return null;
734     }
735     return new SymOp(new double[][] {
736         {parseDouble(tokens[0]), parseDouble(tokens[1]), parseDouble(tokens[2])},
737         {parseDouble(tokens[3]), parseDouble(tokens[4]), parseDouble(tokens[5])},
738         {parseDouble(tokens[6]), parseDouble(tokens[7]), parseDouble(tokens[8])}},
739         new double[] {parseDouble(tokens[9]), parseDouble(tokens[10]), parseDouble(tokens[11])});
740   }
741 
742   /**
743    * Generate a random Cartesian Symmetry Operator.
744    *
745    * @param scalar The range of translations will be from -scalar/2 .. scalar/2.
746    * @return A Cartesian SymOp with a random rotation and translation.
747    */
748   public static SymOp randomSymOpFactory(double scalar) {
749     double[] tr = {scalar * (random() - 0.5), scalar * (random() - 0.5), scalar * (random() - 0.5)};
750     return randomSymOpFactory(tr);
751   }
752 
753   /**
754    * Generate a random Cartesian Symmetry Operator.
755    *
756    * <p>The random rotation matrix is derived from: Arvo, James (1992), "Fast random rotation
757    * matrices", in David Kirk, Graphics Gems III, San Diego: Academic Press Professional, pp.
758    * 117–120, ISBN 978-0-12-409671-4
759    *
760    * @param tr The translations to apply.
761    * @return A Cartesian SymOp with a random rotation and translation.
762    */
763   public static SymOp randomSymOpFactory(double[] tr) {
764     double[][] rot = new double[3][3];
765     double PI2 = 2.0 * PI;
766     double[] x = new double[3];
767     x[0] = random();
768     x[1] = random();
769     x[2] = random();
770     /* Rotation about the pole (Z).      */
771     double theta = x[0] * PI2;
772     /* For direction of pole deflection. */
773     double phi = x[1] * PI2;
774     /* For magnitude of pole deflection. */
775     double z = x[2] * 2.0;
776     /*
777      Compute a vector V used for distributing points over the sphere via
778      the reflection I - V Transpose(V). This formulation of V will
779      guarantee that if x[1] and x[2] are uniformly distributed, the
780      reflected points will be uniform on the sphere. Note that V has
781      length sqrt(2) to eliminate the 2 in the Householder matrix.
782     */
783     double r = sqrt(z);
784     double Vx = sin(phi) * r;
785     double Vy = cos(phi) * r;
786     double Vz = sqrt(2.0 - z);
787     /*
788      Compute the row vector S = Transpose(V) * R, where R is a simple
789      rotation by theta about the z-axis. No need to compute Sz since it's
790      just Vz.
791     */
792     double st = sin(theta);
793     double ct = cos(theta);
794     double Sx = Vx * ct - Vy * st;
795     double Sy = Vx * st + Vy * ct;
796 
797     // Construct the rotation matrix ( V Transpose(V) - I ) R, which is equivalent to V S - R.
798     rot[0][0] = Vx * Sx - ct;
799     rot[0][1] = Vx * Sy - st;
800     rot[0][2] = Vx * Vz;
801     rot[1][0] = Vy * Sx + st;
802     rot[1][1] = Vy * Sy - ct;
803     rot[1][2] = Vy * Vz;
804     rot[2][0] = Vz * Sx;
805     rot[2][1] = Vz * Sy;
806     // This equals Vz * Vz - 1.0
807     rot[2][2] = 1.0 - z;
808 
809     return new SymOp(rot, tr);
810   }
811 
812   /**
813    * trtoString
814    *
815    * @param tr a double.
816    * @return a {@link java.lang.String} object.
817    */
818   private static String trtoString(double tr) {
819     if (tr == f12) {
820       return "+1/2";
821     }
822     if (tr == f13) {
823       return "+1/3";
824     }
825     if (tr == f23) {
826       return "+2/3";
827     }
828     if (tr == f14) {
829       return "+1/4";
830     }
831     if (tr == f34) {
832       return "+3/4";
833     }
834     if (tr == f16) {
835       return "+1/6";
836     }
837     if (tr == f56) {
838       return "+5/6";
839     }
840 
841     return "";
842   }
843 
844   /** {@inheritDoc} */
845   @Override
846   public String toString() {
847     StringBuilder sb = new StringBuilder(" Rotation operator:\n");
848     sb.append(
849         format(
850             " [[%4.1f,%4.1f,%4.1f]\n  [%4.1f,%4.1f,%4.1f]\n  [%4.1f,%4.1f,%4.1f]]\n",
851             rot[0][0], rot[0][1], rot[0][2], rot[1][0], rot[1][1], rot[1][2], rot[2][0], rot[2][1],
852             rot[2][2]));
853     sb.append(" Translation:\n");
854     sb.append(format(" [%4.2f,%4.2f,%4.2f]", tr[0], tr[1], tr[2]));
855     return sb.toString();
856   }
857 
858   /**
859    * Print the symmetry operator with double precision.
860    *
861    * @return String of rotation/translation with double precision.
862    */
863   public String toStringPrecise() {
864     StringBuilder sb = new StringBuilder(" Rotation operator:\n");
865     sb.append(
866         format(
867             " [[%18.16e,%18.16e,%18.16e]\n  [%18.16e,%18.16e,%18.16e]\n  [%18.16e,%18.16e,%18.16e]]\n",
868             rot[0][0], rot[0][1], rot[0][2], rot[1][0], rot[1][1], rot[1][2], rot[2][0], rot[2][1],
869             rot[2][2]));
870     sb.append(" Translation:\n");
871     sb.append(format(" [%18.16e,%18.16e,%18.16e]", tr[0], tr[1], tr[2]));
872     return sb.toString();
873   }
874 
875   /**
876    * toXYZString
877    *
878    * @return a {@link java.lang.String} object.
879    */
880   public String toXYZString() {
881     StringBuilder sb = new StringBuilder();
882     for (int i = 0; i < 3; i++) {
883       boolean s = false;
884       if (rot[i][0] < 0.0) {
885         sb.append("-X");
886         s = true;
887       } else if (rot[i][0] > 0.0) {
888         sb.append("X");
889         s = true;
890       }
891       if (rot[i][1] < 0.0) {
892         sb.append("-Y");
893         s = true;
894       } else if (rot[i][1] > 0.0) {
895         sb.append(s ? "+Y" : "Y");
896         s = true;
897       }
898       if (rot[i][2] < 0.0) {
899         sb.append("-Z");
900       } else if (rot[i][2] > 0.0) {
901         sb.append(s ? "+Z" : "Z");
902       }
903 
904       if (tr[i] > 0.0) {
905         sb.append(trtoString(tr[i]));
906       }
907 
908       if (i < 2) {
909         sb.append(", ");
910       }
911     }
912     return sb.toString();
913   }
914 }