IMP logo
IMP Reference Guide  2.20.2
The Integrative Modeling Platform
pmi1/representation.py
1 #!/usr/bin/env python
2 
3 """@namespace IMP.pmi1.representation
4  Representation of the system.
5 """
6 
7 from __future__ import print_function
8 import IMP
9 import IMP.core
10 import IMP.algebra
11 import IMP.atom
12 import IMP.display
13 import IMP.isd
14 import IMP.pmi1
15 import IMP.pmi1.tools
16 import IMP.pmi1.output
17 import IMP.rmf
18 import IMP.pmi1.topology
19 import RMF
20 from math import pi, sqrt
21 from operator import itemgetter
22 import os
23 import weakref
24 
25 class _Repo(object):
26  def __init__(self, doi, root):
27  self.doi = doi
28  self._root = root
29 
30  def get_fname(self, fname):
31  """Return a path relative to the top of the repository"""
32  return os.path.relpath(fname, self._root)
33 
34 class _StateInfo(object):
35  """Score state-specific information about this representation."""
36  short_name = None
37  long_name = None
38 
39 class Representation(object):
40  # Authors: Peter Cimermancic, Riccardo Pellarin, Charles Greenberg
41 
42  '''
43  Set up the representation of all proteins and nucleic acid macromolecules.
44 
45  Create the molecular hierarchies, representation,
46  sequence connectivity for the various involved proteins and
47  nucleic acid macromolecules:
48 
49  Create a protein, DNA or RNA, represent it as a set of connected balls of appropriate
50  radii and number of residues, PDB at given resolution(s), or ideal helices.
51 
52  How to use the SimplifiedModel class (typical use):
53 
54  see test/test_hierarchy_contruction.py
55 
56  examples:
57 
58  1) Create a chain of helices and flexible parts
59 
60  c_1_119 =self.add_component_necklace("prot1",1,119,20)
61  c_120_131 =self.add_component_ideal_helix("prot1",resolutions=[1,10],resrange=(120,131))
62  c_132_138 =self.add_component_beads("prot1",[(132,138)])
63  c_139_156 =self.add_component_ideal_helix("prot1",resolutions=[1,10],resrange=(139,156))
64  c_157_174 =self.add_component_beads("prot1",[(157,174)])
65  c_175_182 =self.add_component_ideal_helix("prot1",resolutions=[1,10],resrange=(175,182))
66  c_183_194 =self.add_component_beads("prot1",[(183,194)])
67  c_195_216 =self.add_component_ideal_helix("prot1",resolutions=[1,10],resrange=(195,216))
68  c_217_250 =self.add_component_beads("prot1",[(217,250)])
69 
70 
71  self.set_rigid_body_from_hierarchies(c_120_131)
72  self.set_rigid_body_from_hierarchies(c_139_156)
73  self.set_rigid_body_from_hierarchies(c_175_182)
74  self.set_rigid_body_from_hierarchies(c_195_216)
75 
76  clist=[c_1_119,c_120_131,c_132_138,c_139_156,c_157_174,c_175_182,c_183_194,c_195_216,
77  c_217_250]
78 
79  self.set_chain_of_super_rigid_bodies(clist,2,3)
80 
81  self.set_super_rigid_bodies(["prot1"])
82 
83  '''
84 
85  def __init__(self, model, upperharmonic=True, disorderedlength=True):
86  """Constructor.
87  @param model the model
88  @param upperharmonic This flag uses either harmonic (False)
89  or upperharmonic (True) in the intra-pair
90  connectivity restraint.
91  @param disorderedlength This flag uses either disordered length
92  calculated for random coil peptides (True) or zero
93  surface-to-surface distance between beads (False)
94  as optimal distance for the sequence connectivity
95  restraint.
96  """
97 
98  self.state = _StateInfo()
99  self._metadata = []
100  self._file_dataset = {}
101  self._protocol_output = []
102  self._non_modeled_components = {}
103 
104  # this flag uses either harmonic (False) or upperharmonic (True)
105  # in the intra-pair connectivity restraint. Harmonic is used when you want to
106  # remove the intra-ev term from energy calculations, e.g.:
107  # upperharmonic=False
108  # ip=self.get_connected_intra_pairs()
109  # ev.add_excluded_particle_pairs(ip)
110 
111  self.upperharmonic = upperharmonic
112  self.disorderedlength = disorderedlength
113  self.rigid_bodies = []
114  self.fixed_rigid_bodies = []
115  self.fixed_floppy_bodies = []
116  self.floppy_bodies = []
117  # self.super_rigid_bodies is a list of tuples.
118  # each tuple, corresponding to a particular super rigid body
119  # the tuple is (super_rigid_xyzs,super_rigid_rbs)
120  # where super_rigid_xyzs are the flexible xyz particles
121  # and super_rigid_rbs is the list of rigid bodies.
122  self.super_rigid_bodies = []
123  self.rigid_body_symmetries = []
124  self.output_level = "low"
125  self.label = "None"
126 
127  self.maxtrans_rb = 2.0
128  self.maxrot_rb = 0.04
129  self.maxtrans_srb = 2.0
130  self.maxrot_srb = 0.2
131  self.rigidbodiesarefixed = False
132  self.floppybodiesarefixed = False
133  self.maxtrans_fb = 3.0
134  self.resolution = 10.0
135  self.bblenght = 100.0
136  self.kappa = 100.0
137  self.model = model
138 
139  self.representation_is_modified = False
140  self.unmodeledregions_cr_dict = {}
141  self.sortedsegments_cr_dict = {}
142  self.prot = IMP.atom.Hierarchy.setup_particle(IMP.Particle(self.model))
143  self.connected_intra_pairs = []
144  self.hier_dict = {}
145  self.color_dict = {}
146  self.sequence_dict = {}
147  self.hier_geometry_pairs = {}
148  self.hier_db = IMP.pmi1.tools.HierarchyDatabase()
149  # this dictionary stores the hierarchies by component name and representation type
150  # self.hier_representation[name][representation_type]
151  # where representation type is Res:X, Beads, Densities, Representation,
152  # etc...
153  self.hier_representation = {}
154  self.hier_resolution = {}
155  # reference structures is a dictionary that contains the coordinates of
156  # structures that are used to calculate the rmsd
157  self.reference_structures = {}
158  self.elements = {}
159  self.linker_restraints = IMP.RestraintSet(self.model, "linker_restraints")
160  self.linker_restraints.set_was_used(True)
161  self.linker_restraints_dict = {}
162  self.threetoone = {'ALA': 'A', 'ARG': 'R', 'ASN': 'N', 'ASP': 'D',
163  'CYS': 'C', 'GLU': 'E', 'GLN': 'Q', 'GLY': 'G',
164  'HIS': 'H', 'ILE': 'I', 'LEU': 'L', 'LYS': 'K',
165  'MET': 'M', 'PHE': 'F', 'PRO': 'P', 'SER': 'S',
166  'THR': 'T', 'TRP': 'W', 'TYR': 'Y', 'VAL': 'V', 'UNK': 'X'}
167 
168  self.onetothree = dict((v, k) for k, v in self.threetoone.items())
169 
170  self.residuenamekey = IMP.StringKey("ResidueName")
171 
172  @property
173  @IMP.deprecated_method("3.0", "Model should be accessed with `.model`.")
174  def m(self):
175  return self.model
176 
177  def add_metadata(self, m):
178  """Associate some metadata with this modeling.
179  @param m an instance of IMP.pmi1.metadata.Metadata or a subclass.
180  """
181  self._metadata.append(m)
182 
183  def set_file_dataset(self, fname, dataset):
184  """Associate a dataset with a filename.
185  This can be used to identify how the file was produced (in many
186  cases IMP can determine this automatically from a file header or
187  other metadata, but not always). For example, a manually-produced
188  PDB file (not from the PDB database or Modeller) can be
189  identified this way.
190  @param fname filename
191  @dataset the ihm.dataset.Dataset object to associate.
192  """
193  self._file_dataset[os.path.abspath(fname)] = dataset
194 
195  def get_file_dataset(self, fname):
196  """Get the dataset associated with a filename, or None.
197  @param fname filename
198  @return an ihm.dataset.Dataset, or None.
199  """
200  return self._file_dataset.get(os.path.abspath(fname), None)
201 
202  def add_protocol_output(self, p):
203  """Capture details of the modeling protocol.
204  @param p an instance of IMP.pmi1.output.ProtocolOutput or a subclass.
205  """
206  state = p._add_state(self)
207  self._protocol_output.append((p, state))
208  p._each_metadata.append(self._metadata)
209  p._file_datasets.append(self._file_dataset)
210  state.model = self.model
211  state.prot = self.prot
212  protocol_output = property(lambda self:
213  [x[0] for x in self._protocol_output])
214 
215  def set_label(self, label):
216  self.label = label
217 
218  def create_component(self, name, color=0.0):
219  # Note that we set up the component as *both* a Chain and a Molecule.
220  # This is because old PMI1 code expects the top-level particle to be
221  # a Molecule, but we also need it to be a Chain to set the sequence.
222  # This looks a little odd but is a valid IMP atom::Hierarchy.
223  protein_h = IMP.atom.Chain.setup_particle(IMP.Particle(self.model), 'X')
225  protein_h.set_name(name)
226  self.hier_dict[name] = protein_h
227  self.hier_representation[name] = {}
228  self.hier_db.add_name(name)
229  self.prot.add_child(protein_h)
230  self.color_dict[name] = color
231  self.elements[name] = []
232  for p, state in self._protocol_output:
233  p.create_component(state, name, True)
234 
235  def create_transformed_component(self, name, original, transform):
236  """Create a transformed view of a component.
237  This does not copy the coordinates or representation of the
238  component, and the transformed component is not visible to IMP,
239  but when the model is written to a ProtocolOutput, the transformed
240  component is added. This can be used to easily add symmetry-related
241  'copies' of a component without copying the underlying
242  IMP objects."""
243  for p, state in self._protocol_output:
244  p.create_transformed_component(state, name, original, transform)
245 
247  """Create a component that isn't used in the modeling.
248  No coordinates or other structural data for this component will
249  be read or written, but a primary sequence can be assigned. This
250  is useful if the input experimental data is of a system larger
251  than that modeled. Any references to these non-modeled components
252  can then be correctly resolved later."""
253  self._non_modeled_components[name] = None
254  self.elements[name] = []
255  for p, state in self._protocol_output:
256  p.create_component(state, name, False)
257 
258  # Deprecation warning
259 
260  @IMP.deprecated_method("2.5", "Use create_component() instead.")
261  def add_component_name(self, *args, **kwargs):
262  self.create_component(*args, **kwargs)
263 
264  def get_component_names(self):
265  return list(self.hier_dict.keys())
266 
267  def add_component_sequence(self, name, filename, id=None, offs=None,
268  format="FASTA"):
269  '''
270  Add the primary sequence for a single component.
271 
272  @param name Human-readable name of the component
273  @param filename Name of the FASTA file
274  @param id Identifier of the sequence in the FASTA file header
275  (if not provided, use `name` instead)
276  '''
277  record_dict = IMP.pmi1.topology.Sequences(filename)
278  if id is None:
279  id = name
280  if id not in record_dict:
281  raise KeyError("id %s not found in fasta file" % id)
282  length = len(record_dict[id])
283  self.sequence_dict[name] = str(record_dict[id])
284  # No Hierarchy for this component if it is non-modeled
285  if name not in self._non_modeled_components:
286  protein_h = self.hier_dict[name]
287  protein_h.set_sequence(self.sequence_dict[name])
288  if offs is not None:
289  offs_str="-"*offs
290  self.sequence_dict[name]=offs_str+self.sequence_dict[name]
291 
292  self.elements[name].append((length, length, " ", "end"))
293  for p, state in self._protocol_output:
294  p.add_component_sequence(state, name, self.sequence_dict[name])
295 
296  def autobuild_model(self, name, pdbname, chain,
297  resolutions=None, resrange=None,
298  missingbeadsize=20,
299  color=None, pdbresrange=None, offset=0,
300  show=False, isnucleicacid=False,
301  attachbeads=False):
302 
303  self.representation_is_modified = True
304  outhiers = []
305 
306  if color is None:
307  color = self.color_dict[name]
308  else:
309  self.color_dict[name] = color
310 
311  if resolutions is None:
312  resolutions = [1]
313  print("autobuild_model: constructing %s from pdb %s and chain %s" % (name, pdbname, str(chain)))
314 
315  # get the initial and end residues of the pdb
316  t = IMP.atom.read_pdb(pdbname, self.model,
319 
320  # find start and end indexes
321 
322  start = IMP.atom.Residue(
323  t.get_children()[0].get_children()[0]).get_index()
324  end = IMP.atom.Residue(
325  t.get_children()[0].get_children()[-1]).get_index()
326 
327  # check if resrange was defined, otherwise
328  # use the sequence, or the pdb resrange
329 
330  if resrange is None:
331  if name in self.sequence_dict:
332  resrange = (1, len(self.sequence_dict[name]))
333  else:
334  resrange = (start + offset, end + offset)
335  else:
336  if resrange[1] in (-1, 'END'):
337  resrange = (resrange[0],end)
338  start = resrange[0] - offset
339  end = resrange[1] - offset
340 
342  t,
343  resrange[0],
344  resrange[1])
345 
347  t,
348  resrange[0],
349  terminus="N")
351  t,
352  resrange[1],
353  terminus="C")
354  # Done with the PDB
355  IMP.atom.destroy(t)
356 
357  # construct pdb fragments and intervening beads
358  for n, g in enumerate(gaps):
359  first = g[0]
360  last = g[1]
361  if g[2] == "cont":
362  print("autobuild_model: constructing fragment %s from pdb" % (str((first, last))))
363  outhiers += self.add_component_pdb(name, pdbname,
364  chain, resolutions=resolutions,
365  color=color, cacenters=True,
366  resrange=(first, last),
367  offset=offset, isnucleicacid=isnucleicacid)
368  elif g[2] == "gap" and n > 0:
369  print("autobuild_model: constructing fragment %s as a bead" % (str((first, last))))
370  parts = self.hier_db.get_particles_at_closest_resolution(name,
371  first + offset - 1,
372  1)
373  xyz = IMP.core.XYZ(parts[0]).get_coordinates()
374  outhiers += self.add_component_necklace(name,
375  first+offset, last+offset, missingbeadsize, incoord=xyz)
376 
377  elif g[2] == "gap" and n == 0:
378  # add pre-beads
379  print("autobuild_model: constructing fragment %s as a bead" % (str((first, last))))
380  outhiers += self.add_component_necklace(name,
381  first+offset, last+offset, missingbeadsize, incoord=xyznter)
382 
383  return outhiers
384 
385  # Deprecation warning
386 
387  @IMP.deprecated_method("2.5", "Use autobuild_model() instead.")
388  def autobuild_pdb_and_intervening_beads(self, *args, **kwargs):
389  r = self.autobuild_model(*args, **kwargs)
390  return r
391 
392  def add_component_pdb(self, name, pdbname, chain, resolutions, color=None,
393  resrange=None, offset=0, cacenters=True, show=False,
394  isnucleicacid=False, readnonwateratoms=False,
395  read_ca_cb_only=False):
396  '''
397  Add a component that has an associated 3D structure in a PDB file.
398 
399  Reads the PDB, and constructs the fragments corresponding to contiguous
400  sequence stretches.
401 
402  @return a list of hierarchies.
403 
404  @param name (string) the name of the component
405  @param pdbname (string) the name of the PDB file
406  @param chain (string or integer) can be either a string (eg, "A")
407  or an integer (eg, 0 or 1) in case you want
408  to get the corresponding chain number in the PDB.
409  @param resolutions (integers) a list of integers that corresponds
410  to the resolutions that have to be generated
411  @param color (float from 0 to 1) the color applied to the
412  hierarchies generated
413  @param resrange (tuple of integers): the residue range to extract
414  from the PDB. It is a tuple (beg,end). If not specified,
415  all residues belonging to the specified chain are read.
416  @param offset (integer) specifies the residue index offset to be
417  applied when reading the PDB (the FASTA sequence is
418  assumed to start from residue 1, so use this parameter
419  if the PDB file does not start at residue 1)
420  @param cacenters (boolean) if True generates resolution=1 beads
421  centered on C-alpha atoms.
422  @param show (boolean) print out the molecular hierarchy at the end.
423  @param isnucleicacid (boolean) use True if you're reading a PDB
424  with nucleic acids.
425  @param readnonwateratoms (boolean) if True fixes some pathological PDB.
426  @param read_ca_cb_only (boolean) if True, only reads CA/CB
427  '''
428 
429  self.representation_is_modified = True
430  if color is None:
431  # if the color is not passed, then get the stored color
432  color = self.color_dict[name]
433  protein_h = self.hier_dict[name]
434  outhiers = []
435 
436  # determine selector
438  if read_ca_cb_only:
439  cacbsel = IMP.atom.OrPDBSelector(
442  sel = IMP.atom.AndPDBSelector(cacbsel, sel)
443  if type(chain) == str:
444  sel = IMP.atom.ChainPDBSelector([chain]) & sel
445  t = IMP.atom.read_pdb(pdbname, self.model, sel)
446 
447  # get the first and last residue
448  start = IMP.atom.Residue(
449  t.get_children()[0].get_children()[0]).get_index()
450  end = IMP.atom.Residue(
451  t.get_children()[0].get_children()[-1]).get_index()
452  c = IMP.atom.Chain(IMP.atom.get_by_type(t, IMP.atom.CHAIN_TYPE)[0])
453  else:
454  t = IMP.atom.read_pdb(pdbname, self.model, sel)
455  c = IMP.atom.Chain(
456  IMP.atom.get_by_type(t, IMP.atom.CHAIN_TYPE)[chain])
457 
458  # get the first and last residue
459  start = IMP.atom.Residue(c.get_children()[0]).get_index()
460  end = IMP.atom.Residue(c.get_children()[-1]).get_index()
461  chain = c.get_id()
462 
463  if not resrange is None:
464  if resrange[0] > start:
465  start = resrange[0]
466  if resrange[1] < end:
467  end = resrange[1]
468 
469  if not isnucleicacid:
470  # do what you have to do for proteins
471  sel = IMP.atom.Selection(
472  c,
473  residue_indexes=list(range(
474  start,
475  end + 1)),
476  atom_type=IMP.atom.AT_CA)
477 
478  else:
479  # do what you have to do for nucleic-acids
480  # to work, nucleic acids should not be indicated as HETATM in the pdb
481  sel = IMP.atom.Selection(
482  c,
483  residue_indexes=list(range(
484  start,
485  end + 1)),
486  atom_type=IMP.atom.AT_P)
487 
488  ps = sel.get_selected_particles()
489  if len(ps) == 0:
490  raise ValueError("%s no residue found in pdb %s chain %s that overlaps with the queried stretch %s-%s" \
491  % (name, pdbname, str(chain), str(resrange[0]), str(resrange[1])))
492  c0 = IMP.atom.Chain.setup_particle(IMP.Particle(self.model), "X")
493 
494  for p in ps:
495  par = IMP.atom.Atom(p).get_parent()
496  ri = IMP.atom.Residue(par).get_index()
497  # Move residue from original PDB hierarchy to new chain c0
498  IMP.atom.Residue(par).set_index(ri + offset)
499  if par.get_parent() != c0:
500  par.get_parent().remove_child(par)
501  c0.add_child(par)
502  start = start + offset
503  end = end + offset
504 
505  self.elements[name].append(
506  (start, end, pdbname.split("/")[-1] + ":" + chain, "pdb"))
507 
508  hiers = self.coarse_hierarchy(name, start, end,
509  resolutions, isnucleicacid, c0, protein_h, "pdb", color)
510  self._copy_pdb_provenance(t, hiers[0], offset)
511  outhiers += hiers
512  for p, state in self._protocol_output:
513  p.add_pdb_element(state, name, start, end, offset, pdbname, chain,
514  hiers[0])
515 
516  if show:
518 
519  # We cannot simply destroy(c0) since it might not be a well-behaved
520  # hierarchy; in some cases it could contain a given residue more than
521  # once (this is surely a bug but we need to keep this behavior for
522  # backwards compatibility).
523  residues = {}
524  for p in ps:
525  par = IMP.atom.Atom(p).get_parent()
526  residues[par] = None
527  for r in residues.keys():
528  IMP.atom.destroy(r)
529  self.model.remove_particle(c0)
530 
531  IMP.atom.destroy(t)
532 
533  return outhiers
534 
535  def _copy_pdb_provenance(self, pdb, h, offset):
536  """Copy the provenance information from the PDB to our hierarchy"""
537  c = IMP.atom.Chain(IMP.atom.get_by_type(pdb, IMP.atom.CHAIN_TYPE)[0])
538  prov = IMP.core.Provenanced(c).get_provenance()
539  newprov = IMP.core.create_clone(prov)
542  IMP.core.StructureProvenance(newprov).set_residue_offset(offset)
543 
544  def add_component_ideal_helix(
545  self,
546  name,
547  resolutions,
548  resrange,
549  color=None,
550  show=False):
551 
552  self.representation_is_modified = True
553  from math import pi, cos, sin
554 
555  protein_h = self.hier_dict[name]
556  outhiers = []
557  if color is None:
558  color = self.color_dict[name]
559 
560  start = resrange[0]
561  end = resrange[1]
562  self.elements[name].append((start, end, " ", "helix"))
563  c0 = IMP.atom.Chain.setup_particle(IMP.Particle(self.model), "X")
564  for n, res in enumerate(range(start, end + 1)):
565  if name in self.sequence_dict:
566  try:
567  rtstr = self.onetothree[
568  self.sequence_dict[name][res-1]]
569  except:
570  rtstr = "UNK"
571  rt = IMP.atom.ResidueType(rtstr)
572  else:
573  rt = IMP.atom.ResidueType("ALA")
574 
575  # get the residue volume
576  try:
578  # mass=IMP.atom.get_mass_from_residue_type(rt)
579  except IMP.ValueException:
581  IMP.atom.ResidueType("ALA"))
582  # mass=IMP.atom.get_mass_from_residue_type(IMP.atom.ResidueType("ALA"))
584 
585  r = IMP.atom.Residue.setup_particle(IMP.Particle(self.model), rt, res)
586  p = IMP.Particle(self.model)
588  x = 2.3 * cos(n * 2 * pi / 3.6)
589  y = 2.3 * sin(n * 2 * pi / 3.6)
590  z = 6.2 / 3.6 / 2 * n * 2 * pi / 3.6
591  d.set_coordinates(IMP.algebra.Vector3D(x, y, z))
592  d.set_radius(radius)
593  # print d
594  a = IMP.atom.Atom.setup_particle(p, IMP.atom.AT_CA)
595  r.add_child(a)
596  c0.add_child(r)
597 
598  outhiers += self.coarse_hierarchy(name, start, end,
599  resolutions, False, c0, protein_h, "helix", color)
600 
601  if show:
603 
604  IMP.atom.destroy(c0)
605  return outhiers
606 
607  def add_component_beads(self, name, ds, colors=None, incoord=None):
608  """ add beads to the representation
609  @param name the component name
610  @param ds a list of tuples corresponding to the residue ranges
611  of the beads
612  @param colors a list of colors associated to the beads
613  @param incoord the coordinate tuple corresponding to the position
614  of the beads
615  """
616 
617  from math import pi
618  self.representation_is_modified = True
619 
620  protein_h = self.hier_dict[name]
621  outhiers = []
622  if colors is None:
623  colors = [self.color_dict[name]]
624 
625 
626  for n, dss in enumerate(ds):
627  ds_frag = (dss[0], dss[1])
628  self.elements[name].append((dss[0], dss[1], " ", "bead"))
629  prt = IMP.Particle(self.model)
630  if ds_frag[0] == ds_frag[1]:
631  # if the bead represent a single residue
632  if name in self.sequence_dict:
633  try:
634  rtstr = self.onetothree[
635  self.sequence_dict[name][ds_frag[0]-1]]
636  except:
637  rtstr = "UNK"
638  rt = IMP.atom.ResidueType(rtstr)
639  else:
640  rt = IMP.atom.ResidueType("ALA")
641  h = IMP.atom.Residue.setup_particle(prt, rt, ds_frag[0])
642  h.set_name(name + '_%i_bead' % (ds_frag[0]))
643  prt.set_name(name + '_%i_bead' % (ds_frag[0]))
644  resolution = 1
645  else:
647  h.set_name(name + '_%i-%i_bead' % (ds_frag[0], ds_frag[1]))
648  prt.set_name(name + '_%i-%i_bead' % (ds_frag[0], ds_frag[1]))
649  h.set_residue_indexes(list(range(ds_frag[0], ds_frag[1] + 1)))
650  resolution = len(h.get_residue_indexes())
651  if "Beads" not in self.hier_representation[name]:
653  root.set_name("Beads")
654  self.hier_representation[name]["Beads"] = root
655  protein_h.add_child(root)
656  self.hier_representation[name]["Beads"].add_child(h)
657 
658  for kk in range(ds_frag[0], ds_frag[1] + 1):
659  self.hier_db.add_particles(name, kk, resolution, [h])
660 
661  try:
662  clr = IMP.display.get_rgb_color(colors[n])
663  except:
664  clr = IMP.display.get_rgb_color(colors[0])
665 
667 
668  # decorate particles according to their resolution
669  IMP.pmi1.Resolution.setup_particle(prt, resolution)
670 
672  ptem = IMP.core.XYZR(prt)
674  if resolution == 1:
675  try:
677  except IMP.ValueException:
679  IMP.atom.ResidueType("ALA"))
682  ptem.set_radius(radius)
683  else:
684  volume = IMP.atom.get_volume_from_mass(mass)
685  radius = 0.8 * (3.0 / 4.0 / pi * volume) ** (1.0 / 3.0)
687  ptem.set_radius(radius)
688  try:
689  if not tuple(incoord) is None:
690  ptem.set_coordinates(incoord)
691  except TypeError:
692  pass
695  self.floppy_bodies.append(prt)
696  IMP.core.XYZ(prt).set_coordinates_are_optimized(True)
697  outhiers += [h]
698 
699  for p, state in self._protocol_output:
700  p.add_bead_element(state, name, ds[0][0], ds[-1][1], len(ds),
701  outhiers[0])
702 
703  return outhiers
704 
705  def add_component_necklace(self, name, begin, end, length, color=None,
706  incoord=None):
707  '''
708  Generates a string of beads with given length.
709  '''
710  self.representation_is_modified = True
711  outhiers = []
712  if color is None:
713  colors=None
714  else:
715  colors=[color]
716  for chunk in IMP.pmi1.tools.list_chunks_iterator(range(begin, end + 1), length):
717  outhiers += self.add_component_beads(name,
718  [(chunk[0], chunk[-1])], colors=colors,incoord=incoord)
719  return outhiers
720 
722  self, name, hierarchies=None, selection_tuples=None,
723  particles=None,
724  resolution=0.0, num_components=10,
725  inputfile=None, outputfile=None,
726  outputmap=None,
727  kernel_type=None,
728  covariance_type='full', voxel_size=1.0,
729  out_hier_name='',
730  sampled_points=1000000, num_iter=100,
731  simulation_res=1.0,
732  multiply_by_total_mass=True,
733  transform=None,
734  intermediate_map_fn=None,
735  density_ps_to_copy=None,
736  use_precomputed_gaussians=False):
737  '''
738  Sets up a Gaussian Mixture Model for this component.
739  Can specify input GMM file or it will be computed.
740  @param name component name
741  @param hierarchies set up GMM for some hierarchies
742  @param selection_tuples (list of tuples) example (first_residue,last_residue,component_name)
743  @param particles set up GMM for particles directly
744  @param resolution usual PMI resolution for selecting particles from the hierarchies
745  @param inputfile read the GMM from this file
746  @param outputfile fit and write the GMM to this file (text)
747  @param outputmap after fitting, create GMM density file (mrc)
748  @param kernel_type for creating the intermediate density (points are sampled to make GMM). Options are IMP.em.GAUSSIAN, IMP.em.SPHERE, and IMP.em.BINARIZED_SPHERE
749  @param covariance_type for fitting the GMM. options are 'full', 'diagonal' and 'spherical'
750  @param voxel_size for creating the intermediate density map and output map.
751  lower number increases accuracy but also rasterizing time grows
752  @param out_hier_name name of the output density hierarchy
753  @param sampled_points number of points to sample. more will increase accuracy and fitting time
754  @param num_iter num GMM iterations. more will increase accuracy and fitting time
755  @param multiply_by_total_mass multiply the weights of the GMM by this value (only works on creation!)
756  @param transform for input file only, apply a transformation (eg for multiple copies same GMM)
757  @param intermediate_map_fn for debugging, this will write the intermediate (simulated) map
758  @param density_ps_to_copy in case you already created the appropriate GMM (eg, for beads)
759  @param use_precomputed_gaussians Set this flag and pass fragments - will use roughly spherical Gaussian setup
760  '''
761  import numpy as np
762  import sys
763  import IMP.em
764  import IMP.isd.gmm_tools
765 
766  # prepare output
767  self.representation_is_modified = True
768  out_hier = []
769  protein_h = self.hier_dict[name]
770  if "Densities" not in self.hier_representation[name]:
772  root.set_name("Densities")
773  self.hier_representation[name]["Densities"] = root
774  protein_h.add_child(root)
775 
776  # gather passed particles
777  fragment_particles = []
778  if not particles is None:
779  fragment_particles += particles
780  if not hierarchies is None:
781  fragment_particles += IMP.pmi1.tools.select(
782  self, resolution=resolution,
783  hierarchies=hierarchies)
784  if not selection_tuples is None:
785  for st in selection_tuples:
786  fragment_particles += IMP.pmi1.tools.select_by_tuple(
787  self, tupleselection=st,
788  resolution=resolution,
789  name_is_ambiguous=False)
790 
791  # compute or read gaussians
792  density_particles = []
793  if inputfile:
795  inputfile, density_particles,
796  self.model, transform)
797  elif density_ps_to_copy:
798  for ip in density_ps_to_copy:
799  p = IMP.Particle(self.model)
800  shape = IMP.core.Gaussian(ip).get_gaussian()
801  mass = IMP.atom.Mass(ip).get_mass()
804  density_particles.append(p)
805  elif use_precomputed_gaussians:
806  if len(fragment_particles) == 0:
807  print("add_component_density: no particle was selected")
808  return out_hier
809  for p in fragment_particles:
810  if not (IMP.atom.Fragment.get_is_setup(self.model,p.get_particle_index()) and
811  IMP.core.XYZ.get_is_setup(self.model,p.get_particle_index())):
812  raise Exception("The particles you selected must be Fragments and XYZs")
813  nres=len(IMP.atom.Fragment(self.model,p.get_particle_index()).get_residue_indexes())
814  pos=IMP.core.XYZ(self.model,p.get_particle_index()).get_coordinates()
815  density_particles=[]
816  try:
817  IMP.isd.get_data_path("beads/bead_%i.txt"%nres)
818  except:
819  raise Exception("We haven't created a bead file for",nres,"residues")
820  transform = IMP.algebra.Transformation3D(pos)
822  IMP.isd.get_data_path("beads/bead_%i.txt"%nres), density_particles,
823  self.model, transform)
824  else:
825  #compute the gaussians here
826  if len(fragment_particles) == 0:
827  print("add_component_density: no particle was selected")
828  return out_hier
829 
830  density_particles = IMP.isd.gmm_tools.sample_and_fit_to_particles(
831  self.model,
832  fragment_particles,
833  num_components,
834  sampled_points,
835  simulation_res,
836  voxel_size,
837  num_iter,
838  covariance_type,
839  multiply_by_total_mass,
840  outputmap,
841  outputfile)
842 
843  # prepare output hierarchy
845  s0.set_name(out_hier_name)
846  self.hier_representation[name]["Densities"].add_child(s0)
847  out_hier.append(s0)
848  for nps, p in enumerate(density_particles):
849  s0.add_child(p)
850  p.set_name(s0.get_name() + '_gaussian_%i' % nps)
851  return out_hier
852 
853  def get_component_density(self, name):
854  return self.hier_representation[name]["Densities"]
855 
856  def add_all_atom_densities(self, name, hierarchies=None,
857  selection_tuples=None,
858  particles=None,
859  resolution=0,
860  output_map=None,
861  voxel_size=1.0):
862  '''Decorates all specified particles as Gaussians directly.
863  @param name component name
864  @param hierarchies set up GMM for some hierarchies
865  @param selection_tuples (list of tuples) example (first_residue,last_residue,component_name)
866  @param particles set up GMM for particles directly
867  @param resolution usual PMI resolution for selecting particles from the hierarchies
868  '''
869 
870  import IMP.em
871  import numpy as np
872  import sys
873  from math import sqrt
874  self.representation_is_modified = True
875 
876  if particles is None:
877  fragment_particles = []
878  else:
879  fragment_particles = particles
880 
881  if not hierarchies is None:
882  fragment_particles += IMP.pmi1.tools.select(
883  self, resolution=resolution,
884  hierarchies=hierarchies)
885 
886  if not selection_tuples is None:
887  for st in selection_tuples:
888  fragment_particles += IMP.pmi1.tools.select_by_tuple(
889  self, tupleselection=st,
890  resolution=resolution,
891  name_is_ambiguous=False)
892 
893  if len(fragment_particles) == 0:
894  print("add all atom densities: no particle was selected")
895  return
896 
897  # create a spherical gaussian for each particle based on atom type
898  print('setting up all atom gaussians num_particles',len(fragment_particles))
899  for n,p in enumerate(fragment_particles):
900  if IMP.core.Gaussian.get_is_setup(p): continue
901  center=IMP.core.XYZ(p).get_coordinates()
902  rad=IMP.core.XYZR(p).get_radius()
903  mass=IMP.atom.Mass(p).get_mass()
907  print('setting up particle',p.get_name(), " as individual gaussian particle")
908 
909  if not output_map is None:
910  print('writing map to', output_map)
912  fragment_particles,
913  output_map,
914  voxel_size)
915 
916  def add_component_hierarchy_clone(self, name, hierarchy):
917  '''
918  Make a copy of a hierarchy and append it to a component.
919  '''
920  outhier = []
921  self.representation_is_modified = True
922  protein_h = self.hier_dict[name]
923  hierclone = IMP.atom.create_clone(hierarchy)
924  hierclone.set_name(hierclone.get_name() + "_clone")
925  protein_h.add_child(hierclone)
926  outhier.append(hierclone)
927 
928  psmain = IMP.atom.get_leaves(hierarchy)
929  psclone = IMP.atom.get_leaves(hierclone)
930 
931  # copying attributes
932  for n, pmain in enumerate(psmain):
933  pclone = psclone[n]
935  resolution = IMP.pmi1.Resolution(pmain).get_resolution()
936  IMP.pmi1.Resolution.setup_particle(pclone, resolution)
937  for kk in IMP.pmi1.tools.get_residue_indexes(pclone):
938  self.hier_db.add_particles(
939  name,
940  kk,
942  [pclone])
943 
945  uncertainty = IMP.pmi1.Uncertainty(pmain).get_uncertainty()
946  IMP.pmi1.Uncertainty.setup_particle(pclone, uncertainty)
947 
949  symmetric = IMP.pmi1.Symmetric(pmain).get_symmetric()
950  IMP.pmi1.Symmetric.setup_particle(pclone, symmetric)
951 
952  return outhier
953 
954  def dump_particle_descriptors(self):
955  import numpy
956  import pickle
957  import IMP.isd
958  import IMP.isd.gmm_tools
959 
960  particles_attributes={}
961  floppy_body_attributes={}
962  gaussians=[]
963  for h in IMP.atom.get_leaves(self.prot):
964  leaf=h
965  name=h.get_name()
966  hroot=self.prot
967  hparent=h.get_parent()
968  while hparent != hroot:
969  hparent=h.get_parent()
970  name+="|"+hparent.get_name()
971  h=hparent
972  particles_attributes[name]={"COORDINATES":numpy.array(IMP.core.XYZR(leaf.get_particle()).get_coordinates()),
973  "RADIUS":IMP.core.XYZR(leaf.get_particle()).get_radius(),
974  "MASS":IMP.atom.Mass(leaf.get_particle()).get_mass()}
975  if IMP.core.Gaussian.get_is_setup(leaf.get_particle()):
976  gaussians.append(IMP.core.Gaussian(leaf.get_particle()))
977 
978  rigid_body_attributes={}
979  for rb in self.rigid_bodies:
980  name=rb.get_name()
981  rf=rb.get_reference_frame()
982  t=rf.get_transformation_to()
983  trans=t.get_translation()
984  rot=t.get_rotation()
985  rigid_body_attributes[name]={"TRANSLATION":numpy.array(trans),
986  "ROTATION":numpy.array(rot.get_quaternion()),
987  "COORDINATES_NONRIGID_MEMBER":{},
988  "COORDINATES_RIGID_MEMBER":{}}
989  for mi in rb.get_member_indexes():
990  rm=self.model.get_particle(mi)
992  name_part=rm.get_name()
993  xyz=[self.model.get_attribute(fk, rm) for fk in [IMP.FloatKey(4), IMP.FloatKey(5), IMP.FloatKey(6)]]
994  rigid_body_attributes[name]["COORDINATES_NONRIGID_MEMBER"][name_part]=numpy.array(xyz)
995  else:
996  name_part=rm.get_name()
997  xyz=IMP.core.XYZ(rm).get_coordinates()
998  rigid_body_attributes[name]["COORDINATES_RIGID_MEMBER"][name_part]=numpy.array(xyz)
999 
1000 
1001  IMP.isd.gmm_tools.write_gmm_to_text(gaussians,"model_gmm.txt")
1002  pickle.dump(particles_attributes,
1003  open("particles_attributes.pkl", "wb"))
1004  pickle.dump(rigid_body_attributes,
1005  open("rigid_body_attributes.pkl", "wb"))
1006 
1007 
1008 
1009  def load_particle_descriptors(self):
1010  import numpy
1011  import pickle
1012  import IMP.isd
1013  import IMP.isd.gmm_tools
1014 
1015  particles_attributes = pickle.load(open("particles_attributes.pkl",
1016  "rb"))
1017  rigid_body_attributes = pickle.load(open("rigid_body_attributes.pkl",
1018  "rb"))
1019 
1020  particles=[]
1021  hierarchies=[]
1022  gaussians=[]
1023  for h in IMP.atom.get_leaves(self.prot):
1024  leaf=h
1025  name=h.get_name()
1026  hroot=self.prot
1027  hparent=h.get_parent()
1028  while hparent != hroot:
1029  hparent=h.get_parent()
1030  name+="|"+hparent.get_name()
1031  h=hparent
1032 
1033  xyzr=IMP.core.XYZR(leaf.get_particle())
1034  xyzr.set_coordinates(particles_attributes[name]["COORDINATES"])
1035  #xyzr.set_radius(particles_attributes[name]["RADIUS"])
1036  #IMP.atom.Mass(leaf.get_particle()).set_mass(particles_attributes[name]["MASS"])
1037  if IMP.core.Gaussian.get_is_setup(leaf.get_particle()):
1038  gaussians.append(IMP.core.Gaussian(leaf.get_particle()))
1039 
1040  for rb in self.rigid_bodies:
1041  name=rb.get_name()
1042  trans=rigid_body_attributes[name]["TRANSLATION"]
1043  rot=rigid_body_attributes[name]["ROTATION"]
1046  rb.set_reference_frame(rf)
1047  coor_nrm_ref=rigid_body_attributes[name]["COORDINATES_NONRIGID_MEMBER"]
1048  coor_rm_ref_dict=rigid_body_attributes[name]["COORDINATES_RIGID_MEMBER"]
1049  coor_rm_model=[]
1050  coor_rm_ref=[]
1051  for mi in rb.get_member_indexes():
1052  rm=self.model.get_particle(mi)
1054  name_part=rm.get_name()
1055  xyz=coor_nrm_ref[name_part]
1056  for n,fk in enumerate([IMP.FloatKey(4), IMP.FloatKey(5), IMP.FloatKey(6)]):
1057  self.model.set_attribute(fk, rm,xyz[n])
1058  else:
1059  name_part=rm.get_name()
1060  coor_rm_ref.append(IMP.algebra.Vector3D(coor_rm_ref_dict[name_part]))
1061  coor_rm_model.append(IMP.core.XYZ(rm).get_coordinates())
1062  if len(coor_rm_model)==0: continue
1063  t=IMP.algebra.get_transformation_aligning_first_to_second(coor_rm_model,coor_rm_ref)
1064  IMP.core.transform(rb,t)
1065 
1066  IMP.isd.gmm_tools.decorate_gmm_from_text("model_gmm.txt",gaussians,self.model)
1067 
1068  def _compare_rmf_repr_names(self, rmfname, reprname, component_name):
1069  """Print a warning if particle names in RMF and model don't match"""
1070  def match_any_suffix():
1071  # Handle common mismatches like 743 != Nup85_743_pdb
1072  suffixes = ["pdb", "bead_floppy_body_rigid_body_member_floppy_body",
1073  "bead_floppy_body_rigid_body_member",
1074  "bead_floppy_body"]
1075  for s in suffixes:
1076  if "%s_%s_%s" % (component_name, rmfname, s) == reprname:
1077  return True
1078  if rmfname != reprname and not match_any_suffix():
1079  print("set_coordinates_from_rmf: WARNING rmf particle and "
1080  "representation particle names don't match %s %s"
1081  % (rmfname, reprname))
1082 
1083  def set_coordinates_from_rmf(self, component_name, rmf_file_name,
1084  rmf_frame_number,
1085  rmf_component_name=None,
1086  check_number_particles=True,
1087  representation_name_to_rmf_name_map=None,
1088  state_number=0,
1089  skip_gaussian_in_rmf=False,
1090  skip_gaussian_in_representation=False,
1091  save_file=False,
1092  force_rigid_update=False):
1093  '''Read and replace coordinates from an RMF file.
1094  Replace the coordinates of particles with the same name.
1095  It assumes that the RMF and the representation have the particles
1096  in the same order.
1097  @param component_name Component name
1098  @param rmf_component_name Name of the component in the RMF file
1099  (if not specified, use `component_name`)
1100  @param representation_name_to_rmf_name_map a dictionary that map
1101  the original rmf particle name to the recipient particle component name
1102  @param save_file: save a file with the names of particles of the component
1103  @param force_rigid_update: update the coordinates of rigid bodies
1104  (normally this should be called before rigid bodies are set up)
1105  '''
1106  import IMP.pmi1.analysis
1107 
1108  tempm = IMP.Model()
1109  prots = IMP.pmi1.analysis.get_hiers_from_rmf(
1110  tempm,
1111  rmf_frame_number,
1112  rmf_file_name)
1113 
1114  if not prots:
1115  raise ValueError("cannot read hierarchy from rmf")
1116 
1117  prot=prots[0]
1118  # Make sure coordinates of rigid body members in the RMF are correct
1119  if force_rigid_update:
1120  tempm.update()
1121 
1122  # if len(self.rigid_bodies)!=0:
1123  # print "set_coordinates_from_rmf: cannot proceed if rigid bodies were initialized. Use the function before defining the rigid bodies"
1124  # exit()
1125 
1126  allpsrmf = IMP.atom.get_leaves(prot)
1127  psrmf = []
1128  for p in allpsrmf:
1129  (protname, is_a_bead) = IMP.pmi1.tools.get_prot_name_from_particle(
1130  p, self.hier_dict.keys())
1131  if (protname is None) and (rmf_component_name is not None):
1132  (protname, is_a_bead) = IMP.pmi1.tools.get_prot_name_from_particle(
1133  p, rmf_component_name)
1134  if (skip_gaussian_in_rmf):
1136  continue
1137  if (rmf_component_name is not None) and (protname == rmf_component_name):
1138  psrmf.append(p)
1139  elif (rmf_component_name is None) and (protname == component_name):
1140  psrmf.append(p)
1141 
1142  psrepr = IMP.atom.get_leaves(self.hier_dict[component_name])
1143  if (skip_gaussian_in_representation):
1144  allpsrepr = psrepr
1145  psrepr = []
1146  for p in allpsrepr:
1147  #(protname, is_a_bead) = IMP.pmi1.tools.get_prot_name_from_particle(
1148  # p, self.hier_dict.keys())
1150  continue
1151  psrepr.append(p)
1152 
1153  import itertools
1154  reprnames=[p.get_name() for p in psrepr]
1155  rmfnames=[p.get_name() for p in psrmf]
1156 
1157  if save_file:
1158  fl=open(component_name+".txt","w")
1159  for i in itertools.izip_longest(reprnames,rmfnames): fl.write(str(i[0])+","+str(i[1])+"\n")
1160 
1161 
1162  if check_number_particles and not representation_name_to_rmf_name_map:
1163  if len(psrmf) != len(psrepr):
1164  fl=open(component_name+".txt","w")
1165  for i in itertools.izip_longest(reprnames,rmfnames): fl.write(str(i[0])+","+str(i[1])+"\n")
1166  raise ValueError("%s cannot proceed the rmf and the representation don't have the same number of particles; "
1167  "particles in rmf: %s particles in the representation: %s" % (str(component_name), str(len(psrmf)), str(len(psrepr))))
1168 
1169 
1170  if not representation_name_to_rmf_name_map:
1171  for n, prmf in enumerate(psrmf):
1172 
1173  prmfname = prmf.get_name()
1174  preprname = psrepr[n].get_name()
1175  if force_rigid_update:
1176  if IMP.core.RigidBody.get_is_setup(psrepr[n]) \
1177  and not IMP.core.RigidBodyMember.get_is_setup(psrepr[n]):
1178  continue
1179  else:
1181  raise ValueError("component %s cannot proceed if rigid bodies were initialized. Use the function before defining the rigid bodies" % component_name)
1182 
1183  self._compare_rmf_repr_names(prmfname, preprname,
1184  component_name)
1185  if IMP.core.XYZ.get_is_setup(prmf) and IMP.core.XYZ.get_is_setup(psrepr[n]):
1186  xyz = IMP.core.XYZ(prmf).get_coordinates()
1187  IMP.core.XYZ(psrepr[n]).set_coordinates(xyz)
1189  # Set rigid body so that coordinates are preserved
1190  # on future model updates
1191  rbm = IMP.core.RigidBodyMember(psrepr[n])
1192  rbm.set_internal_coordinates(xyz)
1194  rb = rbm.get_rigid_body()
1196  raise ValueError("Cannot handle nested "
1197  "rigid bodies yet")
1198  rb.set_reference_frame_lazy(tr)
1199  else:
1200  print("set_coordinates_from_rmf: WARNING particles are not XYZ decorated %s %s " % (str(IMP.core.XYZ.get_is_setup(prmf)), str(IMP.core.XYZ.get_is_setup(psrepr[n]))))
1201 
1203  gprmf=IMP.core.Gaussian(prmf)
1204  grepr=IMP.core.Gaussian(psrepr[n])
1205  g=gprmf.get_gaussian()
1206  grepr.set_gaussian(g)
1207 
1208  else:
1209  repr_name_particle_map={}
1210  rmf_name_particle_map={}
1211  for p in psrmf:
1212  rmf_name_particle_map[p.get_name()]=p
1213  #for p in psrepr:
1214  # repr_name_particle_map[p.get_name()]=p
1215 
1216  for prepr in psrepr:
1217  try:
1218  prmf=rmf_name_particle_map[representation_name_to_rmf_name_map[prepr.get_name()]]
1219  except KeyError:
1220  print("set_coordinates_from_rmf: WARNING missing particle name in representation_name_to_rmf_name_map, skipping...")
1221  continue
1222  xyz = IMP.core.XYZ(prmf).get_coordinates()
1223  IMP.core.XYZ(prepr).set_coordinates(xyz)
1224 
1225 
1226 
1227  def check_root(self, name, protein_h, resolution):
1228  '''
1229  If the root hierarchy does not exist, construct it.
1230  '''
1231  if "Res:" + str(int(resolution)) not in self.hier_representation[name]:
1233  root.set_name(name + "_Res:" + str(int(resolution)))
1234  self.hier_representation[name][
1235  "Res:" + str(int(resolution))] = root
1236  protein_h.add_child(root)
1237 
1238  def coarse_hierarchy(self, name, start, end, resolutions, isnucleicacid,
1239  input_hierarchy, protein_h, type, color):
1240  '''
1241  Generate all needed coarse grained layers.
1242 
1243  @param name name of the protein
1244  @param resolutions list of resolutions
1245  @param protein_h root hierarchy
1246  @param input_hierarchy hierarchy to coarse grain
1247  @param type a string, typically "pdb" or "helix"
1248  '''
1249  outhiers = []
1250 
1251  if (1 in resolutions) or (0 in resolutions):
1252  # in that case create residues and append atoms
1253 
1254  if 1 in resolutions:
1255  self.check_root(name, protein_h, 1)
1257  s1.set_name('%s_%i-%i_%s' % (name, start, end, type))
1258  # s1.set_residue_indexes(range(start,end+1))
1259  self.hier_representation[name]["Res:1"].add_child(s1)
1260  outhiers += [s1]
1261  if 0 in resolutions:
1262  self.check_root(name, protein_h, 0)
1264  s0.set_name('%s_%i-%i_%s' % (name, start, end, type))
1265  # s0.set_residue_indexes(range(start,end+1))
1266  self.hier_representation[name]["Res:0"].add_child(s0)
1267  outhiers += [s0]
1268 
1269  if not isnucleicacid:
1270  sel = IMP.atom.Selection(
1271  input_hierarchy,
1272  atom_type=IMP.atom.AT_CA)
1273  else:
1274  sel = IMP.atom.Selection(
1275  input_hierarchy,
1276  atom_type=IMP.atom.AT_P)
1277 
1278  for p in sel.get_selected_particles():
1279  resobject = IMP.atom.Residue(IMP.atom.Atom(p).get_parent())
1280  if 0 in resolutions:
1281  # if you ask for atoms
1282  resclone0 = IMP.atom.create_clone(resobject)
1283  resindex = IMP.atom.Residue(resclone0).get_index()
1284  s0.add_child(resclone0)
1285  self.hier_db.add_particles(
1286  name,
1287  resindex,
1288  0,
1289  resclone0.get_children())
1290 
1291  chil = resclone0.get_children()
1292  for ch in chil:
1294  try:
1295  clr = IMP.display.get_rgb_color(color)
1296  except:
1297  clr = IMP.display.get_rgb_color(1.0)
1299 
1300  if 1 in resolutions:
1301  # else clone the residue
1302  resclone1 = IMP.atom.create_clone_one(resobject)
1303  resindex = IMP.atom.Residue(resclone1).get_index()
1304  s1.add_child(resclone1)
1305  self.hier_db.add_particles(name, resindex, 1, [resclone1])
1306 
1307  rt = IMP.atom.Residue(resclone1).get_residue_type()
1308  xyz = IMP.core.XYZ(p).get_coordinates()
1309  prt = resclone1.get_particle()
1310  prt.set_name('%s_%i_%s' % (name, resindex, type))
1311  IMP.core.XYZ.setup_particle(prt).set_coordinates(xyz)
1312 
1313  try:
1315  # mass=IMP.atom.get_mass_from_residue_type(rt)
1316  except IMP.ValueException:
1318  IMP.atom.ResidueType("ALA"))
1319  # mass=IMP.atom.get_mass_from_residue_type(IMP.atom.ResidueType("ALA"))
1321  IMP.core.XYZR.setup_particle(prt).set_radius(radius)
1326  try:
1327  clr = IMP.display.get_rgb_color(color)
1328  except:
1329  clr = IMP.display.get_rgb_color(1.0)
1331 
1332  for r in resolutions:
1333  if r != 0 and r != 1:
1334  self.check_root(name, protein_h, r)
1336  input_hierarchy,
1337  r)
1338 
1339  chil = s.get_children()
1341  s0.set_name('%s_%i-%i_%s' % (name, start, end, type))
1342  # Move all children from s to s0
1343  for ch in chil:
1344  s.remove_child(ch)
1345  s0.add_child(ch)
1346  self.hier_representation[name][
1347  "Res:" + str(int(r))].add_child(s0)
1348  outhiers += [s0]
1349  IMP.atom.destroy(s)
1350 
1351  for prt in IMP.atom.get_leaves(s0):
1352  ri = IMP.atom.Fragment(prt).get_residue_indexes()
1353  first = ri[0]
1354  last = ri[-1]
1355  if first == last:
1356  prt.set_name('%s_%i_%s' % (name, first, type))
1357  else:
1358  prt.set_name('%s_%i_%i_%s' % (name, first, last, type))
1359  for kk in ri:
1360  self.hier_db.add_particles(name, kk, r, [prt])
1361 
1362  radius = IMP.core.XYZR(prt).get_radius()
1366 
1367  # setting up color for each particle in the
1368  # hierarchy, if colors missing in the colors list set it to
1369  # red
1370  try:
1371  clr = IMP.display.get_rgb_color(color)
1372  except:
1373  colors.append(1.0)
1374  clr = IMP.display.get_rgb_color(colors[pdb_part_count])
1376 
1377  return outhiers
1378 
1380  '''
1381  Get the hierarchies at the given resolution.
1382 
1383  The map between resolution and hierarchies is cached to accelerate
1384  the selection algorithm. The cache is invalidated when the
1385  representation was changed by any add_component_xxx.
1386  '''
1387 
1388  if self.representation_is_modified:
1389  rhbr = self.hier_db.get_all_root_hierarchies_by_resolution(
1390  resolution)
1391  self.hier_resolution[resolution] = rhbr
1392  self.representation_is_modified = False
1393  return rhbr
1394  else:
1395  if resolution in self.hier_resolution:
1396  return self.hier_resolution[resolution]
1397  else:
1398  rhbr = self.hier_db.get_all_root_hierarchies_by_resolution(
1399  resolution)
1400  self.hier_resolution[resolution] = rhbr
1401  return rhbr
1402 
1404  self, max_translation=300., max_rotation=2.0 * pi,
1405  avoidcollision=True, cutoff=10.0, niterations=100,
1406  bounding_box=None,
1407  excluded_rigid_bodies=None,
1408  ignore_initial_coordinates=False,
1409  hierarchies_excluded_from_collision=None):
1410  '''
1411  Shuffle configuration; used to restart the optimization.
1412 
1413  The configuration of the system is initialized by placing each
1414  rigid body and each bead randomly in a box with a side of
1415  `max_translation` angstroms, and far enough from each other to
1416  prevent any steric clashes. The rigid bodies are also randomly rotated.
1417 
1418  @param avoidcollision check if the particle/rigid body was
1419  placed close to another particle; uses the optional
1420  arguments cutoff and niterations
1421  @param bounding_box defined by ((x1,y1,z1),(x2,y2,z2))
1422  '''
1423 
1424  if excluded_rigid_bodies is None:
1425  excluded_rigid_bodies = []
1426  if hierarchies_excluded_from_collision is None:
1427  hierarchies_excluded_from_collision = []
1428 
1429  if len(self.rigid_bodies) == 0:
1430  print("shuffle_configuration: rigid bodies were not initialized")
1431 
1433  gcpf.set_distance(cutoff)
1434  ps = []
1435  hierarchies_excluded_from_collision_indexes = []
1436  for p in IMP.atom.get_leaves(self.prot):
1438  ps.append(p)
1440  # remove the densities particles out of the calculation
1441  hierarchies_excluded_from_collision_indexes.extend(
1442  IMP.get_indexes([p]))
1443  allparticleindexes = IMP.get_indexes(ps)
1444 
1445  if bounding_box is not None:
1446  ((x1, y1, z1), (x2, y2, z2)) = bounding_box
1447  lb = IMP.algebra.Vector3D(x1, y1, z1)
1448  ub = IMP.algebra.Vector3D(x2, y2, z2)
1449  bb = IMP.algebra.BoundingBox3D(lb, ub)
1450 
1451  for h in hierarchies_excluded_from_collision:
1452  hierarchies_excluded_from_collision_indexes.extend(
1454 
1455 
1456  allparticleindexes = list(
1457  set(allparticleindexes) - set(hierarchies_excluded_from_collision_indexes))
1458 
1459  print(hierarchies_excluded_from_collision)
1460  print(len(allparticleindexes),len(hierarchies_excluded_from_collision_indexes))
1461 
1462  print('shuffling', len(self.rigid_bodies) - len(excluded_rigid_bodies), 'rigid bodies')
1463  for rb in self.rigid_bodies:
1464  if rb not in excluded_rigid_bodies:
1465  if avoidcollision:
1466  rbindexes = rb.get_member_particle_indexes()
1467  rbindexes = list(
1468  set(rbindexes) - set(hierarchies_excluded_from_collision_indexes))
1469  otherparticleindexes = list(
1470  set(allparticleindexes) - set(rbindexes))
1471 
1472  if len(otherparticleindexes) is None:
1473  continue
1474 
1475  niter = 0
1476  while niter < niterations:
1477  if (ignore_initial_coordinates):
1478  # Move the particle to the origin
1479  transformation = IMP.algebra.Transformation3D(IMP.algebra.get_identity_rotation_3d(), -IMP.core.XYZ(rb).get_coordinates())
1480  IMP.core.transform(rb, transformation)
1481  rbxyz = IMP.core.XYZ(rb).get_coordinates()
1482 
1483  if bounding_box is not None:
1484  # overrides the perturbation
1485  translation = IMP.algebra.get_random_vector_in(bb)
1487  transformation = IMP.algebra.Transformation3D(rotation, translation-rbxyz)
1488  else:
1490  rbxyz,
1491  max_translation,
1492  max_rotation)
1493 
1494  IMP.core.transform(rb, transformation)
1495 
1496  if avoidcollision:
1497  self.model.update()
1498  npairs = len(
1499  gcpf.get_close_pairs(
1500  self.model,
1501  otherparticleindexes,
1502  rbindexes))
1503  if npairs == 0:
1504  niter = niterations
1505  if (ignore_initial_coordinates):
1506  print (rb.get_name(), IMP.core.XYZ(rb).get_coordinates())
1507  else:
1508  niter += 1
1509  print("shuffle_configuration: rigid body placed close to other %d particles, trying again..." % npairs)
1510  print("shuffle_configuration: rigid body name: " + rb.get_name())
1511  if niter == niterations:
1512  raise ValueError("tried the maximum number of iterations to avoid collisions, increase the distance cutoff")
1513  else:
1514  break
1515 
1516  print('shuffling', len(self.floppy_bodies), 'floppy bodies')
1517  for fb in self.floppy_bodies:
1518  if (avoidcollision):
1521  if rm and nrm:
1522  fbindexes = IMP.get_indexes([fb])
1523  otherparticleindexes = list(
1524  set(allparticleindexes) - set(fbindexes))
1525  if len(otherparticleindexes) is None:
1526  continue
1527  else:
1528  continue
1529  else:
1532  if (rm or nrm):
1533  continue
1534 
1536  d=IMP.core.RigidBodyMember(fb).get_rigid_body()
1538  d=IMP.core.RigidBody(fb)
1539  elif IMP.core.XYZ.get_is_setup(fb):
1540  d=IMP.core.XYZ(fb)
1541 
1542  niter = 0
1543  while niter < niterations:
1544  if (ignore_initial_coordinates):
1545  # Move the particle to the origin
1546  transformation = IMP.algebra.Transformation3D(IMP.algebra.get_identity_rotation_3d(), -IMP.core.XYZ(fb).get_coordinates())
1547  IMP.core.transform(d, transformation)
1548  fbxyz = IMP.core.XYZ(fb).get_coordinates()
1549 
1550  if bounding_box is not None:
1551  # overrides the perturbation
1552  translation = IMP.algebra.get_random_vector_in(bb)
1553  transformation = IMP.algebra.Transformation3D(translation-fbxyz)
1554  else:
1556  fbxyz,
1557  max_translation,
1558  max_rotation)
1559 
1560  IMP.core.transform(d, transformation)
1561 
1562  if (avoidcollision):
1563  self.model.update()
1564  npairs = len(
1565  gcpf.get_close_pairs(
1566  self.model,
1567  otherparticleindexes,
1568  fbindexes))
1569  if npairs == 0:
1570  niter = niterations
1571  if (ignore_initial_coordinates):
1572  print (fb.get_name(), IMP.core.XYZ(fb).get_coordinates())
1573  else:
1574  niter += 1
1575  print("shuffle_configuration: floppy body placed close to other %d particles, trying again..." % npairs)
1576  print("shuffle_configuration: floppy body name: " + fb.get_name())
1577  if niter == niterations:
1578  raise ValueError("tried the maximum number of iterations to avoid collisions, increase the distance cutoff")
1579  else:
1580  break
1581 
1582  def set_current_coordinates_as_reference_for_rmsd(self, label="None"):
1583  # getting only coordinates from pdb
1584  ps = IMP.pmi1.tools.select(self, resolution=1.0)
1585  # storing the reference coordinates and the particles
1586  self.reference_structures[label] = (
1587  [IMP.core.XYZ(p).get_coordinates() for p in ps],
1588  ps)
1589 
1590  def get_all_rmsds(self):
1591  rmsds = {}
1592 
1593  for label in self.reference_structures:
1594 
1595  current_coordinates = [IMP.core.XYZ(p).get_coordinates()
1596  for p in self.reference_structures[label][1]]
1597  reference_coordinates = self.reference_structures[label][0]
1598  if len(reference_coordinates) != len(current_coordinates):
1599  print("calculate_all_rmsds: reference and actual coordinates are not the same")
1600  continue
1602  current_coordinates,
1603  reference_coordinates)
1604  rmsd_global = IMP.algebra.get_rmsd(
1605  reference_coordinates,
1606  current_coordinates)
1607  # warning: temporary we are calculating the drms, and not the rmsd,
1608  # for the relative distance
1609  rmsd_relative = IMP.atom.get_drms(
1610  reference_coordinates,
1611  current_coordinates)
1612  rmsds[label + "_GlobalRMSD"] = rmsd_global
1613  rmsds[label + "_RelativeDRMS"] = rmsd_relative
1614  return rmsds
1615 
1616  def setup_component_geometry(self, name, color=None, resolution=1.0):
1617  if color is None:
1618  color = self.color_dict[name]
1619  # this function stores all particle pairs
1620  # ordered by residue number, to be used
1621  # to construct backbone traces
1622  self.hier_geometry_pairs[name] = []
1623  protein_h = self.hier_dict[name]
1624  pbr = IMP.pmi1.tools.select(self, name=name, resolution=resolution)
1625  pbr = IMP.pmi1.tools.sort_by_residues(pbr)
1626 
1627  for n in range(len(pbr) - 1):
1628  self.hier_geometry_pairs[name].append((pbr[n], pbr[n + 1], color))
1629 
1631  self, name, resolution=10, scale=1.0):
1632  '''
1633  Generate restraints between contiguous fragments in the hierarchy.
1634  The linkers are generated at resolution 10 by default.
1635 
1636  '''
1637 
1638  unmodeledregions_cr = IMP.RestraintSet(self.model, "unmodeledregions")
1639  sortedsegments_cr = IMP.RestraintSet(self.model, "sortedsegments")
1640 
1641  protein_h = self.hier_dict[name]
1642  SortedSegments = []
1643  frs = self.hier_db.get_preroot_fragments_by_resolution(
1644  name,
1645  resolution)
1646 
1647  for fr in frs:
1648  try:
1649  start = fr.get_children()[0]
1650  except:
1651  start = fr
1652 
1653  try:
1654  end = fr.get_children()[-1]
1655  except:
1656  end = fr
1657 
1658  startres = IMP.pmi1.tools.get_residue_indexes(start)[0]
1659  endres = IMP.pmi1.tools.get_residue_indexes(end)[-1]
1660  SortedSegments.append((start, end, startres))
1661  SortedSegments = sorted(SortedSegments, key=itemgetter(2))
1662 
1663  # connect the particles
1664  for x in range(len(SortedSegments) - 1):
1665  last = SortedSegments[x][1]
1666  first = SortedSegments[x + 1][0]
1667 
1668  nreslast = len(IMP.pmi1.tools.get_residue_indexes(last))
1669  lastresn = IMP.pmi1.tools.get_residue_indexes(last)[-1]
1670  nresfirst = len(IMP.pmi1.tools.get_residue_indexes(first))
1671  firstresn = IMP.pmi1.tools.get_residue_indexes(first)[0]
1672 
1673  residuegap = firstresn - lastresn - 1
1674  if self.disorderedlength and (nreslast / 2 + nresfirst / 2 + residuegap) > 20.0:
1675  # calculate the distance between the sphere centers using Kohn
1676  # PNAS 2004
1677  optdist = sqrt(5 / 3) * 1.93 * \
1678  (nreslast / 2 + nresfirst / 2 + residuegap) ** 0.6
1679  # optdist2=sqrt(5/3)*1.93*((nreslast)**0.6+(nresfirst)**0.6)/2
1680  if self.upperharmonic:
1681  hu = IMP.core.HarmonicUpperBound(optdist, self.kappa)
1682  else:
1683  hu = IMP.core.Harmonic(optdist, self.kappa)
1684  dps = IMP.core.DistancePairScore(hu)
1685  else: # default
1686  optdist = (0.0 + (float(residuegap) + 1.0) * 3.6) * scale
1687  if self.upperharmonic: # default
1688  hu = IMP.core.HarmonicUpperBound(optdist, self.kappa)
1689  else:
1690  hu = IMP.core.Harmonic(optdist, self.kappa)
1692 
1693  pt0 = last.get_particle()
1694  pt1 = first.get_particle()
1695  r = IMP.core.PairRestraint(self.model, dps, (pt0.get_index(), pt1.get_index()))
1696 
1697  print("Adding sequence connectivity restraint between", pt0.get_name(), " and ", pt1.get_name(), 'of distance', optdist)
1698  sortedsegments_cr.add_restraint(r)
1699  self.linker_restraints_dict[
1700  "LinkerRestraint-" + pt0.get_name() + "-" + pt1.get_name()] = r
1701  self.connected_intra_pairs.append((pt0, pt1))
1702  self.connected_intra_pairs.append((pt1, pt0))
1703 
1704  self.linker_restraints.add_restraint(sortedsegments_cr)
1705  self.linker_restraints.add_restraint(unmodeledregions_cr)
1706  IMP.pmi1.tools.add_restraint_to_model(self.model, self.linker_restraints)
1707  self.sortedsegments_cr_dict[name] = sortedsegments_cr
1708  self.unmodeledregions_cr_dict[name] = unmodeledregions_cr
1709 
1710  def optimize_floppy_bodies(self, nsteps, temperature=1.0):
1711  import IMP.pmi1.samplers
1712  pts = IMP.pmi1.tools.ParticleToSampleList()
1713  for n, fb in enumerate(self.floppy_bodies):
1714  pts.add_particle(fb, "Floppy_Bodies", 1.0, "Floppy_Body_" + str(n))
1715  if len(pts.get_particles_to_sample()) > 0:
1716  mc = IMP.pmi1.samplers.MonteCarlo(self.model, [pts], temperature)
1717  print("optimize_floppy_bodies: optimizing %i floppy bodies" % len(self.floppy_bodies))
1718  mc.optimize(nsteps)
1719  else:
1720  print("optimize_floppy_bodies: no particle to optimize")
1721 
1722  def create_rotational_symmetry(self, maincopy, copies, rotational_axis=IMP.algebra.Vector3D(0, 0, 1.0),
1723  nSymmetry=None, skip_gaussian_in_clones=False):
1724  '''
1725  The copies must not contain rigid bodies.
1726  The symmetry restraints are applied at each leaf.
1727  '''
1728 
1729  from math import pi
1730  self.representation_is_modified = True
1731  ncopies = len(copies) + 1
1732  main_hiers = IMP.atom.get_leaves(self.hier_dict[maincopy])
1733 
1734  for k in range(len(copies)):
1735  if (nSymmetry is None):
1736  rotation_angle = 2.0 * pi / float(ncopies) * float(k + 1)
1737  else:
1738  if ( k % 2 == 0 ):
1739  rotation_angle = 2.0 * pi / float(nSymmetry) * float((k + 2) / 2)
1740  else:
1741  rotation_angle = -2.0 * pi / float(nSymmetry) * float((k + 1) / 2)
1742 
1743  rotation3D = IMP.algebra.get_rotation_about_axis(rotational_axis, rotation_angle)
1744 
1745  sm = IMP.core.TransformationSymmetry(rotation3D)
1746  clone_hiers = IMP.atom.get_leaves(self.hier_dict[copies[k]])
1747 
1748  lc = IMP.container.ListSingletonContainer(self.model)
1749  for n, p in enumerate(main_hiers):
1750  if (skip_gaussian_in_clones):
1752  continue
1753  pc = clone_hiers[n]
1754  #print("setting " + p.get_name() + " as reference for " + pc.get_name())
1755  IMP.core.Reference.setup_particle(pc.get_particle(), p.get_particle())
1756  lc.add(pc.get_particle().get_index())
1757 
1758  c = IMP.container.SingletonsConstraint(sm, None, lc)
1759  self.model.add_score_state(c)
1760  print("Completed setting " + str(maincopy) + " as a reference for "
1761  + str(copies[k]) + " by rotating it by "
1762  + str(rotation_angle / 2.0 / pi * 360)
1763  + " degrees around the " + str(rotational_axis) + " axis.")
1764  self.model.update()
1765 
1766  def create_rigid_body_symmetry(self, particles_reference, particles_copy,label="None",
1767  initial_transformation=IMP.algebra.get_identity_transformation_3d()):
1768  from math import pi
1769  self.representation_is_modified = True
1770 
1771  mainparticles = particles_reference
1772 
1773  t=initial_transformation
1774  p=IMP.Particle(self.model)
1775  p.set_name("RigidBody_Symmetry")
1777 
1779 
1780 
1781  copyparticles = particles_copy
1782 
1783  mainpurged = []
1784  copypurged = []
1785  for n, p in enumerate(mainparticles):
1786  print(p.get_name())
1787  pc = copyparticles[n]
1788 
1789  mainpurged.append(p)
1792  else:
1793  IMP.pmi1.Symmetric(p).set_symmetric(0)
1794 
1795  copypurged.append(pc)
1798  else:
1799  IMP.pmi1.Symmetric(pc).set_symmetric(1)
1800 
1801  lc = IMP.container.ListSingletonContainer(self.model)
1802  for n, p in enumerate(mainpurged):
1803 
1804  pc = copypurged[n]
1805  print("setting " + p.get_name() + " as reference for " + pc.get_name())
1806 
1808  lc.add(pc.get_index())
1809 
1810  c = IMP.container.SingletonsConstraint(sm, None, lc)
1811  self.model.add_score_state(c)
1812 
1813  self.model.update()
1814  self.rigid_bodies.append(rb)
1815  self.rigid_body_symmetries.append(rb)
1816  rb.set_name(label+".rigid_body_symmetry."+str(len(self.rigid_body_symmetries)))
1817 
1818 
1819  def create_amyloid_fibril_symmetry(self, maincopy, axial_copies,
1820  longitudinal_copies, axis=(0, 0, 1), translation_value=4.8):
1821 
1822  from math import pi
1823  self.representation_is_modified = True
1824 
1825  outhiers = []
1826  protein_h = self.hier_dict[maincopy]
1827  mainparts = IMP.atom.get_leaves(protein_h)
1828 
1829  for ilc in range(-longitudinal_copies, longitudinal_copies + 1):
1830  for iac in range(axial_copies):
1831  copyname = maincopy + "_a" + str(ilc) + "_l" + str(iac)
1832  self.create_component(copyname, 0.0)
1833  for hier in protein_h.get_children():
1834  self.add_component_hierarchy_clone(copyname, hier)
1835  copyhier = self.hier_dict[copyname]
1836  outhiers.append(copyhier)
1837  copyparts = IMP.atom.get_leaves(copyhier)
1839  IMP.algebra.Vector3D(axis),
1840  2 * pi / axial_copies * (float(iac)))
1841  translation_vector = tuple(
1842  [translation_value * float(ilc) * x for x in axis])
1843  print(translation_vector)
1844  translation = IMP.algebra.Vector3D(translation_vector)
1846  IMP.algebra.Transformation3D(rotation3D, translation))
1847  lc = IMP.container.ListSingletonContainer(self.model)
1848  for n, p in enumerate(mainparts):
1849  pc = copyparts[n]
1855  lc.add(pc.get_index())
1856  c = IMP.container.SingletonsConstraint(sm, None, lc)
1857  self.model.add_score_state(c)
1858  self.model.update()
1859  return outhiers
1860 
1861  def link_components_to_rmf(self, rmfname, frameindex):
1862  '''
1863  Load coordinates in the current representation.
1864  This should be done only if the hierarchy self.prot is identical
1865  to the one as stored in the rmf i.e. all components were added.
1866  '''
1867  import IMP.rmf
1868  import RMF
1869 
1870  rh = RMF.open_rmf_file_read_only(rmfname)
1871  IMP.rmf.link_hierarchies(rh, [self.prot])
1872  IMP.rmf.load_frame(rh, RMF.FrameID(frameindex))
1873  del rh
1874 
1875  def create_components_from_rmf(self, rmfname, frameindex):
1876  '''
1877  still not working.
1878  create the representation (i.e. hierarchies) from the rmf file.
1879  it will be stored in self.prot, which will be overwritten.
1880  load the coordinates from the rmf file at frameindex.
1881  '''
1882  rh = RMF.open_rmf_file_read_only(rmfname)
1883  self.prot = IMP.rmf.create_hierarchies(rh, self.model)[0]
1885  IMP.rmf.link_hierarchies(rh, [self.prot])
1886  IMP.rmf.load_frame(rh, RMF.FrameID(frameindex))
1887  del rh
1888  for p in self.prot.get_children():
1889  self.create_component(p.get_name())
1890  self.hier_dict[p.get_name()] = p
1891  '''
1892  still missing: save rigid bodies contained in the rmf in self.rigid_bodies
1893  save floppy bodies in self.floppy_bodies
1894  get the connectivity restraints
1895  '''
1896 
1897  def set_rigid_body_from_hierarchies(self, hiers, particles=None):
1898  '''
1899  Construct a rigid body from hierarchies (and optionally particles).
1900 
1901  @param hiers list of hierarchies to make rigid
1902  @param particles list of particles to add to the rigid body
1903  '''
1904 
1905  if particles is None:
1906  rigid_parts = set()
1907  else:
1908  rigid_parts = set(particles)
1909 
1910  name = ""
1911  print("set_rigid_body_from_hierarchies> setting up a new rigid body")
1912  for hier in hiers:
1913  ps = IMP.atom.get_leaves(hier)
1914 
1915  for p in ps:
1917  rb = IMP.core.RigidMember(p).get_rigid_body()
1918  print("set_rigid_body_from_hierarchies> WARNING particle %s already belongs to rigid body %s" % (p.get_name(), rb.get_name()))
1919  else:
1920  rigid_parts.add(p)
1921  name += hier.get_name() + "-"
1922  print("set_rigid_body_from_hierarchies> adding %s to the rigid body" % hier.get_name())
1923 
1924  if len(list(rigid_parts)) != 0:
1925  rb = IMP.atom.create_rigid_body(list(rigid_parts))
1926  rb.set_coordinates_are_optimized(True)
1927  rb.set_name(name + "rigid_body")
1928  self.rigid_bodies.append(rb)
1929  return rb
1930 
1931  else:
1932  print("set_rigid_body_from_hierarchies> rigid body could not be setup")
1933 
1934  def set_rigid_bodies(self, subunits):
1935  '''
1936  Construct a rigid body from a list of subunits.
1937 
1938  Each subunit is a tuple that identifies the residue ranges and the
1939  component name (as used in create_component()).
1940 
1941  subunits: [(name_1,(first_residue_1,last_residue_1)),(name_2,(first_residue_2,last_residue_2)),.....]
1942  or
1943  [name_1,name_2,(name_3,(first_residue_3,last_residue_3)),.....]
1944 
1945  example: ["prot1","prot2",("prot3",(1,10))]
1946 
1947  sometimes, we know about structure of an interaction
1948  and here we make such PPIs rigid
1949  '''
1950 
1951  rigid_parts = set()
1952  for s in subunits:
1953  if type(s) == type(tuple()) and len(s) == 2:
1954  sel = IMP.atom.Selection(
1955  self.prot,
1956  molecule=s[0],
1957  residue_indexes=list(range(s[1][0],
1958  s[1][1] + 1)))
1959  if len(sel.get_selected_particles()) == 0:
1960  print("set_rigid_bodies: selected particle does not exist")
1961  for p in sel.get_selected_particles():
1962  # if not p in self.floppy_bodies:
1964  rb = IMP.core.RigidMember(p).get_rigid_body()
1965  print("set_rigid_body_from_hierarchies> WARNING particle %s already belongs to rigid body %s" % (p.get_name(), rb.get_name()))
1966  else:
1967  rigid_parts.add(p)
1968 
1969  elif type(s) == type(str()):
1970  sel = IMP.atom.Selection(self.prot, molecule=s)
1971  if len(sel.get_selected_particles()) == 0:
1972  print("set_rigid_bodies: selected particle does not exist")
1973  for p in sel.get_selected_particles():
1974  # if not p in self.floppy_bodies:
1976  rb = IMP.core.RigidMember(p).get_rigid_body()
1977  print("set_rigid_body_from_hierarchies> WARNING particle %s already belongs to rigid body %s" % (p.get_name(), rb.get_name()))
1978  else:
1979  rigid_parts.add(p)
1980 
1981  rb = IMP.atom.create_rigid_body(list(rigid_parts))
1982  rb.set_coordinates_are_optimized(True)
1983  rb.set_name(''.join(str(subunits)) + "_rigid_body")
1984  self.rigid_bodies.append(rb)
1985  return rb
1986 
1987  def set_super_rigid_body_from_hierarchies(
1988  self,
1989  hiers,
1990  axis=None,
1991  min_size=1):
1992  # axis is the rotation axis for 2D rotation
1993  super_rigid_xyzs = set()
1994  super_rigid_rbs = set()
1995  name = ""
1996  print("set_super_rigid_body_from_hierarchies> setting up a new SUPER rigid body")
1997 
1998  for hier in hiers:
1999  ps = IMP.atom.get_leaves(hier)
2000  for p in ps:
2002  rb = IMP.core.RigidMember(p).get_rigid_body()
2003  super_rigid_rbs.add(rb)
2004  else:
2005  super_rigid_xyzs.add(p)
2006  print("set_rigid_body_from_hierarchies> adding %s to the rigid body" % hier.get_name())
2007  if len(super_rigid_rbs|super_rigid_xyzs) < min_size:
2008  return
2009  if axis is None:
2010  self.super_rigid_bodies.append((super_rigid_xyzs, super_rigid_rbs))
2011  else:
2012  # these will be 2D rotation SRB or a bond rotamer (axis can be a IMP.algebra.Vector3D or particle Pair)
2013  self.super_rigid_bodies.append(
2014  (super_rigid_xyzs, super_rigid_rbs, axis))
2015 
2016  def fix_rigid_bodies(self, rigid_bodies):
2017  self.fixed_rigid_bodies += rigid_bodies
2018 
2019 
2021  self, hiers, min_length=None, max_length=None, axis=None):
2022  '''
2023  Make a chain of super rigid bodies from a list of hierarchies.
2024 
2025  Takes a linear list of hierarchies (they are supposed to be
2026  sequence-contiguous) and produces a chain of super rigid bodies
2027  with given length range, specified by min_length and max_length.
2028  '''
2029  try:
2030  hiers = IMP.pmi1.tools.flatten_list(hiers)
2031  except:
2032  pass
2033  for hs in IMP.pmi1.tools.sublist_iterator(hiers, min_length, max_length):
2034  self.set_super_rigid_body_from_hierarchies(hs, axis, min_length)
2035 
2036  def set_super_rigid_bodies(self, subunits, coords=None):
2037  super_rigid_xyzs = set()
2038  super_rigid_rbs = set()
2039 
2040  for s in subunits:
2041  if type(s) == type(tuple()) and len(s) == 3:
2042  sel = IMP.atom.Selection(
2043  self.prot,
2044  molecule=s[2],
2045  residue_indexes=list(range(s[0],
2046  s[1] + 1)))
2047  if len(sel.get_selected_particles()) == 0:
2048  print("set_rigid_bodies: selected particle does not exist")
2049  for p in sel.get_selected_particles():
2051  rb = IMP.core.RigidMember(p).get_rigid_body()
2052  super_rigid_rbs.add(rb)
2053  else:
2054  super_rigid_xyzs.add(p)
2055  elif type(s) == type(str()):
2056  sel = IMP.atom.Selection(self.prot, molecule=s)
2057  if len(sel.get_selected_particles()) == 0:
2058  print("set_rigid_bodies: selected particle does not exist")
2059  for p in sel.get_selected_particles():
2060  # if not p in self.floppy_bodies:
2062  rb = IMP.core.RigidMember(p).get_rigid_body()
2063  super_rigid_rbs.add(rb)
2064  else:
2065  super_rigid_xyzs.add(p)
2066  self.super_rigid_bodies.append((super_rigid_xyzs, super_rigid_rbs))
2067 
2068  def remove_floppy_bodies_from_component(self, component_name):
2069  '''
2070  Remove leaves of hierarchies from the floppy bodies list based
2071  on the component name
2072  '''
2073  hs=IMP.atom.get_leaves(self.hier_dict[component_name])
2074  ps=[h.get_particle() for h in hs]
2075  for p in self.floppy_bodies:
2076  try:
2077  if p in ps: self.floppy_bodies.remove(p)
2078  if p in hs: self.floppy_bodies.remove(p)
2079  except:
2080  continue
2081 
2082 
2083  def remove_floppy_bodies(self, hierarchies):
2084  '''
2085  Remove leaves of hierarchies from the floppy bodies list.
2086 
2087  Given a list of hierarchies, get the leaves and remove the
2088  corresponding particles from the floppy bodies list. We need this
2089  function because sometimes
2090  we want to constrain the floppy bodies in a rigid body - for instance
2091  when you want to associate a bead with a density particle.
2092  '''
2093  for h in hierarchies:
2094  ps = IMP.atom.get_leaves(h)
2095  for p in ps:
2096  if p in self.floppy_bodies:
2097  print("remove_floppy_bodies: removing %s from floppy body list" % p.get_name())
2098  self.floppy_bodies.remove(p)
2099 
2100 
2101  def set_floppy_bodies(self):
2102  for p in self.floppy_bodies:
2103  name = p.get_name()
2104  p.set_name(name + "_floppy_body")
2106  print("I'm trying to make this particle flexible although it was assigned to a rigid body", p.get_name())
2107  rb = IMP.core.RigidMember(p).get_rigid_body()
2108  try:
2109  rb.set_is_rigid_member(p.get_index(), False)
2110  except:
2111  # some IMP versions still work with that
2112  rb.set_is_rigid_member(p.get_particle_index(), False)
2113  p.set_name(p.get_name() + "_rigid_body_member")
2114 
2115  def set_floppy_bodies_from_hierarchies(self, hiers):
2116  for hier in hiers:
2117  ps = IMP.atom.get_leaves(hier)
2118  for p in ps:
2119  IMP.core.XYZ(p).set_coordinates_are_optimized(True)
2120  self.floppy_bodies.append(p)
2121 
2123  self,
2124  selection_tuples,
2125  resolution=None):
2126  '''
2127  selection tuples must be [(r1,r2,"name1"),(r1,r2,"name2"),....]
2128  @return the particles
2129  '''
2130  particles = []
2131  print(selection_tuples)
2132  for s in selection_tuples:
2133  ps = IMP.pmi1.tools.select_by_tuple(
2134  representation=self, tupleselection=tuple(s),
2135  resolution=None, name_is_ambiguous=False)
2136  particles += ps
2137  return particles
2138 
2139  def get_connected_intra_pairs(self):
2140  return self.connected_intra_pairs
2141 
2142  def set_rigid_bodies_max_trans(self, maxtrans):
2143  self.maxtrans_rb = maxtrans
2144 
2145  def set_rigid_bodies_max_rot(self, maxrot):
2146  self.maxrot_rb = maxrot
2147 
2148  def set_super_rigid_bodies_max_trans(self, maxtrans):
2149  self.maxtrans_srb = maxtrans
2150 
2151  def set_super_rigid_bodies_max_rot(self, maxrot):
2152  self.maxrot_srb = maxrot
2153 
2154  def set_floppy_bodies_max_trans(self, maxtrans):
2155  self.maxtrans_fb = maxtrans
2156 
2157  def set_rigid_bodies_as_fixed(self, rigidbodiesarefixed=True):
2158  '''
2159  Fix rigid bodies in their actual position.
2160  The get_particles_to_sample() function will return
2161  just the floppy bodies (if they are not fixed).
2162  '''
2163  self.rigidbodiesarefixed = rigidbodiesarefixed
2164 
2165  def set_floppy_bodies_as_fixed(self, floppybodiesarefixed=True):
2166  '''
2167  Fix floppy bodies in their actual position.
2168  The get_particles_to_sample() function will return
2169  just the rigid bodies (if they are not fixed).
2170  '''
2171  self.floppybodiesarefixed=floppybodiesarefixed
2172 
2173  def draw_hierarchy_graph(self):
2174  for c in IMP.atom.Hierarchy(self.prot).get_children():
2175  print("Drawing hierarchy graph for " + c.get_name())
2176  IMP.pmi1.output.get_graph_from_hierarchy(c)
2177 
2178  def get_geometries(self):
2179  # create segments at the lowest resolution
2180  seggeos = []
2181  for name in self.hier_geometry_pairs:
2182  for pt in self.hier_geometry_pairs[name]:
2183  p1 = pt[0]
2184  p2 = pt[1]
2185  color = pt[2]
2186  try:
2187  clr = IMP.display.get_rgb_color(color)
2188  except:
2189  clr = IMP.display.get_rgb_color(1.0)
2190  coor1 = IMP.core.XYZ(p1).get_coordinates()
2191  coor2 = IMP.core.XYZ(p2).get_coordinates()
2192  seg = IMP.algebra.Segment3D(coor1, coor2)
2193  seggeos.append(IMP.display.SegmentGeometry(seg, clr))
2194  return seggeos
2195 
2196  def setup_bonds(self):
2197  # create segments at the lowest resolution
2198  seggeos = []
2199  for name in self.hier_geometry_pairs:
2200  for pt in self.hier_geometry_pairs[name]:
2201 
2202  p1 = pt[0]
2203  p2 = pt[1]
2204  if not IMP.atom.Bonded.get_is_setup(p1):
2206  if not IMP.atom.Bonded.get_is_setup(p2):
2208 
2211  IMP.atom.Bonded(p1),
2212  IMP.atom.Bonded(p2),
2213  1)
2214 
2215  def show_component_table(self, name):
2216  if name in self.sequence_dict:
2217  lastresn = len(self.sequence_dict[name])
2218  firstresn = 1
2219  else:
2220  residues = self.hier_db.get_residue_numbers(name)
2221  firstresn = min(residues)
2222  lastresn = max(residues)
2223 
2224  for nres in range(firstresn, lastresn + 1):
2225  try:
2226  resolutions = self.hier_db.get_residue_resolutions(name, nres)
2227  ps = []
2228  for r in resolutions:
2229  ps += self.hier_db.get_particles(name, nres, r)
2230  print("%20s %7s" % (name, nres), " ".join(["%20s %7s" % (str(p.get_name()),
2231  str(IMP.pmi1.Resolution(p).get_resolution())) for p in ps]))
2232  except:
2233  print("%20s %20s" % (name, nres), "**** not represented ****")
2234 
2235  def draw_hierarchy_composition(self):
2236 
2237  ks = list(self.elements.keys())
2238  ks.sort()
2239 
2240  max = 0
2241  for k in self.elements:
2242  for l in self.elements[k]:
2243  if l[1] > max:
2244  max = l[1]
2245 
2246  for k in ks:
2247  self.draw_component_composition(k, max)
2248 
2249  def draw_component_composition(self, name, max=1000, draw_pdb_names=False):
2250  from matplotlib import pyplot
2251  import matplotlib as mpl
2252  k = name
2253  tmplist = sorted(self.elements[k], key=itemgetter(0))
2254  try:
2255  endres = tmplist[-1][1]
2256  except:
2257  print("draw_component_composition: missing information for component %s" % name)
2258  return
2259  fig = pyplot.figure(figsize=(26.0 * float(endres) / max + 2, 2))
2260  ax = fig.add_axes([0.05, 0.475, 0.9, 0.15])
2261 
2262  # Set the colormap and norm to correspond to the data for which
2263  # the colorbar will be used.
2264  cmap = mpl.cm.cool
2265  norm = mpl.colors.Normalize(vmin=5, vmax=10)
2266  bounds = [1]
2267  colors = []
2268 
2269  for n, l in enumerate(tmplist):
2270  firstres = l[0]
2271  lastres = l[1]
2272  if l[3] != "end":
2273  if bounds[-1] != l[0]:
2274  colors.append("white")
2275  bounds.append(l[0])
2276  if l[3] == "pdb":
2277  colors.append("#99CCFF")
2278  if l[3] == "bead":
2279  colors.append("#FFFF99")
2280  if l[3] == "helix":
2281  colors.append("#33CCCC")
2282  if l[3] != "end":
2283  bounds.append(l[1] + 1)
2284  else:
2285  if l[3] == "pdb":
2286  colors.append("#99CCFF")
2287  if l[3] == "bead":
2288  colors.append("#FFFF99")
2289  if l[3] == "helix":
2290  colors.append("#33CCCC")
2291  if l[3] != "end":
2292  bounds.append(l[1] + 1)
2293  else:
2294  if bounds[-1] - 1 == l[0]:
2295  bounds.pop()
2296  bounds.append(l[0])
2297  else:
2298  colors.append("white")
2299  bounds.append(l[0])
2300 
2301  bounds.append(bounds[-1])
2302  colors.append("white")
2303  cmap = mpl.colors.ListedColormap(colors)
2304  cmap.set_over('0.25')
2305  cmap.set_under('0.75')
2306 
2307  norm = mpl.colors.BoundaryNorm(bounds, cmap.N)
2308  cb2 = mpl.colorbar.ColorbarBase(ax, cmap=cmap,
2309  norm=norm,
2310  # to use 'extend', you must
2311  # specify two extra boundaries:
2312  boundaries=bounds,
2313  ticks=bounds, # optional
2314  spacing='proportional',
2315  orientation='horizontal')
2316 
2317  extra_artists = []
2318  npdb = 0
2319 
2320  if draw_pdb_names:
2321  for l in tmplist:
2322  if l[3] == "pdb":
2323  npdb += 1
2324  mid = 1.0 / endres * float(l[0])
2325  # t =ax.text(mid, float(npdb-1)/2.0+1.5, l[2], ha="left", va="center", rotation=0,
2326  # size=10)
2327  # t=ax.annotate(l[0],2)
2328  t = ax.annotate(
2329  l[2], xy=(mid, 1), xycoords='axes fraction',
2330  xytext=(mid + 0.025, float(npdb - 1) / 2.0 + 1.5), textcoords='axes fraction',
2331  arrowprops=dict(arrowstyle="->",
2332  connectionstyle="angle,angleA=0,angleB=90,rad=10"),
2333  )
2334  extra_artists.append(t)
2335 
2336  # set the title of the bar
2337  title = ax.text(-0.005, 0.5, k, ha="right", va="center", rotation=90,
2338  size=20)
2339 
2340  extra_artists.append(title)
2341  # changing the xticks labels
2342  labels = len(bounds) * [" "]
2343  ax.set_xticklabels(labels)
2344  mid = 1.0 / endres * float(bounds[0])
2345  t = ax.annotate(bounds[0], xy=(mid, 0), xycoords='axes fraction',
2346  xytext=(mid - 0.01, -0.5), textcoords='axes fraction',)
2347  extra_artists.append(t)
2348  offsets = [0, -0.5, -1.0]
2349  nclashes = 0
2350  for n in range(1, len(bounds)):
2351  if bounds[n] == bounds[n - 1]:
2352  continue
2353  mid = 1.0 / endres * float(bounds[n])
2354  if (float(bounds[n]) - float(bounds[n - 1])) / max <= 0.01:
2355  nclashes += 1
2356  offset = offsets[nclashes % 3]
2357  else:
2358  nclashes = 0
2359  offset = offsets[0]
2360  if offset > -0.75:
2361  t = ax.annotate(
2362  bounds[n], xy=(mid, 0), xycoords='axes fraction',
2363  xytext=(mid, -0.5 + offset), textcoords='axes fraction')
2364  else:
2365  t = ax.annotate(
2366  bounds[n], xy=(mid, 0), xycoords='axes fraction',
2367  xytext=(mid, -0.5 + offset), textcoords='axes fraction', arrowprops=dict(arrowstyle="-"))
2368  extra_artists.append(t)
2369 
2370  cb2.add_lines(bounds, ["black"] * len(bounds), [1] * len(bounds))
2371  # cb2.set_label(k)
2372 
2373  pyplot.savefig(
2374  k + "structure.pdf",
2375  dpi=150,
2376  transparent="True",
2377  bbox_extra_artists=(extra_artists),
2378  bbox_inches='tight')
2379  pyplot.show()
2380 
2381  def draw_coordinates_projection(self):
2382  import matplotlib.pyplot as pp
2383  xpairs = []
2384  ypairs = []
2385  for name in self.hier_geometry_pairs:
2386  for pt in self.hier_geometry_pairs[name]:
2387  p1 = pt[0]
2388  p2 = pt[1]
2389  color = pt[2]
2390  coor1 = IMP.core.XYZ(p1).get_coordinates()
2391  coor2 = IMP.core.XYZ(p2).get_coordinates()
2392  x1 = coor1[0]
2393  x2 = coor2[0]
2394  y1 = coor1[1]
2395  y2 = coor2[1]
2396  xpairs.append([x1, x2])
2397  ypairs.append([y1, y2])
2398  xlist = []
2399  ylist = []
2400  for xends, yends in zip(xpairs, ypairs):
2401  xlist.extend(xends)
2402  xlist.append(None)
2403  ylist.extend(yends)
2404  ylist.append(None)
2405  pp.plot(xlist, ylist, 'b-', alpha=0.1)
2406  pp.show()
2407 
2408  def get_prot_name_from_particle(self, particle):
2409  names = self.get_component_names()
2410  particle0 = particle
2411  name = None
2412  while not name in names:
2413  h = IMP.atom.Hierarchy(particle0).get_parent()
2414  name = h.get_name()
2415  particle0 = h.get_particle()
2416  return name
2417 
2418 
2419  def get_particles_to_sample(self):
2420  # get the list of samplable particles with their type
2421  # and the mover displacement. Everything wrapped in a dictionary,
2422  # to be used by samplers modules
2423  ps = {}
2424 
2425  # remove symmetric particles: they are not sampled
2426  rbtmp = []
2427  fbtmp = []
2428  srbtmp = []
2429  if not self.rigidbodiesarefixed:
2430  for rb in self.rigid_bodies:
2432  if IMP.pmi1.Symmetric(rb).get_symmetric() != 1:
2433  rbtmp.append(rb)
2434  else:
2435  if rb not in self.fixed_rigid_bodies:
2436  rbtmp.append(rb)
2437 
2438  if not self.floppybodiesarefixed:
2439  for fb in self.floppy_bodies:
2441  if IMP.pmi1.Symmetric(fb).get_symmetric() != 1:
2442  fbtmp.append(fb)
2443  else:
2444  fbtmp.append(fb)
2445 
2446  for srb in self.super_rigid_bodies:
2447  # going to prune the fixed rigid bodies out
2448  # of the super rigid body list
2449  rigid_bodies = list(srb[1])
2450  filtered_rigid_bodies = []
2451  for rb in rigid_bodies:
2452  if rb not in self.fixed_rigid_bodies:
2453  filtered_rigid_bodies.append(rb)
2454  srbtmp.append((srb[0], filtered_rigid_bodies))
2455 
2456  self.rigid_bodies = rbtmp
2457  self.floppy_bodies = fbtmp
2458  self.super_rigid_bodies = srbtmp
2459 
2460  ps["Rigid_Bodies_SimplifiedModel"] = (
2461  self.rigid_bodies,
2462  self.maxtrans_rb,
2463  self.maxrot_rb)
2464  ps["Floppy_Bodies_SimplifiedModel"] = (
2465  self.floppy_bodies,
2466  self.maxtrans_fb)
2467  ps["SR_Bodies_SimplifiedModel"] = (
2468  self.super_rigid_bodies,
2469  self.maxtrans_srb,
2470  self.maxrot_srb)
2471  return ps
2472 
2473  def set_output_level(self, level):
2474  self.output_level = level
2475 
2476  def _evaluate(self, deriv):
2477  """Evaluate the total score of all added restraints"""
2478  r = IMP.pmi1.tools.get_restraint_set(self.model)
2479  return r.evaluate(deriv)
2480 
2481  def get_output(self):
2482  output = {}
2483  score = 0.0
2484 
2485  output["SimplifiedModel_Total_Score_" +
2486  self.label] = str(self._evaluate(False))
2487  output["SimplifiedModel_Linker_Score_" +
2488  self.label] = str(self.linker_restraints.unprotected_evaluate(None))
2489  for name in self.sortedsegments_cr_dict:
2490  partialscore = self.sortedsegments_cr_dict[name].evaluate(False)
2491  score += partialscore
2492  output[
2493  "SimplifiedModel_Link_SortedSegments_" +
2494  name +
2495  "_" +
2496  self.label] = str(
2497  partialscore)
2498  partialscore = self.unmodeledregions_cr_dict[name].evaluate(False)
2499  score += partialscore
2500  output[
2501  "SimplifiedModel_Link_UnmodeledRegions_" +
2502  name +
2503  "_" +
2504  self.label] = str(
2505  partialscore)
2506  for rb in self.rigid_body_symmetries:
2507  name=rb.get_name()
2508  output[name +"_" +self.label]=str(rb.get_reference_frame().get_transformation_to())
2509  for name in self.linker_restraints_dict:
2510  output[
2511  name +
2512  "_" +
2513  self.label] = str(
2514  self.linker_restraints_dict[
2515  name].unprotected_evaluate(
2516  None))
2517 
2518  if len(self.reference_structures.keys()) != 0:
2519  rmsds = self.get_all_rmsds()
2520  for label in rmsds:
2521  output[
2522  "SimplifiedModel_" +
2523  label +
2524  "_" +
2525  self.label] = rmsds[
2526  label]
2527 
2528  if self.output_level == "high":
2529  for p in IMP.atom.get_leaves(self.prot):
2530  d = IMP.core.XYZR(p)
2531  output["Coordinates_" +
2532  p.get_name() + "_" + self.label] = str(d)
2533 
2534  output["_TotalScore"] = str(score)
2535  return output
2536 
2537  def get_test_output(self):
2538  # this method is called by test functions and return an enriched output
2539  output = self.get_output()
2540  for n, p in enumerate(self.get_particles_to_sample()):
2541  output["Particle_to_sample_" + str(n)] = str(p)
2542  output["Number_of_particles"] = len(IMP.atom.get_leaves(self.prot))
2543  output["Hierarchy_Dictionary"] = list(self.hier_dict.keys())
2544  output["Number_of_floppy_bodies"] = len(self.floppy_bodies)
2545  output["Number_of_rigid_bodies"] = len(self.rigid_bodies)
2546  output["Number_of_super_bodies"] = len(self.super_rigid_bodies)
2547  output["Selection_resolution_1"] = len(
2548  IMP.pmi1.tools.select(self, resolution=1))
2549  output["Selection_resolution_5"] = len(
2550  IMP.pmi1.tools.select(self, resolution=5))
2551  output["Selection_resolution_7"] = len(
2552  IMP.pmi1.tools.select(self, resolution=5))
2553  output["Selection_resolution_10"] = len(
2554  IMP.pmi1.tools.select(self, resolution=10))
2555  output["Selection_resolution_100"] = len(
2556  IMP.pmi1.tools.select(self, resolution=100))
2557  output["Selection_All"] = len(IMP.pmi1.tools.select(self))
2558  output["Selection_resolution=1"] = len(
2559  IMP.pmi1.tools.select(self, resolution=1))
2560  output["Selection_resolution=1,resid=10"] = len(
2561  IMP.pmi1.tools.select(self, resolution=1, residue=10))
2562  for resolution in self.hier_resolution:
2563  output["Hier_resolution_dictionary" +
2564  str(resolution)] = len(self.hier_resolution[resolution])
2565  for name in self.hier_dict:
2566  output[
2567  "Selection_resolution=1,resid=10,name=" +
2568  name] = len(
2570  self,
2571  resolution=1,
2572  name=name,
2573  residue=10))
2574  output[
2575  "Selection_resolution=1,resid=10,name=" +
2576  name +
2577  ",ambiguous"] = len(
2579  self,
2580  resolution=1,
2581  name=name,
2582  name_is_ambiguous=True,
2583  residue=10))
2584  output[
2585  "Selection_resolution=1,resid=10,name=" +
2586  name +
2587  ",ambiguous"] = len(
2589  self,
2590  resolution=1,
2591  name=name,
2592  name_is_ambiguous=True,
2593  residue=10))
2594  output[
2595  "Selection_resolution=1,resrange=(10,20),name=" +
2596  name] = len(
2598  self,
2599  resolution=1,
2600  name=name,
2601  first_residue=10,
2602  last_residue=20))
2603  output[
2604  "Selection_resolution=1,resrange=(10,20),name=" +
2605  name +
2606  ",ambiguous"] = len(
2608  self,
2609  resolution=1,
2610  name=name,
2611  name_is_ambiguous=True,
2612  first_residue=10,
2613  last_residue=20))
2614  output[
2615  "Selection_resolution=10,resrange=(10,20),name=" +
2616  name] = len(
2618  self,
2619  resolution=10,
2620  name=name,
2621  first_residue=10,
2622  last_residue=20))
2623  output[
2624  "Selection_resolution=10,resrange=(10,20),name=" +
2625  name +
2626  ",ambiguous"] = len(
2628  self,
2629  resolution=10,
2630  name=name,
2631  name_is_ambiguous=True,
2632  first_residue=10,
2633  last_residue=20))
2634  output[
2635  "Selection_resolution=100,resrange=(10,20),name=" +
2636  name] = len(
2638  self,
2639  resolution=100,
2640  name=name,
2641  first_residue=10,
2642  last_residue=20))
2643  output[
2644  "Selection_resolution=100,resrange=(10,20),name=" +
2645  name +
2646  ",ambiguous"] = len(
2648  self,
2649  resolution=100,
2650  name=name,
2651  name_is_ambiguous=True,
2652  first_residue=10,
2653  last_residue=20))
2654 
2655  return output
static Symmetric setup_particle(Model *m, ParticleIndex pi, Float symmetric)
Definition: /Symmetric.h:44
def sublist_iterator
Yield all sublists of length >= lmin and <= lmax.
Definition: /tools.py:1168
Select non water and non hydrogen atoms.
Definition: pdb.h:331
Store the representations for a system.
Definition: /tools.py:896
Tools for handling Gaussian Mixture Models.
Definition: gmm_tools.py:1
Add mass to a particle.
Definition: Mass.h:23
double get_volume_from_residue_type(ResidueType rt)
Return an estimate for the volume of a given residue.
static Resolution setup_particle(Model *m, ParticleIndex pi, Float resolution)
Definition: /Resolution.h:45
Simple 3D transformation class.
def add_protocol_output
Capture details of the modeling protocol.
static bool get_is_setup(const IMP::ParticleAdaptor &p)
Definition: Residue.h:158
A decorator to associate a particle with a part of a protein/DNA/RNA.
Definition: Fragment.h:20
def add_component_pdb
Add a component that has an associated 3D structure in a PDB file.
static Gaussian setup_particle(Model *m, ParticleIndex pi)
Definition: core/Gaussian.h:65
void show_molecular_hierarchy(Hierarchy h)
Print out the molecular hierarchy.
double get_drms(const Vector3DsOrXYZs0 &m1, const Vector3DsOrXYZs1 &m2)
Apply a SingletonFunction to a SingletonContainer to maintain an invariant.
static bool get_is_setup(Model *m, ParticleIndex pi)
Definition: /Resolution.h:30
atom::Hierarchies create_hierarchies(RMF::FileConstHandle fh, Model *m)
Set the coordinates of a particle to be a transformed version of a reference.
Definition: core/symmetry.h:76
Sampling of the system.
Definition: /samplers.py:1
A member of a rigid body, it has internal (local) coordinates.
Definition: rigid_bodies.h:634
Add resolution to a particle.
Definition: /Resolution.h:24
static bool get_is_setup(const IMP::ParticleAdaptor &p)
Definition: rigid_bodies.h:635
static Provenanced setup_particle(Model *m, ParticleIndex pi, Provenance p)
Definition: provenance.h:657
static Atom setup_particle(Model *m, ParticleIndex pi, Atom other)
Definition: atom/Atom.h:246
static Fragment setup_particle(Model *m, ParticleIndex pi)
Definition: Fragment.h:67
Upper bound harmonic function (non-zero when feature > mean)
def coarse_hierarchy
Generate all needed coarse grained layers.
static XYZR setup_particle(Model *m, ParticleIndex pi)
Definition: XYZR.h:48
A decorator for a particle which has bonds.
def check_root
If the root hierarchy does not exist, construct it.
Select atoms which are selected by both selectors.
Definition: pdb.h:434
static bool get_is_setup(Model *m, ParticleIndex pi)
Definition: /Symmetric.h:29
Rotation3D get_random_rotation_3d(const Rotation3D &center, double distance)
Pick a rotation at random near the provided one.
double get_mass(ResidueType c)
Get the mass from the residue type.
def add_component_hierarchy_clone
Make a copy of a hierarchy and append it to a component.
Color get_rgb_color(double f)
Return the color for f from the RGB color map.
double get_mass_from_number_of_residues(unsigned int num_aa)
Estimate the mass of a protein from the number of amino acids.
def link_components_to_rmf
Load coordinates in the current representation.
Legacy PMI1 module to represent, score, sample and analyze models.
Provenance create_clone(Provenance p)
Clone provenance (including previous provenance)
def create_transformed_component
Create a transformed view of a component.
double get_ball_radius_from_volume_3d(double volume)
Return the radius of a sphere with a given volume.
Definition: Sphere3D.h:35
Set of python classes to create a multi-state, multi-resolution IMP hierarchy.
def get_particles_from_selection_tuples
selection tuples must be [(r1,r2,"name1"),(r1,r2,"name2"),....
def get_residue_gaps_in_hierarchy
Return the residue index gaps and contiguous segments in the hierarchy.
Definition: /tools.py:656
def get_hierarchies_at_given_resolution
Get the hierarchies at the given resolution.
A score on the distance between the surfaces of two spheres.
static Residue setup_particle(Model *m, ParticleIndex pi, ResidueType t, int index, int insertion_code)
Definition: Residue.h:160
static bool get_is_setup(const IMP::ParticleAdaptor &p)
Definition: rigid_bodies.h:749
static bool get_is_setup(const IMP::ParticleAdaptor &p)
Definition: XYZ.h:49
def shuffle_configuration
Shuffle configuration; used to restart the optimization.
def setup_component_sequence_connectivity
Generate restraints between contiguous fragments in the hierarchy.
def get_residue_indexes
Retrieve the residue indexes for the given particle.
Definition: /tools.py:1064
static XYZ setup_particle(Model *m, ParticleIndex pi)
Definition: XYZ.h:51
void read_pdb(TextInput input, int model, Hierarchy h)
def add_component_beads
add beads to the representation
def list_chunks_iterator
Yield successive length-sized chunks from a list.
Definition: /tools.py:1187
Vector3D get_random_vector_in(const Cylinder3D &c)
Generate a random vector in a cylinder with uniform density.
def get_closest_residue_position
this function works with plain hierarchies, as read from the pdb, no multi-scale hierarchies ...
Definition: /tools.py:572
Sample using Monte Carlo.
Definition: /samplers.py:36
A reference frame in 3D.
Bond create_bond(Bonded a, Bonded b, Bond o)
Connect the two wrapped particles by a custom bond.
Object used to hold a set of restraints.
Definition: RestraintSet.h:41
Tools for clustering and cluster analysis.
Definition: pmi1/Analysis.py:1
static Molecule setup_particle(Model *m, ParticleIndex pi)
Definition: Molecule.h:37
Rotation3D get_rotation_about_axis(const Vector3D &axis, double angle)
Generate a Rotation3D object from a rotation around an axis.
Definition: Rotation3D.h:416
Class for storing model, its restraints, constraints, and particles.
Definition: Model.h:86
def get_prot_name_from_particle
Return the component name provided a particle and a list of names.
Definition: /tools.py:1044
def add_component_density
Sets up a Gaussian Mixture Model for this component.
Select all CB ATOM records.
Definition: pdb.h:156
Add uncertainty to a particle.
Definition: /Uncertainty.h:24
A Gaussian distribution in 3D.
Definition: Gaussian3D.h:25
static bool get_is_setup(Model *m, ParticleIndex pi)
Definition: Fragment.h:46
Miscellaneous utilities.
Definition: /tools.py:1
static Hierarchy setup_particle(Model *m, ParticleIndex pi, ParticleIndexesAdaptor children=ParticleIndexesAdaptor())
Create a Hierarchy of level t by adding the needed attributes.
static bool get_is_setup(Model *m, ParticleIndex pi)
Definition: /Uncertainty.h:30
ParticleIndexPairs get_indexes(const ParticlePairsTemp &ps)
Get the indexes from a list of particle pairs.
A dictionary-like wrapper for reading and storing sequence data.
def add_metadata
Associate some metadata with this modeling.
double get_volume_from_mass(double m, ProteinDensityReference ref=ALBER)
Estimate the volume of a protein from its mass.
The standard decorator for manipulating molecular structures.
Ints get_index(const ParticlesTemp &particles, const Subset &subset, const Subsets &excluded)
Store a list of ParticleIndexes.
def deprecated_method
Python decorator to mark a method as deprecated.
Definition: __init__.py:9538
A decorator for a particle representing an atom.
Definition: atom/Atom.h:238
void transform(XYZ a, const algebra::Transformation3D &tr)
Apply a transformation to the particle.
static Mass setup_particle(Model *m, ParticleIndex pi, Float mass)
Definition: Mass.h:48
static Bonded setup_particle(Model *m, ParticleIndex pi)
The type for a residue.
double get_rmsd(const Vector3DsOrXYZs0 &m1, const Vector3DsOrXYZs1 &m2)
def get_file_dataset
Get the dataset associated with a filename, or None.
Track creation of a system fragment from a PDB file.
Definition: provenance.h:86
Basic utilities for handling cryo-electron microscopy 3D density maps.
void load_frame(RMF::FileConstHandle file, RMF::FrameID frame)
Load the given RMF frame into the state of the linked objects.
def set_coordinates_from_rmf
Read and replace coordinates from an RMF file.
A decorator for a particle with x,y,z coordinates.
Definition: XYZ.h:30
def set_file_dataset
Associate a dataset with a filename.
def add_all_atom_densities
Decorates all specified particles as Gaussians directly.
A base class for Keys.
Definition: Key.h:45
def add_component_sequence
Add the primary sequence for a single component.
static Colored setup_particle(Model *m, ParticleIndex pi, Color color)
Definition: Colored.h:62
def write_gmm_to_map
write density map from GMM.
Definition: gmm_tools.py:118
def set_rigid_body_from_hierarchies
Construct a rigid body from hierarchies (and optionally particles).
static Uncertainty setup_particle(Model *m, ParticleIndex pi, Float uncertainty)
Definition: /Uncertainty.h:45
3D rotation class.
Definition: Rotation3D.h:52
Transformation3D get_identity_transformation_3d()
Return a transformation that does not do anything.
Find all nearby pairs by testing all pairs.
Score a pair of particles based on the distance between their centers.
Simple implementation of segments in 3D.
Definition: Segment3D.h:25
static bool get_is_setup(const IMP::ParticleAdaptor &p)
Definition: rigid_bodies.h:169
def remove_floppy_bodies
Remove leaves of hierarchies from the floppy bodies list.
def select
this function uses representation=SimplifiedModel it returns the corresponding selected particles rep...
Definition: /tools.py:727
A decorator for a residue.
Definition: Residue.h:137
Basic functionality that is expected to be used by a wide variety of IMP users.
General purpose algebraic and geometric methods that are expected to be used by a wide variety of IMP...
def get_restraint_set
Get a RestraintSet containing all PMI restraints added to the model.
Definition: /tools.py:72
def set_chain_of_super_rigid_bodies
Make a chain of super rigid bodies from a list of hierarchies.
static bool get_is_setup(const IMP::ParticleAdaptor &p)
def set_rigid_bodies
Construct a rigid body from a list of subunits.
def add_component_necklace
Generates a string of beads with given length.
The general base class for IMP exceptions.
Definition: exception.h:48
def create_rotational_symmetry
The copies must not contain rigid bodies.
Hierarchy create_simplified_along_backbone(Chain input, const IntRanges &residue_segments, bool keep_detailed=false)
VectorD< 3 > Vector3D
Definition: VectorD.h:425
IMP::core::RigidBody create_rigid_body(Hierarchy h)
static Reference setup_particle(Model *m, ParticleIndex pi, ParticleIndexAdaptor reference)
Definition: core/symmetry.h:32
Rotation3D get_identity_rotation_3d()
Return a rotation that does not do anything.
Definition: Rotation3D.h:351
void link_hierarchies(RMF::FileConstHandle fh, const atom::Hierarchies &hs)
Class to handle individual particles of a Model object.
Definition: Particle.h:43
Bond get_bond(Bonded a, Bonded b)
Get the bond between two particles.
std::string get_data_path(std::string file_name)
Return the full path to one of this module's data files.
Select atoms which are selected by either or both selectors.
Definition: pdb.h:457
Add symmetric attribute to a particle.
Definition: /Symmetric.h:23
Store info for a chain of a protein.
Definition: Chain.h:61
Applies a PairScore to a Pair.
Definition: PairRestraint.h:31
Classes for writing output files and processing them.
Definition: /output.py:1
Select all CA ATOM records.
Definition: pdb.h:142
double get_resolution(Model *m, ParticleIndex pi)
Estimate the resolution of the hierarchy as used by Representation.
A decorator for a rigid body.
Definition: rigid_bodies.h:82
static bool get_is_setup(const IMP::ParticleAdaptor &p)
Definition: core/Gaussian.h:64
def add_restraint_to_model
Add a PMI restraint to the model.
Definition: /tools.py:53
def write_gmm_to_text
write a list of gaussians to text.
Definition: gmm_tools.py:61
Transformation3D get_transformation_aligning_first_to_second(Vector3Ds a, Vector3Ds b)
Output IMP model data in various file formats.
Functionality for loading, creating, manipulating and scoring atomic structures.
Tag part of the system to track how it was created.
Definition: provenance.h:632
static Chain setup_particle(Model *m, ParticleIndex pi, std::string id)
Definition: Chain.h:83
Hierarchies get_leaves(const Selection &h)
An exception for an invalid value being passed to IMP.
Definition: exception.h:136
Select hierarchy particles identified by the biological name.
Definition: Selection.h:70
static RigidBody setup_particle(Model *m, ParticleIndex pi, ParticleIndexesAdaptor ps)
Definition: rigid_bodies.h:180
Select all ATOM and HETATM records with the given chain ids.
Definition: pdb.h:256
def create_non_modeled_component
Create a component that isn't used in the modeling.
Support for the RMF file format for storing hierarchical molecular data and markup.
static bool get_is_setup(Model *m, ParticleIndex pi)
Definition: provenance.h:109
def decorate_gmm_from_text
read the output from write_gmm_to_text, decorate as Gaussian and Mass
Definition: gmm_tools.py:23
def remove_floppy_bodies_from_component
Remove leaves of hierarchies from the floppy bodies list based on the component name.
static bool get_is_setup(const IMP::ParticleAdaptor &p)
Definition: rigid_bodies.h:770
Transformation3D get_random_local_transformation(Vector3D origin, double max_translation=5., double max_angle_in_rad=0.26)
Get a local transformation.
def set_rigid_bodies_as_fixed
Fix rigid bodies in their actual position.
Inferential scoring building on methods developed as part of the Inferential Structure Determination ...
def set_floppy_bodies_as_fixed
Fix floppy bodies in their actual position.
Harmonic function (symmetric about the mean)
Definition: core/Harmonic.h:27
A decorator for a particle with x,y,z coordinates and a radius.
Definition: XYZR.h:27