/*
* 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.Matrix3d;
import javax.vecmath.Matrix4d;
import javax.vecmath.Point3i;
import javax.vecmath.Vector3d;
import java.io.Serializable;
import static java.lang.Math.abs;
/**
* Representation of a transformation in a crystal:
* - a transformation id (each of the transformations in a space group, 0 to m)
* - a crystal translation
* The transformation matrix in crystal basis is stored, representing the basic
* transformation together with the crystal translation.
* Contains methods to check for equivalent transformations.
*
*
* @author duarte_j
*
*/
public class CrystalTransform implements Serializable {
private static final long serialVersionUID = 1L;
public static final Matrix4d IDENTITY = new Matrix4d(1,0,0,0, 0,1,0,0, 0,0,1,0, 0,0,0,1);
/**
* The space group to which this transform belongs
*/
private final SpaceGroup sg;
/**
* The transform id corresponding to the SpaceGroup's transform indices.
* From 0 (identity) to m (m=number of symmetry operations of the space group)
* It is unique within the unit cell but equivalent units of different crystal unit cells
* will have same id
*/
private int transformId;
/**
* The 4-dimensional matrix transformation in crystal basis.
* Note that the translational component of this matrix is not necessarily
* identical to crystalTranslation since some operators have fractional
* translations within the cell
*/
private Matrix4d matTransform;
/**
* The crystal translation (always integer)
*/
private Point3i crystalTranslation;
/**
* Creates a new CrystalTransform representing the identity transform
* in cell (0,0,0)
*/
public CrystalTransform(SpaceGroup sg) {
this.sg = sg;
this.transformId = 0;
this.matTransform = (Matrix4d)IDENTITY.clone();
this.crystalTranslation = new Point3i(0,0,0);
}
/**
* Represents the n-th transform
* @param sg
* @param transformId
*/
public CrystalTransform(SpaceGroup sg, int transformId) {
this.sg = sg;
this.transformId = transformId;
if (sg==null && transformId==0) {
this.matTransform = (Matrix4d)IDENTITY.clone();
} else if (sg==null) {
throw new IllegalArgumentException("Space Group cannot be null if transformId!=0");
} else {
this.matTransform = (Matrix4d)sg.getTransformation(transformId).clone();
}
this.crystalTranslation = new Point3i(0,0,0);
}
/**
* Copy constructor
* @param transform
*/
public CrystalTransform(CrystalTransform transform) {
this.sg = transform.sg;
this.transformId = transform.transformId;
this.matTransform = new Matrix4d(transform.matTransform);
this.crystalTranslation = new Point3i(transform.crystalTranslation);
}
public Matrix4d getMatTransform() {
return matTransform;
}
public void setMatTransform(Matrix4d matTransform) {
this.matTransform = matTransform;
}
public Point3i getCrystalTranslation() {
return crystalTranslation;
}
public void translate(Point3i translation) {
matTransform.m03 = matTransform.m03+translation.x;
matTransform.m13 = matTransform.m13+translation.y;
matTransform.m23 = matTransform.m23+translation.z;
crystalTranslation.add(translation);
}
/**
* Returns true if the given CrystalTransform is equivalent to this one.
* Two crystal transforms are equivalent if one is the inverse of the other, i.e.
* their transformation matrices multiplication is equal to the identity.
* @param other
* @return
*/
public boolean isEquivalent(CrystalTransform other) {
Matrix4d mul = new Matrix4d();
mul.mul(this.matTransform,other.matTransform);
if (mul.epsilonEquals(IDENTITY, 0.0001)) {
return true;
}
return false;
}
/**
* Tells whether this transformation is a pure crystal lattice translation,
* i.e. no rotational component and an integer translation vector.
* @return
*/
public boolean isPureCrystalTranslation() {
return (transformId==0 && (crystalTranslation.x!=0 || crystalTranslation.y!=0 || crystalTranslation.z!=0));
}
/**
* Tells whether this transformation is the identity: no rotation and no translation
* @return
*/
public boolean isIdentity() {
return (transformId==0 && crystalTranslation.x==0 && crystalTranslation.y==0 && crystalTranslation.z==0);
}
/**
* Tells whether this transformation is a pure translation:
* either a pure crystal (lattice) translation or a fractional (within
* unit cell) translation: space groups Ixxx, Cxxx, Fxxx have operators
* with fractional translations within the unit cell.
* @return
*/
public boolean isPureTranslation() {
if (isPureCrystalTranslation()) return true;
if (SpaceGroup.deltaComp(matTransform.m00,1,SpaceGroup.DELTA) &&
SpaceGroup.deltaComp(matTransform.m01,0,SpaceGroup.DELTA) &&
SpaceGroup.deltaComp(matTransform.m02,0,SpaceGroup.DELTA) &&
SpaceGroup.deltaComp(matTransform.m10,0,SpaceGroup.DELTA) &&
SpaceGroup.deltaComp(matTransform.m11,1,SpaceGroup.DELTA) &&
SpaceGroup.deltaComp(matTransform.m12,0,SpaceGroup.DELTA) &&
SpaceGroup.deltaComp(matTransform.m20,0,SpaceGroup.DELTA) &&
SpaceGroup.deltaComp(matTransform.m21,0,SpaceGroup.DELTA) &&
SpaceGroup.deltaComp(matTransform.m22,1,SpaceGroup.DELTA) &&
( Math.abs(matTransform.m03-0.0)>SpaceGroup.DELTA ||
Math.abs(matTransform.m13-0.0)>SpaceGroup.DELTA ||
Math.abs(matTransform.m23-0.0)>SpaceGroup.DELTA)) {
return true;
}
return false;
}
/**
* Tells whether this transformation contains a fractional translational
* component (whatever its rotational component). A fractional translation
* together with a rotation means a screw axis.
* @return
*/
public boolean isFractionalTranslation() {
if ((Math.abs(matTransform.m03-crystalTranslation.x)>SpaceGroup.DELTA) ||
(Math.abs(matTransform.m13-crystalTranslation.y)>SpaceGroup.DELTA) ||
(Math.abs(matTransform.m23-crystalTranslation.z)>SpaceGroup.DELTA)) {
return true;
}
return false;
}
/**
* Tells whether this transformation is a rotation disregarding the translational component,
* i.e. either pure rotation or screw rotation, but not improper rotation.
* @return
*/
public boolean isRotation() {
// if no SG, that means a non-crystallographic entry (e.g. NMR). We return false
if (sg==null) return false;
// this also takes care of case <0 (improper rotations): won't be considered as rotations
if (sg.getAxisFoldType(this.transformId)>1) return true;
return false;
}
/**
* Returns the TransformType of this transformation: AU, crystal translation, fractional translation
* , 2 3 4 6-fold rotations, 2 3 4 6-fold screw rotations, -1 -3 -2 -4 -6 inversions/rotoinversions.
* @return
*/
public TransformType getTransformType() {
// if no SG, that means a non-crystallographic entry (e.g. NMR). We return AU as type
if (sg==null) return TransformType.AU;
int foldType = sg.getAxisFoldType(this.transformId);
boolean isScrewOrGlide = false;
Vector3d translScrewComponent = getTranslScrewComponent();
if (Math.abs(translScrewComponent.x-0.0)>SpaceGroup.DELTA ||
Math.abs(translScrewComponent.y-0.0)>SpaceGroup.DELTA ||
Math.abs(translScrewComponent.z-0.0)>SpaceGroup.DELTA) {
isScrewOrGlide = true;
}
if (foldType>1) {
if (isScrewOrGlide) {
switch (foldType) {
case 2:
return TransformType.TWOFOLDSCREW;
case 3:
return TransformType.THREEFOLDSCREW;
case 4:
return TransformType.FOURFOLDSCREW;
case 6:
return TransformType.SIXFOLDSCREW;
default:
throw new NullPointerException("This transformation did not fall into any of the known types! This is most likely a bug.");
}
} else {
switch (foldType) {
case 2:
return TransformType.TWOFOLD;
case 3:
return TransformType.THREEFOLD;
case 4:
return TransformType.FOURFOLD;
case 6:
return TransformType.SIXFOLD;
default:
throw new NullPointerException("This transformation did not fall into any of the known types! This is most likely a bug.");
}
}
} else if (foldType<0) {
switch (foldType) {
case -1:
return TransformType.ONEBAR;
case -2:
if (isScrewOrGlide) {
return TransformType.GLIDE;
}
return TransformType.TWOBAR;
case -3:
return TransformType.THREEBAR;
case -4:
return TransformType.FOURBAR;
case -6:
return TransformType.SIXBAR;
default:
throw new NullPointerException("This transformation did not fall into any of the known types! This is most likely a bug.");
}
} else {
if (isIdentity()) {
return TransformType.AU;
}
if (isPureCrystalTranslation()) {
return TransformType.XTALTRANSL;
}
if (isFractionalTranslation()) {
return TransformType.CELLTRANSL;
}
throw new NullPointerException("This transformation did not fall into any of the known types! This is most likely a bug.");
}
}
public Vector3d getTranslScrewComponent() {
return getTranslScrewComponent(matTransform);
}
public int getTransformId() {
return transformId;
}
public void setTransformId(int transformId) {
this.transformId = transformId;
}
@Override
public String toString() {
return String.format("[%2d-(%s)]",transformId,toXYZString());
}
/**
* Expresses this transformation in terms of x,y,z fractional coordinates.
*
* Examples:
* @return
*/
public String toXYZString() {
StringBuilder str = new StringBuilder();
for(int i=0;i<3;i++) { //for each row
boolean emptyRow = true;
double coef; // TODO work with rational numbers
// X
coef = matTransform.getElement(i, 0);
// Three cases for coef: zero, one, non-one
if(abs(coef) > 1e-6 ) { // Non-zero
if( abs( abs(coef)-1 ) < 1e-6 ) { // +/- 1
if( coef < 0 ) {
str.append("-");
}
} else {
str.append(formatCoef(coef));
str.append("*");
}
str.append("x");
emptyRow = false;
}
// Y
coef = matTransform.getElement(i, 1);
if(abs(coef) > 1e-6 ) { // Non-zero
if( abs( abs(coef)-1 ) < 1e-6 ) { // +/- 1
if( coef < 0 ) {
str.append("-");
} else if( !emptyRow ) {
str.append("+");
}
} else {
if( !emptyRow && coef > 0) {
str.append("+");
}
str.append(formatCoef(coef));
str.append("*");
}
str.append("y");
emptyRow = false;
}
// Z
coef = matTransform.getElement(i, 2);
if(abs(coef) > 1e-6 ) { // Non-zero
if( abs( abs(coef)-1 ) < 1e-6 ) { // +/- 1
if( coef < 0 ) {
str.append("-");
} else if( !emptyRow ) {
str.append("+");
}
} else {
if( !emptyRow && coef > 0) {
str.append("+");
}
str.append(formatCoef(coef));
str.append("*");
}
str.append("z");
emptyRow = false;
}
// Intercept
coef = matTransform.getElement(i, 3);
if(abs(coef) > 1e-6 ) { // Non-zero
if( !emptyRow && coef > 0) {
str.append("+");
}
str.append(formatCoef(coef));
}
if(i<2) {
str.append(",");
}
}
return str.toString();
}
/**
* helper function to format simple fractions into rationals
* @param coef
* @return
*/
private String formatCoef(double coef) {
double tol = 1e-6; // rounding tolerance
// zero case
if( Math.abs(coef) < tol) {
return "0";
}
// integer case
long num = Math.round(coef);
if( Math.abs(num - coef) < tol) {
return Long.toString(num);
}
// Other small cases
for(int denom = 2; denom < 12; denom++ ) {
num = Math.round(coef*denom);
if( num - coef*denom < tol ) {
return String.format("%d/%d",num, denom);
}
}
// Give up and use floating point;
return String.format("%.3f", coef);
}
/**
* Given a transformation matrix containing a rotation and translation returns the
* screw component of the rotation.
* See http://www.crystallography.fr/mathcryst/pdf/Gargnano/Aroyo_Gargnano_1.pdf
* @param m
* @return
*/
public static Vector3d getTranslScrewComponent(Matrix4d m) {
int foldType = SpaceGroup.getRotAxisType(m);
// For reference see:
// http://www.crystallography.fr/mathcryst/pdf/Gargnano/Aroyo_Gargnano_1.pdf
Vector3d transl = null;
Matrix3d W =
new Matrix3d(m.m00,m.m01,m.m02,
m.m10,m.m11,m.m12,
m.m20,m.m21,m.m22);
if (foldType>=0) {
// the Y matrix: Y = W^k-1 + W^k-2 ... + W + I ; with k the fold type
Matrix3d Y = new Matrix3d(1,0,0, 0,1,0, 0,0,1);
Matrix3d Wk = new Matrix3d(1,0,0, 0,1,0, 0,0,1);
for (int k=0;k<foldType;k++) {
Wk.mul(W); // k=0 Wk=W, k=1 Wk=W^2, k=2 Wk=W^3, ... k=foldType-1, Wk=W^foldType
if (k!=foldType-1) Y.add(Wk);
}
transl = new Vector3d(m.m03, m.m13, m.m23);
Y.transform(transl);
transl.scale(1.0/foldType);
} else {
if (foldType==-2) { // there are glide planes only in -2
Matrix3d Y = new Matrix3d(1,0,0, 0,1,0, 0,0,1);
Y.add(W);
transl = new Vector3d(m.m03, m.m13, m.m23);
Y.transform(transl);
transl.scale(1.0/2.0);
} else { // for -1, -3, -4 and -6 there's nothing to do: fill with 0s
transl = new Vector3d(0,0,0);
}
}
return transl;
}
}