/Bio/PDB/Superimposer.py

http://github.com/biopython/biopython · Python · 57 lines · 29 code · 9 blank · 19 comment · 4 complexity · e4e7af7f6f6677cf232d07e8b0ca24c8 MD5 · raw file

  1. # Copyright (C) 2002, Thomas Hamelryck (thamelry@binf.ku.dk)
  2. #
  3. # This file is part of the Biopython distribution and governed by your
  4. # choice of the "Biopython License Agreement" or the "BSD 3-Clause License".
  5. # Please see the LICENSE file that should have been included as part of this
  6. # package.
  7. """Superimpose two structures."""
  8. import numpy
  9. from Bio.SVDSuperimposer import SVDSuperimposer
  10. from Bio.PDB.PDBExceptions import PDBException
  11. class Superimposer:
  12. """Rotate/translate one set of atoms on top of another to minimize RMSD."""
  13. def __init__(self):
  14. """Initialize the class."""
  15. self.rotran = None
  16. self.rms = None
  17. def set_atoms(self, fixed, moving):
  18. """Prepare translation/rotation to minimize RMSD between atoms.
  19. Put (translate/rotate) the atoms in fixed on the atoms in
  20. moving, in such a way that the RMSD is minimized.
  21. :param fixed: list of (fixed) atoms
  22. :param moving: list of (moving) atoms
  23. :type fixed,moving: [L{Atom}, L{Atom},...]
  24. """
  25. if not len(fixed) == len(moving):
  26. raise PDBException("Fixed and moving atom lists differ in size")
  27. length = len(fixed)
  28. fixed_coord = numpy.zeros((length, 3))
  29. moving_coord = numpy.zeros((length, 3))
  30. for i in range(0, length):
  31. fixed_coord[i] = fixed[i].get_coord()
  32. moving_coord[i] = moving[i].get_coord()
  33. sup = SVDSuperimposer()
  34. sup.set(fixed_coord, moving_coord)
  35. sup.run()
  36. self.rms = sup.get_rms()
  37. self.rotran = sup.get_rotran()
  38. def apply(self, atom_list):
  39. """Rotate/translate a list of atoms."""
  40. if self.rotran is None:
  41. raise PDBException("No transformation has been calculated yet")
  42. rot, tran = self.rotran
  43. rot = rot.astype("f")
  44. tran = tran.astype("f")
  45. for atom in atom_list:
  46. atom.transform(rot, tran)