IMP logo
IMP Reference Guide  2.8.0
The Integrative Modeling Platform
BrownianDynamics.h
Go to the documentation of this file.
1 /**
2  * \file IMP/atom/BrownianDynamics.h
3  * \brief Simple molecular dynamics optimizer.
4  *
5  * Copyright 2007-2017 IMP Inventors. All rights reserved.
6  *
7  */
8 
9 #ifndef IMPATOM_BROWNIAN_DYNAMICS_H
10 #define IMPATOM_BROWNIAN_DYNAMICS_H
11 
12 #include <IMP/atom/atom_config.h>
13 #include "Diffusion.h"
14 #include "Simulator.h"
15 #include "atom_macros.h"
16 #include <IMP/Optimizer.h>
17 #include <IMP/Particle.h>
18 //#include <IMP/utility.h>
19 #include <IMP/internal/units.h>
20 #include <IMP/algebra/Vector3D.h>
21 
22 IMPATOM_BEGIN_NAMESPACE
23 
24 #ifdef IMP_KERNEL_CUDA_LIB
25 #define IMP_ATOM_DEFAULT_BD_RANDOM_POOL_SIZE 1000000
26 #else
27 #define IMP_ATOM_DEFAULT_BD_RANDOM_POOL_SIZE 10000
28 #endif
29 
30 // for swig
31 class SimulationParameters;
32 
33 //! Simple Brownian dynamics simulator.
34 /** This is an implementation of a Brownian Dynamics simulator.
35 
36  _Input particles and score_
37 
38  Each optimized particle must have x,y,z attributes
39  that are optimizable. In addition, each optimized particle must be
40  decorated with the Diffusion decorator. Optionally, the
41  RigidBodyDiffusion decorator can be used to specify a rotational
42  diffusion coefficient for core::RigidBody particles. The
43  optimizer assumes the scoring function to be energy in kcal/mol, and the xyz
44  coordinates to be in angstroms and the diffusion coefficient of
45  each particle be in \f$A^2/fs\f$ (or \f$Radian^2/fs\f$ for rotational
46  diffusion coefficient). Particles without optimized x,y,z
47  and nonoptimized D are skipped.
48 
49  The optimizer can either automatically determine which particles
50  to use from the model or be passed a SingletonContainer for the
51  particles. If such a container is passed, particles added to it
52  during optimization state updates are handled properly.
53 
54  _Simulation_
55 
56  At each simulation time step, each particle is translated in the
57  direction of the sum of a random diffusion vector and the gradient
58  of the scoring function (force field) at the particle
59  coordinates. The translation is proportional to the particle
60  diffusion coefficient, the time step size, and the inverse of kT.
61  Note that particles masses are not considered, only their
62  diffusion coefficients.
63 
64  Similarly, rigid bodies are rotated by the sum of a random torque and a
65  force field torque, proportionally to the rotational diffusion
66  coefficient, the time step size, and inversely proportional kT.
67 
68  If the skt (stochastic Runge Kutta) flag is true, the simulation is
69  altered slightly to apply the SKT scheme.
70 
71  \see Diffusion
72  \see RigidBodyDiffusion
73  */
74 class IMPATOMEXPORT BrownianDynamics : public Simulator {
75  private:
76 
77  double max_step_;
78  bool srk_;
80  IMP::Vector<double> random_pool_; // pool of random doubles ~N(0.0,1.0)
81  unsigned int i_random_pool_; // poistion in pool of random numbers
82 
83  public:
84 
85  //! Create the optimizer
86  /** If sc is not null, that container will be used to find particles
87  to move, otherwise the model will be searched.
88  @param m model associated with bd
89  @param name name of bd object
90  @param wave_factor for wave step function, see Simulator object,
91  if >1.001 or so, creates a wave of time steps
92  that are larger by up to wave_factor from
93  formal maximal time step
94  @param random_pool_size number of random numbers in internal pool
95  used to accelerate random number generation.
96  Memory requirement scales accordingly.
97 
98  @note wave_factor is an advanced feature - if you're not sure, just use
99  its default, see also Simulator::simulate_wave()
100  */
101  BrownianDynamics(Model *m, std::string name = "BrownianDynamics%1%",
102  double wave_factor = 1.0,
103  unsigned int random_pool_size=IMP_ATOM_DEFAULT_BD_RANDOM_POOL_SIZE);
104  void set_maximum_move(double ms) { max_step_ = ms; }
105  void set_use_stochastic_runge_kutta(bool tf) { srk_ = tf; }
106 
108 
109  protected:
110  /** a set of setup operations before a series of simulation steps */
111  virtual void setup(const ParticleIndexes &ps) IMP_OVERRIDE;
112 
113  /** Calls do_advance_chunk() to advance ps in chunks
114 
115  @param sc particles to simulate in this step
116  @param dt maximal step size in femtoseconds
117 
118  @return the time step actually simulated (for this class,
119  it is always equal to the inut dt)
120  */
121  virtual double do_step(const ParticleIndexes &sc,
122  double dt) IMP_OVERRIDE;
123 
124  virtual bool get_is_simulation_particle(ParticleIndex p) const
125  IMP_OVERRIDE;
126 
127  protected:
128  /** advances a chunk of ps from index begin to end
129 
130  @param dtfs time step in femtoseconds
131  @param ikt inverse kT for current chunk step
132  @param ps particle indexes to advance
133  @param begin beginning index of chunk of ps
134  @param end end index of chunk of ps
135  */
136  virtual void do_advance_chunk(double dtfs, double ikt,
137  const ParticleIndexes &ps,
138  unsigned int begin, unsigned int end);
139 
140 
141  protected:
142  //! returns the maximal step size allowed in this simulation
143  double get_max_step() const { return max_step_; }
144 
145  //! returns true if implementing the Stochastic Runga-Kutta
146  //! Brownian Dynamics variant
147  bool get_is_srk() const { return srk_; }
148 
149 
150  //! regenerate internal cached pool of random numbers
151  void reset_random_pool();
152 
153  //! returns a normally distributed sample with
154  //! mean zero and std-dev sigma
155  double get_sample(double sigma)
156  {
157  if(i_random_pool_ >= random_pool_.size()){
158  reset_random_pool();
159  }
160  return random_pool_[i_random_pool_++]*sigma;
161  }
162 
163 
164 #ifndef SWIG
165  //! set the force felt on particle i to f
166  void set_force(unsigned int i, algebra::Vector3D const& f)
167  { forces_[i]=f; }
168 #endif
169 
170  //! get the force vectors felt on each particle in kCal/mol/A
172  { return forces_; }
173 
174  //! get the force felt on particle i in kCal/mol/A
175  algebra::Vector3D const& get_force(unsigned int i) const
176  { return forces_[i]; }
177 
178 private:
179  void advance_coordinates_1(ParticleIndex pi, unsigned int i,
180 
181  double dtfs, double ikT);
182  void advance_coordinates_0(ParticleIndex pi, unsigned int i,
183  double dtfs, double ikT);
184  void advance_orientation_0(ParticleIndex pi, double dtfs, double ikT);
185 
186 };
187 
188 /** Repeatedly run the current model with Brownian dynamics at different time
189  steps to try to find the maximum time step that can be used without
190  the model exploding.
191 */
192 IMPATOMEXPORT double get_maximum_time_step_estimate(BrownianDynamics *bd);
193 
194 #ifndef IMP_DOXYGEN
195 IMPATOMEXPORT double get_harmonic_sigma(double D, double f);
196 #endif
197 IMPATOM_END_NAMESPACE
198 
199 #endif /* IMPATOM_BROWNIAN_DYNAMICS_H */
The base class for simulators.
Definition: Simulator.h:34
A decorator for a diffusing particle.
algebra::Vector3D const & get_force(unsigned int i) const
get the force felt on particle i in kCal/mol/A
#define IMP_OBJECT_METHODS(Name)
Define the basic things needed by any Object.
Definition: object_macros.h:25
algebra::Vector3Ds const & get_forces() const
get the force vectors felt on each particle in kCal/mol/A
virtual bool get_is_simulation_particle(ParticleIndex p) const =0
Return true if the passed particle is appropriate for the simulation.
double get_max_step() const
returns the maximal step size allowed in this simulation
Base class for all optimizers.
virtual double do_step(const ParticleIndexes &sc, double dt)=0
Perform a single time step.
Class for storing model, its restraints, constraints, and particles.
Definition: Model.h:72
virtual void setup(const ParticleIndexes &)
Definition: Simulator.h:169
Simple Brownian dynamics simulator.
Simple molecular dynamics optimizer.
Classes to handle individual model particles. (Note that implementation of inline functions is in int...
Various important macros for implementing decorators.
VectorD< 3 > Vector3D
Definition: VectorD.h:395
Simple 3D vector class.
void set_force(unsigned int i, algebra::Vector3D const &f)
set the force felt on particle i to f
double get_maximum_time_step_estimate(BrownianDynamics *bd)
double get_sample(double sigma)
#define IMP_OVERRIDE
Cause a compile error if this method does not override a parent method.