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 * Created on 08.05.2004 021 * 022 */ 023package org.biojava.nbio.structure ; 024 025import org.biojava.nbio.structure.jama.Matrix; 026import org.slf4j.Logger; 027import org.slf4j.LoggerFactory; 028 029import javax.vecmath.Matrix3d; 030import javax.vecmath.Matrix4d; 031import javax.vecmath.Point3d; 032import javax.vecmath.Vector3d; 033 034/** 035 * Utility operations on Atoms, AminoAcids, etc. 036 * <p> 037 * Currently the coordinates of an Atom are stored as an array 038 * of size 3 (double[3]). It would be more powerful to use Point3D from 039 * javax.vecmath, but unfortunately this is not a part of standard 040 * java installations, since it comes with java3d. So to keep things 041 * simple at the moment biojava does not depend on java3d. 042 * 043 * @author Andreas Prlic 044 * @since 1.4 045 * @version %I% %G% 046 */ 047 048public class Calc { 049 050 private final static Logger logger = LoggerFactory.getLogger(Calc.class); 051 052 /** 053 * calculate distance between two atoms. 054 * 055 * @param a an Atom object 056 * @param b an Atom object 057 * @return a double 058 */ 059 public static final double getDistance(Atom a, Atom b) { 060 double x = a.getX() - b.getX(); 061 double y = a.getY() - b.getY(); 062 double z = a.getZ() - b.getZ(); 063 064 double s = x * x + y * y + z * z; 065 066 return Math.sqrt(s); 067 } 068 069 070 /** 071 * Will calculate the square of distances between two atoms. This will be 072 * faster as it will not perform the final square root to get the actual 073 * distance. Use this if doing large numbers of distance comparisons - it is 074 * marginally faster than getDistance(). 075 * 076 * @param a an Atom object 077 * @param b an Atom object 078 * @return a double 079 */ 080 public static double getDistanceFast(Atom a, Atom b) { 081 double x = a.getX() - b.getX(); 082 double y = a.getY() - b.getY(); 083 double z = a.getZ() - b.getZ(); 084 085 return x * x + y * y + z * z; 086 } 087 088 public static final Atom invert(Atom a) { 089 double[] coords = new double[]{0.0,0.0,0.0} ; 090 Atom zero = new AtomImpl(); 091 zero.setCoords(coords); 092 return subtract(zero, a); 093 } 094 095 096 /** add two atoms ( a + b). 097 * 098 * @param a an Atom object 099 * @param b an Atom object 100 * @return an Atom object 101 */ 102 public static final Atom add(Atom a, Atom b){ 103 104 Atom c = new AtomImpl(); 105 c.setX( a.getX() + b.getX() ); 106 c.setY( a.getY() + b.getY() ); 107 c.setZ( a.getZ() + b.getZ() ); 108 109 return c ; 110 } 111 112 113 114 /** subtract two atoms ( a - b). 115 * 116 * @param a an Atom object 117 * @param b an Atom object 118 * @return n new Atom object representing the difference 119 120 121 */ 122 public static final Atom subtract(Atom a, Atom b) { 123 Atom c = new AtomImpl(); 124 c.setX( a.getX() - b.getX() ); 125 c.setY( a.getY() - b.getY() ); 126 c.setZ( a.getZ() - b.getZ() ); 127 128 return c ; 129 } 130 131 /** Vector product (cross product). 132 * 133 * @param a an Atom object 134 * @param b an Atom object 135 * @return an Atom object 136 */ 137 public static final Atom vectorProduct(Atom a , Atom b){ 138 139 Atom c = new AtomImpl(); 140 c.setX( a.getY() * b.getZ() - a.getZ() * b.getY() ) ; 141 c.setY( a.getZ() * b.getX() - a.getX() * b.getZ() ) ; 142 c.setZ( a.getX() * b.getY() - a.getY() * b.getX() ) ; 143 return c ; 144 145 } 146 147 /** 148 * Scalar product (dot product). 149 * 150 * @param a an Atom object 151 * @param b an Atom object 152 * @return a double 153 */ 154 public static final double scalarProduct(Atom a, Atom b) { 155 return a.getX() * b.getX() + a.getY() * b.getY() + a.getZ() * b.getZ(); 156 } 157 158 /** 159 * Gets the length of the vector (2-norm) 160 * 161 * @param a an Atom object 162 * @return Square root of the sum of the squared elements 163 */ 164 public static final double amount(Atom a){ 165 return Math.sqrt(scalarProduct(a,a)); 166 } 167 168 /** 169 * Gets the angle between two vectors 170 * 171 * @param a an Atom object 172 * @param b an Atom object 173 * @return Angle between a and b in degrees, in range [0,180]. 174 * If either vector has length 0 then angle is not defined and NaN is returned 175 */ 176 public static final double angle(Atom a, Atom b){ 177 178 179 Vector3d va = new Vector3d(a.getCoords()); 180 Vector3d vb = new Vector3d(b.getCoords()); 181 182 return Math.toDegrees(va.angle(vb)); 183 184 } 185 186 /** 187 * Returns the unit vector of vector a . 188 * 189 * @param a an Atom object 190 * @return an Atom object 191 */ 192 public static final Atom unitVector(Atom a) { 193 double amount = amount(a) ; 194 195 double[] coords = new double[3]; 196 197 coords[0] = a.getX() / amount ; 198 coords[1] = a.getY() / amount ; 199 coords[2] = a.getZ() / amount ; 200 201 a.setCoords(coords); 202 return a; 203 204 } 205 206 /** 207 * Calculate the torsion angle, i.e. the angle between the normal vectors of the 208 * two plains a-b-c and b-c-d. 209 * See http://en.wikipedia.org/wiki/Dihedral_angle 210 * @param a an Atom object 211 * @param b an Atom object 212 * @param c an Atom object 213 * @param d an Atom object 214 * @return the torsion angle in degrees, in range +-[0,180]. 215 * If either first 3 or last 3 atoms are colinear then torsion angle is not defined and NaN is returned 216 */ 217 public static final double torsionAngle(Atom a, Atom b, Atom c, Atom d) { 218 219 Atom ab = subtract(a,b); 220 Atom cb = subtract(c,b); 221 Atom bc = subtract(b,c); 222 Atom dc = subtract(d,c); 223 224 Atom abc = vectorProduct(ab,cb); 225 Atom bcd = vectorProduct(bc,dc); 226 227 double angl = angle(abc,bcd) ; 228 229 /* calc the sign: */ 230 Atom vecprod = vectorProduct(abc,bcd); 231 double val = scalarProduct(cb,vecprod); 232 if (val<0.0) angl = -angl ; 233 234 return angl; 235 } 236 237 /** 238 * Calculate the phi angle. 239 * 240 * @param a an AminoAcid object 241 * @param b an AminoAcid object 242 * @return a double 243 * @throws StructureException if aminoacids not connected or if any of the 4 needed atoms missing 244 */ 245 public static final double getPhi(AminoAcid a, AminoAcid b) throws StructureException { 246 247 if ( ! isConnected(a,b)){ 248 throw new StructureException("can not calc Phi - AminoAcids are not connected!") ; 249 } 250 251 Atom a_C = a.getC(); 252 Atom b_N = b.getN(); 253 Atom b_CA = b.getCA(); 254 Atom b_C = b.getC(); 255 256 // C and N were checked in isConnected already 257 if (b_CA==null) throw new StructureException("Can not calculate Phi, CA atom is missing"); 258 259 return torsionAngle(a_C,b_N,b_CA,b_C); 260 } 261 262 /** 263 * Calculate the psi angle. 264 * 265 * @param a an AminoAcid object 266 * @param b an AminoAcid object 267 * @return a double 268 * @throws StructureException if aminoacids not connected or if any of the 4 needed atoms missing 269 */ 270 public static final double getPsi(AminoAcid a, AminoAcid b) throws StructureException { 271 if ( ! isConnected(a,b)) { 272 throw new StructureException("can not calc Psi - AminoAcids are not connected!") ; 273 } 274 275 Atom a_N = a.getN(); 276 Atom a_CA = a.getCA(); 277 Atom a_C = a.getC(); 278 Atom b_N = b.getN(); 279 280 // C and N were checked in isConnected already 281 if (a_CA==null) throw new StructureException("Can not calculate Psi, CA atom is missing"); 282 283 return torsionAngle(a_N,a_CA,a_C,b_N); 284 285 } 286 287 /** 288 * Test if two amino acids are connected, i.e. 289 * if the distance from C to N < 2.5 Angstrom. 290 * 291 * If one of the AminoAcids has an atom missing, returns false. 292 * 293 * @param a an AminoAcid object 294 * @param b an AminoAcid object 295 * @return true if ... 296 */ 297 public static final boolean isConnected(AminoAcid a, AminoAcid b) { 298 Atom C = null ; 299 Atom N = null; 300 301 C = a.getC(); 302 N = b.getN(); 303 304 if ( C == null || N == null) 305 return false; 306 307 // one could also check if the CA atoms are < 4 A... 308 double distance = getDistance(C,N); 309 return distance < 2.5; 310 } 311 312 /** 313 * Rotate a single Atom aroud a rotation matrix. 314 * The rotation Matrix must be a pre-multiplication 3x3 Matrix. 315 * 316 * If the matrix is indexed m[row][col], then the matrix will be 317 * pre-multiplied (y=atom*M) 318 * @param atom atom to be rotated 319 * @param m a rotation matrix represented as a double[3][3] array 320 */ 321 public static final void rotate(Atom atom, double[][] m){ 322 323 double x = atom.getX(); 324 double y = atom.getY() ; 325 double z = atom.getZ(); 326 327 double nx = m[0][0] * x + m[0][1] * y + m[0][2] * z ; 328 double ny = m[1][0] * x + m[1][1] * y + m[1][2] * z ; 329 double nz = m[2][0] * x + m[2][1] * y + m[2][2] * z ; 330 331 332 atom.setX(nx); 333 atom.setY(ny); 334 atom.setZ(nz); 335 } 336 337 /** 338 * Rotate a structure. 339 * The rotation Matrix must be a pre-multiplication Matrix. 340 * 341 * @param structure a Structure object 342 * @param rotationmatrix an array (3x3) of double 343 * representing the rotation matrix. 344 * @throws StructureException ... 345 */ 346 public static final void rotate(Structure structure, 347 double[][] rotationmatrix) throws StructureException { 348 349 if ( rotationmatrix.length != 3 ) { 350 throw new StructureException ("matrix does not have size 3x3 !"); 351 } 352 AtomIterator iter = new AtomIterator(structure) ; 353 while (iter.hasNext()) { 354 Atom atom = iter.next() ; 355 Calc.rotate(atom,rotationmatrix); 356 } 357 } 358 359 /** 360 * Rotate a Group. 361 * The rotation Matrix must be a pre-multiplication Matrix. 362 * 363 * @param group a group object 364 * @param rotationmatrix an array (3x3) of double 365 * representing the rotation matrix. 366 * @throws StructureException ... 367 */ 368 public static final void rotate(Group group, 369 double[][] rotationmatrix) throws StructureException { 370 371 if ( rotationmatrix.length != 3 ) { 372 throw new StructureException ("matrix does not have size 3x3 !"); 373 } 374 AtomIterator iter = new AtomIterator(group) ; 375 while (iter.hasNext()) { 376 Atom atom = null ; 377 378 atom = iter.next() ; 379 rotate(atom,rotationmatrix); 380 381 } 382 } 383 384 /** 385 * Rotate an Atom around a Matrix object. 386 * The rotation Matrix must be a pre-multiplication Matrix. 387 * 388 * @param atom atom to be rotated 389 * @param m rotation matrix to be applied to the atom 390 */ 391 public static final void rotate(Atom atom, Matrix m){ 392 393 double x = atom.getX(); 394 double y = atom.getY(); 395 double z = atom.getZ(); 396 double[][] ad = new double[][]{{x,y,z}}; 397 398 Matrix am = new Matrix(ad); 399 Matrix na = am.times(m); 400 401 atom.setX(na.get(0,0)); 402 atom.setY(na.get(0,1)); 403 atom.setZ(na.get(0,2)); 404 405 } 406 407 /** 408 * Rotate a group object. 409 * The rotation Matrix must be a pre-multiplication Matrix. 410 * 411 * @param group a group to be rotated 412 * @param m a Matrix object representing the rotation matrix 413 */ 414 public static final void rotate(Group group, Matrix m){ 415 416 AtomIterator iter = new AtomIterator(group) ; 417 418 while (iter.hasNext()) { 419 Atom atom = iter.next() ; 420 rotate(atom,m); 421 422 } 423 424 } 425 426 /** 427 * Rotate a structure object. 428 * The rotation Matrix must be a pre-multiplication Matrix. 429 * 430 * @param structure the structure to be rotated 431 * @param m rotation matrix to be applied 432 */ 433 public static final void rotate(Structure structure, Matrix m){ 434 435 AtomIterator iter = new AtomIterator(structure) ; 436 437 while (iter.hasNext()) { 438 Atom atom = iter.next() ; 439 rotate(atom,m); 440 441 } 442 443 } 444 445 /** 446 * Transform an array of atoms at once. 447 * The transformation Matrix must be a post-multiplication Matrix. 448 * 449 * @param ca array of Atoms to shift 450 * @param t transformation Matrix4d 451 */ 452 public static void transform(Atom[] ca, Matrix4d t) { 453 for (Atom atom : ca) Calc.transform(atom, t); 454 } 455 456 /** 457 * Transforms an atom object, given a Matrix4d (i.e. the vecmath library 458 * double-precision 4x4 rotation+translation matrix). 459 * The transformation Matrix must be a post-multiplication Matrix. 460 * 461 * @param atom 462 * @param m 463 */ 464 public static final void transform (Atom atom, Matrix4d m) { 465 466 Point3d p = new Point3d(atom.getX(),atom.getY(),atom.getZ()); 467 m.transform(p); 468 469 atom.setX(p.x); 470 atom.setY(p.y); 471 atom.setZ(p.z); 472 } 473 474 /** 475 * Transforms a group object, given a Matrix4d (i.e. the vecmath library 476 * double-precision 4x4 rotation+translation matrix). 477 * The transformation Matrix must be a post-multiplication Matrix. 478 * 479 * @param group 480 * @param m 481 */ 482 public static final void transform (Group group, Matrix4d m) { 483 AtomIterator iter = new AtomIterator(group) ; 484 485 while (iter.hasNext()) { 486 Atom atom = iter.next() ; 487 transform(atom,m); 488 489 } 490 } 491 492 /** 493 * Transforms a structure object, given a Matrix4d (i.e. the vecmath library 494 * double-precision 4x4 rotation+translation matrix). 495 * The transformation Matrix must be a post-multiplication Matrix. 496 * 497 * @param structure 498 * @param m 499 */ 500 public static final void transform (Structure structure, Matrix4d m) { 501 AtomIterator iter = new AtomIterator(structure) ; 502 503 while (iter.hasNext()) { 504 Atom atom = iter.next() ; 505 transform(atom,m); 506 507 } 508 } 509 510 /** 511 * Transforms a chain object, given a Matrix4d (i.e. the vecmath library 512 * double-precision 4x4 rotation+translation matrix). 513 * The transformation Matrix must be a post-multiplication Matrix. 514 * 515 * @param chain 516 * @param m 517 */ 518 public static final void transform (Chain chain, Matrix4d m) { 519 520 for (Group g:chain.getAtomGroups()) { 521 for (Atom atom: g.getAtoms()) { 522 transform(atom,m); 523 } 524 } 525 } 526 527 /** 528 * Translates an atom object, given a Vector3d (i.e. the vecmath library 529 * double-precision 3-d vector) 530 * @param atom 531 * @param v 532 */ 533 public static final void translate (Atom atom, Vector3d v) { 534 535 atom.setX(atom.getX()+v.x); 536 atom.setY(atom.getY()+v.y); 537 atom.setZ(atom.getZ()+v.z); 538 } 539 540 /** 541 * Translates a group object, given a Vector3d (i.e. the vecmath library 542 * double-precision 3-d vector) 543 * @param group 544 * @param v 545 */ 546 public static final void translate (Group group, Vector3d v) { 547 AtomIterator iter = new AtomIterator(group) ; 548 549 while (iter.hasNext()) { 550 Atom atom = iter.next() ; 551 translate(atom,v); 552 553 } 554 } 555 556 /** 557 * Translates a chain object, given a Vector3d (i.e. the vecmath library 558 * double-precision 3-d vector) 559 * @param chain 560 * @param v 561 */ 562 public static final void translate (Chain chain, Vector3d v) { 563 564 for (Group g:chain.getAtomGroups()) { 565 for (Atom atom: g.getAtoms()) { 566 translate(atom,v); 567 } 568 } 569 } 570 571 /** 572 * Translates a Structure object, given a Vector3d (i.e. the vecmath library 573 * double-precision 3-d vector) 574 * @param structure 575 * @param v 576 */ 577 public static final void translate (Structure structure, Vector3d v) { 578 AtomIterator iter = new AtomIterator(structure) ; 579 580 while (iter.hasNext()) { 581 Atom atom = iter.next() ; 582 translate(atom,v); 583 584 } 585 } 586 587 /** calculate structure + Matrix coodinates ... 588 * 589 * @param s the structure to operate on 590 * @param matrix a Matrix object 591 */ 592 public static final void plus(Structure s, Matrix matrix){ 593 AtomIterator iter = new AtomIterator(s) ; 594 Atom oldAtom = null; 595 Atom rotOldAtom = null; 596 while (iter.hasNext()) { 597 Atom atom = null ; 598 599 atom = iter.next() ; 600 try { 601 if ( oldAtom != null){ 602 logger.debug("before {}", getDistance(oldAtom,atom)); 603 } 604 } catch (Exception e){ 605 logger.error("Exception: ", e); 606 } 607 oldAtom = (Atom)atom.clone(); 608 609 double x = atom.getX(); 610 double y = atom.getY() ; 611 double z = atom.getZ(); 612 double[][] ad = new double[][]{{x,y,z}}; 613 614 Matrix am = new Matrix(ad); 615 Matrix na = am.plus(matrix); 616 617 double[] coords = new double[3] ; 618 coords[0] = na.get(0,0); 619 coords[1] = na.get(0,1); 620 coords[2] = na.get(0,2); 621 atom.setCoords(coords); 622 try { 623 if ( rotOldAtom != null){ 624 logger.debug("after {}", getDistance(rotOldAtom,atom)); 625 } 626 } catch (Exception e){ 627 logger.error("Exception: ", e); 628 } 629 rotOldAtom = (Atom) atom.clone(); 630 } 631 632 } 633 634 635 636 /** shift a structure with a vector. 637 * 638 * @param structure a Structure object 639 * @param a an Atom object representing a shift vector 640 */ 641 public static final void shift(Structure structure, Atom a ){ 642 643 AtomIterator iter = new AtomIterator(structure) ; 644 while (iter.hasNext() ) { 645 Atom atom = null ; 646 647 atom = iter.next() ; 648 649 Atom natom = add(atom,a); 650 double x = natom.getX(); 651 double y = natom.getY() ; 652 double z = natom.getZ(); 653 atom.setX(x); 654 atom.setY(y); 655 atom.setZ(z); 656 657 } 658 } 659 660 /** Shift a vector. 661 * 662 * @param a vector a 663 * @param b vector b 664 */ 665 public static final void shift(Atom a, Atom b){ 666 667 Atom natom = add(a,b); 668 double x = natom.getX(); 669 double y = natom.getY() ; 670 double z = natom.getZ(); 671 a.setX(x); 672 a.setY(y); 673 a.setZ(z); 674 } 675 676 /** Shift a Group with a vector. 677 * 678 * @param group a group object 679 * @param a an Atom object representing a shift vector 680 */ 681 public static final void shift(Group group , Atom a ){ 682 683 AtomIterator iter = new AtomIterator(group) ; 684 while (iter.hasNext() ) { 685 Atom atom = null ; 686 687 atom = iter.next() ; 688 689 Atom natom = add(atom,a); 690 double x = natom.getX(); 691 double y = natom.getY() ; 692 double z = natom.getZ(); 693 atom.setX(x); 694 atom.setY(y); 695 atom.setZ(z); 696 697 } 698 } 699 700 701 702 /** 703 * Returns the center of mass of the set of atoms. 704 * @param atomSet a set of Atoms 705 * @return an Atom representing the Centroid of the set of atoms 706 */ 707 public static final Atom getCentroid(Atom[] atomSet){ 708 709 // if we don't catch this case, the centroid returned is (NaN,NaN,NaN), which can cause lots of problems down the line 710 if (atomSet.length==0) 711 throw new IllegalArgumentException("Atom array has length 0, can't calculate centroid!"); 712 713 double[] coords = new double[3]; 714 715 coords[0] = 0; 716 coords[1] = 0; 717 coords[2] = 0 ; 718 719 for (Atom a : atomSet) { 720 coords[0] += a.getX(); 721 coords[1] += a.getY(); 722 coords[2] += a.getZ(); 723 } 724 725 int n = atomSet.length; 726 coords[0] = coords[0] / n; 727 coords[1] = coords[1] / n; 728 coords[2] = coords[2] / n; 729 730 Atom vec = new AtomImpl(); 731 vec.setCoords(coords); 732 return vec; 733 734 } 735 736 public static Atom centerOfMass(Atom[] points) { 737 Atom center = new AtomImpl(); 738 739 float totalMass = 0.0f; 740 for (Atom a : points) { 741 float mass = a.getElement().getAtomicMass(); 742 totalMass += mass; 743 center = scaleAdd(mass, a, center); 744 } 745 746 center = scaleEquals(center, 1.0f/totalMass); 747 return center; 748 } 749 750 /** 751 * Multiply elements of a by s (in place) 752 * @param a 753 * @param s 754 * @return the modified a 755 */ 756 public static Atom scaleEquals(Atom a, double s) { 757 double x = a.getX(); 758 double y = a.getY(); 759 double z = a.getZ(); 760 761 x *= s; 762 y *= s; 763 z *= s; 764 765 //Atom b = new AtomImpl(); 766 a.setX(x); 767 a.setY(y); 768 a.setZ(z); 769 770 return a; 771 } 772 773 /** 774 * Multiply elements of a by s 775 * @param a 776 * @param s 777 * @return A new Atom with s*a 778 */ 779 public static Atom scale(Atom a, double s) { 780 double x = a.getX(); 781 double y = a.getY(); 782 double z = a.getZ(); 783 784 Atom b = new AtomImpl(); 785 b.setX(x*s); 786 b.setY(y*s); 787 b.setZ(z*s); 788 789 return b; 790 } 791 792 793 /** 794 * Perform linear transformation s*X+B, and store the result in b 795 * @param s Amount to scale x 796 * @param x Input coordinate 797 * @param b Vector to translate (will be modified) 798 * @return b, after modification 799 */ 800 public static Atom scaleAdd(double s, Atom x, Atom b){ 801 802 double xc = s*x.getX() + b.getX(); 803 double yc = s*x.getY() + b.getY(); 804 double zc = s*x.getZ() + b.getZ(); 805 806 //Atom a = new AtomImpl(); 807 b.setX(xc); 808 b.setY(yc); 809 b.setZ(zc); 810 811 return b; 812 } 813 814 /** Returns the Vector that needs to be applied to shift a set of atoms 815 * to the Centroid. 816 * @param atomSet array of Atoms 817 * @return the vector needed to shift the set of atoms to its geometric center 818 */ 819 public static final Atom getCenterVector(Atom[] atomSet){ 820 Atom centroid = getCentroid(atomSet); 821 822 return getCenterVector(atomSet,centroid); 823 824 } 825 826 /** Returns the Vector that needs to be applied to shift a set of atoms 827 * to the Centroid, if the centroid is already known 828 * @param atomSet array of Atoms 829 * @return the vector needed to shift the set of atoms to its geometric center 830 */ 831 public static final Atom getCenterVector(Atom[] atomSet, Atom centroid){ 832 833 834 double[] coords = new double[3]; 835 coords[0] = 0 - centroid.getX(); 836 coords[1] = 0 - centroid.getY(); 837 coords[2] = 0 - centroid.getZ(); 838 839 Atom shiftVec = new AtomImpl(); 840 shiftVec.setCoords(coords); 841 return shiftVec; 842 843 } 844 845 846 /** Center the atoms at the Centroid. 847 * @param atomSet a set of Atoms 848 * @return an Atom representing the Centroid of the set of atoms 849 * @throws StructureException 850 * */ 851 public static final Atom[] centerAtoms(Atom[] atomSet) throws StructureException { 852 853 Atom centroid = getCentroid(atomSet); 854 return centerAtoms(atomSet, centroid); 855 } 856 857 /** Center the atoms at the Centroid, if the centroid is already know. 858 * @param atomSet a set of Atoms 859 * @return an Atom representing the Centroid of the set of atoms 860 * @throws StructureException 861 * */ 862 public static final Atom[] centerAtoms(Atom[] atomSet, Atom centroid) throws StructureException { 863 864 Atom shiftVector = getCenterVector(atomSet, centroid); 865 866 Atom[] newAtoms = new AtomImpl[atomSet.length]; 867 868 for (int i =0 ; i < atomSet.length; i++){ 869 Atom a = atomSet[i]; 870 Atom n = add(a,shiftVector); 871 newAtoms[i] = n ; 872 } 873 return newAtoms; 874 } 875 876 877 878 879 /** creates a virtual C-beta atom. this might be needed when working with GLY 880 * 881 * thanks to Peter Lackner for a python template of this method. 882 * @param amino the amino acid for which a "virtual" CB atom should be calculated 883 * @return a "virtual" CB atom 884 * @throws StructureException 885 */ 886 public static final Atom createVirtualCBAtom(AminoAcid amino) 887 throws StructureException{ 888 889 AminoAcid ala = StandardAminoAcid.getAminoAcid("ALA"); 890 Atom aN = ala.getN(); 891 Atom aCA = ala.getCA(); 892 Atom aC = ala.getC(); 893 Atom aCB = ala.getCB(); 894 895 896 Atom[] arr1 = new Atom[3]; 897 arr1[0] = aN; 898 arr1[1] = aCA; 899 arr1[2] = aC; 900 901 Atom[] arr2 = new Atom[3]; 902 arr2[0] = amino.getN(); 903 arr2[1] = amino.getCA(); 904 arr2[2] = amino.getC(); 905 906 // ok now we got the two arrays, do a SVD: 907 908 SVDSuperimposer svd = new SVDSuperimposer(arr2,arr1); 909 910 Matrix rotMatrix = svd.getRotation(); 911 Atom tranMatrix = svd.getTranslation(); 912 913 Calc.rotate(aCB,rotMatrix); 914 915 Atom virtualCB = Calc.add(aCB,tranMatrix); 916 virtualCB.setName("CB"); 917 918 return virtualCB; 919 } 920 921 922 /** 923 * Gets euler angles for a matrix given in ZYZ convention. 924 * (as e.g. used by Jmol) 925 * 926 * @param m the rotation matrix 927 * @return the euler values for a rotation around Z, Y, Z in degrees... 928 */ 929 public static final double[] getZYZEuler(Matrix m) { 930 double m22 = m.get(2,2); 931 double rY = Math.toDegrees(Math.acos(m22)); 932 double rZ1, rZ2; 933 if (m22 > .999d || m22 < -.999d) { 934 rZ1 = Math.toDegrees(Math.atan2(m.get(1,0), m.get(1,1))); 935 rZ2 = 0; 936 } else { 937 rZ1 = Math.toDegrees(Math.atan2(m.get(2,1), -m.get(2,0))); 938 rZ2 = Math.toDegrees(Math.atan2(m.get(1,2), m.get(0,2))); 939 } 940 return new double[] {rZ1,rY,rZ2}; 941 } 942 943 944 /** Convert a rotation Matrix to Euler angles. 945 * This conversion uses conventions as described on page: 946 * http://www.euclideanspace.com/maths/geometry/rotations/euler/index.htm 947 * Coordinate System: right hand 948 * Positive angle: right hand 949 * Order of euler angles: heading first, then attitude, then bank 950 * 951 * @param m the rotation matrix 952 * @return a array of three doubles containing the three euler angles in radians 953 */ 954 public static final double[] getXYZEuler(Matrix m){ 955 double heading, attitude, bank; 956 957 // Assuming the angles are in radians. 958 if (m.get(1,0) > 0.998) { // singularity at north pole 959 heading = Math.atan2(m.get(0,2),m.get(2,2)); 960 attitude = Math.PI/2; 961 bank = 0; 962 963 } else if (m.get(1,0) < -0.998) { // singularity at south pole 964 heading = Math.atan2(m.get(0,2),m.get(2,2)); 965 attitude = -Math.PI/2; 966 bank = 0; 967 968 } else { 969 heading = Math.atan2(-m.get(2,0),m.get(0,0)); 970 bank = Math.atan2(-m.get(1,2),m.get(1,1)); 971 attitude = Math.asin(m.get(1,0)); 972 } 973 return new double[] { heading, attitude, bank }; 974 } 975 976 977 978 /** This conversion uses NASA standard aeroplane conventions as described on page: 979 * http://www.euclideanspace.com/maths/geometry/rotations/euler/index.htm 980 * Coordinate System: right hand 981 * Positive angle: right hand 982 * Order of euler angles: heading first, then attitude, then bank. 983 * matrix row column ordering: 984 * [m00 m01 m02] 985 * [m10 m11 m12] 986 * [m20 m21 m22] 987 * @param heading in radians 988 * @param attitude in radians 989 * @param bank in radians 990 * @return the rotation matrix */ 991 public static final Matrix matrixFromEuler(double heading, double attitude, double bank) { 992 // Assuming the angles are in radians. 993 double ch = Math.cos(heading); 994 double sh = Math.sin(heading); 995 double ca = Math.cos(attitude); 996 double sa = Math.sin(attitude); 997 double cb = Math.cos(bank); 998 double sb = Math.sin(bank); 999 1000 Matrix m = new Matrix(3,3); 1001 m.set(0,0, ch * ca); 1002 m.set(0,1, sh*sb - ch*sa*cb); 1003 m.set(0,2, ch*sa*sb + sh*cb); 1004 m.set(1,0, sa); 1005 m.set(1,1, ca*cb); 1006 m.set(1,2, -ca*sb); 1007 m.set(2,0, -sh*ca); 1008 m.set(2,1, sh*sa*cb + ch*sb); 1009 m.set(2,2, -sh*sa*sb + ch*cb); 1010 1011 return m; 1012 } 1013 1014 1015 /** 1016 * Calculates the angle from centerPt to targetPt in degrees. 1017 * The return should range from [0,360), rotating CLOCKWISE, 1018 * 0 and 360 degrees represents NORTH, 1019 * 90 degrees represents EAST, etc... 1020 * 1021 * Assumes all points are in the same coordinate space. If they are not, 1022 * you will need to call SwingUtilities.convertPointToScreen or equivalent 1023 * on all arguments before passing them to this function. 1024 * 1025 * @param centerPt Point we are rotating around. 1026 * @param targetPt Point we want to calculate the angle to. 1027 * @return angle in degrees. This is the angle from centerPt to targetPt. 1028 */ 1029 public static double calcRotationAngleInDegrees(Atom centerPt, Atom targetPt) 1030 { 1031 // calculate the angle theta from the deltaY and deltaX values 1032 // (atan2 returns radians values from [-PI,PI]) 1033 // 0 currently points EAST. 1034 // NOTE: By preserving Y and X param order to atan2, we are expecting 1035 // a CLOCKWISE angle direction. 1036 double theta = Math.atan2(targetPt.getY() - centerPt.getY(), targetPt.getX() - centerPt.getX()); 1037 1038 // rotate the theta angle clockwise by 90 degrees 1039 // (this makes 0 point NORTH) 1040 // NOTE: adding to an angle rotates it clockwise. 1041 // subtracting would rotate it counter-clockwise 1042 theta += Math.PI/2.0; 1043 1044 // convert from radians to degrees 1045 // this will give you an angle from [0->270],[-180,0] 1046 double angle = Math.toDegrees(theta); 1047 1048 // convert to positive range [0-360) 1049 // since we want to prevent negative angles, adjust them now. 1050 // we can assume that atan2 will not return a negative value 1051 // greater than one partial rotation 1052 if (angle < 0) { 1053 angle += 360; 1054 } 1055 1056 return angle; 1057 } 1058 1059 1060 public static void main(String[] args){ 1061 Atom a =new AtomImpl(); 1062 a.setX(0); 1063 a.setY(0); 1064 a.setZ(0); 1065 1066 Atom b = new AtomImpl(); 1067 b.setX(1); 1068 b.setY(1); 1069 b.setZ(0); 1070 1071 logger.info("Angle between atoms: ", calcRotationAngleInDegrees(a, b)); 1072 } 1073 1074 1075 public static void rotate(Atom[] ca, Matrix matrix) { 1076 for (Atom atom : ca) Calc.rotate(atom, matrix); 1077 } 1078 1079 /** 1080 * Shift an array of atoms at once. 1081 * @param ca array of Atoms to shift 1082 * @param b reference Atom vector 1083 */ 1084 public static void shift(Atom[] ca, Atom b) { 1085 for (Atom atom : ca) Calc.shift(atom, b); 1086 } 1087 1088 /** 1089 * Convert JAMA rotation and translation to a Vecmath transformation matrix. 1090 * Because the JAMA matrix is a pre-multiplication matrix and the Vecmath 1091 * matrix is a post-multiplication one, the rotation matrix is transposed to 1092 * ensure that the transformation they produce is the same. 1093 * 1094 * @param rot 3x3 Rotation matrix 1095 * @param trans 3x1 Translation matrix 1096 * @return 4x4 transformation matrix 1097 */ 1098 public static Matrix4d getTransformation(Matrix rot, Matrix trans) { 1099 return new Matrix4d( new Matrix3d(rot.getColumnPackedCopy()), 1100 new Vector3d(trans.getColumnPackedCopy()), 1101 1.0); 1102 } 1103 1104 /** 1105 * Convert JAMA rotation and translation to a Vecmath transformation matrix. 1106 * Because the JAMA matrix is a pre-multiplication matrix and the Vecmath 1107 * matrix is a post-multiplication one, the rotation matrix is transposed to 1108 * ensure that the transformation they produce is the same. 1109 * 1110 * @param rot 3x3 Rotation matrix 1111 * @param trans 3x1 translation vector in Atom coordinates 1112 * @return 4x4 transformation matrix 1113 */ 1114 public static Matrix4d getTransformation(Matrix rot, Atom trans) { 1115 return new Matrix4d( new Matrix3d(rot.getColumnPackedCopy()), 1116 new Vector3d(trans.getCoords()), 1117 1.0); 1118 } 1119 1120 /** 1121 * Convert Vecmath transformation into a JAMA rotation matrix. 1122 * Because the JAMA matrix is a pre-multiplication matrix and the Vecmath 1123 * matrix is a post-multiplication one, the rotation matrix is transposed to 1124 * ensure that the transformation they produce is the same. 1125 * 1126 * @param transform Matrix4d with transposed rotation matrix 1127 * @return 1128 */ 1129 public static Matrix getRotationMatrix(Matrix4d transform){ 1130 1131 Matrix rot = new Matrix(3,3); 1132 for (int i=0;i<3;i++) { 1133 for (int j=0;j<3;j++) { 1134 rot.set(j, i, transform.getElement(i, j)); //transposed 1135 } 1136 } 1137 return rot; 1138 } 1139 1140 /** 1141 * Extract the translational vector of a Vecmath transformation. 1142 * 1143 * @param transform Matrix4d 1144 * @return Atom shift vector 1145 */ 1146 public static Atom getTranslationVector(Matrix4d transform){ 1147 1148 Atom transl = new AtomImpl(); 1149 double[] coords = {transform.m03, transform.m13, transform.m23}; 1150 transl.setCoords(coords); 1151 return transl; 1152 } 1153 1154 /** 1155 * Convert an array of atoms into an array of vecmath points 1156 * @param atoms list of atoms 1157 * @return list of Point3ds storing the x,y,z coordinates of each atom 1158 */ 1159 public static Point3d[] atomsToPoints(Atom[] atoms) { 1160 Point3d[] points = new Point3d[atoms.length]; 1161 for(int i = 0; i< atoms.length;i++) { 1162 points[i] = new Point3d(atoms[i].getCoords()); 1163 } 1164 return points; 1165 } 1166}