org.biojava.nbio.structure.xtal.CrystalCell Maven / Gradle / Ivy
Show all versions of biojava-structure Show documentation
/*
* BioJava development code
*
* This code may be freely distributed and modified under the
* terms of the GNU Lesser General Public Licence. This should
* be distributed with the code. If you do not have a copy,
* see:
*
* http://www.gnu.org/copyleft/lesser.html
*
* Copyright for this code is held jointly by the individual
* authors. These should be listed in @author doc comments.
*
* For more information on the BioJava project and its aims,
* or to join the biojava-l mailing list, visit the home page
* at:
*
* http://www.biojava.org/
*
*/
package org.biojava.nbio.structure.xtal;
import javax.vecmath.*;
import java.io.Serializable;
import java.util.ArrayList;
import java.util.Collections;
import java.util.Locale;
//import org.slf4j.Logger;
//import org.slf4j.LoggerFactory;
/**
* A crystal cell's parameters.
*
* @author duarte_j
*
*/
public class CrystalCell implements Serializable {
private static final long serialVersionUID = 1L;
//private static final Logger logger = LoggerFactory.getLogger(CrystalCell.class);
public static final double MIN_VALID_CELL_SIZE = 10.0; // the minimum admitted for a crystal cell
private double a;
private double b;
private double c;
private double alpha;
private double beta;
private double gamma;
private double alphaRad;
private double betaRad;
private double gammaRad;
private double volume; // cached volume
private double maxDimension; // cached max dimension
private Matrix3d M; // cached basis change transformation matrix
private Matrix3d Minv; // cached basis change transformation matrix
private Matrix3d Mtransp;
private Matrix3d MtranspInv;
public CrystalCell() {
}
public CrystalCell(double a, double b, double c, double alpha, double beta, double gamma){
this.a = a;
this.b = b;
this.c = c;
this.setAlpha(alpha); // initialises both alpha and alphaRad
this.setBeta(beta);
this.setGamma(gamma);
}
public double getA() {
return a;
}
public void setA(double a) {
this.a = a;
}
public double getB() {
return b;
}
public void setB(double b) {
this.b = b;
}
public double getC() {
return c;
}
public void setC(double c) {
this.c = c;
}
public double getAlpha() {
return alpha;
}
public void setAlpha(double alpha) {
this.alpha = alpha;
this.alphaRad = Math.toRadians(alpha);
}
public double getBeta() {
return beta;
}
public void setBeta(double beta) {
this.beta = beta;
this.betaRad = Math.toRadians(beta);
}
public double getGamma() {
return gamma;
}
public void setGamma(double gamma) {
this.gamma = gamma;
this.gammaRad = Math.toRadians(gamma);
}
/**
* Returns the volume of this unit cell.
* See http://en.wikipedia.org/wiki/Parallelepiped
* @return
*/
public double getVolume() {
if (volume!=0) {
return volume;
}
volume = a*b*c*
Math.sqrt(1-Math.cos(alphaRad)*Math.cos(alphaRad)-Math.cos(betaRad)*Math.cos(betaRad)-Math.cos(gammaRad)*Math.cos(gammaRad)
+2.0*Math.cos(alphaRad)*Math.cos(betaRad)*Math.cos(gammaRad));
return volume;
}
/**
* Get the index of a unit cell to which the query point belongs.
*
* For instance, all points in the unit cell at the origin will return (0,0,0);
* Points in the unit cell one unit further along the `a` axis will return (1,0,0),
* etc.
* @param pt Input point (in orthonormal coordinates)
* @return A new point with the three indices of the cell containing pt
*/
public Point3i getCellIndices(Tuple3d pt) {
Point3d p = new Point3d(pt);
this.transfToCrystal(p);
int x = (int)Math.floor(p.x);
int y = (int)Math.floor(p.y);
int z = (int)Math.floor(p.z);
return new Point3i(x,y,z);
}
/**
* Converts the coordinates in pt so that they occur within the (0,0,0) unit cell
* @param pt
*/
public void transfToOriginCell(Tuple3d pt) {
transfToCrystal(pt);
// convert all coordinates to [0,1) interval
pt.x = pt.x<0 ? (pt.x%1.0 + 1.0)%1.0 : pt.x%1.0;
pt.y = pt.y<0 ? (pt.y%1.0 + 1.0)%1.0 : pt.y%1.0;
pt.z = pt.z<0 ? (pt.z%1.0 + 1.0)%1.0 : pt.z%1.0;
transfToOrthonormal(pt);
}
/**
* Converts a set of points so that the reference point falls in the unit cell.
*
* This is useful to transform a whole chain at once, allowing some of the
* atoms to be outside the unit cell, but forcing the centroid to be within it.
*
* @param points A set of points to transform (in orthonormal coordinates)
* @param reference The reference point, which is unmodified but which would
* be in the unit cell were it to have been transformed. It is safe to
* use a member of the points array here.
*/
public void transfToOriginCell(Tuple3d[] points, Tuple3d reference) {
reference = new Point3d(reference);//clone
transfToCrystal(reference);
int x = (int)Math.floor(reference.x);
int y = (int)Math.floor(reference.y);
int z = (int)Math.floor(reference.z);
for( Tuple3d point: points ) {
transfToCrystal(point);
point.x -= x;
point.y -= y;
point.z -= z;
transfToOrthonormal(point);
}
}
/**
*
* @param ops Set of operations in orthonormal coordinates
* @param reference Reference point, which should be in the unit cell after
* each operation (also in orthonormal coordinates)
* @return A set of orthonormal operators with equivalent rotation to the
* inputs, but with translation such that the reference point would fall
* within the unit cell
*/
public Matrix4d[] transfToOriginCellOrthonormal(Matrix4d[] ops, Tuple3d reference) {
// reference in crystal coords
Point3d refXtal = new Point3d(reference);
transfToCrystal(refXtal);
// Convert to crystal coords
Matrix4d[] opsXtal = new Matrix4d[ops.length];
for(int i=0;i vertDists = new ArrayList<>();
vertDists.add(vert0.distance(vert7));
vertDists.add(vert3.distance(vert4));
vertDists.add(vert1.distance(vert6));
vertDists.add(vert2.distance(vert5));
maxDimension = Collections.max(vertDists);
return maxDimension;
}
/**
* Given a scale matrix parsed from the PDB entry (SCALE1,2,3 records),
* checks that the matrix is a consistent scale matrix by comparing the
* cell volume to the inverse of the scale matrix determinant (tolerance of 1/100).
* If they don't match false is returned.
* See the PDB documentation for the SCALE record.
* See also last equation of section 2.5 of "Fundamentals of Crystallography" C. Giacovazzo
* @param scaleMatrix
* @return
*/
public boolean checkScaleMatrixConsistency(Matrix4d scaleMatrix) {
double vol = getVolume();
Matrix3d m = new Matrix3d();
scaleMatrix.getRotationScale(m);
// note we need to have a relaxed tolerance here as the PDB scale matrix is given with not such high precision
// plus we don't want to have false positives, so we stay conservative
double tolerance = vol/100.0;
if ((Math.abs(vol - 1.0/m.determinant() )>tolerance)) {
//System.err.println("Warning! SCALE matrix from PDB does not match 1/determinat == cell volume: "+
// String.format("vol=%6.3f 1/det=%6.3f",vol,1.0/m.determinant()));
return false;
}
// this would be to check our own matrix, must always match!
//if (!deltaComp(vol,1.0/getMTranspose().determinant())) {
// System.err.println("Our calculated SCALE matrix does not match 1/det=cell volume");
//}
return true;
}
/**
* Given a scale matrix parsed from a PDB entry (SCALE1,2,3 records),
* compares it to our calculated Mtranspose matrix to see if they coincide and
* returns true if they do.
* If they don't that means that the PDB entry is not in the standard
* orthogonalisation (NCODE=1 in ccp4).
* In 2011's remediation only 148 PDB entries were found not to be in
* a non-standard orthogonalisation. See:
* http://www.wwpdb.org/documentation/2011remediation_overview-061711.pdf
* For normal cases the scale matrix is diagonal without a translation component.
* Additionally the translation component of the SCALE matrix is also checked to
* make sure it is (0,0,0), if not false is return
* @param scaleMatrix
* @return
*/
public boolean checkScaleMatrix(Matrix4d scaleMatrix) {
Matrix3d mtranspose = getMTranspose();
for (int i=0;i<3;i++) {
for (int j=0;j<3;j++) {
if (!deltaComp(mtranspose.getElement(i, j),scaleMatrix.getElement(i, j))) {
//System.out.println("Our value ("+i+","+j+"): "+getM().getElement(i,j));
//System.out.println("Their value ("+i+","+j+"): "+scaleMatrix.getElement(i,j));
return false;
}
}
}
for (int i=0;i<3;i++) {
if (!deltaComp(scaleMatrix.getElement(i, 3),0)) {
return false;
}
}
return true;
}
private boolean deltaComp(double d1, double d2) {
if (Math.abs(d1-d2)<0.0001) return true;
return false;
}
/**
* Checks whether the dimensions of this crystal cell are reasonable for protein
* crystallography: if all 3 dimensions are below {@value #MIN_VALID_CELL_SIZE} the cell
* is considered unrealistic and false returned
* @return
*/
public boolean isCellReasonable() {
// this check is necessary mostly when reading PDB files that can contain the default 1 1 1 crystal cell
// if we use that further things can go wrong, for instance for interface calculation
// For instance programs like coot produce by default a 1 1 1 cell
if (this.getA()