Source code

001/*
002 *                    BioJava development code
003 *
004 * This code may be freely distributed and modified under the
005 * terms of the GNU Lesser General Public Licence.  This should
006 * be distributed with the code.  If you do not have a copy,
007 * see:
008 *
009 *      http://www.gnu.org/copyleft/lesser.html
010 *
011 * Copyright for this code is held jointly by the individual
012 * authors.  These should be listed in @author doc comments.
013 *
014 * For more information on the BioJava project and its aims,
015 * or to join the biojava-l mailing list, visit the home page
016 * at:
017 *
018 *      http://www.biojava.org/
019 *
020 */
021/**
022 *
023 */
024package org.biojava.nbio.structure.symmetry.jmolScript;
025
026import org.biojava.nbio.structure.symmetry.axis.RotationAxisAligner;
027import org.biojava.nbio.structure.symmetry.core.Rotation;
028import org.biojava.nbio.structure.symmetry.core.RotationGroup;
029import org.biojava.nbio.structure.symmetry.core.QuatSymmetrySubunits;
030import org.biojava.nbio.structure.symmetry.geometry.Polyhedron;
031import org.jcolorbrewer.ColorBrewer;
032
033import javax.vecmath.*;
034
035import java.awt.*;
036import java.util.*;
037import java.util.List;
038
039/**
040 * @author Peter
041 *
042 */
043public abstract class JmolSymmetryScriptGeneratorPointGroup extends JmolSymmetryScriptGenerator {
044        private static String N_FOLD_AXIS_COLOR = "red";
045        private static String TWO_FOLD_AXIS_COLOR = "deepskyblue";
046        private static String THREE_FOLD_AXIS_COLOR = "lime";
047        private static double AXIS_SCALE_FACTOR = 1.2;
048
049        private RotationAxisAligner rotationAxisAligner = null;
050        private RotationGroup rotationGroup = null;
051        private Polyhedron polyhedron = null;
052        private String name = "";
053        private String defaultColoring = "";
054        private boolean onTheFly = true;
055
056        public JmolSymmetryScriptGeneratorPointGroup(RotationAxisAligner rotationAxisAligner, String name) {
057                this.rotationAxisAligner = rotationAxisAligner;
058                this.rotationGroup = this.rotationAxisAligner.getRotationGroup();
059                this.name = name;
060        }
061
062        /* (non-Javadoc)
063         * @see org.biojava.nbio.structure.quaternary.jmolScript.JMolSymmetryScriptInterface#getZoom()
064         */
065        @Override
066        abstract public int getZoom();
067
068        @Override
069        public void setOnTheFly(boolean onTheFly) {
070                this.onTheFly = onTheFly;
071        }
072
073        /* (non-Javadoc)
074         * @see org.biojava.nbio.structure.quaternary.jmolScript.JMolSymmetryScriptInterface#getDefaultOrientation()
075         */
076        @Override
077        public String getDefaultOrientation() {
078                StringBuilder s = new StringBuilder();
079                s.append(setCentroid());
080
081                // calculate  orientation
082                Quat4d q = new Quat4d();
083                q.set(polyhedron.getViewMatrix(0));
084                q.normalize();
085                Quat4d r = new Quat4d();
086                r.set(rotationAxisAligner.getRotationMatrix());
087                r.normalize();
088                q.mul(r);
089                q.normalize();
090
091                // set orientation
092                s.append("moveto 0 quaternion{");
093                s.append(jMolFloat(q.x));
094                s.append(",");
095                s.append(jMolFloat(q.y));
096                s.append(",");
097                s.append(jMolFloat(q.z));
098                s.append(",");
099                s.append(jMolFloat(q.w));
100                s.append("};");
101                return s.toString();
102        }
103
104        /* (non-Javadoc)
105         * @see org.biojava.nbio.structure.quaternary.jmolScript.JMolSymmetryScriptInterface#getOrientationCount()
106         */
107        @Override
108        public int getOrientationCount() {
109                return polyhedron.getViewCount();
110        }
111
112        /* (non-Javadoc)
113         * @see org.biojava.nbio.structure.quaternary.jmolScript.JMolSymmetryScriptInterface#getOrientation(int)
114         */
115        @Override
116        public String getOrientation(int index) {
117                StringBuilder s = new StringBuilder();
118                s.append(setCentroid());
119
120                // calculate  orientation
121                Quat4d q = new Quat4d();
122                q.set(polyhedron.getViewMatrix(index));
123                q.normalize();
124                Quat4d r = new Quat4d();
125                r.set(rotationAxisAligner.getRotationMatrix());
126                r.normalize();
127                q.mul(r);
128                q.normalize();
129
130                // set orientation
131                s.append("moveto 4 quaternion{");
132                s.append(jMolFloat(q.x));
133                s.append(",");
134                s.append(jMolFloat(q.y));
135                s.append(",");
136                s.append(jMolFloat(q.z));
137                s.append(",");
138                s.append(jMolFloat(q.w));
139                s.append("}");
140                s.append(";");
141                return s.toString();
142        }
143
144        /* (non-Javadoc)
145         * @see org.biojava.nbio.structure.quaternary.jmolScript.JMolSymmetryScriptInterface#getOrientationWithZoom(int)
146         */
147        @Override
148        public String getOrientationWithZoom(int index) {
149                StringBuilder s = new StringBuilder();
150                s.append(getOrientation(index));
151                s.insert(s.length()-1, getZoom());
152                return s.toString();
153
154        }
155        /* (non-Javadoc)
156         * @see org.biojava.nbio.structure.quaternary.jmolScript.JMolSymmetryScriptInterface#getOrientationName(int)
157         */
158        @Override
159        public String getOrientationName(int index) {
160                return polyhedron.getViewName(index);
161        }
162
163        /* (non-Javadoc)
164         * @see org.biojava.nbio.structure.quaternary.jmolScript.JMolSymmetryScriptInterface#getTransformation()
165         */
166        @Override
167        public Matrix4d getTransformation() {
168                return rotationAxisAligner.getTransformation();
169        }
170
171        @Override
172        public void setDefaultColoring(String colorScript) {
173                this.defaultColoring = colorScript;
174        }
175
176        /* (non-Javadoc)
177         * @see org.biojava.nbio.structure.quaternary.jmolScript.JMolSymmetryScriptInterface#drawPolyhedron()
178         */
179        @Override
180        public String drawPolyhedron() {
181                StringBuilder s = new StringBuilder();
182
183                Point3d[] vertices = getPolyhedronVertices();
184
185                int index = 0;
186                double width = getMaxExtension()*0.015;
187
188                for (int[] lineLoop: polyhedron.getLineLoops()) {
189                        s.append("draw polyhedron");
190                        s.append(name);
191                        s.append(index++);
192                        s.append(" line");
193                        for (int i: lineLoop) {
194                                s.append(getJmolPoint(vertices[i]));
195                        }
196                        s.append("width ");
197                        s.append(fDot2(width));
198                        s.append(" color");
199                        Color4f c = getPolyhedronColor();
200                        s.append(getJmolColor(c));
201                        s.append(" off;");
202                }
203
204                return s.toString();
205        }
206
207        /* (non-Javadoc)
208         * @see org.biojava.nbio.structure.quaternary.jmolScript.JMolSymmetryScriptInterface#hidePolyhedron()
209         */
210        @Override
211        public String hidePolyhedron() {
212                return "draw polyhedron* off;";
213        }
214
215        /* (non-Javadoc)
216         * @see org.biojava.nbio.structure.quaternary.jmolScript.JMolSymmetryScriptInterface#showPolyhedron()
217         */
218        @Override
219        public String showPolyhedron() {
220                return "draw polyhedron* on;";
221        }
222
223        /* (non-Javadoc)
224         * @see org.biojava.nbio.structure.quaternary.jmolScript.JMolSymmetryScriptInterface#drawAxes()
225         */
226        @Override
227        public String drawAxes() {
228                if ("C1".equals(rotationGroup.getPointGroup())) {
229                        return drawInertiaAxes();
230                } else {
231                        return drawSymmetryAxes();
232                }
233        }
234
235        /* (non-Javadoc)
236         * @see org.biojava.nbio.structure.quaternary.jmolScript.JMolSymmetryScriptInterface#hideAxes()
237         */
238        @Override
239        public String hideAxes() {
240                return "draw axes* off;";
241        }
242
243        /* (non-Javadoc)
244         * @see org.biojava.nbio.structure.quaternary.jmolScript.JMolSymmetryScriptInterface#showAxes()
245         */
246        @Override
247        public String showAxes() {
248                return "draw axes* on;";
249        }
250
251        /* (non-Javadoc)
252         * @see org.biojava.nbio.structure.quaternary.jmolScript.JMolSymmetryScriptInterface#playOrientations()
253         */
254        @Override
255        public String playOrientations() {
256                StringBuilder s = new StringBuilder();
257
258                // draw point group
259
260                if ( "C1".equals(rotationGroup.getPointGroup())) {
261                        s.append(drawFooter("Asymmetric", "white"));
262                } else {
263                        s.append(drawFooter("Point group " + rotationGroup.getPointGroup(), "white"));
264                }
265
266                // draw polygon
267                s.append(drawPolyhedron()); // draw invisibly
268                s.append(showPolyhedron());
269
270                // draw axes
271                s.append(drawAxes());
272                s.append(showAxes());
273
274                // loop over all orientations with 4 sec. delay
275                for (int i = 0; i < getOrientationCount(); i++) {
276                        s.append(deleteHeader());
277                        s.append(getOrientationWithZoom(i));
278                        s.append(drawHeader(polyhedron.getViewName(i), "white"));
279                        s.append("delay 4;");
280                }
281
282                // go back to first orientation
283                s.append(deleteHeader());
284                s.append(getOrientationWithZoom(0));
285                s.append(drawHeader(polyhedron.getViewName(0), "white"));
286
287                return s.toString();
288        }
289
290        /* (non-Javadoc)
291         * @see org.biojava.nbio.structure.quaternary.jmolScript.JMolSymmetryScriptInterface#colorBySubunit()
292         */
293        @Override
294        public String colorBySubunit() {
295                QuatSymmetrySubunits subunits = rotationAxisAligner.getSubunits();
296                List<Integer> modelNumbers = subunits.getModelNumbers();
297                List<String> chainIds = subunits.getChainIds();
298                List<List<Integer>> orbits = rotationAxisAligner.getOrbits();
299                int fold = rotationGroup.getRotation(0).getFold();
300
301                Color[] col = null;
302                Color4f[] colors = null;
303                if (fold > 1) {
304                        col = ColorBrewer.Spectral.getColorPalette(2*fold);
305                        colors = ColorConverter.convertColor4f(col);
306                } else {
307                        col = ColorBrewer.Spectral.getColorPalette(orbits.size());
308                        colors = ColorConverter.convertColor4f(col);
309                }
310                int half = colors.length/2;
311                for (int i = 0; i < half; i++) {
312                        if (i % 2 != 0) {
313                           Color4f temp = colors[i];
314                           colors[i] = colors[half+i];
315                           colors[half+i] = temp;
316                        }
317                }
318                Map<Color4f, List<String>> colorMap = new HashMap<Color4f, List<String>>();
319
320                for (int i = 0; i < orbits.size(); i++) {
321                        for (int j = 0; j < fold; j++) {
322                                // assign alternating color sets to adjacent orbits
323                                int colorIndex = i;
324                                if (fold > 1) {
325                                        if (i % 2 == 0) {
326                                                colorIndex = j;
327                                        } else {
328                                                colorIndex = fold + j;
329                                        }
330                                }
331                                int subunit = orbits.get(i).get(j);
332                                Color4f c = colors[colorIndex];
333                                List<String> ids = colorMap.get(c);
334                                if (ids == null) {
335                                        ids = new ArrayList<>();
336                                        colorMap.put(c,  ids);
337                                }
338                                String id = getChainSpecification(modelNumbers, chainIds, subunit);
339                                ids.add(id);
340                        }
341                }
342                return defaultColoring + getJmolColorScript(colorMap) + getJmolLigandScript();
343        }
344
345
346        /* (non-Javadoc)
347         * @see org.biojava.nbio.structure.quaternary.jmolScript.JMolSymmetryScriptInterface#colorBySequenceCluster()
348         */
349        @Override
350        public String colorBySequenceCluster() {
351                QuatSymmetrySubunits subunits = rotationAxisAligner.getSubunits();
352                int n = subunits.getSubunitCount();
353                List<Integer> modelNumbers = subunits.getModelNumbers();
354                List<String> chainIds = subunits.getChainIds();
355                List<Integer> seqClusterIds = subunits.getClusterIds();
356                int clusters = Collections.max(seqClusterIds) + 1;
357                Color[] col = ColorBrewer.BrBG.getColorPalette(clusters);
358                Color4f[] colors = ColorConverter.convertColor4f(col);
359                Map<Color4f, List<String>> colorMap = new HashMap<Color4f, List<String>>();
360
361                for (int i = 0; i < n; i++) {
362                        Color4f c = colors[seqClusterIds.get(i)];
363                        List<String> ids = colorMap.get(c);
364                        if (ids == null) {
365                                ids = new ArrayList<>();
366                                colorMap.put(c,  ids);
367                        }
368                        String id = getChainSpecification(modelNumbers, chainIds, i);
369                        ids.add(id);
370
371                }
372
373                return defaultColoring + getJmolColorScript(colorMap) + getJmolLigandScript();
374        }
375
376        /* (non-Javadoc)
377         * @see org.biojava.nbio.structure.quaternary.jmolScript.JMolSymmetryScriptInterface#colorBySymmetry()
378         */
379        @Override
380        public String colorBySymmetry() {
381                // TODO needs some refactoring
382                String pointGroup = rotationGroup.getPointGroup();
383                QuatSymmetrySubunits subunits = rotationAxisAligner.getSubunits();
384                List<Integer> modelNumbers = subunits.getModelNumbers();
385                List<String> chainIds = subunits.getChainIds();
386                List<List<Integer>> orbits = rotationAxisAligner.getOrbits();
387
388                int n = subunits.getSubunitCount();
389                int fold = rotationGroup.getRotation(0).getFold();
390
391                Map<Color4f, List<String>> colorMap = new HashMap<Color4f, List<String>>();
392
393                // Simple Cn symmetry
394                if (pointGroup.startsWith("C") && n == fold) {
395                        colorMap = getCnColorMap();
396                        // complex cases
397                } else if ((pointGroup.startsWith("D") && orbits.size() > 2) ||
398                                "T".equals(pointGroup)|| "O".equals(pointGroup) || "I".equals(pointGroup)) {
399                        int nColor = 0;
400                        if (orbits.size() % 2 == 0) {
401                                nColor = orbits.size()/2;
402                        } else {
403                                nColor = (orbits.size() + 1)/2;
404                        }
405                        Color4f[] colors = getSymmetryColors(nColor);
406
407                        for (int i = 0; i < orbits.size(); i++) {
408                                int colorIndex = i;
409
410                                // reverse colors once the center of the structure has been reached
411                                if (i >= nColor) {
412                                        colorIndex = orbits.size() - 1 - i;
413                                }
414                                Color4f c = colors[colorIndex];
415                                List<String> ids = colorMap.get(c);
416                                if (ids == null) {
417                                        ids = new ArrayList<>();
418                                        colorMap.put(c,  ids);
419                                }
420                                for (int subunit: orbits.get(i)) {
421                                        String id = getChainSpecification(modelNumbers, chainIds, subunit);
422                                        ids.add(id);
423                                }
424                        }
425
426                        // Simple Dn symmetry
427                } else {
428                        Color4f[] colors = getSymmetryColors(orbits.size());
429
430                        for (int i = 0; i < orbits.size(); i++) {
431                                Color4f c = new Color4f(colors[i]);
432                                List<String> ids = colorMap.get(c);
433                                if (ids == null) {
434                                        ids = new ArrayList<>();
435                                        colorMap.put(c,  ids);
436                                }
437                                List<Integer> orbit = orbits.get(i);
438                                for (int j = 0; j < orbit.size(); j++) {
439                                        String id = getChainSpecification(modelNumbers, chainIds, orbit.get(j));
440                                        ids.add(id);
441                                }
442                        }
443                }
444
445                return defaultColoring + getJmolColorScript(colorMap) + getJmolLigandScript();
446        }
447
448        // --- protected methods ---
449        /**
450         * Returns the maximum extension (length) of structure
451         * @return
452         */
453        protected double getMaxExtension() {
454                Vector3d dimension = rotationAxisAligner.getDimension();
455                double maxExtension = Math.max(dimension.x, dimension.y);
456                maxExtension = Math.max(maxExtension, dimension.z);
457                return maxExtension;
458        }
459
460        /**
461         * Returns the mean extension (length) of structure
462         * @return
463         */
464        protected double getMeanExtension() {
465                Vector3d dimension = rotationAxisAligner.getDimension();
466                return (dimension.x+dimension.y+dimension.z)/3;
467        }
468
469        /**
470         * @return the axisTransformation
471         */
472        protected RotationAxisAligner getAxisTransformation() {
473                return rotationAxisAligner;
474        }
475
476        /**
477         * @param axisTransformation the axisTransformation to set
478         */
479        protected void setAxisTransformation(RotationAxisAligner axisTransformation) {
480                this.rotationAxisAligner = axisTransformation;
481        }
482
483        /**
484         * @return the rotationGroup
485         */
486        protected RotationGroup getRotationGroup() {
487                return rotationGroup;
488        }
489
490        /**
491         * @param rotationGroup the rotationGroup to set
492         */
493        protected void setRotationGroup(RotationGroup rotationGroup) {
494                this.rotationGroup = rotationGroup;
495        }
496
497        /**
498         * @return the polyhedron
499         */
500        protected Polyhedron getPolyhedron() {
501                return polyhedron;
502        }
503
504        /**
505         * @param polyhedron the polyhedron to set
506         */
507        protected void setPolyhedron(Polyhedron polyhedron) {
508                this.polyhedron = polyhedron;
509        }
510
511//  --- private methods ---
512
513        private String getChainSpecification(List<Integer> modelNumbers, List<String> chainIds, int subunit) {
514                if (onTheFly) {
515                        if (Collections.max(modelNumbers) > 1) {
516                                return chainIds.get(subunit) + "&symop=" + (modelNumbers.get(subunit)+1);
517                        } else {
518                                // if there is only a single symop, Jmol does not accept the symop syntax
519                                return chainIds.get(subunit);
520                        }
521                } else {
522                        return chainIds.get(subunit) + "/" + (modelNumbers.get(subunit)+1);
523                }
524        }
525
526        private Map<Color4f, List<String>> getCnColorMap() {
527                QuatSymmetrySubunits subunits = rotationAxisAligner.getSubunits();
528                List<Integer> modelNumbers = subunits.getModelNumbers();
529                List<String> chainIds = subunits.getChainIds();
530                List<List<Integer>> orbits = rotationAxisAligner.getOrbits();
531
532                int fold = rotationGroup.getRotation(0).getFold();
533
534                Map<Color4f, List<String>> colorMap = new HashMap<Color4f, List<String>>();
535                Color4f[] colors = getSymmetryColors(fold);
536
537                for (List<Integer> orbit: orbits) {
538                        for (int i = 0; i < fold; i++) {
539                                int subunit = orbit.get(i);
540                                String id = null;
541                                id = getChainSpecification(modelNumbers, chainIds, subunit);
542                                Color4f c = colors[i];
543                                List<String> ids = colorMap.get(c);
544                                if (ids == null) {
545                                        ids = new ArrayList<>();
546                                        colorMap.put(c, ids);
547                                }
548                                ids.add(id);
549                        }
550                }
551
552                return colorMap;
553        }
554
555        private Point3d[] getPolyhedronVertices() {
556                Point3d[] vertices = polyhedron.getVertices();
557                Matrix4d reverseTransformation = rotationAxisAligner.getGeometicCenterTransformation();
558                for (int i = 0; i < vertices.length; i++) {
559                        reverseTransformation.transform(vertices[i]);
560                }
561                return vertices;
562        }
563
564        /**
565         * Return a color that is complementary to the symmetry color
566         * @return
567         */
568        private Color4f getPolyhedronColor() {
569                Color4f[] colors = getSymmetryColors(5);
570                Color4f strongestColor = colors[4];
571                Color4f complement = new Color4f(Color.WHITE);
572                complement.sub(strongestColor);
573                return complement;
574        }
575
576        /**
577         * Returns a unique color palette based on point group
578         * @param nColors
579         * @return
580         */
581        private Color4f[] getSymmetryColors(int nColors) {
582                String pointGroup = rotationGroup.getPointGroup();
583                Color[] col = null;
584                Color4f[] colors = null;
585                if ("C1".equals(pointGroup)) {
586                        col = ColorBrewer.Greys.getColorPalette(nColors);
587                        colors = ColorConverter.convertColor4f(col);
588                } else if (pointGroup.startsWith("C")) {
589                        col = ColorBrewer.YlGnBu.getColorPalette(nColors);
590                        colors = ColorConverter.convertColor4f(col);
591                } else if (pointGroup.startsWith("D")) {
592                        col = ColorBrewer.YlOrRd.getColorPalette(nColors);
593                        colors = ColorConverter.convertColor4f(col);
594                } else if ("T".equals(pointGroup)) {
595                        col = ColorBrewer.Greens.getColorPalette(nColors);
596                        colors = ColorConverter.convertColor4f(col);
597                } else if ("O".equals(pointGroup)) {
598                        col = ColorBrewer.Blues.getColorPalette(nColors);
599                        colors = ColorConverter.convertColor4f(col);
600                } else if ("I".equals(pointGroup)) {
601                        col = ColorBrewer.BuPu.getColorPalette(nColors);
602                        colors = ColorConverter.convertColor4f(col);
603                } else {
604                        col = ColorBrewer.Greys.getColorPalette(nColors);
605                        colors = ColorConverter.convertColor4f(col);
606                }
607                System.arraycopy(colors, 0, colors, 0, nColors);
608                return colors;
609        }
610
611        private String drawInertiaAxes() {
612                StringBuilder s = new StringBuilder();
613                Point3d centroid = rotationAxisAligner.getGeometricCenter();
614                Vector3d[] axes = rotationAxisAligner.getPrincipalAxesOfInertia();
615
616                for (int i = 0; i < axes.length; i++) {
617                        s.append("draw axesInertia");
618                        s.append(name);
619                        s.append(i);
620                        s.append(" ");
621                        s.append("line");
622                        Point3d v1 = new Point3d(axes[i]);
623                        if (i == 0) {
624                                v1.scale(AXIS_SCALE_FACTOR*rotationAxisAligner.getDimension().y);
625                        } else if (i == 1) {
626                                v1.scale(AXIS_SCALE_FACTOR*rotationAxisAligner.getDimension().x);
627                        } else if (i == 2) {
628                                v1.scale(AXIS_SCALE_FACTOR*rotationAxisAligner.getDimension().z);
629                        }
630                        Point3d v2 = new Point3d(v1);
631                        v2.negate();
632                        v1.add(centroid);
633                        v2.add(centroid);
634                        s.append(getJmolPoint(v1));
635                        s.append(getJmolPoint(v2));
636                        s.append("width 0.5 ");
637                        s.append(" color white");
638                        s.append(" off;");
639                }
640                return s.toString();
641        };
642
643        private String drawSymmetryAxes() {
644                StringBuilder s = new StringBuilder();
645
646                int n = rotationGroup.getOrder();
647                if (n == 0) {
648                        return s.toString();
649                }
650
651                float diameter = 0.5f;
652                double radius = 0;
653                String color = "";
654
655                List<Rotation> axes = getUniqueAxes();
656
657                int i = 0;
658                for (Rotation r: axes) {
659                        if (rotationGroup.getPointGroup().startsWith("C") || (rotationGroup.getPointGroup().startsWith("D") && r.getDirection() == 0)) {
660                                radius =  rotationAxisAligner.getDimension().z; // principal axis uses z-dimension
661                                color = N_FOLD_AXIS_COLOR;
662                        } else {
663                                radius = polyhedron.getCirumscribedRadius();
664
665                                if (r.getFold() == 2) {
666                                        color = TWO_FOLD_AXIS_COLOR;
667                                } else if (r.getFold() == 3) {
668                                        color = THREE_FOLD_AXIS_COLOR;
669                                } else {
670                                        color = N_FOLD_AXIS_COLOR;
671                                }
672                        }
673
674
675                        Point3d center = rotationAxisAligner.getGeometricCenter();
676                        AxisAngle4d axisAngle = r.getAxisAngle();
677                        Vector3d axis = new Vector3d(axisAngle.x, axisAngle.y, axisAngle.z);
678                        Vector3d refAxis = rotationAxisAligner.getRotationReferenceAxis();
679
680                        s.append(getSymmetryAxis(i, i+axes.size(), rotationGroup.getPointGroup(), r.getFold(), refAxis, radius, diameter, color, center, axis));
681                        i++;
682                }
683
684                return s.toString();
685        }
686
687
688        private Vector3d getAlignmentVector(Point3d point, Vector3d axis) {
689                // for system with a single Cn axis
690                if (rotationGroup.getPointGroup().startsWith("C") || rotationGroup.getPointGroup().startsWith("D")) {
691                        // if axis is orthogonal to principal axis, use principal axis as reference axis
692                        if (axis.dot(rotationAxisAligner.getPrincipalRotationAxis()) < 0.1) {
693                                return rotationAxisAligner.getPrincipalRotationAxis();
694                        } else {
695                                return rotationAxisAligner.getRotationReferenceAxis();
696                        }
697                }
698
699                // for T, O, and I point groups find reference axis with respect to
700                // nearest polyhedron vertex
701                Vector3d ref = new Vector3d();
702                double dSqThreshold = 25;
703                double dSqMin = Double.MAX_VALUE;
704                Point3d refPoint = null;
705                // find nearest point on polyhedron as reference point,
706                // but do not choose a point on the same axis (therefore, we
707                // apply a distance threshold squared 5A*5A = 25A^2
708                for (Point3d v: getPolyhedronVertices()) {
709                        double dSq = point.distanceSquared(v);
710                        if (dSq > dSqThreshold) {
711                                if (dSq < dSqMin) {
712                                        dSqMin = dSq;
713                                        refPoint = v;
714                                }
715                        }
716                }
717
718
719                ref.sub(point, refPoint);
720
721                // this ref vector is usually not orthogonal to the
722                // axis. The following double-cross product makes it
723                // orthogonal.
724                ref.cross(axis, ref);
725                ref.cross(axis, ref); // note, done twice on purpose
726                ref.normalize();
727                return ref;
728        }
729
730        private String getSymmetryAxis(int i, int j, String pointGroup, int n, Vector3d referenceAxis, double radius, float diameter, String color, Point3d center, Vector3d axis) {
731                boolean drawPolygon = true;
732
733                Point3d p1 = new Point3d(axis);
734                p1.scaleAdd(-AXIS_SCALE_FACTOR * radius, center);
735
736                Point3d p2 = new Point3d(axis);
737                p2.scaleAdd(AXIS_SCALE_FACTOR * radius, center);
738
739                StringBuilder s = new StringBuilder();
740                s.append("draw");
741                s.append(" axesSymmetry");
742                s.append(name);
743                s.append(i);
744                s.append(" cylinder");
745                s.append(getJmolPoint(p1));
746                s.append(getJmolPoint(p2));
747                s.append("diameter ");
748                s.append(diameter);
749                s.append(" color ");
750                s.append(color);
751                s.append(" off;");
752
753                // calc. point to center symmetry symbols. They are offset by 0.01
754                // to avoid overlap with the polyhedron
755                p1 = new Point3d(axis);
756                p1.scaleAdd(-1.01*radius, center);
757
758                p2 = new Point3d(axis);
759                p2.scaleAdd(1.01*radius, center);
760
761                if (drawPolygon) {
762                        double polygonRadius = getMeanExtension() * 0.06;
763                        if (n == 2) {
764                                referenceAxis = getAlignmentVector(p1, axis);
765                                s.append(getC2PolygonJmol(i, p1, referenceAxis, axis, color, polygonRadius, name));
766                                referenceAxis = getAlignmentVector(p2, axis);
767                                s.append(getC2PolygonJmol(j, p2,  referenceAxis, axis, color, polygonRadius, name));
768                        } else if (n > 2) {
769                                referenceAxis = getAlignmentVector(p1, axis);
770                                s.append(getPolygonJmol(i, p1, referenceAxis, axis, n, color, polygonRadius, name));
771                                referenceAxis = getAlignmentVector(p2, axis);
772                                s.append(getPolygonJmol(j, p2, referenceAxis, axis, n, color, polygonRadius, name));
773                        }
774                }
775
776                return s.toString();
777        }
778
779        private static String getPolygonJmol(int index, Point3d center, Vector3d referenceAxis, Vector3d axis, int n, String color, double radius, String name) {
780                StringBuilder s = new StringBuilder();
781                s.append("draw axesSymbol");
782                s.append(name);
783                s.append(index);
784                s.append(" ");
785                s.append("polygon");
786                s.append(" ");
787                s.append(n+1);
788                s.append(getJmolPoint(center));
789
790                Vector3d[] vertexes = getPolygonVertices(axis, referenceAxis, center, n, radius);
791                // create vertex list
792                for (Vector3d v: vertexes) {
793                        s.append(getJmolPoint(v));
794                }
795
796                // create face list
797                s.append(n);
798                for (int i = 1; i <= n; i++) {
799                        s.append("[");
800                        s.append(0);
801                        s.append(" ");
802                        s.append(i);
803                        s.append(" ");
804                        if (i < n) {
805                                s.append(i+1);
806                        } else {
807                                s.append(1);
808                        }
809                        s.append(" ");
810                        s.append(7);
811                        s.append("]");
812                }
813
814                if (n == 2) {
815                        s.append("mesh off");
816                }
817                s.append(" color ");
818                s.append(color);
819                s.append(" off;");
820
821                return s.toString();
822        }
823
824        private static Vector3d[] getPolygonVertices(Vector3d axis, Vector3d referenceAxis, Point3d center, int n, double radius) {
825                Vector3d ref = new Vector3d(referenceAxis);
826                ref.scale(radius);
827
828                AxisAngle4d axisAngle = new AxisAngle4d(axis, 0);
829                Vector3d[] vectors = new Vector3d[n];
830                Matrix4d m = new Matrix4d();
831
832                for (int i = 0; i < n; i++) {
833                        axisAngle.angle = i * 2 * Math.PI/n;
834                        vectors[i] = new Vector3d(ref);
835                        m.set(axisAngle);
836                        // make sure matrix element m33 is 1.0. It's 0 on Linux.
837                        m.setElement(3, 3, 1.0);
838                        m.transform(vectors[i]);
839                        vectors[i].add(center);
840                }
841                return vectors;
842        }
843
844        private static String getC2PolygonJmol(int index, Point3d center, Vector3d referenceAxis, Vector3d axis, String color, double radius, String name) {
845                StringBuilder s = new StringBuilder();
846                int n = 10;
847                s.append("draw axesSymbol");
848                s.append(name);
849                s.append(index);
850                s.append(" ");
851                s.append("polygon");
852                s.append(" ");
853                s.append(n-1);
854                s.append(getJmolPoint(center));
855
856                Vector3d[] vertexes = getC2PolygonVertices(axis, referenceAxis, center, n, radius);
857                // create vertex list
858                for (Vector3d v: vertexes) {
859                        s.append(getJmolPoint(v));
860                }
861
862                // create face list
863                s.append(n-2);
864
865                for (int i = 1; i < n-1; i++) {
866                        s.append("[");
867                        s.append(0);
868                        s.append(" ");
869                        s.append(i);
870                        s.append(" ");
871                        if (i < n-2) {
872                                s.append(i+1);
873                        } else {
874                                s.append(1);
875                        }
876                        s.append(" ");
877                        s.append(7);
878                        s.append("]");
879                }
880
881                s.append("color ");
882                s.append(color);
883                s.append(" off;");
884
885                return s.toString();
886        }
887        private static Vector3d[] getC2PolygonVertices(Vector3d axis, Vector3d referenceAxis, Point3d center, int n, double radius) {
888                Vector3d ref = new Vector3d(referenceAxis);
889                ref.scale(4*radius);
890
891                AxisAngle4d axisAngle = new AxisAngle4d(axis, 0);
892                int k = n / 2;
893                // draw arc 1/6 of a full circle
894                int f = 6;
895                Vector3d[] vectors = new Vector3d[n-2];
896                Matrix4d m = new Matrix4d();
897
898                // first point of arc
899                axisAngle.angle = (k+0.5) * 2 * Math.PI/(f*k);
900                Vector3d begin = new Vector3d(ref);
901                m.set(axisAngle);
902                // make sure matrix element m33 is 1.0. It's 0 on Linux.
903                m.setElement(3, 3, 1.0);
904                m.transform(begin);
905
906                // last point of arc
907                axisAngle.angle = (2*k-1+0.5) * 2 * Math.PI/(f*k);
908                Vector3d end = new Vector3d(ref);
909                m.set(axisAngle);
910                // make sure matrix element m33 is 1.0. It's 0 on Linux.
911                m.setElement(3, 3, 1.0);
912                m.transform(end);
913
914                // center of arc
915                Vector3d arcCenter = new Vector3d();
916                arcCenter.interpolate(begin, end, 0.5);
917                arcCenter.negate();
918
919                // add translation component
920                Vector3d trans =  new Vector3d();
921                trans.sub(center, arcCenter);
922
923                // draw arc (1/6 of a full circle)
924                for (int i = 0; i < k; i++) {
925                        axisAngle.angle = (k + i + 0.5) * 2 * Math.PI/(f*k);
926                        vectors[i] = new Vector3d(ref);
927                        m.set(axisAngle);
928                        // make sure matrix element m33 is 1.0. It's 0 on Linux.
929                        m.setElement(3, 3, 1.0);
930                        m.transform(vectors[i]);
931                        vectors[i].add(arcCenter);
932                        vectors[i].add(center);
933                }
934                // in reverse order, draw reflected half of arc (1/6 of full circle)
935                // don't draw first and last element, since the are already part of the previous arc
936                for (int i = k; i < 2*k-2; i++) {
937                        axisAngle.angle = (f/2*k + i + 1.5) * 2 * Math.PI/(f*k);
938                        vectors[i] = new Vector3d(ref);
939                        m.set(axisAngle);
940                        // make sure matrix element m33 is 1.0. It's 0 on Linux.
941                        m.setElement(3, 3, 1.0);
942                        m.transform(vectors[i]);
943                        vectors[i].sub(arcCenter);
944                        vectors[i].add(center);
945                }
946                return vectors;
947        }
948
949
950        private List<Rotation> getUniqueAxes() {
951                List<Rotation> uniqueRotations = new ArrayList<>();
952
953                for (int i = 0, n = rotationGroup.getOrder(); i < n; i++) {
954                        Rotation rotationI = rotationGroup.getRotation(i);
955                        AxisAngle4d axisAngleI = rotationI.getAxisAngle();
956                        Vector3d axisI = new Vector3d(axisAngleI.x, axisAngleI.y, axisAngleI.z);
957
958                        boolean redundant = false;
959                        for (Rotation r: uniqueRotations) {
960                                AxisAngle4d axisAngleJ = r.getAxisAngle();
961                                Vector3d axisJ = new Vector3d(axisAngleJ.x, axisAngleJ.y, axisAngleJ.z);
962                                if (Math.abs(axisI.dot(axisJ)) > 0.99) {
963                                        redundant = true;
964                                        break;
965                                }
966                        }
967
968                        if (! redundant) {
969                                uniqueRotations.add(rotationI);
970                        }
971                }
972                return uniqueRotations;
973        }
974
975        private String drawHeader(String text, String color) {
976                StringBuilder s = new StringBuilder();
977                s.append("set echo top center;");
978                s.append("color echo ");
979                s.append(color);
980                s.append(";");
981                s.append("font echo 24 sanserif;");
982                s.append("echo ");
983                s.append(text);
984                s.append(";");
985                return s.toString();
986        }
987
988        private String deleteHeader() {
989                return "set echo top center;echo ;";
990        }
991
992        private String drawFooter(String text, String color) {
993                StringBuilder s = new StringBuilder();
994                s.append("set echo bottom center;");
995                s.append("color echo ");
996                s.append(color);
997                s.append(";");
998                s.append("font echo 24 sanserif;");
999                s.append("echo ").append(text);
1000                //s.append("echo Point group ");
1001                //s.append(rotationGroup.getPointGroup());
1002                s.append(";");
1003                return s.toString();
1004        }
1005
1006        private String setCentroid() {
1007                // calculate center of rotation
1008        //      Point3d centroid = axisTransformation.getGeometricCenter();
1009                Point3d centroid = rotationAxisAligner.getCentroid();
1010
1011                // set centroid
1012                StringBuilder s = new StringBuilder();
1013                s.append("center");
1014                s.append(getJmolPoint(centroid));
1015                s.append(";");
1016                return s.toString();
1017        }
1018}