IMP logo
IMP Reference Guide  2.9.0
The Integrative Modeling Platform
SimulationData.h
Go to the documentation of this file.
1 /**
2  * \file npctransport/SimulationData.h
3  * \brief description
4  *
5  * Copyright 2007-2018 IMP Inventors. All rights reserved.
6  */
7 
8 #ifndef IMPNPCTRANSPORT_SIMULATION_DATA_H
9 #define IMPNPCTRANSPORT_SIMULATION_DATA_H
10 
11 
12 #include "npctransport_config.h"
13 #include <IMP/Model.h>
14 #include <IMP/PairContainer.h>
16 #include <IMP/atom/Hierarchy.h>
17 #include <IMP/SingletonContainer.h>
19 //#include <IMP/container/DynamicListSingletonContainer.h>
24 #include <IMP/core/Typed.h>
27 #include <IMP/Pointer.h>
28 #include <IMP/algebra/Vector3D.h>
29 #include <IMP/algebra/Sphere3D.h>
30 #include <boost/unordered_map.hpp>
31 #include <boost/unordered_set.hpp>
32 #include "io.h"
33 #include "Parameter.h"
34 #include "Scoring.h"
35 #include "Statistics.h"
36 #include "npctransport_proto.fwd.h"
37 #include <string>
38 
39 IMPNPCTRANSPORT_BEGIN_NAMESPACE
40 
41 // Version 2.5 - turned slab into a particle
42 // Version 3.0 - added harmonic bond with proper tau
43 // Version 4.0 - more properly handle various suffixes of FG chains (=different types of beads)
44 #define IMPNPCTRANSPORT_VERSION 4.0
45 
46 //! Store all parameters for a simulation
47 class IMPNPCTRANSPORTEXPORT SimulationData : public Object {
48  private:
49  // params
50  Parameter<double> output_npctransport_version_;
51  Parameter<double> box_side_;
52  Parameter<double> tunnel_radius_; // note this is the initial pore radius (major radius if toroidal pore)
53  Parameter<double> tunnel_radius_k_; // k for harmonic restraint on pore radius
54  Parameter<double> pore_anchored_beads_k_; // k for harmonic restraint on beads anchored to pore
55  Parameter<double> slab_thickness_; // note this is the initial thickness (also minor vertical radius if toroidal pore)
56  Parameter<bool> box_is_on_;
57  Parameter<int> slab_is_on_;
58  Parameter<int> number_of_trials_;
59  Parameter<int> number_of_frames_;
60  Parameter<int> dump_interval_frames_;
61  Parameter<bool> is_backbone_harmonic_;
62  Parameter<double> backbone_tau_ns_;
63  Parameter<double> angular_d_factor_;
64  Parameter<double> range_;
65  Parameter<double> statistics_fraction_;
66  Parameter<int> statistics_interval_frames_;
67  Parameter<int> output_statistics_interval_frames_;
68  Parameter<double> time_step_;
69  Parameter<double> time_step_wave_factor_;
70  Parameter<double> maximum_number_of_minutes_;
71  Parameter<double> fg_anchor_inflate_factor_; // By how much to inflate static
72  // anchors of FG nups.
73  Parameter<int> is_exclude_floaters_from_slab_initially_;
74  Parameter<double> are_floaters_on_one_slab_side_;
75  Parameter<int> is_xyz_hist_stats_;
76 
77  // time when simulation has started for this process
78  Parameter<double> initial_simulation_time_ns_;
79  Parameter<double> temperature_k_;
80 
81  public:
82  double get_output_npctransport_version() const { return output_npctransport_version_; }
83 
84  /** returns the maximal interaction range between particles */
85  double get_range() const { return range_; }
86 
87  int get_output_statistics_interval_frames()
88  { return output_statistics_interval_frames_; }
89 
90  /** returns whether should exclude floaters from slab during initialization */
92  { return is_exclude_floaters_from_slab_initially_; }
93 
94  /** returns whether should exclude floaters from one slab side,
95  if excluded at all
96  */
98  { return are_floaters_on_one_slab_side_;}
99 
100  bool get_is_xyz_hist_stats()
101  { return is_xyz_hist_stats_; }
102 
103  /** returns the simulation angular d factor */
104  double get_angular_d_factor() const
105  { return angular_d_factor_; }
106 
107 
108  private:
109  // the model on which the simulation is run
111 
112  // The BrownianDynamic simulator
114 
115  // The scoring function wrapper for the simulation
117  scoring_;
118 
119  // The statistics manger for this simulation
121  statistics_;
122 
123  // all beads in the simulation (=fine-level particles)
124  Particles beads_;
125 
126  // a writer to an RMF (Rich Molecular Format) type file
128 
129  // the root of the model hierarchy
131 
132  // Membrane slab, if exists (mutable but is expected to be accessed only from get_slab_particle()
133  mutable PointerMember<Particle> slab_particle_;
134 
135  // fg bead types - a set of all fg/floater/obstacle types that were
136  // added via create_fgs/floaters/obstacles(), including
137  // the suffixes of individual beads in a chain added to their chain root name
138  ParticleTypeSet fg_bead_types_;
139  // fg chain types - same as fg_bead_types but just for the root particles
140  // without the suffixes
141  ParticleTypeSet fg_chain_types_;
142  ParticleTypeSet floater_types_;
143  ParticleTypeSet obstacle_types_;
144 
146 
147  boost::unordered_map<core::ParticleType, algebra::Sphere3Ds> sites_;
148 
149  boost::unordered_map<core::ParticleType, double> ranges_;
150 
151  // the RMF format file to which simulation output is dumped:
152  std::string rmf_file_name_;
153 
154  // whether to save restraints to rmf when
155  // calling get_rmf_sos_writer()
156  bool is_save_restraints_to_rmf_;
157 
158  private:
159 
160  //! initialize slab_particle_ based on current parameters
161  void create_slab_particle();
162 
163  /**
164  Adds the FG Nup chains to the model hierarchy,
165  based on the settings in data
166 
167  @param fg_data data for fgs of this type in protobuf format as specified
168  in data/npctransport.proto
169  */
170  void create_fgs
171  ( const ::npctransport_proto::Assignment_FGAssignment &fg_data );
172 
173  /**
174  Adds the 'floaters' (free diffusing particles) to the model hierarchy,
175  based on the settings in data
176 
177  @param f_data data for floater in protobuf format as specified
178  in data/npctransport.proto
179  @param color a color for all floaters of this type
180  */
181  void create_floaters(
182  const ::npctransport_proto::Assignment_FloaterAssignment &f_data,
183  display::Color color);
184 
185  /**
186  Adds the 'obstacles' (possibly static e.g. nups that make the pore)
187  to the model hierarchy, based on the settings in data
188 
189  @param o_data data for obstacles in protobuf format as specified in
190  data/npctransport.proto
191  */
192  void create_obstacles
193  ( const ::npctransport_proto::Assignment_ObstacleAssignment &o_data);
194 
195  /** Initializes the simulation based on the specified assignment file
196 
197  @param prev_output_file protobuf file with the current assignment, and
198  possibly previous progress and statistics info
199  @param[in] new_output_file new protobuf file for subsequent output and
200  statistics.
201  If it is "", then prev_output_file is being
202  rewritten.
203  @param quick whether to perform a quick simulation with a small number
204  of iterations
205  */
206  void initialize(std::string prev_output_file,
207  std::string new_output_file,
208  bool quick);
209 
210  public:
211  /**
212  @param[in] prev_output_file name of protobuf file that encapsulates the
213  assignment data for initializing the
214  simulation params and the output
215  statistics file. The output file is
216  assumed to already exist at this
217  point and contain the assignment.
218  If it contains an rmf_conformation
219  or conformation field, it will be
220  used to initialize particle
221  positions, in that priority.
222  @param[in] quick if true, perform a very short simulation,
223  typically for calibration or for initial testing
224  @param[out] rmf_file_name RMF file to which simulation trajectory
225  is recorded (if empty, no RMF file is created
226  upon construction)
227  \see set_rmf_file()
228  @param[in] new_output_file new protobuf file for subsequent output and statistics
229  If it is "", then prev_output_file is being rewritten
230  (= default).
231  */
232  SimulationData(std::string prev_output_file, bool quick,
233  std::string rmf_file_name = std::string(),
234  std::string new_output_file = "");
235 
236  Model *get_model();
237 #ifndef SWIG
238  Model * get_model() const {
239  IMP_USAGE_CHECK(m_, "model not initialized in get_model()");
240  return m_;
241  }
242 #endif
243 
244  /** returns a scoring object that is updated with the current particles
245  hierarch and associated with this sd
246  */
247  Scoring * get_scoring();
248 
249  /** returns a Statistics object that is updated by this simulation data
250  */
251  Statistics* get_statistics();
252 
253 #ifndef SWIG
254  /** returns a scoring object that is updated with the current particles
255  hierarchy and associated with this sd
256  */
257  Scoring const* get_scoring() const;
258 
259  /** returns a Statistics object that is updated by this simulation data
260  */
261  Statistics const* get_statistics() const;
262 
263 #endif
264 
265  /** gets the Brownian Dynamics object that is capable of simulating
266  this data. Statistics are not yet active for the simulation.
267 
268  @param recreate if true, forces recreation of the bd object
269  */
270  atom::BrownianDynamics *get_bd(bool recreate = false);
271 
272  //! activates Brownian Dynamics statistics tracking
273  //! by adding all appropriate optimizer states, if they weren't already
274  void activate_statistics();
275 
276  /** returns the requested fraction of time for taking statistics */
277  double get_statistics_fraction() const { return statistics_fraction_; }
278 
279  /** returns true if bead has an fg type
280  that is, particle was added within create_fgs()
281  */
282  bool get_is_fg_bead(ParticleIndex pi) const;
283 
284 
285  /** returns true if particle type is of fg type
286  (that is, it is one of the types added via create_fgs())
287  */
289  return fg_bead_types_.find(pt) != fg_bead_types_.end();
290  }
291 
292  /** returns true if particle has an fg chain type
293  that is, particle has the same type as the root of an
294  fg chain that was added via create_fgs()
295  */
296  bool get_is_fg_chain(ParticleIndex pi) const;
297 
298 
299  /** returns true if particle type is an fg chain type
300  (that is, it is one of the types of a chain added
301  via create_fgs())
302  */
304  return fg_chain_types_.find(pt) != fg_chain_types_.end();
305  }
306 
307 
308  /** return all the types of fg beads that were added
309  via create_fgs() including the suffix appended to their
310  chain type name */
311  ParticleTypeSet const& get_fg_bead_types() const {
312  return fg_bead_types_;
313  }
314 
315  /** return all the types of fg chains that were added
316  via create_fgs(), i.e. the types of chain roots */
317  ParticleTypeSet const& get_fg_chain_types() const {
318  return fg_chain_types_;
319  }
320 
321 
322  /** return all the types of floaters that were added
323  via create_floaters() */
324  ParticleTypeSet const& get_floater_types() const {
325  return floater_types_;
326  }
327 
328  /** return all the types of obstacles that were added
329  via create_obstacles() */
330  ParticleTypeSet const& get_obstacle_types() const {
331  return obstacle_types_;
332  }
333 
334  /**
335  retrieve fg chain hierarchies
336 
337  @return all the fg hierarchies in the simulation data object
338  that stand for individual FG chains
339  */
340  atom::Hierarchies get_fg_chain_roots() const;
341 
342  /** \deprecated_at{2.2}, use get_fg_chain_roots() instead */
343  IMPCORE_DEPRECATED_METHOD_DECL(2.2)
344  atom::Hierarchies get_fg_chains() const
345  { return get_fg_chain_roots(); }
346 
347 
348  /** return all the obstacle particles */
349  ParticlesTemp get_obstacle_particles() const;
350 
351 
352 #ifndef SWIG
353 
354  /**
355  returns all diffusing fine-level beads in this
356  simulation data object (or an empty list if none exists)
357  @note efficient - returns an existing container by ref
358  */
359  Particles& get_beads_byref(){ return beads_; }
360 #endif
361 
362  /**
363  returns all diffusing fine-level beads in this
364  simulation data object (or an empty list if none exists)
365  */
366  ParticlesTemp get_beads(){ return beads_; }
367 
368 
369  /**
370  returns all fine-level beads that currently exist in
371  this simulation data object get_sd(), which are also optimizable
372  (or an empty list if no such bead exists)
373  */
374  ParticlesTemp get_optimizable_beads();
375 
376  /**
377  returns all fine-level beads that currently exist in
378  this simulation data object get_sd(), which are also non-optimizable
379  (or an empty list if no such bead exists)
380  */
381  ParticlesTemp get_non_optimizable_beads();
382 
383  bool get_is_backbone_harmonic() const
384  { return is_backbone_harmonic_; }
385 
386  double get_backbone_tau_ns() const {
387  return backbone_tau_ns_;
388  }
389 
390  double get_temperature_k() const
391  { return temperature_k_; }
392 
393  /** get time from which simulation begins */
395  {return initial_simulation_time_ns_; }
396 
397  /**
398  Create n interaction sites spread around the surface of a ball of
399  radius r, and associate these sites with all particles of type t0
400 
401  If n==-1, creates a single site, centered at the bead, and r is ignored
402 
403  @param t0 type of particle to associate with sites
404  @param n number of sites
405  @param r radius at which to position sites relative to t0 origin
406  @param sr radius from site center within which site
407  has maximal interaction energy
408  */
409  void set_sites(core::ParticleType t0, int n, double r, double sr);
410 
411  /**
412  returns a list of 3D coordinates for the interaction sites associated
413  with particles of type t0. The site coordinates are in the
414  particles local reference frame
415  */
417  algebra::Vector3Ds ret;
418  if (sites_.find(t0) != sites_.end()) {
419  algebra::Sphere3Ds sites = sites_.find(t0)->second;
420  for(unsigned int i = 0; i < sites.size(); i++){
421  ret.push_back(sites[i].get_center());
422  }
423  return ret;
424  } else {
425  return algebra::Vector3Ds();
426  }
427  }
428 
429  /**
430  returns a list of Spheres3D for the interaction sites associated
431  with particles of type t0. The site coordinates are in the
432  particles local reference frame.
433  */
435  if (sites_.find(t0) != sites_.end()) {
436  return sites_.find(t0)->second;
437  } else {
438  return algebra::Sphere3Ds();
439  }
440  }
441 
442  /** Set the sites explicitly. */
444  IMP_USAGE_CHECK(sites.size()>0, "trying to set zero sites for particle type"
445  << t0.get_string() );
446  sites_[t0] = sites;
447  }
448 
449  /**
450  Returns the effective interaction range radius of a
451  site on a floater */
452  // TODO: this is not the true value - the true one might depend on the
453  // specific range of each interaction type, so range_ is more
454  // like an upper bound
455  double get_site_display_radius(core::ParticleType) const { return std::max(range_ / 2, 6.0); }
456 
457  //! Return the maximum number of minutes the simulation can run
458  /** Or 0 for no limit. */
460  return maximum_number_of_minutes_;
461  }
462 
463 // swig doesn't equate the two protobuf types
464 #ifndef SWIG
465  /**
466  add the specified interaction to the scoring.
467  In additions, adds statistics optimizer state about pairs of this type.
468 
469  @param idata the protobuf data about the interaction (particle types,
470  interaction coefficients, etc.)
471  */
472  void add_interaction(
473  const ::npctransport_proto::Assignment_InteractionAssignment &idata);
474 #endif
475 
476  // get the bounding box for this simulation
477  algebra::BoundingBox3D get_box() const;
478 
479  double get_box_size() const;
480 
481  void set_box_size(double box_size);
482 
483  //* returns true if a slab is defnied */
484  bool get_has_slab() const { return slab_is_on_!=0; }
485 
486  /** returns true if a slab is defined and has a cylindrical pore */
488  { return slab_is_on_==1; }
489 
490  /** returns true if a slab is defined and has a toroidal pore */
492  { return slab_is_on_==2; }
493 
494  /** returns the slab particle, initializing it based on current parameters
495  if needed
496  */
498  if(!slab_particle_ && get_has_slab()){
499  // Const cast below - cause the initialization of slab_particle_ is transparent when
500  // using get_slab_particle(), and it should only be accessed from this method
501  const_cast<SimulationData*>(this)->create_slab_particle();
502  }
503  return slab_particle_.get();
504  }
505 
506  // get the cylinder in the slab for this simulation
507  algebra::Cylinder3D get_cylinder() const;
508 
509  bool get_has_bounding_box() const { return box_is_on_; }
510 
511  /**
512  Open the specified RMF file, links it to the hierarchies of this object, and
513  copy
514  the XYZ coordinates from its last frame into the particles of this object.
515 
516  @note this method assumes that the hierarchies that are stored in the RMF
517  file
518  was constructed in the same way as the hierarchies within this SimulationData
519  object. The only major difference between the two is assumed to be the
520  coordinates of
521  the various particles. Otherwise, results might be unexpected, and it is not
522  guaranteed
523  that a mismatch would be detected at runtime.
524  \see IMP::rmf::link_hierarchies()
525  \see IMP::rmf::load_frame()
526 
527  The frame number is the frame in the file to use, if, eg, there are several
528  initial frames to use. If it is -1, the last frame is used.
529 
530  @exception RMF::IOException if couldn't open RMF file, or unsupported file
531  format
532  */
533  void initialize_positions_from_rmf(RMF::FileConstHandle fh,
534  int frame_number = -1);
535 
536  /** Links a handle to an open rmf file with the
537  hierarchy and constraints of this simulation data,
538  so that writing frames into this handle will write
539  the updated state of the simulation
540 
541  @param fh the handle to be linked
542  @param is_restraints if true, save restraints to RMF
543 
544  \exception RMF::IOException IO error with RMF file handle
545  */
546  void link_rmf_file_handle(RMF::FileHandle fh, bool is_restraints = true);
547 
548  /**
549  Returns the internal periodic SaveOptimizerState writer that
550  periodically outputs the particles hierarchy and restraints,
551  using the file name returned by SimulationData::get_rmf_file_name(), which is
552  assumed to be non empty. If the writer does not exist, then
553  it is constructed first.
554 
555  \exception RMF::IOException couldn't create RMF file
556  */
557  rmf::SaveOptimizerState *get_rmf_sos_writer();
558 
559  /**
560  Resets the RMF file to which the internal periodic SaveOptimizerState
561  writer dumps the output. If it does not exist, create it.
562 
563  \exception RMF::IOException couldn't create RMF file or other IO related
564  problems with RMF
565  */
566  void reset_rmf();
567 
568 
569  /*
570  temporarily suspend output to RMF file
571 
572  @param suspend if true, suspend an existing rmf save optimizer state,
573  associated with the bd, otherwise restores an old one, if it
574  exists or was created.
575  @note the suspension will expire when calling get_bd(true)
576  */
577  void switch_suspend_rmf(bool suspend);
578 
579  /**
580  Write the geometry to the file path 'out'
581  */
582  void write_geometry(std::string out);
583 
584  void dump_geometry();
585 
586  /**
587  Returns the root of all chains of type 'type'
588  */
589  atom::Hierarchy get_root_of_type(core::ParticleType type) const;
590 
591  unsigned int get_number_of_frames() const { return number_of_frames_; }
592 
593  unsigned int get_number_of_trials() const { return number_of_trials_; }
594 
595  atom::Hierarchy get_root() const { return atom::Hierarchy(root_); }
596 
597  double get_slab_thickness() const;
598 
599  //! returns te current tunnel radius
600  double get_tunnel_radius() const;
601 
602  //! alias to get_tunnel_radius
603  double get_pore_radius() const{
604  return get_tunnel_radius();
605  }
606 
607  //! returns the force coefficient for keeping the tunnel radius in
608  //! equilibrium (radius is dynamic only if this coefficient is positive)
609  double get_tunnel_radius_k() const {
610  return tunnel_radius_k_;
611  }
612 
613  //! alias to get_tunnel_radius_k
614  double get_pore_radius_k() const {
615  return get_tunnel_radius_k();
616  }
617 
618  //! returns true if pore radius can change dynamically
620  return get_has_slab() && get_tunnel_radius_k()>0.0;
621  }
622 
623  double get_pore_anchored_beads_k() const {
624  return pore_anchored_beads_k_;
625  }
626 
627  display::Geometry *get_static_geometry();
628 
629  int get_rmf_dump_interval_frames() const { return dump_interval_frames_; }
630 
631  std::string get_rmf_file_name() const { return rmf_file_name_; }
632 
633  /**
634  resets the name of the RMF file that records the simulation. If
635  an old one exists from a previous call, then the previous RMF
636  writer is invalidated by this action (closed and flushed). If
637  the Brownian Dynamics object has already been initialized, a new
638  writer with the new name is added as an optimizer state to it
639  instead of the existing one.
640  TODO: make sure the old writer is indeed closed and flushed
641 
642  @param new_name the new name of the rmf file
643  @param is_save_restraints_to_rmf whether to save restraints to this
644  rmf file
645  */
646  void set_rmf_file(const std::string &new_name,
647  bool is_save_restraints_to_rmf = true);
648 
650 
651 };
652 
653 
654 inline WeakObjectKey get_simulation_data_key() {
655  static WeakObjectKey simdata("simulation data");
656  return simdata;
657 }
658 
659 IMPNPCTRANSPORT_END_NAMESPACE
660 
661 #endif /* IMPNPCTRANSPORT_SIMULATION_DATA_H */
Apply a PairScore to each Pair in a list.
Represent an RGB color.
Definition: Color.h:24
Particle * get_slab_particle() const
Simple molecular dynamics optimizer.
void set_sites(core::ParticleType t0, const algebra::Sphere3Ds &sites)
algebra::Vector3Ds get_site_centers(core::ParticleType t0) const
scoring associated with a SimulationData object
A container for Singletons.
Vector< VectorD< 3 > > Vector3Ds
Definition: VectorD.h:397
The base class for geometry.
#define IMP_OBJECT_METHODS(Name)
Define the basic things needed by any Object.
Definition: object_macros.h:25
Store a set of PairContainers.
Storage of a model, its restraints, constraints and particles.
double get_maximum_number_of_minutes() const
Return the maximum number of minutes the simulation can run.
statistics and order parameters about the simulations that is associated with a SimulationData object...
Represent a cylinder in 3D.
Definition: Cylinder3D.h:26
bool get_is_pore_radius_dynamic() const
returns true if pore radius can change dynamically
Score particles based on a bounding box.
Dump the state of all associated objects into the RMF file.
description
ParticleTypeSet const & get_fg_bead_types() const
Implement geometry for the basic shapes from IMP.algebra.
A particle with a user-defined type.
const std::string get_string() const
Turn a key into a pretty string.
Definition: Key.h:165
double get_site_display_radius(core::ParticleType) const
double get_initial_simulation_time_ns() const
Class for storing model, its restraints, constraints, and particles.
Definition: Model.h:72
Decorator for helping deal with a hierarchy of molecules.
A container for Pairs.
Simple Brownian dynamics simulator.
bool get_is_slab_with_cylindrical_pore() const
The standard decorator for manipulating molecular structures.
Common base class for heavy weight IMP objects.
Definition: Object.h:106
int get_number_of_frames(const ::npctransport_proto::Assignment &config, double time_step)
Return all pairs from a SingletonContainer.
A smart pointer to a ref-counted Object that is a class member.
Definition: Pointer.h:146
algebra::Sphere3Ds get_sites(core::ParticleType t0) const
IMP::Vector< Sphere3D > Sphere3Ds
Definition: SphereD.h:88
ParticleTypeSet const & get_floater_types() const
Store all parameters for a simulation.
ParticleTypeSet const & get_obstacle_types() const
double get_pore_radius() const
alias to get_tunnel_radius
bool get_is_fg_chain(core::ParticleType pt) const
double get_pore_radius_k() const
alias to get_tunnel_radius_k
A nullptr-initialized pointer to an IMP Object.
Periodically dump the state of all associated objects into the RMF file.
Define some predicates.
Hierarchy get_root(Hierarchy h)
Return the root of the hierarchy.
Simple 3D vector class.
Class to handle individual particles of a Model object.
Definition: Particle.h:41
ParticleTypeSet const & get_fg_chain_types() const
#define IMP_USAGE_CHECK(expr, message)
A runtime test for incorrect usage of a class or method.
Definition: check_macros.h:168
description
Simple 3D sphere class.
bool get_is_fg_bead(core::ParticleType pt) const
forward declaration of incomplete protobuf files for the main data structures used ...