IMP logo
IMP Reference Guide  2.17.0
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 whe 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,
318 
319  # find start and end indexes
320 
321  start = IMP.atom.Residue(
322  t.get_children()[0].get_children()[0]).get_index()
323  end = IMP.atom.Residue(
324  t.get_children()[0].get_children()[-1]).get_index()
325 
326  # check if resrange was defined, otherwise
327  # use the sequence, or the pdb resrange
328 
329  if resrange is None:
330  if name in self.sequence_dict:
331  resrange = (1, len(self.sequence_dict[name]))
332  else:
333  resrange = (start + offset, end + offset)
334  else:
335  if resrange[1] in (-1, 'END'):
336  resrange = (resrange[0],end)
337  start = resrange[0] - offset
338  end = resrange[1] - offset
339 
341  t,
342  resrange[0],
343  resrange[1])
344 
346  t,
347  resrange[0],
348  terminus="N")
350  t,
351  resrange[1],
352  terminus="C")
353  # Done with the PDB
354  IMP.atom.destroy(t)
355 
356  # construct pdb fragments and intervening beads
357  for n, g in enumerate(gaps):
358  first = g[0]
359  last = g[1]
360  if g[2] == "cont":
361  print("autobuild_model: constructing fragment %s from pdb" % (str((first, last))))
362  outhiers += self.add_component_pdb(name, pdbname,
363  chain, resolutions=resolutions,
364  color=color, cacenters=True,
365  resrange=(first, last),
366  offset=offset, isnucleicacid=isnucleicacid)
367  elif g[2] == "gap" and n > 0:
368  print("autobuild_model: constructing fragment %s as a bead" % (str((first, last))))
369  parts = self.hier_db.get_particles_at_closest_resolution(name,
370  first + offset - 1,
371  1)
372  xyz = IMP.core.XYZ(parts[0]).get_coordinates()
373  outhiers += self.add_component_necklace(name,
374  first+offset, last+offset, missingbeadsize, incoord=xyz)
375 
376  elif g[2] == "gap" and n == 0:
377  # add pre-beads
378  print("autobuild_model: constructing fragment %s as a bead" % (str((first, last))))
379  outhiers += self.add_component_necklace(name,
380  first+offset, last+offset, missingbeadsize, incoord=xyznter)
381 
382  return outhiers
383 
384  # Deprecation warning
385 
386  @IMP.deprecated_method("2.5", "Use autobuild_model() instead.")
387  def autobuild_pdb_and_intervening_beads(self, *args, **kwargs):
388  r = self.autobuild_model(*args, **kwargs)
389  return r
390 
391  def add_component_pdb(self, name, pdbname, chain, resolutions, color=None,
392  resrange=None, offset=0, cacenters=True, show=False,
393  isnucleicacid=False, readnonwateratoms=False,
394  read_ca_cb_only=False):
395  '''
396  Add a component that has an associated 3D structure in a PDB file.
397 
398  Reads the PDB, and constructs the fragments corresponding to contiguous
399  sequence stretches.
400 
401  @return a list of hierarchies.
402 
403  @param name (string) the name of the component
404  @param pdbname (string) the name of the PDB file
405  @param chain (string or integer) can be either a string (eg, "A")
406  or an integer (eg, 0 or 1) in case you want
407  to get the corresponding chain number in the PDB.
408  @param resolutions (integers) a list of integers that corresponds
409  to the resolutions that have to be generated
410  @param color (float from 0 to 1) the color applied to the
411  hierarchies generated
412  @param resrange (tuple of integers): the residue range to extract
413  from the PDB. It is a tuple (beg,end). If not specified,
414  all residues belonging to the specified chain are read.
415  @param offset (integer) specifies the residue index offset to be
416  applied when reading the PDB (the FASTA sequence is
417  assumed to start from residue 1, so use this parameter
418  if the PDB file does not start at residue 1)
419  @param cacenters (boolean) if True generates resolution=1 beads
420  centered on C-alpha atoms.
421  @param show (boolean) print out the molecular hierarchy at the end.
422  @param isnucleicacid (boolean) use True if you're reading a PDB
423  with nucleic acids.
424  @param readnonwateratoms (boolean) if True fixes some pathological PDB.
425  @param read_ca_cb_only (boolean) if True, only reads CA/CB
426  '''
427 
428  self.representation_is_modified = True
429  if color is None:
430  # if the color is not passed, then get the stored color
431  color = self.color_dict[name]
432  protein_h = self.hier_dict[name]
433  outhiers = []
434 
435  # determine selector
437  if read_ca_cb_only:
438  cacbsel = IMP.atom.OrPDBSelector(
441  sel = IMP.atom.AndPDBSelector(cacbsel, sel)
442  if type(chain) == str:
445  sel)
446  t = IMP.atom.read_pdb(pdbname, self.model, sel)
447 
448  # get the first and last residue
449  start = IMP.atom.Residue(
450  t.get_children()[0].get_children()[0]).get_index()
451  end = IMP.atom.Residue(
452  t.get_children()[0].get_children()[-1]).get_index()
453  c = IMP.atom.Chain(IMP.atom.get_by_type(t, IMP.atom.CHAIN_TYPE)[0])
454  else:
455  t = IMP.atom.read_pdb(pdbname, self.model, sel)
456  c = IMP.atom.Chain(
457  IMP.atom.get_by_type(t, IMP.atom.CHAIN_TYPE)[chain])
458 
459  # get the first and last residue
460  start = IMP.atom.Residue(c.get_children()[0]).get_index()
461  end = IMP.atom.Residue(c.get_children()[-1]).get_index()
462  chain = c.get_id()
463 
464  if not resrange is None:
465  if resrange[0] > start:
466  start = resrange[0]
467  if resrange[1] < end:
468  end = resrange[1]
469 
470  if not isnucleicacid:
471  # do what you have to do for proteins
472  sel = IMP.atom.Selection(
473  c,
474  residue_indexes=list(range(
475  start,
476  end + 1)),
477  atom_type=IMP.atom.AT_CA)
478 
479  else:
480  # do what you have to do for nucleic-acids
481  # to work, nucleic acids should not be indicated as HETATM in the pdb
482  sel = IMP.atom.Selection(
483  c,
484  residue_indexes=list(range(
485  start,
486  end + 1)),
487  atom_type=IMP.atom.AT_P)
488 
489  ps = sel.get_selected_particles()
490  if len(ps) == 0:
491  raise ValueError("%s no residue found in pdb %s chain %s that overlaps with the queried stretch %s-%s" \
492  % (name, pdbname, str(chain), str(resrange[0]), str(resrange[1])))
493  c0 = IMP.atom.Chain.setup_particle(IMP.Particle(self.model), "X")
494 
495  for p in ps:
496  par = IMP.atom.Atom(p).get_parent()
497  ri = IMP.atom.Residue(par).get_index()
498  # Move residue from original PDB hierarchy to new chain c0
499  IMP.atom.Residue(par).set_index(ri + offset)
500  if par.get_parent() != c0:
501  par.get_parent().remove_child(par)
502  c0.add_child(par)
503  start = start + offset
504  end = end + offset
505 
506  self.elements[name].append(
507  (start, end, pdbname.split("/")[-1] + ":" + chain, "pdb"))
508 
509  hiers = self.coarse_hierarchy(name, start, end,
510  resolutions, isnucleicacid, c0, protein_h, "pdb", color)
511  self._copy_pdb_provenance(t, hiers[0], offset)
512  outhiers += hiers
513  for p, state in self._protocol_output:
514  p.add_pdb_element(state, name, start, end, offset, pdbname, chain,
515  hiers[0])
516 
517  if show:
519 
520  # We cannot simply destroy(c0) since it might not be a well-behaved
521  # hierarchy; in some cases it could contain a given residue more than
522  # once (this is surely a bug but we need to keep this behavior for
523  # backwards compatibility).
524  residues = {}
525  for p in ps:
526  par = IMP.atom.Atom(p).get_parent()
527  residues[par] = None
528  for r in residues.keys():
529  IMP.atom.destroy(r)
530  self.model.remove_particle(c0)
531 
532  IMP.atom.destroy(t)
533 
534  return outhiers
535 
536  def _copy_pdb_provenance(self, pdb, h, offset):
537  """Copy the provenance information from the PDB to our hierarchy"""
538  c = IMP.atom.Chain(IMP.atom.get_by_type(pdb, IMP.atom.CHAIN_TYPE)[0])
539  prov = IMP.core.Provenanced(c).get_provenance()
540  newprov = IMP.core.create_clone(prov)
543  IMP.core.StructureProvenance(newprov).set_residue_offset(offset)
544 
545  def add_component_ideal_helix(
546  self,
547  name,
548  resolutions,
549  resrange,
550  color=None,
551  show=False):
552 
553  self.representation_is_modified = True
554  from math import pi, cos, sin
555 
556  protein_h = self.hier_dict[name]
557  outhiers = []
558  if color is None:
559  color = self.color_dict[name]
560 
561  start = resrange[0]
562  end = resrange[1]
563  self.elements[name].append((start, end, " ", "helix"))
564  c0 = IMP.atom.Chain.setup_particle(IMP.Particle(self.model), "X")
565  for n, res in enumerate(range(start, end + 1)):
566  if name in self.sequence_dict:
567  try:
568  rtstr = self.onetothree[
569  self.sequence_dict[name][res-1]]
570  except:
571  rtstr = "UNK"
572  rt = IMP.atom.ResidueType(rtstr)
573  else:
574  rt = IMP.atom.ResidueType("ALA")
575 
576  # get the residue volume
577  try:
579  # mass=IMP.atom.get_mass_from_residue_type(rt)
580  except IMP.ValueException:
582  IMP.atom.ResidueType("ALA"))
583  # mass=IMP.atom.get_mass_from_residue_type(IMP.atom.ResidueType("ALA"))
585 
586  r = IMP.atom.Residue.setup_particle(IMP.Particle(self.model), rt, res)
587  p = IMP.Particle(self.model)
589  x = 2.3 * cos(n * 2 * pi / 3.6)
590  y = 2.3 * sin(n * 2 * pi / 3.6)
591  z = 6.2 / 3.6 / 2 * n * 2 * pi / 3.6
592  d.set_coordinates(IMP.algebra.Vector3D(x, y, z))
593  d.set_radius(radius)
594  # print d
595  a = IMP.atom.Atom.setup_particle(p, IMP.atom.AT_CA)
596  r.add_child(a)
597  c0.add_child(r)
598 
599  outhiers += self.coarse_hierarchy(name, start, end,
600  resolutions, False, c0, protein_h, "helix", color)
601 
602  if show:
604 
605  IMP.atom.destroy(c0)
606  return outhiers
607 
608  def add_component_beads(self, name, ds, colors=None, incoord=None):
609  """ add beads to the representation
610  @param name the component name
611  @param ds a list of tuples corresponding to the residue ranges
612  of the beads
613  @param colors a list of colors associated to the beads
614  @param incoord the coordinate tuple corresponding to the position
615  of the beads
616  """
617 
618  from math import pi
619  self.representation_is_modified = True
620 
621  protein_h = self.hier_dict[name]
622  outhiers = []
623  if colors is None:
624  colors = [self.color_dict[name]]
625 
626 
627  for n, dss in enumerate(ds):
628  ds_frag = (dss[0], dss[1])
629  self.elements[name].append((dss[0], dss[1], " ", "bead"))
630  prt = IMP.Particle(self.model)
631  if ds_frag[0] == ds_frag[1]:
632  # if the bead represent a single residue
633  if name in self.sequence_dict:
634  try:
635  rtstr = self.onetothree[
636  self.sequence_dict[name][ds_frag[0]-1]]
637  except:
638  rtstr = "UNK"
639  rt = IMP.atom.ResidueType(rtstr)
640  else:
641  rt = IMP.atom.ResidueType("ALA")
642  h = IMP.atom.Residue.setup_particle(prt, rt, ds_frag[0])
643  h.set_name(name + '_%i_bead' % (ds_frag[0]))
644  prt.set_name(name + '_%i_bead' % (ds_frag[0]))
645  resolution = 1
646  else:
648  h.set_name(name + '_%i-%i_bead' % (ds_frag[0], ds_frag[1]))
649  prt.set_name(name + '_%i-%i_bead' % (ds_frag[0], ds_frag[1]))
650  h.set_residue_indexes(list(range(ds_frag[0], ds_frag[1] + 1)))
651  resolution = len(h.get_residue_indexes())
652  if "Beads" not in self.hier_representation[name]:
654  root.set_name("Beads")
655  self.hier_representation[name]["Beads"] = root
656  protein_h.add_child(root)
657  self.hier_representation[name]["Beads"].add_child(h)
658 
659  for kk in range(ds_frag[0], ds_frag[1] + 1):
660  self.hier_db.add_particles(name, kk, resolution, [h])
661 
662  try:
663  clr = IMP.display.get_rgb_color(colors[n])
664  except:
665  clr = IMP.display.get_rgb_color(colors[0])
666 
668 
669  # decorate particles according to their resolution
670  IMP.pmi1.Resolution.setup_particle(prt, resolution)
671 
673  ptem = IMP.core.XYZR(prt)
675  if resolution == 1:
676  try:
678  except IMP.ValueException:
680  IMP.atom.ResidueType("ALA"))
683  ptem.set_radius(radius)
684  else:
685  volume = IMP.atom.get_volume_from_mass(mass)
686  radius = 0.8 * (3.0 / 4.0 / pi * volume) ** (1.0 / 3.0)
688  ptem.set_radius(radius)
689  try:
690  if not tuple(incoord) is None:
691  ptem.set_coordinates(incoord)
692  except TypeError:
693  pass
696  self.floppy_bodies.append(prt)
697  IMP.core.XYZ(prt).set_coordinates_are_optimized(True)
698  outhiers += [h]
699 
700  for p, state in self._protocol_output:
701  p.add_bead_element(state, name, ds[0][0], ds[-1][1], len(ds),
702  outhiers[0])
703 
704  return outhiers
705 
706  def add_component_necklace(self, name, begin, end, length, color=None,
707  incoord=None):
708  '''
709  Generates a string of beads with given length.
710  '''
711  self.representation_is_modified = True
712  outhiers = []
713  if color is None:
714  colors=None
715  else:
716  colors=[color]
717  for chunk in IMP.pmi1.tools.list_chunks_iterator(range(begin, end + 1), length):
718  outhiers += self.add_component_beads(name,
719  [(chunk[0], chunk[-1])], colors=colors,incoord=incoord)
720  return outhiers
721 
723  self, name, hierarchies=None, selection_tuples=None,
724  particles=None,
725  resolution=0.0, num_components=10,
726  inputfile=None, outputfile=None,
727  outputmap=None,
728  kernel_type=None,
729  covariance_type='full', voxel_size=1.0,
730  out_hier_name='',
731  sampled_points=1000000, num_iter=100,
732  simulation_res=1.0,
733  multiply_by_total_mass=True,
734  transform=None,
735  intermediate_map_fn=None,
736  density_ps_to_copy=None,
737  use_precomputed_gaussians=False):
738  '''
739  Sets up a Gaussian Mixture Model for this component.
740  Can specify input GMM file or it will be computed.
741  @param name component name
742  @param hierarchies set up GMM for some hierarchies
743  @param selection_tuples (list of tuples) example (first_residue,last_residue,component_name)
744  @param particles set up GMM for particles directly
745  @param resolution usual PMI resolution for selecting particles from the hierarchies
746  @param inputfile read the GMM from this file
747  @param outputfile fit and write the GMM to this file (text)
748  @param outputmap after fitting, create GMM density file (mrc)
749  @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
750  @param covariance_type for fitting the GMM. options are 'full', 'diagonal' and 'spherical'
751  @param voxel_size for creating the intermediate density map and output map.
752  lower number increases accuracy but also rasterizing time grows
753  @param out_hier_name name of the output density hierarchy
754  @param sampled_points number of points to sample. more will increase accuracy and fitting time
755  @param num_iter num GMM iterations. more will increase accuracy and fitting time
756  @param multiply_by_total_mass multiply the weights of the GMM by this value (only works on creation!)
757  @param transform for input file only, apply a transformation (eg for multiple copies same GMM)
758  @param intermediate_map_fn for debugging, this will write the intermediate (simulated) map
759  @param density_ps_to_copy in case you already created the appropriate GMM (eg, for beads)
760  @param use_precomputed_gaussians Set this flag and pass fragments - will use roughly spherical Gaussian setup
761  '''
762  import numpy as np
763  import sys
764  import IMP.em
765  import IMP.isd.gmm_tools
766 
767  # prepare output
768  self.representation_is_modified = True
769  out_hier = []
770  protein_h = self.hier_dict[name]
771  if "Densities" not in self.hier_representation[name]:
773  root.set_name("Densities")
774  self.hier_representation[name]["Densities"] = root
775  protein_h.add_child(root)
776 
777  # gather passed particles
778  fragment_particles = []
779  if not particles is None:
780  fragment_particles += particles
781  if not hierarchies is None:
782  fragment_particles += IMP.pmi1.tools.select(
783  self, resolution=resolution,
784  hierarchies=hierarchies)
785  if not selection_tuples is None:
786  for st in selection_tuples:
787  fragment_particles += IMP.pmi1.tools.select_by_tuple(
788  self, tupleselection=st,
789  resolution=resolution,
790  name_is_ambiguous=False)
791 
792  # compute or read gaussians
793  density_particles = []
794  if inputfile:
796  inputfile, density_particles,
797  self.model, transform)
798  elif density_ps_to_copy:
799  for ip in density_ps_to_copy:
800  p = IMP.Particle(self.model)
801  shape = IMP.core.Gaussian(ip).get_gaussian()
802  mass = IMP.atom.Mass(ip).get_mass()
805  density_particles.append(p)
806  elif use_precomputed_gaussians:
807  if len(fragment_particles) == 0:
808  print("add_component_density: no particle was selected")
809  return out_hier
810  for p in fragment_particles:
811  if not (IMP.atom.Fragment.get_is_setup(self.model,p.get_particle_index()) and
812  IMP.core.XYZ.get_is_setup(self.model,p.get_particle_index())):
813  raise Exception("The particles you selected must be Fragments and XYZs")
814  nres=len(IMP.atom.Fragment(self.model,p.get_particle_index()).get_residue_indexes())
815  pos=IMP.core.XYZ(self.model,p.get_particle_index()).get_coordinates()
816  density_particles=[]
817  try:
818  IMP.isd.get_data_path("beads/bead_%i.txt"%nres)
819  except:
820  raise Exception("We haven't created a bead file for",nres,"residues")
821  transform = IMP.algebra.Transformation3D(pos)
823  IMP.isd.get_data_path("beads/bead_%i.txt"%nres), density_particles,
824  self.model, transform)
825  else:
826  #compute the gaussians here
827  if len(fragment_particles) == 0:
828  print("add_component_density: no particle was selected")
829  return out_hier
830 
831  density_particles = IMP.isd.gmm_tools.sample_and_fit_to_particles(
832  self.model,
833  fragment_particles,
834  num_components,
835  sampled_points,
836  simulation_res,
837  voxel_size,
838  num_iter,
839  covariance_type,
840  multiply_by_total_mass,
841  outputmap,
842  outputfile)
843 
844  # prepare output hierarchy
846  s0.set_name(out_hier_name)
847  self.hier_representation[name]["Densities"].add_child(s0)
848  out_hier.append(s0)
849  for nps, p in enumerate(density_particles):
850  s0.add_child(p)
851  p.set_name(s0.get_name() + '_gaussian_%i' % nps)
852  return out_hier
853 
854  def get_component_density(self, name):
855  return self.hier_representation[name]["Densities"]
856 
857  def add_all_atom_densities(self, name, hierarchies=None,
858  selection_tuples=None,
859  particles=None,
860  resolution=0,
861  output_map=None,
862  voxel_size=1.0):
863  '''Decorates all specified particles as Gaussians directly.
864  @param name component name
865  @param hierarchies set up GMM for some hierarchies
866  @param selection_tuples (list of tuples) example (first_residue,last_residue,component_name)
867  @param particles set up GMM for particles directly
868  @param resolution usual PMI resolution for selecting particles from the hierarchies
869  '''
870 
871  import IMP.em
872  import numpy as np
873  import sys
874  from math import sqrt
875  self.representation_is_modified = True
876 
877  if particles is None:
878  fragment_particles = []
879  else:
880  fragment_particles = particles
881 
882  if not hierarchies is None:
883  fragment_particles += IMP.pmi1.tools.select(
884  self, resolution=resolution,
885  hierarchies=hierarchies)
886 
887  if not selection_tuples is None:
888  for st in selection_tuples:
889  fragment_particles += IMP.pmi1.tools.select_by_tuple(
890  self, tupleselection=st,
891  resolution=resolution,
892  name_is_ambiguous=False)
893 
894  if len(fragment_particles) == 0:
895  print("add all atom densities: no particle was selected")
896  return
897 
898  # create a spherical gaussian for each particle based on atom type
899  print('setting up all atom gaussians num_particles',len(fragment_particles))
900  for n,p in enumerate(fragment_particles):
901  if IMP.core.Gaussian.get_is_setup(p): continue
902  center=IMP.core.XYZ(p).get_coordinates()
903  rad=IMP.core.XYZR(p).get_radius()
904  mass=IMP.atom.Mass(p).get_mass()
908  print('setting up particle',p.get_name(), " as individual gaussian particle")
909 
910  if not output_map is None:
911  print('writing map to', output_map)
913  fragment_particles,
914  output_map,
915  voxel_size)
916 
917  def add_component_hierarchy_clone(self, name, hierarchy):
918  '''
919  Make a copy of a hierarchy and append it to a component.
920  '''
921  outhier = []
922  self.representation_is_modified = True
923  protein_h = self.hier_dict[name]
924  hierclone = IMP.atom.create_clone(hierarchy)
925  hierclone.set_name(hierclone.get_name() + "_clone")
926  protein_h.add_child(hierclone)
927  outhier.append(hierclone)
928 
929  psmain = IMP.atom.get_leaves(hierarchy)
930  psclone = IMP.atom.get_leaves(hierclone)
931 
932  # copying attributes
933  for n, pmain in enumerate(psmain):
934  pclone = psclone[n]
936  resolution = IMP.pmi1.Resolution(pmain).get_resolution()
937  IMP.pmi1.Resolution.setup_particle(pclone, resolution)
938  for kk in IMP.pmi1.tools.get_residue_indexes(pclone):
939  self.hier_db.add_particles(
940  name,
941  kk,
943  [pclone])
944 
946  uncertainty = IMP.pmi1.Uncertainty(pmain).get_uncertainty()
947  IMP.pmi1.Uncertainty.setup_particle(pclone, uncertainty)
948 
950  symmetric = IMP.pmi1.Symmetric(pmain).get_symmetric()
951  IMP.pmi1.Symmetric.setup_particle(pclone, symmetric)
952 
953  return outhier
954 
955  def dump_particle_descriptors(self):
956  import numpy
957  import pickle
958  import IMP.isd
959  import IMP.isd.gmm_tools
960 
961  particles_attributes={}
962  floppy_body_attributes={}
963  gaussians=[]
964  for h in IMP.atom.get_leaves(self.prot):
965  leaf=h
966  name=h.get_name()
967  hroot=self.prot
968  hparent=h.get_parent()
969  while hparent != hroot:
970  hparent=h.get_parent()
971  name+="|"+hparent.get_name()
972  h=hparent
973  particles_attributes[name]={"COORDINATES":numpy.array(IMP.core.XYZR(leaf.get_particle()).get_coordinates()),
974  "RADIUS":IMP.core.XYZR(leaf.get_particle()).get_radius(),
975  "MASS":IMP.atom.Mass(leaf.get_particle()).get_mass()}
976  if IMP.core.Gaussian.get_is_setup(leaf.get_particle()):
977  gaussians.append(IMP.core.Gaussian(leaf.get_particle()))
978 
979  rigid_body_attributes={}
980  for rb in self.rigid_bodies:
981  name=rb.get_name()
982  rf=rb.get_reference_frame()
983  t=rf.get_transformation_to()
984  trans=t.get_translation()
985  rot=t.get_rotation()
986  rigid_body_attributes[name]={"TRANSLATION":numpy.array(trans),
987  "ROTATION":numpy.array(rot.get_quaternion()),
988  "COORDINATES_NONRIGID_MEMBER":{},
989  "COORDINATES_RIGID_MEMBER":{}}
990  for mi in rb.get_member_indexes():
991  rm=self.model.get_particle(mi)
993  name_part=rm.get_name()
994  xyz=[self.model.get_attribute(fk, rm) for fk in [IMP.FloatKey(4), IMP.FloatKey(5), IMP.FloatKey(6)]]
995  rigid_body_attributes[name]["COORDINATES_NONRIGID_MEMBER"][name_part]=numpy.array(xyz)
996  else:
997  name_part=rm.get_name()
998  xyz=IMP.core.XYZ(rm).get_coordinates()
999  rigid_body_attributes[name]["COORDINATES_RIGID_MEMBER"][name_part]=numpy.array(xyz)
1000 
1001 
1002  IMP.isd.gmm_tools.write_gmm_to_text(gaussians,"model_gmm.txt")
1003  pickle.dump(particles_attributes,
1004  open("particles_attributes.pkl", "wb"))
1005  pickle.dump(rigid_body_attributes,
1006  open("rigid_body_attributes.pkl", "wb"))
1007 
1008 
1009 
1010  def load_particle_descriptors(self):
1011  import numpy
1012  import pickle
1013  import IMP.isd
1014  import IMP.isd.gmm_tools
1015 
1016  particles_attributes = pickle.load(open("particles_attributes.pkl",
1017  "rb"))
1018  rigid_body_attributes = pickle.load(open("rigid_body_attributes.pkl",
1019  "rb"))
1020 
1021  particles=[]
1022  hierarchies=[]
1023  gaussians=[]
1024  for h in IMP.atom.get_leaves(self.prot):
1025  leaf=h
1026  name=h.get_name()
1027  hroot=self.prot
1028  hparent=h.get_parent()
1029  while hparent != hroot:
1030  hparent=h.get_parent()
1031  name+="|"+hparent.get_name()
1032  h=hparent
1033 
1034  xyzr=IMP.core.XYZR(leaf.get_particle())
1035  xyzr.set_coordinates(particles_attributes[name]["COORDINATES"])
1036  #xyzr.set_radius(particles_attributes[name]["RADIUS"])
1037  #IMP.atom.Mass(leaf.get_particle()).set_mass(particles_attributes[name]["MASS"])
1038  if IMP.core.Gaussian.get_is_setup(leaf.get_particle()):
1039  gaussians.append(IMP.core.Gaussian(leaf.get_particle()))
1040 
1041  for rb in self.rigid_bodies:
1042  name=rb.get_name()
1043  trans=rigid_body_attributes[name]["TRANSLATION"]
1044  rot=rigid_body_attributes[name]["ROTATION"]
1047  rb.set_reference_frame(rf)
1048  coor_nrm_ref=rigid_body_attributes[name]["COORDINATES_NONRIGID_MEMBER"]
1049  coor_rm_ref_dict=rigid_body_attributes[name]["COORDINATES_RIGID_MEMBER"]
1050  coor_rm_model=[]
1051  coor_rm_ref=[]
1052  for mi in rb.get_member_indexes():
1053  rm=self.model.get_particle(mi)
1055  name_part=rm.get_name()
1056  xyz=coor_nrm_ref[name_part]
1057  for n,fk in enumerate([IMP.FloatKey(4), IMP.FloatKey(5), IMP.FloatKey(6)]):
1058  self.model.set_attribute(fk, rm,xyz[n])
1059  else:
1060  name_part=rm.get_name()
1061  coor_rm_ref.append(IMP.algebra.Vector3D(coor_rm_ref_dict[name_part]))
1062  coor_rm_model.append(IMP.core.XYZ(rm).get_coordinates())
1063  if len(coor_rm_model)==0: continue
1064  t=IMP.algebra.get_transformation_aligning_first_to_second(coor_rm_model,coor_rm_ref)
1065  IMP.core.transform(rb,t)
1066 
1067  IMP.isd.gmm_tools.decorate_gmm_from_text("model_gmm.txt",gaussians,self.model)
1068 
1069  def _compare_rmf_repr_names(self, rmfname, reprname, component_name):
1070  """Print a warning if particle names in RMF and model don't match"""
1071  def match_any_suffix():
1072  # Handle common mismatches like 743 != Nup85_743_pdb
1073  suffixes = ["pdb", "bead_floppy_body_rigid_body_member_floppy_body",
1074  "bead_floppy_body_rigid_body_member",
1075  "bead_floppy_body"]
1076  for s in suffixes:
1077  if "%s_%s_%s" % (component_name, rmfname, s) == reprname:
1078  return True
1079  if rmfname != reprname and not match_any_suffix():
1080  print("set_coordinates_from_rmf: WARNING rmf particle and "
1081  "representation particle names don't match %s %s"
1082  % (rmfname, reprname))
1083 
1084  def set_coordinates_from_rmf(self, component_name, rmf_file_name,
1085  rmf_frame_number,
1086  rmf_component_name=None,
1087  check_number_particles=True,
1088  representation_name_to_rmf_name_map=None,
1089  state_number=0,
1090  skip_gaussian_in_rmf=False,
1091  skip_gaussian_in_representation=False,
1092  save_file=False,
1093  force_rigid_update=False):
1094  '''Read and replace coordinates from an RMF file.
1095  Replace the coordinates of particles with the same name.
1096  It assumes that the RMF and the representation have the particles
1097  in the same order.
1098  @param component_name Component name
1099  @param rmf_component_name Name of the component in the RMF file
1100  (if not specified, use `component_name`)
1101  @param representation_name_to_rmf_name_map a dictionary that map
1102  the original rmf particle name to the recipient particle component name
1103  @param save_file: save a file with the names of particles of the component
1104  @param force_rigid_update: update the coordinates of rigid bodies
1105  (normally this should be called before rigid bodies are set up)
1106  '''
1107  import IMP.pmi1.analysis
1108 
1109  tempm = IMP.Model()
1110  prots = IMP.pmi1.analysis.get_hiers_from_rmf(
1111  tempm,
1112  rmf_frame_number,
1113  rmf_file_name)
1114 
1115  if not prots:
1116  raise ValueError("cannot read hierarchy from rmf")
1117 
1118  prot=prots[0]
1119  # Make sure coordinates of rigid body members in the RMF are correct
1120  if force_rigid_update:
1121  tempm.update()
1122 
1123  # if len(self.rigid_bodies)!=0:
1124  # print "set_coordinates_from_rmf: cannot proceed if rigid bodies were initialized. Use the function before defining the rigid bodies"
1125  # exit()
1126 
1127  allpsrmf = IMP.atom.get_leaves(prot)
1128  psrmf = []
1129  for p in allpsrmf:
1130  (protname, is_a_bead) = IMP.pmi1.tools.get_prot_name_from_particle(
1131  p, self.hier_dict.keys())
1132  if (protname is None) and (rmf_component_name is not None):
1133  (protname, is_a_bead) = IMP.pmi1.tools.get_prot_name_from_particle(
1134  p, rmf_component_name)
1135  if (skip_gaussian_in_rmf):
1137  continue
1138  if (rmf_component_name is not None) and (protname == rmf_component_name):
1139  psrmf.append(p)
1140  elif (rmf_component_name is None) and (protname == component_name):
1141  psrmf.append(p)
1142 
1143  psrepr = IMP.atom.get_leaves(self.hier_dict[component_name])
1144  if (skip_gaussian_in_representation):
1145  allpsrepr = psrepr
1146  psrepr = []
1147  for p in allpsrepr:
1148  #(protname, is_a_bead) = IMP.pmi1.tools.get_prot_name_from_particle(
1149  # p, self.hier_dict.keys())
1151  continue
1152  psrepr.append(p)
1153 
1154  import itertools
1155  reprnames=[p.get_name() for p in psrepr]
1156  rmfnames=[p.get_name() for p in psrmf]
1157 
1158  if save_file:
1159  fl=open(component_name+".txt","w")
1160  for i in itertools.izip_longest(reprnames,rmfnames): fl.write(str(i[0])+","+str(i[1])+"\n")
1161 
1162 
1163  if check_number_particles and not representation_name_to_rmf_name_map:
1164  if len(psrmf) != len(psrepr):
1165  fl=open(component_name+".txt","w")
1166  for i in itertools.izip_longest(reprnames,rmfnames): fl.write(str(i[0])+","+str(i[1])+"\n")
1167  raise ValueError("%s cannot proceed the rmf and the representation don't have the same number of particles; "
1168  "particles in rmf: %s particles in the representation: %s" % (str(component_name), str(len(psrmf)), str(len(psrepr))))
1169 
1170 
1171  if not representation_name_to_rmf_name_map:
1172  for n, prmf in enumerate(psrmf):
1173 
1174  prmfname = prmf.get_name()
1175  preprname = psrepr[n].get_name()
1176  if force_rigid_update:
1177  if IMP.core.RigidBody.get_is_setup(psrepr[n]) \
1178  and not IMP.core.RigidBodyMember.get_is_setup(psrepr[n]):
1179  continue
1180  else:
1182  raise ValueError("component %s cannot proceed if rigid bodies were initialized. Use the function before defining the rigid bodies" % component_name)
1183 
1184  self._compare_rmf_repr_names(prmfname, preprname,
1185  component_name)
1186  if IMP.core.XYZ.get_is_setup(prmf) and IMP.core.XYZ.get_is_setup(psrepr[n]):
1187  xyz = IMP.core.XYZ(prmf).get_coordinates()
1188  IMP.core.XYZ(psrepr[n]).set_coordinates(xyz)
1190  # Set rigid body so that coordinates are preserved
1191  # on future model updates
1192  rbm = IMP.core.RigidBodyMember(psrepr[n])
1193  rbm.set_internal_coordinates(xyz)
1195  rb = rbm.get_rigid_body()
1197  raise ValueError("Cannot handle nested "
1198  "rigid bodies yet")
1199  rb.set_reference_frame_lazy(tr)
1200  else:
1201  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]))))
1202 
1204  gprmf=IMP.core.Gaussian(prmf)
1205  grepr=IMP.core.Gaussian(psrepr[n])
1206  g=gprmf.get_gaussian()
1207  grepr.set_gaussian(g)
1208 
1209  else:
1210  repr_name_particle_map={}
1211  rmf_name_particle_map={}
1212  for p in psrmf:
1213  rmf_name_particle_map[p.get_name()]=p
1214  #for p in psrepr:
1215  # repr_name_particle_map[p.get_name()]=p
1216 
1217  for prepr in psrepr:
1218  try:
1219  prmf=rmf_name_particle_map[representation_name_to_rmf_name_map[prepr.get_name()]]
1220  except KeyError:
1221  print("set_coordinates_from_rmf: WARNING missing particle name in representation_name_to_rmf_name_map, skipping...")
1222  continue
1223  xyz = IMP.core.XYZ(prmf).get_coordinates()
1224  IMP.core.XYZ(prepr).set_coordinates(xyz)
1225 
1226 
1227 
1228  def check_root(self, name, protein_h, resolution):
1229  '''
1230  If the root hierarchy does not exist, construct it.
1231  '''
1232  if "Res:" + str(int(resolution)) not in self.hier_representation[name]:
1234  root.set_name(name + "_Res:" + str(int(resolution)))
1235  self.hier_representation[name][
1236  "Res:" + str(int(resolution))] = root
1237  protein_h.add_child(root)
1238 
1239  def coarse_hierarchy(self, name, start, end, resolutions, isnucleicacid,
1240  input_hierarchy, protein_h, type, color):
1241  '''
1242  Generate all needed coarse grained layers.
1243 
1244  @param name name of the protein
1245  @param resolutions list of resolutions
1246  @param protein_h root hierarchy
1247  @param input_hierarchy hierarchy to coarse grain
1248  @param type a string, typically "pdb" or "helix"
1249  '''
1250  outhiers = []
1251 
1252  if (1 in resolutions) or (0 in resolutions):
1253  # in that case create residues and append atoms
1254 
1255  if 1 in resolutions:
1256  self.check_root(name, protein_h, 1)
1258  s1.set_name('%s_%i-%i_%s' % (name, start, end, type))
1259  # s1.set_residue_indexes(range(start,end+1))
1260  self.hier_representation[name]["Res:1"].add_child(s1)
1261  outhiers += [s1]
1262  if 0 in resolutions:
1263  self.check_root(name, protein_h, 0)
1265  s0.set_name('%s_%i-%i_%s' % (name, start, end, type))
1266  # s0.set_residue_indexes(range(start,end+1))
1267  self.hier_representation[name]["Res:0"].add_child(s0)
1268  outhiers += [s0]
1269 
1270  if not isnucleicacid:
1271  sel = IMP.atom.Selection(
1272  input_hierarchy,
1273  atom_type=IMP.atom.AT_CA)
1274  else:
1275  sel = IMP.atom.Selection(
1276  input_hierarchy,
1277  atom_type=IMP.atom.AT_P)
1278 
1279  for p in sel.get_selected_particles():
1280  resobject = IMP.atom.Residue(IMP.atom.Atom(p).get_parent())
1281  if 0 in resolutions:
1282  # if you ask for atoms
1283  resclone0 = IMP.atom.create_clone(resobject)
1284  resindex = IMP.atom.Residue(resclone0).get_index()
1285  s0.add_child(resclone0)
1286  self.hier_db.add_particles(
1287  name,
1288  resindex,
1289  0,
1290  resclone0.get_children())
1291 
1292  chil = resclone0.get_children()
1293  for ch in chil:
1295  try:
1296  clr = IMP.display.get_rgb_color(color)
1297  except:
1298  clr = IMP.display.get_rgb_color(1.0)
1300 
1301  if 1 in resolutions:
1302  # else clone the residue
1303  resclone1 = IMP.atom.create_clone_one(resobject)
1304  resindex = IMP.atom.Residue(resclone1).get_index()
1305  s1.add_child(resclone1)
1306  self.hier_db.add_particles(name, resindex, 1, [resclone1])
1307 
1308  rt = IMP.atom.Residue(resclone1).get_residue_type()
1309  xyz = IMP.core.XYZ(p).get_coordinates()
1310  prt = resclone1.get_particle()
1311  prt.set_name('%s_%i_%s' % (name, resindex, type))
1312  IMP.core.XYZ.setup_particle(prt).set_coordinates(xyz)
1313 
1314  try:
1316  # mass=IMP.atom.get_mass_from_residue_type(rt)
1317  except IMP.ValueException:
1319  IMP.atom.ResidueType("ALA"))
1320  # mass=IMP.atom.get_mass_from_residue_type(IMP.atom.ResidueType("ALA"))
1322  IMP.core.XYZR.setup_particle(prt).set_radius(radius)
1327  try:
1328  clr = IMP.display.get_rgb_color(color)
1329  except:
1330  clr = IMP.display.get_rgb_color(1.0)
1332 
1333  for r in resolutions:
1334  if r != 0 and r != 1:
1335  self.check_root(name, protein_h, r)
1337  input_hierarchy,
1338  r)
1339 
1340  chil = s.get_children()
1342  s0.set_name('%s_%i-%i_%s' % (name, start, end, type))
1343  # Move all children from s to s0
1344  for ch in chil:
1345  s.remove_child(ch)
1346  s0.add_child(ch)
1347  self.hier_representation[name][
1348  "Res:" + str(int(r))].add_child(s0)
1349  outhiers += [s0]
1350  IMP.atom.destroy(s)
1351 
1352  for prt in IMP.atom.get_leaves(s0):
1353  ri = IMP.atom.Fragment(prt).get_residue_indexes()
1354  first = ri[0]
1355  last = ri[-1]
1356  if first == last:
1357  prt.set_name('%s_%i_%s' % (name, first, type))
1358  else:
1359  prt.set_name('%s_%i_%i_%s' % (name, first, last, type))
1360  for kk in ri:
1361  self.hier_db.add_particles(name, kk, r, [prt])
1362 
1363  radius = IMP.core.XYZR(prt).get_radius()
1367 
1368  # setting up color for each particle in the
1369  # hierarchy, if colors missing in the colors list set it to
1370  # red
1371  try:
1372  clr = IMP.display.get_rgb_color(color)
1373  except:
1374  colors.append(1.0)
1375  clr = IMP.display.get_rgb_color(colors[pdb_part_count])
1377 
1378  return outhiers
1379 
1381  '''
1382  Get the hierarchies at the given resolution.
1383 
1384  The map between resolution and hierarchies is cached to accelerate
1385  the selection algorithm. The cache is invalidated when the
1386  representation was changed by any add_component_xxx.
1387  '''
1388 
1389  if self.representation_is_modified:
1390  rhbr = self.hier_db.get_all_root_hierarchies_by_resolution(
1391  resolution)
1392  self.hier_resolution[resolution] = rhbr
1393  self.representation_is_modified = False
1394  return rhbr
1395  else:
1396  if resolution in self.hier_resolution:
1397  return self.hier_resolution[resolution]
1398  else:
1399  rhbr = self.hier_db.get_all_root_hierarchies_by_resolution(
1400  resolution)
1401  self.hier_resolution[resolution] = rhbr
1402  return rhbr
1403 
1405  self, max_translation=300., max_rotation=2.0 * pi,
1406  avoidcollision=True, cutoff=10.0, niterations=100,
1407  bounding_box=None,
1408  excluded_rigid_bodies=None,
1409  ignore_initial_coordinates=False,
1410  hierarchies_excluded_from_collision=None):
1411  '''
1412  Shuffle configuration; used to restart the optimization.
1413 
1414  The configuration of the system is initialized by placing each
1415  rigid body and each bead randomly in a box with a side of
1416  `max_translation` angstroms, and far enough from each other to
1417  prevent any steric clashes. The rigid bodies are also randomly rotated.
1418 
1419  @param avoidcollision check if the particle/rigid body was
1420  placed close to another particle; uses the optional
1421  arguments cutoff and niterations
1422  @param bounding_box defined by ((x1,y1,z1),(x2,y2,z2))
1423  '''
1424 
1425  if excluded_rigid_bodies is None:
1426  excluded_rigid_bodies = []
1427  if hierarchies_excluded_from_collision is None:
1428  hierarchies_excluded_from_collision = []
1429 
1430  if len(self.rigid_bodies) == 0:
1431  print("shuffle_configuration: rigid bodies were not intialized")
1432 
1434  gcpf.set_distance(cutoff)
1435  ps = []
1436  hierarchies_excluded_from_collision_indexes = []
1437  for p in IMP.atom.get_leaves(self.prot):
1439  ps.append(p)
1441  # remove the densities particles out of the calculation
1442  hierarchies_excluded_from_collision_indexes.extend(
1443  IMP.get_indexes([p]))
1444  allparticleindexes = IMP.get_indexes(ps)
1445 
1446  if bounding_box is not None:
1447  ((x1, y1, z1), (x2, y2, z2)) = bounding_box
1448  lb = IMP.algebra.Vector3D(x1, y1, z1)
1449  ub = IMP.algebra.Vector3D(x2, y2, z2)
1450  bb = IMP.algebra.BoundingBox3D(lb, ub)
1451 
1452  for h in hierarchies_excluded_from_collision:
1453  hierarchies_excluded_from_collision_indexes.extend(
1455 
1456 
1457  allparticleindexes = list(
1458  set(allparticleindexes) - set(hierarchies_excluded_from_collision_indexes))
1459 
1460  print(hierarchies_excluded_from_collision)
1461  print(len(allparticleindexes),len(hierarchies_excluded_from_collision_indexes))
1462 
1463  print('shuffling', len(self.rigid_bodies) - len(excluded_rigid_bodies), 'rigid bodies')
1464  for rb in self.rigid_bodies:
1465  if rb not in excluded_rigid_bodies:
1466  if avoidcollision:
1467  rbindexes = rb.get_member_particle_indexes()
1468  rbindexes = list(
1469  set(rbindexes) - set(hierarchies_excluded_from_collision_indexes))
1470  otherparticleindexes = list(
1471  set(allparticleindexes) - set(rbindexes))
1472 
1473  if len(otherparticleindexes) is None:
1474  continue
1475 
1476  niter = 0
1477  while niter < niterations:
1478  if (ignore_initial_coordinates):
1479  # Move the particle to the origin
1480  transformation = IMP.algebra.Transformation3D(IMP.algebra.get_identity_rotation_3d(), -IMP.core.XYZ(rb).get_coordinates())
1481  IMP.core.transform(rb, transformation)
1482  rbxyz = IMP.core.XYZ(rb).get_coordinates()
1483 
1484  if bounding_box is not None:
1485  # overrides the perturbation
1486  translation = IMP.algebra.get_random_vector_in(bb)
1488  transformation = IMP.algebra.Transformation3D(rotation, translation-rbxyz)
1489  else:
1491  rbxyz,
1492  max_translation,
1493  max_rotation)
1494 
1495  IMP.core.transform(rb, transformation)
1496 
1497  if avoidcollision:
1498  self.model.update()
1499  npairs = len(
1500  gcpf.get_close_pairs(
1501  self.model,
1502  otherparticleindexes,
1503  rbindexes))
1504  if npairs == 0:
1505  niter = niterations
1506  if (ignore_initial_coordinates):
1507  print (rb.get_name(), IMP.core.XYZ(rb).get_coordinates())
1508  else:
1509  niter += 1
1510  print("shuffle_configuration: rigid body placed close to other %d particles, trying again..." % npairs)
1511  print("shuffle_configuration: rigid body name: " + rb.get_name())
1512  if niter == niterations:
1513  raise ValueError("tried the maximum number of iterations to avoid collisions, increase the distance cutoff")
1514  else:
1515  break
1516 
1517  print('shuffling', len(self.floppy_bodies), 'floppy bodies')
1518  for fb in self.floppy_bodies:
1519  if (avoidcollision):
1522  if rm and nrm:
1523  fbindexes = IMP.get_indexes([fb])
1524  otherparticleindexes = list(
1525  set(allparticleindexes) - set(fbindexes))
1526  if len(otherparticleindexes) is None:
1527  continue
1528  else:
1529  continue
1530  else:
1533  if (rm or nrm):
1534  continue
1535 
1537  d=IMP.core.RigidBodyMember(fb).get_rigid_body()
1539  d=IMP.core.RigidBody(fb)
1540  elif IMP.core.XYZ.get_is_setup(fb):
1541  d=IMP.core.XYZ(fb)
1542 
1543  niter = 0
1544  while niter < niterations:
1545  if (ignore_initial_coordinates):
1546  # Move the particle to the origin
1547  transformation = IMP.algebra.Transformation3D(IMP.algebra.get_identity_rotation_3d(), -IMP.core.XYZ(fb).get_coordinates())
1548  IMP.core.transform(d, transformation)
1549  fbxyz = IMP.core.XYZ(fb).get_coordinates()
1550 
1551  if bounding_box is not None:
1552  # overrides the perturbation
1553  translation = IMP.algebra.get_random_vector_in(bb)
1554  transformation = IMP.algebra.Transformation3D(translation-fbxyz)
1555  else:
1557  fbxyz,
1558  max_translation,
1559  max_rotation)
1560 
1561  IMP.core.transform(d, transformation)
1562 
1563  if (avoidcollision):
1564  self.model.update()
1565  npairs = len(
1566  gcpf.get_close_pairs(
1567  self.model,
1568  otherparticleindexes,
1569  fbindexes))
1570  if npairs == 0:
1571  niter = niterations
1572  if (ignore_initial_coordinates):
1573  print (fb.get_name(), IMP.core.XYZ(fb).get_coordinates())
1574  else:
1575  niter += 1
1576  print("shuffle_configuration: floppy body placed close to other %d particles, trying again..." % npairs)
1577  print("shuffle_configuration: floppy body name: " + fb.get_name())
1578  if niter == niterations:
1579  raise ValueError("tried the maximum number of iterations to avoid collisions, increase the distance cutoff")
1580  else:
1581  break
1582 
1583  def set_current_coordinates_as_reference_for_rmsd(self, label="None"):
1584  # getting only coordinates from pdb
1585  ps = IMP.pmi1.tools.select(self, resolution=1.0)
1586  # storing the reference coordinates and the particles
1587  self.reference_structures[label] = (
1588  [IMP.core.XYZ(p).get_coordinates() for p in ps],
1589  ps)
1590 
1591  def get_all_rmsds(self):
1592  rmsds = {}
1593 
1594  for label in self.reference_structures:
1595 
1596  current_coordinates = [IMP.core.XYZ(p).get_coordinates()
1597  for p in self.reference_structures[label][1]]
1598  reference_coordinates = self.reference_structures[label][0]
1599  if len(reference_coordinates) != len(current_coordinates):
1600  print("calculate_all_rmsds: reference and actual coordinates are not the same")
1601  continue
1603  current_coordinates,
1604  reference_coordinates)
1605  rmsd_global = IMP.algebra.get_rmsd(
1606  reference_coordinates,
1607  current_coordinates)
1608  # warning: temporary we are calculating the drms, and not the rmsd,
1609  # for the relative distance
1610  rmsd_relative = IMP.atom.get_drms(
1611  reference_coordinates,
1612  current_coordinates)
1613  rmsds[label + "_GlobalRMSD"] = rmsd_global
1614  rmsds[label + "_RelativeDRMS"] = rmsd_relative
1615  return rmsds
1616 
1617  def setup_component_geometry(self, name, color=None, resolution=1.0):
1618  if color is None:
1619  color = self.color_dict[name]
1620  # this function stores all particle pairs
1621  # ordered by residue number, to be used
1622  # to construct backbone traces
1623  self.hier_geometry_pairs[name] = []
1624  protein_h = self.hier_dict[name]
1625  pbr = IMP.pmi1.tools.select(self, name=name, resolution=resolution)
1626  pbr = IMP.pmi1.tools.sort_by_residues(pbr)
1627 
1628  for n in range(len(pbr) - 1):
1629  self.hier_geometry_pairs[name].append((pbr[n], pbr[n + 1], color))
1630 
1632  self, name, resolution=10, scale=1.0):
1633  '''
1634  Generate restraints between contiguous fragments in the hierarchy.
1635  The linkers are generated at resolution 10 by default.
1636 
1637  '''
1638 
1639  unmodeledregions_cr = IMP.RestraintSet(self.model, "unmodeledregions")
1640  sortedsegments_cr = IMP.RestraintSet(self.model, "sortedsegments")
1641 
1642  protein_h = self.hier_dict[name]
1643  SortedSegments = []
1644  frs = self.hier_db.get_preroot_fragments_by_resolution(
1645  name,
1646  resolution)
1647 
1648  for fr in frs:
1649  try:
1650  start = fr.get_children()[0]
1651  except:
1652  start = fr
1653 
1654  try:
1655  end = fr.get_children()[-1]
1656  except:
1657  end = fr
1658 
1659  startres = IMP.pmi1.tools.get_residue_indexes(start)[0]
1660  endres = IMP.pmi1.tools.get_residue_indexes(end)[-1]
1661  SortedSegments.append((start, end, startres))
1662  SortedSegments = sorted(SortedSegments, key=itemgetter(2))
1663 
1664  # connect the particles
1665  for x in range(len(SortedSegments) - 1):
1666  last = SortedSegments[x][1]
1667  first = SortedSegments[x + 1][0]
1668 
1669  nreslast = len(IMP.pmi1.tools.get_residue_indexes(last))
1670  lastresn = IMP.pmi1.tools.get_residue_indexes(last)[-1]
1671  nresfirst = len(IMP.pmi1.tools.get_residue_indexes(first))
1672  firstresn = IMP.pmi1.tools.get_residue_indexes(first)[0]
1673 
1674  residuegap = firstresn - lastresn - 1
1675  if self.disorderedlength and (nreslast / 2 + nresfirst / 2 + residuegap) > 20.0:
1676  # calculate the distance between the sphere centers using Kohn
1677  # PNAS 2004
1678  optdist = sqrt(5 / 3) * 1.93 * \
1679  (nreslast / 2 + nresfirst / 2 + residuegap) ** 0.6
1680  # optdist2=sqrt(5/3)*1.93*((nreslast)**0.6+(nresfirst)**0.6)/2
1681  if self.upperharmonic:
1682  hu = IMP.core.HarmonicUpperBound(optdist, self.kappa)
1683  else:
1684  hu = IMP.core.Harmonic(optdist, self.kappa)
1685  dps = IMP.core.DistancePairScore(hu)
1686  else: # default
1687  optdist = (0.0 + (float(residuegap) + 1.0) * 3.6) * scale
1688  if self.upperharmonic: # default
1689  hu = IMP.core.HarmonicUpperBound(optdist, self.kappa)
1690  else:
1691  hu = IMP.core.Harmonic(optdist, self.kappa)
1693 
1694  pt0 = last.get_particle()
1695  pt1 = first.get_particle()
1696  r = IMP.core.PairRestraint(self.model, dps, (pt0.get_index(), pt1.get_index()))
1697 
1698  print("Adding sequence connectivity restraint between", pt0.get_name(), " and ", pt1.get_name(), 'of distance', optdist)
1699  sortedsegments_cr.add_restraint(r)
1700  self.linker_restraints_dict[
1701  "LinkerRestraint-" + pt0.get_name() + "-" + pt1.get_name()] = r
1702  self.connected_intra_pairs.append((pt0, pt1))
1703  self.connected_intra_pairs.append((pt1, pt0))
1704 
1705  self.linker_restraints.add_restraint(sortedsegments_cr)
1706  self.linker_restraints.add_restraint(unmodeledregions_cr)
1707  IMP.pmi1.tools.add_restraint_to_model(self.model, self.linker_restraints)
1708  self.sortedsegments_cr_dict[name] = sortedsegments_cr
1709  self.unmodeledregions_cr_dict[name] = unmodeledregions_cr
1710 
1711  def optimize_floppy_bodies(self, nsteps, temperature=1.0):
1712  import IMP.pmi1.samplers
1713  pts = IMP.pmi1.tools.ParticleToSampleList()
1714  for n, fb in enumerate(self.floppy_bodies):
1715  pts.add_particle(fb, "Floppy_Bodies", 1.0, "Floppy_Body_" + str(n))
1716  if len(pts.get_particles_to_sample()) > 0:
1717  mc = IMP.pmi1.samplers.MonteCarlo(self.model, [pts], temperature)
1718  print("optimize_floppy_bodies: optimizing %i floppy bodies" % len(self.floppy_bodies))
1719  mc.optimize(nsteps)
1720  else:
1721  print("optimize_floppy_bodies: no particle to optimize")
1722 
1723  def create_rotational_symmetry(self, maincopy, copies, rotational_axis=IMP.algebra.Vector3D(0, 0, 1.0),
1724  nSymmetry=None, skip_gaussian_in_clones=False):
1725  '''
1726  The copies must not contain rigid bodies.
1727  The symmetry restraints are applied at each leaf.
1728  '''
1729 
1730  from math import pi
1731  self.representation_is_modified = True
1732  ncopies = len(copies) + 1
1733  main_hiers = IMP.atom.get_leaves(self.hier_dict[maincopy])
1734 
1735  for k in range(len(copies)):
1736  if (nSymmetry is None):
1737  rotation_angle = 2.0 * pi / float(ncopies) * float(k + 1)
1738  else:
1739  if ( k % 2 == 0 ):
1740  rotation_angle = 2.0 * pi / float(nSymmetry) * float((k + 2) / 2)
1741  else:
1742  rotation_angle = -2.0 * pi / float(nSymmetry) * float((k + 1) / 2)
1743 
1744  rotation3D = IMP.algebra.get_rotation_about_axis(rotational_axis, rotation_angle)
1745 
1746  sm = IMP.core.TransformationSymmetry(rotation3D)
1747  clone_hiers = IMP.atom.get_leaves(self.hier_dict[copies[k]])
1748 
1749  lc = IMP.container.ListSingletonContainer(self.model)
1750  for n, p in enumerate(main_hiers):
1751  if (skip_gaussian_in_clones):
1753  continue
1754  pc = clone_hiers[n]
1755  #print("setting " + p.get_name() + " as reference for " + pc.get_name())
1756  IMP.core.Reference.setup_particle(pc.get_particle(), p.get_particle())
1757  lc.add(pc.get_particle().get_index())
1758 
1759  c = IMP.container.SingletonsConstraint(sm, None, lc)
1760  self.model.add_score_state(c)
1761  print("Completed setting " + str(maincopy) + " as a reference for "
1762  + str(copies[k]) + " by rotating it by "
1763  + str(rotation_angle / 2.0 / pi * 360)
1764  + " degrees around the " + str(rotational_axis) + " axis.")
1765  self.model.update()
1766 
1767  def create_rigid_body_symmetry(self, particles_reference, particles_copy,label="None",
1768  initial_transformation=IMP.algebra.get_identity_transformation_3d()):
1769  from math import pi
1770  self.representation_is_modified = True
1771 
1772  mainparticles = particles_reference
1773 
1774  t=initial_transformation
1775  p=IMP.Particle(self.model)
1776  p.set_name("RigidBody_Symmetry")
1778 
1780 
1781 
1782  copyparticles = particles_copy
1783 
1784  mainpurged = []
1785  copypurged = []
1786  for n, p in enumerate(mainparticles):
1787  print(p.get_name())
1788  pc = copyparticles[n]
1789 
1790  mainpurged.append(p)
1793  else:
1794  IMP.pmi1.Symmetric(p).set_symmetric(0)
1795 
1796  copypurged.append(pc)
1799  else:
1800  IMP.pmi1.Symmetric(pc).set_symmetric(1)
1801 
1802  lc = IMP.container.ListSingletonContainer(self.model)
1803  for n, p in enumerate(mainpurged):
1804 
1805  pc = copypurged[n]
1806  print("setting " + p.get_name() + " as reference for " + pc.get_name())
1807 
1809  lc.add(pc.get_index())
1810 
1811  c = IMP.container.SingletonsConstraint(sm, None, lc)
1812  self.model.add_score_state(c)
1813 
1814  self.model.update()
1815  self.rigid_bodies.append(rb)
1816  self.rigid_body_symmetries.append(rb)
1817  rb.set_name(label+".rigid_body_symmetry."+str(len(self.rigid_body_symmetries)))
1818 
1819 
1820  def create_amyloid_fibril_symmetry(self, maincopy, axial_copies,
1821  longitudinal_copies, axis=(0, 0, 1), translation_value=4.8):
1822 
1823  from math import pi
1824  self.representation_is_modified = True
1825 
1826  outhiers = []
1827  protein_h = self.hier_dict[maincopy]
1828  mainparts = IMP.atom.get_leaves(protein_h)
1829 
1830  for ilc in range(-longitudinal_copies, longitudinal_copies + 1):
1831  for iac in range(axial_copies):
1832  copyname = maincopy + "_a" + str(ilc) + "_l" + str(iac)
1833  self.create_component(copyname, 0.0)
1834  for hier in protein_h.get_children():
1835  self.add_component_hierarchy_clone(copyname, hier)
1836  copyhier = self.hier_dict[copyname]
1837  outhiers.append(copyhier)
1838  copyparts = IMP.atom.get_leaves(copyhier)
1840  IMP.algebra.Vector3D(axis),
1841  2 * pi / axial_copies * (float(iac)))
1842  translation_vector = tuple(
1843  [translation_value * float(ilc) * x for x in axis])
1844  print(translation_vector)
1845  translation = IMP.algebra.Vector3D(translation_vector)
1847  IMP.algebra.Transformation3D(rotation3D, translation))
1848  lc = IMP.container.ListSingletonContainer(self.model)
1849  for n, p in enumerate(mainparts):
1850  pc = copyparts[n]
1856  lc.add(pc.get_index())
1857  c = IMP.container.SingletonsConstraint(sm, None, lc)
1858  self.model.add_score_state(c)
1859  self.model.update()
1860  return outhiers
1861 
1862  def link_components_to_rmf(self, rmfname, frameindex):
1863  '''
1864  Load coordinates in the current representation.
1865  This should be done only if the hierarchy self.prot is identical
1866  to the one as stored in the rmf i.e. all components were added.
1867  '''
1868  import IMP.rmf
1869  import RMF
1870 
1871  rh = RMF.open_rmf_file_read_only(rmfname)
1872  IMP.rmf.link_hierarchies(rh, [self.prot])
1873  IMP.rmf.load_frame(rh, RMF.FrameID(frameindex))
1874  del rh
1875 
1876  def create_components_from_rmf(self, rmfname, frameindex):
1877  '''
1878  still not working.
1879  create the representation (i.e. hierarchies) from the rmf file.
1880  it will be stored in self.prot, which will be overwritten.
1881  load the coordinates from the rmf file at frameindex.
1882  '''
1883  rh = RMF.open_rmf_file_read_only(rmfname)
1884  self.prot = IMP.rmf.create_hierarchies(rh, self.model)[0]
1886  IMP.rmf.link_hierarchies(rh, [self.prot])
1887  IMP.rmf.load_frame(rh, RMF.FrameID(frameindex))
1888  del rh
1889  for p in self.prot.get_children():
1890  self.create_component(p.get_name())
1891  self.hier_dict[p.get_name()] = p
1892  '''
1893  still missing: save rigid bodies contained in the rmf in self.rigid_bodies
1894  save floppy bodies in self.floppy_bodies
1895  get the connectivity restraints
1896  '''
1897 
1898  def set_rigid_body_from_hierarchies(self, hiers, particles=None):
1899  '''
1900  Construct a rigid body from hierarchies (and optionally particles).
1901 
1902  @param hiers list of hierarchies to make rigid
1903  @param particles list of particles to add to the rigid body
1904  '''
1905 
1906  if particles is None:
1907  rigid_parts = set()
1908  else:
1909  rigid_parts = set(particles)
1910 
1911  name = ""
1912  print("set_rigid_body_from_hierarchies> setting up a new rigid body")
1913  for hier in hiers:
1914  ps = IMP.atom.get_leaves(hier)
1915 
1916  for p in ps:
1918  rb = IMP.core.RigidMember(p).get_rigid_body()
1919  print("set_rigid_body_from_hierarchies> WARNING particle %s already belongs to rigid body %s" % (p.get_name(), rb.get_name()))
1920  else:
1921  rigid_parts.add(p)
1922  name += hier.get_name() + "-"
1923  print("set_rigid_body_from_hierarchies> adding %s to the rigid body" % hier.get_name())
1924 
1925  if len(list(rigid_parts)) != 0:
1926  rb = IMP.atom.create_rigid_body(list(rigid_parts))
1927  rb.set_coordinates_are_optimized(True)
1928  rb.set_name(name + "rigid_body")
1929  self.rigid_bodies.append(rb)
1930  return rb
1931 
1932  else:
1933  print("set_rigid_body_from_hierarchies> rigid body could not be setup")
1934 
1935  def set_rigid_bodies(self, subunits):
1936  '''
1937  Construct a rigid body from a list of subunits.
1938 
1939  Each subunit is a tuple that identifies the residue ranges and the
1940  component name (as used in create_component()).
1941 
1942  subunits: [(name_1,(first_residue_1,last_residue_1)),(name_2,(first_residue_2,last_residue_2)),.....]
1943  or
1944  [name_1,name_2,(name_3,(first_residue_3,last_residue_3)),.....]
1945 
1946  example: ["prot1","prot2",("prot3",(1,10))]
1947 
1948  sometimes, we know about structure of an interaction
1949  and here we make such PPIs rigid
1950  '''
1951 
1952  rigid_parts = set()
1953  for s in subunits:
1954  if type(s) == type(tuple()) and len(s) == 2:
1955  sel = IMP.atom.Selection(
1956  self.prot,
1957  molecule=s[0],
1958  residue_indexes=list(range(s[1][0],
1959  s[1][1] + 1)))
1960  if len(sel.get_selected_particles()) == 0:
1961  print("set_rigid_bodies: selected particle does not exist")
1962  for p in sel.get_selected_particles():
1963  # if not p in self.floppy_bodies:
1965  rb = IMP.core.RigidMember(p).get_rigid_body()
1966  print("set_rigid_body_from_hierarchies> WARNING particle %s already belongs to rigid body %s" % (p.get_name(), rb.get_name()))
1967  else:
1968  rigid_parts.add(p)
1969 
1970  elif type(s) == type(str()):
1971  sel = IMP.atom.Selection(self.prot, molecule=s)
1972  if len(sel.get_selected_particles()) == 0:
1973  print("set_rigid_bodies: selected particle does not exist")
1974  for p in sel.get_selected_particles():
1975  # if not p in self.floppy_bodies:
1977  rb = IMP.core.RigidMember(p).get_rigid_body()
1978  print("set_rigid_body_from_hierarchies> WARNING particle %s already belongs to rigid body %s" % (p.get_name(), rb.get_name()))
1979  else:
1980  rigid_parts.add(p)
1981 
1982  rb = IMP.atom.create_rigid_body(list(rigid_parts))
1983  rb.set_coordinates_are_optimized(True)
1984  rb.set_name(''.join(str(subunits)) + "_rigid_body")
1985  self.rigid_bodies.append(rb)
1986  return rb
1987 
1988  def set_super_rigid_body_from_hierarchies(
1989  self,
1990  hiers,
1991  axis=None,
1992  min_size=1):
1993  # axis is the rotation axis for 2D rotation
1994  super_rigid_xyzs = set()
1995  super_rigid_rbs = set()
1996  name = ""
1997  print("set_super_rigid_body_from_hierarchies> setting up a new SUPER rigid body")
1998 
1999  for hier in hiers:
2000  ps = IMP.atom.get_leaves(hier)
2001  for p in ps:
2003  rb = IMP.core.RigidMember(p).get_rigid_body()
2004  super_rigid_rbs.add(rb)
2005  else:
2006  super_rigid_xyzs.add(p)
2007  print("set_rigid_body_from_hierarchies> adding %s to the rigid body" % hier.get_name())
2008  if len(super_rigid_rbs|super_rigid_xyzs) < min_size:
2009  return
2010  if axis is None:
2011  self.super_rigid_bodies.append((super_rigid_xyzs, super_rigid_rbs))
2012  else:
2013  # these will be 2D rotation SRB or a bond rotamer (axis can be a IMP.algebra.Vector3D or particle Pair)
2014  self.super_rigid_bodies.append(
2015  (super_rigid_xyzs, super_rigid_rbs, axis))
2016 
2017  def fix_rigid_bodies(self, rigid_bodies):
2018  self.fixed_rigid_bodies += rigid_bodies
2019 
2020 
2022  self, hiers, min_length=None, max_length=None, axis=None):
2023  '''
2024  Make a chain of super rigid bodies from a list of hierarchies.
2025 
2026  Takes a linear list of hierarchies (they are supposed to be
2027  sequence-contiguous) and produces a chain of super rigid bodies
2028  with given length range, specified by min_length and max_length.
2029  '''
2030  try:
2031  hiers = IMP.pmi1.tools.flatten_list(hiers)
2032  except:
2033  pass
2034  for hs in IMP.pmi1.tools.sublist_iterator(hiers, min_length, max_length):
2035  self.set_super_rigid_body_from_hierarchies(hs, axis, min_length)
2036 
2037  def set_super_rigid_bodies(self, subunits, coords=None):
2038  super_rigid_xyzs = set()
2039  super_rigid_rbs = set()
2040 
2041  for s in subunits:
2042  if type(s) == type(tuple()) and len(s) == 3:
2043  sel = IMP.atom.Selection(
2044  self.prot,
2045  molecule=s[2],
2046  residue_indexes=list(range(s[0],
2047  s[1] + 1)))
2048  if len(sel.get_selected_particles()) == 0:
2049  print("set_rigid_bodies: selected particle does not exist")
2050  for p in sel.get_selected_particles():
2052  rb = IMP.core.RigidMember(p).get_rigid_body()
2053  super_rigid_rbs.add(rb)
2054  else:
2055  super_rigid_xyzs.add(p)
2056  elif type(s) == type(str()):
2057  sel = IMP.atom.Selection(self.prot, molecule=s)
2058  if len(sel.get_selected_particles()) == 0:
2059  print("set_rigid_bodies: selected particle does not exist")
2060  for p in sel.get_selected_particles():
2061  # if not p in self.floppy_bodies:
2063  rb = IMP.core.RigidMember(p).get_rigid_body()
2064  super_rigid_rbs.add(rb)
2065  else:
2066  super_rigid_xyzs.add(p)
2067  self.super_rigid_bodies.append((super_rigid_xyzs, super_rigid_rbs))
2068 
2069  def remove_floppy_bodies_from_component(self, component_name):
2070  '''
2071  Remove leaves of hierarchies from the floppy bodies list based
2072  on the component name
2073  '''
2074  hs=IMP.atom.get_leaves(self.hier_dict[component_name])
2075  ps=[h.get_particle() for h in hs]
2076  for p in self.floppy_bodies:
2077  try:
2078  if p in ps: self.floppy_bodies.remove(p)
2079  if p in hs: self.floppy_bodies.remove(p)
2080  except:
2081  continue
2082 
2083 
2084  def remove_floppy_bodies(self, hierarchies):
2085  '''
2086  Remove leaves of hierarchies from the floppy bodies list.
2087 
2088  Given a list of hierarchies, get the leaves and remove the
2089  corresponding particles from the floppy bodies list. We need this
2090  function because sometimes
2091  we want to constrain the floppy bodies in a rigid body - for instance
2092  when you want to associate a bead with a density particle.
2093  '''
2094  for h in hierarchies:
2095  ps = IMP.atom.get_leaves(h)
2096  for p in ps:
2097  if p in self.floppy_bodies:
2098  print("remove_floppy_bodies: removing %s from floppy body list" % p.get_name())
2099  self.floppy_bodies.remove(p)
2100 
2101 
2102  def set_floppy_bodies(self):
2103  for p in self.floppy_bodies:
2104  name = p.get_name()
2105  p.set_name(name + "_floppy_body")
2107  print("I'm trying to make this particle flexible although it was assigned to a rigid body", p.get_name())
2108  rb = IMP.core.RigidMember(p).get_rigid_body()
2109  try:
2110  rb.set_is_rigid_member(p.get_index(), False)
2111  except:
2112  # some IMP versions still work with that
2113  rb.set_is_rigid_member(p.get_particle_index(), False)
2114  p.set_name(p.get_name() + "_rigid_body_member")
2115 
2116  def set_floppy_bodies_from_hierarchies(self, hiers):
2117  for hier in hiers:
2118  ps = IMP.atom.get_leaves(hier)
2119  for p in ps:
2120  IMP.core.XYZ(p).set_coordinates_are_optimized(True)
2121  self.floppy_bodies.append(p)
2122 
2124  self,
2125  selection_tuples,
2126  resolution=None):
2127  '''
2128  selection tuples must be [(r1,r2,"name1"),(r1,r2,"name2"),....]
2129  @return the particles
2130  '''
2131  particles = []
2132  print(selection_tuples)
2133  for s in selection_tuples:
2134  ps = IMP.pmi1.tools.select_by_tuple(
2135  representation=self, tupleselection=tuple(s),
2136  resolution=None, name_is_ambiguous=False)
2137  particles += ps
2138  return particles
2139 
2140  def get_connected_intra_pairs(self):
2141  return self.connected_intra_pairs
2142 
2143  def set_rigid_bodies_max_trans(self, maxtrans):
2144  self.maxtrans_rb = maxtrans
2145 
2146  def set_rigid_bodies_max_rot(self, maxrot):
2147  self.maxrot_rb = maxrot
2148 
2149  def set_super_rigid_bodies_max_trans(self, maxtrans):
2150  self.maxtrans_srb = maxtrans
2151 
2152  def set_super_rigid_bodies_max_rot(self, maxrot):
2153  self.maxrot_srb = maxrot
2154 
2155  def set_floppy_bodies_max_trans(self, maxtrans):
2156  self.maxtrans_fb = maxtrans
2157 
2158  def set_rigid_bodies_as_fixed(self, rigidbodiesarefixed=True):
2159  '''
2160  Fix rigid bodies in their actual position.
2161  The get_particles_to_sample() function will return
2162  just the floppy bodies (if they are not fixed).
2163  '''
2164  self.rigidbodiesarefixed = rigidbodiesarefixed
2165 
2166  def set_floppy_bodies_as_fixed(self, floppybodiesarefixed=True):
2167  '''
2168  Fix floppy bodies in their actual position.
2169  The get_particles_to_sample() function will return
2170  just the rigid bodies (if they are not fixed).
2171  '''
2172  self.floppybodiesarefixed=floppybodiesarefixed
2173 
2174  def draw_hierarchy_graph(self):
2175  for c in IMP.atom.Hierarchy(self.prot).get_children():
2176  print("Drawing hierarchy graph for " + c.get_name())
2177  IMP.pmi1.output.get_graph_from_hierarchy(c)
2178 
2179  def get_geometries(self):
2180  # create segments at the lowest resolution
2181  seggeos = []
2182  for name in self.hier_geometry_pairs:
2183  for pt in self.hier_geometry_pairs[name]:
2184  p1 = pt[0]
2185  p2 = pt[1]
2186  color = pt[2]
2187  try:
2188  clr = IMP.display.get_rgb_color(color)
2189  except:
2190  clr = IMP.display.get_rgb_color(1.0)
2191  coor1 = IMP.core.XYZ(p1).get_coordinates()
2192  coor2 = IMP.core.XYZ(p2).get_coordinates()
2193  seg = IMP.algebra.Segment3D(coor1, coor2)
2194  seggeos.append(IMP.display.SegmentGeometry(seg, clr))
2195  return seggeos
2196 
2197  def setup_bonds(self):
2198  # create segments at the lowest resolution
2199  seggeos = []
2200  for name in self.hier_geometry_pairs:
2201  for pt in self.hier_geometry_pairs[name]:
2202 
2203  p1 = pt[0]
2204  p2 = pt[1]
2205  if not IMP.atom.Bonded.get_is_setup(p1):
2207  if not IMP.atom.Bonded.get_is_setup(p2):
2209 
2212  IMP.atom.Bonded(p1),
2213  IMP.atom.Bonded(p2),
2214  1)
2215 
2216  def show_component_table(self, name):
2217  if name in self.sequence_dict:
2218  lastresn = len(self.sequence_dict[name])
2219  firstresn = 1
2220  else:
2221  residues = self.hier_db.get_residue_numbers(name)
2222  firstresn = min(residues)
2223  lastresn = max(residues)
2224 
2225  for nres in range(firstresn, lastresn + 1):
2226  try:
2227  resolutions = self.hier_db.get_residue_resolutions(name, nres)
2228  ps = []
2229  for r in resolutions:
2230  ps += self.hier_db.get_particles(name, nres, r)
2231  print("%20s %7s" % (name, nres), " ".join(["%20s %7s" % (str(p.get_name()),
2232  str(IMP.pmi1.Resolution(p).get_resolution())) for p in ps]))
2233  except:
2234  print("%20s %20s" % (name, nres), "**** not represented ****")
2235 
2236  def draw_hierarchy_composition(self):
2237 
2238  ks = list(self.elements.keys())
2239  ks.sort()
2240 
2241  max = 0
2242  for k in self.elements:
2243  for l in self.elements[k]:
2244  if l[1] > max:
2245  max = l[1]
2246 
2247  for k in ks:
2248  self.draw_component_composition(k, max)
2249 
2250  def draw_component_composition(self, name, max=1000, draw_pdb_names=False):
2251  from matplotlib import pyplot
2252  import matplotlib as mpl
2253  k = name
2254  tmplist = sorted(self.elements[k], key=itemgetter(0))
2255  try:
2256  endres = tmplist[-1][1]
2257  except:
2258  print("draw_component_composition: missing information for component %s" % name)
2259  return
2260  fig = pyplot.figure(figsize=(26.0 * float(endres) / max + 2, 2))
2261  ax = fig.add_axes([0.05, 0.475, 0.9, 0.15])
2262 
2263  # Set the colormap and norm to correspond to the data for which
2264  # the colorbar will be used.
2265  cmap = mpl.cm.cool
2266  norm = mpl.colors.Normalize(vmin=5, vmax=10)
2267  bounds = [1]
2268  colors = []
2269 
2270  for n, l in enumerate(tmplist):
2271  firstres = l[0]
2272  lastres = l[1]
2273  if l[3] != "end":
2274  if bounds[-1] != l[0]:
2275  colors.append("white")
2276  bounds.append(l[0])
2277  if l[3] == "pdb":
2278  colors.append("#99CCFF")
2279  if l[3] == "bead":
2280  colors.append("#FFFF99")
2281  if l[3] == "helix":
2282  colors.append("#33CCCC")
2283  if l[3] != "end":
2284  bounds.append(l[1] + 1)
2285  else:
2286  if l[3] == "pdb":
2287  colors.append("#99CCFF")
2288  if l[3] == "bead":
2289  colors.append("#FFFF99")
2290  if l[3] == "helix":
2291  colors.append("#33CCCC")
2292  if l[3] != "end":
2293  bounds.append(l[1] + 1)
2294  else:
2295  if bounds[-1] - 1 == l[0]:
2296  bounds.pop()
2297  bounds.append(l[0])
2298  else:
2299  colors.append("white")
2300  bounds.append(l[0])
2301 
2302  bounds.append(bounds[-1])
2303  colors.append("white")
2304  cmap = mpl.colors.ListedColormap(colors)
2305  cmap.set_over('0.25')
2306  cmap.set_under('0.75')
2307 
2308  norm = mpl.colors.BoundaryNorm(bounds, cmap.N)
2309  cb2 = mpl.colorbar.ColorbarBase(ax, cmap=cmap,
2310  norm=norm,
2311  # to use 'extend', you must
2312  # specify two extra boundaries:
2313  boundaries=bounds,
2314  ticks=bounds, # optional
2315  spacing='proportional',
2316  orientation='horizontal')
2317 
2318  extra_artists = []
2319  npdb = 0
2320 
2321  if draw_pdb_names:
2322  for l in tmplist:
2323  if l[3] == "pdb":
2324  npdb += 1
2325  mid = 1.0 / endres * float(l[0])
2326  # t =ax.text(mid, float(npdb-1)/2.0+1.5, l[2], ha="left", va="center", rotation=0,
2327  # size=10)
2328  # t=ax.annotate(l[0],2)
2329  t = ax.annotate(
2330  l[2], xy=(mid, 1), xycoords='axes fraction',
2331  xytext=(mid + 0.025, float(npdb - 1) / 2.0 + 1.5), textcoords='axes fraction',
2332  arrowprops=dict(arrowstyle="->",
2333  connectionstyle="angle,angleA=0,angleB=90,rad=10"),
2334  )
2335  extra_artists.append(t)
2336 
2337  # set the title of the bar
2338  title = ax.text(-0.005, 0.5, k, ha="right", va="center", rotation=90,
2339  size=20)
2340 
2341  extra_artists.append(title)
2342  # changing the xticks labels
2343  labels = len(bounds) * [" "]
2344  ax.set_xticklabels(labels)
2345  mid = 1.0 / endres * float(bounds[0])
2346  t = ax.annotate(bounds[0], xy=(mid, 0), xycoords='axes fraction',
2347  xytext=(mid - 0.01, -0.5), textcoords='axes fraction',)
2348  extra_artists.append(t)
2349  offsets = [0, -0.5, -1.0]
2350  nclashes = 0
2351  for n in range(1, len(bounds)):
2352  if bounds[n] == bounds[n - 1]:
2353  continue
2354  mid = 1.0 / endres * float(bounds[n])
2355  if (float(bounds[n]) - float(bounds[n - 1])) / max <= 0.01:
2356  nclashes += 1
2357  offset = offsets[nclashes % 3]
2358  else:
2359  nclashes = 0
2360  offset = offsets[0]
2361  if offset > -0.75:
2362  t = ax.annotate(
2363  bounds[n], xy=(mid, 0), xycoords='axes fraction',
2364  xytext=(mid, -0.5 + offset), textcoords='axes fraction')
2365  else:
2366  t = ax.annotate(
2367  bounds[n], xy=(mid, 0), xycoords='axes fraction',
2368  xytext=(mid, -0.5 + offset), textcoords='axes fraction', arrowprops=dict(arrowstyle="-"))
2369  extra_artists.append(t)
2370 
2371  cb2.add_lines(bounds, ["black"] * len(bounds), [1] * len(bounds))
2372  # cb2.set_label(k)
2373 
2374  pyplot.savefig(
2375  k + "structure.pdf",
2376  dpi=150,
2377  transparent="True",
2378  bbox_extra_artists=(extra_artists),
2379  bbox_inches='tight')
2380  pyplot.show()
2381 
2382  def draw_coordinates_projection(self):
2383  import matplotlib.pyplot as pp
2384  xpairs = []
2385  ypairs = []
2386  for name in self.hier_geometry_pairs:
2387  for pt in self.hier_geometry_pairs[name]:
2388  p1 = pt[0]
2389  p2 = pt[1]
2390  color = pt[2]
2391  coor1 = IMP.core.XYZ(p1).get_coordinates()
2392  coor2 = IMP.core.XYZ(p2).get_coordinates()
2393  x1 = coor1[0]
2394  x2 = coor2[0]
2395  y1 = coor1[1]
2396  y2 = coor2[1]
2397  xpairs.append([x1, x2])
2398  ypairs.append([y1, y2])
2399  xlist = []
2400  ylist = []
2401  for xends, yends in zip(xpairs, ypairs):
2402  xlist.extend(xends)
2403  xlist.append(None)
2404  ylist.extend(yends)
2405  ylist.append(None)
2406  pp.plot(xlist, ylist, 'b-', alpha=0.1)
2407  pp.show()
2408 
2409  def get_prot_name_from_particle(self, particle):
2410  names = self.get_component_names()
2411  particle0 = particle
2412  name = None
2413  while not name in names:
2414  h = IMP.atom.Hierarchy(particle0).get_parent()
2415  name = h.get_name()
2416  particle0 = h.get_particle()
2417  return name
2418 
2419 
2420  def get_particles_to_sample(self):
2421  # get the list of samplable particles with their type
2422  # and the mover displacement. Everything wrapped in a dictionary,
2423  # to be used by samplers modules
2424  ps = {}
2425 
2426  # remove symmetric particles: they are not sampled
2427  rbtmp = []
2428  fbtmp = []
2429  srbtmp = []
2430  if not self.rigidbodiesarefixed:
2431  for rb in self.rigid_bodies:
2433  if IMP.pmi1.Symmetric(rb).get_symmetric() != 1:
2434  rbtmp.append(rb)
2435  else:
2436  if rb not in self.fixed_rigid_bodies:
2437  rbtmp.append(rb)
2438 
2439  if not self.floppybodiesarefixed:
2440  for fb in self.floppy_bodies:
2442  if IMP.pmi1.Symmetric(fb).get_symmetric() != 1:
2443  fbtmp.append(fb)
2444  else:
2445  fbtmp.append(fb)
2446 
2447  for srb in self.super_rigid_bodies:
2448  # going to prune the fixed rigid bodies out
2449  # of the super rigid body list
2450  rigid_bodies = list(srb[1])
2451  filtered_rigid_bodies = []
2452  for rb in rigid_bodies:
2453  if rb not in self.fixed_rigid_bodies:
2454  filtered_rigid_bodies.append(rb)
2455  srbtmp.append((srb[0], filtered_rigid_bodies))
2456 
2457  self.rigid_bodies = rbtmp
2458  self.floppy_bodies = fbtmp
2459  self.super_rigid_bodies = srbtmp
2460 
2461  ps["Rigid_Bodies_SimplifiedModel"] = (
2462  self.rigid_bodies,
2463  self.maxtrans_rb,
2464  self.maxrot_rb)
2465  ps["Floppy_Bodies_SimplifiedModel"] = (
2466  self.floppy_bodies,
2467  self.maxtrans_fb)
2468  ps["SR_Bodies_SimplifiedModel"] = (
2469  self.super_rigid_bodies,
2470  self.maxtrans_srb,
2471  self.maxrot_srb)
2472  return ps
2473 
2474  def set_output_level(self, level):
2475  self.output_level = level
2476 
2477  def _evaluate(self, deriv):
2478  """Evaluate the total score of all added restraints"""
2479  r = IMP.pmi1.tools.get_restraint_set(self.model)
2480  return r.evaluate(deriv)
2481 
2482  def get_output(self):
2483  output = {}
2484  score = 0.0
2485 
2486  output["SimplifiedModel_Total_Score_" +
2487  self.label] = str(self._evaluate(False))
2488  output["SimplifiedModel_Linker_Score_" +
2489  self.label] = str(self.linker_restraints.unprotected_evaluate(None))
2490  for name in self.sortedsegments_cr_dict:
2491  partialscore = self.sortedsegments_cr_dict[name].evaluate(False)
2492  score += partialscore
2493  output[
2494  "SimplifiedModel_Link_SortedSegments_" +
2495  name +
2496  "_" +
2497  self.label] = str(
2498  partialscore)
2499  partialscore = self.unmodeledregions_cr_dict[name].evaluate(False)
2500  score += partialscore
2501  output[
2502  "SimplifiedModel_Link_UnmodeledRegions_" +
2503  name +
2504  "_" +
2505  self.label] = str(
2506  partialscore)
2507  for rb in self.rigid_body_symmetries:
2508  name=rb.get_name()
2509  output[name +"_" +self.label]=str(rb.get_reference_frame().get_transformation_to())
2510  for name in self.linker_restraints_dict:
2511  output[
2512  name +
2513  "_" +
2514  self.label] = str(
2515  self.linker_restraints_dict[
2516  name].unprotected_evaluate(
2517  None))
2518 
2519  if len(self.reference_structures.keys()) != 0:
2520  rmsds = self.get_all_rmsds()
2521  for label in rmsds:
2522  output[
2523  "SimplifiedModel_" +
2524  label +
2525  "_" +
2526  self.label] = rmsds[
2527  label]
2528 
2529  if self.output_level == "high":
2530  for p in IMP.atom.get_leaves(self.prot):
2531  d = IMP.core.XYZR(p)
2532  output["Coordinates_" +
2533  p.get_name() + "_" + self.label] = str(d)
2534 
2535  output["_TotalScore"] = str(score)
2536  return output
2537 
2538  def get_test_output(self):
2539  # this method is called by test functions and return an enriched output
2540  output = self.get_output()
2541  for n, p in enumerate(self.get_particles_to_sample()):
2542  output["Particle_to_sample_" + str(n)] = str(p)
2543  output["Number_of_particles"] = len(IMP.atom.get_leaves(self.prot))
2544  output["Hierarchy_Dictionary"] = list(self.hier_dict.keys())
2545  output["Number_of_floppy_bodies"] = len(self.floppy_bodies)
2546  output["Number_of_rigid_bodies"] = len(self.rigid_bodies)
2547  output["Number_of_super_bodies"] = len(self.super_rigid_bodies)
2548  output["Selection_resolution_1"] = len(
2549  IMP.pmi1.tools.select(self, resolution=1))
2550  output["Selection_resolution_5"] = len(
2551  IMP.pmi1.tools.select(self, resolution=5))
2552  output["Selection_resolution_7"] = len(
2553  IMP.pmi1.tools.select(self, resolution=5))
2554  output["Selection_resolution_10"] = len(
2555  IMP.pmi1.tools.select(self, resolution=10))
2556  output["Selection_resolution_100"] = len(
2557  IMP.pmi1.tools.select(self, resolution=100))
2558  output["Selection_All"] = len(IMP.pmi1.tools.select(self))
2559  output["Selection_resolution=1"] = len(
2560  IMP.pmi1.tools.select(self, resolution=1))
2561  output["Selection_resolution=1,resid=10"] = len(
2562  IMP.pmi1.tools.select(self, resolution=1, residue=10))
2563  for resolution in self.hier_resolution:
2564  output["Hier_resolution_dictionary" +
2565  str(resolution)] = len(self.hier_resolution[resolution])
2566  for name in self.hier_dict:
2567  output[
2568  "Selection_resolution=1,resid=10,name=" +
2569  name] = len(
2571  self,
2572  resolution=1,
2573  name=name,
2574  residue=10))
2575  output[
2576  "Selection_resolution=1,resid=10,name=" +
2577  name +
2578  ",ambiguous"] = len(
2580  self,
2581  resolution=1,
2582  name=name,
2583  name_is_ambiguous=True,
2584  residue=10))
2585  output[
2586  "Selection_resolution=1,resid=10,name=" +
2587  name +
2588  ",ambiguous"] = len(
2590  self,
2591  resolution=1,
2592  name=name,
2593  name_is_ambiguous=True,
2594  residue=10))
2595  output[
2596  "Selection_resolution=1,resrange=(10,20),name=" +
2597  name] = len(
2599  self,
2600  resolution=1,
2601  name=name,
2602  first_residue=10,
2603  last_residue=20))
2604  output[
2605  "Selection_resolution=1,resrange=(10,20),name=" +
2606  name +
2607  ",ambiguous"] = len(
2609  self,
2610  resolution=1,
2611  name=name,
2612  name_is_ambiguous=True,
2613  first_residue=10,
2614  last_residue=20))
2615  output[
2616  "Selection_resolution=10,resrange=(10,20),name=" +
2617  name] = len(
2619  self,
2620  resolution=10,
2621  name=name,
2622  first_residue=10,
2623  last_residue=20))
2624  output[
2625  "Selection_resolution=10,resrange=(10,20),name=" +
2626  name +
2627  ",ambiguous"] = len(
2629  self,
2630  resolution=10,
2631  name=name,
2632  name_is_ambiguous=True,
2633  first_residue=10,
2634  last_residue=20))
2635  output[
2636  "Selection_resolution=100,resrange=(10,20),name=" +
2637  name] = len(
2639  self,
2640  resolution=100,
2641  name=name,
2642  first_residue=10,
2643  last_residue=20))
2644  output[
2645  "Selection_resolution=100,resrange=(10,20),name=" +
2646  name +
2647  ",ambiguous"] = len(
2649  self,
2650  resolution=100,
2651  name=name,
2652  name_is_ambiguous=True,
2653  first_residue=10,
2654  last_residue=20))
2655 
2656  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:245
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:48
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:616
Add resolution to a particle.
Definition: /Resolution.h:24
static bool get_is_setup(const IMP::ParticleAdaptor &p)
Definition: rigid_bodies.h:617
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:242
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:348
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:731
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:37
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:421
Class for storing model, its restraints, constraints, and particles.
Definition: Model.h:73
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:91
Add uncertainty to a particle.
Definition: /Uncertainty.h:24
A Gaussian distribution in 3D.
Definition: Gaussian3D.h:24
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:10610
A decorator for a particle representing an atom.
Definition: atom/Atom.h:234
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:44
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:115
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:51
Transformation3D get_identity_transformation_3d()
Return a transformation that does not do anything.
Find all nearby pairs by testing all pairs.
Simple implementation of segments in 3D.
Definition: Segment3D.h:24
static bool get_is_setup(const IMP::ParticleAdaptor &p)
Definition: rigid_bodies.h:163
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:421
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:356
void link_hierarchies(RMF::FileConstHandle fh, const atom::Hierarchies &hs)
Class to handle individual particles of a Model object.
Definition: Particle.h:41
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:371
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:29
Classes for writing output files and processing them.
Definition: /output.py:1
Select all CA ATOM records.
Definition: pdb.h:77
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:47
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:64
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:81
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:174
Select all ATOM and HETATM records with the given chain ids.
Definition: pdb.h:189
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:752
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:24
A decorator for a particle with x,y,z coordinates and a radius.
Definition: XYZR.h:27