IonicStepper.h 3.03 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
// IonicStepper.h:
//
////////////////////////////////////////////////////////////////////////////////
18
// $Id: IonicStepper.h,v 1.15 2010-04-16 21:43:11 fgygi Exp $
19 20 21 22 23

#ifndef IONICSTEPPER_H
#define IONICSTEPPER_H

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

class IonicStepper
{
  protected:
30

31
  Sample& s_;
32
  AtomSet& atoms_;
33
  ConstraintSet& constraints_;
Francois Gygi committed
34 35 36
  double dt_;
  int    nsp_;
  int    natoms_;
Francois Gygi committed
37
  // ndofs_ is the total number of degrees of freedom after
Francois Gygi committed
38
  // constraints are considered
39
  int                            ndofs_;
40
  std::vector<int>               na_;  // number of atoms per species na_[nsp_]
41 42 43 44
  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_]
45
  std::vector<double>            pmass_;   // pmass_[nsp_]
46 47

  public:
48 49

  IonicStepper (Sample& s) : s_(s), atoms_(s.atoms),
50
    constraints_(s.constraints), dt_(s.ctrl.dt)
51
  {
52
    ndofs_ = 3 * atoms_.size() - constraints_.ndofs();
53 54
    // if there are more constraints than dofs, set ndofs_ to zero
    if ( ndofs_ < 0 ) ndofs_ = 0;
55 56
    nsp_ = atoms_.nsp();
    na_.resize(nsp_);
57 58 59
    r0_.resize(nsp_);
    rp_.resize(nsp_);
    v0_.resize(nsp_);
60 61 62 63 64
    pmass_.resize(nsp_);
    for ( int is = 0; is < nsp_; is++ )
    {
      const int nais = atoms_.na(is);
      na_[is] = nais;
65 66 67
      r0_[is].resize(3*nais);
      rp_[is].resize(3*nais);
      v0_[is].resize(3*nais);
68 69
      pmass_[is] = atoms_.species_list[is]->mass() * 1822.89;
    }
Francois Gygi committed
70
    natoms_ = atoms_.size();
71
    atoms_.get_positions(r0_);
72
    atoms_.get_velocities(v0_);
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 84
  void setup_constraints(void)
  {
    constraints_.setup(atoms_);
  }
85
  virtual void compute_r(double e0,
86
    const std::vector<std::vector< double> >& f0) = 0;
87
  virtual void compute_v(double e0,
88
    const std::vector<std::vector< double> >& f0) = 0;
89
  virtual void reset(void) {}
90
  virtual double ekin(void) const { return 0.0; }
91
  virtual double ekin_stepper(void) const { return 0.0; }
92
  virtual double temp(void) const { return 0.0; }
93

94
  virtual ~IonicStepper() {}
95

96 97
};
#endif