IMP logo
IMP Reference Guide  2.17.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-2022 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_velocities_key(), pi, v.get_coordinates());
27  }
28  static FloatsKey get_velocities_key() {
29  static const FloatsKey key("linvel");
30  return key;
31  }
32 
33 public:
34  static bool get_is_setup(Model *m, ParticleIndex pi) {
35  return m->get_has_attribute(get_velocities_key(), pi);
36  }
37 
41 
42  void set_velocity(const algebra::Vector3D &v) {
43  Model *m = get_model();
45  m->set_attribute(get_velocities_key(), pi, v.get_coordinates());
46  }
47 
48  algebra::Vector3D get_velocity() const {
49  Model *m = get_model();
51  return algebra::Vector3D(m->get_attribute(get_velocities_key(), pi));
52  }
53 };
54 
55 //! A particle with angular velocity
56 /** Typically this is used for RigidBody particles in combination
57  with the MolecularDynamics optimizer. The velocity is stored as
58  a quaternion.
59  */
60 class IMPATOMEXPORT AngularVelocity : public Decorator {
61  static void do_setup_particle(Model *m, ParticleIndex pi,
62  const algebra::Vector4D v = algebra::Vector4D(0, 0, 0, 0)) {
63  m->add_attribute(get_velocities_key(), pi, v.get_coordinates());
64  }
65  static FloatsKey get_velocities_key() {
66  static const FloatsKey key("angvel");
67  return key;
68  }
69 
70 public:
71 
72  static bool get_is_setup(Model *m, ParticleIndex pi) {
73  return m->get_has_attribute(get_velocities_key(), pi);
74  }
75 
79 
80  void set_velocity(const algebra::Vector4D &v) {
81  Model *m = get_model();
83  m->set_attribute(get_velocities_key(), pi, v.get_coordinates());
84  }
85 
86  algebra::Vector4D get_velocity() const {
87  Model *m = get_model();
89  return algebra::Vector4D(m->get_attribute(get_velocities_key(), pi));
90  }
91 };
92 
93 //! Simple molecular dynamics optimizer.
94 /** The particles to be optimized must have optimizable x,y,z attributes
95  and a non-optimizable mass attribute; this optimizer assumes the score
96  to be energy in kcal/mol, the xyz coordinates to be in angstroms, and
97  the mass to be in AMU (g/mol).
98 
99  \note RigidBody particles are not handled properly.
100 
101  Particles without optimized x,y,z and nonoptimized mass are skipped.
102  \see VelocityScalingOptimizerState
103  \see LangevinThermostatOptimizerState
104  \see BerendsenThermostatOptimizerState
105  \see RemoveRigidMotionOptimizerState
106  */
107 class IMPATOMEXPORT MolecularDynamics : public Simulator {
108  public:
109  /** Score based on the provided model */
111 
112  //! Return the current kinetic energy of the system, in kcal/mol
113  virtual Float get_kinetic_energy() const;
114 
115  //! Return the current kinetic temperature of the system
116  /** \param[in] ekinetic kinetic energy, e.g. from get_kinetic_energy()
117  */
118  Float get_kinetic_temperature(Float ekinetic) const;
119 
120  //! Set maximum velocity in A/fs
121  /** At each dynamics time step, the absolute value of each velocity
122  component is capped at this value. This prevents spurious strong forces
123  (occasionally encountered with frustrated conformations) from causing
124  large oscillations in the system.
125  By default, velocities are not capped.
126 
127  \note The actual velocities that are capped are the half-step velocities
128  in the velocity Verlet algorithm.
129  */
130  void set_velocity_cap(Float velocity_cap) { velocity_cap_ = velocity_cap; }
131 
132  //! Assign velocities representative of the given temperature
133  virtual void assign_velocities(Float temperature);
134  virtual void setup(const ParticleIndexes &ps) override;
135  virtual double do_step(const ParticleIndexes &sc,
136  double dt) override;
137  virtual bool get_is_simulation_particle(ParticleIndex p) const
138  override;
139 
141 
142  protected:
143  void initialize();
144 
145  virtual void setup_degrees_of_freedom(const ParticleIndexes &ps);
146 
147  //! First part of velocity Verlet (update coordinates and half-step velocity)
148  virtual void propagate_coordinates(const ParticleIndexes &ps,
149  double step_size);
150 
151  //! Second part of velocity Verlet (update velocity)
152  virtual void propagate_velocities(const ParticleIndexes &ps,
153  double step_size);
154 
155  //! Cap a velocity component to the maximum value.
156  inline void cap_velocity_component(Float &vel) {
157  if (vel >= 0.0) {
158  vel = std::min(vel, velocity_cap_);
159  } else {
160  vel = std::max(vel, -velocity_cap_);
161  }
162  }
163 
164  //! Number of degrees of freedom in the system
166 
167  //! Maximum absolute value of a single velocity component
169 };
170 
171 IMPATOM_END_NAMESPACE
172 
173 #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:190
#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:193
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:425
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
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
A base class for Keys.
Definition: Key.h:44
Base class for "simulators", such as molecular dynamics.
#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.
Macros for maintaining molecular hierarchies.
Float velocity_cap_
Maximum absolute value of a single velocity component.
#define IMP_DECORATOR_METHODS(Name, Parent)
VectorD< 3 > Vector3D
Definition: VectorD.h:421
double Float
Basic floating-point value (could be float, double...)
Definition: types.h:19
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