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