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