colvaratoms.h

// -*- c++ -*-
 
// This file is part of the Collective Variables module (Colvars).
// The original version of Colvars and its updates are located at:
// https://github.com/Colvars/colvars
// Please update all Colvars source files before making any changes.
// If you wish to distribute your changes, please submit them to the
// Colvars repository at GitHub.
 
#ifndef COLVARATOMS_H
#define COLVARATOMS_H
 
#include "colvarmodule.h"
#include "colvarproxy.h"
#include "colvarparse.h"
#include "colvardeps.h"
 
template <typename T1, typename T2>
struct rotation_derivative;
 
 
 
class colvarmodule::atom {
 
protected:
 
  int index;
 
public:
 
  int             id;
 
  cvm::real       mass;
 
  cvm::real       charge;
 
  cvm::atom_pos   pos;
 
  cvm::rvector    vel;
 
  cvm::rvector    total_force;
 
  cvm::rvector   grad;
 
  atom();
 
  atom(int atom_number);
 
  atom(cvm::residue_id const &residue,
       std::string const     &atom_name,
       std::string const     &segment_id);
 
  atom(atom const &a);
 
  ~atom();
 
  atom & operator = (atom const &a);
 
  inline void reset_data()
  {
    pos = cvm::atom_pos(0.0);
    vel = grad = total_force = cvm::rvector(0.0);
  }
 
  inline void update_mass()
  {
    colvarproxy *p = cvm::proxy;
    mass = p->get_atom_mass(index);
  }
 
  inline void update_charge()
  {
    colvarproxy *p = cvm::proxy;
    charge = p->get_atom_charge(index);
  }
 
  inline void read_position()
  {
    pos = (cvm::proxy)->get_atom_position(index);
  }
 
  inline void read_velocity()
  {
    vel = (cvm::proxy)->get_atom_velocity(index);
  }
 
  inline void read_total_force()
  {
    total_force = (cvm::proxy)->get_atom_total_force(index);
  }
 
  inline void apply_force(cvm::rvector const &new_force) const
  {
    (cvm::proxy)->apply_atom_force(index, new_force);
  }
};
 
 
 
