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 * Created on 08.05.2004
021 *
022 */
023package org.biojava.nbio.structure;
024
025import java.util.ArrayList;
026import java.util.Collection;
027import java.util.List;
028
029import javax.vecmath.Matrix3d;
030import javax.vecmath.Matrix4d;
031import javax.vecmath.Point3d;
032import javax.vecmath.Vector3d;
033
034import org.biojava.nbio.structure.geometry.CalcPoint;
035import org.biojava.nbio.structure.geometry.Matrices;
036import org.biojava.nbio.structure.geometry.SuperPositionSVD;
037import org.biojava.nbio.structure.jama.Matrix;
038import org.slf4j.Logger;
039import org.slf4j.LoggerFactory;
040
041/**
042 * Utility operations on Atoms, AminoAcids, Matrices, Point3d, etc.
043 * <p>
044 * Currently the coordinates of an Atom are stored as an array of size 3
045 * (double[3]). It would be more powerful to use Point3d from javax.vecmath.
046 * Overloaded methods for Point3d operations exist in the {@link CalcPoint}
047 * Class.
048 *
049 * @author Andreas Prlic
050 * @author Aleix Lafita
051 * @since 1.4
052 * @version %I% %G%
053 */
054
055public class Calc {
056
057        private final static Logger logger = LoggerFactory.getLogger(Calc.class);
058
059        /**
060         * calculate distance between two atoms.
061         *
062         * @param a
063         *            an Atom object
064         * @param b
065         *            an Atom object
066         * @return a double
067         */
068        public static final double getDistance(Atom a, Atom b) {
069                double x = a.getX() - b.getX();
070                double y = a.getY() - b.getY();
071                double z = a.getZ() - b.getZ();
072
073                double s = x * x + y * y + z * z;
074
075                return Math.sqrt(s);
076        }
077
078        /**
079         * Will calculate the square of distances between two atoms. This will be
080         * faster as it will not perform the final square root to get the actual
081         * distance. Use this if doing large numbers of distance comparisons - it is
082         * marginally faster than getDistance().
083         *
084         * @param a
085         *            an Atom object
086         * @param b
087         *            an Atom object
088         * @return a double
089         */
090        public static double getDistanceFast(Atom a, Atom b) {
091                double x = a.getX() - b.getX();
092                double y = a.getY() - b.getY();
093                double z = a.getZ() - b.getZ();
094
095                return x * x + y * y + z * z;
096        }
097
098        public static final Atom invert(Atom a) {
099                double[] coords = new double[] { 0.0, 0.0, 0.0 };
100                Atom zero = new AtomImpl();
101                zero.setCoords(coords);
102                return subtract(zero, a);
103        }
104
105        /**
106         * add two atoms ( a + b).
107         *
108         * @param a
109         *            an Atom object
110         * @param b
111         *            an Atom object
112         * @return an Atom object
113         */
114        public static final Atom add(Atom a, Atom b) {
115
116                Atom c = new AtomImpl();
117                c.setX(a.getX() + b.getX());
118                c.setY(a.getY() + b.getY());
119                c.setZ(a.getZ() + b.getZ());
120
121                return c;
122        }
123
124        /**
125         * subtract two atoms ( a - b).
126         *
127         * @param a
128         *            an Atom object
129         * @param b
130         *            an Atom object
131         * @return n new Atom object representing the difference
132         */
133        public static final Atom subtract(Atom a, Atom b) {
134                Atom c = new AtomImpl();
135                c.setX(a.getX() - b.getX());
136                c.setY(a.getY() - b.getY());
137                c.setZ(a.getZ() - b.getZ());
138
139                return c;
140        }
141
142        /**
143         * Vector product (cross product).
144         *
145         * @param a
146         *            an Atom object
147         * @param b
148         *            an Atom object
149         * @return an Atom object
150         */
151        public static final Atom vectorProduct(Atom a, Atom b) {
152
153                Atom c = new AtomImpl();
154                c.setX(a.getY() * b.getZ() - a.getZ() * b.getY());
155                c.setY(a.getZ() * b.getX() - a.getX() * b.getZ());
156                c.setZ(a.getX() * b.getY() - a.getY() * b.getX());
157                return c;
158
159        }
160
161        /**
162         * Scalar product (dot product).
163         *
164         * @param a
165         *            an Atom object
166         * @param b
167         *            an Atom object
168         * @return a double
169         */
170        public static final double scalarProduct(Atom a, Atom b) {
171                return a.getX() * b.getX() + a.getY() * b.getY() + a.getZ() * b.getZ();
172        }
173
174        /**
175         * Gets the length of the vector (2-norm)
176         *
177         * @param a
178         *            an Atom object
179         * @return Square root of the sum of the squared elements
180         */
181        public static final double amount(Atom a) {
182                return Math.sqrt(scalarProduct(a, a));
183        }
184
185        /**
186         * Gets the angle between two vectors
187         *
188         * @param a
189         *            an Atom object
190         * @param b
191         *            an Atom object
192         * @return Angle between a and b in degrees, in range [0,180]. If either
193         *         vector has length 0 then angle is not defined and NaN is returned
194         */
195        public static final double angle(Atom a, Atom b) {
196
197                Vector3d va = new Vector3d(a.getCoordsAsPoint3d());
198                Vector3d vb = new Vector3d(b.getCoordsAsPoint3d());
199
200                return Math.toDegrees(va.angle(vb));
201
202        }
203
204        /**
205         * Returns the unit vector of vector a .
206         *
207         * @param a
208         *            an Atom object
209         * @return an Atom object
210         */
211        public static final Atom unitVector(Atom a) {
212                double amount = amount(a);
213
214                double[] coords = new double[3];
215
216                coords[0] = a.getX() / amount;
217                coords[1] = a.getY() / amount;
218                coords[2] = a.getZ() / amount;
219
220                a.setCoords(coords);
221                return a;
222
223        }
224
225        /**
226         * Calculate the torsion angle, i.e. the angle between the normal vectors of
227         * the two plains a-b-c and b-c-d. See
228         * http://en.wikipedia.org/wiki/Dihedral_angle
229         * 
230         * @param a
231         *            an Atom object
232         * @param b
233         *            an Atom object
234         * @param c
235         *            an Atom object
236         * @param d
237         *            an Atom object
238         * @return the torsion angle in degrees, in range +-[0,180]. If either first
239         *         3 or last 3 atoms are colinear then torsion angle is not defined
240         *         and NaN is returned
241         */
242        public static final double torsionAngle(Atom a, Atom b, Atom c, Atom d) {
243
244                Atom ab = subtract(a, b);
245                Atom cb = subtract(c, b);
246                Atom bc = subtract(b, c);
247                Atom dc = subtract(d, c);
248
249                Atom abc = vectorProduct(ab, cb);
250                Atom bcd = vectorProduct(bc, dc);
251
252                double angl = angle(abc, bcd);
253
254                /* calc the sign: */
255                Atom vecprod = vectorProduct(abc, bcd);
256                double val = scalarProduct(cb, vecprod);
257                if (val < 0.0)
258                        angl = -angl;
259
260                return angl;
261        }
262
263        /**
264         * Calculate the phi angle.
265         *
266         * @param a
267         *            an AminoAcid object
268         * @param b
269         *            an AminoAcid object
270         * @return a double
271         * @throws StructureException
272         *             if aminoacids not connected or if any of the 4 needed atoms
273         *             missing
274         */
275        public static final double getPhi(AminoAcid a, AminoAcid b)
276                        throws StructureException {
277
278                if (!isConnected(a, b)) {
279                        throw new StructureException(
280                                        "can not calc Phi - AminoAcids are not connected!");
281                }
282
283                Atom a_C = a.getC();
284                Atom b_N = b.getN();
285                Atom b_CA = b.getCA();
286                Atom b_C = b.getC();
287
288                // C and N were checked in isConnected already
289                if (b_CA == null)
290                        throw new StructureException(
291                                        "Can not calculate Phi, CA atom is missing");
292
293                return torsionAngle(a_C, b_N, b_CA, b_C);
294        }
295
296        /**
297         * Calculate the psi angle.
298         *
299         * @param a
300         *            an AminoAcid object
301         * @param b
302         *            an AminoAcid object
303         * @return a double
304         * @throws StructureException
305         *             if aminoacids not connected or if any of the 4 needed atoms
306         *             missing
307         */
308        public static final double getPsi(AminoAcid a, AminoAcid b)
309                        throws StructureException {
310                if (!isConnected(a, b)) {
311                        throw new StructureException(
312                                        "can not calc Psi - AminoAcids are not connected!");
313                }
314
315                Atom a_N = a.getN();
316                Atom a_CA = a.getCA();
317                Atom a_C = a.getC();
318                Atom b_N = b.getN();
319
320                // C and N were checked in isConnected already
321                if (a_CA == null)
322                        throw new StructureException(
323                                        "Can not calculate Psi, CA atom is missing");
324
325                return torsionAngle(a_N, a_CA, a_C, b_N);
326
327        }
328
329        /**
330         * Test if two amino acids are connected, i.e. if the distance from C to N <
331         * 2.5 Angstrom.
332         *
333         * If one of the AminoAcids has an atom missing, returns false.
334         *
335         * @param a
336         *            an AminoAcid object
337         * @param b
338         *            an AminoAcid object
339         * @return true if ...
340         */
341        public static final boolean isConnected(AminoAcid a, AminoAcid b) {
342                Atom C = null;
343                Atom N = null;
344
345                C = a.getC();
346                N = b.getN();
347
348                if (C == null || N == null)
349                        return false;
350
351                // one could also check if the CA atoms are < 4 A...
352                double distance = getDistance(C, N);
353                return distance < 2.5;
354        }
355
356        /**
357         * Rotate a single Atom aroud a rotation matrix. The rotation Matrix must be
358         * a pre-multiplication 3x3 Matrix.
359         *
360         * If the matrix is indexed m[row][col], then the matrix will be
361         * pre-multiplied (y=atom*M)
362         * 
363         * @param atom
364         *            atom to be rotated
365         * @param m
366         *            a rotation matrix represented as a double[3][3] array
367         */
368        public static final void rotate(Atom atom, double[][] m) {
369
370                double x = atom.getX();
371                double y = atom.getY();
372                double z = atom.getZ();
373
374                double nx = m[0][0] * x + m[0][1] * y + m[0][2] * z;
375                double ny = m[1][0] * x + m[1][1] * y + m[1][2] * z;
376                double nz = m[2][0] * x + m[2][1] * y + m[2][2] * z;
377
378                atom.setX(nx);
379                atom.setY(ny);
380                atom.setZ(nz);
381        }
382
383        /**
384         * Rotate a structure. The rotation Matrix must be a pre-multiplication
385         * Matrix.
386         *
387         * @param structure
388         *            a Structure object
389         * @param rotationmatrix
390         *            an array (3x3) of double representing the rotation matrix.
391         * @throws StructureException
392         *             ...
393         */
394        public static final void rotate(Structure structure,
395                        double[][] rotationmatrix) throws StructureException {
396
397                if (rotationmatrix.length != 3) {
398                        throw new StructureException("matrix does not have size 3x3 !");
399                }
400                AtomIterator iter = new AtomIterator(structure);
401                while (iter.hasNext()) {
402                        Atom atom = iter.next();
403                        Calc.rotate(atom, rotationmatrix);
404                }
405        }
406
407        /**
408         * Rotate a Group. The rotation Matrix must be a pre-multiplication Matrix.
409         *
410         * @param group
411         *            a group object
412         * @param rotationmatrix
413         *            an array (3x3) of double representing the rotation matrix.
414         * @throws StructureException
415         *             ...
416         */
417        public static final void rotate(Group group, double[][] rotationmatrix)
418                        throws StructureException {
419
420                if (rotationmatrix.length != 3) {
421                        throw new StructureException("matrix does not have size 3x3 !");
422                }
423                AtomIterator iter = new AtomIterator(group);
424                while (iter.hasNext()) {
425                        Atom atom = null;
426
427                        atom = iter.next();
428                        rotate(atom, rotationmatrix);
429
430                }
431        }
432
433        /**
434         * Rotate an Atom around a Matrix object. The rotation Matrix must be a
435         * pre-multiplication Matrix.
436         *
437         * @param atom
438         *            atom to be rotated
439         * @param m
440         *            rotation matrix to be applied to the atom
441         */
442        public static final void rotate(Atom atom, Matrix m) {
443
444                double x = atom.getX();
445                double y = atom.getY();
446                double z = atom.getZ();
447                double[][] ad = new double[][] { { x, y, z } };
448
449                Matrix am = new Matrix(ad);
450                Matrix na = am.times(m);
451
452                atom.setX(na.get(0, 0));
453                atom.setY(na.get(0, 1));
454                atom.setZ(na.get(0, 2));
455
456        }
457
458        /**
459         * Rotate a group object. The rotation Matrix must be a pre-multiplication
460         * Matrix.
461         *
462         * @param group
463         *            a group to be rotated
464         * @param m
465         *            a Matrix object representing the rotation matrix
466         */
467        public static final void rotate(Group group, Matrix m) {
468
469                AtomIterator iter = new AtomIterator(group);
470
471                while (iter.hasNext()) {
472                        Atom atom = iter.next();
473                        rotate(atom, m);
474
475                }
476
477        }
478
479        /**
480         * Rotate a structure object. The rotation Matrix must be a
481         * pre-multiplication Matrix.
482         *
483         * @param structure
484         *            the structure to be rotated
485         * @param m
486         *            rotation matrix to be applied
487         */
488        public static final void rotate(Structure structure, Matrix m) {
489
490                AtomIterator iter = new AtomIterator(structure);
491
492                while (iter.hasNext()) {
493                        Atom atom = iter.next();
494                        rotate(atom, m);
495
496                }
497
498        }
499
500        /**
501         * Transform an array of atoms at once. The transformation Matrix must be a
502         * post-multiplication Matrix.
503         *
504         * @param ca
505         *            array of Atoms to shift
506         * @param t
507         *            transformation Matrix4d
508         */
509        public static void transform(Atom[] ca, Matrix4d t) {
510                for (Atom atom : ca)
511                        Calc.transform(atom, t);
512        }
513
514        /**
515         * Transforms an atom object, given a Matrix4d (i.e. the vecmath library
516         * double-precision 4x4 rotation+translation matrix). The transformation
517         * Matrix must be a post-multiplication Matrix.
518         *
519         * @param atom
520         * @param m
521         */
522        public static final void transform(Atom atom, Matrix4d m) {
523
524                Point3d p = new Point3d(atom.getX(), atom.getY(), atom.getZ());
525                m.transform(p);
526
527                atom.setX(p.x);
528                atom.setY(p.y);
529                atom.setZ(p.z);
530        }
531
532        /**
533         * Transforms a group object, given a Matrix4d (i.e. the vecmath library
534         * double-precision 4x4 rotation+translation matrix). The transformation
535         * Matrix must be a post-multiplication Matrix.
536         *
537         * @param group
538         * @param m
539         */
540        public static final void transform(Group group, Matrix4d m) {
541                for (Atom atom : group.getAtoms()) {
542                        transform(atom, m);
543                }
544                for (Group altG : group.getAltLocs()) {
545                        transform(altG, m);
546                }
547        }
548
549        /**
550         * Transforms a structure object, given a Matrix4d (i.e. the vecmath library
551         * double-precision 4x4 rotation+translation matrix). The transformation
552         * Matrix must be a post-multiplication Matrix.
553         *
554         * @param structure
555         * @param m
556         */
557        public static final void transform(Structure structure, Matrix4d m) {
558                for (int n=0; n<structure.nrModels();n++) {
559                        for (Chain c : structure.getChains(n)) {
560                                transform(c, m);
561                        }
562                }
563        }
564
565        /**
566         * Transforms a chain object, given a Matrix4d (i.e. the vecmath library
567         * double-precision 4x4 rotation+translation matrix). The transformation
568         * Matrix must be a post-multiplication Matrix.
569         *
570         * @param chain
571         * @param m
572         */
573        public static final void transform(Chain chain, Matrix4d m) {
574
575                for (Group g : chain.getAtomGroups()) {
576                        transform(g, m);
577                }
578        }
579
580        /**
581         * Translates an atom object, given a Vector3d (i.e. the vecmath library
582         * double-precision 3-d vector)
583         * 
584         * @param atom
585         * @param v
586         */
587        public static final void translate(Atom atom, Vector3d v) {
588
589                atom.setX(atom.getX() + v.x);
590                atom.setY(atom.getY() + v.y);
591                atom.setZ(atom.getZ() + v.z);
592        }
593
594        /**
595         * Translates a group object, given a Vector3d (i.e. the vecmath library
596         * double-precision 3-d vector)
597         * 
598         * @param group
599         * @param v
600         */
601        public static final void translate(Group group, Vector3d v) {
602
603                for (Atom atom : group.getAtoms()) {
604                        translate(atom, v);
605                }
606                for (Group altG : group.getAltLocs()) {
607                        translate(altG, v);
608                }
609        }
610
611        /**
612         * Translates a chain object, given a Vector3d (i.e. the vecmath library
613         * double-precision 3-d vector)
614         * 
615         * @param chain
616         * @param v
617         */
618        public static final void translate(Chain chain, Vector3d v) {
619
620                for (Group g : chain.getAtomGroups()) {
621                        translate(g, v);
622                }
623        }
624
625        /**
626         * Translates a Structure object, given a Vector3d (i.e. the vecmath library
627         * double-precision 3-d vector)
628         * 
629         * @param structure
630         * @param v
631         */
632        public static final void translate(Structure structure, Vector3d v) {
633                
634                for (int n=0; n<structure.nrModels();n++) {
635                        for (Chain c : structure.getChains(n)) {
636                                translate(c, v);
637                        }
638                }
639        }
640
641        /**
642         * calculate structure + Matrix coodinates ...
643         *
644         * @param s
645         *            the structure to operate on
646         * @param matrix
647         *            a Matrix object
648         */
649        public static final void plus(Structure s, Matrix matrix) {
650                AtomIterator iter = new AtomIterator(s);
651                Atom oldAtom = null;
652                Atom rotOldAtom = null;
653                while (iter.hasNext()) {
654                        Atom atom = null;
655
656                        atom = iter.next();
657                        try {
658                                if (oldAtom != null) {
659                                        logger.debug("before {}", getDistance(oldAtom, atom));
660                                }
661                        } catch (Exception e) {
662                                logger.error("Exception: ", e);
663                        }
664                        oldAtom = (Atom) atom.clone();
665
666                        double x = atom.getX();
667                        double y = atom.getY();
668                        double z = atom.getZ();
669                        double[][] ad = new double[][] { { x, y, z } };
670
671                        Matrix am = new Matrix(ad);
672                        Matrix na = am.plus(matrix);
673
674                        double[] coords = new double[3];
675                        coords[0] = na.get(0, 0);
676                        coords[1] = na.get(0, 1);
677                        coords[2] = na.get(0, 2);
678                        atom.setCoords(coords);
679                        try {
680                                if (rotOldAtom != null) {
681                                        logger.debug("after {}", getDistance(rotOldAtom, atom));
682                                }
683                        } catch (Exception e) {
684                                logger.error("Exception: ", e);
685                        }
686                        rotOldAtom = (Atom) atom.clone();
687                }
688
689        }
690
691        /**
692         * shift a structure with a vector.
693         *
694         * @param structure
695         *            a Structure object
696         * @param a
697         *            an Atom object representing a shift vector
698         */
699        public static final void shift(Structure structure, Atom a) {
700
701                AtomIterator iter = new AtomIterator(structure);
702                while (iter.hasNext()) {
703                        Atom atom = null;
704
705                        atom = iter.next();
706
707                        Atom natom = add(atom, a);
708                        double x = natom.getX();
709                        double y = natom.getY();
710                        double z = natom.getZ();
711                        atom.setX(x);
712                        atom.setY(y);
713                        atom.setZ(z);
714
715                }
716        }
717
718        /**
719         * Shift a vector.
720         *
721         * @param a
722         *            vector a
723         * @param b
724         *            vector b
725         */
726        public static final void shift(Atom a, Atom b) {
727
728                Atom natom = add(a, b);
729                double x = natom.getX();
730                double y = natom.getY();
731                double z = natom.getZ();
732                a.setX(x);
733                a.setY(y);
734                a.setZ(z);
735        }
736
737        /**
738         * Shift a Group with a vector.
739         *
740         * @param group
741         *            a group object
742         * @param a
743         *            an Atom object representing a shift vector
744         */
745        public static final void shift(Group group, Atom a) {
746
747                AtomIterator iter = new AtomIterator(group);
748                while (iter.hasNext()) {
749                        Atom atom = null;
750
751                        atom = iter.next();
752
753                        Atom natom = add(atom, a);
754                        double x = natom.getX();
755                        double y = natom.getY();
756                        double z = natom.getZ();
757                        atom.setX(x);
758                        atom.setY(y);
759                        atom.setZ(z);
760
761                }
762        }
763
764        /**
765         * Returns the centroid of the set of atoms.
766         * 
767         * @param atomSet
768         *            a set of Atoms
769         * @return an Atom representing the Centroid of the set of atoms
770         */
771        public static final Atom getCentroid(Atom[] atomSet) {
772                
773                // if we don't catch this case, the centroid returned is (NaN,NaN,NaN), which can cause lots of problems down the line
774                if (atomSet.length==0) 
775                        throw new IllegalArgumentException("Atom array has length 0, can't calculate centroid!");
776
777
778                double[] coords = new double[3];
779
780                coords[0] = 0;
781                coords[1] = 0;
782                coords[2] = 0;
783
784                for (Atom a : atomSet) {
785                        coords[0] += a.getX();
786                        coords[1] += a.getY();
787                        coords[2] += a.getZ();
788                }
789
790                int n = atomSet.length;
791                coords[0] = coords[0] / n;
792                coords[1] = coords[1] / n;
793                coords[2] = coords[2] / n;
794
795                Atom vec = new AtomImpl();
796                vec.setCoords(coords);
797                return vec;
798
799        }
800
801        /**
802         * Returns the center of mass of the set of atoms. Atomic masses of the
803         * Atoms are used.
804         * 
805         * @param points
806         *            a set of Atoms
807         * @return an Atom representing the center of mass
808         */
809        public static Atom centerOfMass(Atom[] points) {
810                Atom center = new AtomImpl();
811
812                float totalMass = 0.0f;
813                for (Atom a : points) {
814                        float mass = a.getElement().getAtomicMass();
815                        totalMass += mass;
816                        center = scaleAdd(mass, a, center);
817                }
818
819                center = scaleEquals(center, 1.0f / totalMass);
820                return center;
821        }
822
823        /**
824         * Multiply elements of a by s (in place)
825         * 
826         * @param a
827         * @param s
828         * @return the modified a
829         */
830        public static Atom scaleEquals(Atom a, double s) {
831                double x = a.getX();
832                double y = a.getY();
833                double z = a.getZ();
834
835                x *= s;
836                y *= s;
837                z *= s;
838
839                // Atom b = new AtomImpl();
840                a.setX(x);
841                a.setY(y);
842                a.setZ(z);
843
844                return a;
845        }
846
847        /**
848         * Multiply elements of a by s
849         * 
850         * @param a
851         * @param s
852         * @return A new Atom with s*a
853         */
854        public static Atom scale(Atom a, double s) {
855                double x = a.getX();
856                double y = a.getY();
857                double z = a.getZ();
858
859                Atom b = new AtomImpl();
860                b.setX(x * s);
861                b.setY(y * s);
862                b.setZ(z * s);
863
864                return b;
865        }
866
867        /**
868         * Perform linear transformation s*X+B, and store the result in b
869         * 
870         * @param s
871         *            Amount to scale x
872         * @param x
873         *            Input coordinate
874         * @param b
875         *            Vector to translate (will be modified)
876         * @return b, after modification
877         */
878        public static Atom scaleAdd(double s, Atom x, Atom b) {
879
880                double xc = s * x.getX() + b.getX();
881                double yc = s * x.getY() + b.getY();
882                double zc = s * x.getZ() + b.getZ();
883
884                // Atom a = new AtomImpl();
885                b.setX(xc);
886                b.setY(yc);
887                b.setZ(zc);
888
889                return b;
890        }
891
892        /**
893         * Returns the Vector that needs to be applied to shift a set of atoms to
894         * the Centroid.
895         * 
896         * @param atomSet
897         *            array of Atoms
898         * @return the vector needed to shift the set of atoms to its geometric
899         *         center
900         */
901        public static final Atom getCenterVector(Atom[] atomSet) {
902                Atom centroid = getCentroid(atomSet);
903
904                return getCenterVector(atomSet, centroid);
905
906        }
907
908        /**
909         * Returns the Vector that needs to be applied to shift a set of atoms to
910         * the Centroid, if the centroid is already known
911         * 
912         * @param atomSet
913         *            array of Atoms
914         * @return the vector needed to shift the set of atoms to its geometric
915         *         center
916         */
917        public static final Atom getCenterVector(Atom[] atomSet, Atom centroid) {
918
919                double[] coords = new double[3];
920                coords[0] = 0 - centroid.getX();
921                coords[1] = 0 - centroid.getY();
922                coords[2] = 0 - centroid.getZ();
923
924                Atom shiftVec = new AtomImpl();
925                shiftVec.setCoords(coords);
926                return shiftVec;
927
928        }
929
930        /**
931         * Center the atoms at the Centroid.
932         * 
933         * @param atomSet
934         *            a set of Atoms
935         * @return an Atom representing the Centroid of the set of atoms
936         * @throws StructureException
937         * */
938        public static final Atom[] centerAtoms(Atom[] atomSet)
939                        throws StructureException {
940
941                Atom centroid = getCentroid(atomSet);
942                return centerAtoms(atomSet, centroid);
943        }
944
945        /**
946         * Center the atoms at the Centroid, if the centroid is already know.
947         * 
948         * @param atomSet
949         *            a set of Atoms
950         * @return an Atom representing the Centroid of the set of atoms
951         * @throws StructureException
952         * */
953        public static final Atom[] centerAtoms(Atom[] atomSet, Atom centroid)
954                        throws StructureException {
955
956                Atom shiftVector = getCenterVector(atomSet, centroid);
957
958                Atom[] newAtoms = new AtomImpl[atomSet.length];
959
960                for (int i = 0; i < atomSet.length; i++) {
961                        Atom a = atomSet[i];
962                        Atom n = add(a, shiftVector);
963                        newAtoms[i] = n;
964                }
965                return newAtoms;
966        }
967
968        /**
969         * creates a virtual C-beta atom. this might be needed when working with GLY
970         *
971         * thanks to Peter Lackner for a python template of this method.
972         * 
973         * @param amino
974         *            the amino acid for which a "virtual" CB atom should be
975         *            calculated
976         * @return a "virtual" CB atom
977         * @throws StructureException
978         */
979        public static final Atom createVirtualCBAtom(AminoAcid amino)
980                        throws StructureException {
981
982                AminoAcid ala = StandardAminoAcid.getAminoAcid("ALA");
983                Atom aN = ala.getN();
984                Atom aCA = ala.getCA();
985                Atom aC = ala.getC();
986                Atom aCB = ala.getCB();
987
988                Atom[] arr1 = new Atom[3];
989                arr1[0] = aN;
990                arr1[1] = aCA;
991                arr1[2] = aC;
992
993                Atom[] arr2 = new Atom[3];
994                arr2[0] = amino.getN();
995                arr2[1] = amino.getCA();
996                arr2[2] = amino.getC();
997
998                // ok now we got the two arrays, do a Superposition:
999
1000                SuperPositionSVD svd = new SuperPositionSVD(false);
1001
1002                Matrix4d transform = svd.superpose(Calc.atomsToPoints(arr1), Calc.atomsToPoints(arr2));
1003                Matrix rotMatrix = Matrices.getRotationJAMA(transform);
1004                Atom tranMatrix = getTranslationVector(transform);
1005                
1006                Calc.rotate(aCB, rotMatrix);
1007
1008                Atom virtualCB = Calc.add(aCB, tranMatrix);
1009                virtualCB.setName("CB");
1010
1011                return virtualCB;
1012        }
1013
1014        /**
1015         * Gets euler angles for a matrix given in ZYZ convention. (as e.g. used by
1016         * Jmol)
1017         *
1018         * @param m
1019         *            the rotation matrix
1020         * @return the euler values for a rotation around Z, Y, Z in degrees...
1021         */
1022        public static final double[] getZYZEuler(Matrix m) {
1023                double m22 = m.get(2, 2);
1024                double rY = Math.toDegrees(Math.acos(m22));
1025                double rZ1, rZ2;
1026                if (m22 > .999d || m22 < -.999d) {
1027                        rZ1 = Math.toDegrees(Math.atan2(m.get(1, 0), m.get(1, 1)));
1028                        rZ2 = 0;
1029                } else {
1030                        rZ1 = Math.toDegrees(Math.atan2(m.get(2, 1), -m.get(2, 0)));
1031                        rZ2 = Math.toDegrees(Math.atan2(m.get(1, 2), m.get(0, 2)));
1032                }
1033                return new double[] { rZ1, rY, rZ2 };
1034        }
1035
1036        /**
1037         * Convert a rotation Matrix to Euler angles. This conversion uses
1038         * conventions as described on page:
1039         * http://www.euclideanspace.com/maths/geometry/rotations/euler/index.htm
1040         * Coordinate System: right hand Positive angle: right hand Order of euler
1041         * angles: heading first, then attitude, then bank
1042         *
1043         * @param m
1044         *            the rotation matrix
1045         * @return a array of three doubles containing the three euler angles in
1046         *         radians
1047         */
1048        public static final double[] getXYZEuler(Matrix m) {
1049                double heading, attitude, bank;
1050
1051                // Assuming the angles are in radians.
1052                if (m.get(1, 0) > 0.998) { // singularity at north pole
1053                        heading = Math.atan2(m.get(0, 2), m.get(2, 2));
1054                        attitude = Math.PI / 2;
1055                        bank = 0;
1056
1057                } else if (m.get(1, 0) < -0.998) { // singularity at south pole
1058                        heading = Math.atan2(m.get(0, 2), m.get(2, 2));
1059                        attitude = -Math.PI / 2;
1060                        bank = 0;
1061
1062                } else {
1063                        heading = Math.atan2(-m.get(2, 0), m.get(0, 0));
1064                        bank = Math.atan2(-m.get(1, 2), m.get(1, 1));
1065                        attitude = Math.asin(m.get(1, 0));
1066                }
1067                return new double[] { heading, attitude, bank };
1068        }
1069
1070        /**
1071         * This conversion uses NASA standard aeroplane conventions as described on
1072         * page:
1073         * http://www.euclideanspace.com/maths/geometry/rotations/euler/index.htm
1074         * Coordinate System: right hand Positive angle: right hand Order of euler
1075         * angles: heading first, then attitude, then bank. matrix row column
1076         * ordering: [m00 m01 m02] [m10 m11 m12] [m20 m21 m22]
1077         * 
1078         * @param heading
1079         *            in radians
1080         * @param attitude
1081         *            in radians
1082         * @param bank
1083         *            in radians
1084         * @return the rotation matrix
1085         */
1086        public static final Matrix matrixFromEuler(double heading, double attitude,
1087                        double bank) {
1088                // Assuming the angles are in radians.
1089                double ch = Math.cos(heading);
1090                double sh = Math.sin(heading);
1091                double ca = Math.cos(attitude);
1092                double sa = Math.sin(attitude);
1093                double cb = Math.cos(bank);
1094                double sb = Math.sin(bank);
1095
1096                Matrix m = new Matrix(3, 3);
1097                m.set(0, 0, ch * ca);
1098                m.set(0, 1, sh * sb - ch * sa * cb);
1099                m.set(0, 2, ch * sa * sb + sh * cb);
1100                m.set(1, 0, sa);
1101                m.set(1, 1, ca * cb);
1102                m.set(1, 2, -ca * sb);
1103                m.set(2, 0, -sh * ca);
1104                m.set(2, 1, sh * sa * cb + ch * sb);
1105                m.set(2, 2, -sh * sa * sb + ch * cb);
1106
1107                return m;
1108        }
1109
1110        /**
1111         * Calculates the angle from centerPt to targetPt in degrees. The return
1112         * should range from [0,360), rotating CLOCKWISE, 0 and 360 degrees
1113         * represents NORTH, 90 degrees represents EAST, etc...
1114         *
1115         * Assumes all points are in the same coordinate space. If they are not, you
1116         * will need to call SwingUtilities.convertPointToScreen or equivalent on
1117         * all arguments before passing them to this function.
1118         *
1119         * @param centerPt
1120         *            Point we are rotating around.
1121         * @param targetPt
1122         *            Point we want to calculate the angle to.
1123         * @return angle in degrees. This is the angle from centerPt to targetPt.
1124         */
1125        public static double calcRotationAngleInDegrees(Atom centerPt, Atom targetPt) {
1126                // calculate the angle theta from the deltaY and deltaX values
1127                // (atan2 returns radians values from [-PI,PI])
1128                // 0 currently points EAST.
1129                // NOTE: By preserving Y and X param order to atan2, we are expecting
1130                // a CLOCKWISE angle direction.
1131                double theta = Math.atan2(targetPt.getY() - centerPt.getY(),
1132                                targetPt.getX() - centerPt.getX());
1133
1134                // rotate the theta angle clockwise by 90 degrees
1135                // (this makes 0 point NORTH)
1136                // NOTE: adding to an angle rotates it clockwise.
1137                // subtracting would rotate it counter-clockwise
1138                theta += Math.PI / 2.0;
1139
1140                // convert from radians to degrees
1141                // this will give you an angle from [0->270],[-180,0]
1142                double angle = Math.toDegrees(theta);
1143
1144                // convert to positive range [0-360)
1145                // since we want to prevent negative angles, adjust them now.
1146                // we can assume that atan2 will not return a negative value
1147                // greater than one partial rotation
1148                if (angle < 0) {
1149                        angle += 360;
1150                }
1151
1152                return angle;
1153        }
1154
1155        public static void main(String[] args) {
1156                Atom a = new AtomImpl();
1157                a.setX(0);
1158                a.setY(0);
1159                a.setZ(0);
1160
1161                Atom b = new AtomImpl();
1162                b.setX(1);
1163                b.setY(1);
1164                b.setZ(0);
1165
1166                logger.info("Angle between atoms: ", calcRotationAngleInDegrees(a, b));
1167        }
1168
1169        public static void rotate(Atom[] ca, Matrix matrix) {
1170                for (Atom atom : ca)
1171                        Calc.rotate(atom, matrix);
1172        }
1173
1174        /**
1175         * Shift an array of atoms at once.
1176         * 
1177         * @param ca
1178         *            array of Atoms to shift
1179         * @param b
1180         *            reference Atom vector
1181         */
1182        public static void shift(Atom[] ca, Atom b) {
1183                for (Atom atom : ca)
1184                        Calc.shift(atom, b);
1185        }
1186
1187        /**
1188         * Convert JAMA rotation and translation to a Vecmath transformation matrix.
1189         * Because the JAMA matrix is a pre-multiplication matrix and the Vecmath
1190         * matrix is a post-multiplication one, the rotation matrix is transposed to
1191         * ensure that the transformation they produce is the same.
1192         *
1193         * @param rot
1194         *            3x3 Rotation matrix
1195         * @param trans
1196         *            3x1 translation vector in Atom coordinates
1197         * @return 4x4 transformation matrix
1198         */
1199        public static Matrix4d getTransformation(Matrix rot, Atom trans) {
1200                return new Matrix4d(new Matrix3d(rot.getColumnPackedCopy()),
1201                                new Vector3d(trans.getCoordsAsPoint3d()), 1.0);
1202        }
1203        
1204        /**
1205         * Extract the translational vector as an Atom of a transformation matrix.
1206         *
1207         * @param transform
1208         *            Matrix4d
1209         * @return Atom shift vector
1210         */
1211        public static Atom getTranslationVector(Matrix4d transform) {
1212
1213                Atom transl = new AtomImpl();
1214                double[] coords = { transform.m03, transform.m13, transform.m23 };
1215                transl.setCoords(coords);
1216                return transl;
1217        }
1218
1219        /**
1220         * Convert an array of atoms into an array of vecmath points
1221         * 
1222         * @param atoms
1223         *            list of atoms
1224         * @return list of Point3ds storing the x,y,z coordinates of each atom
1225         */
1226        public static Point3d[] atomsToPoints(Atom[] atoms) {
1227                Point3d[] points = new Point3d[atoms.length];
1228                for (int i = 0; i < atoms.length; i++) {
1229                        points[i] = atoms[i].getCoordsAsPoint3d();
1230                }
1231                return points;
1232        }
1233        /**
1234         * Convert an array of atoms into an array of vecmath points
1235         * 
1236         * @param atoms
1237         *            list of atoms
1238         * @return list of Point3ds storing the x,y,z coordinates of each atom
1239         */
1240        public static List<Point3d> atomsToPoints(Collection<Atom> atoms) {
1241                ArrayList<Point3d> points = new ArrayList<>(atoms.size());
1242                for (Atom atom : atoms) {
1243                        points.add(atom.getCoordsAsPoint3d());
1244                }
1245                return points;
1246        }
1247
1248        /**
1249         * Calculate the RMSD of two Atom arrays, already superposed.
1250         * 
1251         * @param x
1252         *            array of Atoms superposed to y
1253         * @param y
1254         *            array of Atoms superposed to x
1255         * @return RMSD
1256         */
1257        public static double rmsd(Atom[] x, Atom[] y) {
1258                return CalcPoint.rmsd(atomsToPoints(x), atomsToPoints(y));
1259        }
1260
1261        /**
1262         * Calculate the TM-Score for the superposition.
1263         *
1264         * <em>Normalizes by the <strong>minimum</strong>-length structure (that is, {@code min\{len1,len2\}}).</em>
1265         *
1266         * Atom sets must be pre-rotated.
1267         *
1268         * <p>
1269         * Citation:<br/>
1270         * <i>Zhang Y and Skolnick J (2004). "Scoring function for automated
1271         * assessment of protein structure template quality". Proteins 57: 702 -
1272         * 710.</i>
1273         *
1274         * @param atomSet1
1275         *            atom array 1
1276         * @param atomSet2
1277         *            atom array 2
1278         * @param len1
1279         *            The full length of the protein supplying atomSet1
1280         * @param len2
1281         *            The full length of the protein supplying atomSet2
1282         * @return The TM-Score
1283         * @throws StructureException
1284         */
1285        public static double getTMScore(Atom[] atomSet1, Atom[] atomSet2, int len1, int len2) throws StructureException {
1286                return getTMScore(atomSet1, atomSet2, len1, len2,true);
1287        }
1288
1289        /**
1290         * Calculate the TM-Score for the superposition.
1291         *
1292         * Atom sets must be pre-rotated.
1293         *
1294         * <p>
1295         * Citation:<br/>
1296         * <i>Zhang Y and Skolnick J (2004). "Scoring function for automated
1297         * assessment of protein structure template quality". Proteins 57: 702 -
1298         * 710.</i>
1299         *
1300         * @param atomSet1
1301         *            atom array 1
1302         * @param atomSet2
1303         *            atom array 2
1304         * @param len1
1305         *            The full length of the protein supplying atomSet1
1306         * @param len2
1307         *            The full length of the protein supplying atomSet2
1308         * @param normalizeMin
1309         *            Whether to normalize by the <strong>minimum</strong>-length structure,
1310         *            that is, {@code min\{len1,len2\}}. If false, normalized by the {@code max\{len1,len2\}}).
1311         *
1312         * @return The TM-Score
1313         * @throws StructureException
1314         */
1315        public static double getTMScore(Atom[] atomSet1, Atom[] atomSet2, int len1,
1316                        int len2, boolean normalizeMin) throws StructureException {
1317                if (atomSet1.length != atomSet2.length) {
1318                        throw new StructureException(
1319                                        "The two atom sets are not of same length!");
1320                }
1321                if (atomSet1.length > len1) {
1322                        throw new StructureException(
1323                                        "len1 must be greater or equal to the alignment length!");
1324                }
1325                if (atomSet2.length > len2) {
1326                        throw new StructureException(
1327                                        "len2 must be greater or equal to the alignment length!");
1328                }
1329
1330                int Lnorm;
1331                if (normalizeMin) {
1332                        Lnorm = Math.min(len1, len2);
1333                } else {
1334                        Lnorm = Math.max(len1, len2);
1335                }
1336
1337                int Laln = atomSet1.length;
1338
1339                double d0 = 1.24 * Math.cbrt(Lnorm - 15.) - 1.8;
1340                double d0sq = d0 * d0;
1341
1342                double sum = 0;
1343                for (int i = 0; i < Laln; i++) {
1344                        double d = Calc.getDistance(atomSet1[i], atomSet2[i]);
1345                        sum += 1. / (1 + d * d / d0sq);
1346                }
1347
1348                return sum / Lnorm;
1349        }
1350}