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.Matrix3d; 024import javax.vecmath.Matrix4d; 025import javax.vecmath.Point3i; 026import javax.vecmath.Vector3d; 027import java.io.Serializable; 028 029import static java.lang.Math.abs; 030 031 032/** 033 * Representation of a transformation in a crystal: 034 * - a transformation id (each of the transformations in a space group, 0 to m) 035 * - a crystal translation 036 * The transformation matrix in crystal basis is stored, representing the basic 037 * transformation together with the crystal translation. 038 * Contains methods to check for equivalent transformations. 039 * 040 * 041 * @author duarte_j 042 * 043 */ 044public class CrystalTransform implements Serializable { 045 046 private static final long serialVersionUID = 1L; 047 048 049 public static final Matrix4d IDENTITY = new Matrix4d(1,0,0,0, 0,1,0,0, 0,0,1,0, 0,0,0,1); 050 051 /** 052 * The space group to which this transform belongs 053 */ 054 private final SpaceGroup sg; 055 056 /** 057 * The transform id corresponding to the SpaceGroup's transform indices. 058 * From 0 (identity) to m (m=number of symmetry operations of the space group) 059 * It is unique within the unit cell but equivalent units of different crystal unit cells 060 * will have same id 061 */ 062 private int transformId; 063 064 /** 065 * The 4-dimensional matrix transformation in crystal basis. 066 * Note that the translational component of this matrix is not necessarily 067 * identical to crystalTranslation since some operators have fractional 068 * translations within the cell 069 */ 070 private Matrix4d matTransform; 071 072 /** 073 * The crystal translation (always integer) 074 */ 075 private Point3i crystalTranslation; 076 077 078 /** 079 * Creates a new CrystalTransform representing the identity transform 080 * in cell (0,0,0) 081 */ 082 public CrystalTransform(SpaceGroup sg) { 083 this.sg = sg; 084 this.transformId = 0; 085 this.matTransform = (Matrix4d)IDENTITY.clone(); 086 this.crystalTranslation = new Point3i(0,0,0); 087 } 088 089 /** 090 * Represents the n-th transform 091 * @param sg 092 * @param transformId 093 */ 094 public CrystalTransform(SpaceGroup sg, int transformId) { 095 this.sg = sg; 096 this.transformId = transformId; 097 if (sg==null && transformId==0) { 098 this.matTransform = (Matrix4d)IDENTITY.clone(); 099 } else if (sg==null) { 100 throw new IllegalArgumentException("Space Group cannot be null if transformId!=0"); 101 } else { 102 this.matTransform = (Matrix4d)sg.getTransformation(transformId).clone(); 103 } 104 this.crystalTranslation = new Point3i(0,0,0); 105 } 106 107 /** 108 * Copy constructor 109 * @param transform 110 */ 111 public CrystalTransform(CrystalTransform transform) { 112 this.sg = transform.sg; 113 this.transformId = transform.transformId; 114 this.matTransform = new Matrix4d(transform.matTransform); 115 this.crystalTranslation = new Point3i(transform.crystalTranslation); 116 } 117 118 public Matrix4d getMatTransform() { 119 return matTransform; 120 } 121 122 public void setMatTransform(Matrix4d matTransform) { 123 this.matTransform = matTransform; 124 } 125 126 public Point3i getCrystalTranslation() { 127 return crystalTranslation; 128 } 129 130 public void translate(Point3i translation) { 131 matTransform.m03 = matTransform.m03+translation.x; 132 matTransform.m13 = matTransform.m13+translation.y; 133 matTransform.m23 = matTransform.m23+translation.z; 134 135 crystalTranslation.add(translation); 136 137 } 138 139 /** 140 * Returns true if the given CrystalTransform is equivalent to this one. 141 * Two crystal transforms are equivalent if one is the inverse of the other, i.e. 142 * their transformation matrices multiplication is equal to the identity. 143 * @param other 144 * @return 145 */ 146 public boolean isEquivalent(CrystalTransform other) { 147 Matrix4d mul = new Matrix4d(); 148 mul.mul(this.matTransform,other.matTransform); 149 150 if (mul.epsilonEquals(IDENTITY, 0.0001)) { 151 return true; 152 } 153 return false; 154 } 155 156 /** 157 * Tells whether this transformation is a pure crystal lattice translation, 158 * i.e. no rotational component and an integer translation vector. 159 * @return 160 */ 161 public boolean isPureCrystalTranslation() { 162 return (transformId==0 && (crystalTranslation.x!=0 || crystalTranslation.y!=0 || crystalTranslation.z!=0)); 163 } 164 165 /** 166 * Tells whether this transformation is the identity: no rotation and no translation 167 * @return 168 */ 169 public boolean isIdentity() { 170 return (transformId==0 && crystalTranslation.x==0 && crystalTranslation.y==0 && crystalTranslation.z==0); 171 } 172 173 /** 174 * Tells whether this transformation is a pure translation: 175 * either a pure crystal (lattice) translation or a fractional (within 176 * unit cell) translation: space groups Ixxx, Cxxx, Fxxx have operators 177 * with fractional translations within the unit cell. 178 * @return 179 */ 180 public boolean isPureTranslation() { 181 if (isPureCrystalTranslation()) return true; 182 if (SpaceGroup.deltaComp(matTransform.m00,1,SpaceGroup.DELTA) && 183 SpaceGroup.deltaComp(matTransform.m01,0,SpaceGroup.DELTA) && 184 SpaceGroup.deltaComp(matTransform.m02,0,SpaceGroup.DELTA) && 185 186 SpaceGroup.deltaComp(matTransform.m10,0,SpaceGroup.DELTA) && 187 SpaceGroup.deltaComp(matTransform.m11,1,SpaceGroup.DELTA) && 188 SpaceGroup.deltaComp(matTransform.m12,0,SpaceGroup.DELTA) && 189 190 SpaceGroup.deltaComp(matTransform.m20,0,SpaceGroup.DELTA) && 191 SpaceGroup.deltaComp(matTransform.m21,0,SpaceGroup.DELTA) && 192 SpaceGroup.deltaComp(matTransform.m22,1,SpaceGroup.DELTA) && 193 ( Math.abs(matTransform.m03-0.0)>SpaceGroup.DELTA || 194 Math.abs(matTransform.m13-0.0)>SpaceGroup.DELTA || 195 Math.abs(matTransform.m23-0.0)>SpaceGroup.DELTA)) { 196 return true; 197 } 198 199 return false; 200 } 201 202 /** 203 * Tells whether this transformation contains a fractional translational 204 * component (whatever its rotational component). A fractional translation 205 * together with a rotation means a screw axis. 206 * @return 207 */ 208 public boolean isFractionalTranslation() { 209 if ((Math.abs(matTransform.m03-crystalTranslation.x)>SpaceGroup.DELTA) || 210 (Math.abs(matTransform.m13-crystalTranslation.y)>SpaceGroup.DELTA) || 211 (Math.abs(matTransform.m23-crystalTranslation.z)>SpaceGroup.DELTA)) { 212 return true; 213 } 214 return false; 215 } 216 217 /** 218 * Tells whether this transformation is a rotation disregarding the translational component, 219 * i.e. either pure rotation or screw rotation, but not improper rotation. 220 * @return 221 */ 222 public boolean isRotation() { 223 // if no SG, that means a non-crystallographic entry (e.g. NMR). We return false 224 if (sg==null) return false; 225 226 // this also takes care of case <0 (improper rotations): won't be considered as rotations 227 if (sg.getAxisFoldType(this.transformId)>1) return true; 228 229 return false; 230 } 231 232 /** 233 * Returns the TransformType of this transformation: AU, crystal translation, fractional translation 234 * , 2 3 4 6-fold rotations, 2 3 4 6-fold screw rotations, -1 -3 -2 -4 -6 inversions/rotoinversions. 235 * @return 236 */ 237 public TransformType getTransformType() { 238 239 // if no SG, that means a non-crystallographic entry (e.g. NMR). We return AU as type 240 if (sg==null) return TransformType.AU; 241 242 int foldType = sg.getAxisFoldType(this.transformId); 243 boolean isScrewOrGlide = false; 244 Vector3d translScrewComponent = getTranslScrewComponent(); 245 if (Math.abs(translScrewComponent.x-0.0)>SpaceGroup.DELTA || 246 Math.abs(translScrewComponent.y-0.0)>SpaceGroup.DELTA || 247 Math.abs(translScrewComponent.z-0.0)>SpaceGroup.DELTA) { 248 249 isScrewOrGlide = true; 250 } 251 252 if (foldType>1) { 253 254 if (isScrewOrGlide) { 255 switch (foldType) { 256 case 2: 257 return TransformType.TWOFOLDSCREW; 258 case 3: 259 return TransformType.THREEFOLDSCREW; 260 case 4: 261 return TransformType.FOURFOLDSCREW; 262 case 6: 263 return TransformType.SIXFOLDSCREW; 264 default: 265 throw new NullPointerException("This transformation did not fall into any of the known types! This is most likely a bug."); 266 } 267 } else { 268 switch (foldType) { 269 case 2: 270 return TransformType.TWOFOLD; 271 case 3: 272 return TransformType.THREEFOLD; 273 case 4: 274 return TransformType.FOURFOLD; 275 case 6: 276 return TransformType.SIXFOLD; 277 default: 278 throw new NullPointerException("This transformation did not fall into any of the known types! This is most likely a bug."); 279 } 280 } 281 282 } else if (foldType<0) { 283 switch (foldType) { 284 case -1: 285 return TransformType.ONEBAR; 286 case -2: 287 if (isScrewOrGlide) { 288 return TransformType.GLIDE; 289 } 290 return TransformType.TWOBAR; 291 case -3: 292 return TransformType.THREEBAR; 293 case -4: 294 return TransformType.FOURBAR; 295 case -6: 296 return TransformType.SIXBAR; 297 default: 298 throw new NullPointerException("This transformation did not fall into any of the known types! This is most likely a bug."); 299 } 300 } else { 301 if (isIdentity()) { 302 return TransformType.AU; 303 } 304 if (isPureCrystalTranslation()) { 305 return TransformType.XTALTRANSL; 306 } 307 if (isFractionalTranslation()) { 308 return TransformType.CELLTRANSL; 309 } 310 throw new NullPointerException("This transformation did not fall into any of the known types! This is most likely a bug."); 311 } 312 313 } 314 315 public Vector3d getTranslScrewComponent() { 316 317 return getTranslScrewComponent(matTransform); 318 319 } 320 321 public int getTransformId() { 322 return transformId; 323 } 324 325 public void setTransformId(int transformId) { 326 this.transformId = transformId; 327 } 328 329 @Override 330 public String toString() { 331 return String.format("[%2d-(%s)]",transformId,toXYZString()); 332 } 333 334 /** 335 * Expresses this transformation in terms of x,y,z fractional coordinates. 336 * 337 * Examples: 338 * @return 339 */ 340 public String toXYZString() { 341 StringBuilder str = new StringBuilder(); 342 343 for(int i=0;i<3;i++) { //for each row 344 boolean emptyRow = true; 345 346 double coef; // TODO work with rational numbers 347 348 349 // X 350 coef = matTransform.getElement(i, 0); 351 352 // Three cases for coef: zero, one, non-one 353 if(abs(coef) > 1e-6 ) { // Non-zero 354 if( abs( abs(coef)-1 ) < 1e-6 ) { // +/- 1 355 if( coef < 0 ) { 356 str.append("-"); 357 } 358 } else { 359 str.append(formatCoef(coef)); 360 str.append("*"); 361 } 362 str.append("x"); 363 emptyRow = false; 364 } 365 366 // Y 367 coef = matTransform.getElement(i, 1); 368 369 if(abs(coef) > 1e-6 ) { // Non-zero 370 if( abs( abs(coef)-1 ) < 1e-6 ) { // +/- 1 371 if( coef < 0 ) { 372 str.append("-"); 373 } else if( !emptyRow ) { 374 str.append("+"); 375 } 376 } else { 377 if( !emptyRow && coef > 0) { 378 str.append("+"); 379 } 380 str.append(formatCoef(coef)); 381 str.append("*"); 382 } 383 str.append("y"); 384 emptyRow = false; 385 } 386 387 // Z 388 coef = matTransform.getElement(i, 2); 389 390 if(abs(coef) > 1e-6 ) { // Non-zero 391 if( abs( abs(coef)-1 ) < 1e-6 ) { // +/- 1 392 if( coef < 0 ) { 393 str.append("-"); 394 } else if( !emptyRow ) { 395 str.append("+"); 396 } 397 } else { 398 if( !emptyRow && coef > 0) { 399 str.append("+"); 400 } 401 str.append(formatCoef(coef)); 402 str.append("*"); 403 } 404 str.append("z"); 405 emptyRow = false; 406 } 407 408 // Intercept 409 coef = matTransform.getElement(i, 3); 410 411 if(abs(coef) > 1e-6 ) { // Non-zero 412 if( !emptyRow && coef > 0) { 413 str.append("+"); 414 } 415 str.append(formatCoef(coef)); 416 } 417 418 if(i<2) { 419 str.append(","); 420 } 421 } 422 423 424 return str.toString(); 425 } 426 /** 427 * helper function to format simple fractions into rationals 428 * @param coef 429 * @return 430 */ 431 private String formatCoef(double coef) { 432 double tol = 1e-6; // rounding tolerance 433 434 // zero case 435 if( Math.abs(coef) < tol) { 436 return "0"; 437 } 438 439 // integer case 440 long num = Math.round(coef); 441 if( Math.abs(num - coef) < tol) { 442 return Long.toString(num); 443 } 444 445 // Other small cases 446 for(int denom = 2; denom < 12; denom++ ) { 447 num = Math.round(coef*denom); 448 if( num - coef*denom < tol ) { 449 return String.format("%d/%d",num, denom); 450 } 451 } 452 453 // Give up and use floating point; 454 return String.format("%.3f", coef); 455 } 456 457 /** 458 * Given a transformation matrix containing a rotation and translation returns the 459 * screw component of the rotation. 460 * See http://www.crystallography.fr/mathcryst/pdf/Gargnano/Aroyo_Gargnano_1.pdf 461 * @param m 462 * @return 463 */ 464 public static Vector3d getTranslScrewComponent(Matrix4d m) { 465 466 int foldType = SpaceGroup.getRotAxisType(m); 467 // For reference see: 468 // http://www.crystallography.fr/mathcryst/pdf/Gargnano/Aroyo_Gargnano_1.pdf 469 470 Vector3d transl = null; 471 472 Matrix3d W = 473 new Matrix3d(m.m00,m.m01,m.m02, 474 m.m10,m.m11,m.m12, 475 m.m20,m.m21,m.m22); 476 477 if (foldType>=0) { 478 479 // the Y matrix: Y = W^k-1 + W^k-2 ... + W + I ; with k the fold type 480 Matrix3d Y = new Matrix3d(1,0,0, 0,1,0, 0,0,1); 481 Matrix3d Wk = new Matrix3d(1,0,0, 0,1,0, 0,0,1); 482 483 for (int k=0;k<foldType;k++) { 484 Wk.mul(W); // k=0 Wk=W, k=1 Wk=W^2, k=2 Wk=W^3, ... k=foldType-1, Wk=W^foldType 485 if (k!=foldType-1) Y.add(Wk); 486 } 487 488 transl = new Vector3d(m.m03, m.m13, m.m23); 489 Y.transform(transl); 490 491 transl.scale(1.0/foldType); 492 493 } else { 494 495 if (foldType==-2) { // there are glide planes only in -2 496 Matrix3d Y = new Matrix3d(1,0,0, 0,1,0, 0,0,1); 497 Y.add(W); 498 499 transl = new Vector3d(m.m03, m.m13, m.m23); 500 Y.transform(transl); 501 502 transl.scale(1.0/2.0); 503 } else { // for -1, -3, -4 and -6 there's nothing to do: fill with 0s 504 transl = new Vector3d(0,0,0); 505 } 506 } 507 508 return transl; 509 } 510}