text: #ifndef RPMD_NECKLACE_H
  #define RPMD_NECKLACE_H
  
  #include <cassert>
  #include <cstddef>
  
  #include <armadillo>
  
  namespace parts {
  
  //
  // rpmd_necklace for RPMD; manages memory for cartesian/normal mode
  // positions and momenta and handles the transforms
  //
  
  struct rpmd_necklace {
  
    rpmd_necklace();
  
    inline size_t ndim() const { return m_ndim; };
    inline size_t natoms() const { return m_natoms; };
    inline size_t ndofs() const { return m_ndofs; };
    inline size_t nbeads() const { return m_nbeads; };
    inline size_t beta() const { return m_beta; };
  
  private:
    rpmd_necklace(const rpmd_necklace &);
    rpmd_necklace &operator=(const rpmd_necklace &);
  
    size_t m_ndim;
    size_t m_natoms;
    size_t m_ndofs;
    size_t m_nbeads;
  
  public:
    double m_dt;
  
    double m_beta;
    double m_beta_n;
    double m_omega_n;
  
    void setup(size_t ndim, size_t natom, size_t nbead, double, double, double);
  
    void pos_c2n();
    void pos_n2c();
  
    void mom_c2n();
    void mom_n2c();
  
    void mom2vel();
  
    // layout is bead1, bead2, ..., beadN, where
    // each bead consists of ndofs elements
  
    arma::mat m_pos_cart;
    arma::mat m_pos_nmode;
  
    arma::mat m_mom_cart;
    arma::mat m_mom_nmode;
    arma::mat m_vel_cart;
  
    arma::mat m_frc_cart;
  
    arma::mat m_cart_to_nm;
    arma::cube m_freerp_propagator;
  
    arma::vec m_omega_k;
  
    arma::vec m_mass;
    arma::vec m_sqrt_mass;
  
    constexpr static double engunit =
        1.0; // conversion from internal units to kcal/mol
    constexpr static double kB = 1.0; // Boltzmann constant in internal units
    constexpr static double hbar = 1.0;
  };
  
  // inline const double& rpmd_necklace::lambda(size_t b) const
  //{
  //    assert(b < nbeads());
  //
  //    return m_lambda[b];
  //}
  
  } // namespace parts
  
  #endif // RPMD_NECKLACE_H
  