1 /* $RCSfile$
2  * $Author: hansonr $
3  * $Date: 2007-05-18 08:19:45 -0500 (Fri, 18 May 2007) $
4  * $Revision: 7742 $
5 
6  *
7  * Copyright (C) 2003-2005  The Jmol Development Team
8  *
9  * Contact: jmol-developers@lists.sf.net
10  *
11  *  This library is free software; you can redistribute it and/or
12  *  modify it under the terms of the GNU Lesser General Public
13  *  License as published by the Free Software Foundation; either
14  *  version 2.1 of the License, or (at your option) any later version.
15  *
16  *  This library is distributed in the hope that it will be useful,
17  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
18  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
19  *  Lesser General Public License for more details.
20  *
21  *  You should have received a copy of the GNU Lesser General Public
22  *  License along with this library; if not, write to the Free Software
23  *  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA.
24  */
25 
26 package org.jmol.symmetry;
27 
28 import java.util.Hashtable;
29 import java.util.Map;
30 
31 import javajs.util.Lst;
32 import javajs.util.M4;
33 import javajs.util.Measure;
34 import javajs.util.P3;
35 import javajs.util.P4;
36 import javajs.util.PT;
37 import javajs.util.Quat;
38 import javajs.util.SB;
39 import javajs.util.T3;
40 import javajs.util.V3;
41 
42 import org.jmol.api.SymmetryInterface;
43 import javajs.util.BS;
44 import org.jmol.modelset.Atom;
45 import org.jmol.modelset.ModelSet;
46 import org.jmol.script.T;
47 import org.jmol.util.Escape;
48 import org.jmol.util.Logger;
49 
50 /**
51  * A class to handle requests for information about space groups and symmetry
52  * operations.
53  *
54  * Two entry points, both from Symmetry:
55  *
56  * getSymopInfo
57  *
58  * getSpaceGroupInfo
59  *
60  *
61  *
62  */
63 public class SymmetryDesc {
64 
65   private ModelSet modelSet;
66 
SymmetryDesc()67   public SymmetryDesc() {
68     // for reflection
69   }
70 
set(ModelSet modelSet)71   public SymmetryDesc set(ModelSet modelSet) {
72     this.modelSet = modelSet;
73     return this;
74   }
75 
76   private final static String[] keys = { "xyz", "xyzOriginal", "label",
77       null /*draw*/, "fractionalTranslation", "cartesianTranslation",
78       "inversionCenter", null /*point*/, "axisVector", "rotationAngle",
79       "matrix", "unitTranslation", "centeringVector", "timeReversal", "plane",
80       "_type", "id", "cif2", "xyzCanonical" };
81 
82   private final static int KEY_DRAW = 3;
83   private final static int KEY_POINT = 7;
84 
85   ////// "public" methods ////////
86 
87   /**
88    *
89    * @param iAtom
90    * @param xyz
91    * @param op
92    * @param translation TODO
93    * @param pt
94    * @param pt2
95    * @param id
96    * @param type
97    * @param scaleFactor
98    * @param nth
99    * @param options 0 or T.offset
100    * @return "" or a bitset of matching atoms, or
101    */
getSymopInfo(int iAtom, String xyz, int op, P3 translation, P3 pt, P3 pt2, String id, int type, float scaleFactor, int nth, int options)102   Object getSymopInfo(int iAtom, String xyz, int op, P3 translation, P3 pt, P3 pt2,
103                       String id, int type, float scaleFactor, int nth, int options) {
104     if (type == 0)
105       type = getType(id);
106     Object ret = (type == T.atoms ? new BS() : "");
107     if (iAtom < 0)
108       return ret;
109 
110     // get model symmetry
111 
112     int iModel = modelSet.at[iAtom].mi;
113     SymmetryInterface uc = modelSet.am[iModel].biosymmetry;
114     if (uc == null && (uc = modelSet.getUnitCell(iModel)) == null)
115       return ret;
116 
117     // generally get the result from getSymmetryInfo
118 
119     if (type != T.draw || op != Integer.MAX_VALUE) {
120       return getSymmetryInfo((Symmetry) uc, iModel, iAtom, (Symmetry) uc, xyz,
121           op, translation, pt, pt2, id, type, scaleFactor, nth, options);
122     }
123 
124     // draw SPACEGROUP
125 
126     String s = "";
127     M4[] ops = uc.getSymmetryOperations();
128     if (ops != null) {
129       if (id == null)
130         id = "sg";
131       int n = ops.length;
132       for (op = 0; op < n; op++)
133         s += (String) getSymmetryInfo((Symmetry) uc, iModel, iAtom,
134             (Symmetry) uc, xyz, op, translation, pt, pt2, id + op, T.draw, scaleFactor, nth, options);
135     }
136     return s;
137   }
138 
139   @SuppressWarnings("unchecked")
getSpaceGroupInfo(Symmetry sym, int modelIndex, String sgName, int symOp, P3 pt1, P3 pt2, String drawID, float scaleFactor, int nth, boolean isFull, boolean isForModel, int options, SymmetryInterface cellInfo)140   Map<String, Object> getSpaceGroupInfo(Symmetry sym, int modelIndex,
141                                         String sgName, int symOp, P3 pt1,
142                                         P3 pt2, String drawID,
143                                         float scaleFactor, int nth,
144                                         boolean isFull, boolean isForModel, int options, SymmetryInterface cellInfo) {
145     Map<String, Object> info = null;
146     boolean isStandard = (pt1 == null && drawID == null && nth <= 0);
147     boolean isBio = false;
148     String sgNote = null;
149     boolean haveName = (sgName != null && sgName.length() > 0);
150     boolean haveRawName = (haveName && sgName.indexOf("[--]") >= 0);
151     if (isForModel || !haveName) {
152       boolean saveModelInfo = (isStandard && symOp == 0);
153       if (modelIndex < 0)
154         modelIndex = (pt1 instanceof Atom ? ((Atom) pt1).mi
155             : modelSet.vwr.am.cmi);
156       if (modelIndex < 0)
157         sgNote = "no single current model";
158       else if (cellInfo == null && !(isBio = (cellInfo = modelSet.am[modelIndex].biosymmetry) != null)
159           && (cellInfo = modelSet.getUnitCell(modelIndex)) == null)
160         sgNote = "not applicable";
161       if (sgNote != null) {
162         info = new Hashtable<String, Object>();
163         info.put("spaceGroupInfo", "");
164         info.put("spaceGroupNote", sgNote);
165         info.put("symmetryInfo", "");
166       } else if (isStandard) {
167         info = (Map<String, Object>) modelSet.getInfo(modelIndex,
168             "spaceGroupInfo");
169       }
170       // created once
171       if (info != null)
172         return info;
173 
174       // full check
175       info = new Hashtable<String, Object>();
176       sgName = cellInfo.getSpaceGroupName();
177       SymmetryOperation[] ops = (SymmetryOperation[]) cellInfo
178           .getSymmetryOperations();
179       SpaceGroup sg = (isBio ? ((Symmetry) cellInfo).spaceGroup : null);
180       String slist = (haveRawName ? "" : null);
181       int opCount = 0;
182       if (ops != null) {
183         if (isBio)
184           sym.spaceGroup = SpaceGroup.getNull(false, false, false);
185         else
186           sym.setSpaceGroup(false);
187         // check to make sure that new group has been created magnetic or not
188         if (ops[0].timeReversal != 0)
189           ((SymmetryOperation) sym.getSpaceGroupOperation(0)).timeReversal = 1;
190         Object[][] infolist = new Object[ops.length][];
191         String sops = "";
192         for (int i = 0, nop = 0; i < ops.length && nop != nth; i++) {
193 
194           SymmetryOperation op = ops[i];
195           boolean isNewIncomm = (i == 0 && op.xyz.indexOf("x4") >= 0);
196           int iop = (!isNewIncomm && sym.getSpaceGroupOperation(i) != null ? i : isBio ? sym
197               .addBioMoleculeOperation(sg.finalOperations[i], false) : sym
198               .addSpaceGroupOperation("=" + op.xyz, i + 1));
199           if (iop < 0)
200             continue;
201           String xyzOriginal = op.xyzOriginal;
202           op = (SymmetryOperation) sym.getSpaceGroupOperation(i);
203           if (op == null)
204             continue;
205           op.xyzOriginal = xyzOriginal;
206           if (op.timeReversal != 0 || op.modDim > 0)
207             isStandard = false;
208           if (slist != null)
209             slist += ";" + op.xyz;
210           Object[] ret = (symOp > 0 && symOp - 1 != iop ? null
211               : createInfoArray(op, cellInfo, pt1, pt2, drawID, scaleFactor, options, false));
212           if (ret != null) {
213             if (nth > 0 && ++nop != nth)
214               continue;
215             infolist[i] = ret;
216             sops += "\n" + (i + 1) + "\t" + ret[0] + "\t" + ret[2];
217             opCount++;
218           }
219         }
220         info.put("operations", infolist);
221         info.put("symmetryInfo", (sops.length() == 0 ? "" : sops.substring(1)));
222       }
223       sgNote = (opCount == 0 ? "\n no symmetry operations"
224           : nth <= 0 && symOp <= 0 ? "\n" + opCount
225             + " symmetry operation" + (opCount == 1 ? ":\n" : "s:\n") : "");
226       if (slist != null)
227         sgName = slist.substring(slist.indexOf(";") + 1);
228       if (saveModelInfo)
229         modelSet.setInfo(modelIndex, "spaceGroupInfo", info);
230     } else {
231       info = new Hashtable<String, Object>();
232     }
233     info.put("spaceGroupName", sgName);
234     info.put("spaceGroupNote", sgNote == null ? "" : sgNote);
235     Object data;
236     if (isBio) {
237       data = sgName;
238     } else {
239       if (haveName && !haveRawName)
240         sym.setSpaceGroupName(sgName);
241       data = sym.getSpaceGroupInfoObj(sgName, cellInfo, isFull);
242       if (data == null || data.equals("?")) {
243         data = "?";
244         info.put("spaceGroupNote", "could not identify space group from name: " + sgName
245             + "\nformat: show spacegroup \"2\" or \"P 2c\" "
246             + "or \"C m m m\" or \"x, y, z;-x ,-y, -z\"");
247       }
248     }
249     info.put("spaceGroupInfo", data);
250     return info;
251   }
252 
253   //////////// private methods ///////////
254 
255   /** Determine the type of this request.
256    * Note that label and xyz will be returned as T.xys and T.label
257    *
258    * @param id
259    * @return a code that identifies this request.
260    */
getType(String id)261   private static int getType(String id) {
262     int type;
263     if (id == null)
264       return T.list;
265     if (id.equalsIgnoreCase("matrix"))
266       return T.matrix4f;
267     if (id.equalsIgnoreCase("description"))
268       return T.label;
269     if (id.equalsIgnoreCase("axispoint"))
270       return T.point;
271     if (id.equalsIgnoreCase("time"))
272       return T.times;
273     if (id.equalsIgnoreCase("info"))
274       return T.array;
275     // center, draw, plane, axis, atom, translation, angle, array, list
276     type = T.getTokFromName(id);
277     if (type != 0)
278       return type;
279     type = getKeyType(id);
280     return (type < 0 ? type : T.all);
281   }
282 
getKeyType(String id)283   private static int getKeyType(String id) {
284     if ("type".equals(id))
285       id = "_type";
286     for (int type = 0; type < keys.length; type++)
287       if (id.equalsIgnoreCase(keys[type]))
288         return -1 - type;
289     return 0;
290   }
291 
nullReturn(int type)292   private static Object nullReturn(int type) {
293     switch(type) {
294     case T.draw:
295       return "draw ID sym_* delete";
296     case T.full:
297     case T.label:
298     case T.id:
299     case T.xyz:
300     case T.matrix3f:
301     case T.origin:
302       return "";
303     case T.atoms:
304       return new BS();
305     default:
306       return null;
307     }
308   }
309 
310   /**
311    * Return information about a symmetry operator by type:
312    *
313    * array, angle, axis, center, draw, full, info, label, matrix4f, point, time,
314    * plane, translation, unitcell, xyz, all,
315    *
316    * or a negative number (-length, -1]:
317    *
318    * { "xyz", etc. }
319    *
320    * where "all" is the info array itself,
321    *
322    * @param info
323    * @param type
324    * @return object specified
325    *
326    */
327 
getInfo(Object[] info, int type)328   private static Object getInfo(Object[] info, int type) {
329     if (info.length == 0)
330       return "";
331     if (type < 0 && -type <= keys.length && -type <= info.length)
332       return info[-1 - type];
333     switch (type) {
334     case T.all:
335     case T.info:
336       return info; // 'info' is internal only; if user selects "info", it is turned into "array"
337     case T.array:
338       Map<String, Object> lst = new Hashtable<String, Object>();
339       for (int j = 0, n = info.length; j < n; j++) {
340         String key = (j == KEY_DRAW ? "draw" : j == KEY_POINT ? "axispoint" : keys[j]);
341         if (info[j] != null)
342           lst.put(key, info[j]);
343       }
344       return lst;
345     case T.full:
346       return info[0] + "  \t" + info[2];
347     case T.xyz:
348       return info[0];
349     case T.origin:
350       return info[1];
351     default:
352     case T.label:
353       return info[2];
354     case T.draw:
355       return info[3] + "\nprint " + PT.esc(info[0] + " " + info[2]);
356     case T.fracxyz:
357       return info[4]; // fractional translation "fxyz"
358     case T.translation:
359       return info[5]; // cartesian translation
360     case T.center:
361       return info[6];
362     case T.point:
363       return info[7];
364     case T.axis:
365       return info[8];
366     case T.angle:
367       return info[9];
368     case T.matrix4f:
369       return info[10];
370     case T.unitcell:
371       return info[11];
372     case T.translate:
373       // centering
374       return info[12];
375     case T.times:
376       return info[13];
377     case T.plane:
378       return info[14];
379     case T.type:
380       return info[15];
381     case T.id:
382       return info[16];
383     }
384   }
385 
386   /**
387    *
388    * @param op
389    * @param uc
390    * @param pta00
391    *        optional initial atom point
392    * @param ptTarget
393    *        optional target atom point
394    * @param id
395    * @param scaleFactor
396    *        scale for rotation vector only
397    * @param options
398    *        0 or T.offset
399    * @param haveTranslation TODO
400    * @return Object[] containing:
401    *
402    *         [0] xyz (Jones-Faithful calculated from matrix)
403    *
404    *         [1] xyzOriginal (Provided by calling method)
405    *
406    *         [2] info ("C2 axis", for example)
407    *
408    *         [3] draw commands
409    *
410    *         [4] translation vector (fractional)
411    *
412    *         [5] translation vector (Cartesian)
413    *
414    *         [6] inversion point
415    *
416    *         [7] axis point
417    *
418    *         [8] axis vector (defines plane if angle = 0
419    *
420    *         [9] angle of rotation
421    *
422    *         [10] matrix representation
423    *
424    *         [11] lattice translation
425    *
426    *         [12] centering
427    *
428    *         [13] time reversal
429    *
430    *         [14] plane
431    *
432    *         [15] _type
433    *
434    *         [16] id
435    *
436    */
437 
createInfoArray(SymmetryOperation op, SymmetryInterface uc, P3 pta00, P3 ptTarget, String id, float scaleFactor, int options, boolean haveTranslation)438   private Object[] createInfoArray(SymmetryOperation op, SymmetryInterface uc,
439                                    P3 pta00, P3 ptTarget, String id,
440                                    float scaleFactor, int options, boolean haveTranslation) {
441     if (!op.isFinalized)
442       op.doFinalize();
443     boolean isTimeReversed = (op.timeReversal == -1);
444     if (scaleFactor == 0)
445       scaleFactor = 1f;
446     V3 vtemp = new V3();
447     P3 ptemp = new P3();
448     P3 ptemp2 = new P3();
449     P3 pta01 = new P3();
450     P3 pta02 = new P3();
451     V3 ftrans = new V3();
452     V3 vtrans = new V3();
453     P4 plane = null;
454 
455     if (pta00 == null || Float.isNaN(pta00.x))
456       pta00 = new P3();
457     if (ptTarget != null) {
458 
459       // Check to see that this is the correct operator
460 
461       setFractional(uc, pta00, pta01, ptemp);
462       op.rotTrans(pta01);
463       uc.toCartesian(pta01, false);
464       uc.toUnitCell(pta01, ptemp);
465       pta02.setT(ptTarget);
466       uc.toUnitCell(pta02, ptemp);
467       if (pta01.distance(pta02) > 0.1f)
468         return null;
469 
470       // Check to see if the two points only differ by
471       // a translation after transformation.
472       // If so, add that difference to the matrix transformation
473 
474       setFractional(uc, pta00, pta01, null);
475       op.rotTrans(pta01);
476       setFractional(uc, ptTarget, pta02, null);
477       vtrans.sub2(pta02, pta01);
478     }
479 
480     // get the frame vectors and points
481 
482     pta01.set(1, 0, 0);
483     pta02.set(0, 1, 0);
484     P3 pta03 = P3.new3(0, 0, 1);
485     pta01.add(pta00);
486     pta02.add(pta00);
487     pta03.add(pta00);
488 
489     // target point, rotated, inverted, and translated
490 
491     P3 pt0 = rotTransCart(op, uc, pta00, vtrans);
492     P3 pt1 = rotTransCart(op, uc, pta01, vtrans);
493     P3 pt2 = rotTransCart(op, uc, pta02, vtrans);
494     P3 pt3 = rotTransCart(op, uc, pta03, vtrans);
495 
496     V3 vt1 = V3.newVsub(pt1, pt0);
497     V3 vt2 = V3.newVsub(pt2, pt0);
498     V3 vt3 = V3.newVsub(pt3, pt0);
499 
500     approx(vtrans);
501 
502     // check for inversion
503 
504     vtemp.cross(vt1, vt2);
505     boolean haveInversion = (vtemp.dot(vt3) < 0);
506 
507     // The first trick is to check cross products to see if we still have a
508     // right-hand axis.
509 
510     if (haveInversion) {
511 
512       // undo inversion for quaternion analysis (requires proper rotations only)
513 
514       pt1.sub2(pt0, vt1);
515       pt2.sub2(pt0, vt2);
516       pt3.sub2(pt0, vt3);
517 
518     }
519 
520     // The second trick is to use quaternions. Each of the three faces of the
521     // frame (xy, yz, and zx)
522     // is checked. The helix() function will return data about the local helical
523     // axis, and the
524     // symop(sym,{0 0 0}) function will return the overall translation.
525 
526     T3[] info = Measure.computeHelicalAxis(pta00, pt0,
527         Quat.getQuaternionFrame(pt0, pt1, pt2)
528             .div(Quat.getQuaternionFrame(pta00, pta01, pta02)));
529     // new T3[] { pt_a_prime, n, r, P3.new3(theta, pitch, residuesPerTurn), pt_b_prime };
530     P3 pa1 = (P3) info[0];
531     V3 ax1 = (V3) info[1];
532     int ang1 = (int) Math.abs(PT.approx(((P3) info[3]).x, 1));
533     float pitch1 = SymmetryOperation.approxF(((P3) info[3]).y);
534 
535     if (haveInversion) {
536 
537       // redo inversion
538 
539       pt1.add2(pt0, vt1);
540       pt2.add2(pt0, vt2);
541       pt3.add2(pt0, vt3);
542 
543     }
544 
545     V3 trans = V3.newVsub(pt0, pta00);
546     if (trans.length() < 0.1f)
547       trans = null;
548 
549     // ////////// determination of type of operation from first principles
550 
551     P3 ptinv = null; // inverted point for translucent frame
552     P3 ipt = null; // inversion center
553     P3 ptref = null; // reflection center
554 
555     boolean isTranslation = (ang1 == 0);
556     boolean isRotation = !isTranslation;
557     boolean isInversionOnly = false;
558     boolean isMirrorPlane = false;
559 
560     if (isRotation || haveInversion) {
561       // the translation will be taken care of other ways
562       trans = null;
563     }
564 
565     // handle inversion
566 
567     if (haveInversion && isTranslation) {
568 
569       // simple inversion operation
570 
571       ipt = P3.newP(pta00);
572       ipt.add(pt0);
573       ipt.scale(0.5f);
574       ptinv = pt0;
575       isInversionOnly = true;
576 
577     } else if (haveInversion) {
578 
579       /*
580        *
581        * We must convert simple rotations to rotation-inversions; 2-screws to
582        * planes and glide planes.
583        *
584        * The idea here is that there is a relationship between the axis for a
585        * simple or screw rotation of an inverted frame and one for a
586        * rotation-inversion. The relationship involves two adjacent equilateral
587        * triangles:
588        *
589        *
590        *       o
591        *      / \
592        *     /   \    i'
593        *    /     \
594        *   /   i   \
595        * A/_________\A'
596        *  \         /
597        *   \   j   /
598        *    \     /
599        *     \   /
600        *      \ /
601        *       x
602        *
603        * Points i and j are at the centers of the triangles. Points A and A' are
604        * the frame centers; an operation at point i, j, x, or o is taking A to
605        * A'. Point i is 2/3 of the way from x to o. In addition, point j is half
606        * way between i and x.
607        *
608        * The question is this: Say you have an rotation/inversion taking A to
609        * A'. The relationships are:
610        *
611        * 6-fold screw x for inverted frame corresponds to 6-bar at i for actual
612        * frame 3-fold screw i for inverted frame corresponds to 3-bar at x for
613        * actual frame
614        *
615        * The proof of this follows. Consider point x. Point x can transform A to
616        * A' as a clockwise 6-fold screw axis. So, say we get that result for the
617        * inverted frame. What we need for the real frame is a 6-bar axis
618        * instead. Remember, though, that we inverted the frame at A to get this
619        * result. The real frame isn't inverted. The 6-bar must do that inversion
620        * AND also get the frame to point A' with the same (clockwise) rotation.
621        * The key is to see that there is another axis -- at point i -- that does
622        * the trick.
623        *
624        * Take a look at the angles and distances that arise when you project A
625        * through point i. The result is a frame at i'. Since the distance i-i'
626        * is the same as i-A (and thus i-A') and the angle i'-i-A' is 60 degrees,
627        * point i is also a 6-bar axis transforming A to A'.
628        *
629        * Note that both the 6-fold screw axis at x and the 6-bar axis at i are
630        * both clockwise.
631        *
632        * Similar analysis shows that the 3-fold screw i corresponds to the 3-bar
633        * axis at x.
634        *
635        * So in each case we just calculate the vector i-j or x-o and then factor
636        * appropriately.
637        *
638        * The 4-fold case is simpler -- just a parallelogram.
639        */
640 
641       V3 d = (pitch1 == 0 ? new V3() : ax1);
642       float f = 0;
643       switch (ang1) {
644       case 60: // 6_1 at x to 6-bar at i
645         f = 2f / 3f;
646         break;
647       case 120: // 3_1 at i to 3-bar at x
648         f = 2;
649         break;
650       case 90: // 4_1 to 4-bar at opposite corner
651         f = 1;
652         break;
653       case 180: // 2_1 to mirror plane
654         // C2 with inversion is a mirror plane -- but could have a glide
655         // component.
656         ptref = P3.newP(pta00);
657         ptref.add(d);
658         pa1.scaleAdd2(0.5f, d, pta00);
659         if (ptref.distance(pt0) > 0.1f) {
660           trans = V3.newVsub(pt0, ptref);
661           setFractional(uc, trans, ptemp, null);
662           ftrans.setT(ptemp);
663         } else {
664           trans = null;
665         }
666         isRotation = false;
667         haveInversion = false;
668         isMirrorPlane = true;
669         break;
670       default:
671         haveInversion = false;
672         break;
673       }
674       if (f != 0) {
675         vtemp.sub2(pta00, pa1);
676         vtemp.add(pt0);
677         vtemp.sub(pa1);
678         vtemp.sub(d);
679         vtemp.scale(f);
680         pa1.add(vtemp);
681         ipt = new P3();
682         ipt.scaleAdd2(0.5f, d, pa1);
683         ptinv = new P3();
684         ptinv.scaleAdd2(-2, ipt, pta00);
685         ptinv.scale(-1);
686       }
687 
688     } else if (trans != null) {
689 
690       // get rid of unnecessary translations added to keep most operations
691       // within cell 555
692 
693       ptemp.setT(trans);
694       uc.toFractional(ptemp, false);
695       //      if (SymmetryOperation.approxF(ptemp.x) == 1)
696       //        ptemp.x = 0;
697       //      if (SymmetryOperation.approxF(ptemp.y) == 1)
698       //        ptemp.y = 0;
699       //      if (SymmetryOperation.approxF(ptemp.z) == 1)
700       //        ptemp.z = 0;
701 
702       ftrans.setT(ptemp);
703       uc.toCartesian(ptemp, false);
704       trans.setT(ptemp);
705     }
706 
707     // fix angle based on direction of axis
708 
709     int ang = ang1;
710     approx0(ax1);
711 
712     if (isRotation) {
713 
714       P3 ptr = new P3();
715 
716       vtemp.setT(ax1);
717 
718       // draw the lines associated with a rotation
719 
720       int ang2 = ang1;
721       if (haveInversion) {
722         ptr.add2(pa1, vtemp);
723         ang2 = Math.round(Measure.computeTorsion(ptinv, pa1, ptr, pt0, true));
724       } else if (pitch1 == 0) {
725         ptr.setT(pa1);
726         ptemp.scaleAdd2(1, ptr, vtemp);
727         ang2 = Math.round(Measure.computeTorsion(pta00, pa1, ptemp, pt0, true));
728       } else {
729         ptemp.add2(pa1, vtemp);
730         ptr.scaleAdd2(0.5f, vtemp, pa1);
731         ang2 = Math.round(Measure.computeTorsion(pta00, pa1, ptemp, pt0, true));
732       }
733 
734       if (ang2 != 0)
735         ang1 = ang2;
736 
737       if (!haveInversion && pitch1 == 0 && (ax1.z < 0
738           || ax1.z == 0 && (ax1.y < 0 || ax1.y == 0 && ax1.x < 0))) {
739         ax1.scale(-1);
740         ang1 = -ang1;
741       }
742 
743     }
744 
745     // time to get the description
746 
747     String info1 = "identity";
748     String type = info1;
749 
750     if (isInversionOnly) {
751       ptemp.setT(ipt);
752       uc.toFractional(ptemp, false);
753       info1 = "Ci: " + strCoord(op, ptemp, op.isBio);
754       type = "inversion center";
755     } else if (isRotation) {
756       if (haveInversion) {
757         type = info1 = (360 / ang) + "-bar axis";
758       } else if (pitch1 != 0) {
759         type = info1 = (360 / ang) + "-fold screw axis";
760         ptemp.setT(ax1);
761         uc.toFractional(ptemp, false);
762         info1 += "|translation: " + strCoord(op, ptemp, op.isBio);
763       } else {
764         type = info1 = "C" + (360 / ang) + " axis";
765       }
766     } else if (trans != null) {
767       String s = " " + strCoord(op, ftrans, op.isBio);
768       if (isTranslation) {
769         type = info1 = "translation";
770         info1 += ":" + s;
771       } else if (isMirrorPlane) {
772         float fx = Math.abs(SymmetryOperation.approxF(ftrans.x));
773         float fy = Math.abs(SymmetryOperation.approxF(ftrans.y));
774         float fz = Math.abs(SymmetryOperation.approxF(ftrans.z));
775         s = " " + strCoord(op, ftrans, op.isBio);
776         // set ITA Table 2.1.2.1
777         if (fx != 0 && fy != 0 && fz != 0) {
778           if (fx == 1 / 4f && fy == 1 / 4f && fz == 1 / 4f) {
779             // diamond
780             info1 = "d-";
781           } else if (fx == 1 / 2f && fy == 1 / 2f && fz == 1 / 2f) {
782             info1 = "n-";
783           } else {
784             info1 = "g-";
785           }
786         } else if (fx != 0 && fy != 0 || fy != 0 && fz != 0
787             || fz != 0 && fx != 0) {
788           // any two
789           if (fx == 1 / 4f && fy == 1 / 4f || fx == 1 / 4f && fz == 1 / 4f
790               || fy == 1 / 4f && fz == 1 / 4f) {
791             info1 = "d-";
792           } else if (fx == 1 / 2f && fy == 1 / 2f
793               || fx == 1 / 2f && fz == 1 / 2f || fy == 1 / 2f && fz == 1 / 2f) {
794             // making sure here that this is truly a diagonal in the plane, not just
795             // a glide parallel to a face on a diagonal plane! Mois Aroyo 2018
796             if (fx == 0 && ax1.x == 0 || fy == 0 && ax1.y == 0
797                 || fz == 0 && ax1.z == 0) {
798               info1 = "g-";
799             } else {
800               info1 = "n-";
801             }
802           } else {
803             info1 = "g-";
804           }
805         } else if (fx != 0)
806           info1 = "a-";
807         else if (fy != 0)
808           info1 = "b-";
809         else
810           info1 = "c-";
811         type = info1 = info1 + "glide plane";
812         info1 += "|translation:" + s;
813       }
814     } else if (isMirrorPlane) {
815       type = info1 = "mirror plane";
816     }
817 
818     if (haveInversion && !isInversionOnly) {
819       ptemp.setT(ipt);
820       uc.toFractional(ptemp, false);
821       info1 += "|inversion center at " + strCoord(op, ptemp, op.isBio);
822     }
823 
824     if (isTimeReversed) {
825       info1 += "|time-reversed";
826       type += " (time-reversed)";
827     }
828 
829     String cmds = null;
830     String xyzNew = (op.isBio ? op.xyzOriginal
831         : SymmetryOperation.getXYZFromMatrix(op, false, false, false));
832 
833     // check for drawing
834 
835     if (id != null) {
836 
837       String opType = null;
838       String drawid = "\ndraw ID " + id + "_";
839 
840       // delete previous elements of this user-settable ID
841 
842       SB draw1 = new SB();
843 
844       draw1.append(drawid).append("* delete");
845       //    .append(
846       //    ("print " + PT.esc(
847       //        id + " " + (op.index + 1) + " " + op.fixMagneticXYZ(op, op.xyzOriginal, false) + "|"
848       //            + op.fixMagneticXYZ(op, xyzNew, true) + "|" + info1).replace(
849       //        '\n', ' '))).append("\n")
850 
851       // draw the initial frame
852 
853       drawLine(draw1, drawid + "frame1X", 0.15f, pta00, pta01, "red");
854       drawLine(draw1, drawid + "frame1Y", 0.15f, pta00, pta02, "green");
855       drawLine(draw1, drawid + "frame1Z", 0.15f, pta00, pta03, "blue");
856 
857       String color;
858 
859       if (isRotation) {
860 
861         P3 ptr = new P3();
862 
863         color = "red";
864 
865         ang = ang1;
866         float scale = 1f;
867         vtemp.setT(ax1);
868 
869         // draw the lines associated with a rotation
870 
871         if (haveInversion) {
872           opType = drawid + "rotinv";
873           ptr.add2(pa1, vtemp);
874           if (pitch1 == 0) {
875             ptr.setT(ipt);
876             vtemp.scale(3 * scaleFactor);
877             ptemp.scaleAdd2(-1, vtemp, pa1);
878             drawVector(draw1, drawid, "rotVector2", "", pa1, ptemp, "red");
879           }
880           scale = pt0.distance(ptr);
881           draw1.append(drawid).append("rotLine1 ").append(Escape.eP(ptr))
882               .append(Escape.eP(ptinv)).append(" color red");
883           draw1.append(drawid).append("rotLine2 ").append(Escape.eP(ptr))
884               .append(Escape.eP(pt0)).append(" color red");
885         } else if (pitch1 == 0) {
886           opType = drawid + "rot";
887           boolean isSpecial = (pta00.distance(pt0) < 0.2f);
888           if (!isSpecial) {
889             draw1.append(drawid).append("rotLine1 ").append(Escape.eP(pta00))
890                 .append(Escape.eP(pa1)).append(" color red");
891             draw1.append(drawid).append("rotLine2 ").append(Escape.eP(pt0))
892                 .append(Escape.eP(pa1)).append(" color red");
893           }
894           vtemp.scale(3 * scaleFactor);
895           ptemp.scaleAdd2(-1, vtemp, pa1);
896           drawVector(draw1, drawid, "rotVector2", "", pa1, ptemp,
897               isTimeReversed ? "gray" : "red");
898           ptr.setT(pa1);
899           if (pitch1 == 0 && pta00.distance(pt0) < 0.2)
900             ptr.scaleAdd2(0.5f, vtemp, ptr);
901         } else {
902           opType = drawid + "screw";
903           color = "orange";
904           draw1.append(drawid).append("rotLine1 ").append(Escape.eP(pta00))
905               .append(Escape.eP(pa1)).append(" color red");
906           ptemp.add2(pa1, vtemp);
907           draw1.append(drawid).append("rotLine2 ").append(Escape.eP(pt0))
908               .append(Escape.eP(ptemp)).append(" color red");
909           ptr.scaleAdd2(0.5f, vtemp, pa1);
910         }
911 
912         // draw arc arrow
913 
914         ptemp.add2(ptr, vtemp);
915         if (haveInversion && pitch1 != 0) {
916           draw1.append(drawid).append("rotRotLine1").append(Escape.eP(ptr))
917               .append(Escape.eP(ptinv)).append(" color red");
918           draw1.append(drawid).append("rotRotLine2").append(Escape.eP(ptr))
919               .append(Escape.eP(pt0)).append(" color red");
920         }
921         draw1.append(drawid)
922             .append(
923                 "rotRotArrow arrow width 0.1 scale " + PT.escF(scale) + " arc ")
924             .append(Escape.eP(ptr)).append(Escape.eP(ptemp));
925         ptemp.setT(haveInversion ? ptinv : pta00);
926         if (ptemp.distance(pt0) < 0.1f)
927           ptemp.set((float) Math.random(), (float) Math.random(),
928               (float) Math.random());
929         draw1.append(Escape.eP(ptemp));
930         ptemp.set(0, ang, 0);
931         draw1.append(Escape.eP(ptemp)).append(" color red");
932 
933         // draw the main vector
934 
935         drawVector(draw1, drawid, "rotVector1", "vector", pa1, vtemp,
936             isTimeReversed ? "Gray" : color);
937 
938       } else if (isMirrorPlane) {
939 
940         // lavender arrow across plane from pt00 to pt0
941 
942         ptemp.sub2(ptref, pta00);
943         if (pta00.distance(ptref) > 0.2)
944           drawVector(draw1, drawid, "planeVector", "vector", pta00, ptemp,
945               isTimeReversed ? "Gray" : "cyan");
946 
947         // faint inverted frame if mirror trans is not null
948 
949         opType = drawid + "plane";
950         if (trans != null) {
951           opType = drawid + "glide";
952           drawFrameLine("X", ptref, vt1, 0.15f, ptemp, draw1, opType, "red");
953           drawFrameLine("Y", ptref, vt2, 0.15f, ptemp, draw1, opType, "green");
954           drawFrameLine("Z", ptref, vt3, 0.15f, ptemp, draw1, opType, "blue");
955         }
956 
957         color = (trans == null ? "green" : "blue");
958 
959         // ok, now HERE's a good trick. We use the Marching Cubes
960         // algorithm to find the intersection points of a plane and the unit
961         // cell.
962         // We expand the unit cell by 5% in all directions just so we are
963         // guaranteed to get cutoffs.
964 
965         vtemp.setT(ax1);
966         vtemp.normalize();
967         // ax + by + cz + d = 0
968         // so if a point is in the plane, then N dot X = -d
969         float w = -vtemp.x * pa1.x - vtemp.y * pa1.y - vtemp.z * pa1.z;
970         plane = P4.new4(vtemp.x, vtemp.y, vtemp.z, w);
971         float margin = (Math.abs(w) < 0.01f && vtemp.x * vtemp.y > 0.4 ? 1.30f
972             : 1.05f);
973         // returns triangles and lines
974         Lst<Object> v = modelSet.vwr.getTriangulator().intersectPlane(plane,
975             uc.getCanonicalCopy(margin, true), 3);
976         if (v != null)
977           for (int i = v.size(); --i >= 0;) {
978             P3[] pts = (P3[]) v.get(i);
979             draw1.append(drawid).append("planep").appendI(i).append(" ")
980                 .append(Escape.eP(pts[0])).append(Escape.eP(pts[1]));
981             if (pts.length == 3)
982               draw1.append(Escape.eP(pts[2]));
983             draw1.append(" color translucent ").append(color);
984           }
985 
986         // and JUST in case that does not work, at least draw a circle
987 
988         if (v == null || v.size() == 0) {
989           ptemp.add2(pa1, ax1);
990           draw1.append(drawid).append("planeCircle scale 2.0 circle ")
991               .append(Escape.eP(pa1)).append(Escape.eP(ptemp))
992               .append(" color translucent ").append(color).append(" mesh fill");
993         }
994       }
995 
996       if (haveInversion) {
997         opType = drawid + "inv";
998         draw1.append(drawid).append("invPoint diameter 0.4 ")
999             .append(Escape.eP(ipt));
1000         ptemp.sub2(ptinv, pta00);
1001         drawVector(draw1, drawid, "invArrow", "vector", pta00, ptemp,
1002             isTimeReversed ? "gray" : "cyan");
1003         if (!isInversionOnly && options != T.offset) {
1004           // n-bar: draw a faint frame showing the inversion
1005           drawFrameLine("X", ptinv, vt1, 0.15f, ptemp, draw1, opType, "red");
1006           drawFrameLine("Y", ptinv, vt2, 0.15f, ptemp, draw1, opType, "green");
1007           drawFrameLine("Z", ptinv, vt3, 0.15f, ptemp, draw1, opType, "blue");
1008         }
1009       }
1010 
1011       // and display translation if still not {0 0 0}
1012 
1013       if (trans != null) {
1014         if (ptref == null)
1015           ptref = P3.newP(pta00);
1016         drawVector(draw1, drawid, "transVector", "vector", ptref, trans,
1017             isTimeReversed && !haveInversion && !isMirrorPlane && !isRotation
1018                 ? "darkGray"
1019                 : "gold");
1020       }
1021 
1022       // draw the final frame just a bit fatter and shorter, in case they
1023       // overlap
1024 
1025       ptemp2.setT(pt0);
1026       ptemp.sub2(pt1, pt0);
1027       ptemp.scaleAdd2(0.9f, ptemp, ptemp2);
1028       drawLine(draw1, drawid + "frame2X", 0.2f, ptemp2, ptemp, "red");
1029       ptemp.sub2(pt2, pt0);
1030       ptemp.scaleAdd2(0.9f, ptemp, ptemp2);
1031       drawLine(draw1, drawid + "frame2Y", 0.2f, ptemp2, ptemp, "green");
1032       ptemp.sub2(pt3, pt0);
1033       ptemp.scaleAdd2(0.9f, ptemp, ptemp2);
1034       drawLine(draw1, drawid + "frame2Z", 0.2f, ptemp2, ptemp, "purple");
1035 
1036       // color the targeted atoms opaque and add another frame if necessary
1037 
1038       draw1.append("\nsym_point = " + Escape.eP(pta00));
1039       draw1.append("\nvar p0 = " + Escape.eP(ptemp2));
1040       draw1.append(
1041           "\nvar set2 = within(0.2,p0);if(!set2){set2 = within(0.2,p0.uxyz.xyz)}");
1042       if (pta00 instanceof Atom)
1043         draw1.append("\n set2 &= {_" + ((Atom) pta00).getElementSymbol() + "}");
1044       draw1.append("\nsym_target = set2;if (set2) {");
1045       //      if (haveCentering)
1046       //      draw1.append(drawid).append(
1047       //        "cellOffsetVector arrow @p0 @set2 color grey");
1048       if (options != T.offset && ptTarget == null && !haveTranslation) {
1049         draw1.append(drawid)
1050             .append("offsetFrameX diameter 0.20 @{set2.xyz} @{set2.xyz + ")
1051             .append(Escape.eP(vt1)).append("*0.9} color red");
1052         draw1.append(drawid)
1053             .append("offsetFrameY diameter 0.20 @{set2.xyz} @{set2.xyz + ")
1054             .append(Escape.eP(vt2)).append("*0.9} color green");
1055         draw1.append(drawid)
1056             .append("offsetFrameZ diameter 0.20 @{set2.xyz} @{set2.xyz + ")
1057             .append(Escape.eP(vt3)).append("*0.9} color purple");
1058       }
1059       draw1.append("\n}\n");
1060       cmds = draw1.toString();
1061       if (Logger.debugging)
1062         Logger.info(cmds);
1063       draw1 = null;
1064       drawid = null;
1065     }
1066     if (trans == null)
1067       ftrans = null;
1068     if (isRotation) {
1069       if (haveInversion) {
1070       } else if (pitch1 == 0) {
1071       } else {
1072         // screw
1073         trans = V3.newV(ax1);
1074         ptemp.setT(trans);
1075         uc.toFractional(ptemp, false);
1076         ftrans = V3.newV(ptemp);
1077       }
1078     }
1079     if (isMirrorPlane) {
1080       ang1 = 0;
1081     }
1082     if (haveInversion) {
1083       if (isInversionOnly) {
1084         pa1 = null;
1085         ax1 = null;
1086         trans = null;
1087         ftrans = null;
1088       }
1089     } else if (isTranslation) {
1090       pa1 = null;
1091       ax1 = null;
1092     }
1093 
1094     // and display translation if still not {0 0 0}
1095     if (ax1 != null)
1096       ax1.normalize();
1097     M4 m2 = null;
1098     m2 = M4.newM4(op);
1099     if (vtrans.length() != 0) {
1100       m2.m03 += vtrans.x;
1101       m2.m13 += vtrans.y;
1102       m2.m23 += vtrans.z;
1103     }
1104     // TODO: here we need magnetic
1105     xyzNew = (op.isBio ? m2.toString()
1106         : op.modDim > 0 ? op.xyzOriginal
1107             : SymmetryOperation.getXYZFromMatrix(m2, false, false, false));
1108     if (op.timeReversal != 0)
1109       xyzNew = op.fixMagneticXYZ(m2, xyzNew, true);
1110     T3 f0 = approx0(ftrans);
1111     T3 cift = null;
1112     int cifi = (op.number < 0 ? 0 : op.number);
1113     if (!xyzNew.equals(op.xyzOriginal)) {
1114       if (op.number > 0) {
1115         M4 orig = SymmetryOperation.getMatrixFromXYZ(op.xyzOriginal);
1116         orig.sub(m2);
1117         cift = new P3();
1118         orig.getTranslation(cift);
1119       }
1120     }
1121     String cif2 = cifi + (cift == null ? " [0 0 0]"
1122         : " [" + (int) -cift.x + " " + (int) -cift.y + " " + (int) -cift.z + "]");
1123     return new Object[] { xyzNew, op.xyzOriginal, info1, cmds, f0,
1124         approx0(trans), approx0(ipt), approx0(pa1),
1125         (plane == null ? approx0(ax1) : null),
1126         (ang1 != 0 ? Integer.valueOf(ang1) : null), m2,
1127         (vtrans.lengthSquared() > 0 ? vtrans : null), op.getCentering(),
1128         Integer.valueOf(op.timeReversal), plane, type,
1129         Integer.valueOf(op.number), cif2, op.xyzCanonical };
1130   }
1131 
drawLine(SB s, String id, float diameter, P3 pt0, P3 pt1, String color)1132   private static void drawLine(SB s, String id, float diameter, P3 pt0, P3 pt1,
1133                                String color) {
1134     s.append(id).append(" diameter ").appendF(diameter).append(Escape.eP(pt0))
1135         .append(Escape.eP(pt1)).append(" color ").append(color);
1136   }
1137 
drawFrameLine(String xyz, P3 pt, V3 v, float width, P3 ptemp, SB draw1, String key, String color)1138   private static void drawFrameLine(String xyz, P3 pt, V3 v, float width,
1139                                     P3 ptemp, SB draw1, String key, String color) {
1140     ptemp.setT(pt);
1141     ptemp.add(v);
1142     drawLine(draw1, key + "Pt" + xyz, width, pt, ptemp, "translucent " + color);
1143   }
1144 
drawVector(SB draw1, String drawid, String label, String type, T3 pt1, T3 v, String color)1145   private static void drawVector(SB draw1, String drawid, String label,
1146                                  String type, T3 pt1, T3 v, String color) {
1147     draw1.append(drawid).append(label).append(" diameter 0.1 ").append(type)
1148         .append(Escape.eP(pt1)).append(Escape.eP(v)).append(" color ")
1149         .append(color);
1150   }
1151 
1152   /**
1153    * Set pt01 to pt00, possibly adding offset into unit cell
1154    *
1155    * @param uc
1156    * @param pt00
1157    * @param pt01
1158    * @param offset
1159    */
setFractional(SymmetryInterface uc, T3 pt00, P3 pt01, P3 offset)1160   private static void setFractional(SymmetryInterface uc, T3 pt00, P3 pt01,
1161                                     P3 offset) {
1162     pt01.setT(pt00);
1163     if (offset != null)
1164       uc.toUnitCell(pt01, offset);
1165     uc.toFractional(pt01, false);
1166   }
1167 
rotTransCart(SymmetryOperation op, SymmetryInterface uc, P3 pt00, V3 vtrans)1168   private static P3 rotTransCart(SymmetryOperation op, SymmetryInterface uc,
1169                                  P3 pt00, V3 vtrans) {
1170     P3 p0 = P3.newP(pt00);
1171     uc.toFractional(p0, false);
1172     op.rotTrans(p0);
1173     p0.add(vtrans);
1174     uc.toCartesian(p0, false);
1175     return p0;
1176   }
1177 
strCoord(SymmetryOperation op, T3 p, boolean isBio)1178   private static String strCoord(SymmetryOperation op, T3 p, boolean isBio) {
1179     approx0(p);
1180     return (isBio ? p.x + " " + p.y + " " + p.z : op.fcoord2(p));
1181   }
1182 
approx0(T3 pt)1183   private static T3 approx0(T3 pt) {
1184     if (pt != null) {
1185       if (Math.abs(pt.x) < 0.0001f)
1186         pt.x = 0;
1187       if (Math.abs(pt.y) < 0.0001f)
1188         pt.y = 0;
1189       if (Math.abs(pt.z) < 0.0001f)
1190         pt.z = 0;
1191     }
1192     return pt;
1193   }
1194 
approx(T3 pt)1195   private static T3 approx(T3 pt) {
1196     if (pt != null) {
1197       pt.x = SymmetryOperation.approxF(pt.x);
1198       pt.y = SymmetryOperation.approxF(pt.y);
1199       pt.z = SymmetryOperation.approxF(pt.z);
1200     }
1201     return pt;
1202   }
1203 
1204   /**
1205    * multipurpose function handling a variety of tasks, including:
1206    *
1207    * processing of "lattice", "list", "atom", "point", and some "draw" output
1208    * types
1209    *
1210    * finding the operator in the given space group
1211    *
1212    * creating a temporary space group for an xyz operator
1213    *
1214    *
1215    * @param sym
1216    * @param iModel
1217    * @param iatom
1218    * @param uc
1219    * @param xyz
1220    * @param op
1221    * @param translation [i j k] to be added to operator
1222    * @param pt
1223    * @param pt2 second point or offset
1224    * @param id
1225    * @param type
1226    * @param scaleFactor
1227    * @param nth
1228    * @param options 0 or T.offset
1229    * @return a string or an Object[] containing information
1230    */
getSymmetryInfo(Symmetry sym, int iModel, int iatom, Symmetry uc, String xyz, int op, P3 translation, P3 pt, P3 pt2, String id, int type, float scaleFactor, int nth, int options)1231   private Object getSymmetryInfo(Symmetry sym, int iModel, int iatom,
1232                                  Symmetry uc, String xyz, int op, P3 translation,
1233                                  P3 pt, P3 pt2, String id,
1234                                  int type, float scaleFactor, int nth, int options) {
1235     int returnType = 0;
1236     Object nullRet = nullReturn(type);
1237     switch (type) {
1238     case T.lattice:
1239       return uc.getLatticeType();
1240     case T.array:
1241       returnType = getType(id);
1242       switch (returnType) {
1243       case T.atoms:
1244       case T.draw:
1245       case T.full:
1246       case T.list:
1247       case T.point:
1248         type = returnType;
1249         break;
1250       default:
1251         returnType = getKeyType(id);
1252         break;
1253       }
1254       break;
1255     }
1256     int iop = op;
1257     P3 offset = (options == T.offset  && (type == T.atoms || type == T.point)? pt2 : null);
1258     if (offset != null)
1259       pt2 = null;
1260     Object[] info = null;
1261     String xyzOriginal = null;
1262     if (pt2 == null) {
1263       if (xyz == null) {
1264         SymmetryOperation[] ops = (SymmetryOperation[]) uc
1265             .getSymmetryOperations();
1266         if (ops == null || op == 0 || Math.abs(op) > ops.length)
1267           return nullRet;
1268         iop = Math.abs(op) - 1;
1269         xyz = (translation == null ? ops[iop].xyz : ops[iop].getxyzTrans(translation));
1270         xyzOriginal = ops[iop].xyzOriginal;
1271       } else {
1272         iop = op = 0;
1273       }
1274       SymmetryInterface symTemp = new Symmetry();
1275       symTemp.setSpaceGroup(false);
1276       boolean isBio = uc.isBio();
1277       int i = (isBio ? symTemp.addBioMoleculeOperation(
1278           uc.spaceGroup.finalOperations[iop], op < 0) : symTemp
1279           .addSpaceGroupOperation((op < 0 ? "!" : "=") + xyz, Math.abs(op)));
1280 
1281       if (i < 0)
1282         return nullRet;
1283       SymmetryOperation opTemp = (SymmetryOperation) symTemp
1284           .getSpaceGroupOperation(i);
1285       if (xyzOriginal != null)
1286         opTemp.xyzOriginal = xyzOriginal;
1287       opTemp.number = op;
1288       if (!isBio)
1289         opTemp.getCentering();
1290       if (pt == null && iatom >= 0)
1291         pt = modelSet.at[iatom];
1292       if (type == T.point || type == T.atoms) {
1293         if (isBio)
1294           return nullRet;
1295         symTemp.setUnitCell(uc);
1296         pt = P3.newP(pt);
1297         uc.toFractional(pt, false);
1298         if (Float.isNaN(pt.x))
1299           return nullRet;
1300         P3 sympt = new P3();
1301         symTemp.newSpaceGroupPoint(i, pt, sympt, 0, 0, 0, null);
1302 
1303         if (options == T.offset) {
1304           uc.unitize(sympt);//uc.unitize01(sympt);
1305           if (options == T.offset)
1306             sympt.add(offset);
1307         }
1308         symTemp.toCartesian(sympt, false);
1309         return (type == T.atoms ? getAtom(uc, iModel, iatom, sympt) : sympt);
1310       }
1311       info = createInfoArray(opTemp, uc, pt, null, (id == null ? "sym" : id),
1312           scaleFactor, options, (translation != null));
1313       if (type == T.array && id != null) {
1314         returnType = getKeyType(id);
1315       }
1316     } else {
1317       // pt1, pt2
1318       String stype = "info";
1319       boolean asString = false;
1320       switch (type) {
1321       case T.array: // new Jmol 14.29.45
1322         returnType = getKeyType(id);
1323         id = stype = null;
1324         asString = false;
1325         if (nth == 0)
1326           nth = -1;
1327         break;
1328       case T.list:
1329         id = stype = null;
1330         asString = true;
1331         if (nth == 0)
1332           nth = -1;
1333         break;
1334       case T.draw:
1335         if (id == null)
1336           id = "sym";
1337         stype = "all";
1338         asString = true;
1339         break;
1340       case T.atoms:
1341         id = stype = null;
1342         //$FALL-THROUGH$
1343       default:
1344         if (nth == 0)
1345           nth = 1;
1346       }
1347       Object ret1 = getSymopInfoForPoints(sym, iModel, op, translation, pt, pt2, id,
1348           stype, scaleFactor, nth, asString, options);
1349       if (asString)
1350         return ret1;
1351       if (ret1 instanceof String)
1352         return nullRet; // two atoms are not connected, no such oper
1353       info = (Object[]) ret1;
1354       if (type == T.atoms) {
1355         if (!(pt instanceof Atom) && !(pt2 instanceof Atom))
1356           iatom = -1;
1357         return (info == null ? nullRet : getAtom(uc, iModel, iatom,
1358             (T3) info[7]));
1359       }
1360     }
1361     if (info == null)
1362       return nullRet;
1363     boolean isList = (info.length > 0 && info[0] instanceof Object[]);
1364     if (nth < 0 && op <= 0 && (type == T.array || isList)) {
1365      if (type == T.array && info.length > 0 && !(info[0] instanceof Object[]))
1366        info = new Object[] { info };
1367       Lst<Object> lst = new Lst<Object>();
1368       for (int i = 0; i < info.length; i++)
1369         lst.addLast(getInfo((Object[])info[i], returnType < 0 ? returnType : type));
1370       return lst;
1371     } else if (returnType < 0 && (nth >= 0 || op > 0)) {
1372       type = returnType;
1373     }
1374     if (nth > 0 && isList)
1375       info = (Object[]) ((Object[]) info)[0];
1376     return getInfo(info, type);
1377   }
1378 
getAtom(Symmetry uc, int iModel, int iAtom, T3 sympt)1379   private BS getAtom(Symmetry uc, int iModel, int iAtom, T3 sympt) {
1380     BS bsElement = null;
1381     if (iAtom >= 0)
1382       modelSet.getAtomBitsMDa(T.elemno,
1383           Integer.valueOf(modelSet.at[iAtom].getElementNumber()),
1384           bsElement = new BS());
1385     BS bsResult = new BS();
1386     modelSet.getAtomsWithin(0.02f, sympt, bsResult, iModel);
1387     if (bsElement != null)
1388       bsResult.and(bsElement);
1389     if (bsResult.isEmpty()) {
1390       sympt = P3.newP(sympt);
1391       uc.toUnitCell(sympt, null);
1392       uc.toCartesian(sympt, false);
1393       modelSet.getAtomsWithin(0.02f, sympt, bsResult, iModel);
1394       if (bsElement != null)
1395         bsResult.and(bsElement);
1396     }
1397     return bsResult;
1398   }
1399 
1400   /**
1401    * get information about a symmetry operation relating two specific points or atoms
1402    *
1403    * @param sym
1404    * @param modelIndex
1405    * @param symOp
1406    * @param translation TODO
1407    * @param pt1
1408    * @param pt2
1409    * @param drawID
1410    * @param stype
1411    * @param scaleFactor
1412    * @param nth
1413    * @param asString
1414    * @param options 0 or T.offset
1415    * @return Object[] or String or Object[Object[]] (nth = 0, "array")
1416    *
1417    */
getSymopInfoForPoints(Symmetry sym, int modelIndex, int symOp, P3 translation, P3 pt1, P3 pt2, String drawID, String stype, float scaleFactor, int nth, boolean asString, int options)1418   private Object getSymopInfoForPoints(Symmetry sym, int modelIndex, int symOp,
1419                                        P3 translation, P3 pt1, P3 pt2,
1420                                        String drawID, String stype, float scaleFactor,
1421                                        int nth, boolean asString, int options) {
1422     Object ret = (asString ? "" : null);
1423     Map<String, Object> sginfo = getSpaceGroupInfo(sym, modelIndex, null,
1424         symOp, pt1, pt2, drawID, scaleFactor, nth, false, true, options, null);
1425     if (sginfo == null)
1426       return ret;
1427     Object[][] infolist = (Object[][]) sginfo.get("operations");
1428     // at this point, if we have two points, we have a full list of operations, but
1429     // some are null.
1430     if (infolist == null)
1431       return ret;
1432     SB sb = (asString ? new SB() : null);
1433     symOp--;
1434     boolean isAll = (!asString && symOp < 0);
1435     String strOperations = (String) sginfo.get("symmetryInfo");
1436     boolean labelOnly = "label".equals(stype);
1437     int n = 0;
1438     for (int i = 0; i < infolist.length; i++) {
1439       if (infolist[i] == null || symOp >= 0 && symOp != i)
1440         continue;
1441       if (!asString) {
1442         if (!isAll)
1443           return infolist[i];
1444         infolist[n++] = infolist[i];
1445         continue;
1446       }
1447       if (drawID != null)
1448         return ((String) infolist[i][3]) + "\nprint " + PT.esc(strOperations);
1449       if (sb.length() > 0)
1450         sb.appendC('\n');
1451       if (!labelOnly) {
1452         if (symOp < 0)
1453           sb.appendI(i + 1).appendC('\t');
1454         sb.append((String) infolist[i][0]).appendC('\t'); //xyz
1455       }
1456       sb.append((String) infolist[i][2]); //desc
1457     }
1458     if (!asString) {
1459       Object[] a = new Object[n];
1460       for (int i = 0; i < n; i++)
1461         a[i] = infolist[i];
1462       return a;
1463     }
1464     if (sb.length() == 0)
1465       return (drawID != null ? "draw " + drawID + "* delete" : ret);
1466     return sb.toString();
1467   }
1468 }
1469