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.*;
024import java.io.Serializable;
025import java.util.ArrayList;
026import java.util.Collections;
027
028//import org.slf4j.Logger;
029//import org.slf4j.LoggerFactory;
030
031/**
032 * A crystal cell's parameters.
033 *
034 * @author duarte_j
035 *
036 */
037public class CrystalCell implements Serializable {
038
039        private static final long serialVersionUID = 1L;
040        //private static final Logger logger = LoggerFactory.getLogger(CrystalCell.class);
041
042        public static final double MIN_VALID_CELL_SIZE = 10.0; // the minimum admitted for a crystal cell
043
044
045        private double a;
046        private double b;
047        private double c;
048
049        private double alpha;
050        private double beta;
051        private double gamma;
052
053        private double alphaRad;
054        private double betaRad;
055        private double gammaRad;
056
057        private double volume; // cached volume
058
059        private double maxDimension; // cached max dimension
060
061        private Matrix3d M;     // cached basis change transformation matrix
062        private Matrix3d Minv;  // cached basis change transformation matrix
063        private Matrix3d Mtransp;
064        private Matrix3d MtranspInv;
065
066        public CrystalCell() {
067
068        }
069
070        public CrystalCell(double a, double b, double c, double alpha, double beta, double gamma){
071                this.a = a;
072                this.b = b;
073                this.c = c;
074                this.setAlpha(alpha); // initialises both alpha and alphaRad
075                this.setBeta(beta);
076                this.setGamma(gamma);
077        }
078
079        public double getA() {
080                return a;
081        }
082
083        public void setA(double a) {
084                this.a = a;
085        }
086
087        public double getB() {
088                return b;
089        }
090
091        public void setB(double b) {
092                this.b = b;
093        }
094
095        public double getC() {
096                return c;
097        }
098
099        public void setC(double c) {
100                this.c = c;
101        }
102
103        public double getAlpha() {
104                return alpha;
105        }
106
107        public void setAlpha(double alpha) {
108                this.alpha = alpha;
109                this.alphaRad = Math.toRadians(alpha);
110        }
111
112        public double getBeta() {
113                return beta;
114        }
115
116        public void setBeta(double beta) {
117                this.beta = beta;
118                this.betaRad = Math.toRadians(beta);
119        }
120
121        public double getGamma() {
122                return gamma;
123        }
124
125        public void setGamma(double gamma) {
126                this.gamma = gamma;
127                this.gammaRad = Math.toRadians(gamma);
128        }
129
130        /**
131         * Returns the volume of this unit cell.
132         * See http://en.wikipedia.org/wiki/Parallelepiped
133         * @return
134         */
135        public double getVolume() {
136                if (volume!=0) {
137                        return volume;
138                }
139                volume =  a*b*c*
140                Math.sqrt(1-Math.cos(alphaRad)*Math.cos(alphaRad)-Math.cos(betaRad)*Math.cos(betaRad)-Math.cos(gammaRad)*Math.cos(gammaRad)
141                                +2.0*Math.cos(alphaRad)*Math.cos(betaRad)*Math.cos(gammaRad));
142
143                return volume;
144        }
145
146        /**
147         * Get the index of a unit cell to which the query point belongs.
148         *
149         * <p>For instance, all points in the unit cell at the origin will return (0,0,0);
150         * Points in the unit cell one unit further along the `a` axis will return (1,0,0),
151         * etc.
152         * @param pt Input point (in orthonormal coordinates)
153         * @return A new point with the three indices of the cell containing pt
154         */
155        public Point3i getCellIndices(Tuple3d pt) {
156                Point3d p = new Point3d(pt);
157                this.transfToCrystal(p);
158
159                int x = (int)Math.floor(p.x);
160                int y = (int)Math.floor(p.y);
161                int z = (int)Math.floor(p.z);
162                return new Point3i(x,y,z);
163        }
164
165        /**
166         * Converts the coordinates in pt so that they occur within the (0,0,0) unit cell
167         * @param pt
168         */
169        public void transfToOriginCell(Tuple3d pt) {
170                transfToCrystal(pt);
171
172                // convert all coordinates to [0,1) interval
173                pt.x = pt.x<0 ? (pt.x%1.0 + 1.0)%1.0 : pt.x%1.0;
174                pt.y = pt.y<0 ? (pt.y%1.0 + 1.0)%1.0 : pt.y%1.0;
175                pt.z = pt.z<0 ? (pt.z%1.0 + 1.0)%1.0 : pt.z%1.0;
176
177                transfToOrthonormal(pt);
178        }
179
180        /**
181         * Converts a set of points so that the reference point falls in the unit cell.
182         *
183         * This is useful to transform a whole chain at once, allowing some of the
184         * atoms to be outside the unit cell, but forcing the centroid to be within it.
185         *
186         * @param points A set of points to transform (in orthonormal coordinates)
187         * @param reference The reference point, which is unmodified but which would
188         *    be in the unit cell were it to have been transformed. It is safe to
189         *    use a member of the points array here.
190         */
191        public void transfToOriginCell(Tuple3d[] points, Tuple3d reference) {
192                reference = new Point3d(reference);//clone
193                transfToCrystal(reference);
194
195                int x = (int)Math.floor(reference.x);
196                int y = (int)Math.floor(reference.y);
197                int z = (int)Math.floor(reference.z);
198
199                for( Tuple3d point: points ) {
200                        transfToCrystal(point);
201                        point.x -= x;
202                        point.y -= y;
203                        point.z -= z;
204                        transfToOrthonormal(point);
205                }
206        }
207
208        /**
209         *
210         * @param ops Set of operations in orthonormal coordinates
211         * @param reference Reference point, which should be in the unit cell after
212         *  each operation (also in orthonormal coordinates)
213         * @return A set of orthonormal operators with equivalent rotation to the
214         *  inputs, but with translation such that the reference point would fall
215         *  within the unit cell
216         */
217        public Matrix4d[] transfToOriginCellOrthonormal(Matrix4d[] ops, Tuple3d reference) {
218                // reference in crystal coords
219                Point3d refXtal = new Point3d(reference);
220                transfToCrystal(refXtal);
221
222                // Convert to crystal coords
223                Matrix4d[] opsXtal = new Matrix4d[ops.length];
224                for(int i=0;i<ops.length;i++) {
225                        opsXtal[i] = transfToCrystal(ops[i]);
226                }
227
228                Matrix4d[] transformed = transfToOriginCellCrystal(opsXtal, refXtal);
229
230                // Convert back to orthonormal
231                for(int i=0;i<ops.length;i++) {
232                        transformed[i] = transfToOrthonormal(transformed[i]);
233                }
234                return transformed;
235        }
236        /**
237         *
238         * @param ops Set of operations in crystal coordinates
239         * @param reference Reference point, which should be in the unit cell after
240         *  each operation (also in crystal coordinates)
241         * @return A set of crystal operators with equivalent rotation to the
242         *  inputs, but with translation such that the reference point would fall
243         *  within the unit cell
244         */
245        public Matrix4d[] transfToOriginCellCrystal(Matrix4d[] ops, Tuple3d reference) {
246                Matrix4d[] transformed = new Matrix4d[ops.length];
247                for(int j=0;j<ops.length;j++) {
248                        Matrix4d op = ops[j];
249
250                        // transform the reference point
251                        Point3d xXtal = new Point3d(reference);
252                        op.transform(xXtal);
253
254                        // Calculate unit cell of transformed reference
255                        int x = (int)Math.floor(xXtal.x);
256                        int y = (int)Math.floor(xXtal.y);
257                        int z = (int)Math.floor(xXtal.z);
258
259                        Matrix4d translation = new Matrix4d();
260                        translation.set(new Vector3d(-x,-y,-z));
261
262                        // Compose op with an additional translation operator
263                        translation.mul(op);
264
265                        Point3d ref2 = new Point3d(reference);
266                        translation.transform(ref2);
267
268                        transformed[j] = translation;
269                }
270
271                return transformed;
272        }
273
274        /**
275         * Transform given Matrix4d in crystal basis to the orthonormal basis using
276         * the PDB axes convention (NCODE=1)
277         * @param m
278         * @return
279         */
280        public Matrix4d transfToOrthonormal(Matrix4d m) {
281                Vector3d trans = new Vector3d(m.m03,m.m13,m.m23);
282                transfToOrthonormal(trans);
283
284                Matrix3d rot = new Matrix3d();
285                m.getRotationScale(rot);
286                // see Giacovazzo section 2.E, eq. 2.E.1
287                // Rprime = MT-1 * R * MT
288                rot.mul(this.getMTranspose());
289                rot.mul(this.getMTransposeInv(),rot);
290
291                return new Matrix4d(rot,trans,1.0);
292        }
293
294        /**
295         * Transforms the given crystal basis coordinates into orthonormal coordinates.
296         * e.g. transfToOrthonormal(new Point3d(1,1,1)) returns the orthonormal coordinates of the
297         * vertex of the unit cell.
298         * See Giacovazzo section 2.E, eq. 2.E.1 (or any linear algebra manual)
299         * @param v
300         */
301        public void transfToOrthonormal(Tuple3d v) {
302                // see Giacovazzo section 2.E, eq. 2.E.1
303                getMTransposeInv().transform(v);
304        }
305
306        /**
307         * Transform given Matrix4d in orthonormal basis to the crystal basis using
308         * the PDB axes convention (NCODE=1)
309         * @param m
310         * @return
311         */
312        public Matrix4d transfToCrystal(Matrix4d m) {
313                Vector3d trans = new Vector3d(m.m03,m.m13,m.m23);
314                transfToCrystal(trans);
315
316                Matrix3d rot = new Matrix3d();
317                m.getRotationScale(rot);
318                // see Giacovazzo section 2.E, eq. 2.E.1 (the inverse equation)
319                // R = MT * Rprime * MT-1
320                rot.mul(this.getMTransposeInv());
321                rot.mul(this.getMTranspose(),rot);
322
323                return new Matrix4d(rot,trans,1.0);
324        }
325
326        /**
327         * Transforms the given orthonormal basis coordinates into crystal coordinates.
328         * See Giacovazzo eq 2.20 (or any linear algebra manual)
329         * @param v
330         */
331        public void transfToCrystal(Tuple3d v) {
332                getMTranspose().transform(v);
333        }
334
335        /**
336         * Returns the change of basis (crystal to orthonormal) transform matrix, that is
337         * M inverse in the notation of Giacovazzo.
338         * Using the PDB axes convention
339         * (CCP4 uses NCODE identifiers to distinguish the different conventions, the PDB one is called NCODE=1)
340         * The matrix is only calculated upon first call of this method, thereafter it is cached.
341         * See "Fundamentals of Crystallography" C. Giacovazzo, section 2.5 (eq 2.30)
342         *
343         * The non-standard orthogonalisation codes (NCODE for ccp4) are flagged in REMARK 285 after 2011's remediation
344         * with text: "THE ENTRY COORDINATES ARE NOT PRESENTED IN THE STANDARD CRYSTAL FRAME". There were only 148 PDB
345         * entries with non-standard code in 2011. See:
346         * http://www.wwpdb.org/documentation/2011remediation_overview-061711.pdf
347         * The SCALE1,2,3 records contain the correct transformation matrix (what Giacovazzo calls M matrix).
348         * In those cases if we calculate the M matrix following Giacovazzo's equations here, we get an entirely wrong one.
349         * Examples of PDB with non-standard orthogonalisation are 1bab and 1bbb.
350         * @return
351         */
352        private Matrix3d getMInv() {
353                if (Minv!=null) {
354                        return Minv;
355                }
356
357                // see eq. 2.30 Giacovazzo
358                Minv =  new Matrix3d(                    this.a,                                           0,              0,
359                                                          this.b*Math.cos(gammaRad),                   this.b*Math.sin(gammaRad),              0,
360                                                          this.c*Math.cos(betaRad) , -this.c*Math.sin(betaRad)*getCosAlphaStar(),  1.0/getCstar());
361                return Minv;
362        }
363
364        // another axes convention (from Giacovazzo) eq. 2.31b, not sure what would be the NCODE of this
365//      private Matrix3d getMInv() {
366//              if (Minv!=null) {
367//                      return Minv;
368//              }
369//              Minv = new Matrix3d( 1/getAstar(),  -(getCosGammaStar()/getSinGammaStar())/getAstar(), this.a*Math.cos(betaRad),
370//                                                      0,                   1/(getBstar()*getSinGammaStar()), this.b*Math.sin(alphaRad),
371//                                                      0,                                                  0,                  this.c);
372//              return Minv;
373//      }
374
375        // and yet another axes convention (from Giacovazzo) eq. 2.31c, not sure what would be the NCODE of this
376//      private Matrix3d getMInv() {
377//              if (Minv!=null) {
378//                      return Minv;
379//              }
380//              Minv = new Matrix3d( alphaRad*Math.sin(gammaRad)*getSinBetaStar(), this.a*Math.cos(gammaRad), this.a*Math.sin(gammaRad)*getCosBetaStar(),
381//                                                                                      0,                    this.b,                                          0,
382//                                                                                      0, this.c*Math.cos(alphaRad),                   this.c*Math.sin(alphaRad));
383//              return Minv;
384//      }
385
386        // relationships among direct and reciprocal lattice parameters
387        // see Table 2.1 of chapter 2 of Giacovazzo
388        @SuppressWarnings("unused")
389        private double getAstar() {
390                return (this.b*this.c*Math.sin(alphaRad))/getVolume();
391        }
392
393        @SuppressWarnings("unused")
394        private double getBstar() {
395                return (this.a*this.c*Math.sin(betaRad))/getVolume();
396        }
397
398        private double getCstar() {
399                return (this.a*this.b*Math.sin(gammaRad))/getVolume();
400        }
401
402        private double getCosAlphaStar() {
403                return (Math.cos(betaRad)*Math.cos(gammaRad)-Math.cos(alphaRad))/(Math.sin(betaRad)*Math.sin(gammaRad));
404        }
405
406        @SuppressWarnings("unused")
407        private double getCosBetaStar() {
408                return (Math.cos(alphaRad)*Math.cos(gammaRad)-Math.cos(betaRad))/(Math.sin(alphaRad)*Math.sin(gammaRad));
409        }
410
411        @SuppressWarnings("unused")
412        private double getCosGammaStar() {
413                return (Math.cos(alphaRad)*Math.cos(betaRad)-Math.cos(gammaRad))/(Math.sin(alphaRad)*Math.sin(betaRad));
414        }
415
416        @SuppressWarnings("unused")
417        private double getSinAlphaStar() {
418                return getVolume()/(this.a*this.b*this.c*Math.sin(betaRad)*Math.sin(gammaRad));
419        }
420
421        @SuppressWarnings("unused")
422        private double getSinBetaStar() {
423                return getVolume()/(this.a*this.b*this.c*Math.sin(alphaRad)*Math.sin(gammaRad));
424        }
425
426        @SuppressWarnings("unused")
427        private double getSinGammaStar() {
428                return getVolume()/(this.a*this.b*this.c*Math.sin(alphaRad)*Math.sin(betaRad));
429        }
430
431        /**
432         * Returns the change of basis (orthonormal to crystal) transform matrix, that is
433         * M in the notation of Giacovazzo.
434         * Using the PDB convention (NCODE=1).
435         * The matrix is only calculated upon first call of this method, thereafter it is cached.
436         * See "Fundamentals of Crystallography" C. Giacovazzo, section 2.5
437         * @return
438         */
439        private Matrix3d getM() {
440                if (M!=null){
441                        return M;
442                }
443                M = new Matrix3d();
444                M.invert(getMInv());
445                return M;
446        }
447
448        public Matrix3d getMTranspose() {
449                if (Mtransp!=null){
450                        return Mtransp;
451                }
452                Matrix3d M = getM();
453                Mtransp = new Matrix3d();
454                Mtransp.transpose(M);
455                return Mtransp;
456        }
457
458        private Matrix3d getMTransposeInv() {
459                if (MtranspInv!=null){
460                        return MtranspInv;
461                }
462                MtranspInv = new Matrix3d();
463                MtranspInv.invert(getMTranspose());
464                return MtranspInv;
465        }
466
467        /**
468         * Gets the maximum dimension of the unit cell.
469         * @return
470         */
471        public double getMaxDimension() {
472                if (maxDimension!=0) {
473                        return maxDimension;
474                }
475                Point3d vert0 = new Point3d(0,0,0);
476                Point3d vert1 = new Point3d(1,0,0);
477                transfToOrthonormal(vert1);
478                Point3d vert2 = new Point3d(0,1,0);
479                transfToOrthonormal(vert2);
480                Point3d vert3 = new Point3d(0,0,1);
481                transfToOrthonormal(vert3);
482                Point3d vert4 = new Point3d(1,1,0);
483                transfToOrthonormal(vert4);
484                Point3d vert5 = new Point3d(1,0,1);
485                transfToOrthonormal(vert5);
486                Point3d vert6 = new Point3d(0,1,1);
487                transfToOrthonormal(vert6);
488                Point3d vert7 = new Point3d(1,1,1);
489                transfToOrthonormal(vert7);
490
491                ArrayList<Double> vertDists = new ArrayList<Double>();
492                vertDists.add(vert0.distance(vert7));
493                vertDists.add(vert3.distance(vert4));
494                vertDists.add(vert1.distance(vert6));
495                vertDists.add(vert2.distance(vert5));
496                maxDimension = Collections.max(vertDists);
497                return maxDimension;
498        }
499
500        /**
501         * Given a scale matrix parsed from the PDB entry (SCALE1,2,3 records),
502         * checks that the matrix is a consistent scale matrix by comparing the
503         * cell volume to the inverse of the scale matrix determinant (tolerance of 1/100).
504         * If they don't match false is returned.
505         * See the PDB documentation for the SCALE record.
506         * See also last equation of section 2.5 of "Fundamentals of Crystallography" C. Giacovazzo
507         * @param scaleMatrix
508         * @return
509         */
510        public boolean checkScaleMatrixConsistency(Matrix4d scaleMatrix) {
511
512                double vol = getVolume();
513                Matrix3d m = new Matrix3d();
514                scaleMatrix.getRotationScale(m);
515
516                // note we need to have a relaxed tolerance here as the PDB scale matrix is given with not such high precision
517                // plus we don't want to have false positives, so we stay conservative
518                double tolerance = vol/100.0;
519                if ((Math.abs(vol - 1.0/m.determinant() )>tolerance)) {
520                        //System.err.println("Warning! SCALE matrix from PDB does not match 1/determinat == cell volume: "+
521                        //              String.format("vol=%6.3f  1/det=%6.3f",vol,1.0/m.determinant()));
522                        return false;
523                }
524                // this would be to check our own matrix, must always match!
525                //if (!deltaComp(vol,1.0/getMTranspose().determinant())) {
526                //      System.err.println("Our calculated SCALE matrix does not match 1/det=cell volume");
527                //}
528
529                return true;
530
531        }
532
533        /**
534         * Given a scale matrix parsed from a PDB entry (SCALE1,2,3 records),
535         * compares it to our calculated Mtranspose matrix to see if they coincide and
536         * returns true if they do.
537         * If they don't that means that the PDB entry is not in the standard
538         * orthogonalisation (NCODE=1 in ccp4).
539         * In 2011's remediation only 148 PDB entries were found not to be in
540         * a non-standard orthogonalisation. See:
541         * http://www.wwpdb.org/documentation/2011remediation_overview-061711.pdf
542         * For normal cases the scale matrix is diagonal without a translation component.
543         * Additionally the translation component of the SCALE matrix is also checked to
544         * make sure it is (0,0,0), if not false is return
545         * @param scaleMatrix
546         * @return
547         */
548        public boolean checkScaleMatrix(Matrix4d scaleMatrix) {
549
550                Matrix3d mtranspose = getMTranspose();
551                for (int i=0;i<3;i++) {
552                        for (int j=0;j<3;j++) {
553                                if (!deltaComp(mtranspose.getElement(i, j),scaleMatrix.getElement(i, j))) {
554                                        //System.out.println("Our value   ("+i+","+j+"): "+getM().getElement(i,j));
555                                        //System.out.println("Their value ("+i+","+j+"): "+scaleMatrix.getElement(i,j));
556                                        return false;
557                                }
558                        }
559                }
560                for (int i=0;i<3;i++) {
561                        if (!deltaComp(scaleMatrix.getElement(i, 3),0)) {
562                                return false;
563                        }
564                }
565                return true;
566        }
567
568        private boolean deltaComp(double d1, double d2) {
569                if (Math.abs(d1-d2)<0.0001) return true;
570                return false;
571        }
572
573        /**
574         * Checks whether the dimensions of this crystal cell are reasonable for protein
575         * crystallography: if all 3 dimensions are below {@value #MIN_VALID_CELL_SIZE} the cell
576         * is considered unrealistic and false returned
577         * @return
578         */
579        public boolean isCellReasonable() {
580                // this check is necessary mostly when reading PDB files that can contain the default 1 1 1 crystal cell
581                // if we use that further things can go wrong, for instance for interface calculation
582                // For instance programs like coot produce by default a 1 1 1 cell
583
584                if (this.getA()<MIN_VALID_CELL_SIZE &&
585                                this.getB()<MIN_VALID_CELL_SIZE &&
586                                this.getC()<MIN_VALID_CELL_SIZE) {
587                        return false;
588                }
589
590                return true;
591
592        }
593
594        @Override
595        public String toString() {
596                return String.format("a%7.2f b%7.2f c%7.2f alpha%6.2f beta%6.2f gamma%6.2f", a, b, c, alpha, beta, gamma);
597        }
598}