IMP logo
IMP Reference Guide  2.6.2
The Integrative Modeling Platform
macros.py
1 """@namespace IMP.pmi.macros
2 Protocols for sampling structures and analyzing them.
3 """
4 
5 from __future__ import print_function, division
6 import IMP
8 import IMP.pmi.tools
9 import IMP.pmi.samplers
10 import IMP.pmi.output
11 import IMP.pmi.analysis
12 import IMP.pmi.io
13 import IMP.rmf
14 import IMP.isd
15 import IMP.pmi.dof
16 import RMF
17 import os
18 import glob
19 from operator import itemgetter
20 from collections import defaultdict
21 import numpy as np
22 
23 class ReplicaExchange0(object):
24  """A macro to help setup and run replica exchange.
25  Supports Monte Carlo and molecular dynamics.
26  Produces trajectory RMF files, best PDB structures,
27  and output stat files.
28  """
29  def __init__(self, model,
30  representation=None,
31  root_hier=None,
32  sample_objects=None, # DEPRECATED
33  monte_carlo_sample_objects=None,
34  molecular_dynamics_sample_objects=None,
35  output_objects=[],
36  crosslink_restraints=None,
37  monte_carlo_temperature=1.0,
38  simulated_annealing=False,
39  simulated_annealing_minimum_temperature=1.0,
40  simulated_annealing_maximum_temperature=2.5,
41  simulated_annealing_minimum_temperature_nframes=100,
42  simulated_annealing_maximum_temperature_nframes=100,
43  replica_exchange_minimum_temperature=1.0,
44  replica_exchange_maximum_temperature=2.5,
45  num_sample_rounds=1,
46  number_of_best_scoring_models=500,
47  monte_carlo_steps=10,
48  molecular_dynamics_steps=10,
49  molecular_dynamics_max_time_step=1.0,
50  number_of_frames=1000,
51  nframes_write_coordinates=1,
52  write_initial_rmf=True,
53  initial_rmf_name_suffix="initial",
54  stat_file_name_suffix="stat",
55  best_pdb_name_suffix="model",
56  do_clean_first=True,
57  do_create_directories=True,
58  global_output_directory="./",
59  rmf_dir="rmfs/",
60  best_pdb_dir="pdbs/",
61  replica_stat_file_suffix="stat_replica",
62  em_object_for_rmf=None,
63  atomistic=False,
64  replica_exchange_object=None,
65  test_mode=False):
66  """Constructor.
67  @param model The IMP model
68  @param representation PMI.representation.Representation object
69  (or list of them, for multi-state modeling)
70  @param root_hier Instead of passing Representation, just pass
71  a hierarchy
72  @param monte_carlo_sample_objects Objects for MC sampling; a list of
73  structural components (generally the representation) that will
74  be moved and restraints with parameters that need to
75  be sampled.
76  For PMI2: just pass flat list of movers
77  @param molecular_dynamics_sample_objects Objects for MD sampling
78  For PMI2: just pass flat list of particles
79  @param output_objects A list of structural objects and restraints
80  that will be included in output (statistics files). Any object
81  that provides a get_output() method can be used here.
82  @param crosslink_restraints List of cross-link restraints that will
83  be included in output RMF files (for visualization).
84  @param monte_carlo_temperature MC temp (may need to be optimized
85  based on post-sampling analysis)
86  @param simulated_annealing If True, perform simulated annealing
87  @param simulated_annealing_minimum_temperature Should generally be
88  the same as monte_carlo_temperature.
89  @param simulated_annealing_minimum_temperature_nframes Number of
90  frames to compute at minimum temperature.
91  @param simulated_annealing_maximum_temperature_nframes Number of
92  frames to compute at
93  temps > simulated_annealing_maximum_temperature.
94  @param replica_exchange_minimum_temperature Low temp for REX; should
95  generally be the same as monte_carlo_temperature.
96  @param replica_exchange_maximum_temperature High temp for REX
97  @param num_sample_rounds Number of rounds of MC/MD per cycle
98  @param number_of_best_scoring_models Number of top-scoring PDB models
99  to keep around for analysis
100  @param monte_carlo_steps Number of MC steps per round
101  @param molecular_dynamics_steps Number of MD steps per round
102  @param molecular_dynamics_max_time_step Max time step for MD
103  @param number_of_frames Number of REX frames to run
104  @param nframes_write_coordinates How often to write the coordinates
105  of a frame
106  @param write_initial_rmf Write the initial configuration
107  @param global_output_directory Folder that will be created to house
108  output.
109  @param test_mode Set to True to avoid writing any files, just test one frame.
110  """
111  self.model = model
112  self.vars = {}
113 
114  ### add check hierarchy is multistate
115  if representation:
116  if type(representation) == list:
117  self.is_multi_state = True
118  self.root_hiers = [r.prot for r in representation]
119  self.vars["number_of_states"] = len(representation)
120  else:
121  self.is_multi_state = False
122  self.root_hier = representation.prot
123  self.vars["number_of_states"] = 1
124  elif root_hier:
125  states = IMP.atom.get_by_type(root_hier,IMP.atom.STATE_TYPE)
126  self.vars["number_of_states"] = len(states)
127  if len(states)>1:
128  self.root_hiers = states
129  self.is_multi_state = True
130  else:
131  self.root_hier = root_hier
132  self.is_multi_state = False
133  else:
134  print("ERROR: Must provide representation or root_hier")
135  return
136 
137  self.crosslink_restraints = crosslink_restraints
138  self.em_object_for_rmf = em_object_for_rmf
139  self.monte_carlo_sample_objects = monte_carlo_sample_objects
140  if sample_objects is not None:
141  self.monte_carlo_sample_objects+=sample_objects
142  self.molecular_dynamics_sample_objects=molecular_dynamics_sample_objects
143  self.output_objects = output_objects
144  self.replica_exchange_object = replica_exchange_object
145  self.molecular_dynamics_max_time_step = molecular_dynamics_max_time_step
146  self.vars["monte_carlo_temperature"] = monte_carlo_temperature
147  self.vars[
148  "replica_exchange_minimum_temperature"] = replica_exchange_minimum_temperature
149  self.vars[
150  "replica_exchange_maximum_temperature"] = replica_exchange_maximum_temperature
151 
152  self.vars["simulated_annealing"]=\
153  simulated_annealing
154  self.vars["simulated_annealing_minimum_temperature"]=\
155  simulated_annealing_minimum_temperature
156  self.vars["simulated_annealing_maximum_temperature"]=\
157  simulated_annealing_maximum_temperature
158  self.vars["simulated_annealing_minimum_temperature_nframes"]=\
159  simulated_annealing_minimum_temperature_nframes
160  self.vars["simulated_annealing_maximum_temperature_nframes"]=\
161  simulated_annealing_maximum_temperature_nframes
162 
163  self.vars["num_sample_rounds"] = num_sample_rounds
164  self.vars[
165  "number_of_best_scoring_models"] = number_of_best_scoring_models
166  self.vars["monte_carlo_steps"] = monte_carlo_steps
167  self.vars["molecular_dynamics_steps"]=molecular_dynamics_steps
168  self.vars["number_of_frames"] = number_of_frames
169  self.vars["nframes_write_coordinates"] = nframes_write_coordinates
170  self.vars["write_initial_rmf"] = write_initial_rmf
171  self.vars["initial_rmf_name_suffix"] = initial_rmf_name_suffix
172  self.vars["best_pdb_name_suffix"] = best_pdb_name_suffix
173  self.vars["stat_file_name_suffix"] = stat_file_name_suffix
174  self.vars["do_clean_first"] = do_clean_first
175  self.vars["do_create_directories"] = do_create_directories
176  self.vars["global_output_directory"] = global_output_directory
177  self.vars["rmf_dir"] = rmf_dir
178  self.vars["best_pdb_dir"] = best_pdb_dir
179  self.vars["atomistic"] = atomistic
180  self.vars["replica_stat_file_suffix"] = replica_stat_file_suffix
181  self.vars["geometries"] = None
182  self.test_mode = test_mode
183  if root_hier:
184  self.output_objects.append(IMP.pmi.io.TotalScoreOutput(self.model))
185  def add_geometries(self, geometries):
186  if self.vars["geometries"] is None:
187  self.vars["geometries"] = list(geometries)
188  else:
189  self.vars["geometries"].extend(geometries)
190 
191  def show_info(self):
192  print("ReplicaExchange0: it generates initial.*.rmf3, stat.*.out, rmfs/*.rmf3 for each replica ")
193  print("--- it stores the best scoring pdb models in pdbs/")
194  print("--- the stat.*.out and rmfs/*.rmf3 are saved only at the lowest temperature")
195  print("--- variables:")
196  keys = list(self.vars.keys())
197  keys.sort()
198  for v in keys:
199  print("------", v.ljust(30), self.vars[v])
200 
201  def get_replica_exchange_object(self):
202  return self.replica_exchange_object
203 
204  def execute_macro(self):
205  temp_index_factor = 100000.0
206  samplers=[]
207  sampler_mc=None
208  sampler_md=None
209  if self.monte_carlo_sample_objects is not None:
210  print("Setting up MonteCarlo")
211  sampler_mc = IMP.pmi.samplers.MonteCarlo(self.model,
212  self.monte_carlo_sample_objects,
213  self.vars["monte_carlo_temperature"])
214  if self.vars["simulated_annealing"]:
215  tmin=self.vars["simulated_annealing_minimum_temperature"]
216  tmax=self.vars["simulated_annealing_maximum_temperature"]
217  nfmin=self.vars["simulated_annealing_minimum_temperature_nframes"]
218  nfmax=self.vars["simulated_annealing_maximum_temperature_nframes"]
219  sampler_mc.set_simulated_annealing(tmin,tmax,nfmin,nfmax)
220  self.output_objects.append(sampler_mc)
221  samplers.append(sampler_mc)
222 
223 
224  if self.molecular_dynamics_sample_objects is not None:
225  print("Setting up MolecularDynamics")
226  sampler_md = IMP.pmi.samplers.MolecularDynamics(self.model,
227  self.molecular_dynamics_sample_objects,
228  self.vars["monte_carlo_temperature"],
229  maximum_time_step=self.molecular_dynamics_max_time_step)
230  if self.vars["simulated_annealing"]:
231  tmin=self.vars["simulated_annealing_minimum_temperature"]
232  tmax=self.vars["simulated_annealing_maximum_temperature"]
233  nfmin=self.vars["simulated_annealing_minimum_temperature_nframes"]
234  nfmax=self.vars["simulated_annealing_maximum_temperature_nframes"]
235  sampler_md.set_simulated_annealing(tmin,tmax,nfmin,nfmax)
236  self.output_objects.append(sampler_md)
237  samplers.append(sampler_md)
238 # -------------------------------------------------------------------------
239 
240  print("Setting up ReplicaExchange")
241  rex = IMP.pmi.samplers.ReplicaExchange(self.model,
242  self.vars[
243  "replica_exchange_minimum_temperature"],
244  self.vars[
245  "replica_exchange_maximum_temperature"],
246  samplers,
247  replica_exchange_object=self.replica_exchange_object)
248  self.replica_exchange_object = rex.rem
249 
250  myindex = rex.get_my_index()
251  self.output_objects.append(rex)
252 
253  # must reset the minimum temperature due to the
254  # different binary length of rem.get_my_parameter double and python
255  # float
256  min_temp_index = int(min(rex.get_temperatures()) * temp_index_factor)
257 
258 # -------------------------------------------------------------------------
259 
260  globaldir = self.vars["global_output_directory"] + "/"
261  rmf_dir = globaldir + self.vars["rmf_dir"]
262  pdb_dir = globaldir + self.vars["best_pdb_dir"]
263 
264  if not self.test_mode:
265  if self.vars["do_clean_first"]:
266  pass
267 
268  if self.vars["do_create_directories"]:
269 
270  try:
271  os.makedirs(globaldir)
272  except:
273  pass
274  try:
275  os.makedirs(rmf_dir)
276  except:
277  pass
278 
279  if not self.is_multi_state:
280  try:
281  os.makedirs(pdb_dir)
282  except:
283  pass
284  else:
285  for n in range(self.vars["number_of_states"]):
286  try:
287  os.makedirs(pdb_dir + "/" + str(n))
288  except:
289  pass
290 
291 # -------------------------------------------------------------------------
292 
293  sw = IMP.pmi.tools.Stopwatch()
294  self.output_objects.append(sw)
295 
296  print("Setting up stat file")
297  output = IMP.pmi.output.Output(atomistic=self.vars["atomistic"])
298  low_temp_stat_file = globaldir + \
299  self.vars["stat_file_name_suffix"] + "." + str(myindex) + ".out"
300  if not self.test_mode:
301  output.init_stat2(low_temp_stat_file,
302  self.output_objects,
303  extralabels=["rmf_file", "rmf_frame_index"])
304 
305  print("Setting up replica stat file")
306  replica_stat_file = globaldir + \
307  self.vars["replica_stat_file_suffix"] + "." + str(myindex) + ".out"
308  if not self.test_mode:
309  output.init_stat2(replica_stat_file, [rex], extralabels=["score"])
310 
311  if not self.test_mode:
312  print("Setting up best pdb files")
313  if not self.is_multi_state:
314  if self.vars["number_of_best_scoring_models"] > 0:
315  output.init_pdb_best_scoring(pdb_dir + "/" +
316  self.vars["best_pdb_name_suffix"],
317  self.root_hier,
318  self.vars[
319  "number_of_best_scoring_models"],
320  replica_exchange=True)
321  output.write_psf(pdb_dir + "/" +"model.psf",pdb_dir + "/" +
322  self.vars["best_pdb_name_suffix"]+".0.pdb")
323  else:
324  if self.vars["number_of_best_scoring_models"] > 0:
325  for n in range(self.vars["number_of_states"]):
326  output.init_pdb_best_scoring(pdb_dir + "/" + str(n) + "/" +
327  self.vars["best_pdb_name_suffix"],
328  self.root_hiers[n],
329  self.vars[
330  "number_of_best_scoring_models"],
331  replica_exchange=True)
332  output.write_psf(pdb_dir + "/" + str(n) + "/" +"model.psf",pdb_dir + "/" + str(n) + "/" +
333  self.vars["best_pdb_name_suffix"]+".0.pdb")
334 # ---------------------------------------------
335 
336  if not self.em_object_for_rmf is None:
337  if not self.is_multi_state:
338  output_hierarchies = [
339  self.root_hier,
340  self.em_object_for_rmf.get_density_as_hierarchy(
341  )]
342  else:
343  output_hierarchies = self.root_hiers
344  output_hierarchies.append(
345  self.em_object_for_rmf.get_density_as_hierarchy())
346  else:
347  if not self.is_multi_state:
348  output_hierarchies = [self.root_hier]
349  else:
350  output_hierarchies = self.root_hiers
351 
352 #----------------------------------------------
353  if not self.test_mode:
354  print("Setting up and writing initial rmf coordinate file")
355  init_suffix = globaldir + self.vars["initial_rmf_name_suffix"]
356  output.init_rmf(init_suffix + "." + str(myindex) + ".rmf3",
357  output_hierarchies)
358  if self.crosslink_restraints:
359  output.add_restraints_to_rmf(
360  init_suffix + "." + str(myindex) + ".rmf3",
361  self.crosslink_restraints)
362  output.write_rmf(init_suffix + "." + str(myindex) + ".rmf3")
363  output.close_rmf(init_suffix + "." + str(myindex) + ".rmf3")
364 
365 #----------------------------------------------
366 
367  if not self.test_mode:
368  print("Setting up production rmf files")
369  rmfname = rmf_dir + "/" + str(myindex) + ".rmf3"
370  output.init_rmf(rmfname, output_hierarchies, geometries=self.vars["geometries"])
371 
372  if self.crosslink_restraints:
373  output.add_restraints_to_rmf(rmfname, self.crosslink_restraints)
374 
375  ntimes_at_low_temp = 0
376 
377  if myindex == 0:
378  self.show_info()
379  self.replica_exchange_object.set_was_used(True)
380  for i in range(self.vars["number_of_frames"]):
381  for nr in range(self.vars["num_sample_rounds"]):
382  if sampler_md is not None:
383  sampler_md.optimize(self.vars["molecular_dynamics_steps"])
384  if sampler_mc is not None:
385  sampler_mc.optimize(self.vars["monte_carlo_steps"])
386  score = IMP.pmi.tools.get_restraint_set(self.model).evaluate(False)
387  output.set_output_entry("score", score)
388 
389  my_temp_index = int(rex.get_my_temp() * temp_index_factor)
390 
391  if min_temp_index == my_temp_index:
392  print("--- frame %s score %s " % (str(i), str(score)))
393 
394  if not self.test_mode:
395  if i % self.vars["nframes_write_coordinates"]==0:
396  print('--- writing coordinates')
397  if self.vars["number_of_best_scoring_models"] > 0:
398  output.write_pdb_best_scoring(score)
399  output.write_rmf(rmfname)
400  output.set_output_entry("rmf_file", rmfname)
401  output.set_output_entry("rmf_frame_index", ntimes_at_low_temp)
402  else:
403  output.set_output_entry("rmf_file", rmfname)
404  output.set_output_entry("rmf_frame_index", '-1')
405  output.write_stat2(low_temp_stat_file)
406  ntimes_at_low_temp += 1
407 
408  if not self.test_mode:
409  output.write_stat2(replica_stat_file)
410  rex.swap_temp(i, score)
411 
412 # ----------------------------------------------------------------------
413 class BuildSystem(object):
414  """A macro to build a IMP::pmi::topology::System based on a TopologyReader object.
415  Easily create multi-state systems by calling this macro
416  repeatedly with different TopologyReader objects!
417  A useful function is get_molecules() which returns the PMI Molecules grouped by state
418  as a dictionary with key = (molecule name), value = IMP.pmi.topology.Molecule
419 
420  Quick multi-state system:
421  @code{.python}
422  mdl = IMP.Model()
423  reader1 = IMP.pmi.topology.TopologyReader(tfile1)
424  reader2 = IMP.pmi.topology.TopologyReader(tfile2)
425 
426  bs = IMP.pmi.macros.BuildSystem(mdl)
427  bs.add_state(reader1)
428  bs.add_state(reader2)
429  bs.execute_macro() # build everything including degrees of freedom
430 
431  IMP.atom.show_molecular_hierarchy(bs.get_hierarchy())
432  ### now you have a two state system, you add restraints etc
433  @endcode
434 
435  \note The "domain name" entry of the topology reader is not used.
436  All molecules are set up by the component name, but split into rigid bodies
437  as requested.
438  """
439  def __init__(self,
440  mdl,
441  sequence_connectivity_scale=4.0,
442  force_create_gmm_files=False):
443  """Constructor
444  @param mdl An IMP Model
445  @param sequence_connectivity_scale For scaling the connectivity restraint
446  @param force_create_gmm_files If True, will sample and create GMMs
447  no matter what. If False, will only only sample if the
448  files don't exist. If number of Gaussians is zero, won't
449  do anything.
450  """
451  self.mdl = mdl
452  self.system = IMP.pmi.topology.System(self.mdl)
453  self._readers = [] #internal storage of the TopologyReaders (one per state)
454  self._domains = [] #internal storage of key=domain name, value=(atomic_res,non_atomic_res). (one per state)
455  self.force_create_gmm_files = force_create_gmm_files
456  def add_state(self,reader):
457  """Add a state using the topology info in a IMP::pmi::topology::TopologyReader object.
458  When you are done adding states, call execute_macro()
459  @param reader The TopologyReader object
460  """
461  state = self.system.create_state()
462  self._readers.append(reader)
463  these_domains = {}
464 
465  ### setup representation
466  # setup each component once, then each domain gets seperate structure, representation
467  for molname in reader.get_unique_molecules():
468  mlist = reader.get_unique_molecules()[molname]
469  seq = IMP.pmi.topology.Sequences(mlist[0].fasta_file)
470  mol = state.create_molecule(molname,seq[mlist[0].fasta_id],mlist[0].chain)
471  for domain in mlist:
472  # we build everything in the residue range, even if it extends beyond what's in the actual PDB file
473  if domain.residue_range==[] or domain.residue_range is None:
474  domain_res = mol.get_residues()
475  else:
476  start = domain.residue_range[0]+domain.pdb_offset
477  if domain.residue_range[1]==-1:
478  end = -1
479  else:
480  end = domain.residue_range[1]+domain.pdb_offset
481  domain_res = mol.residue_range(start,end)
482  if domain.pdb_file=="BEADS":
483  mol.add_representation(domain_res,
484  resolutions=[domain.bead_size],
485  setup_particles_as_densities=(domain.em_residues_per_gaussian!=0))
486  these_domains[domain.domain_name] = (set(),domain_res)
487  elif domain.pdb_file=="IDEAL_HELIX":
488  mol.add_representation(domain_res,
489  resolutions=reader.resolutions,
490  ideal_helix=True,
491  density_residues_per_component=domain.em_residues_per_gaussian,
492  density_prefix=domain.density_prefix,
493  density_force_compute=self.force_create_gmm_files)
494  these_domains[domain.domain_name] = (domain_res,set())
495  else:
496  domain_atomic = mol.add_structure(domain.pdb_file,
497  domain.chain,
498  domain.residue_range,
499  domain.pdb_offset,
500  soft_check=True)
501  domain_non_atomic = domain_res - domain_atomic
502  if domain.em_residues_per_gaussian==0:
503  mol.add_representation(domain_atomic,
504  resolutions=reader.resolutions)
505  if len(domain_non_atomic)>0:
506  mol.add_representation(domain_non_atomic,
507  resolutions=[domain.bead_size])
508  else:
509  mol.add_representation(domain_atomic,
510  resolutions=reader.resolutions,
511  density_residues_per_component=domain.em_residues_per_gaussian,
512  density_prefix=domain.density_prefix,
513  density_force_compute=self.force_create_gmm_files)
514  if len(domain_non_atomic)>0:
515  mol.add_representation(domain_non_atomic,
516  resolutions=[domain.bead_size],
517  setup_particles_as_densities=True)
518  these_domains[domain.domain_name] = (domain_atomic,domain_non_atomic)
519  self._domains.append(these_domains)
520  print('State',len(self.system.states),'added')
521 
522  def get_molecules(self):
523  """Return list of all molecules grouped by state.
524  For each state, it's a dictionary of Molecules where key is the molecule name
525  """
526  return [s.get_molecules() for s in self.system.get_states()]
527 
528  def execute_macro(self):
529  """Builds representations and sets up degrees of freedom"""
530  print("BuildSystem: building representations")
531  self.root_hier = self.system.build()
532 
533  print("BuildSystem: setting up degrees of freedom")
534  self.dof = IMP.pmi.dof.DegreesOfFreedom(self.mdl)
535  for nstate,reader in enumerate(self._readers):
536  rbs = reader.get_rigid_bodies()
537  srbs = reader.get_super_rigid_bodies()
538  csrbs = reader.get_chains_of_super_rigid_bodies()
539 
540  # add rigid bodies
541  domains_in_rbs = set()
542  for rblist in rbs:
543  all_res = IMP.pmi.tools.OrderedSet()
544  bead_res = IMP.pmi.tools.OrderedSet()
545  for domain in rblist:
546  all_res|=self._domains[nstate][domain][0]
547  bead_res|=self._domains[nstate][domain][1]
548  domains_in_rbs.add(domain)
549  all_res|=bead_res
550  self.dof.create_rigid_body(all_res,
551  nonrigid_parts=bead_res)
552 
553  # if you have any BEAD domains not in an RB, set them as flexible beads
554  for molname in reader.get_unique_molecules():
555  mlist = reader.get_unique_molecules()[molname]
556  for domain in mlist:
557  if domain.pdb_file=="BEADS" and domain not in domains_in_rbs:
558  self.dof.create_flexible_beads(self._domains[nstate][domain.domain_name][1])
559 
560  # add super rigid bodies
561  for srblist in srbs:
562  all_res = IMP.pmi.tools.OrderedSet()
563  for domain in srblist:
564  all_res|=self._domains[nstate][domain][0]
565  self.dof.create_super_rigid_body(all_res)
566 
567  # add chains of super rigid bodies
568  for csrblist in csrbs:
569  all_res = IMP.pmi.tools.OrderedSet()
570  for domain in csrblist:
571  all_res|=self._domains[nstate][domain][0]
572  all_res = list(all_res)
573  all_res.sort(key=lambda r:r.get_index())
574  self.dof.create_super_rigid_body(all_res,chain_min_length=2,chain_max_length=3)
575  return self.root_hier,self.dof
576 
577 @IMP.deprecated_object("2.5", "Use BuildSystem instead")
578 class BuildModel(object):
579  """A macro to build a Representation based on a Topology and lists of movers
580  DEPRECATED - Use BuildSystem instead.
581  """
582  def __init__(self,
583  model,
584  component_topologies,
585  list_of_rigid_bodies=[],
586  list_of_super_rigid_bodies=[],
587  chain_of_super_rigid_bodies=[],
588  sequence_connectivity_scale=4.0,
589  add_each_domain_as_rigid_body=False,
590  force_create_gmm_files=False):
591  """Constructor.
592  @param model The IMP model
593  @param component_topologies List of
594  IMP.pmi.topology.ComponentTopology items
595  @param list_of_rigid_bodies List of lists of domain names that will
596  be moved as rigid bodies.
597  @param list_of_super_rigid_bodies List of lists of domain names
598  that will move together in an additional Monte Carlo move.
599  @param chain_of_super_rigid_bodies List of lists of domain names
600  (choices can only be from the same molecule). Each of these
601  groups will be moved rigidly. This helps to sample more
602  efficiently complex topologies, made of several rigid bodies,
603  connected by flexible linkers.
604  @param sequence_connectivity_scale For scaling the connectivity
605  restraint
606  @param add_each_domain_as_rigid_body That way you don't have to
607  put all of them in the list
608  @param force_create_gmm_files If True, will sample and create GMMs
609  no matter what. If False, will only only sample if the
610  files don't exist. If number of Gaussians is zero, won't
611  do anything.
612  """
613  self.m = model
614  self.simo = IMP.pmi.representation.Representation(self.m,
615  upperharmonic=True,
616  disorderedlength=False)
617 
618  data=component_topologies
619  if list_of_rigid_bodies==[]:
620  print("WARNING: No list of rigid bodies inputted to build_model()")
621  if list_of_super_rigid_bodies==[]:
622  print("WARNING: No list of super rigid bodies inputted to build_model()")
623  if chain_of_super_rigid_bodies==[]:
624  print("WARNING: No chain of super rigid bodies inputted to build_model()")
625  all_dnames = set([d for sublist in list_of_rigid_bodies+list_of_super_rigid_bodies\
626  +chain_of_super_rigid_bodies for d in sublist])
627  all_available = set([c.domain_name for c in component_topologies])
628  if not all_dnames <= all_available:
629  raise ValueError("All requested movers must reference domain "
630  "names in the component topologies")
631 
632  self.domain_dict={}
633  self.resdensities={}
634  super_rigid_bodies={}
635  chain_super_rigid_bodies={}
636  rigid_bodies={}
637 
638  for c in data:
639  comp_name = c.name
640  hier_name = c.domain_name
641  color = c.color
642  fasta_file = c.fasta_file
643  fasta_id = c.fasta_id
644  pdb_name = c.pdb_file
645  chain_id = c.chain
646  res_range = c.residue_range
647  offset = c.pdb_offset
648  bead_size = c.bead_size
649  em_num_components = c.em_residues_per_gaussian
650  em_txt_file_name = c.gmm_file
651  em_mrc_file_name = c.mrc_file
652 
653  if comp_name not in self.simo.get_component_names():
654  self.simo.create_component(comp_name,color=0.0)
655  self.simo.add_component_sequence(comp_name,fasta_file,fasta_id)
656 
657  # create hierarchy (adds resolutions, beads) with autobuild and optionally add EM data
658  if em_num_components==0:
659  read_em_files=False
660  include_res0=False
661  else:
662  if (not os.path.isfile(em_txt_file_name)) or force_create_gmm_files:
663  read_em_files=False
664  include_res0=True
665  else:
666  read_em_files=True
667  include_res0=False
668 
669  outhier=self.autobuild(self.simo,comp_name,pdb_name,chain_id,
670  res_range,include_res0,beadsize=bead_size,
671  color=color,offset=offset)
672  if em_num_components!=0:
673  if read_em_files:
674  print("will read GMM files")
675  else:
676  print("will calculate GMMs")
677 
678  dens_hier,beads=self.create_density(self.simo,comp_name,outhier,em_txt_file_name,
679  em_mrc_file_name,em_num_components,read_em_files)
680  self.simo.add_all_atom_densities(comp_name, particles=beads)
681  dens_hier+=beads
682  else:
683  dens_hier=[]
684 
685  self.resdensities[hier_name]=dens_hier
686  self.domain_dict[hier_name]=outhier+dens_hier
687 
688  # setup basic restraints
689  for c in self.simo.get_component_names():
690  self.simo.setup_component_sequence_connectivity(c,scale=sequence_connectivity_scale)
691  self.simo.setup_component_geometry(c)
692 
693  # create movers
694  for rblist in list_of_rigid_bodies:
695  rb=[]
696  for rbmember in rblist:
697  rb+=[h for h in self.domain_dict[rbmember]]
698  self.simo.set_rigid_body_from_hierarchies(rb)
699  for srblist in list_of_super_rigid_bodies:
700  srb=[]
701  for srbmember in rblist:
702  srb+=[h for h in self.domain_dict[srbmember]]
703  self.simo.set_super_rigid_body_from_hierarchies(srb)
704  for clist in chain_of_super_rigid_bodies:
705  crb=[]
706  for crbmember in rblist:
707  crb+=[h for h in self.domain_dict[crbmember]]
708  self.simo.set_chain_of_super_rigid_bodies(crb,2,3)
709 
710  self.simo.set_floppy_bodies()
711  self.simo.setup_bonds()
712 
714  '''Return the Representation object'''
715  return self.simo
716 
717  def get_density_hierarchies(self,hier_name_list):
718  # return a list of density hierarchies
719  # specify the list of hierarchy names
720  dens_hier_list=[]
721  for hn in hier_name_list:
722  print(hn)
723  dens_hier_list+=self.resdensities[hn]
724  return dens_hier_list
725 
726  def set_gmm_models_directory(self,directory_name):
727  self.gmm_models_directory=directory_name
728 
729  def get_pdb_bead_bits(self,hierarchy):
730  pdbbits=[]
731  beadbits=[]
732  helixbits=[]
733  for h in hierarchy:
734  if "_pdb" in h.get_name():pdbbits.append(h)
735  if "_bead" in h.get_name():beadbits.append(h)
736  if "_helix" in h.get_name():helixbits.append(h)
737  return (pdbbits,beadbits,helixbits)
738 
739  def scale_bead_radii(self,nresidues,scale):
740  scaled_beads=set()
741  for h in self.domain_dict:
742  (pdbbits,beadbits,helixbits)=self.get_pdb_bead_bits(self.domain_dict[h])
743  slope=(1.0-scale)/(1.0-float(nresidues))
744 
745  for b in beadbits:
746  # I have to do the following
747  # because otherwise we'll scale more than once
748  if b not in scaled_beads:
749  scaled_beads.add(b)
750  else:
751  continue
752  radius=IMP.core.XYZR(b).get_radius()
753  num_residues=len(IMP.pmi.tools.get_residue_indexes(b))
754  scale_factor=slope*float(num_residues)+1.0
755  print(scale_factor)
756  new_radius=scale_factor*radius
757  IMP.core.XYZR(b).set_radius(new_radius)
758  print(b.get_name())
759  print("particle with radius "+str(radius)+" and "+str(num_residues)+" residues scaled to a new radius "+str(new_radius))
760 
761 
762  def create_density(self,simo,compname,comphier,txtfilename,mrcfilename,num_components,read=True):
763  #density generation for the EM restraint
764  (pdbbits,beadbits,helixbits)=self.get_pdb_bead_bits(comphier)
765  #get the number of residues from the pdb bits
766  res_ind=[]
767  for pb in pdbbits+helixbits:
768  for p in IMP.core.get_leaves(pb):
770 
771  number_of_residues=len(set(res_ind))
772  outhier=[]
773  if read:
774  if len(pdbbits)!=0:
775  outhier+=simo.add_component_density(compname,
776  pdbbits,
777  num_components=num_components,
778  resolution=0,
779  inputfile=txtfilename)
780  if len(helixbits)!=0:
781  outhier+=simo.add_component_density(compname,
782  helixbits,
783  num_components=num_components,
784  resolution=1,
785  inputfile=txtfilename)
786 
787 
788  else:
789  if len(pdbbits)!=0:
790  num_components=number_of_residues//abs(num_components)+1
791  outhier+=simo.add_component_density(compname,
792  pdbbits,
793  num_components=num_components,
794  resolution=0,
795  outputfile=txtfilename,
796  outputmap=mrcfilename,
797  multiply_by_total_mass=True)
798 
799  if len(helixbits)!=0:
800  num_components=number_of_residues//abs(num_components)+1
801  outhier+=simo.add_component_density(compname,
802  helixbits,
803  num_components=num_components,
804  resolution=1,
805  outputfile=txtfilename,
806  outputmap=mrcfilename,
807  multiply_by_total_mass=True)
808 
809  return outhier,beadbits
810 
811  def autobuild(self,simo,comname,pdbname,chain,resrange,include_res0=False,
812  beadsize=5,color=0.0,offset=0):
813  if pdbname is not None and pdbname is not "IDEAL_HELIX" and pdbname is not "BEADS" :
814  if include_res0:
815  outhier=simo.autobuild_model(comname,
816  pdbname=pdbname,
817  chain=chain,
818  resrange=resrange,
819  resolutions=[0,1,10],
820  offset=offset,
821  color=color,
822  missingbeadsize=beadsize)
823  else:
824  outhier=simo.autobuild_model(comname,
825  pdbname=pdbname,
826  chain=chain,
827  resrange=resrange,
828  resolutions=[1,10],
829  offset=offset,
830  color=color,
831  missingbeadsize=beadsize)
832 
833 
834  elif pdbname is not None and pdbname is "IDEAL_HELIX" and pdbname is not "BEADS" :
835  outhier=simo.add_component_ideal_helix(comname,
836  resolutions=[1,10],
837  resrange=resrange,
838  color=color,
839  show=False)
840 
841  elif pdbname is not None and pdbname is not "IDEAL_HELIX" and pdbname is "BEADS" :
842  outhier=simo.add_component_necklace(comname,resrange[0],resrange[1],beadsize,color=color)
843 
844  else:
845 
846  seq_len=len(simo.sequence_dict[comname])
847  outhier=simo.add_component_necklace(comname,
848  begin=1,
849  end=seq_len,
850  length=beadsize)
851 
852  return outhier
853 
854 
855 @IMP.deprecated_object("2.5", "Use BuildSystem instead")
856 class BuildModel1(object):
857  """Deprecated building macro - use BuildSystem()"""
858 
859  def __init__(self, representation):
860  """Constructor.
861  @param representation The PMI representation
862  """
863  self.simo=representation
864  self.gmm_models_directory="."
865  self.rmf_file={}
866  self.rmf_frame_number={}
867  self.rmf_names_map={}
868 
869  def set_gmm_models_directory(self,directory_name):
870  self.gmm_models_directory=directory_name
871 
872  def build_model(self,data_structure,sequence_connectivity_scale=4.0,
873  sequence_connectivity_resolution=10,rmf_file=None,rmf_frame_number=0,rmf_file_map=None):
874  """Create model.
875  @param data_structure List of lists containing these entries:
876  comp_name, hier_name, color, fasta_file, fasta_id, pdb_name, chain_id,
877  res_range, read_em_files, bead_size, rb, super_rb,
878  em_num_components, em_txt_file_name, em_mrc_file_name
879  @param sequence_connectivity_scale
880  @param rmf_file
881  @param rmf_frame_number
882  @param rmf_file_map : a dictionary that map key=component_name:value=(rmf_file_name,
883  rmf_frame_number,
884  rmf_component_name)
885  """
886  self.domain_dict={}
887  self.resdensities={}
888  super_rigid_bodies={}
889  chain_super_rigid_bodies={}
890  rigid_bodies={}
891 
892  for d in data_structure:
893  comp_name = d[0]
894  hier_name = d[1]
895  color = d[2]
896  fasta_file = d[3]
897  fasta_id = d[4]
898  pdb_name = d[5]
899  chain_id = d[6]
900  res_range = d[7][0:2]
901  try:
902  offset = d[7][2]
903  except:
904  offset = 0
905  read_em_files = d[8]
906  bead_size = d[9]
907  rb = d[10]
908  super_rb = d[11]
909  em_num_components = d[12]
910  em_txt_file_name = d[13]
911  em_mrc_file_name = d[14]
912  chain_of_super_rb = d[15]
913 
914  if comp_name not in self.simo.get_component_names():
915  self.simo.create_component(comp_name,color=0.0)
916  self.simo.add_component_sequence(comp_name,fasta_file,fasta_id)
917  outhier=self.autobuild(self.simo,comp_name,pdb_name,chain_id,res_range,read=read_em_files,beadsize=bead_size,color=color,offset=offset)
918 
919 
920  if not read_em_files is None:
921  if em_txt_file_name is " ": em_txt_file_name=self.gmm_models_directory+"/"+hier_name+".txt"
922  if em_mrc_file_name is " ": em_mrc_file_name=self.gmm_models_directory+"/"+hier_name+".mrc"
923 
924 
925  dens_hier,beads=self.create_density(self.simo,comp_name,outhier,em_txt_file_name,em_mrc_file_name,em_num_components,read_em_files)
926 
927  self.simo.add_all_atom_densities(comp_name, particles=beads)
928  dens_hier+=beads
929 
930  else:
931  dens_hier=[]
932 
933  self.resdensities[hier_name]=dens_hier
934  self.domain_dict[hier_name]=outhier+dens_hier
935 
936  if rb is not None:
937  if rb not in rigid_bodies:
938  rigid_bodies[rb]=[h for h in self.domain_dict[hier_name]]
939  else:
940  rigid_bodies[rb]+=[h for h in self.domain_dict[hier_name]]
941 
942 
943  if super_rb is not None:
944  for k in super_rb:
945  if k not in super_rigid_bodies:
946  super_rigid_bodies[k]=[h for h in self.domain_dict[hier_name]]
947  else:
948  super_rigid_bodies[k]+=[h for h in self.domain_dict[hier_name]]
949 
950  if chain_of_super_rb is not None:
951  for k in chain_of_super_rb:
952  if k not in chain_super_rigid_bodies:
953  chain_super_rigid_bodies[k]=[h for h in self.domain_dict[hier_name]]
954  else:
955  chain_super_rigid_bodies[k]+=[h for h in self.domain_dict[hier_name]]
956 
957 
958 
959  self.rigid_bodies=rigid_bodies
960 
961  for c in self.simo.get_component_names():
962  if rmf_file is not None:
963  rf=rmf_file
964  rfn=rmf_frame_number
965  self.simo.set_coordinates_from_rmf(c, rf,rfn)
966  elif rmf_file_map:
967  for k in rmf_file_map:
968  cname=k
969  rf=rmf_file_map[k][0]
970  rfn=rmf_file_map[k][1]
971  rcname=rmf_file_map[k][2]
972  self.simo.set_coordinates_from_rmf(cname, rf,rfn,rcname)
973  else:
974  if c in self.rmf_file:
975  rf=self.rmf_file[c]
976  rfn=self.rmf_frame_number[c]
977  rfm=self.rmf_names_map[c]
978  self.simo.set_coordinates_from_rmf(c, rf,rfn,representation_name_to_rmf_name_map=rfm)
979  self.simo.setup_component_sequence_connectivity(c,resolution=sequence_connectivity_resolution,scale=sequence_connectivity_scale)
980  self.simo.setup_component_geometry(c)
981 
982  for rb in rigid_bodies:
983  self.simo.set_rigid_body_from_hierarchies(rigid_bodies[rb])
984 
985  for k in super_rigid_bodies:
986  self.simo.set_super_rigid_body_from_hierarchies(super_rigid_bodies[k])
987 
988  for k in chain_super_rigid_bodies:
989  self.simo.set_chain_of_super_rigid_bodies(chain_super_rigid_bodies[k],2,3)
990 
991  self.simo.set_floppy_bodies()
992  self.simo.setup_bonds()
993 
994 
995 
996  def set_rmf_file(self,component_name,rmf_file,rmf_frame_number,rmf_names_map=None):
997  self.rmf_file[component_name]=rmf_file
998  self.rmf_frame_number[component_name]=rmf_frame_number
999  self.rmf_names_map[component_name]=rmf_names_map
1000 
1001  def get_density_hierarchies(self,hier_name_list):
1002  # return a list of density hierarchies
1003  # specify the list of hierarchy names
1004  dens_hier_list=[]
1005  for hn in hier_name_list:
1006  print(hn)
1007  dens_hier_list+=self.resdensities[hn]
1008  return dens_hier_list
1009 
1010  def get_pdb_bead_bits(self,hierarchy):
1011  pdbbits=[]
1012  beadbits=[]
1013  helixbits=[]
1014  for h in hierarchy:
1015  if "_pdb" in h.get_name():pdbbits.append(h)
1016  if "_bead" in h.get_name():beadbits.append(h)
1017  if "_helix" in h.get_name():helixbits.append(h)
1018  return (pdbbits,beadbits,helixbits)
1019 
1020  def scale_bead_radii(self,nresidues,scale):
1021  scaled_beads=set()
1022  for h in self.domain_dict:
1023  (pdbbits,beadbits,helixbits)=self.get_pdb_bead_bits(self.domain_dict[h])
1024  slope=(1.0-scale)/(1.0-float(nresidues))
1025 
1026  for b in beadbits:
1027  # I have to do the following
1028  # because otherwise we'll scale more than once
1029  if b not in scaled_beads:
1030  scaled_beads.add(b)
1031  else:
1032  continue
1033  radius=IMP.core.XYZR(b).get_radius()
1034  num_residues=len(IMP.pmi.tools.get_residue_indexes(b))
1035  scale_factor=slope*float(num_residues)+1.0
1036  print(scale_factor)
1037  new_radius=scale_factor*radius
1038  IMP.core.XYZR(b).set_radius(new_radius)
1039  print(b.get_name())
1040  print("particle with radius "+str(radius)+" and "+str(num_residues)+" residues scaled to a new radius "+str(new_radius))
1041 
1042 
1043  def create_density(self,simo,compname,comphier,txtfilename,mrcfilename,num_components,read=True):
1044  #density generation for the EM restraint
1045  (pdbbits,beadbits,helixbits)=self.get_pdb_bead_bits(comphier)
1046 
1047  #get the number of residues from the pdb bits
1048  res_ind=[]
1049  for pb in pdbbits+helixbits:
1050  for p in IMP.core.get_leaves(pb):
1052 
1053  number_of_residues=len(set(res_ind))
1054  outhier=[]
1055  if read:
1056  if len(pdbbits)!=0:
1057  outhier+=simo.add_component_density(compname,
1058  pdbbits,
1059  num_components=num_components, # number of gaussian into which the simulated density is approximated
1060  resolution=0, # resolution that you want to calculate the simulated density
1061  inputfile=txtfilename) # read what it was calculated before
1062  if len(helixbits)!=0:
1063  outhier+=simo.add_component_density(compname,
1064  helixbits,
1065  num_components=num_components, # number of gaussian into which the simulated density is approximated
1066  resolution=1, # resolution that you want to calculate the simulated density
1067  inputfile=txtfilename) # read what it was calculated before
1068 
1069 
1070  else:
1071  if len(pdbbits)!=0:
1072  if num_components<0:
1073  #if negative calculate the number of gmm components automatically
1074  # from the number of residues
1075  num_components=number_of_residues/abs(num_components)
1076  outhier+=simo.add_component_density(compname,
1077  pdbbits,
1078  num_components=num_components, # number of gaussian into which the simulated density is approximated
1079  resolution=0, # resolution that you want to calculate the simulated density
1080  outputfile=txtfilename, # do the calculation
1081  outputmap=mrcfilename,
1082  multiply_by_total_mass=True) # do the calculation and output the mrc
1083 
1084  if len(helixbits)!=0:
1085  if num_components<0:
1086  #if negative calculate the number of gmm components automatically
1087  # from the number of residues
1088  num_components=number_of_residues/abs(num_components)
1089  outhier+=simo.add_component_density(compname,
1090  helixbits,
1091  num_components=num_components, # number of gaussian into which the simulated density is approximated
1092  resolution=1, # resolution that you want to calculate the simulated density
1093  outputfile=txtfilename, # do the calculation
1094  outputmap=mrcfilename,
1095  multiply_by_total_mass=True) # do the calculation and output the mrc
1096 
1097  return outhier,beadbits
1098 
1099  def autobuild(self,simo,comname,pdbname,chain,resrange,read=True,beadsize=5,color=0.0,offset=0):
1100 
1101  if pdbname is not None and pdbname is not "IDEAL_HELIX" and pdbname is not "BEADS" and pdbname is not "DENSITY" :
1102  if resrange[-1]==-1: resrange=(resrange[0],len(simo.sequence_dict[comname]))
1103  if read==False:
1104  outhier=simo.autobuild_model(comname,
1105  pdbname=pdbname,
1106  chain=chain,
1107  resrange=resrange,
1108  resolutions=[0,1,10],
1109  offset=offset,
1110  color=color,
1111  missingbeadsize=beadsize)
1112  else:
1113  outhier=simo.autobuild_model(comname,
1114  pdbname=pdbname,
1115  chain=chain,
1116  resrange=resrange,
1117  resolutions=[1,10],
1118  offset=offset,
1119  color=color,
1120  missingbeadsize=beadsize)
1121 
1122  elif pdbname is not None and pdbname is "IDEAL_HELIX" and pdbname is not "BEADS" and pdbname is not "DENSITY" :
1123  outhier=simo.add_component_ideal_helix(comname,
1124  resolutions=[1,10],
1125  resrange=resrange,
1126  color=color,
1127  show=False)
1128 
1129  elif pdbname is not None and pdbname is not "IDEAL_HELIX" and pdbname is "BEADS" and pdbname is not "DENSITY" :
1130  seq_len=resrange[1]
1131  if resrange[1]==-1:
1132  seq_len=len(simo.sequence_dict[comname])
1133  outhier=simo.add_component_necklace(comname,resrange[0],seq_len,beadsize,color=color)
1134 
1135  elif pdbname is not None and pdbname is not "IDEAL_HELIX" and pdbname is not "BEADS" and pdbname is "DENSITY" :
1136  outhier=[]
1137 
1138  else:
1139 
1140  seq_len=len(simo.sequence_dict[comname])
1141  outhier=simo.add_component_necklace(comname,
1142  begin=1,
1143  end=seq_len,
1144  length=beadsize)
1145 
1146  return outhier
1147 
1148  def set_coordinates(self,hier_name,xyz_tuple):
1149  hier=self.domain_dict[hier_name]
1150  for h in IMP.atom.get_leaves(hier):
1151  p=h.get_particle()
1153  pass
1154  else:
1155  IMP.core.XYZ(p).set_coordinates(xyz_tuple)
1156 
1157  def save_rmf(self,rmfname):
1158 
1160  o.init_rmf(rmfname,[self.simo.prot])
1161  o.write_rmf(rmfname)
1162  o.close_rmf(rmfname)
1163 # -----------------------------------------------------------------------
1164 
1165 @IMP.deprecated_object("2.5", "Use BuildSystem instead")
1166 def BuildModel0(
1167  m,
1168  data,
1169  resolutions=[1,
1170  10],
1171  missing_bead_size=20,
1172  residue_per_gaussian=None):
1173  '''
1174  Construct a component for each subunit (no splitting, nothing fancy).
1175 
1176  You can pass the resolutions and the bead size for the missing residue regions.
1177  To use this macro, you must provide the following data structure:
1178 
1179  Component pdbfile chainid rgb color fastafile sequence id
1180  in fastafile
1181 
1182 data = [("Rpb1", pdbfile, "A", 0.00000000, (fastafile, 0)),
1183  ("Rpb2", pdbfile, "B", 0.09090909, (fastafile, 1)),
1184  ("Rpb3", pdbfile, "C", 0.18181818, (fastafile, 2)),
1185  ("Rpb4", pdbfile, "D", 0.27272727, (fastafile, 3)),
1186  ("Rpb5", pdbfile, "E", 0.36363636, (fastafile, 4)),
1187  ("Rpb6", pdbfile, "F", 0.45454545, (fastafile, 5)),
1188  ("Rpb7", pdbfile, "G", 0.54545455, (fastafile, 6)),
1189  ("Rpb8", pdbfile, "H", 0.63636364, (fastafile, 7)),
1190  ("Rpb9", pdbfile, "I", 0.72727273, (fastafile, 8)),
1191  ("Rpb10", pdbfile, "L", 0.81818182, (fastafile, 9)),
1192  ("Rpb11", pdbfile, "J", 0.90909091, (fastafile, 10)),
1193  ("Rpb12", pdbfile, "K", 1.00000000, (fastafile, 11))]
1194 
1195  '''
1196 
1198 
1199  # the dictionary for the hierarchies,
1200  hierarchies = {}
1201 
1202  for d in data:
1203  # retrieve the information from the data structure
1204  component_name = d[0]
1205  pdb_file = d[1]
1206  chain_id = d[2]
1207  color_id = d[3]
1208  fasta_file = d[4][0]
1209  # this function
1210  fastids = IMP.pmi.tools.get_ids_from_fasta_file(fasta_file)
1211  fasta_file_id = d[4][1]
1212  # avoid to add a component with the same name
1213  r.create_component(component_name,
1214  color=color_id)
1215 
1216  r.add_component_sequence(component_name,
1217  fasta_file,
1218  id=fastids[fasta_file_id])
1219 
1220  hierarchies = r.autobuild_model(component_name,
1221  pdb_file,
1222  chain_id,
1223  resolutions=resolutions,
1224  missingbeadsize=missing_bead_size)
1225 
1226  r.show_component_table(component_name)
1227 
1228  r.set_rigid_bodies([component_name])
1229 
1230  r.set_chain_of_super_rigid_bodies(
1231  hierarchies,
1232  min_length=2,
1233  max_length=2)
1234 
1235  r.setup_component_sequence_connectivity(component_name, resolution=1)
1236  r.setup_component_geometry(component_name)
1237 
1238  r.setup_bonds()
1239  # put it at the end of rigid bodies
1240  r.set_floppy_bodies()
1241 
1242  # set current coordinates as reference for RMSD calculation
1243  r.set_current_coordinates_as_reference_for_rmsd("Reference")
1244 
1245  return r
1246 
1247 # ----------------------------------------------------------------------
1248 
1250  """A macro for running all the basic operations of analysis.
1251  Includes clustering, precision analysis, and making ensemble density maps.
1252  A number of plots are also supported.
1253  """
1254  def __init__(self, model,
1255  merge_directories=["./"],
1256  stat_file_name_suffix="stat",
1257  best_pdb_name_suffix="model",
1258  do_clean_first=True,
1259  do_create_directories=True,
1260  global_output_directory="output/",
1261  replica_stat_file_suffix="stat_replica",
1262  global_analysis_result_directory="./analysis/"):
1263  """Constructor.
1264  @param model The IMP model
1265  @param stat_file_name_suffix
1266  @param merge_directories The directories containing output files
1267  @param best_pdb_name_suffix
1268  @param do_clean_first
1269  @param do_create_directories
1270  @param global_output_directory Where everything is
1271  @param replica_stat_file_suffix
1272  @param global_analysis_result_directory
1273  """
1274 
1275  try:
1276  from mpi4py import MPI
1277  self.comm = MPI.COMM_WORLD
1278  self.rank = self.comm.Get_rank()
1279  self.number_of_processes = self.comm.size
1280  except ImportError:
1281  self.rank = 0
1282  self.number_of_processes = 1
1283 
1284  self.model = model
1285  stat_dir = global_output_directory
1286  self.stat_files = []
1287  # it contains the position of the root directories
1288  for rd in merge_directories:
1289  stat_files = glob.glob(os.path.join(rd,stat_dir,"stat.*.out"))
1290  if len(stat_files)==0:
1291  print("WARNING: no stat files found in",os.path.join(rd,stat_dir))
1292  self.stat_files += stat_files
1293 
1294 
1296  score_key="SimplifiedModel_Total_Score_None",
1297  rmf_file_key="rmf_file",
1298  rmf_file_frame_key="rmf_frame_index",
1299  outputdir="./",
1300  get_every=1,
1301  nframes_trajectory=10000):
1302  """ Get a trajectory of the modeling run, for generating demonstrative movies
1303  @param score_key The score for ranking models
1304  @param rmf_file_key Key pointing to RMF filename
1305  @param rmf_file_frame_key Key pointing to RMF frame number
1306  @param outputdir The local output directory used in the run
1307  @param get_every Extract every nth frame
1308  @param nframes_trajectory Total number of frames of the trajectory
1309  """
1310  from operator import itemgetter
1311  import math
1312 
1313  trajectory_models = IMP.pmi.io.get_trajectory_models(self.stat_files,
1314  score_key,
1315  rmf_file_key,
1316  rmf_file_frame_key,
1317  get_every)
1318  rmf_file_list=trajectory_models[0]
1319  rmf_file_frame_list=trajectory_models[1]
1320  score_list=list(map(float, trajectory_models[2]))
1321 
1322  max_score=max(score_list)
1323  min_score=min(score_list)
1324 
1325 
1326  bins=[(max_score-min_score)*math.exp(-float(i))+min_score for i in range(nframes_trajectory)]
1327  binned_scores=[None]*nframes_trajectory
1328  binned_model_indexes=[-1]*nframes_trajectory
1329 
1330  for model_index,s in enumerate(score_list):
1331  bins_score_diffs=[abs(s-b) for b in bins]
1332  bin_index=min(enumerate(bins_score_diffs), key=itemgetter(1))[0]
1333  if binned_scores[bin_index]==None:
1334  binned_scores[bin_index]=s
1335  binned_model_indexes[bin_index]=model_index
1336  else:
1337  old_diff=abs(binned_scores[bin_index]-bins[bin_index])
1338  new_diff=abs(s-bins[bin_index])
1339  if new_diff < old_diff:
1340  binned_scores[bin_index]=s
1341  binned_model_indexes[bin_index]=model_index
1342 
1343  print(binned_scores)
1344  print(binned_model_indexes)
1345 
1346 
1347 
1348  def clustering(self,
1349  score_key="SimplifiedModel_Total_Score_None",
1350  rmf_file_key="rmf_file",
1351  rmf_file_frame_key="rmf_frame_index",
1352  state_number=0,
1353  prefiltervalue=None,
1354  feature_keys=[],
1355  outputdir="./",
1356  alignment_components=None,
1357  number_of_best_scoring_models=10,
1358  rmsd_calculation_components=None,
1359  distance_matrix_file='distances.mat',
1360  load_distance_matrix_file=False,
1361  skip_clustering=False,
1362  number_of_clusters=1,
1363  display_plot=False,
1364  exit_after_display=True,
1365  get_every=1,
1366  first_and_last_frames=None,
1367  density_custom_ranges=None,
1368  write_pdb_with_centered_coordinates=False,
1369  voxel_size=5.0):
1370  """ Get the best scoring models, compute a distance matrix, cluster them, and create density maps
1371  @param score_key The score for ranking models
1372  @param rmf_file_key Key pointing to RMF filename
1373  @param rmf_file_frame_key Key pointing to RMF frame number
1374  @param state_number State number to analyze
1375  @param prefiltervalue Only include frames where the score key is below this value
1376  @param feature_keys Keywords for which you want to calculate average,
1377  medians, etc,
1378  @param outputdir The local output directory used in the run
1379  @param alignment_components List of tuples for aligning the structures
1380  e.g. ["Rpb1", (20,100,"Rpb2"), .....]
1381  @param number_of_best_scoring_models Num models to keep per run
1382  @param rmsd_calculation_components List of tuples for calculating RMSD
1383  e.g. ["Rpb1", (20,100,"Rpb2"), .....]
1384  @param distance_matrix_file Where to store/read the distance matrix
1385  @param load_distance_matrix_file Try to load the distance matrix file
1386  @param skip_clustering Just extract the best scoring models and save the pdbs
1387  @param number_of_clusters Number of k-means clusters
1388  @param display_plot Display the distance matrix
1389  @param exit_after_display Exit after displaying distance matrix
1390  @param get_every Extract every nth frame
1391  @param first_and_last_frames A tuple with the first and last frames to be
1392  analyzed. Values are percentages!
1393  Default: get all frames
1394  @param density_custom_ranges List of tuples or strings for density calculation
1395  e.g. ["Rpb1", (20,100,"Rpb2"), .....]
1396  @param write_pdb_with_centered_coordinates
1397  @param voxel_size Used for the density output
1398  """
1399 
1400 
1401  if self.rank==0:
1402  try:
1403  os.mkdir(outputdir)
1404  except:
1405  pass
1406 
1407 
1408  if not load_distance_matrix_file:
1409  if len(self.stat_files)==0: print("ERROR: no stat file found in the given path"); return
1410  my_stat_files=IMP.pmi.tools.chunk_list_into_segments(
1411  self.stat_files,self.number_of_processes)[self.rank]
1412  best_models = IMP.pmi.io.get_best_models(my_stat_files,
1413  score_key,
1414  feature_keys,
1415  rmf_file_key,
1416  rmf_file_frame_key,
1417  prefiltervalue,
1418  get_every)
1419  rmf_file_list=best_models[0]
1420  rmf_file_frame_list=best_models[1]
1421  score_list=best_models[2]
1422  feature_keyword_list_dict=best_models[3]
1423 
1424 # ------------------------------------------------------------------------
1425 # collect all the files and scores
1426 # ------------------------------------------------------------------------
1427 
1428  if self.number_of_processes > 1:
1429  score_list = IMP.pmi.tools.scatter_and_gather(score_list)
1430  rmf_file_list = IMP.pmi.tools.scatter_and_gather(rmf_file_list)
1431  rmf_file_frame_list = IMP.pmi.tools.scatter_and_gather(
1432  rmf_file_frame_list)
1433  for k in feature_keyword_list_dict:
1434  feature_keyword_list_dict[k] = IMP.pmi.tools.scatter_and_gather(
1435  feature_keyword_list_dict[k])
1436 
1437  # sort by score and get the best scoring ones
1438  score_rmf_tuples = list(zip(score_list,
1439  rmf_file_list,
1440  rmf_file_frame_list,
1441  list(range(len(score_list)))))
1442 
1443  if density_custom_ranges:
1444  for k in density_custom_ranges:
1445  if type(density_custom_ranges[k]) is not list:
1446  raise Exception("Density custom ranges: values must be lists of tuples")
1447 
1448  # keep subset of frames if requested
1449  if first_and_last_frames is not None:
1450  nframes = len(score_rmf_tuples)
1451  first_frame = int(first_and_last_frames[0] * nframes)
1452  last_frame = int(first_and_last_frames[1] * nframes)
1453  if last_frame > len(score_rmf_tuples):
1454  last_frame = -1
1455  score_rmf_tuples = score_rmf_tuples[first_frame:last_frame]
1456 
1457  # sort RMFs by the score_key in ascending order, and store the rank
1458  best_score_rmf_tuples = sorted(score_rmf_tuples,
1459  key=lambda x: float(x[0]))[:number_of_best_scoring_models]
1460  best_score_rmf_tuples=[t+(n,) for n,t in enumerate(best_score_rmf_tuples)]
1461 
1462  # sort the feature scores in the same way
1463  best_score_feature_keyword_list_dict = defaultdict(list)
1464  for tpl in best_score_rmf_tuples:
1465  index = tpl[3]
1466  for f in feature_keyword_list_dict:
1467  best_score_feature_keyword_list_dict[f].append(
1468  feature_keyword_list_dict[f][index])
1469 
1470  my_best_score_rmf_tuples = IMP.pmi.tools.chunk_list_into_segments(
1471  best_score_rmf_tuples,
1472  self.number_of_processes)[self.rank]
1473 
1474 
1475 #-------------------------------------------------------------
1476 # read the coordinates
1477 # ------------------------------------------------------------
1478  rmsd_weights = IMP.pmi.io.get_bead_sizes(self.model,
1479  my_best_score_rmf_tuples[0],
1480  rmsd_calculation_components,
1481  state_number=state_number)
1482  got_coords = IMP.pmi.io.read_coordinates_of_rmfs(self.model,
1483  my_best_score_rmf_tuples,
1484  alignment_components,
1485  rmsd_calculation_components,
1486  state_number=state_number)
1487 
1488  # note! the coordinates are simply float tuples, NOT decorators, NOT Vector3D,
1489  # NOR particles, because these object cannot be serialized. We need serialization
1490  # for the parallel computation based on mpi.
1491  all_coordinates=got_coords[0] # dict:key=component name,val=coords per hit
1492  alignment_coordinates=got_coords[1] # same as above, limited to alignment bits
1493  rmsd_coordinates=got_coords[2] # same as above, limited to RMSD bits
1494  rmf_file_name_index_dict=got_coords[3] # dictionary with key=RMF, value=score rank
1495  all_rmf_file_names=got_coords[4] # RMF file per hit
1496 
1497 # ------------------------------------------------------------------------
1498 # optionally don't compute distance matrix or cluster, just write top files
1499 # ------------------------------------------------------------------------
1500  if skip_clustering:
1501  if density_custom_ranges:
1502  DensModule = IMP.pmi.analysis.GetModelDensity(
1503  density_custom_ranges,
1504  voxel=voxel_size)
1505 
1506  dircluster=os.path.join(outputdir,"all_models."+str(n))
1507  try:
1508  os.mkdir(outputdir)
1509  except:
1510  pass
1511  try:
1512  os.mkdir(dircluster)
1513  except:
1514  pass
1515  clusstat=open(os.path.join(dircluster,"stat."+str(self.rank)+".out"),"w")
1516  for cnt,tpl in enumerate(my_best_score_rmf_tuples):
1517  rmf_name=tpl[1]
1518  rmf_frame_number=tpl[2]
1519  tmp_dict={}
1520  index=tpl[4]
1521  for key in best_score_feature_keyword_list_dict:
1522  tmp_dict[key]=best_score_feature_keyword_list_dict[key][index]
1523 
1524  if cnt==0:
1525  prots,rs = IMP.pmi.analysis.get_hiers_and_restraints_from_rmf(
1526  self.model,
1527  rmf_frame_number,
1528  rmf_name)
1529  else:
1530  linking_successful=IMP.pmi.analysis.link_hiers_and_restraints_to_rmf(
1531  self.model,
1532  prots,
1533  rs,
1534  rmf_frame_number,
1535  rmf_name)
1536  if not linking_successful:
1537  continue
1538 
1539  if not prots:
1540  continue
1541 
1542  if IMP.pmi.get_is_canonical(prots[0]):
1543  states = IMP.atom.get_by_type(prots[0],IMP.atom.STATE_TYPE)
1544  prot = states[state_number]
1545  else:
1546  prot = prots[state_number]
1547 
1548 
1549  if cnt==0:
1550  coords_f1=alignment_coordinates[cnt]
1551  if cnt > 0:
1552  coords_f2=alignment_coordinates[cnt]
1553  if coords_f2:
1554  Ali = IMP.pmi.analysis.Alignment(coords_f1, coords_f2)
1555  transformation = Ali.align()[1]
1556  else:
1558 
1559  rbs = set()
1560  for p in IMP.atom.get_leaves(prot):
1561  if not IMP.core.XYZR.get_is_setup(p):
1563  IMP.core.XYZR(p).set_radius(0.0001)
1564  IMP.core.XYZR(p).set_coordinates((0, 0, 0))
1565 
1567  rb = IMP.core.RigidBodyMember(p).get_rigid_body()
1568  rbs.add(rb)
1569  else:
1571  transformation)
1572  for rb in rbs:
1573  IMP.core.transform(rb,transformation)
1574 
1576  out_pdb_fn=os.path.join(dircluster,str(cnt)+"."+str(self.rank)+".pdb")
1577  out_rmf_fn=os.path.join(dircluster,str(cnt)+"."+str(self.rank)+".rmf3")
1578  o.init_pdb(out_pdb_fn,prot)
1579  o.write_pdb(out_pdb_fn,
1580  translate_to_geometric_center=write_pdb_with_centered_coordinates)
1581 
1582  tmp_dict["local_pdb_file_name"]=os.path.basename(out_pdb_fn)
1583  tmp_dict["rmf_file_full_path"]=rmf_name
1584  tmp_dict["local_rmf_file_name"]=os.path.basename(out_rmf_fn)
1585  tmp_dict["local_rmf_frame_number"]=0
1586 
1587  clusstat.write(str(tmp_dict)+"\n")
1588  o.init_rmf(out_rmf_fn,[prot],rs)
1589  o.write_rmf(out_rmf_fn)
1590  o.close_rmf(out_rmf_fn)
1591  # add the density
1592  if density_custom_ranges:
1593  DensModule.add_subunits_density(prot)
1594 
1595  if density_custom_ranges:
1596  DensModule.write_mrc(path=dircluster)
1597  del DensModule
1598  return
1599 
1600 
1601 
1602  # broadcast the coordinates
1603  if self.number_of_processes > 1:
1604  all_coordinates = IMP.pmi.tools.scatter_and_gather(
1605  all_coordinates)
1606  all_rmf_file_names = IMP.pmi.tools.scatter_and_gather(
1607  all_rmf_file_names)
1608  rmf_file_name_index_dict = IMP.pmi.tools.scatter_and_gather(
1609  rmf_file_name_index_dict)
1610  alignment_coordinates=IMP.pmi.tools.scatter_and_gather(
1611  alignment_coordinates)
1612  rmsd_coordinates=IMP.pmi.tools.scatter_and_gather(
1613  rmsd_coordinates)
1614 
1615  if self.rank == 0:
1616  # save needed informations in external files
1617  self.save_objects(
1618  [best_score_feature_keyword_list_dict,
1619  rmf_file_name_index_dict],
1620  ".macro.pkl")
1621 
1622 # ------------------------------------------------------------------------
1623 # Calculate distance matrix and cluster
1624 # ------------------------------------------------------------------------
1625  print("setup clustering class")
1626  Clusters = IMP.pmi.analysis.Clustering(rmsd_weights)
1627 
1628  for n, model_coordinate_dict in enumerate(all_coordinates):
1629  template_coordinate_dict = {}
1630  # let's try to align
1631  if alignment_components is not None and len(Clusters.all_coords) == 0:
1632  # set the first model as template coordinates
1633  Clusters.set_template(alignment_coordinates[n])
1634  Clusters.fill(all_rmf_file_names[n], rmsd_coordinates[n])
1635  print("Global calculating the distance matrix")
1636 
1637  # calculate distance matrix, all against all
1638  Clusters.dist_matrix()
1639 
1640  # perform clustering and optionally display
1641  if self.rank == 0:
1642  Clusters.do_cluster(number_of_clusters)
1643  if display_plot:
1644  if self.rank == 0:
1645  Clusters.plot_matrix(figurename=os.path.join(outputdir,'dist_matrix.pdf'))
1646  if exit_after_display:
1647  exit()
1648  Clusters.save_distance_matrix_file(file_name=distance_matrix_file)
1649 
1650 # ------------------------------------------------------------------------
1651 # Alteratively, load the distance matrix from file and cluster that
1652 # ------------------------------------------------------------------------
1653  else:
1654  if self.rank==0:
1655  print("setup clustering class")
1656  Clusters = IMP.pmi.analysis.Clustering()
1657  Clusters.load_distance_matrix_file(file_name=distance_matrix_file)
1658  print("clustering with %s clusters" % str(number_of_clusters))
1659  Clusters.do_cluster(number_of_clusters)
1660  [best_score_feature_keyword_list_dict,
1661  rmf_file_name_index_dict] = self.load_objects(".macro.pkl")
1662  if display_plot:
1663  if self.rank == 0:
1664  Clusters.plot_matrix(figurename=os.path.join(outputdir,'dist_matrix.pdf'))
1665  if exit_after_display:
1666  exit()
1667  if self.number_of_processes > 1:
1668  self.comm.Barrier()
1669 
1670 # ------------------------------------------------------------------------
1671 # now save all informations about the clusters
1672 # ------------------------------------------------------------------------
1673 
1674  if self.rank == 0:
1675  print(Clusters.get_cluster_labels())
1676  for n, cl in enumerate(Clusters.get_cluster_labels()):
1677  print("rank %s " % str(self.rank))
1678  print("cluster %s " % str(n))
1679  print("cluster label %s " % str(cl))
1680  print(Clusters.get_cluster_label_names(cl))
1681 
1682  # first initialize the Density class if requested
1683  if density_custom_ranges:
1684  DensModule = IMP.pmi.analysis.GetModelDensity(
1685  density_custom_ranges,
1686  voxel=voxel_size)
1687 
1688  dircluster = outputdir + "/cluster." + str(n) + "/"
1689  try:
1690  os.mkdir(dircluster)
1691  except:
1692  pass
1693 
1694  rmsd_dict = {"AVERAGE_RMSD":
1695  str(Clusters.get_cluster_label_average_rmsd(cl))}
1696  clusstat = open(dircluster + "stat.out", "w")
1697  for k, structure_name in enumerate(Clusters.get_cluster_label_names(cl)):
1698  # extract the features
1699  tmp_dict = {}
1700  tmp_dict.update(rmsd_dict)
1701  index = rmf_file_name_index_dict[structure_name]
1702  for key in best_score_feature_keyword_list_dict:
1703  tmp_dict[
1704  key] = best_score_feature_keyword_list_dict[
1705  key][
1706  index]
1707 
1708  # get the rmf name and the frame number from the list of
1709  # frame names
1710  rmf_name = structure_name.split("|")[0]
1711  rmf_frame_number = int(structure_name.split("|")[1])
1712  clusstat.write(str(tmp_dict) + "\n")
1713 
1714  # extract frame (open or link to existing)
1715  if k==0:
1716  prots,rs = IMP.pmi.analysis.get_hiers_and_restraints_from_rmf(
1717  self.model,
1718  rmf_frame_number,
1719  rmf_name)
1720  else:
1721  linking_successful = IMP.pmi.analysis.link_hiers_and_restraints_to_rmf(
1722  self.model,
1723  prots,
1724  rs,
1725  rmf_frame_number,
1726  rmf_name)
1727  if not linking_successful:
1728  continue
1729  if not prots:
1730  continue
1731 
1732  if IMP.pmi.get_is_canonical(prots[0]):
1733  states = IMP.atom.get_by_type(prots[0],IMP.atom.STATE_TYPE)
1734  prot = states[state_number]
1735  else:
1736  prot = prots[state_number]
1737 
1738  # transform clusters onto first
1739  if k > 0:
1740  model_index = Clusters.get_model_index_from_name(
1741  structure_name)
1742  transformation = Clusters.get_transformation_to_first_member(
1743  cl,
1744  model_index)
1745  rbs = set()
1746  for p in IMP.atom.get_leaves(prot):
1747  if not IMP.core.XYZR.get_is_setup(p):
1749  IMP.core.XYZR(p).set_radius(0.0001)
1750  IMP.core.XYZR(p).set_coordinates((0, 0, 0))
1751 
1753  rb = IMP.core.RigidBodyMember(p).get_rigid_body()
1754  rbs.add(rb)
1755  else:
1757  transformation)
1758  for rb in rbs:
1759  IMP.core.transform(rb,transformation)
1760 
1761  # add the density
1762  if density_custom_ranges:
1763  DensModule.add_subunits_density(prot)
1764 
1765  # pdb writing should be optimized!
1766  o = IMP.pmi.output.Output()
1767  o.init_pdb(dircluster + str(k) + ".pdb", prot)
1768  o.write_pdb(dircluster + str(k) + ".pdb")
1769 
1770  if IMP.pmi.get_is_canonical(prot):
1771  # create a single-state System and write that
1773  h.set_name("System")
1774  h.add_child(prot)
1775  o.init_rmf(dircluster + str(k) + ".rmf3", [h], rs)
1776  else:
1777  o.init_rmf(dircluster + str(k) + ".rmf3", [prot],rs)
1778  o.write_rmf(dircluster + str(k) + ".rmf3")
1779  o.close_rmf(dircluster + str(k) + ".rmf3")
1780 
1781  del o
1782  # IMP.atom.destroy(prot)
1783 
1784  if density_custom_ranges:
1785  DensModule.write_mrc(path=dircluster)
1786  del DensModule
1787 
1788  if self.number_of_processes>1:
1789  self.comm.Barrier()
1790 
1791  def save_objects(self, objects, file_name):
1792  import pickle
1793  outf = open(file_name, 'wb')
1794  pickle.dump(objects, outf)
1795  outf.close()
1796 
1797  def load_objects(self, file_name):
1798  import pickle
1799  inputf = open(file_name, 'rb')
1800  objects = pickle.load(inputf)
1801  inputf.close()
1802  return objects
A class to simplify create of constraints and movers for an IMP Hierarchy.
A macro for running all the basic operations of analysis.
Definition: macros.py:1249
def get_restraint_set
Get a RestraintSet containing all PMI restraints added to the model.
Definition: tools.py:49
Sample using molecular dynamics.
Definition: samplers.py:349
A member of a rigid body, it has internal (local) coordinates.
Definition: rigid_bodies.h:358
static bool get_is_setup(const IMP::ParticleAdaptor &p)
Definition: rigid_bodies.h:359
def __init__
Constructor.
Definition: macros.py:599
static XYZR setup_particle(Model *m, ParticleIndex pi)
Definition: XYZR.h:48
Utility classes and functions for reading and storing PMI files.
def get_best_models
Given a list of stat files, read them all and find the best models.
A helper output for model evaluation.
Miscellaneous utilities.
Definition: tools.py:1
def get_molecules
Return list of all molecules grouped by state.
Definition: macros.py:522
def __init__
Constructor.
Definition: macros.py:68
def clustering
Get the best scoring models, compute a distance matrix, cluster them, and create density maps...
Definition: macros.py:1348
A macro to build a IMP::pmi::topology::System based on a TopologyReader object.
Definition: macros.py:413
GenericHierarchies get_leaves(Hierarchy mhd)
Get all the leaves of the bit of hierarchy.
This class initializes the root node of the global IMP.atom.Hierarchy.
A class to cluster structures.
Representation of the system.
static bool get_is_setup(const IMP::ParticleAdaptor &p)
Definition: XYZR.h:47
def get_modeling_trajectory
Get a trajectory of the modeling run, for generating demonstrative movies.
Definition: macros.py:1295
static Hierarchy setup_particle(Model *m, ParticleIndex pi, ParticleIndexesAdaptor children=ParticleIndexesAdaptor())
Create a Hierarchy of level t by adding the needed attributes.
def get_trajectory_models
Given a list of stat files, read them all and find a trajectory of models.
def build_model
Create model.
Definition: macros.py:872
Performs alignment and RMSD calculation for two sets of coordinates.
Definition: pmi/Analysis.py:23
def scatter_and_gather
Synchronize data over a parallel run.
Definition: tools.py:1042
void transform(XYZ a, const algebra::Transformation3D &tr)
Apply a transformation to the particle.
A decorator for a particle with x,y,z coordinates.
Definition: XYZ.h:30
def BuildModel0
Construct a component for each subunit (no splitting, nothing fancy).
Definition: macros.py:1165
Class for easy writing of PDBs, RMFs, and stat files.
Definition: output.py:32
def __init__
Constructor.
Definition: macros.py:859
A macro to build a Representation based on a Topology and lists of movers DEPRECATED - Use BuildSyste...
Definition: macros.py:577
Tools for clustering and cluster analysis.
Definition: pmi/Analysis.py:1
bool get_is_canonical(atom::Hierarchy h)
Walk up a PMI2 hierarchy/representations and check if the root is named System.
Definition: utilities.h:98
Transformation3D get_identity_transformation_3d()
Return a transformation that does not do anything.
Classes for writing output files and processing them.
Definition: output.py:1
def deprecated_object
Python decorator to mark a class as deprecated.
Definition: __init__.py:9614
Sampling of the system.
Definition: samplers.py:1
Sample using Monte Carlo.
Definition: samplers.py:36
Create movers and setup constraints for PMI objects.
def add_state
Add a state using the topology info in a IMP::pmi::topology::TopologyReader object.
Definition: macros.py:456
Deprecated building macro - use BuildSystem()
Definition: macros.py:855
The general base class for IMP exceptions.
Definition: exception.h:49
Class to handle individual model particles.
Definition: Particle.h:37
def execute_macro
Builds representations and sets up degrees of freedom.
Definition: macros.py:528
def read_coordinates_of_rmfs
Read in coordinates of a set of RMF tuples.
def __init__
Constructor.
Definition: macros.py:439
void add_geometries(RMF::FileHandle file, const display::GeometriesTemp &r)
Add geometries to the file.
Set up the representation of all proteins and nucleic acid macromolecules.
A macro to help setup and run replica exchange.
Definition: macros.py:23
A dictionary-like wrapper for reading and storing sequence data.
Hierarchies get_leaves(const Selection &h)
Compute mean density maps from structures.
Support for the RMF file format for storing hierarchical molecular data and markup.
def get_residue_indexes
Retrieve the residue indexes for the given particle.
Definition: tools.py:1000
static bool get_is_setup(const IMP::ParticleAdaptor &p)
Definition: rigid_bodies.h:490
Sample using replica exchange.
Definition: samplers.py:452
Inferential scoring building on methods developed as part of the Inferential Structure Determination ...
def get_representation
Return the Representation object.
Definition: macros.py:713
A decorator for a particle with x,y,z coordinates and a radius.
Definition: XYZR.h:27