IMP logo
IMP Reference Guide  2.9.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-2018 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  _Time step_
72  The time step is always equal precisely to Simulater::get_maximum_time_step()
73  when using either Simulator::simulate() or Optimizer::optimize()
74 
75  \see Diffusion
76  \see RigidBodyDiffusion
77  */
78 class IMPATOMEXPORT BrownianDynamics : public Simulator {
79  private:
80 
81  double max_step_in_A_;
82  bool srk_;
84  IMP::Vector<double> random_pool_; // pool of random doubles ~N(0.0,1.0)
85  unsigned int i_random_pool_; // poistion in pool of random numbers
86 
87  public:
88 
89  //! Create the optimizer
90  /** If sc is not null, that container will be used to find particles
91  to move, otherwise the model will be searched.
92  @param m model associated with bd
93  @param name name of bd object
94  @param wave_factor for wave step function, see Simulator object,
95  if >1.001 or so, creates a wave of time steps
96  that are larger by up to wave_factor from
97  formal maximal time step
98  @param random_pool_size number of random numbers in internal pool
99  used to accelerate random number generation.
100  Memory requirement scales accordingly.
101 
102  @note wave_factor is an advanced feature - if you're not sure, just use
103  its default, see also Simulator::simulate_wave()
104  */
105  BrownianDynamics(Model *m, std::string name = "BrownianDynamics%1%",
106  double wave_factor = 1.0,
107  unsigned int random_pool_size=IMP_ATOM_DEFAULT_BD_RANDOM_POOL_SIZE);
108  //! sets the maximum move in A along either x,y or z axes
109  void set_maximum_move(double ms_in_A) { max_step_in_A_ = ms_in_A; }
110  void set_use_stochastic_runge_kutta(bool tf) { srk_ = tf; }
111 
113 
114  protected:
115  /** a set of setup operations before a series of simulation steps */
116  virtual void setup(const ParticleIndexes &ps) IMP_OVERRIDE;
117 
118  /** Calls do_advance_chunk() to advance ps in chunks
119 
120  @param sc particles to simulate in this step
121  @param dt_fs step size in femtoseconds
122 
123  @return the time step actually simulated (for this class,
124  it is always equal to the input dt_fs)
125  */
126  virtual double do_step(const ParticleIndexes &sc,
127  double dt_fs) IMP_OVERRIDE;
128 
129  virtual bool get_is_simulation_particle(ParticleIndex p) const
130  IMP_OVERRIDE;
131 
132  protected:
133  /** advances a chunk of ps from index begin to end
134 
135  @param dtfs time step in femtoseconds
136  @param ikt inverse kT for current chunk step
137  @param ps particle indexes to advance
138  @param begin beginning index of chunk of ps
139  @param end end index of chunk of ps
140  */
141  virtual void do_advance_chunk(double dtfs, double ikt,
142  const ParticleIndexes &ps,
143  unsigned int begin, unsigned int end);
144 
145 
146  protected:
147  //! returns the maximal step size allowed in this simulation
148  //! in A along x, y or z axes
149  double get_max_step() const { return max_step_in_A_; }
150 
151  //! returns true if implementing the Stochastic Runga-Kutta
152  //! Brownian Dynamics variant
153  bool get_is_srk() const { return srk_; }
154 
155 
156  //! regenerate internal cached pool of random numbers
157  void reset_random_pool();
158 
159  //! returns a normally distributed sample with
160  //! mean zero and std-dev sigma
161  double get_sample(double sigma)
162  {
163  if(i_random_pool_ >= random_pool_.size()){
164  reset_random_pool();
165  }
166  return random_pool_[i_random_pool_++]*sigma;
167  }
168 
169 
170 #ifndef SWIG
171  //! set the force felt on particle i to f
172  void set_force(unsigned int i, algebra::Vector3D const& f)
173  { forces_[i]=f; }
174 #endif
175 
176  //! get the force vectors felt on each particle in kCal/mol/A
178  { return forces_; }
179 
180  //! get the force felt on particle i in kCal/mol/A
181  algebra::Vector3D const& get_force(unsigned int i) const
182  { return forces_[i]; }
183 
184 private:
185  void advance_coordinates_1(ParticleIndex pi, unsigned int i,
186 
187  double dtfs, double ikT);
188  void advance_coordinates_0(ParticleIndex pi, unsigned int i,
189  double dtfs, double ikT);
190  void advance_orientation_0(ParticleIndex pi, double dtfs, double ikT);
191 
192 };
193 
194 /** Repeatedly run the current model with Brownian dynamics at different time
195  steps to try to find the maximum time step that can be used without
196  the model exploding.
197 */
198 IMPATOMEXPORT double get_maximum_time_step_estimate(BrownianDynamics *bd);
199 
200 #ifndef IMP_DOXYGEN
201 IMPATOMEXPORT double get_harmonic_sigma(double D, double f);
202 #endif
203 IMPATOM_END_NAMESPACE
204 
205 #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
void set_maximum_move(double ms_in_A)
sets the maximum move in A along either x,y or z axes
#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.
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.