IMP logo
IMP Reference Guide  2.19.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-2022 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) override;
121 
122  //! Calls do_advance_chunk() to advance ps in chunks
123  /** @param sc particles to simulate in this step
124  @param dt_fs step size in femtoseconds
125 
126  @return the time step actually simulated (for this class,
127  it is always equal to the input dt_fs)
128  */
129  virtual double do_step(const ParticleIndexes &sc,
130  double dt_fs) override;
131 
132  virtual bool get_is_simulation_particle(ParticleIndex p) const
133  override;
134 
135  protected:
136  //! advances a chunk of ps from index begin to end
137  /** @param dtfs time step in femtoseconds
138  @param ikt inverse kT for current chunk step
139  @param ps particle indexes to advance
140  @param begin beginning index of chunk of ps
141  @param end end index of chunk of ps
142  */
143  virtual void do_advance_chunk(double dtfs, double ikt,
144  const ParticleIndexes &ps,
145  unsigned int begin, unsigned int end);
146 
147 
148  protected:
149  //! returns the maximal step size allowed in this simulation
150  //! in A along x, y or z axes
151  double get_max_step() const { return max_step_in_A_; }
152 
153  //! returns true if implementing the Stochastic Runga-Kutta
154  //! Brownian Dynamics variant
155  bool get_is_srk() const { return srk_; }
156 
157 
158  //! regenerate internal cached pool of random numbers
159  void reset_random_pool();
160 
161  //! returns a normally distributed sample with
162  //! mean zero and std-dev sigma
163  double get_sample(double sigma)
164  {
165  if(i_random_pool_ >= random_pool_.size()){
166  reset_random_pool();
167  }
168  return random_pool_[i_random_pool_++]*sigma;
169  }
170 
171 
172 #ifndef SWIG
173  //! set the force felt on particle i to f
174  void set_force(unsigned int i, algebra::Vector3D const& f)
175  { forces_[i]=f; }
176 #endif
177 
178  //! get the force vectors felt on each particle in kCal/mol/A
180  { return forces_; }
181 
182  //! get the force felt on particle i in kCal/mol/A
183  algebra::Vector3D const& get_force(unsigned int i) const
184  { return forces_[i]; }
185 
186 private:
187  void advance_coordinates_1(ParticleIndex pi, unsigned int i,
188 
189  double dtfs, double ikT);
190  void advance_coordinates_0(ParticleIndex pi, unsigned int i,
191  double dtfs, double ikT);
192  void advance_orientation_0(ParticleIndex pi, double dtfs, double ikT);
193 
194 };
195 
196 /** Repeatedly run the current model with Brownian dynamics at different time
197  steps to try to find the maximum time step that can be used without
198  the model exploding.
199 */
200 IMPATOMEXPORT double get_maximum_time_step_estimate(BrownianDynamics *bd);
201 
202 #ifndef IMP_DOXYGEN
203 IMPATOMEXPORT double get_harmonic_sigma(double D, double f);
204 #endif
205 IMPATOM_END_NAMESPACE
206 
207 #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:86
virtual void setup(const ParticleIndexes &)
Definition: Simulator.h:178
Simple Brownian dynamics simulator.
Base class for "simulators", such as molecular dynamics.
Classes to handle individual model particles. (Note that implementation of inline functions is in int...
Macros for maintaining molecular hierarchies.
VectorD< 3 > Vector3D
Definition: VectorD.h:425
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)