Commit 75b770f5 by Francois Gygi

new IonicStepper logic including constraints.


git-svn-id: http://qboxcode.org/svn/qb/trunk@406 cba15fb0-1239-40c8-b417-11db7ca47a34
parent 1bcfb529
......@@ -3,12 +3,13 @@
// IonicStepper.h:
//
////////////////////////////////////////////////////////////////////////////////
// $Id: IonicStepper.h,v 1.6 2004-12-17 23:31:00 fgygi Exp $
// $Id: IonicStepper.h,v 1.7 2005-06-27 22:23:51 fgygi Exp $
#ifndef IONICSTEPPER_H
#define IONICSTEPPER_H
#include "Sample.h"
#include "Species.h"
#include <vector>
using namespace std;
......@@ -18,20 +19,25 @@ class IonicStepper
Sample& s_;
AtomSet& atoms_;
ConstraintSet& constraints_;
double dt_;
int nsp_;
int ndofs_;
vector<int> na_; // number of atoms per species na_[nsp_]
vector<vector< double> > r0_; // r0_[nsp_][3*na_]
vector<vector< double> > rp_; // rp_[nsp_][3*na_]
vector<vector< double> > rm_; // rm_[nsp_][3*na_]
vector<vector< double> > v0_; // v0_[nsp_][3*na_]
vector<double> pmass_; // pmass_[nsp_]
public:
IonicStepper (Sample& s) : s_(s), atoms_(s.atoms), dt_(s.ctrl.dt)
IonicStepper (Sample& s) : s_(s), atoms_(s.atoms),
constraints_(s.constraints), dt_(s.ctrl.dt)
{
ndofs_ = 3 * s.atoms.size();
ndofs_ = 3 * atoms_.size() - constraints_.size();
// if there are more constraints than dofs, set ndofs_ to zero
if ( ndofs_ < 0 ) ndofs_ = 0;
nsp_ = atoms_.nsp();
na_.resize(nsp_);
r0_.resize(nsp_);
......@@ -48,19 +54,21 @@ class IonicStepper
pmass_[is] = atoms_.species_list[is]->mass() * 1822.89;
}
atoms_.get_positions(r0_);
atoms_.get_velocities(v0_);
}
double r0(int is, int i) const { return r0_[is][i]; }
double rp(int is, int i) const { return rp_[is][i]; }
double v0(int is, int i) const { return v0_[is][i]; }
const vector<vector<double> >& r0(void) const { return r0_; }
const vector<vector<double> >& rp(void) const { return rp_; }
const vector<vector<double> >& v0(void) const { return v0_; }
const vector<double>& pmass(void) const { return pmass_; }
virtual void compute_rp(const vector<vector< double> >& f0) = 0;
virtual void compute_v0(const vector<vector< double> >& f0) = 0;
virtual void update_r(void) = 0;
virtual void update_v(void) = 0;
void setup_constraints(void)
{
constraints_.setup(atoms_);
}
virtual void compute_r(double e0, const vector<vector< double> >& f0) = 0;
virtual void compute_v(double e0, const vector<vector< double> >& f0) = 0;
virtual void reset(void) {}
virtual double ekin(void) const { return 0.0; }
virtual double temp(void) const { return 0.0; }
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment