IMP logo
IMP Reference Guide  2.7.0
The Integrative Modeling Platform
atom/MolecularDynamics.h
Go to the documentation of this file.
1 /**
2  * \file IMP/atom/MolecularDynamics.h
3  * \brief Simple molecular dynamics optimizer.
4  *
5  * Copyright 2007-2017 IMP Inventors. All rights reserved.
6  *
7  */
8 
9 #ifndef IMPATOM_MOLECULAR_DYNAMICS_H
10 #define IMPATOM_MOLECULAR_DYNAMICS_H
11 
12 #include <IMP/atom/atom_config.h>
13 #include "Simulator.h"
14 #include "atom_macros.h"
15 #include <IMP/Particle.h>
16 #include <IMP/Optimizer.h>
17 
18 IMPATOM_BEGIN_NAMESPACE
19 
20 //! A particle with linear (XYZ) velocity
21 /** Typically this is used in combination with the MolecularDynamics optimizer.
22  */
23 class IMPATOMEXPORT LinearVelocity : public Decorator {
24  static void do_setup_particle(Model *m, ParticleIndex pi,
25  const algebra::Vector3D v = algebra::Vector3D(0, 0, 0)) {
26  m->add_attribute(get_velocity_key(0), pi, v[0]);
27  m->add_attribute(get_velocity_key(1), pi, v[1]);
28  m->add_attribute(get_velocity_key(2), pi, v[2]);
29  }
30 
31 public:
32  static FloatKey get_velocity_key(unsigned int i) {
33  IMP_USAGE_CHECK(i < 3, "Out of range coordinate");
34  static const FloatKey keys[] = {FloatKey("vx"), FloatKey("vy"),
35  FloatKey("vz")};
36  return keys[i];
37  }
38 
39  static bool get_is_setup(Model *m, ParticleIndex pi) {
40  return m->get_has_attribute(get_velocity_key(0), pi)
41  && m->get_has_attribute(get_velocity_key(1), pi)
42  && m->get_has_attribute(get_velocity_key(2), pi);
43  }
44 
48 
49  void set_velocity(const algebra::Vector3D &v) {
50  Model *m = get_model();
52  m->set_attribute(get_velocity_key(0), pi, v[0]);
53  m->set_attribute(get_velocity_key(1), pi, v[1]);
54  m->set_attribute(get_velocity_key(2), pi, v[2]);
55  }
56 
57  algebra::Vector3D get_velocity() const {
58  Model *m = get_model();
60  return algebra::Vector3D(m->get_attribute(get_velocity_key(0), pi),
61  m->get_attribute(get_velocity_key(1), pi),
62  m->get_attribute(get_velocity_key(2), pi));
63  }
64 };
65 
66 //! A particle with angular velocity
67 /** Typically this is used for RigidBody particles in combination
68  with the MolecularDynamics optimizer. The velocity is stored as
69  a quaternion.
70  */
71 class IMPATOMEXPORT AngularVelocity : public Decorator {
72  static void do_setup_particle(Model *m, ParticleIndex pi,
73  const algebra::Vector4D v = algebra::Vector4D(0, 0, 0, 0)) {
74  m->add_attribute(get_velocity_key(0), pi, v[0]);
75  m->add_attribute(get_velocity_key(1), pi, v[1]);
76  m->add_attribute(get_velocity_key(2), pi, v[2]);
77  m->add_attribute(get_velocity_key(3), pi, v[3]);
78  }
79 
80 public:
81  static FloatKey get_velocity_key(unsigned int i) {
82  IMP_USAGE_CHECK(i < 4, "Out of range coordinate");
83  static const FloatKey keys[] = {FloatKey("angvel0"), FloatKey("angvel1"),
84  FloatKey("angvel2"), FloatKey("angvel3")};
85  return keys[i];
86  }
87 
88  static bool get_is_setup(Model *m, ParticleIndex pi) {
89  return m->get_has_attribute(get_velocity_key(0), pi)
90  && m->get_has_attribute(get_velocity_key(1), pi)
91  && m->get_has_attribute(get_velocity_key(2), pi)
92  && m->get_has_attribute(get_velocity_key(3), pi);
93  }
94 
98 
99  void set_velocity(const algebra::Vector4D &v) {
100  Model *m = get_model();
102  m->set_attribute(get_velocity_key(0), pi, v[0]);
103  m->set_attribute(get_velocity_key(1), pi, v[1]);
104  m->set_attribute(get_velocity_key(2), pi, v[2]);
105  m->set_attribute(get_velocity_key(3), pi, v[3]);
106  }
107 
108  algebra::Vector4D get_velocity() const {
109  Model *m = get_model();
111  return algebra::Vector4D(m->get_attribute(get_velocity_key(0), pi),
112  m->get_attribute(get_velocity_key(1), pi),
113  m->get_attribute(get_velocity_key(2), pi),
114  m->get_attribute(get_velocity_key(3), pi));
115  }
116 };
117 
118 //! Simple molecular dynamics optimizer.
119 /** The particles to be optimized must have optimizable x,y,z attributes
120  and a non-optimizable mass attribute; this optimizer assumes the score
121  to be energy in kcal/mol, the xyz coordinates to be in angstroms, and
122  the mass to be in AMU (g/mol).
123 
124  \note RigidBody particles are not handled properly.
125 
126  Particles without optimized x,y,z and nonoptimized mass are skipped.
127  \see VelocityScalingOptimizerState
128  \see LangevinThermostatOptimizerState
129  \see BerendsenThermostatOptimizerState
130  \see RemoveRigidMotionOptimizerState
131  */
132 class IMPATOMEXPORT MolecularDynamics : public Simulator {
133  public:
134  /** Score based on the provided model */
136 
137  //! Return the current kinetic energy of the system, in kcal/mol
138  virtual Float get_kinetic_energy() const;
139 
140  //! Return the current kinetic temperature of the system
141  /** \param[in] ekinetic kinetic energy, e.g. from get_kinetic_energy()
142  */
143  Float get_kinetic_temperature(Float ekinetic) const;
144 
145  //! Set maximum velocity in A/fs
146  /** At each dynamics time step, the absolute value of each velocity
147  component is capped at this value. This prevents spurious strong forces
148  (occasionally encountered with frustrated conformations) from causing
149  large oscillations in the system.
150  By default, velocities are not capped.
151 
152  \note The actual velocities that are capped are the half-step velocities
153  in the velocity Verlet algorithm.
154  */
155  void set_velocity_cap(Float velocity_cap) { velocity_cap_ = velocity_cap; }
156 
157  //! Assign velocities representative of the given temperature
158  virtual void assign_velocities(Float temperature);
159  virtual void setup(const ParticleIndexes &ps) IMP_OVERRIDE;
160  virtual double do_step(const ParticleIndexes &sc,
161  double dt) IMP_OVERRIDE;
162  virtual bool get_is_simulation_particle(ParticleIndex p) const
163  IMP_OVERRIDE;
164 
166 
167  protected:
168  void initialize();
169 
170  virtual void setup_degrees_of_freedom(const ParticleIndexes &ps);
171 
172  //! First part of velocity Verlet (update coordinates and half-step velocity)
173  virtual void propagate_coordinates(const ParticleIndexes &ps,
174  double step_size);
175 
176  //! Second part of velocity Verlet (update velocity)
177  virtual void propagate_velocities(const ParticleIndexes &ps,
178  double step_size);
179 
180  //! Cap a velocity component to the maximum value.
181  inline void cap_velocity_component(Float &vel) {
182  if (vel >= 0.0) {
183  vel = std::min(vel, velocity_cap_);
184  } else {
185  vel = std::max(vel, -velocity_cap_);
186  }
187  }
188 
189  //! Number of degrees of freedom in the system
191 
192  //! Maximum absolute value of a single velocity component
194 };
195 
196 IMPATOM_END_NAMESPACE
197 
198 #endif /* IMPATOM_MOLECULAR_DYNAMICS_H */
A particle with angular velocity.
The base class for simulators.
Definition: Simulator.h:34
ParticleIndex get_particle_index() const
Returns the particle index decorated by this decorator.
Definition: Decorator.h:188
Key< 0 > FloatKey
The type used to identify float attributes in the Particles.
Definition: base_types.h:32
#define IMP_DECORATOR_SETUP_1(Name, FirstArgumentType, first_argument_name)
#define IMP_OBJECT_METHODS(Name)
Define the basic things needed by any Object.
Definition: object_macros.h:25
Model * get_model() const
Returns the Model containing the particle.
Definition: Decorator.h:191
virtual bool get_is_simulation_particle(ParticleIndex p) const =0
Return true if the passed particle is appropriate for the simulation.
A particle with linear (XYZ) velocity.
Base class for all optimizers.
VectorD< 4 > Vector4D
Definition: VectorD.h:399
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
void set_velocity_cap(Float velocity_cap)
Set maximum velocity in A/fs.
void add_attribute(TypeKey attribute_key, ParticleIndex particle, Type value)
add particle atribute with the specied key and initial value
Simple molecular dynamics optimizer.
int degrees_of_freedom_
Number of degrees of freedom in the system.
void set_attribute(TypeKey attribute_key, ParticleIndex particle, Type value)
set the value of particle attribute with the specified key
Simple molecular dynamics optimizer.
#define IMP_DECORATOR_SETUP_0(Name)
Interface to specialized Particle types (e.g. atoms)
Definition: Decorator.h:118
Classes to handle individual model particles. (Note that implementation of inline functions is in int...
void cap_velocity_component(Float &vel)
Cap a velocity component to the maximum value.
Various important macros for implementing decorators.
Float velocity_cap_
Maximum absolute value of a single velocity component.
#define IMP_DECORATOR_METHODS(Name, Parent)
VectorD< 3 > Vector3D
Definition: VectorD.h:395
double Float
Basic floating-point value (could be float, double...)
Definition: types.h:20
#define IMP_USAGE_CHECK(expr, message)
A runtime test for incorrect usage of a class or method.
Definition: check_macros.h:168
bool get_has_attribute(TypeKey attribute_key, ParticleIndex particle) const
return true if particle has attribute with the specified key
Type get_attribute(TypeKey attribute_key, ParticleIndex particle)
get the value of the particle attribute with the specified key
#define IMP_OVERRIDE
Cause a compile error if this method does not override a parent method.