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
022package org.biojava.nbio.structure.symmetry.core;
023
024import javax.vecmath.AxisAngle4d;
025import javax.vecmath.Matrix4d;
026import javax.vecmath.Vector3d;
027import java.util.ArrayList;
028import java.util.Collections;
029import java.util.Comparator;
030import java.util.Iterator;
031import java.util.List;
032
033/**
034 * @see http://en.wikipedia.org/wiki/Rotation_group_SO(3)
035 * @author Peter
036 */
037public class RotationGroup implements Iterable<Rotation> {
038        private List<Rotation> rotations = new ArrayList<Rotation>();
039        private int principalAxisIndex = 0;
040        private int higherOrderRotationAxis = 0;
041        private int twoFoldsPerpendicular = 0;
042        private int highestOrder = 0;
043        private String pointGroup = "C1";
044        private double symmetryDeviation = 0;
045        private boolean complete = true;
046        private boolean modified = true;
047
048
049        public int getOrder() {
050                return rotations.size();
051        }
052
053        public Rotation getRotation(int index) {
054                return rotations.get(index);
055        }
056
057        public void addRotation(Rotation rotation) {
058                rotations.add(rotation);
059                modified = true;
060        }
061
062        public void setC1(int n) {
063                Rotation r = new Rotation();
064                List<Integer> permutation = new ArrayList<Integer>(n);
065                for (int i = 0; i < n; i++) {
066                        permutation.add(i);
067                }
068                r.setPermutation(permutation);
069                Matrix4d m = new Matrix4d();
070                m.setIdentity();
071                r.setTransformation(m);
072                r.setAxisAngle(new AxisAngle4d());
073                r.setFold(1);
074                r.setScores(new QuatSymmetryScores());
075                rotations.add(r);
076                pointGroup = "C1";
077        }
078
079        public void removeRotation(int index) {
080                rotations.remove(index);
081                modified = true;
082        }
083
084        public void complete() {
085                if (modified) {
086                        if (rotations.size() > 0) {
087                                findHighestOrderAxis();
088                                setEAxis();
089                                calcAxesDirections();
090                                findHigherOrderAxes();
091                                findTwoFoldsPerpendicular();
092                                calcPointGroup();
093                                sortByFoldDecending();
094                        }
095                        modified = false;
096                }
097        }
098
099        public String getPointGroup() {
100                if (modified) {
101                        if (rotations.size() == 0) {
102                                return "C1";
103                        }
104                        complete();
105                }
106                return pointGroup;
107        }
108
109        /**
110         * Returns QuatSymmetryScores averaged over all rotations
111         * (except the first rotation, which is the unit operation E)
112         * @return mean scores average over rotations
113         */
114        public QuatSymmetryScores getScores() {
115                QuatSymmetryScores scores = new QuatSymmetryScores();
116
117                int n = rotations.size()-1;
118
119                if (n > 0) {
120                        double[] values = new double[n];
121
122                        // minRmsd
123                        for (int i = 1; i < rotations.size(); i++) {
124                                values[i-1] = rotations.get(i).getScores().getMinRmsd();
125                        }
126                        scores.setMinRmsd(minScores(values));
127
128                        // maxRmsd
129                        for (int i = 1; i < rotations.size(); i++) {
130                                values[i-1] = rotations.get(i).getScores().getMaxRmsd();
131                        }
132                        scores.setMaxRmsd(maxScores(values));
133
134                        // Rmsd
135                        for (int i = 1; i < rotations.size(); i++) {
136                                values[i-1] = rotations.get(i).getScores().getRmsd();
137                        }
138                        scores.setRmsd(averageScores(values));
139
140                        // minTm
141                        for (int i = 1; i < rotations.size(); i++) {
142                                values[i-1] = rotations.get(i).getScores().getMinTm();
143                        }
144                        scores.setMinTm(minScores(values));
145
146                        // maxTm
147                        for (int i = 1; i < rotations.size(); i++) {
148                                values[i-1] = rotations.get(i).getScores().getMaxTm();
149                        }
150                        scores.setMaxTm(maxScores(values));
151
152                        // Tm
153                        for (int i = 1; i < rotations.size(); i++) {
154                                values[i-1] = rotations.get(i).getScores().getTm();
155                        }
156                        scores.setTm(averageScores(values));
157
158                        // Rmsd subunit centers
159                        for (int i = 1; i < rotations.size(); i++) {
160                                values[i-1] = rotations.get(i).getScores().getRmsdCenters();
161                        }
162                        scores.setRmsdCenters(averageScores(values));
163                        // TmIntra
164                        for (int i = 1; i < rotations.size(); i++) {
165                                values[i-1] = rotations.get(i).getScores().getTmIntra();
166                        }
167                        scores.setTmIntra(averageScores(values));
168
169                        // RmsdIntra
170                        for (int i = 1; i < rotations.size(); i++) {
171                                values[i-1] = rotations.get(i).getScores().getRmsdIntra();
172                        }
173                        scores.setRmsdIntra(averageScores(values));
174
175                        // SymDeviation
176                        scores.setSymDeviation(symmetryDeviation);
177                }
178                return scores;
179        }
180
181        /**
182         * @param symmetryDeviation the symmetryDeviation to set
183         */
184        public void setSymmetryDeviation(double symmetryDeviation) {
185                this.symmetryDeviation = symmetryDeviation;
186        }
187
188        public boolean isComplete() {
189                return complete;
190        }
191
192        @Override
193        public String toString() {
194                StringBuilder sb = new StringBuilder();
195                sb.append("Rotations: " + rotations.size() + "\n");
196                for (Rotation s: rotations) {
197                        sb.append(s.toString()).append("\n");
198                }
199                return sb.toString();
200        }
201
202        private double averageScores(double[] scores) {
203                double sum = 0;
204                for (double s: scores) {
205                        sum += s;
206                }
207                return sum/scores.length;
208        }
209
210        private double minScores(double[] scores) {
211                double score = Double.MAX_VALUE;
212                for (double s: scores) {
213                        score = Math.min(score, s);
214                }
215                return score;
216        }
217
218        private double maxScores(double[] scores) {
219                double score = Double.MIN_VALUE;
220                for (double s: scores) {
221                        score = Math.max(score, s);
222                }
223                return score;
224        }
225
226        private void findHighestOrderAxis() {
227                highestOrder = 1;
228                principalAxisIndex = 0;
229                double rmsd  = Double.MAX_VALUE;
230
231                for (int i = 0; i < rotations.size(); i++) {
232                        Rotation s = rotations.get(i);
233                        if (s.getFold() > highestOrder) {
234                                highestOrder = s.getFold();
235                                principalAxisIndex = i;
236                                rmsd = s.getTraceRmsd();
237                        } else if (s.getFold() >= highestOrder && s.getTraceRmsd() < rmsd) {
238                                highestOrder = s.getFold();
239                                principalAxisIndex = i;
240                                rmsd = s.getTraceRmsd();
241                        }
242                }
243        }
244
245        /**
246         * Add E operation to the highest order rotation axis. By definition
247         * E belongs to the highest order axis.
248         */
249        private void setEAxis() {
250                Rotation e = rotations.get(0);
251                Rotation h = rotations.get(principalAxisIndex);
252                e.setAxisAngle(new AxisAngle4d(h.getAxisAngle()));
253                e.getAxisAngle().angle = 0.0;
254                e.setFold(h.getFold());
255        }
256
257        private void findHigherOrderAxes() {
258                higherOrderRotationAxis = 0;
259                for (Rotation s: rotations) {
260                        if (s.getFold() > 2 && s.getDirection() == 1) {
261                                higherOrderRotationAxis++;
262                        }
263                }
264        }
265
266        private void calcAxesDirections() {
267                if (highestOrder == 1) {
268                        for (Rotation s: rotations) {
269                                s.setDirection(0);
270                        }
271                        return;
272                }
273
274                AxisAngle4d pa = rotations.get(principalAxisIndex).getAxisAngle();
275                Vector3d pv = new Vector3d(pa.x, pa.y, pa.z);
276
277                for (Rotation s: rotations) {
278                        AxisAngle4d axis = s.getAxisAngle();
279                        Vector3d av = new Vector3d(axis.x, axis.y, axis.z);
280                        if (Math.abs(pv.dot(av)) > 0.9f) {
281                                // co-linear with principal axis
282                                s.setDirection(0);
283                        } else {
284                                // not co-linear or perpendicular to principal axis
285                                s.setDirection(1);
286                        }
287                }
288                rotations.get(0).setDirection(0); // set the E axis to the principal axis (by definition)
289        }
290
291        private void findTwoFoldsPerpendicular() {
292                twoFoldsPerpendicular = 0;
293                for (Rotation s: rotations) {
294                        if (s.getFold() == 2 && s.getDirection() == 1) {
295                                twoFoldsPerpendicular++;
296                        }
297                }
298        }
299
300
301        public int getHigherOrderRotationAxis(){
302                return higherOrderRotationAxis;
303        }
304
305        public int getTwoFoldsPerpendicular(){
306                return twoFoldsPerpendicular;
307        }
308
309        public int getPrincipalAxisIndex(){
310                return principalAxisIndex;
311        }
312
313        private void calcPointGroup() {
314                complete = false;
315                if (higherOrderRotationAxis > 1) {
316                        // cubic groups
317                        if (highestOrder == 5) {
318                                // rotational icosahedral symmetry or chiral icosahedral symmetry
319                                pointGroup = "I";
320                                complete = rotations.size() == 60;
321                        } else if (highestOrder == 4) {
322                                // rotational octahedral symmetry or chiral octahedral symmetry
323                                pointGroup = "O";
324                                complete = rotations.size() == 24;
325                        } else if (highestOrder == 3) {
326                                // rotational tetrahedral symmetry or chiral tetrahedral symmetry
327                                pointGroup = "T";
328                                complete = rotations.size() == 12;
329                        }
330                } else {
331                        // Cn and Dn groups
332                        // if E is not counted, subtract 1
333                        if (Math.abs(twoFoldsPerpendicular - highestOrder) <= 1 && highestOrder > 1) {
334                                pointGroup = "D" + highestOrder;
335                                complete = rotations.size() == 2 * highestOrder;
336                        } else {
337                                pointGroup = "C" + highestOrder;
338                                complete = rotations.size() == highestOrder;
339                        }
340                }
341
342                if (!complete) {
343                        // the rotation group is incomplete, remove partial results. This happens
344                        // when a structure is symmetric, some subunits are below the rmsd threshold,
345                        // and some are just above the rmsd threshold
346                        int n = 0;
347                        if (rotations.size() > 0) {
348                                n = rotations.get(0).getPermutation().size();
349                                rotations.clear();
350                        }
351                        setC1(n);
352                }
353        }
354
355        public void sortByFoldDecending() {
356                Collections.sort(rotations, new Comparator<Rotation>() {
357                        @Override
358                        public int compare(Rotation o1, Rotation o2) {
359                                int delta = o1.getDirection() - o2.getDirection();
360                                if (delta != 0) {
361                                        return delta;
362                                }
363                                delta = Math.round(Math.signum(o2.getFold() - o1.getFold()));
364                                if (delta != 0) {
365                                        return delta;
366                                }
367
368                                delta = (int)(Math.signum(o1.getAxisAngle().angle - o2.getAxisAngle().angle));
369                                return delta;
370                        }
371                });
372        }
373
374        @Override
375        public Iterator<Rotation> iterator() {
376                return rotations.iterator();
377        }
378}