IMP logo
IMP Reference Guide  2.20.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 simulator.
94 /**
95  The particles to be optimized must have optimizable x,y,z attributes
96  and a non-optimizable mass attribute; this optimizer assumes the score
97  to be energy in kcal/mol, the xyz coordinates to be in angstroms, and
98  the mass to be in AMU (g/mol).
99 
100  \note RigidBody particles are not handled properly.
101 
102  Particles without optimized x,y,z and nonoptimized mass are skipped.
103  \see VelocityScalingOptimizerState
104  \see LangevinThermostatOptimizerState
105  \see BerendsenThermostatOptimizerState
106  \see RemoveRigidMotionOptimizerState
107  */
108 class IMPATOMEXPORT MolecularDynamics : public Simulator {
109  public:
110  /** Score based on the provided model */
112 
113  //! Return the current kinetic energy of the system, in kcal/mol
114  virtual Float get_kinetic_energy() const;
115 
116  //! Return the current kinetic temperature of the system
117  /** \param[in] ekinetic kinetic energy, e.g. from get_kinetic_energy()
118  */
119  Float get_kinetic_temperature(Float ekinetic) const;
120 
121  //! Set maximum velocity in A/fs
122  /** At each dynamics time step, the absolute value of each velocity
123  component is capped at this value. This prevents spurious strong forces
124  (occasionally encountered with frustrated conformations) from causing
125  large oscillations in the system.
126  By default, velocities are not capped.
127 
128  \note The actual velocities that are capped are the half-step velocities
129  in the velocity Verlet algorithm.
130  */
131  void set_velocity_cap(Float velocity_cap) { velocity_cap_ = velocity_cap; }
132 
133  //! Assign velocities representative of the given temperature
134  virtual void assign_velocities(Float temperature);
135  virtual void setup(const ParticleIndexes &ps) override;
136  virtual double do_step(const ParticleIndexes &sc,
137  double dt) override;
138  virtual bool get_is_simulation_particle(ParticleIndex p) const
139  override;
140 
142 
143  protected:
144  void initialize();
145 
146  virtual void setup_degrees_of_freedom(const ParticleIndexes &ps);
147 
148  //! First part of velocity Verlet (update coordinates and half-step velocity)
149  virtual void propagate_coordinates(const ParticleIndexes &ps,
150  double step_size);
151 
152  //! Second part of velocity Verlet (update velocity)
153  virtual void propagate_velocities(const ParticleIndexes &ps,
154  double step_size);
155 
156  //! Cap a velocity component to the maximum value.
157  inline void cap_velocity_component(Float &vel) {
158  if (vel >= 0.0) {
159  vel = std::min(vel, velocity_cap_);
160  } else {
161  vel = std::max(vel, -velocity_cap_);
162  }
163  }
164 
165  //! Number of degrees of freedom in the system
167 
168  //! Maximum absolute value of a single velocity component
170 };
171 
172 IMPATOM_END_NAMESPACE
173 
174 #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:211
#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:214
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:429
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
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 attribute with the specified key and initial value
Simple molecular dynamics simulator.
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:45
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:119
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:425
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