class colvarmodule::atom_group
  : public colvarparse, public colvardeps
{
public:
 
 
  atom_group();
 
  atom_group(char const *key);
 
  atom_group(std::vector<cvm::atom> const &atoms_in);
 
  ~atom_group() override;
 
  std::string name;
 
  // TODO Make this field part of the data structures that link a group to a CVC
  std::string key;
 
  int init();
 
  int init_dependencies() override;
 
  int setup();
 
  int parse(std::string const &conf);
 
  int add_atom_numbers(std::string const &numbers_conf);
  int add_atoms_of_group(atom_group const * ag);
  int add_index_group(std::string const &index_group_name, bool silent = false);
  int add_atom_numbers_range(std::string const &range_conf);
  int add_atom_name_residue_range(std::string const &psf_segid,
                                  std::string const &range_conf);
  int parse_fitting_options(std::string const &group_conf);
 
  int add_atom(cvm::atom const &a);
 
  int add_atom_id(int aid);
 
  int remove_atom(cvm::atom_iter ai);
 
  int set_dummy();
 
  int set_dummy_pos(cvm::atom_pos const &pos);
 
  void print_properties(std::string const &colvar_name, int i, int j);
 
  static std::vector<feature *> ag_features;
 
  const std::vector<feature *> &features() const override { return ag_features; }
 
  std::vector<feature *> &modify_features() override { return ag_features; }
 
  static void delete_features()
  {
    for (size_t i = 0; i < ag_features.size(); i++) {
      delete ag_features[i];
    }
    ag_features.clear();
  }
 
protected:
 
  std::vector<cvm::atom> atoms;
 
  std::vector<int> atoms_ids;
 
  std::vector<int> sorted_atoms_ids;
 
  std::vector<int> sorted_atoms_ids_map;
 
  cvm::atom_pos dummy_atom_pos;
 
  int index;
 
  std::vector<cvm::rvector> group_forces;
 
public:
 
  class group_force_object {
  public:
    group_force_object(cvm::atom_group* ag);
    ~group_force_object();
    void add_atom_force(size_t i, const cvm::rvector& force);
  private:
    cvm::atom_group* m_ag;
    cvm::atom_group* m_group_for_fit;
    bool m_has_fitting_force;
    void apply_force_with_fitting_group();
  };
 
  group_force_object get_group_force_object();
 
  inline cvm::atom & operator [] (size_t const i)
  {
    return atoms[i];
  }
 
  inline cvm::atom const & operator [] (size_t const i) const
  {
    return atoms[i];
  }
 
  inline cvm::atom_iter begin()
  {
    return atoms.begin();
  }
 
  inline cvm::atom_const_iter begin() const
  {
    return atoms.begin();
  }
 
  inline cvm::atom_iter end()
  {
    return atoms.end();
  }
 
  inline cvm::atom_const_iter end() const
  {
    return atoms.end();
  }
 
  inline size_t size() const
  {
    return atoms.size();
  }
 
  bool b_dummy;
 
  inline std::vector<int> const &ids() const
  {
    return atoms_ids;
  }
 
  std::string const print_atom_ids() const;
 
  int create_sorted_ids();
 
  inline std::vector<int> const &sorted_ids() const
  {
    return sorted_atoms_ids;
  }
 
  inline std::vector<int> const &sorted_ids_map() const
  {
    return sorted_atoms_ids_map;
  }
 
  static int overlap(const atom_group &g1, const atom_group &g2);
 
  cvm::rotation rot;
 
  rotation_derivative<cvm::atom, cvm::atom_pos>* rot_deriv;
 
  bool b_user_defined_fit;
 
  std::vector<cvm::atom_pos> ref_pos;
 
  cvm::atom_pos              ref_pos_cog;
 
  atom_group                *fitting_group;
 
  cvm::real total_mass;
 
  void update_total_mass();
 
  cvm::real total_charge;
 
  void update_total_charge();
 
  bool noforce;
 
  void read_positions();
 
  void calc_apply_roto_translation();
 
  void setup_rotation_derivative();
 
  void center_ref_pos();
 
  void apply_translation(cvm::rvector const &t);
 
  void read_velocities();
 
  void read_total_forces();
 
  inline void reset_atoms_data()
  {
    for (cvm::atom_iter ai = atoms.begin(); ai != atoms.end(); ai++)
      ai->reset_data();
    if (fitting_group)
      fitting_group->reset_atoms_data();
  }
 
  int calc_required_properties();
 
  std::vector<cvm::atom_pos> positions() const;
 
  int calc_center_of_geometry();
 
private:
 
  cvm::atom_pos cog;
 
  cvm::atom_pos cog_orig;
 
  std::vector<cvm::atom_pos> pos_unrotated;
 
public:
 
  inline cvm::atom_pos center_of_geometry() const
  {
    return cog;
  }
 
  int calc_center_of_mass();
 
private:
 
  cvm::atom_pos com;
 
  // TODO for scalable calculations of more complex variables (e.g. rotation),
  // use a colvarvalue of vectors to hold the entire derivative
  cvm::rvector scalar_com_gradient;
 
public:
 
  inline cvm::atom_pos center_of_mass() const
  {
    return com;
  }
 
  inline cvm::rvector center_of_mass_scalar_gradient() const
  {
    return scalar_com_gradient;
  }
 
  std::vector<cvm::atom_pos> positions_shifted(cvm::rvector const &shift) const;
 
  std::vector<cvm::rvector> velocities() const;
 
  int calc_dipole(cvm::atom_pos const &dipole_center);
 
private:
 
  cvm::rvector dip;
 
public:
 
  inline cvm::rvector dipole() const
  {
    return dip;
  }
 
  std::vector<cvm::rvector> total_forces() const;
 
  cvm::rvector total_force() const;
 
 
  void set_weighted_gradient(cvm::rvector const &grad);
 
  void calc_fit_gradients();
 
  template <bool B_ag_center, bool B_ag_rotate,
            typename main_force_accessor_T, typename fitting_force_accessor_T>
  void calc_fit_forces_impl(
    main_force_accessor_T accessor_main,
    fitting_force_accessor_T accessor_fitting) const;
 
  template <typename main_force_accessor_T, typename fitting_force_accessor_T>
  void calc_fit_forces(
    main_force_accessor_T accessor_main,
    fitting_force_accessor_T accessor_fitting) const;
 
  std::vector<cvm::atom_pos> fit_gradients;
 
  void apply_colvar_force(cvm::real const &force);
 
  void apply_force(cvm::rvector const &force);
 
  void do_feature_side_effects(int id) override;
};
 
 
#endif