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}