IonicStepper.h 3.22 KB
Newer Older
1 2
////////////////////////////////////////////////////////////////////////////////
//
Francois Gygi committed
3 4 5 6
// Copyright (c) 2008 The Regents of the University of California
//
// This file is part of Qbox
//
Francois Gygi committed
7 8
// Qbox is distributed under the terms of the GNU General Public License
// as published by the Free Software Foundation, either version 2 of
Francois Gygi committed
9 10 11 12 13 14
// the License, or (at your option) any later version.
// See the file COPYING in the root directory of this distribution
// or <http://www.gnu.org/licenses/>.
//
////////////////////////////////////////////////////////////////////////////////
//
15 16 17 18 19 20 21 22
// IonicStepper.h:
//
////////////////////////////////////////////////////////////////////////////////

#ifndef IONICSTEPPER_H
#define IONICSTEPPER_H

#include "Sample.h"
23
#include "Species.h"
24 25 26 27 28
#include <vector>

class IonicStepper
{
  protected:
29

30
  Sample& s_;
31
  AtomSet& atoms_;
32
  ConstraintSet& constraints_;
Francois Gygi committed
33 34 35
  double dt_;
  int    nsp_;
  int    natoms_;
Francois Gygi committed
36
  // ndofs_ is the total number of degrees of freedom after
Francois Gygi committed
37
  // constraints are considered
38
  int                            ndofs_;
39
  std::vector<int>               na_;  // number of atoms per species na_[nsp_]
40 41 42 43
  std::vector<std::vector< double> >  r0_; // r0_[nsp_][3*na_]
  std::vector<std::vector< double> >  rp_; // rp_[nsp_][3*na_]
  std::vector<std::vector< double> >  rm_; // rm_[nsp_][3*na_]
  std::vector<std::vector< double> >  v0_; // v0_[nsp_][3*na_]
44
  std::vector<double>            pmass_;   // pmass_[nsp_]
45 46

  public:
47 48

  IonicStepper (Sample& s) : s_(s), atoms_(s.atoms),
49
    constraints_(s.constraints), dt_(s.ctrl.dt)
50
  {
51
    ndofs_ = 3 * atoms_.size() - constraints_.ndofs();
52 53
    // if there are more constraints than dofs, set ndofs_ to zero
    if ( ndofs_ < 0 ) ndofs_ = 0;
54 55
    nsp_ = atoms_.nsp();
    na_.resize(nsp_);
56 57 58
    r0_.resize(nsp_);
    rp_.resize(nsp_);
    v0_.resize(nsp_);
59 60 61 62 63
    pmass_.resize(nsp_);
    for ( int is = 0; is < nsp_; is++ )
    {
      const int nais = atoms_.na(is);
      na_[is] = nais;
64 65 66
      r0_[is].resize(3*nais);
      rp_[is].resize(3*nais);
      v0_[is].resize(3*nais);
67 68
      pmass_[is] = atoms_.species_list[is]->mass() * 1822.89;
    }
Francois Gygi committed
69
    natoms_ = atoms_.size();
70 71 72
    // get positions and velocities from atoms_
    get_positions();
    get_velocities();
73
  }
74

75 76
  double r0(int is, int i) const { return r0_[is][i]; }
  double v0(int is, int i) const { return v0_[is][i]; }
77 78 79
  const std::vector<std::vector<double> >& r0(void) const { return r0_; }
  const std::vector<std::vector<double> >& v0(void) const { return v0_; }
  const std::vector<double>& pmass(void) const { return pmass_; }
80 81 82 83
  void get_positions(void) { atoms_.get_positions(r0_); }
  void get_velocities(void) { atoms_.get_velocities(v0_); }
  void set_positions(void) { atoms_.set_positions(r0_); }
  void set_velocities(void) { atoms_.set_velocities(v0_); }
84

85 86 87 88
  void setup_constraints(void)
  {
    constraints_.setup(atoms_);
  }
89
  virtual void compute_r(double e0,
90
    const std::vector<std::vector< double> >& f0) = 0;
91
  virtual void compute_v(double e0,
92
    const std::vector<std::vector< double> >& f0) = 0;
93
  virtual void reset(void) {}
94
  virtual double ekin(void) const { return 0.0; }
95
  virtual double ekin_stepper(void) const { return 0.0; }
96
  virtual double temp(void) const { return 0.0; }
97

98
  virtual ~IonicStepper() {}
99

100 101
};
#endif