1 //////////////////////////////////////////////////////////////////////////////// 2 // 3 // Copyright (c) 2008 The Regents of the University of California 4 // 5 // This file is part of Qbox 6 // 7 // Qbox is distributed under the terms of the GNU General Public License 8 // as published by the Free Software Foundation, either version 2 of 9 // the License, or (at your option) any later version. 10 // See the file COPYING in the root directory of this distribution 11 // or <http://www.gnu.org/licenses/>. 12 // 13 //////////////////////////////////////////////////////////////////////////////// 14 // 15 // IonicStepper.h: 16 // 17 //////////////////////////////////////////////////////////////////////////////// 18 19 #ifndef IONICSTEPPER_H 20 #define IONICSTEPPER_H 21 22 #include "Sample.h" 23 #include "Species.h" 24 #include <vector> 25 26 class IonicStepper 27 { 28 protected: 29 30 Sample& s_; 31 AtomSet& atoms_; 32 ConstraintSet& constraints_; 33 double dt_; 34 int nsp_; 35 int natoms_; 36 // ndofs_ is the total number of degrees of freedom after 37 // constraints are considered 38 int ndofs_; 39 std::vector<int> na_; // number of atoms per species na_[nsp_] 40 std::vector<std::vector< double> > r0_; // r0_[nsp_][3*na_] 41 std::vector<std::vector< double> > rp_; // rp_[nsp_][3*na_] 42 std::vector<std::vector< double> > rm_; // rm_[nsp_][3*na_] 43 std::vector<std::vector< double> > v0_; // v0_[nsp_][3*na_] 44 std::vector<double> pmass_; // pmass_[nsp_] 45 46 public: 47 IonicStepper(Sample & s)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 // if there are more constraints than dofs, set ndofs_ to zero 53 if ( ndofs_ < 0 ) ndofs_ = 0; 54 nsp_ = atoms_.nsp(); 55 na_.resize(nsp_); 56 r0_.resize(nsp_); 57 rp_.resize(nsp_); 58 v0_.resize(nsp_); 59 pmass_.resize(nsp_); 60 for ( int is = 0; is < nsp_; is++ ) 61 { 62 const int nais = atoms_.na(is); 63 na_[is] = nais; 64 r0_[is].resize(3*nais); 65 rp_[is].resize(3*nais); 66 v0_[is].resize(3*nais); 67 pmass_[is] = atoms_.species_list[is]->mass() * 1822.89; 68 } 69 natoms_ = atoms_.size(); 70 // get positions and velocities from atoms_ 71 get_positions(); 72 get_velocities(); 73 } 74 r0(int is,int i)75 double r0(int is, int i) const { return r0_[is][i]; } v0(int is,int i)76 double v0(int is, int i) const { return v0_[is][i]; } r0(void)77 const std::vector<std::vector<double> >& r0(void) const { return r0_; } v0(void)78 const std::vector<std::vector<double> >& v0(void) const { return v0_; } pmass(void)79 const std::vector<double>& pmass(void) const { return pmass_; } get_positions(void)80 void get_positions(void) { atoms_.get_positions(r0_); } get_velocities(void)81 void get_velocities(void) { atoms_.get_velocities(v0_); } set_positions(void)82 void set_positions(void) { atoms_.set_positions(r0_); } set_velocities(void)83 void set_velocities(void) { atoms_.set_velocities(v0_); } 84 setup_constraints(void)85 void setup_constraints(void) 86 { 87 constraints_.setup(atoms_); 88 } 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; reset(void)93 virtual void reset(void) {} ekin(void)94 virtual double ekin(void) const { return 0.0; } ekin_stepper(void)95 virtual double ekin_stepper(void) const { return 0.0; } temp(void)96 virtual double temp(void) const { return 0.0; } 97 ~IonicStepper()98 virtual ~IonicStepper() {} 99 100 }; 101 #endif 102