/*
* 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 org.biojava.nbio.structure.jama.EigenvalueDecomposition;
import org.biojava.nbio.structure.jama.Matrix;
import org.biojava.nbio.structure.xtal.io.TransfAlgebraicAdapter;
import org.slf4j.Logger;
import org.slf4j.LoggerFactory;
import javax.vecmath.AxisAngle4d;
import javax.vecmath.Matrix3d;
import javax.vecmath.Matrix4d;
import javax.vecmath.Vector3d;
import javax.xml.bind.JAXBContext;
import javax.xml.bind.JAXBException;
import javax.xml.bind.Marshaller;
import javax.xml.bind.annotation.XmlAccessType;
import javax.xml.bind.annotation.XmlAccessorType;
import javax.xml.bind.annotation.XmlRootElement;
import javax.xml.bind.annotation.adapters.XmlJavaTypeAdapter;
import java.io.ByteArrayOutputStream;
import java.io.PrintStream;
import java.io.Serializable;
import java.util.ArrayList;
import java.util.List;
import java.util.regex.Matcher;
import java.util.regex.Pattern;
/**
* A crystallographic space group. We store the standard numeric identifier,
* the international short symbol and the transformations corresponding to
* each space group (as Matrix4ds and in algebraic notation).
* The information for all (protein crystallography) space groups can be
* parsed from the XML file in the resource directory.
*
* See: http://en.wikipedia.org/wiki/Space_group
*
* @author duarte_j
* @see SymoplibParser
*/
@XmlRootElement(name = "SpaceGroup", namespace ="http://www.biojava.org")
@XmlAccessorType(XmlAccessType.PUBLIC_MEMBER)
public class SpaceGroup implements Serializable {
private static final long serialVersionUID = 1L;
private static final Logger logger = LoggerFactory.getLogger(SpaceGroup.class);
private static final Pattern splitPat1 = Pattern.compile("((?:[+-]?[XYZ])+)([+-][0-9/.]+)");
private static final Pattern splitPat2 = Pattern.compile("([+-]?[0-9/.]+)((?:[+-][XYZ])+)");
private static final Pattern coordPat = Pattern.compile("(?:([+-])?([XYZ]))+?"); // the last +? is for ungreedy matching
private static final Pattern transCoefPat = Pattern.compile("([-+]?[0-9.]+)(?:/([0-9.]+))?");
private static final Pattern nonEnantPat = Pattern.compile("[-abcmnd]");
protected static final double DELTA=0.0000001;
private int id;
private int multiplicity;
private int primitiveMultiplicity;
private String shortSymbol;
private String altShortSymbol;
private List<Matrix4d> transformations;
private List<String> transfAlgebraic;
private Vector3d[] cellTranslations; // in space groups I, C, F or H there are pure cell translations corresponding to recenterings
private AxisAngle4d[] axisAngles;
private int[] axisTypes; // indices of array are transformIds
private BravaisLattice bravLattice;
@SuppressWarnings("unused")
private SpaceGroup(){
// required by JAXB
}
public SpaceGroup(int id, int multiplicity, int primitiveMultiplicity, String shortSymbol, String altShortSymbol, BravaisLattice bravLattice) {
this.id = id;
this.multiplicity = multiplicity;
this.primitiveMultiplicity = primitiveMultiplicity;
this.shortSymbol = shortSymbol;
this.altShortSymbol = altShortSymbol;
transformations = new ArrayList<Matrix4d>(multiplicity);
transfAlgebraic = new ArrayList<String>(multiplicity);
cellTranslations = new Vector3d[multiplicity/primitiveMultiplicity];
this.bravLattice = bravLattice;
}
/**
* Get the space group for the given international short name, using
* the PDB format, e.g. 'P 21 21 21' or 'C 1 c 1'
* @param shortName
* @return the SpaceGroup or null if the shortName is not valid
* @see SymoplibParser#getSpaceGroup(String)
*/
public static SpaceGroup parseSpaceGroup(String shortName) {
return SymoplibParser.getSpaceGroup(shortName);
}
public void addTransformation(String transfAlgebraic) {
this.transfAlgebraic.add(transfAlgebraic);
this.transformations.add(getMatrixFromAlgebraic(transfAlgebraic));
}
protected void initializeCellTranslations() {
if ( cellTranslations != null && cellTranslations.length >0) {
// we already initialized this
return;
}
cellTranslations = new Vector3d[multiplicity/primitiveMultiplicity];
cellTranslations[0] = new Vector3d(0,0,0);
if ( transformations == null){
logger.warn("transformations == null" + this.toXML());
}
if (multiplicity==primitiveMultiplicity) {
return;
}
int fold = multiplicity/primitiveMultiplicity;
for (int n=1;n<fold;n++) {
if ( transformations.size() < (n* primitiveMultiplicity)){
logger.warn("WARNING number of transformations < " +(n*primitiveMultiplicity));
logger.warn(this.toXML());
}
Matrix4d t = transformations.get(n*primitiveMultiplicity);
cellTranslations[n] = new Vector3d(t.m03,t.m13,t.m23);
}
}
public int getMultiplicity() {
return multiplicity;
}
public int getPrimitiveMultiplicity() {
return primitiveMultiplicity;
}
public Vector3d[] getCellTranslations() {
return cellTranslations;
}
public Vector3d getCellTranslation(int i) {
return cellTranslations[i];
}
public static Matrix4d getMatrixFromAlgebraic(String transfAlgebraic) {
String[] parts = transfAlgebraic.toUpperCase().split(",");
double[] xCoef = convertAlgebraicStrToCoefficients(parts[0].trim());
double[] yCoef = convertAlgebraicStrToCoefficients(parts[1].trim());
double[] zCoef = convertAlgebraicStrToCoefficients(parts[2].trim());
Matrix4d mat = new Matrix4d();
mat.setIdentity();
mat.setRotation(new Matrix3d(xCoef[0],xCoef[1],xCoef[2],yCoef[0],yCoef[1],yCoef[2],zCoef[0],zCoef[1],zCoef[2]));
mat.setTranslation(new Vector3d(xCoef[3],yCoef[3],zCoef[3]));
return mat;
//return new Matrix4d(xCoef[0],xCoef[1],xCoef[2],xCoef[3],
// yCoef[0],yCoef[1],yCoef[2],yCoef[3],
// zCoef[0],zCoef[1],zCoef[2],zCoef[3],
// 0,0,0,1);
}
private static double[] convertAlgebraicStrToCoefficients(String algString) {
String letters = null;
String noLetters = null;
Matcher m = splitPat1.matcher(algString);
if (m.matches()) {
letters = m.group(1);
noLetters = m.group(2);
} else {
m = splitPat2.matcher(algString);
if (m.matches()) {
letters = m.group(2);
noLetters = m.group(1);
} else {
letters = algString;
}
}
double[] coefficients = new double[4];
m = coordPat.matcher(letters);
while(m.find()){
String sign = "";
if (m.group(1)!=null) {
sign = m.group(1);
}
double s = 1.0;
if (sign.equals("-")){
s = -1.0;
}
String coord = m.group(2);
if (coord.equals("X")) {
coefficients[0] = s;
} else if (coord.equals("Y")) {
coefficients[1] = s;
} else if (coord.equals("Z")) {
coefficients[2] = s;
}
}
if (noLetters!=null) {
m = transCoefPat.matcher(noLetters);
if (m.matches()) {
double num = Double.parseDouble(m.group(1));
double den = 1;
if (m.group(2)!=null) {
den = Double.parseDouble(m.group(2));
}
coefficients[3] = num/den;
}
} else {
coefficients[3]=0;
}
return coefficients;
}
/**
* Gets the standard numeric identifier for the space group.
* See for example http://en.wikipedia.org/wiki/Space_group
* or the IUCr crystallographic tables
* @return
*/
public int getId() {
return id;
}
/**
* Gets the international short name (as used in PDB),
* e.g. "P 21 21 21" or "C 1 c 1"
* @return
*/
public String getShortSymbol() {
return shortSymbol;
}
/**
* Gets the alternative international short name (as sometimes used in PDB),
* e.g. "I 1 2 1" instead of "I 2"
* @return
*/
public String getAltShortSymbol() {
return altShortSymbol;
}
/**
* Gets all transformations except for the identity in crystal axes basis.
* @return
*/
public List<Matrix4d> getTransformations() {
List<Matrix4d> transfs = new ArrayList<Matrix4d>();
for (int i=1;i<this.transformations.size();i++){
transfs.add(transformations.get(i));
}
return transfs;
}
private void calcRotAxesAndAngles() {
axisAngles = new AxisAngle4d[multiplicity];
// identity operator (transformId==0)
axisAngles[0] = new AxisAngle4d(new Vector3d(0,0,0), 0.0);
for (int i=1;i<this.transformations.size();i++){
Matrix3d r = new Matrix3d(transformations.get(i).m00,transformations.get(i).m01,transformations.get(i).m02,
transformations.get(i).m10,transformations.get(i).m11,transformations.get(i).m12,
transformations.get(i).m20,transformations.get(i).m21,transformations.get(i).m22);
axisAngles[i] = getRotAxisAndAngle(r);
}
}
/**
* Calculates the axis fold type (1, 2, 3, 4, 5, 6 for rotations or -1, -2, -3, -4, -6 improper rotations)
* from the trace of the rotation matrix, see for instance
* http://www.crystallography.fr/mathcryst/pdf/Gargnano/Aroyo_Gargnano_1.pdf
*/
private void calcAxisFoldTypes() {
axisTypes = new int[multiplicity];
for (int i=0;i<this.transformations.size();i++){
axisTypes[i] = getRotAxisType(transformations.get(i));
}
}
public AxisAngle4d getRotAxisAngle(int transformId) {
if (this.axisAngles == null) calcRotAxesAndAngles();
return this.axisAngles[transformId];
}
/**
* Returns true if both given transform ids belong to the same crystallographic axis (a, b or c)
* For two non-rotation transformations (i.e. identity operators) it returns true
* @param tId1
* @param tId2
* @return
*/
public boolean areInSameAxis(int tId1, int tId2) {
if (tId1==tId2) return true;
if (axisAngles== null) calcRotAxesAndAngles();
if (getAxisFoldType(tId1)==1 && getAxisFoldType(tId2)==1) return true;
// we can't deal yet with improper rotations: we return false whenever either of them is improper
if (getAxisFoldType(tId1)<0 || getAxisFoldType(tId2)<0) return false;
Vector3d axis1 = new Vector3d(axisAngles[tId1].x, axisAngles[tId1].y, axisAngles[tId1].z);
Vector3d axis2 = new Vector3d(axisAngles[tId2].x, axisAngles[tId2].y, axisAngles[tId2].z);
// TODO revise: we might need to consider that the 2 are in same direction but opposite senses
// the method is not used at the moment anyway
if (deltaComp(axis1.angle(axis2), 0.0, DELTA)) return true;
return false;
}
/**
* Given a transformId returns the type of axis of rotation: 1 (no rotation), 2, 3, 4 or 6 -fold
* and for improper rotations: -1, -2, -3, -4 and -6
*
* @param transformId
* @return
*/
public int getAxisFoldType(int transformId) {
if (axisTypes== null) calcAxisFoldTypes();
return axisTypes[transformId];
}
/**
* Gets a transformation by index expressed in crystal axes basis.
* Index 0 corresponds always to the identity transformation.
* Beware the returned Matrix4d is not a copy but it stays linked
* to the one stored in this SpaceGroup object
* @param i
* @return
*/
public Matrix4d getTransformation(int i) {
return transformations.get(i);
}
/**
* Gets a transformation algebraic string given its index.
* Index 0 corresponds always to the identity transformation.
* @param i
* @return
*/
public String getTransfAlgebraic(int i) {
return transfAlgebraic.get(i);
}
@Override
public int hashCode() {
final int prime = 31;
int result = 1;
result = prime * result + id;
return result;
}
@Override
public boolean equals(Object o) {
if (! (o instanceof SpaceGroup)) {
return false;
}
SpaceGroup other = (SpaceGroup) o;
if (other.getId()==this.getId()) {
return true;
}
return false;
}
/**
* Gets the number of symmetry operators corresponding to this SpaceGroup (counting
* the identity operator)
* @return
*/
public int getNumOperators() {
return this.transformations.size();
}
public static String getAlgebraicFromMatrix(Matrix4d m) {
String x = formatAlg(m.m00,m.m01,m.m02,m.m03);
String y = formatAlg(m.m10,m.m11,m.m12,m.m13);
String z = formatAlg(m.m20,m.m21,m.m22,m.m23);
String alg = x+","+y+","+z;
return alg;
}
private static String formatAlg(double xcoef, double ycoef, double zcoef, double trans) {
boolean[] leading = {false,false,false};
if (xcoef!=0) {
leading[0] = true;
} else if (ycoef!=0) {
leading[1] = true;
} else if (zcoef!=0) {
leading[2] = true;
}
String x = deltaComp(xcoef,0,DELTA)?"":formatCoef(xcoef,leading[0])+"X";
String y = deltaComp(ycoef,0,DELTA)?"":formatCoef(ycoef,leading[1])+"Y";
String z = deltaComp(zcoef,0,DELTA)?"":formatCoef(zcoef, leading[2])+"Z";
String t = deltaComp(trans,0,DELTA)?"":formatTransCoef(trans);
return x+y+z+t;
}
private static String formatCoef(double c, boolean leading) {
if (leading) {
return (deltaComp(Math.abs(c),1,DELTA)?(c>0?"":"-"):String.format("%4.2f",c));
} else {
return (deltaComp(Math.abs(c),1,DELTA)?(c>0?"+":"-"):String.format("%+4.2f",c));
}
}
private static String formatTransCoef(double c) {
if (Math.abs((Math.rint(c)-c))<DELTA) { // this is an integer
return String.format("%+d",(int)Math.rint(c));
} else { // it is a fraction
int num,den;
int floor = (int)Math.floor(c);
double decPart = c - floor;
if (deltaComp(decPart,0.3333333,DELTA)) {
num=1;den=3;
} else if (deltaComp(decPart,0.6666667,DELTA)) {
num=2;den=3;
} else if (deltaComp(decPart,0.2500000,DELTA)) {
num=1;den=4;
} else if (deltaComp(decPart,0.5000000,DELTA)) {
num=1;den=2;
} else if (deltaComp(decPart,0.7500000,DELTA)) {
num=3;den=4;
} else if (deltaComp(decPart,0.1666667,DELTA)) {
num=1;den=6;
} else if (deltaComp(decPart,0.8333333,DELTA)) {
num=5;den=6;
} else {
num=0;den=0; // this in an error
}
num = floor*den+num;
return String.format("%+d/%d", num,den);
//return String.format("%+4.2f",c);
}
}
protected static boolean deltaComp(double d1, double d2, double delta) {
return Math.abs(d1-d2)<delta;
}
public BravaisLattice getBravLattice() {
return bravLattice;
}
public boolean isEnantiomorphic() {
Matcher m = nonEnantPat.matcher(shortSymbol);
if (m.find()) {
return false;
}
return true;
}
/**
* Given a rotation matrix calculates the rotation axis and angle for it.
* The angle is calculated from the trace, the axis from the eigenvalue
* decomposition.
* If given matrix is improper rotation or identity matrix then
* axis (0,0,0) and angle 0 are returned.
* @param m
* @return
* @throws IllegalArgumentException if given matrix is not a rotation matrix (determinant not 1 or -1)
*/
public static AxisAngle4d getRotAxisAndAngle(Matrix3d m) {
double determinant = m.determinant();
if (!(Math.abs(determinant)-1.0<DELTA)) throw new IllegalArgumentException("Given matrix is not a rotation matrix");
AxisAngle4d axisAndAngle = new AxisAngle4d(new Vector3d(0,0,0),0);
double[] d = {m.m00,m.m10,m.m20,
m.m01,m.m11,m.m21,
m.m02,m.m12,m.m22};
Matrix r = new Matrix(d,3);
if (!deltaComp(r.det(), 1.0, DELTA)) {
// improper rotation: we return axis 0,0,0 and angle 0
return axisAndAngle;
}
EigenvalueDecomposition evd = new EigenvalueDecomposition(r);
Matrix eval = evd.getD();
if (deltaComp(eval.get(0, 0),1.0,DELTA) && deltaComp(eval.get(1, 1),1.0,DELTA) && deltaComp(eval.get(2, 2),1.0,DELTA)) {
// the rotation is an identity: we return axis 0,0,0 and angle 0
return axisAndAngle;
}
int indexOfEv1;
for (indexOfEv1=0;indexOfEv1<3;indexOfEv1++) {
if (deltaComp(eval.get(indexOfEv1, indexOfEv1),1,DELTA)) break;
}
Matrix evec = evd.getV();
axisAndAngle.set(new Vector3d(evec.get(0,indexOfEv1), evec.get(1, indexOfEv1), evec.get(2, indexOfEv1)),
Math.acos((eval.trace()-1.0)/2.0));
return axisAndAngle;
}
/**
* Given a transformation matrix containing a rotation returns the type of rotation:
* 1 for identity, 2 for 2-fold rotation, 3 for 3-fold rotation, 4 for 4-fold rotation,
* 6 for 6-fold rotation,
* -1 for inversions, -2 for mirror planes, -3 for 3-fold improper rotation,
* -4 for 4-fold improper rotation and -6 for 6-fold improper rotation
* @param m
* @return
*/
public static int getRotAxisType(Matrix4d m) {
int axisType = 0;
Matrix3d rot = new Matrix3d(m.m00,m.m01,m.m02,
m.m10,m.m11,m.m12,
m.m20,m.m21,m.m22);
double determinant = rot.determinant();
if (!deltaComp(determinant,1.0,DELTA) && !deltaComp(determinant, -1.0, DELTA)) {
throw new IllegalArgumentException("Given matrix does not seem to be a rotation matrix.");
}
int trace = (int)(rot.m00+rot.m11+rot.m22);
if (determinant>0) {
switch (trace) {
case 3:
axisType=1;
break;
case -1:
axisType=2;
break;
case 0:
axisType=3;
break;
case 1:
axisType=4;
break;
case 2:
axisType=6;
break;
default:
throw new RuntimeException("Trace of transform does not correspond to one of the expected types. This is most likely a bug");
}
} else {
switch (trace) {
case -3:
axisType=-1;
break;
case 1:
axisType=-2;
break;
case 0:
axisType=-3;
break;
case -1:
axisType=-4;
break;
case -2:
axisType=-6;
break;
default:
throw new RuntimeException("Trace of transform does not correspond to one of the expected types. This is most likely a bug");
}
}
return axisType;
}
@Override
public String toString() {
return getShortSymbol();
}
public String toXML(){
JAXBContext jaxbContextStringSortedSet = null;
try {
jaxbContextStringSortedSet= JAXBContext.newInstance(SpaceGroup.class);
} catch (JAXBException e){
logger.error("Error converting to XML",e);
return null;
}
ByteArrayOutputStream baos = new ByteArrayOutputStream();
PrintStream ps = new PrintStream(baos);
try {
Marshaller m = jaxbContextStringSortedSet.createMarshaller();
m.setProperty(Marshaller.JAXB_FORMATTED_OUTPUT, Boolean.TRUE);
m.marshal( this, ps);
} catch (JAXBException e){
logger.error("Error converting to XML",e);
}
return baos.toString();
}
//@XmlElementWrapper(name="transfAlgebraicList", namespace="http://www.biojava.org")
//@XmlElement
@XmlJavaTypeAdapter(TransfAlgebraicAdapter.class)
public List<String> getTransfAlgebraic() {
return transfAlgebraic;
}
public void setTransfAlgebraic(List<String> transfAlgebraic) {
//System.out.println("setting transfAlgebraic " + transfAlgebraic);
if ( transformations == null || transformations.size() == 0)
transformations = new ArrayList<Matrix4d>(transfAlgebraic.size());
if ( this.transfAlgebraic == null || this.transfAlgebraic.size() == 0)
this.transfAlgebraic = new ArrayList<String>(transfAlgebraic.size());
for ( String transf : transfAlgebraic){
addTransformation(transf);
}
}
public int[] getAxisTypes() {
return axisTypes;
}
public void setAxisTypes(int[] axisTypes) {
this.axisTypes = axisTypes;
}
public static long getSerialversionuid() {
return serialVersionUID;
}
public static Pattern getSplitpat1() {
return splitPat1;
}
public static Pattern getSplitpat2() {
return splitPat2;
}
public static Pattern getCoordpat() {
return coordPat;
}
public static Pattern getTranscoefpat() {
return transCoefPat;
}
public static Pattern getNonenantpat() {
return nonEnantPat;
}
public static double getDelta() {
return DELTA;
}
public void setId(int id) {
this.id = id;
}
public void setMultiplicity(int multiplicity) {
this.multiplicity = multiplicity;
}
public void setPrimitiveMultiplicity(int primitiveMultiplicity) {
this.primitiveMultiplicity = primitiveMultiplicity;
}
public void setShortSymbol(String shortSymbol) {
this.shortSymbol = shortSymbol;
}
public void setAltShortSymbol(String altShortSymbol) {
this.altShortSymbol = altShortSymbol;
}
public void setBravLattice(BravaisLattice bravLattice) {
this.bravLattice = bravLattice;
}
}