Update Colvars to version 2018-05-15
[namd.git] / colvars / src / colvaratoms.h
1 // -*- c++ -*-
2
3 // This file is part of the Collective Variables module (Colvars).
4 // The original version of Colvars and its updates are located at:
5 // https://github.com/colvars/colvars
6 // Please update all Colvars source files before making any changes.
7 // If you wish to distribute your changes, please submit them to the
8 // Colvars repository at GitHub.
9
10 #ifndef COLVARATOMS_H
11 #define COLVARATOMS_H
12
13 #include "colvarmodule.h"
14 #include "colvarproxy.h"
15 #include "colvarparse.h"
16 #include "colvardeps.h"
17
18
19 /// \brief Stores numeric id, mass and all mutable data for an atom,
20 /// mostly used by a \link cvc \endlink
21 ///
22 /// This class may be used to keep atomic data such as id, mass,
23 /// position and collective variable derivatives) altogether.
24 /// There may be multiple instances with identical
25 /// numeric id, all acting independently: forces communicated through
26 /// these instances will be summed together.
27
28 class colvarmodule::atom {
29
30 protected:
31
32   /// Index in the colvarproxy arrays (\b NOT in the global topology!)
33   int index;
34
35 public:
36
37   /// Identifier for the MD program (0-based)
38   int             id;
39
40   /// Mass
41   cvm::real       mass;
42
43   /// Charge
44   cvm::real       charge;
45
46   /// \brief Current position (copied from the program, can be
47   /// modified if necessary)
48   cvm::atom_pos   pos;
49
50   /// \brief Current velocity (copied from the program, can be
51   /// modified if necessary)
52   cvm::rvector    vel;
53
54   /// \brief System force at the previous step (copied from the
55   /// program, can be modified if necessary)
56   cvm::rvector    total_force;
57
58   /// \brief Gradient of a scalar collective variable with respect
59   /// to this atom
60   ///
61   /// This can only handle a scalar collective variable (i.e. when
62   /// the \link colvarvalue::real_value \endlink member is used
63   /// from the \link colvarvalue \endlink class), which is also the
64   /// most frequent case. For more complex types of \link
65   /// colvarvalue \endlink objects, atomic gradients should be
66   /// defined within the specific \link cvc \endlink
67   /// implementation
68   cvm::rvector   grad;
69
70   /// \brief Default constructor (sets index and id both to -1)
71   atom();
72
73   /// \brief Initialize an atom for collective variable calculation
74   /// and get its internal identifier \param atom_number Atom index in
75   /// the system topology (1-based)
76   atom(int atom_number);
77
78   /// \brief Initialize an atom for collective variable calculation
79   /// and get its internal identifier \param residue Residue number
80   /// \param atom_name Name of the atom in the residue \param
81   /// segment_id For PSF topologies, the segment identifier; for other
82   /// type of topologies, may not be required
83   atom(cvm::residue_id const &residue,
84        std::string const     &atom_name,
85        std::string const     &segment_id);
86
87   /// Copy constructor
88   atom(atom const &a);
89
90   /// Destructor
91   ~atom();
92
93   /// Set mutable data (everything except id and mass) to zero
94   inline void reset_data()
95   {
96     pos = cvm::atom_pos(0.0);
97     vel = grad = total_force = cvm::rvector(0.0);
98   }
99
100   /// Get the latest value of the mass
101   inline void update_mass()
102   {
103     mass = (cvm::proxy)->get_atom_mass(index);
104   }
105
106   /// Get the latest value of the charge
107   inline void update_charge()
108   {
109     charge = (cvm::proxy)->get_atom_charge(index);
110   }
111
112   /// Get the current position
113   inline void read_position()
114   {
115     pos = (cvm::proxy)->get_atom_position(index);
116   }
117
118   /// Get the current velocity
119   inline void read_velocity()
120   {
121     vel = (cvm::proxy)->get_atom_velocity(index);
122   }
123
124   /// Get the total force
125   inline void read_total_force()
126   {
127     total_force = (cvm::proxy)->get_atom_total_force(index);
128   }
129
130   /// \brief Apply a force to the atom
131   ///
132   /// Note: the force is not applied instantly, but will be used later
133   /// by the MD integrator (the colvars module does not integrate
134   /// equations of motion.
135   ///
136   /// Multiple calls to this function by either the same
137   /// \link atom \endlink object or different objects with identical
138   /// \link id \endlink will all be added together.
139   inline void apply_force(cvm::rvector const &new_force) const
140   {
141     (cvm::proxy)->apply_atom_force(index, new_force);
142   }
143 };
144
145
146
147 /// \brief Group of \link atom \endlink objects, mostly used by a
148 /// \link cvc \endlink object to gather all atomic data
149 class colvarmodule::atom_group
150   : public colvarparse, public colvardeps
151 {
152 public:
153
154
155   /// \brief Default constructor
156   atom_group();
157
158   /// \brief Create a group object, assign a name to it
159   atom_group(char const *key);
160
161   /// \brief Initialize the group after a (temporary) vector of atoms
162   atom_group(std::vector<cvm::atom> const &atoms_in);
163
164   /// \brief Destructor
165   ~atom_group();
166
167   /// \brief Optional name to reuse properties of this in other groups
168   std::string name;
169
170   /// \brief Keyword used to define the group
171   // TODO Make this field part of the data structures that link a group to a CVC
172   std::string key;
173
174   /// \brief Set default values for common flags
175   int init();
176
177   /// \brief Update data required to calculate cvc's
178   int setup();
179
180   /// \brief Initialize the group by looking up its configuration
181   /// string in conf and parsing it
182   int parse(std::string const &conf);
183
184   int add_atom_numbers(std::string const &numbers_conf);
185   int add_atoms_of_group(atom_group const * ag);
186   int add_index_group(std::string const &index_group_name);
187   int add_atom_numbers_range(std::string const &range_conf);
188   int add_atom_name_residue_range(std::string const &psf_segid,
189                                   std::string const &range_conf);
190   int parse_fitting_options(std::string const &group_conf);
191
192   /// \brief Add an atom object to this group
193   int add_atom(cvm::atom const &a);
194
195   /// \brief Add an atom ID to this group (the actual atomicdata will be not be handled by the group)
196   int add_atom_id(int aid);
197
198   /// \brief Remove an atom object from this group
199   int remove_atom(cvm::atom_iter ai);
200
201   /// \brief Re-initialize the total mass of a group.
202   /// This is needed in case the hosting MD code has an option to
203   /// change atom masses after their initialization.
204   void reset_mass(std::string &name, int i, int j);
205
206   /// \brief Implementation of the feature list for atom group
207   static std::vector<feature *> ag_features;
208
209   /// \brief Implementation of the feature list accessor for atom group
210   virtual const std::vector<feature *> &features()
211   {
212     return ag_features;
213   }
214   virtual std::vector<feature *> &modify_features()
215   {
216     return ag_features;
217   }
218   static void delete_features() {
219     for (size_t i=0; i < ag_features.size(); i++) {
220       delete ag_features[i];
221     }
222     ag_features.clear();
223   }
224
225 protected:
226
227   /// \brief Array of atom objects
228   std::vector<cvm::atom> atoms;
229
230   /// \brief Internal atom IDs for host code
231   std::vector<int> atoms_ids;
232
233   /// Sorted list of internal atom IDs (populated on-demand by
234   /// create_sorted_ids); used to read coordinate files
235   std::vector<int> sorted_atoms_ids;
236
237   /// Map entries of sorted_atoms_ids onto the original positions in the group
238   std::vector<int> sorted_atoms_ids_map;
239
240   /// \brief Dummy atom position
241   cvm::atom_pos dummy_atom_pos;
242
243   /// \brief Index in the colvarproxy arrays (if the group is scalable)
244   int index;
245
246 public:
247
248   inline cvm::atom & operator [] (size_t const i)
249   {
250     return atoms[i];
251   }
252
253   inline cvm::atom const & operator [] (size_t const i) const
254   {
255     return atoms[i];
256   }
257
258   inline cvm::atom_iter begin()
259   {
260     return atoms.begin();
261   }
262
263   inline cvm::atom_const_iter begin() const
264   {
265     return atoms.begin();
266   }
267
268   inline cvm::atom_iter end()
269   {
270     return atoms.end();
271   }
272
273   inline cvm::atom_const_iter end() const
274   {
275     return atoms.end();
276   }
277
278   inline size_t size() const
279   {
280     return atoms.size();
281   }
282
283   /// \brief If this option is on, this group merely acts as a wrapper
284   /// for a fixed position; any calls to atoms within or to
285   /// functions that return disaggregated data will fail
286   bool b_dummy;
287
288   /// Internal atom IDs (populated during initialization)
289   inline std::vector<int> const &ids() const
290   {
291     return atoms_ids;
292   }
293
294   std::string const print_atom_ids() const;
295
296   /// Allocates and populates sorted_ids and sorted_ids_map
297   int create_sorted_ids();
298
299   /// Sorted internal atom IDs (populated on-demand by create_sorted_ids);
300   /// used to read coordinate files
301   inline std::vector<int> const &sorted_ids() const
302   {
303     return sorted_atoms_ids;
304   }
305
306   /// Map entries of sorted_atoms_ids onto the original positions in the group
307   inline std::vector<int> const &sorted_ids_map() const
308   {
309     return sorted_atoms_ids_map;
310   }
311
312   /// Detect whether two groups share atoms
313   /// If yes, returns 1-based number of a common atom; else, returns 0
314   static int overlap(const atom_group &g1, const atom_group &g2);
315
316   /// \brief When updating atomic coordinates, translate them to align with the
317   /// center of mass of the reference coordinates
318   bool b_center;
319
320   /// \brief When updating atom coordinates (and after
321   /// centering them if b_center is set), rotate the group to
322   /// align with the reference coordinates.
323   ///
324   /// Note: gradients will be calculated in the rotated frame: when
325   /// forces will be applied, they will rotated back to the original
326   /// frame
327   bool b_rotate;
328   /// The rotation calculated automatically if b_rotate is defined
329   cvm::rotation rot;
330
331   /// \brief Indicates that the user has explicitly set centerReference or
332   /// rotateReference, and the corresponding reference:
333   /// cvc's (eg rmsd, eigenvector) will not override the user's choice
334   bool b_user_defined_fit;
335
336   /// \brief use reference coordinates for b_center or b_rotate
337   std::vector<cvm::atom_pos> ref_pos;
338
339   /// \brief Center of geometry of the reference coordinates; regardless
340   /// of whether b_center is true, ref_pos is centered to zero at
341   /// initialization, and ref_pos_cog serves to center the positions
342   cvm::atom_pos              ref_pos_cog;
343
344   /// \brief If b_center or b_rotate is true, use this group to
345   /// define the transformation (default: this group itself)
346   atom_group                *fitting_group;
347
348   /// Total mass of the atom group
349   cvm::real total_mass;
350   void update_total_mass();
351
352   /// Total charge of the atom group
353   cvm::real total_charge;
354   void update_total_charge();
355
356   /// \brief Don't apply any force on this group (use its coordinates
357   /// only to calculate a colvar)
358   bool        noforce;
359
360   /// \brief Get the current positions
361   void read_positions();
362
363   /// \brief (Re)calculate the optimal roto-translation
364   void calc_apply_roto_translation();
365
366   /// \brief Save aside the center of geometry of the reference positions,
367   /// then subtract it from them
368   ///
369   /// In this way it will be possible to use ref_pos also for the
370   /// rotational fit.
371   /// This is called either by atom_group::parse or by CVCs that assign
372   /// reference positions (eg. RMSD, eigenvector).
373   void center_ref_pos();
374
375   /// \brief Move all positions
376   void apply_translation(cvm::rvector const &t);
377
378   /// \brief Get the current velocities; this must be called always
379   /// *after* read_positions(); if b_rotate is defined, the same
380   /// rotation applied to the coordinates will be used
381   void read_velocities();
382
383   /// \brief Get the current total_forces; this must be called always
384   /// *after* read_positions(); if b_rotate is defined, the same
385   /// rotation applied to the coordinates will be used
386   void read_total_forces();
387
388   /// Call reset_data() for each atom
389   inline void reset_atoms_data()
390   {
391     for (cvm::atom_iter ai = atoms.begin(); ai != atoms.end(); ai++)
392       ai->reset_data();
393     if (fitting_group)
394       fitting_group->reset_atoms_data();
395   }
396
397   /// \brief Recompute all mutable quantities that are required to compute CVCs
398   int calc_required_properties();
399
400   /// \brief Return a copy of the current atom positions
401   std::vector<cvm::atom_pos> positions() const;
402
403   /// \brief Calculate the center of geometry of the atomic positions, assuming
404   /// that they are already pbc-wrapped
405   int calc_center_of_geometry();
406
407 private:
408
409   /// \brief Center of geometry
410   cvm::atom_pos cog;
411
412   /// \brief Center of geometry before any fitting
413   cvm::atom_pos cog_orig;
414
415 public:
416
417   /// \brief Return the center of geometry of the atomic positions
418   inline cvm::atom_pos center_of_geometry() const
419   {
420     return cog;
421   }
422
423   /// \brief Calculate the center of mass of the atomic positions, assuming that
424   /// they are already pbc-wrapped
425   int calc_center_of_mass();
426 private:
427   /// \brief Center of mass
428   cvm::atom_pos com;
429   /// \brief The derivative of a scalar variable with respect to the COM
430   // TODO for scalable calculations of more complex variables (e.g. rotation),
431   // use a colvarvalue of vectors to hold the entire derivative
432   cvm::rvector scalar_com_gradient;
433 public:
434   /// \brief Return the center of mass of the atomic positions
435   inline cvm::atom_pos center_of_mass() const
436   {
437     return com;
438   }
439
440   /// \brief Return a copy of the current atom positions, shifted by a constant vector
441   std::vector<cvm::atom_pos> positions_shifted(cvm::rvector const &shift) const;
442
443   /// \brief Return a copy of the current atom velocities
444   std::vector<cvm::rvector> velocities() const;
445
446   ///\brief Calculate the dipole of the atom group around the specified center
447   int calc_dipole(cvm::atom_pos const &com);
448 private:
449   cvm::rvector dip;
450 public:
451   ///\brief Return the (previously calculated) dipole of the atom group
452   inline cvm::rvector dipole() const
453   {
454     return dip;
455   }
456
457   /// \brief Return a copy of the total forces
458   std::vector<cvm::rvector> total_forces() const;
459
460   /// \brief Return a copy of the aggregated total force on the group
461   cvm::rvector total_force() const;
462
463
464   /// \brief Shorthand: save the specified gradient on each atom,
465   /// weighting with the atom mass (mostly used in combination with
466   /// \link center_of_mass() \endlink)
467   void set_weighted_gradient(cvm::rvector const &grad);
468
469   /// \brief Calculate the derivatives of the fitting transformation
470   void calc_fit_gradients();
471
472   /// \brief Derivatives of the fitting transformation
473   std::vector<cvm::atom_pos> fit_gradients;
474
475   /// \brief Used by a (scalar) colvar to apply its force on its \link
476   /// atom_group \endlink members
477   ///
478   /// The (scalar) force is multiplied by the colvar gradient for each
479   /// atom; this should be used when a colvar with scalar \link
480   /// colvarvalue \endlink type is used (this is the most frequent
481   /// case: for colvars with a non-scalar type, the most convenient
482   /// solution is to sum together the Cartesian forces from all the
483   /// colvar components, and use apply_force() or apply_forces()).  If
484   /// the group is being rotated to a reference frame (e.g. to express
485   /// the colvar independently from the solute rotation), the
486   /// gradients are temporarily rotated to the original frame.
487   void apply_colvar_force(cvm::real const &force);
488
489   /// \brief Apply a force "to the center of mass", i.e. the force is
490   /// distributed on each atom according to its mass
491   ///
492   /// If the group is being rotated to a reference frame (e.g. to
493   /// express the colvar independently from the solute rotation), the
494   /// force is rotated back to the original frame.  Colvar gradients
495   /// are not used, either because they were not defined (e.g because
496   /// the colvar has not a scalar value) or the biases require to
497   /// micromanage the force.
498   /// This function will be phased out eventually, in favor of
499   /// apply_colvar_force() once that is implemented for non-scalar values
500   void apply_force(cvm::rvector const &force);
501
502   /// Implements possible actions to be carried out
503   /// when a given feature is enabled
504   /// This overloads the base function in colvardeps
505   void do_feature_side_effects(int id);
506 };
507
508
509 #endif