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.align.util; 022 023import org.biojava.nbio.structure.Atom; 024import org.biojava.nbio.structure.AtomImpl; 025import org.biojava.nbio.structure.Calc; 026import org.biojava.nbio.structure.StructureException; 027import org.biojava.nbio.structure.align.model.AFPChain; 028import org.biojava.nbio.structure.contact.Pair; 029import org.biojava.nbio.structure.jama.Matrix; 030 031import javax.vecmath.AxisAngle4d; 032import javax.vecmath.Matrix4d; 033import javax.vecmath.Vector3d; 034 035import java.io.StringWriter; 036 037/** 038 * Calculates the rotation axis for an alignment 039 * 040 * <p>A superposition of two structures is generally represented as a rotation 041 * matrix plus a translation vector. However, it can also be represented as an 042 * axis of rotation plus some translation. 043 * 044 * <p>This class calculates the rotation axis and stores it as four properties: 045 * <ul><li>A unit vector parallel to the rotation axis ({@link #getRotationAxis()}) 046 * <li>The angle of rotation ({@link #getAngle()}) 047 * <li>A point on the rotation axis ({@link #getRotationPos()}) 048 * <li>Some translation parallel to the axis ({@link #getScrewTranslation()}) 049 * </ul> 050 * 051 * <p>The axis of rotation is poorly defined and numerically unstable for small 052 * angles. Therefore it's direction is left as null for angles less than 053 * {@link #MIN_ANGLE}. 054 * 055 * @author Spencer Bliven 056 * 057 */ 058public final class RotationAxis { 059 060 /** 061 * Minimum angle to calculate rotation axes. 5 degrees. 062 */ 063 static final double MIN_ANGLE = Math.toRadians(5.); 064 065 private double theta; 066 private Atom rotationAxis; // axis of rotation 067 private Atom rotationPos; // a point on the axis of rotation 068 private Atom screwTranslation; //translation parallel to the axis of rotation 069 private Atom otherTranslation; // translation perpendicular to the axis of rotation 070 071 /** 072 * The rotation angle 073 * @return the angle, in radians 074 */ 075 public double getAngle() { 076 return theta; 077 } 078 079 /** 080 * Get a unit vector along the rotation axis 081 * @return rotationAxis 082 */ 083 public Atom getRotationAxis() { 084 return rotationAxis; 085 } 086 087 /** 088 * Returns the rotation axis and angle in a single javax.vecmath.AxisAngle4d object 089 * @return 090 */ 091 public AxisAngle4d getAxisAngle4d() { 092 return new AxisAngle4d(rotationAxis.getX(),rotationAxis.getY(),rotationAxis.getZ(),theta); 093 } 094 095 /** 096 * Get a position on the rotation axis. 097 * 098 * Specifically, project the origin onto the rotation axis 099 * @return rotationPos 100 */ 101 public Atom getRotationPos() { 102 return rotationPos; 103 } 104 105 /** 106 * Get the component of translation parallel to the axis of rotation 107 * @return 108 */ 109 public Atom getScrewTranslation() { 110 return screwTranslation; 111 } 112 113 public Vector3d getVector3dScrewTranslation() { 114 return new Vector3d(screwTranslation.getX(),screwTranslation.getY(),screwTranslation.getZ()); 115 } 116 117 public double getTranslation() { 118 return Calc.amount(screwTranslation); 119 } 120 121 /** 122 * Get the component of translation perpendicular to the axis of rotation. 123 * This isn't particularly meaningful, but is calculated internally and 124 * was useful for debugging. 125 * @return 126 */ 127 @Deprecated 128 public Atom getOtherTranslation() { 129 return otherTranslation; 130 } 131 132 /** 133 * Calculate the rotation axis for the first block of an AFPChain 134 * @param afpChain 135 * @throws StructureException 136 * @throws NullPointerException if afpChain does not contain a valid rotation matrix and shift vector 137 */ 138 public RotationAxis(AFPChain afpChain) throws StructureException { 139 if(afpChain.getAlnLength() < 1) { 140 throw new StructureException("No aligned residues"); 141 } 142 init(afpChain.getBlockRotationMatrix()[0],afpChain.getBlockShiftVector()[0]); 143 } 144 145 /** 146 * Create a rotation axis from a vector, a point, and an angle. 147 * 148 * The result will be a pure rotation, with no screw component. 149 * @param axis A vector parallel to the axis of rotation 150 * @param pos A point on the axis of rotation 151 * @param theta The angle to rotate (radians) 152 */ 153 public RotationAxis(Atom axis, Atom pos, double theta) { 154 this.rotationAxis = Calc.unitVector(axis); 155 this.rotationPos = (Atom) pos.clone(); 156 this.theta = theta; 157 this.screwTranslation = new AtomImpl(); //zero 158 this.otherTranslation = null; //deprecated 159 } 160 161 /** 162 * Determine the location of the rotation axis based on a rotation matrix and a translation vector 163 * @param rotation 164 * @param translation 165 */ 166 public RotationAxis(Matrix rotation, Atom translation) { 167 init(rotation, translation); 168 } 169 170 /** 171 * Create a rotation axis from a Matrix4d containing a rotational 172 * component and a translational component. 173 * 174 * @param transform 175 */ 176 public RotationAxis(Matrix4d transform) { 177 178 Matrix rot = Calc.getRotationMatrix(transform); 179 Atom transl = Calc.getTranslationVector(transform); 180 init(rot,transl); 181 } 182 183 /** 184 * Get the rotation matrix corresponding to this axis 185 * @return A 3x3 rotation matrix 186 */ 187 public Matrix getRotationMatrix() { 188 return getRotationMatrix(theta); 189 } 190 191 /** 192 * Get the rotation matrix corresponding to a rotation about this axis 193 * @param theta The amount to rotate 194 * @return A 3x3 rotation matrix 195 */ 196 public Matrix getRotationMatrix(double theta) { 197 if( rotationAxis == null) { 198 // special case for pure translational axes 199 return Matrix.identity(3, 3); 200 } 201 double x = rotationAxis.getX(); 202 double y = rotationAxis.getY(); 203 double z = rotationAxis.getZ(); 204 double cos = Math.cos(theta); 205 double sin = Math.sin(theta); 206 double com = 1 - cos; 207 return new Matrix(new double[][] { 208 {com*x*x + cos, com*x*y+sin*z, com*x*z+-sin*y}, 209 {com*x*y-sin*z, com*y*y+cos, com*y*z+sin*x}, 210 {com*x*z+sin*y, com*y*z-sin*x, com*z*z+cos}, 211 }); 212 } 213 214 /** 215 * Returns the rotation order o that gives the lowest value of {@code |2PI / o - theta}, 216 * given that the value is strictly lower than {@code threshold}, for orders {@code o=1,...,maxOrder}. 217 */ 218 public int guessOrderFromAngle(double threshold, int maxOrder) { 219 double bestDelta = threshold; 220 int bestOrder = 1; 221 for (int order = 2; order < maxOrder; order++) { 222 double delta = Math.abs(2 * Math.PI / order - theta); 223 if (delta < bestDelta) { 224 bestOrder = order; 225 bestDelta = delta; 226 } 227 } 228 return bestOrder; 229 } 230 231 232 /** 233 * Returns a matrix that describes both rotation and translation. 234 */ 235 public Matrix getFullMatrix() { 236 return null; // TODO, easy 237 } 238 239 /** 240 * Initialize variables 241 * 242 * @param rotation 243 * @param translation 244 */ 245 private void init(Matrix rotation, Atom translation) { 246 if(rotation.getColumnDimension() != 3 || rotation.getRowDimension() != 3) { 247 throw new IllegalArgumentException("Expected 3x3 rotation matrix"); 248 } 249 250 251 // Calculate angle 252 double c = (rotation.trace()-1)/2.0; //=cos(theta) 253 // c is sometimes slightly out of the [-1,1] range due to numerical instabilities 254 if( -1-1e-8 < c && c < -1 ) c = -1; 255 if( 1+1e-8 > c && c > 1 ) c = 1; 256 if( -1 > c || c > 1 ) { 257 throw new IllegalArgumentException("Input matrix is not a valid rotation matrix."); 258 } 259 this.theta = Math.acos(c); 260 261 if(theta < MIN_ANGLE) { 262 calculateTranslationalAxis(rotation,translation); 263 } else { 264 calculateRotationalAxis(rotation, translation, c); 265 } 266 } 267 268 /** 269 * Calculate the rotation axis for the normal case, where there is a 270 * significant rotation angle 271 * @param rotation 272 * @param translation 273 * @param c 274 */ 275 private void calculateRotationalAxis(Matrix rotation, Atom translation, 276 double c) { 277 // Calculate magnitude of rotationAxis components, but not signs 278 double sum=0; 279 double[] rotAx = new double[3]; 280 for(int i=0;i<3;i++) { 281 rotAx[i] = Math.sqrt(rotation.get(i, i)-c); 282 sum+=rotAx[i]*rotAx[i]; 283 } 284 for(int i=0;i<3;i++) { 285 rotAx[i] /= Math.sqrt(sum); 286 } 287 288 // Now determine signs 289 double d0 = rotation.get(2,1)-rotation.get(1,2); //=2u[0]*sin(theta) 290 double d1 = rotation.get(0,2)-rotation.get(2,0); //=2u[1]*sin(theta) 291 double d2 = rotation.get(1,0)-rotation.get(0,1); //=2u[2]*sin(theta) 292 293 double s12 = rotation.get(2,1)+rotation.get(1,2); //=2*u[1]*u[2]*(1-cos(theta)) 294 double s02 = rotation.get(0,2)+rotation.get(2,0); //=2*u[0]*u[2]*(1-cos(theta)) 295 double s01 = rotation.get(1,0)+rotation.get(0,1); //=2*u[0]*u[1]*(1-cos(theta)) 296 297 //Only use biggest d for the sign directly, for numerical stability around 180deg 298 if( Math.abs(d0) < Math.abs(d1) ) { // not d0 299 if( Math.abs(d1) < Math.abs(d2) ) { //d2 300 if(d2>=0){ //u[2] positive 301 if( s02 < 0 ) rotAx[0] = -rotAx[0]; 302 if( s12 < 0 ) rotAx[1] = -rotAx[1]; 303 } else { //u[2] negative 304 rotAx[2] = -rotAx[2]; 305 if( s02 >= 0 ) rotAx[0] = -rotAx[0]; 306 if( s12 >= 0 ) rotAx[1] = -rotAx[1]; 307 } 308 } else { //d1 309 if(d1>=0) {//u[1] positive 310 if( s01 < 0) rotAx[0] = -rotAx[0]; 311 if( s12 < 0) rotAx[2] = -rotAx[2]; 312 } else { //u[1] negative 313 rotAx[1] = -rotAx[1]; 314 if( s01 >= 0) rotAx[0] = -rotAx[0]; 315 if( s12 >= 0) rotAx[2] = -rotAx[2]; 316 } 317 } 318 } else { // not d1 319 if( Math.abs(d0) < Math.abs(d2) ) { //d2 320 if(d2>=0){ //u[2] positive 321 if( s02 < 0 ) rotAx[0] = -rotAx[0]; 322 if( s12 < 0 ) rotAx[1] = -rotAx[1]; 323 } else { //u[2] negative 324 rotAx[2] = -rotAx[2]; 325 if( s02 >= 0 ) rotAx[0] = -rotAx[0]; 326 if( s12 >= 0 ) rotAx[1] = -rotAx[1]; 327 } 328 } else { //d0 329 if(d0>=0) { //u[0] positive 330 if( s01 < 0 ) rotAx[1] = -rotAx[1]; 331 if( s02 < 0 ) rotAx[2] = -rotAx[2]; 332 } else { //u[0] negative 333 rotAx[0] = -rotAx[0]; 334 if( s01 >= 0 ) rotAx[1] = -rotAx[1]; 335 if( s02 >= 0 ) rotAx[2] = -rotAx[2]; 336 } 337 } 338 } 339 340 rotationAxis = new AtomImpl(); 341 rotationAxis.setCoords(rotAx); 342 343 // Calculate screw = (rotationAxis dot translation)*u 344 double dotProduct = Calc.scalarProduct(rotationAxis, translation); 345 346 screwTranslation = Calc.scale(rotationAxis, dotProduct); 347 otherTranslation = Calc.subtract(translation, screwTranslation); 348 349 Atom hypot = Calc.vectorProduct(otherTranslation,rotationAxis); 350 Calc.scaleEquals(hypot,.5/Math.tan(theta/2.0)); 351 352 // Calculate rotation axis position 353 rotationPos = Calc.scaleAdd(.5,otherTranslation, hypot); 354 } 355 356 /** 357 * Handle cases with small angles of rotation 358 * @param rotation 359 * @param translation 360 */ 361 private void calculateTranslationalAxis(Matrix rotation, Atom translation) { 362 // set axis parallel to translation 363 rotationAxis = Calc.scale(translation, 1./Calc.amount(translation)); 364 365 // position is undefined 366 rotationPos = null; 367 368 screwTranslation = translation; 369 otherTranslation = new AtomImpl(); 370 otherTranslation.setCoords(new double[] {0,0,0}); 371 } 372 373 /** 374 * Returns a Jmol script which will display the axis of rotation. This 375 * consists of a cyan arrow along the axis, plus an arc showing the angle 376 * of rotation. 377 * <p> 378 * As the rotation angle gets smaller, the axis of rotation becomes poorly 379 * defined and would need to get farther and farther away from the protein. 380 * This is not particularly useful, so we arbitrarily draw it parallel to 381 * the translation and omit the arc. 382 * @param atoms Some atoms from the protein, used for determining the bounds 383 * of the axis. 384 * 385 * @return The Jmol script, suitable for calls to 386 * {@link org.biojava.nbio.structure.align.gui.jmol.StructureAlignmentJmol#evalString() jmol.evalString()} 387 */ 388 public String getJmolScript(Atom[] atoms){ 389 return getJmolScript(atoms, 0); 390 } 391 392 /** 393 * Find a segment of the axis that covers the specified set of atoms. 394 * <p> 395 * Projects the input atoms onto the rotation axis and returns the bounding 396 * points. 397 * <p> 398 * In the case of a pure translational axis, the axis location is undefined 399 * so the center of mass will be used instead. 400 * @param atoms 401 * @return two points defining the axis segment 402 */ 403 public Pair<Atom> getAxisEnds(Atom[] atoms) { 404 // Project each Atom onto the rotation axis to determine limits 405 double min, max; 406 min = max = Calc.scalarProduct(rotationAxis,atoms[0]); 407 for(int i=1;i<atoms.length;i++) { 408 double prod = Calc.scalarProduct(rotationAxis,atoms[i]); 409 if(prod<min) min = prod; 410 if(prod>max) max = prod; 411 } 412 double uLen = Calc.scalarProduct(rotationAxis,rotationAxis);// Should be 1, but double check 413 min/=uLen; 414 max/=uLen; 415 416 // Project the origin onto the axis. If the axis is undefined, use the center of mass 417 Atom axialPt; 418 if(rotationPos == null) { 419 Atom center = Calc.centerOfMass(atoms); 420 421 // Project center onto the axis 422 Atom centerOnAxis = Calc.scale(rotationAxis, Calc.scalarProduct(center, rotationAxis)); 423 424 // Remainder is projection of origin onto axis 425 axialPt = Calc.subtract(center, centerOnAxis); 426 427 } else { 428 axialPt = rotationPos; 429 } 430 431 // Find end points of the rotation axis to display 432 Atom axisMin = (Atom) axialPt.clone(); 433 Calc.scaleAdd(min, rotationAxis, axisMin); 434 Atom axisMax = (Atom) axialPt.clone(); 435 Calc.scaleAdd(max, rotationAxis, axisMax); 436 437 return new Pair<>(axisMin, axisMax); 438 } 439 /** 440 * Returns a Jmol script which will display the axis of rotation. This 441 * consists of a cyan arrow along the axis, plus an arc showing the angle 442 * of rotation. 443 * <p> 444 * As the rotation angle gets smaller, the axis of rotation becomes poorly 445 * defined and would need to get farther and farther away from the protein. 446 * This is not particularly useful, so we arbitrarily draw it parallel to 447 * the translation and omit the arc. 448 * @param atoms Some atoms from the protein, used for determining the bounds 449 * of the axis. 450 * @param axisID in case of representing more than one axis in the same jmol 451 * panel, indicate the ID number. 452 * 453 * @return The Jmol script, suitable for calls to 454 * {@link org.biojava.nbio.structure.align.gui.jmol.StructureAlignmentJmol#evalString() jmol.evalString()} 455 */ 456 public String getJmolScript(Atom[] atoms, int axisID){ 457 final double width=.5;// width of JMol object 458 final String axisColor = "yellow"; //axis color 459 final String screwColor = "orange"; //screw translation color 460 461 Pair<Atom> endPoints = getAxisEnds(atoms); 462 Atom axisMin = endPoints.getFirst(); 463 Atom axisMax = endPoints.getSecond(); 464 465 StringWriter result = new StringWriter(); 466 467 // set arrow heads to a reasonable length 468 result.append("set defaultDrawArrowScale 2.0;"); 469 470 // draw axis of rotation 471 result.append( 472 String.format("draw ID rot"+axisID+" CYLINDER {%f,%f,%f} {%f,%f,%f} WIDTH %f COLOR %s ;", 473 axisMin.getX(),axisMin.getY(),axisMin.getZ(), 474 axisMax.getX(),axisMax.getY(),axisMax.getZ(), width, axisColor )); 475 476 // draw screw component 477 boolean positiveScrew = Math.signum(rotationAxis.getX()) == Math.signum(screwTranslation.getX()); 478 if( positiveScrew ) { 479 // screw is in the same direction as the axis 480 result.append( String.format( 481 "draw ID screw"+axisID+" VECTOR {%f,%f,%f} {%f,%f,%f} WIDTH %f COLOR %s ;", 482 axisMax.getX(),axisMax.getY(),axisMax.getZ(), 483 screwTranslation.getX(),screwTranslation.getY(),screwTranslation.getZ(), 484 width, screwColor )); 485 } else { 486 // screw is in the opposite direction as the axis 487 result.append( String.format( 488 "draw ID screw"+axisID+" VECTOR {%f,%f,%f} {%f,%f,%f} WIDTH %f COLOR %s ;", 489 axisMin.getX(),axisMin.getY(),axisMin.getZ(), 490 screwTranslation.getX(),screwTranslation.getY(),screwTranslation.getZ(), 491 width, screwColor )); 492 } 493 494 // draw angle of rotation 495 if(rotationPos != null) { 496 result.append(System.getProperty("line.separator")); 497 result.append(String.format("draw ID rotArc"+axisID+" ARC {%f,%f,%f} {%f,%f,%f} {0,0,0} {0,%f,%d} SCALE 500 DIAMETER %f COLOR %s;", 498 axisMin.getX(),axisMin.getY(),axisMin.getZ(), 499 axisMax.getX(),axisMax.getY(),axisMax.getZ(), 500 Math.toDegrees(theta), 501 positiveScrew ? 0 : 1 , // draw at the opposite end from the screw arrow 502 width, axisColor )); 503 } 504 505 return result.toString(); 506 } 507 508 /** 509 * Projects a given point onto the axis of rotation 510 * @param point 511 * @return An atom which lies on the axis, or null if the RotationAxis is purely translational 512 */ 513 public Atom getProjectedPoint(Atom point) { 514 if(rotationPos == null) { 515 // translation only 516 return null; 517 } 518 519 Atom localPoint = Calc.subtract(point, rotationPos); 520 double dot = Calc.scalarProduct(localPoint, rotationAxis); 521 522 Atom localProjected = Calc.scale(rotationAxis, dot); 523 Atom projected = Calc.add(localProjected, rotationPos); 524 return projected; 525 } 526 527 /** 528 * Get the distance from a point to the axis of rotation 529 * @param point 530 * @return The distance to the axis, or NaN if the RotationAxis is purely translational 531 */ 532 public double getProjectedDistance(Atom point) { 533 Atom projected = getProjectedPoint(point); 534 if( projected == null) { 535 // translation only 536 return Double.NaN; 537 } 538 539 540 return Calc.getDistance(point, projected); 541 542 } 543 544 public void rotate(Atom[] atoms, double theta) { 545 Matrix rot = getRotationMatrix(theta); 546 if(rotationPos == null) { 547 // Undefined rotation axis; do nothing 548 return; 549 } 550 Atom negPos; 551 552 negPos = Calc.invert(rotationPos); 553 554 for(Atom a: atoms) { 555 Calc.shift(a, negPos); 556 } 557 Calc.rotate(atoms, rot); 558 for(Atom a: atoms) { 559 Calc.shift(a, rotationPos); 560 } 561 } 562 563 /** 564 * Calculate the rotation angle for a structure 565 * @param afpChain 566 * @return The rotation angle, in radians 567 * @throws StructureException If the alignment doesn't contain any blocks 568 * @throws NullPointerException If the alignment doesn't have a rotation matrix set 569 */ 570 public static double getAngle(AFPChain afpChain) throws StructureException { 571 if(afpChain.getBlockNum() < 1) { 572 throw new StructureException("No aligned residues"); 573 } 574 Matrix rotation = afpChain.getBlockRotationMatrix()[0]; 575 576 if(rotation == null) { 577 throw new NullPointerException("AFPChain does not contain a rotation matrix"); 578 } 579 return getAngle(rotation); 580 } 581 /** 582 * Calculate the rotation angle for a given matrix 583 * @param rotation Rotation matrix 584 * @return The angle, in radians 585 */ 586 public static double getAngle(Matrix rotation) { 587 double c = (rotation.trace()-1)/2.0; //=cos(theta) 588 return Math.acos(c); 589 } 590 591 /** 592 * 593 * @return If the rotation axis is well defined, rather than purely translational 594 */ 595 public boolean isDefined() { 596 return rotationPos != null; 597 } 598}