Source code

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