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 */
021package org.biojava.nbio.structure.xtal;
022
023import javax.vecmath.Matrix3d;
024import javax.vecmath.Matrix4d;
025import javax.vecmath.Point3i;
026import javax.vecmath.Vector3d;
027import java.io.Serializable;
028
029import static java.lang.Math.abs;
030
031
032/**
033 * Representation of a transformation in a crystal:
034 * - a transformation id (each of the transformations in a space group, 0 to m)
035 * - a crystal translation
036 * The transformation matrix in crystal basis is stored, representing the basic
037 * transformation together with the crystal translation.
038 * Contains methods to check for equivalent transformations.
039 *
040 *
041 * @author duarte_j
042 *
043 */
044public class CrystalTransform implements Serializable {
045
046        private static final long serialVersionUID = 1L;
047
048
049        public static final Matrix4d IDENTITY = new Matrix4d(1,0,0,0, 0,1,0,0, 0,0,1,0, 0,0,0,1);
050
051        /**
052         * The space group to which this transform belongs
053         */
054        private final SpaceGroup sg;
055
056        /**
057         * The transform id corresponding to the SpaceGroup's transform indices.
058         * From 0 (identity) to m (m=number of symmetry operations of the space group)
059         * It is unique within the unit cell but equivalent units of different crystal unit cells
060         * will have same id
061         */
062        private int transformId;
063
064        /**
065         * The 4-dimensional matrix transformation in crystal basis.
066         * Note that the translational component of this matrix is not necessarily
067         * identical to crystalTranslation since some operators have fractional
068         * translations within the cell
069         */
070        private Matrix4d matTransform;
071
072        /**
073         * The crystal translation (always integer)
074         */
075        private Point3i crystalTranslation;
076
077
078        /**
079         * Creates a new CrystalTransform representing the identity transform
080         * in cell (0,0,0)
081         */
082        public CrystalTransform(SpaceGroup sg) {
083                this.sg = sg;
084                this.transformId = 0;
085                this.matTransform = (Matrix4d)IDENTITY.clone();
086                this.crystalTranslation = new Point3i(0,0,0);
087        }
088
089        /**
090         * Represents the n-th transform
091         * @param sg
092         * @param transformId
093         */
094        public CrystalTransform(SpaceGroup sg, int transformId) {
095                this.sg = sg;
096                this.transformId = transformId;
097                if (sg==null && transformId==0) {
098                        this.matTransform = (Matrix4d)IDENTITY.clone();
099                } else if (sg==null) {
100                        throw new IllegalArgumentException("Space Group cannot be null if transformId!=0");
101                } else {
102                        this.matTransform = (Matrix4d)sg.getTransformation(transformId).clone();
103                }
104                this.crystalTranslation = new Point3i(0,0,0);
105        }
106
107        /**
108         * Copy constructor
109         * @param transform
110         */
111        public CrystalTransform(CrystalTransform transform) {
112                this.sg = transform.sg;
113                this.transformId = transform.transformId;
114                this.matTransform = new Matrix4d(transform.matTransform);
115                this.crystalTranslation = new Point3i(transform.crystalTranslation);
116        }
117
118        public Matrix4d getMatTransform() {
119                return matTransform;
120        }
121
122        public void setMatTransform(Matrix4d matTransform) {
123                this.matTransform = matTransform;
124        }
125
126        public Point3i getCrystalTranslation() {
127                return crystalTranslation;
128        }
129
130        public void translate(Point3i translation) {
131                matTransform.m03 = matTransform.m03+translation.x;
132                matTransform.m13 = matTransform.m13+translation.y;
133                matTransform.m23 = matTransform.m23+translation.z;
134
135                crystalTranslation.add(translation);
136
137        }
138
139        /**
140         * Returns true if the given CrystalTransform is equivalent to this one.
141         * Two crystal transforms are equivalent if one is the inverse of the other, i.e.
142         * their transformation matrices multiplication is equal to the identity.
143         * @param other
144         * @return
145         */
146        public boolean isEquivalent(CrystalTransform other) {
147                Matrix4d mul = new Matrix4d();
148                mul.mul(this.matTransform,other.matTransform);
149
150                if (mul.epsilonEquals(IDENTITY, 0.0001)) {
151                        return true;
152                }
153                return false;
154        }
155
156        /**
157         * Tells whether this transformation is a pure crystal lattice translation,
158         * i.e. no rotational component and an integer translation vector.
159         * @return
160         */
161        public boolean isPureCrystalTranslation() {
162                return (transformId==0 && (crystalTranslation.x!=0 || crystalTranslation.y!=0 || crystalTranslation.z!=0));
163        }
164
165        /**
166         * Tells whether this transformation is the identity: no rotation and no translation
167         * @return
168         */
169        public boolean isIdentity() {
170                return (transformId==0 && crystalTranslation.x==0 && crystalTranslation.y==0 && crystalTranslation.z==0);
171        }
172
173        /**
174         * Tells whether this transformation is a pure translation:
175         * either a pure crystal (lattice) translation or a fractional (within
176         * unit cell) translation: space groups Ixxx, Cxxx, Fxxx have operators
177         * with fractional translations within the unit cell.
178         * @return
179         */
180        public boolean isPureTranslation() {
181                if (isPureCrystalTranslation()) return true;
182                if (SpaceGroup.deltaComp(matTransform.m00,1,SpaceGroup.DELTA) &&
183                        SpaceGroup.deltaComp(matTransform.m01,0,SpaceGroup.DELTA) &&
184                        SpaceGroup.deltaComp(matTransform.m02,0,SpaceGroup.DELTA) &&
185
186                        SpaceGroup.deltaComp(matTransform.m10,0,SpaceGroup.DELTA) &&
187                        SpaceGroup.deltaComp(matTransform.m11,1,SpaceGroup.DELTA) &&
188                        SpaceGroup.deltaComp(matTransform.m12,0,SpaceGroup.DELTA) &&
189
190                        SpaceGroup.deltaComp(matTransform.m20,0,SpaceGroup.DELTA) &&
191                        SpaceGroup.deltaComp(matTransform.m21,0,SpaceGroup.DELTA) &&
192                        SpaceGroup.deltaComp(matTransform.m22,1,SpaceGroup.DELTA) &&
193                        (       Math.abs(matTransform.m03-0.0)>SpaceGroup.DELTA ||
194                                Math.abs(matTransform.m13-0.0)>SpaceGroup.DELTA ||
195                                Math.abs(matTransform.m23-0.0)>SpaceGroup.DELTA)) {
196                        return true;
197                }
198
199                return false;
200        }
201
202        /**
203         * Tells whether this transformation contains a fractional translational
204         * component (whatever its rotational component). A fractional translation
205         * together with a rotation means a screw axis.
206         * @return
207         */
208        public boolean isFractionalTranslation() {
209                if ((Math.abs(matTransform.m03-crystalTranslation.x)>SpaceGroup.DELTA) ||
210                        (Math.abs(matTransform.m13-crystalTranslation.y)>SpaceGroup.DELTA) ||
211                        (Math.abs(matTransform.m23-crystalTranslation.z)>SpaceGroup.DELTA)) {
212                        return true;
213                }
214                return false;
215        }
216
217        /**
218         * Tells whether this transformation is a rotation disregarding the translational component,
219         * i.e. either pure rotation or screw rotation, but not improper rotation.
220         * @return
221         */
222        public boolean isRotation() {
223                // if no SG, that means a non-crystallographic entry (e.g. NMR). We return false
224                if (sg==null) return false;
225
226                // this also takes care of case <0 (improper rotations): won't be considered as rotations
227                if (sg.getAxisFoldType(this.transformId)>1) return true;
228
229                return false;
230        }
231
232        /**
233         * Returns the TransformType of this transformation: AU, crystal translation, fractional translation
234         * , 2 3 4 6-fold rotations, 2 3 4 6-fold screw rotations, -1 -3 -2 -4 -6 inversions/rotoinversions.
235         * @return
236         */
237        public TransformType getTransformType() {
238
239                // if no SG, that means a non-crystallographic entry (e.g. NMR). We return AU as type
240                if (sg==null) return TransformType.AU;
241
242                int foldType = sg.getAxisFoldType(this.transformId);
243                boolean isScrewOrGlide = false;
244                Vector3d translScrewComponent = getTranslScrewComponent();
245                if (Math.abs(translScrewComponent.x-0.0)>SpaceGroup.DELTA ||
246                        Math.abs(translScrewComponent.y-0.0)>SpaceGroup.DELTA ||
247                        Math.abs(translScrewComponent.z-0.0)>SpaceGroup.DELTA) {
248
249                        isScrewOrGlide = true;
250                }
251
252                if (foldType>1) {
253
254                        if (isScrewOrGlide) {
255                                switch (foldType) {
256                                case 2:
257                                        return TransformType.TWOFOLDSCREW;
258                                case 3:
259                                        return TransformType.THREEFOLDSCREW;
260                                case 4:
261                                        return TransformType.FOURFOLDSCREW;
262                                case 6:
263                                        return TransformType.SIXFOLDSCREW;
264                                default:
265                                        throw new NullPointerException("This transformation did not fall into any of the known types! This is most likely a bug.");
266                                }
267                        } else {
268                                switch (foldType) {
269                                case 2:
270                                        return TransformType.TWOFOLD;
271                                case 3:
272                                        return TransformType.THREEFOLD;
273                                case 4:
274                                        return TransformType.FOURFOLD;
275                                case 6:
276                                        return TransformType.SIXFOLD;
277                                default:
278                                        throw new NullPointerException("This transformation did not fall into any of the known types! This is most likely a bug.");
279                                }
280                        }
281
282                } else if (foldType<0) {
283                        switch (foldType) {
284                        case -1:
285                                return TransformType.ONEBAR;
286                        case -2:
287                                if (isScrewOrGlide) {
288                                        return TransformType.GLIDE;
289                                }
290                                return TransformType.TWOBAR;
291                        case -3:
292                                return TransformType.THREEBAR;
293                        case -4:
294                                return TransformType.FOURBAR;
295                        case -6:
296                                return TransformType.SIXBAR;
297                        default:
298                                throw new NullPointerException("This transformation did not fall into any of the known types! This is most likely a bug.");
299                        }
300                } else {
301                        if (isIdentity()) {
302                                return TransformType.AU;
303                        }
304                        if (isPureCrystalTranslation()) {
305                                return TransformType.XTALTRANSL;
306                        }
307                        if (isFractionalTranslation()) {
308                                return TransformType.CELLTRANSL;
309                        }
310                        throw new NullPointerException("This transformation did not fall into any of the known types! This is most likely a bug.");
311                }
312
313        }
314
315        public Vector3d getTranslScrewComponent() {
316
317                return getTranslScrewComponent(matTransform);
318
319        }
320
321        public int getTransformId() {
322                return transformId;
323        }
324
325        public void setTransformId(int transformId) {
326                this.transformId = transformId;
327        }
328
329        @Override
330        public String toString() {
331                return String.format("[%2d-(%s)]",transformId,toXYZString());
332        }
333
334        /**
335         * Expresses this transformation in terms of x,y,z fractional coordinates.
336         *
337         * Examples:
338         * @return
339         */
340        public String toXYZString() {
341                StringBuilder str = new StringBuilder();
342
343                for(int i=0;i<3;i++) { //for each row
344                        boolean emptyRow = true;
345
346                        double coef; // TODO work with rational numbers
347
348
349                        // X
350                        coef = matTransform.getElement(i, 0);
351
352                        // Three cases for coef: zero, one, non-one
353                        if(abs(coef) > 1e-6 ) { // Non-zero
354                                if( abs( abs(coef)-1 ) < 1e-6 ) { // +/- 1
355                                        if( coef < 0 ) {
356                                                str.append("-");
357                                        }
358                                } else {
359                                        str.append(formatCoef(coef));
360                                        str.append("*");
361                                }
362                                str.append("x");
363                                emptyRow = false;
364                        }
365
366                        // Y
367                        coef = matTransform.getElement(i, 1);
368
369                        if(abs(coef) > 1e-6 ) { // Non-zero
370                                if( abs( abs(coef)-1 ) < 1e-6 ) { // +/- 1
371                                        if( coef < 0 ) {
372                                                str.append("-");
373                                        } else if( !emptyRow ) {
374                                                str.append("+");
375                                        }
376                                } else {
377                                        if( !emptyRow && coef > 0) {
378                                                str.append("+");
379                                        }
380                                        str.append(formatCoef(coef));
381                                        str.append("*");
382                                }
383                                str.append("y");
384                                emptyRow = false;
385                        }
386
387                        // Z
388                        coef = matTransform.getElement(i, 2);
389
390                        if(abs(coef) > 1e-6 ) { // Non-zero
391                                if( abs( abs(coef)-1 ) < 1e-6 ) { // +/- 1
392                                        if( coef < 0 ) {
393                                                str.append("-");
394                                        } else if( !emptyRow ) {
395                                                str.append("+");
396                                        }
397                                } else {
398                                        if( !emptyRow && coef > 0) {
399                                                str.append("+");
400                                        }
401                                        str.append(formatCoef(coef));
402                                        str.append("*");
403                                }
404                                str.append("z");
405                                emptyRow = false;
406                        }
407
408                        // Intercept
409                        coef = matTransform.getElement(i, 3);
410
411                        if(abs(coef) > 1e-6 ) { // Non-zero
412                                if( !emptyRow && coef > 0) {
413                                        str.append("+");
414                                }
415                                str.append(formatCoef(coef));
416                        }
417
418                        if(i<2) {
419                                str.append(",");
420                        }
421                }
422
423
424                return str.toString();
425        }
426        /**
427         * helper function to format simple fractions into rationals
428         * @param coef
429         * @return
430         */
431        private String formatCoef(double coef) {
432                double tol = 1e-6; // rounding tolerance
433
434                // zero case
435                if( Math.abs(coef) < tol) {
436                        return "0";
437                }
438
439                // integer case
440                long num = Math.round(coef);
441                if( Math.abs(num - coef) < tol) {
442                        return Long.toString(num);
443                }
444
445                // Other small cases
446                for(int denom = 2; denom < 12; denom++ ) {
447                        num = Math.round(coef*denom);
448                        if( num - coef*denom < tol ) {
449                                return String.format("%d/%d",num, denom);
450                        }
451                }
452
453                // Give up and use floating point;
454                return String.format("%.3f", coef);
455        }
456
457        /**
458         * Given a transformation matrix containing a rotation and translation returns the
459         * screw component of the rotation.
460         * See http://www.crystallography.fr/mathcryst/pdf/Gargnano/Aroyo_Gargnano_1.pdf
461         * @param m
462         * @return
463         */
464        public static Vector3d getTranslScrewComponent(Matrix4d m) {
465
466                int foldType = SpaceGroup.getRotAxisType(m);
467                // For reference see:
468                // http://www.crystallography.fr/mathcryst/pdf/Gargnano/Aroyo_Gargnano_1.pdf
469
470                Vector3d transl = null;
471
472                Matrix3d W =
473                                new Matrix3d(m.m00,m.m01,m.m02,
474                                                m.m10,m.m11,m.m12,
475                                                m.m20,m.m21,m.m22);
476
477                if (foldType>=0) {
478
479                        // the Y matrix: Y = W^k-1 + W^k-2 ... + W + I  ; with k the fold type
480                        Matrix3d Y = new Matrix3d(1,0,0, 0,1,0, 0,0,1);
481                        Matrix3d Wk = new Matrix3d(1,0,0, 0,1,0, 0,0,1);
482
483                        for (int k=0;k<foldType;k++) {
484                                Wk.mul(W); // k=0 Wk=W, k=1 Wk=W^2, k=2 Wk=W^3, ... k=foldType-1, Wk=W^foldType
485                                if (k!=foldType-1) Y.add(Wk);
486                        }
487
488                        transl = new Vector3d(m.m03, m.m13, m.m23);
489                        Y.transform(transl);
490
491                        transl.scale(1.0/foldType);
492
493                } else {
494
495                        if (foldType==-2) { // there are glide planes only in -2
496                                Matrix3d Y = new Matrix3d(1,0,0, 0,1,0, 0,0,1);
497                                Y.add(W);
498
499                                transl = new Vector3d(m.m03, m.m13, m.m23);
500                                Y.transform(transl);
501
502                                transl.scale(1.0/2.0);
503                        } else { // for -1, -3, -4 and -6 there's nothing to do: fill with 0s
504                                transl = new Vector3d(0,0,0);
505                        }
506                }
507
508                return transl;
509        }
510}