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