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.symmetry.axis; 022 023import org.biojava.nbio.structure.geometry.CalcPoint; 024import org.biojava.nbio.structure.geometry.MomentsOfInertia; 025import org.biojava.nbio.structure.symmetry.core.QuatSymmetryResults; 026import org.biojava.nbio.structure.symmetry.core.Rotation; 027import org.biojava.nbio.structure.symmetry.core.RotationGroup; 028import org.biojava.nbio.structure.symmetry.core.QuatSymmetrySubunits; 029import org.slf4j.Logger; 030import org.slf4j.LoggerFactory; 031 032import javax.vecmath.*; 033 034import java.util.*; 035 036public class RotationAxisAligner extends AxisAligner{ 037 038 private static final Logger logger = LoggerFactory 039 .getLogger(RotationAxisAligner.class); 040 041 private static final Vector3d X_AXIS = new Vector3d(1,0,0); 042 private static final Vector3d Y_AXIS = new Vector3d(0,1,0); 043 private static final Vector3d Z_AXIS = new Vector3d(0,0,1); 044 045 private QuatSymmetrySubunits subunits = null; 046 private RotationGroup rotationGroup = null; 047 048 private Matrix4d transformationMatrix = new Matrix4d(); 049 private Matrix4d reverseTransformationMatrix = new Matrix4d(); 050 private Vector3d referenceVector = new Vector3d(); 051 private Vector3d principalRotationVector = new Vector3d(); 052 private Vector3d[] principalAxesOfInertia = null; 053 List<List<Integer>> alignedOrbits = null; 054 055 private Vector3d minBoundary = new Vector3d(); 056 private Vector3d maxBoundary = new Vector3d(); 057 private double xyRadiusMax = Double.MIN_VALUE; 058 059 boolean modified = true; 060 061 public RotationAxisAligner(QuatSymmetryResults results) { 062 this.subunits = new QuatSymmetrySubunits(results.getSubunitClusters()); 063 this.rotationGroup = results.getRotationGroup(); 064 065 if (subunits == null) { 066 throw new IllegalArgumentException("AxisTransformation: Subunits are null"); 067 } else if (rotationGroup == null) { 068 throw new IllegalArgumentException("AxisTransformation: RotationGroup is null"); 069 } else if (subunits.getSubunitCount() == 0) { 070 throw new IllegalArgumentException("AxisTransformation: Subunits is empty"); 071 } else if (rotationGroup.getOrder() == 0) { 072 throw new IllegalArgumentException("AxisTransformation: RotationGroup is empty"); 073 } 074 } 075 076 /* (non-Javadoc) 077 * @see org.biojava.nbio.structure.quaternary.core.AxisAligner#getTransformation() 078 */ 079 @Override 080 public String getSymmetry() { 081 run(); 082 return rotationGroup.getPointGroup(); 083 } 084 085 @Override 086 public Matrix4d getTransformation() { 087 run(); 088 return transformationMatrix; 089 } 090 091 @Override 092 public Matrix3d getRotationMatrix() { 093 run(); 094 Matrix3d m = new Matrix3d(); 095 transformationMatrix.getRotationScale(m); 096 return m; 097 } 098 099 @Override 100 public Matrix4d getReverseTransformation() { 101 run(); 102 return reverseTransformationMatrix; 103 } 104 105 @Override 106 public Vector3d getPrincipalRotationAxis() { 107 run(); 108 return principalRotationVector; 109 } 110 111 @Override 112 public Vector3d getRotationReferenceAxis() { 113 run(); 114 return referenceVector; 115 } 116 117 @Override 118 public Vector3d[] getPrincipalAxesOfInertia() { 119 run(); 120 return principalAxesOfInertia; 121 } 122 123 @Override 124 public Vector3d getDimension() { 125 run(); 126 Vector3d dimension = new Vector3d(); 127 dimension.sub(maxBoundary, minBoundary); 128 dimension.scale(0.5); 129 return dimension; 130 } 131 132 /** 133 * Returns the radius for drawing the minor rotation axis in the xy-plane 134 * @return double radius in xy-plane 135 */ 136 @Override 137 public double getRadius() { 138 run(); 139 return xyRadiusMax; 140 } 141 142 /** 143 * Returns a transformation matrix transform polyhedra for Cn structures. 144 * The center in this matrix is the geometric center, rather then the centroid. 145 * In Cn structures those are usually not the same. 146 * @return 147 */ 148 @Override 149 public Matrix4d getGeometicCenterTransformation() { 150 run(); 151 152 Matrix4d geometricCentered = new Matrix4d(reverseTransformationMatrix); 153 geometricCentered.setTranslation(new Vector3d(getGeometricCenter())); 154 155 return geometricCentered; 156 } 157 158 /** 159 * Returns the geometric center of polyhedron. In the case of the Cn 160 * point group, the centroid and geometric center are usually not 161 * identical. 162 * @return 163 */ 164 @Override 165 public Point3d getGeometricCenter() { 166 run(); 167 168 Point3d geometricCenter = new Point3d(); 169 Vector3d translation = new Vector3d(); 170 reverseTransformationMatrix.get(translation); 171 172 // calculate adjustment around z-axis and transform adjustment to 173 // original coordinate frame with the reverse transformation 174 if (rotationGroup.getPointGroup().startsWith("C")) { 175 Vector3d corr = new Vector3d(0,0, minBoundary.z+getDimension().z); 176 reverseTransformationMatrix.transform(corr); 177 geometricCenter.set(corr); 178 } 179 180 geometricCenter.add(translation); 181 return geometricCenter; 182 } 183 184 @Override 185 public Point3d getCentroid() { 186 return new Point3d(subunits.getCentroid()); 187 } 188 189 @Override 190 public QuatSymmetrySubunits getSubunits() { 191 return subunits; 192 } 193 194 public RotationGroup getRotationGroup() { 195 return rotationGroup; 196 } 197 198 @Override 199 public List<List<Integer>> getOrbits() { 200 return alignedOrbits; 201 } 202 203 /** 204 * @return 205 */ 206 207 private void run () { 208 if (modified) { 209 calcPrincipalRotationVector(); 210 calcPrincipalAxes(); 211 // initial alignment with draft reference axis 212 calcReferenceVector(); 213 calcTransformation(); 214 215 // refine ref. axis for cyclic and dihedral systems 216 if ((rotationGroup.getPointGroup().startsWith("C") && 217 !rotationGroup.getPointGroup().startsWith("C2")) || 218 (rotationGroup.getPointGroup().startsWith("D") && 219 !rotationGroup.getPointGroup().startsWith("D2")) 220 ) { 221 refineReferenceVector(); 222 calcTransformation(); 223 } 224 calcReverseTransformation(); 225 calcBoundaries(); 226 if (! rotationGroup.getPointGroup().equals("Helical")) { 227 calcAlignedOrbits(); 228 } 229 modified = false; 230 } 231 } 232 /** 233 * Returns a list of orbits (an orbit is a cyclic permutation of subunit indices that are related 234 * by a rotation around the principal rotation axis) ordered from the +z direction to the -z direction (z-depth). 235 * Within an orbit, subunit indices are ordered with a clockwise rotation around the z-axis. 236 * @return list of orbits ordered by z-depth 237 */ 238 private void calcAlignedOrbits() { 239 Map<Double, List<Integer>> depthMap = new TreeMap<Double, List<Integer>>(); 240 double[] depth = getSubunitZDepth(); 241 alignedOrbits = calcOrbits(); 242 243 // calculate the mean depth of orbit along z-axis 244 for (List<Integer> orbit: alignedOrbits) { 245 // calculate the mean depth along z-axis for each orbit 246 double meanDepth = 0; 247 for (int subunit: orbit) { 248 meanDepth += depth[subunit]; 249 } 250 meanDepth /= orbit.size(); 251 252 if (depthMap.get(meanDepth) != null) { 253 // System.out.println("Conflict in depthMap"); 254 meanDepth += 0.01; 255 } 256 depthMap.put(meanDepth, orbit); 257 } 258 259 // now fill orbits back into list ordered by depth 260 alignedOrbits.clear(); 261 for (List<Integer> orbit: depthMap.values()) { 262 // order subunit in a clockwise rotation around the z-axis 263 /// starting at the 12 O-clock position (+y position) 264 alignWithReferenceAxis(orbit); 265 alignedOrbits.add(orbit); 266 } 267 } 268 269 /** 270 * Returns an ordered list of subunit ids (orbit) in such a way that the subunit 271 * indices start at the 12 o-clock (+y axis) and proceed in a clockwise direction 272 * to the 11 o-clock position to close the "orbit". 273 * 274 * @param orbit list of subunit indices that are transformed into each other by a rotation 275 * @return list of subunit indices ordered in a clockwise direction 276 */ 277 278 private List<Integer> alignWithReferenceAxis(List<Integer> orbit) { 279 int n = subunits.getSubunitCount(); 280 int fold = rotationGroup.getRotation(0).getFold(); 281 if (fold < 2) { 282 return orbit; 283 } 284 Vector3d probe = new Vector3d(); 285 286 double dotMin1 = Double.MIN_VALUE; 287 double dotMin2 = Double.MIN_VALUE; 288 int index1 = 0; 289 int index2 = 0; 290 Vector3d Y1 = new Vector3d(0,1,0); 291 Vector3d Y2 = new Vector3d(0,1,0); 292 Matrix3d m = new Matrix3d(); 293 double angle = -2*Math.PI/fold; 294 m.rotZ(0.1*angle); // add small offset, since two subunits may be equidistant to the y-axis 295 m.transform(Y1); 296 m.rotZ(1.1*angle); 297 m.transform(Y2); 298 // transform subunit centers into z-aligned position and calculate 299 // width in xy direction. 300 for (int i: orbit) { 301 Point3d p = subunits.getCenters().get(i); 302 probe.set(p); 303 transformationMatrix.transform(probe); 304 // find subunit that lines up with y-axis 305 double dot1 = Y1.dot(probe); 306 if (dot1 > dotMin1) { 307 dotMin1 = dot1; 308 index1 = i; 309 } 310 // find next subunit (rotated by one fold around z-axis - clockwise) 311 double dot2 = Y2.dot(probe); 312 if (dot2 > dotMin2) { 313 dotMin2 = dot2; 314 index2 = i; 315 } 316 } 317// System.out.println("Index1/2: " + index1 + " - " + index2); 318// System.out.println("Orbit0: " + orbit); 319 // order subunit indices in a clockwise orientation around the z-axis 320 // bring subunit into position 0 321 for (int i = 0; i < n; i++) { 322 if (orbit.get(0) == index1) { 323 break; 324 } 325 Collections.rotate(orbit,1); 326 } 327// System.out.println("Orbit1: " + orbit); 328 // bring second subunit onto position 1 329 if (orbit.get(1) == index2) { 330 return orbit; 331 } 332 Collections.reverse(orbit.subList(1, orbit.size())); 333 if (orbit.get(1) != index2) { 334 logger.warn("alignWithReferenceAxis failed"); 335 } 336// System.out.println("Orbit2: " + orbit); 337 return orbit; 338 } 339 340 341 342 private void calcTransformation() { 343 if (rotationGroup.getPointGroup().equals("C1")) { 344 calcTransformationByInertiaAxes(); 345 } else { 346 calcTransformationBySymmetryAxes(); 347 } 348 // make sure this value is zero. On Linux this value is 0. 349 transformationMatrix.setElement(3, 3, 1.0); 350 } 351 352 private void calcReverseTransformation() { 353 reverseTransformationMatrix.invert(transformationMatrix); 354 } 355 356 private void calcTransformationBySymmetryAxes() { 357 Vector3d[] axisVectors = new Vector3d[2]; 358 axisVectors[0] = new Vector3d(principalRotationVector); 359 axisVectors[1] = new Vector3d(referenceVector); 360 361 // y,z axis centered at the centroid of the subunits 362 Vector3d[] referenceVectors = new Vector3d[2]; 363 referenceVectors[0] = new Vector3d(Z_AXIS); 364 referenceVectors[1] = new Vector3d(Y_AXIS); 365 366 transformationMatrix = alignAxes(axisVectors, referenceVectors); 367 368 // combine with translation 369 Matrix4d combined = new Matrix4d(); 370 combined.setIdentity(); 371 Vector3d trans = new Vector3d(subunits.getCentroid()); 372 trans.negate(); 373 combined.setTranslation(trans); 374 transformationMatrix.mul(combined); 375 376 // for cyclic geometry, set a canonical view for the Z direction 377 if (rotationGroup.getPointGroup().startsWith("C")) { 378 calcZDirection(); 379 } 380 } 381 382 private void calcTransformationByInertiaAxes() { 383 Vector3d[] axisVectors = new Vector3d[2]; 384 axisVectors[0] = new Vector3d(principalAxesOfInertia[0]); 385 axisVectors[1] = new Vector3d(principalAxesOfInertia[1]); 386 387 Vector3d[] referenceVectors = new Vector3d[2]; 388 referenceVectors[0] = new Vector3d(Y_AXIS); 389 referenceVectors[1] = new Vector3d(X_AXIS); 390 391 // align inertia axes with y-x plane 392 transformationMatrix = alignAxes(axisVectors, referenceVectors); 393 394 // combine with translation 395 Matrix4d translation = new Matrix4d(); 396 translation.setIdentity(); 397 Vector3d trans = new Vector3d(subunits.getCentroid()); 398 trans.negate(); 399 translation.setTranslation(trans); 400 transformationMatrix.mul(translation); 401 } 402 403 /** 404 * Returns a transformation matrix that rotates refPoints to match 405 * coordPoints 406 * @param refPoints the points to be aligned 407 * @param referenceVectors 408 * @return 409 */ 410 private Matrix4d alignAxes(Vector3d[] axisVectors, Vector3d[] referenceVectors) { 411 Matrix4d m1 = new Matrix4d(); 412 AxisAngle4d a = new AxisAngle4d(); 413 Vector3d axis = new Vector3d(); 414 415 // calculate rotation matrix to rotate refPoints[0] into coordPoints[0] 416 Vector3d v1 = new Vector3d(axisVectors[0]); 417 Vector3d v2 = new Vector3d(referenceVectors[0]); 418 double dot = v1.dot(v2); 419 if (Math.abs(dot) < 0.999) { 420 axis.cross(v1,v2); 421 axis.normalize(); 422 a.set(axis, v1.angle(v2)); 423 m1.set(a); 424 // make sure matrix element m33 is 1.0. It's 0 on Linux. 425 m1.setElement(3, 3, 1.0); 426 } else if (dot > 0) { 427 // parallel axis, nothing to do -> identity matrix 428 m1.setIdentity(); 429 } else if (dot < 0) { 430 // anti-parallel axis, flip around x-axis 431 m1.set(flipX()); 432 } 433 434 // apply transformation matrix to all refPoints 435 m1.transform(axisVectors[0]); 436 m1.transform(axisVectors[1]); 437 438 // calculate rotation matrix to rotate refPoints[1] into coordPoints[1] 439 v1 = new Vector3d(axisVectors[1]); 440 v2 = new Vector3d(referenceVectors[1]); 441 Matrix4d m2 = new Matrix4d(); 442 dot = v1.dot(v2); 443 if (Math.abs(dot) < 0.999) { 444 axis.cross(v1,v2); 445 axis.normalize(); 446 a.set(axis, v1.angle(v2)); 447 m2.set(a); 448 // make sure matrix element m33 is 1.0. It's 0 on Linux. 449 m2.setElement(3, 3, 1.0); 450 } else if (dot > 0) { 451 // parallel axis, nothing to do -> identity matrix 452 m2.setIdentity(); 453 } else if (dot < 0) { 454 // anti-parallel axis, flip around z-axis 455 m2.set(flipZ()); 456 } 457 458 // apply transformation matrix to all refPoints 459 m2.transform(axisVectors[0]); 460 m2.transform(axisVectors[1]); 461 462 // combine the two rotation matrices 463 m2.mul(m1); 464 465 // the RMSD should be close to zero 466 Point3d[] axes = new Point3d[2]; 467 axes[0] = new Point3d(axisVectors[0]); 468 axes[1] = new Point3d(axisVectors[1]); 469 Point3d[] ref = new Point3d[2]; 470 ref[0] = new Point3d(referenceVectors[0]); 471 ref[1] = new Point3d(referenceVectors[1]); 472 if (CalcPoint.rmsd(axes, ref) > 0.1) { 473 logger.warn("AxisTransformation: axes alignment is off. RMSD: " 474 + CalcPoint.rmsd(axes, ref)); 475 } 476 477 return m2; 478 } 479 480 private void calcPrincipalAxes() { 481 MomentsOfInertia moi = new MomentsOfInertia(); 482 483 for (Point3d[] list: subunits.getTraces()) { 484 for (Point3d p: list) { 485 moi.addPoint(p, 1.0); 486 } 487 } 488 principalAxesOfInertia = moi.getPrincipalAxes(); 489 } 490 491 /** 492 * Calculates the min and max boundaries of the structure after it has been 493 * transformed into its canonical orientation. 494 */ 495 private void calcBoundaries() { 496 minBoundary.x = Double.MAX_VALUE; 497 maxBoundary.x = Double.MIN_VALUE; 498 minBoundary.y = Double.MAX_VALUE; 499 maxBoundary.x = Double.MIN_VALUE; 500 minBoundary.z = Double.MAX_VALUE; 501 maxBoundary.z = Double.MIN_VALUE; 502 503 Point3d probe = new Point3d(); 504 505 for (Point3d[] list: subunits.getTraces()) { 506 for (Point3d p: list) { 507 probe.set(p); 508 transformationMatrix.transform(probe); 509 510 minBoundary.x = Math.min(minBoundary.x, probe.x); 511 maxBoundary.x = Math.max(maxBoundary.x, probe.x); 512 minBoundary.y = Math.min(minBoundary.y, probe.y); 513 maxBoundary.y = Math.max(maxBoundary.y, probe.y); 514 minBoundary.z = Math.min(minBoundary.z, probe.z); 515 maxBoundary.z = Math.max(maxBoundary.z, probe.z); 516 xyRadiusMax = Math.max(xyRadiusMax, Math.sqrt(probe.x*probe.x + probe.y * probe.y)); 517 } 518 } 519 } 520 521 /* 522 * Modifies the rotation part of the transformation axis for 523 * a Cn symmetric complex, so that the narrower end faces the 524 * viewer, and the wider end faces away from the viewer. Example: 3LSV 525 */ 526 private void calcZDirection() { 527 calcBoundaries(); 528 529 // if the longer part of the structure faces towards the back (-z direction), 530 // rotate around y-axis so the longer part faces the viewer (+z direction) 531 if (Math.abs(minBoundary.z) > Math.abs(maxBoundary.z)) { 532 Matrix4d rot = flipY(); 533 rot.mul(transformationMatrix); 534 transformationMatrix.set(rot); 535 } 536 } 537 538 /** 539 * 540 */ 541 private List<List<Integer>> getOrbitsByXYWidth() { 542 Map<Double, List<Integer>> widthMap = new TreeMap<Double, List<Integer>>(); 543 double[] width = getSubunitXYWidth(); 544 List<List<Integer>> orbits = calcOrbits(); 545 546 // calculate the mean width of orbits in XY-plane 547 for (List<Integer> orbit: orbits) { 548 double meanWidth = 0; 549 for (int subunit: orbit) { 550 meanWidth += width[subunit]; 551 } 552 meanWidth /= orbit.size(); 553 554 if (widthMap.get(meanWidth) != null) { 555 meanWidth += 0.01; 556 } 557 widthMap.put(meanWidth, orbit); 558 } 559 560 // now fill orbits back into list ordered by width 561 orbits.clear(); 562 for (List<Integer> orbit: widthMap.values()) { 563 orbits.add(orbit); 564 } 565 return orbits; 566 } 567 568 private double[] getSubunitXYWidth() { 569 int n = subunits.getSubunitCount(); 570 double[] width = new double[n]; 571 Point3d probe = new Point3d(); 572 573 // transform subunit centers into z-aligned position and calculate 574 // width in xy direction. 575 for (int i = 0; i < n; i++) { 576 width[i] = Double.MIN_VALUE; 577 for (Point3d p: subunits.getTraces().get(i)) { 578 probe.set(p); 579 transformationMatrix.transform(probe); 580 width[i] = Math.max(width[i], Math.sqrt(probe.x*probe.x + probe.y*probe.y)); 581 } 582 } 583 return width; 584 } 585 586 private double[] getSubunitZDepth() { 587 int n = subunits.getSubunitCount(); 588 double[] depth = new double[n]; 589 Point3d probe = new Point3d(); 590 591 // transform subunit centers into z-aligned position and calculate 592 // z-coordinates (depth) along the z-axis. 593 for (int i = 0; i < n; i++) { 594 Point3d p= subunits.getCenters().get(i); 595 probe.set(p); 596 transformationMatrix.transform(probe); 597 depth[i] = probe.z; 598 } 599 return depth; 600 } 601 602 /** 603 * Returns a list of list of subunit ids that form an "orbit", i.e. they 604 * are transformed into each other during a rotation around the principal symmetry axis (z-axis) 605 * @return 606 */ 607 private List<List<Integer>> calcOrbits() { 608 int n = subunits.getSubunitCount(); 609 int fold = rotationGroup.getRotation(0).getFold(); 610 611 List<List<Integer>> orbits = new ArrayList<List<Integer>>(); 612 boolean[] used = new boolean[n]; 613 Arrays.fill(used, false); 614 615 for (int i = 0; i < n; i++) { 616 if (! used[i]) { 617 // determine the equivalent subunits 618 List<Integer> orbit = new ArrayList<Integer>(fold); 619 for (int j = 0; j < fold; j++) { 620 List<Integer> permutation = rotationGroup.getRotation(j).getPermutation(); 621 orbit.add(permutation.get(i)); 622 used[permutation.get(i)] = true; 623 } 624 orbits.add(deconvolute(orbit)); 625 } 626 } 627 return orbits; 628 } 629 630 private List<Integer> deconvolute(List<Integer> orbit) { 631 if (rotationGroup.getOrder() < 2) { 632 return orbit; 633 } 634 List<Integer> p0 = rotationGroup.getRotation(0).getPermutation(); 635 List<Integer> p1 = rotationGroup.getRotation(1).getPermutation(); 636// System.out.println("deconvolute"); 637// System.out.println("Permutation0: " + p0); 638// System.out.println("Permutation1: " + p1); 639 640 List<Integer> inRotationOrder = new ArrayList<Integer>(orbit.size()); 641 inRotationOrder.add(orbit.get(0)); 642 for (int i = 1; i < orbit.size(); i++) { 643 int current = inRotationOrder.get(i-1); 644 int index = p0.indexOf(current); 645 int next = p1.get(index); 646 if (!orbit.contains(next)) { 647 logger.warn("deconvolute: inconsistency in permuation. Returning original order"); 648 return orbit; 649 } 650 inRotationOrder.add(next); 651 } 652// System.out.println("In order: " + inRotationOrder); 653 return inRotationOrder; 654 } 655 656 /** 657 * Returns a vector along the principal rotation axis for the 658 * alignment of structures along the z-axis 659 * @return principal rotation vector 660 */ 661 private void calcPrincipalRotationVector() { 662 Rotation rotation = rotationGroup.getRotation(0); // the rotation around the principal axis is the first rotation 663 AxisAngle4d axisAngle = rotation.getAxisAngle(); 664 principalRotationVector = new Vector3d(axisAngle.x, axisAngle.y, axisAngle.z); 665 } 666 667 /** 668 * Returns a vector perpendicular to the principal rotation vector 669 * for the alignment of structures in the xy-plane 670 * @return reference vector 671 */ 672 private void calcReferenceVector() { 673 referenceVector = null; 674 if (rotationGroup.getPointGroup().startsWith("C")) { 675 referenceVector = getReferenceAxisCylic(); 676 } else if (rotationGroup.getPointGroup().startsWith("D")) { 677 referenceVector = getReferenceAxisDihedral(); 678 } else if (rotationGroup.getPointGroup().equals("T")) { 679 referenceVector = getReferenceAxisTetrahedral(); 680 } else if (rotationGroup.getPointGroup().equals("O")) { 681 referenceVector = getReferenceAxisOctahedral(); 682 } else if (rotationGroup.getPointGroup().equals("I")) { 683 referenceVector = getReferenceAxisIcosahedral(); 684 } else if (rotationGroup.getPointGroup().equals("Helical")) { 685 // TODO what should the reference vector be?? 686 referenceVector = getReferenceAxisCylic(); 687 } 688 689 if (referenceVector == null) { 690 logger.warn("no reference vector found. Using y-axis."); 691 referenceVector = new Vector3d(Y_AXIS); 692 } 693 // make sure reference vector is perpendicular principal roation vector 694 referenceVector = orthogonalize(principalRotationVector, referenceVector); 695 } 696 697 /** 698 * Returns a normalized vector that represents a minor rotation axis, except 699 * for Cn, this represents an axis orthogonal to the principal axis. 700 * @return minor rotation axis 701 */ 702 private void refineReferenceVector() { 703 referenceVector = new Vector3d(Y_AXIS); 704 if (rotationGroup.getPointGroup().startsWith("C")) { 705 referenceVector = getReferenceAxisCylicWithSubunitAlignment(); 706 } else if (rotationGroup.getPointGroup().startsWith("D")) { 707 referenceVector = getReferenceAxisDihedralWithSubunitAlignment(); 708 } 709 710 referenceVector = orthogonalize(principalRotationVector, referenceVector); 711 } 712 713 private Vector3d orthogonalize(Vector3d vector1, Vector3d vector2) { 714 double dot = vector1.dot(vector2); 715 Vector3d ref = new Vector3d(vector2); 716// System.out.println("p.r: " + dot); 717// System.out.println("Orig refVector: " + referenceVector); 718 if (dot < 0) { 719 vector2.negate(); 720 } 721 vector2.cross(vector1, vector2); 722// System.out.println("Intermed. refVector: " + vector2); 723 vector2.normalize(); 724// referenceVector.cross(referenceVector, principalRotationVector); 725 vector2.cross(vector1, vector2); 726 vector2.normalize(); 727 if (ref.dot(vector2) < 0) { 728 vector2.negate(); 729 } 730// System.out.println("Mod. refVector: " + vector2); 731 return vector2; 732 } 733 /** 734 * Returns the default reference vector for the alignment of Cn structures 735 * @return 736 */ 737 private Vector3d getReferenceAxisCylic() { 738 if (rotationGroup.getPointGroup().equals("C2")) { 739 Vector3d vr = new Vector3d(subunits.getOriginalCenters().get(0)); 740 vr.sub(subunits.getCentroid()); 741 vr.normalize(); 742 return vr; 743 } 744 745 // get principal axis vector that is perpendicular to the principal 746 // rotation vector 747 Vector3d vmin = null; 748 double dotMin = 1.0; 749 for (Vector3d v: principalAxesOfInertia) { 750 if (Math.abs(principalRotationVector.dot(v)) < dotMin) { 751 dotMin = Math.abs(principalRotationVector.dot(v)); 752 vmin = new Vector3d(v); 753 } 754 } 755 if (principalRotationVector.dot(vmin) < 0) { 756 vmin.negate(); 757 } 758 759 return vmin; 760 } 761 762 763 /** 764 * Returns a reference vector for the alignment of Cn structures. 765 * @return reference vector 766 */ 767 private Vector3d getReferenceAxisCylicWithSubunitAlignment() { 768 if (rotationGroup.getPointGroup().equals("C2")) { 769 return referenceVector; 770 } 771 772 // find subunit that extends the most in the xy-plane 773 List<List<Integer>> orbits = getOrbitsByXYWidth(); 774 // get the last orbit which is the widest 775 List<Integer> widestOrbit = orbits.get(orbits.size()-1); 776 List<Point3d> centers = subunits.getCenters(); 777 int subunit = widestOrbit.get(0); 778 779 // calculate reference vector 780 Vector3d refAxis = new Vector3d(); 781 refAxis.sub(centers.get(subunit), subunits.getCentroid()); 782 refAxis.normalize(); 783 return refAxis; 784 } 785 786 /** 787 * 788 */ 789 private Vector3d getReferenceAxisDihedralWithSubunitAlignment() { 790 int maxFold = rotationGroup.getRotation(0).getFold(); 791 792 double minAngle = Double.MAX_VALUE; 793 Vector3d refVec = null; 794 795 Vector3d ref = getSubunitReferenceVector(); 796 797 for (int i = 0; i < rotationGroup.getOrder(); i++) { 798 if (rotationGroup.getRotation(i).getDirection() == 1 && 799 (rotationGroup.getRotation(i).getFold() < maxFold) || 800 rotationGroup.getPointGroup().equals("D2")) { 801 802 AxisAngle4d axisAngle = rotationGroup.getRotation(i).getAxisAngle(); 803 Vector3d v = new Vector3d(axisAngle.x, axisAngle.y, axisAngle.z); 804 v.normalize(); 805 806// System.out.println("Ref axis angle(+): " + Math.toDegrees(v.angle(ref))); 807 double angle = v.angle(ref); 808 if (angle < minAngle) { 809 minAngle = angle; 810 refVec = v; 811 } 812 Vector3d vn = new Vector3d(v); 813 vn.negate(); 814// System.out.println("Ref axis angle(-): " + Math.toDegrees(vn.angle(ref))); 815 angle = vn.angle(ref); 816 if (angle < minAngle) { 817 minAngle = angle; 818 refVec = vn; 819 } 820 } 821 } 822 refVec.normalize(); 823 return refVec; 824 } 825 826 /** 827 * 828 */ 829 private Vector3d getReferenceAxisDihedral() { 830 int maxFold = rotationGroup.getRotation(0).getFold(); 831 // one exception: D2 832 if (maxFold == 2) { 833 maxFold = 3; 834 } 835 // TODO how about D2, where minor axis = 2 = principal axis?? 836 for (int i = 0; i < rotationGroup.getOrder(); i++) { 837 if (rotationGroup.getRotation(i).getDirection() == 1 && rotationGroup.getRotation(i).getFold() < maxFold) { 838 AxisAngle4d axisAngle = rotationGroup.getRotation(i).getAxisAngle(); 839 Vector3d v = new Vector3d(axisAngle.x, axisAngle.y, axisAngle.z); 840 v.normalize(); 841 return v; 842 } 843 } 844 return null; 845 } 846 847 private Vector3d getReferenceAxisTetrahedral() { 848 for (int i = 0; i < rotationGroup.getOrder(); i++) { 849 AxisAngle4d axisAngle = rotationGroup.getRotation(i).getAxisAngle(); 850 Vector3d v = new Vector3d(axisAngle.x, axisAngle.y, axisAngle.z); 851 double d = v.dot(principalRotationVector); 852 if (rotationGroup.getRotation(i).getFold() == 3) { 853 // the dot product 0 is between to adjacent 3-fold axes 854 if (d > 0.3 && d < 0.9) { 855 return v; 856 } 857 } 858 } 859 return null; 860 } 861 862 private Vector3d getReferenceAxisOctahedral() { 863 for (int i = 0; i < rotationGroup.getOrder(); i++) { 864 AxisAngle4d axisAngle = rotationGroup.getRotation(i).getAxisAngle(); 865 Vector3d v = new Vector3d(axisAngle.x, axisAngle.y, axisAngle.z); 866 double d = v.dot(principalRotationVector); 867 if (rotationGroup.getRotation(i).getFold() == 4) { 868 // the dot product 0 is between to adjacent 4-fold axes 869 if (d > -0.1 && d < 0.1 ) { 870 return v; 871 } 872 } 873 } 874 return null; 875 } 876 877 private Vector3d getReferenceAxisIcosahedral() { 878 for (int i = 0; i < rotationGroup.getOrder(); i++) { 879 AxisAngle4d axisAngle = rotationGroup.getRotation(i).getAxisAngle(); 880 Vector3d v = new Vector3d(axisAngle.x, axisAngle.y, axisAngle.z); 881 double d = v.dot(principalRotationVector); 882 if (rotationGroup.getRotation(i).getFold() == 5) { 883 // the dot product of 0.447.. is between to adjacent 5-fold axes 884// if (d > 0.447 && d < 0.448) { 885 if (d > 0.4 && d < 0.5) { 886 return v; 887 } 888 } 889 } 890 return null; 891 } 892 893 private Vector3d getSubunitReferenceVector() { 894 int n = subunits.getSubunitCount(); 895 Point3d probe = new Point3d(); 896 897 // transform subunit centers into z-aligned position and calculate 898 // width in xy direction. 899 double maxWidthSq = 0; 900 Point3d ref = null; 901 for (int i = 0; i < n; i++) { 902 for (Point3d p: subunits.getTraces().get(i)) { 903 probe.set(p); 904 transformationMatrix.transform(probe); 905 double widthSq = probe.x*probe.x + probe.y*probe.y; 906 if (widthSq > maxWidthSq) { 907 maxWidthSq = widthSq; 908 ref = p; 909 } 910 } 911 } 912 // System.out.println("width: " + maxWidthSq); 913 Vector3d refVector = new Vector3d(); 914 refVector.sub(ref, subunits.getCentroid()); 915 refVector.normalize(); 916 return refVector; 917 } 918 919 private static Matrix4d flipX() { 920 Matrix4d rot = new Matrix4d(); 921 rot.m00 = 1; 922 rot.m11 = -1; 923 rot.m22 = -1; 924 rot.m33 = 1; 925 return rot; 926 } 927 928 private static Matrix4d flipY() { 929 Matrix4d rot = new Matrix4d(); 930 rot.m00 = -1; 931 rot.m11 = 1; 932 rot.m22 = -1; 933 rot.m33 = 1; 934 return rot; 935 } 936 937 private static Matrix4d flipZ() { 938 Matrix4d rot = new Matrix4d(); 939 rot.m00 = -1; 940 rot.m11 = -1; 941 rot.m22 = 1; 942 rot.m33 = 1; 943 return rot; 944 } 945 946}