Commit e229a784 by Francois Gygi

release 1.13.0 Implemented stress calculation


git-svn-id: http://qboxcode.org/svn/qb/trunk@181 cba15fb0-1239-40c8-b417-11db7ca47a34
parent 2d8665b2
......@@ -3,7 +3,7 @@
// AtomSet.h
//
////////////////////////////////////////////////////////////////////////////////
// $Id: AtomSet.h,v 1.9 2003-11-20 20:42:24 fgygi Exp $
// $Id: AtomSet.h,v 1.10 2004-03-11 21:52:32 fgygi Exp $
#ifndef ATOMSET_H
#define ATOMSET_H
......@@ -53,6 +53,7 @@ class AtomSet
void set_positions(const vector<vector<double> >& tau);
void get_velocities(vector<vector<double> >& vel) const;
void set_velocities(const vector<vector<double> >& vel);
void reset_velocities(void);
int size(void);
};
ostream& operator << ( ostream &os, AtomSet &as );
......
......@@ -3,7 +3,7 @@
// BOSampleStepper.h
//
////////////////////////////////////////////////////////////////////////////////
// $Id: BOSampleStepper.h,v 1.1 2003-11-21 20:01:06 fgygi Exp $
// $Id: BOSampleStepper.h,v 1.2 2004-03-11 21:52:32 fgygi Exp $
#ifndef BOSAMPLESTEPPER_H
#define BOSAMPLESTEPPER_H
......@@ -22,8 +22,8 @@ class BOSampleStepper : public SampleStepper
Wavefunction dwf;
Wavefunction* wfv;
vector<vector<double> > tau0,taum,vel,fion;
int nite_;
EnergyFunctional& ef_;
WavefunctionStepper* wf_stepper;
IonicStepper* ionic_stepper;
......@@ -35,9 +35,9 @@ class BOSampleStepper : public SampleStepper
mutable TimerMap tmap;
void step(EnergyFunctional& e, int niter);
void step(int niter);
BOSampleStepper(Sample& s, int nite);
BOSampleStepper(Sample& s, EnergyFunctional& ef, int nite);
//~BOSampleStepper();
};
#endif
////////////////////////////////////////////////////////////////////////////////
//
// SampleStepper.C
// CPSampleStepper.C
//
////////////////////////////////////////////////////////////////////////////////
// $Id: CPSampleStepper.C,v 1.2 2004-02-04 19:55:16 fgygi Exp $
// $Id: CPSampleStepper.C,v 1.3 2004-03-11 21:52:32 fgygi Exp $
#include "CPSampleStepper.h"
#include "EnergyFunctional.h"
#include "SlaterDet.h"
#include "MDWavefunctionStepper.h"
#include "MDIonicStepper.h"
#include "SDCellStepper.h"
#include "Basis.h"
#include <iostream>
......@@ -17,25 +18,28 @@
using namespace std;
////////////////////////////////////////////////////////////////////////////////
CPSampleStepper::CPSampleStepper(Sample& s) : SampleStepper(s), dwf(s.wf),
wfv(s.wfv)
CPSampleStepper::CPSampleStepper(Sample& s, EnergyFunctional& ef) :
SampleStepper(s), ef_(ef), dwf(s.wf), wfv(s.wfv)
{
mdwf_stepper = new MDWavefunctionStepper(s,tmap);
assert(mdwf_stepper!=0);
mdionic_stepper = 0;
if ( s.ctrl.atoms_dyn != "LOCKED" )
mdionic_stepper = new MDIonicStepper(s);
fion.resize(s_.atoms.nsp());
for ( int is = 0; is < fion.size(); is++ )
fion[is].resize(3*s_.atoms.na(is));
}
////////////////////////////////////////////////////////////////////////////////
void CPSampleStepper::step(EnergyFunctional& e, int niter)
CPSampleStepper::~CPSampleStepper(void)
{
delete mdwf_stepper;
if ( mdionic_stepper != 0 ) delete mdionic_stepper;
}
////////////////////////////////////////////////////////////////////////////////
void CPSampleStepper::step(int niter)
{
AtomSet& atoms = s_.atoms;
Wavefunction& wf = s_.wf;
valarray<double> sigma(6);
const double dt = s_.ctrl.dt;
double ekin_ion=0.0,ekin_e, temp_ion, eta;
......@@ -43,16 +47,29 @@ void CPSampleStepper::step(EnergyFunctional& e, int niter)
const string wf_dyn = s_.ctrl.wf_dyn;
assert(wf_dyn=="MD");
const string atoms_dyn = s_.ctrl.atoms_dyn;
const string cell_dyn = s_.ctrl.cell_dyn;
const bool compute_hpsi = true;
const bool compute_forces = ( atoms_dyn != "LOCKED" );
const bool compute_stress = false;
const bool compute_stress = ( s_.ctrl.stress == "ON" );
CellStepper* cell_stepper = 0;
if ( cell_dyn == "SD" )
cell_stepper = new SDCellStepper(s_);
if ( s_.wfv == 0 )
{
s_.wfv = new Wavefunction(wf);
s_.wfv->clear();
}
Timer tm_iter;
double energy =
e.energy(compute_hpsi,dwf,compute_forces,fion,compute_stress,sigma);
ef_.energy(compute_hpsi,dwf,compute_forces,fion,compute_stress,sigma_eks);
mdwf_stepper->compute_wfm(dwf);
for ( int iter = 0; iter < niter; iter++ )
{
tm_iter.reset();
......@@ -60,10 +77,7 @@ void CPSampleStepper::step(EnergyFunctional& e, int niter)
if ( s_.ctxt_.mype() == 0 )
cout << " <iteration count=\"" << iter+1 << "\">\n";
if ( iter == 0 )
mdwf_stepper->stoermer_start(dwf);
else
mdwf_stepper->update(dwf);
mdwf_stepper->update(dwf);
ekin_e = mdwf_stepper->ekin();
......@@ -72,20 +86,30 @@ void CPSampleStepper::step(EnergyFunctional& e, int niter)
cout.setf(ios::fixed,ios::floatfield);
cout.setf(ios::right,ios::adjustfield);
cout << " <ekin> " << setprecision(8)
<< setw(15) << e.ekin() << " </ekin>\n"
<< " <eps> " << setw(15) << e.eps() << " </eps>\n"
<< " <enl> " << setw(15) << e.enl() << " </enl>\n"
<< " <ecoul> " << setw(15) << e.ecoul() << " </ecoul>\n"
<< " <exc> " << setw(15) << e.exc() << " </exc>\n"
<< " <esr> " << setw(15) << e.esr() << " </esr>\n"
<< " <eself> " << setw(15) << e.eself() << " </eself>\n"
<< " <etotal> " << setw(15) << e.etotal() << " </etotal>\n"
<< setw(15) << ef_.ekin() << " </ekin>\n";
if ( compute_stress )
{
cout << " <econf_int> " << setw(15) << ef_.econf()
<< " </econf_int>\n";
}
cout << " <eps> " << setw(15) << ef_.eps() << " </eps>\n"
<< " <enl> " << setw(15) << ef_.enl() << " </enl>\n"
<< " <ecoul> " << setw(15) << ef_.ecoul() << " </ecoul>\n"
<< " <exc> " << setw(15) << ef_.exc() << " </exc>\n"
<< " <esr> " << setw(15) << ef_.esr() << " </esr>\n"
<< " <eself> " << setw(15) << ef_.eself() << " </eself>\n"
<< " <etotal> " << setw(15) << ef_.etotal() << " </etotal>\n"
<< flush;
}
if ( compute_forces )
{
if ( s_.wf.context().onpe0() )
mdionic_stepper->compute_v0(fion);
mdionic_stepper->update_v();
mdionic_stepper->compute_rp(fion);
ekin_ion = mdionic_stepper->ekin();
if ( s_.ctxt_.onpe0() )
{
for ( int is = 0; is < atoms.atom_list.size(); is++ )
{
......@@ -97,28 +121,26 @@ void CPSampleStepper::step(EnergyFunctional& e, int niter)
<< " species=\"" << pa->species()
<< "\">\n"
<< " <position> "
<< mdionic_stepper->tau0(is,i) << " "
<< mdionic_stepper->tau0(is,i+1) << " "
<< mdionic_stepper->tau0(is,i+2) << " </position>\n"
<< " <force> " << fion[is][i] << " "
<< fion[is][i+1] << " " << fion[is][i+2]
<< mdionic_stepper->r0(is,i) << " "
<< mdionic_stepper->r0(is,i+1) << " "
<< mdionic_stepper->r0(is,i+2) << " </position>\n"
<< " <velocity> "
<< mdionic_stepper->v0(is,i) << " "
<< mdionic_stepper->v0(is,i+1) << " "
<< mdionic_stepper->v0(is,i+2) << " </velocity>\n"
<< " <force> "
<< fion[is][i] << " "
<< fion[is][i+1] << " "
<< fion[is][i+2]
<< " </force>\n </atom>" << endl;
i += 3;
}
}
}
if ( iter == 0 )
mdionic_stepper->preprocess(fion);
mdionic_stepper->update(fion);
ekin_ion = mdionic_stepper->ekin();
e.atoms_moved();
}
if ( s_.wf.context().onpe0() )
if ( s_.ctxt_.onpe0() )
{
cout << " <ekin_e> " << ekin_e << " </ekin_e>\n";
......@@ -132,8 +154,31 @@ void CPSampleStepper::step(EnergyFunctional& e, int niter)
cout << " <ekin_ec> " << energy+ekin_ion+2*ekin_e << " </ekin_ec>\n";
}
if ( compute_stress )
{
compute_sigma();
print_stress();
if ( cell_dyn != "LOCKED" )
{
cell_stepper->compute_new_cell(sigma);
// Update cell
cell_stepper->update_cell();
ef_.cell_moved();
ef_.atoms_moved();
}
}
if ( compute_forces )
{
mdionic_stepper->update_r();
ef_.atoms_moved();
}
energy =
e.energy(compute_hpsi,dwf,compute_forces,fion,compute_stress,sigma);
ef_.energy(compute_hpsi,dwf,compute_forces,fion,compute_stress,sigma_eks);
if ( s_.ctxt_.mype() == 0 )
cout << " </iteration>" << endl;
......@@ -157,10 +202,14 @@ void CPSampleStepper::step(EnergyFunctional& e, int niter)
// dwf and fion now contain the forces on wavefunctions and ions at the
// endpoint
mdwf_stepper->stoermer_end(dwf); // replace wfm by wfv
mdwf_stepper->compute_wfv(dwf); // replace wfm by wfv
if ( compute_forces )
{
mdionic_stepper->postprocess(fion);
// Note: next line function call updates velocities in the AtomSet
mdionic_stepper->compute_v0(fion);
mdionic_stepper->update_v();
}
if ( cell_stepper != 0 ) delete cell_stepper;
}
......@@ -3,7 +3,7 @@
// CPSampleStepper.h
//
////////////////////////////////////////////////////////////////////////////////
// $Id: CPSampleStepper.h,v 1.1 2003-11-21 20:01:06 fgygi Exp $
// $Id: CPSampleStepper.h,v 1.2 2004-03-11 21:52:31 fgygi Exp $
#ifndef CPSAMPLESTEPPER_H
#define CPSAMPLESTEPPER_H
......@@ -20,9 +20,9 @@ class CPSampleStepper : public SampleStepper
{
private:
EnergyFunctional& ef_;
Wavefunction dwf;
Wavefunction* wfv;
vector<vector<double> > fion;
MDWavefunctionStepper* mdwf_stepper;
MDIonicStepper* mdionic_stepper;
......@@ -34,9 +34,9 @@ class CPSampleStepper : public SampleStepper
mutable TimerMap tmap;
void step(EnergyFunctional& e, int niter);
void step(int niter);
CPSampleStepper(Sample& s);
//~CPSampleStepper();
CPSampleStepper(Sample& s, EnergyFunctional& ef);
~CPSampleStepper();
};
#endif
......@@ -3,7 +3,7 @@
// CellDyn.h
//
////////////////////////////////////////////////////////////////////////////////
// $Id: CellDyn.h,v 1.1 2004-02-04 19:52:27 fgygi Exp $
// $Id: CellDyn.h,v 1.2 2004-03-11 21:52:31 fgygi Exp $
#ifndef CELLDYN_H
#define CELLDYN_H
......@@ -35,10 +35,10 @@ class CellDyn : public Var
}
string v = argv[1];
if ( !( v == "LOCKED" || v == "SD" || v == "MD" ) )
if ( !( v == "LOCKED" || v == "SD" ) )
{
if ( ui->onpe0() )
cout << " cell_dyn must be in [LOCKED,SD,MD]" << endl;
cout << " cell_dyn must be in [LOCKED,SD]" << endl;
return 1;
}
......
......@@ -3,7 +3,7 @@
// Control.h:
//
////////////////////////////////////////////////////////////////////////////////
// $Id: Control.h,v 1.6 2004-02-04 19:55:17 fgygi Exp $
// $Id: Control.h,v 1.7 2004-03-11 21:52:32 fgygi Exp $
#ifndef CONTROL_H
#define CONTROL_H
......@@ -14,6 +14,7 @@
struct Control
{
// control variables
string debug, timing;
string wf_dyn, atoms_dyn; // dynamics string flags
int nite;
double emass; // electron mass
......@@ -31,7 +32,7 @@ struct Control
double gms_mix; // mixing factor for generalized minimum spread functions
string thermostat;
double th_temp,th_time; // thermostat control
double th_temp,th_time, th_width; // thermostat control
string stress;
string cell_dyn;
......
......@@ -3,7 +3,7 @@
// EnergyFunctional.h
//
////////////////////////////////////////////////////////////////////////////////
// $Id: EnergyFunctional.h,v 1.10 2004-02-04 19:55:17 fgygi Exp $
// $Id: EnergyFunctional.h,v 1.11 2004-03-11 21:52:32 fgygi Exp $
#ifndef ENERGYFUNCTIONAL_H
#define ENERGYFUNCTIONAL_H
......@@ -26,6 +26,7 @@ class UnitCell;
class FourierTransform;
class XCPotential;
class NonLocalPotential;
class ConfinementPotential;
typedef map<string,Timer> TimerMap;
......@@ -41,6 +42,7 @@ class EnergyFunctional
StructureFactor sf;
XCPotential* xcp;
NonLocalPotential* nlp;
vector<ConfinementPotential*> cfp; // cfp[ikp]
vector<vector<double> > vps, dvps, rhops;
vector<complex<double> > tmp_r, vion_local_g, dvion_local_g, vlocal_g,
......@@ -55,12 +57,9 @@ class EnergyFunctional
int nsp_;
double ekin_, econf_, eps_, enl_, ehart_,
ecoul_, exc_, esr_, eself_, etotal_;
valarray<double> fstress_, dfstress_;
valarray<double> sigma_ekin,sigma_econf,sigma_eps,sigma_ehart,sigma_exc,
sigma_enl, sigma_esr, sigma;
void init(void);
public:
mutable TimerMap tmap;
......@@ -80,10 +79,15 @@ class EnergyFunctional
double esr(void) const { return esr_; }
double eself(void) const { return eself_; }
const ConfinementPotential *confpot(int ikp) const { return cfp[ikp]; }
void atoms_moved(void);
void cell_moved(void);
void print(ostream& os) const;
EnergyFunctional(const Sample& s);
~EnergyFunctional();
};
ostream& operator << ( ostream& os, const EnergyFunctional& e );
#endif
......@@ -3,7 +3,7 @@
// IonicStepper.h:
//
////////////////////////////////////////////////////////////////////////////////
// $Id: IonicStepper.h,v 1.2 2003-11-27 01:18:41 fgygi Exp $
// $Id: IonicStepper.h,v 1.3 2004-03-11 21:52:31 fgygi Exp $
#ifndef IONICSTEPPER_H
#define IONICSTEPPER_H
......@@ -21,8 +21,9 @@ class IonicStepper
int nsp_;
int ndofs_;
vector<int> na_; // number of atoms per species na_[nsp_]
vector<vector< double> > tau0_; // tau0_[nsp_][3*na_]
vector<vector< double> > vel_; // vel_[nsp_][3*na_]
vector<vector< double> > r0_; // r0_[nsp_][3*na_]
vector<vector< double> > rp_; // rp_[nsp_][3*na_]
vector<vector< double> > v0_; // v0_[nsp_][3*na_]
vector<double> pmass_; // pmass_[nsp_]
double ekin_; // kinetic energy
......@@ -33,27 +34,33 @@ class IonicStepper
ndofs_ = 3 * s.atoms.size();
nsp_ = atoms_.nsp();
na_.resize(nsp_);
tau0_.resize(nsp_);
vel_.resize(nsp_);
r0_.resize(nsp_);
rp_.resize(nsp_);
v0_.resize(nsp_);
pmass_.resize(nsp_);
for ( int is = 0; is < nsp_; is++ )
{
const int nais = atoms_.na(is);
na_[is] = nais;
tau0_[is].resize(3*nais);
vel_[is].resize(3*nais);
r0_[is].resize(3*nais);
rp_[is].resize(3*nais);
v0_[is].resize(3*nais);
pmass_[is] = atoms_.species_list[is]->mass() * 1822.89;
}
atoms_.get_positions(tau0_);
atoms_.get_velocities(vel_);
atoms_.get_positions(r0_);
}
double tau0(int is, int i) const { return tau0_[is][i]; }
double vel(int is, int i) const { return vel_[is][i]; }
const vector<vector<double> >& tau0(void) const { return tau0_; }
virtual void preprocess(const vector<vector<double> >&fion) {}
virtual void postprocess(const vector<vector<double> >&fion) {}
virtual void update(const vector<vector< double> >& fion) = 0;
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_; }
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;
double ekin(void) const { return ekin_; }
double temp(void) const
{
......
......@@ -3,111 +3,84 @@
// MDIonicStepper.C
//
////////////////////////////////////////////////////////////////////////////////
// $Id: MDIonicStepper.C,v 1.3 2003-12-02 20:25:28 fgygi Exp $
// $Id: MDIonicStepper.C,v 1.4 2004-03-11 21:52:32 fgygi Exp $
#include "MDIonicStepper.h"
////////////////////////////////////////////////////////////////////////////////
void MDIonicStepper::update(const vector<vector< double> >& fion)
void MDIonicStepper::compute_rp(const vector<vector< double> >& f0)
{
eta_ = 0.0;
// Verlet update
if ( thermostat_ )
// f0 contains forces at r0
// Compute new positions rp using the Verlet algorithm
// compute velocity at t+1/2 vhalf_
ekin_ = 0.0;
atoms_.get_velocities(v0_);
for ( int is = 0; is < v0_.size(); is++ )
{
double ekin_ion = 0.0;
// compute damping factor eta
// compute ekin_ion and temp_ion before step using a first order
// approximation. Note: ekin_ is recomputed after the step using
// a second-order approximation.
for ( int is = 0; is < tau0_.size(); is++ )
const double dt2bym = dt_ * dt_ / pmass_[is];
const double dtby2m = dt_ / ( 2.0 * pmass_[is] );
if ( dt_ != 0.0 )
{
const double dt2bym = dt_ * dt_ / pmass_[is];
if ( dt_ != 0.0 )
for ( int i = 0; i < v0_[is].size(); i++ )
{
for ( int i = 0; i < tau0_[is].size(); i++ )
{
const double v = ( tau0_[is][i] - taum_[is][i] ) / dt_;
ekin_ion += 0.5 * pmass_[is] * v * v;
}
const double v = v0_[is][i] + dtby2m * f0[is][i];
ekin_ += 0.5 * pmass_[is] * v0_[is][i] * v0_[is][i];
vhalf_[is][i] = v;
}
}
// Next line: linear thermostat, width of tanh = 100.0
eta_ = tanh ( ( temp() - th_temp_ ) / 100. ) / th_time_;
}
// ekin_ is the kinetic energy computed from v0
if ( thermostat_ )
{
eta_ = tanh ( ( temp() - th_temp_ ) / th_width_ ) / th_time_;
}
// make Verlet step with damping eta, and recompute ekin_ to second order
ekin_ = 0.0;
for ( int is = 0; is < tau0_.size(); is++ )
// compute rp
atoms_.get_positions(r0_);
for ( int is = 0; is < r0_.size(); is++ )
{
const double dt2bym = dt_ * dt_ / pmass_[is];
for ( int i = 0; i < tau0_[is].size(); i++ )
const double dt2by2m = dt_ * dt_ / ( 2.0 * pmass_[is] );
for ( int i = 0; i < r0_[is].size(); i++ )
{
const double taup = tau0_[is][i] +
(1.0 - dt_*eta_) * (tau0_[is][i] - taum_[is][i]) +
dt2bym * fion[is][i];
if ( dt_ != 0.0 )
{
const double v = 0.5 * ( taup - taum_[is][i] ) / dt_;
ekin_ += 0.5 * pmass_[is] * v * v;
vel_[is][i] = v;
}
taum_[is][i] = tau0_[is][i];
tau0_[is][i] = taup;
rp_[is][i] = r0_[is][i] +
(1.0 - dt_*eta_) * v0_[is][i] * dt_ + dt2by2m * f0[is][i];
}
}
atoms_.set_positions(tau0_);
atoms_.set_velocities(vel_);
}
////////////////////////////////////////////////////////////////////////////////
void MDIonicStepper::stoermer_start(const vector<vector< double> >& fion)
void MDIonicStepper::compute_v0(const vector<vector< double> >& f0)
{
// compute taum from tau0,vel,fion and dt before first step
// First step of Stoermer's rule
// x1 = x0 + h * v0 + 0.5 * dt2/m * f0
// compute velocities v0_ using vhalf_ and force f0(r0)
// Note: vhalf contains the velocity at t-1/2 since compute_v0 is called
// after a call to update_r()
// x-1 = x0 - h * v0 + 0.5 * dt2/m * f0
atoms_.get_velocities(vel_);
ekin_ = 0.0;
for ( int is = 0; is < tau0_.size(); is++ )
for ( int is = 0; is < vhalf_.size(); is++ )
{
const double dt2bym = dt_ * dt_ / pmass_[is];
for ( int i = 0; i < tau0_[is].size(); i++ )
assert(pmass_[is] > 0.0);
const double dtby2m = dt_ / ( 2.0 * pmass_[is] );
for ( int i = 0; i < vhalf_[is].size(); i++ )
{
const double taum = tau0_[is][i] - dt_ * vel_[is][i] +
0.5 * dt2bym * fion[is][i];
taum_[is][i] = taum;
const double v = vel_[is][i];
const double v = vhalf_[is][i] + dtby2m * f0[is][i];
ekin_ += 0.5 * pmass_[is] * v * v;
v0_[is][i] = v;
}
}
// ekin_ contains the kinetic energy computed from v0
}
////////////////////////////////////////////////////////////////////////////////
void MDIonicStepper::stoermer_end(const vector<vector< double> >& fion)
void MDIonicStepper::update_r(void)
{
// compute vel from tau0, taum, fion and dt at last step
// Last step of Stoermer's rule
// v = (xn - x(n-1) )/dt + 0.5 * dt/m * fn
ekin_ = 0.0;
if ( dt_ != 0.0 )
{
for ( int is = 0; is < tau0_.size(); is++ )
{
const double dtbym = dt_ / pmass_[is];
for ( int i = 0; i < tau0_[is].size(); i++ )
{
vel_[is][i] = ( tau0_[is][i] - taum_[is][i] ) / dt_ +
0.5 * dtbym * fion[is][i];
const double v = vel_[is][i];
ekin_ += 0.5 * pmass_[is] * v * v;
}
}
}
atoms_.set_velocities(vel_);
r0_ = rp_;
atoms_.set_positions(r0_);
}
////////////////////////////////////////////////////////////////////////////////
void MDIonicStepper::update_v(void)
{
atoms_.set_velocities(v0_);
}
......@@ -3,7 +3,31 @@
// MDIonicStepper.h:
//
////////////////////////////////////////////////////////////////////////////////
// $Id: MDIonicStepper.h,v 1.3 2003-12-02 20:25:28 fgygi Exp $
// $Id: MDIonicStepper.h,v 1.4 2004-03-11 21:52:31 fgygi Exp $
//
// IonicStepper is used in the following way
//
// input: r0,v0
//
// for ( k=0; k<n; k++ )
// {
// compute forces f0(r0)
// if ( k > 0 ) stepper->compute_v0(f0);
// optional: modify velocities
// stepper->update_v();
// // new r0,v0,f0 known
// stepper->compute_rp(f0); // using r0, v0 and f0
// optional: restore constraints using rp and r0
// stepper->update_r(); // r0 <- rp
// // Note: r0 and v0 are not consistent at this point
// }
// compute forces f0(r0)
// stepper->compute_v0(f0);
// stepper->update_v();
//
// // r0,v0,f0 consistent at this point
//
#ifndef MDIONICSTEPPER_H
#define MDIONICSTEPPER_H
......@@ -14,13 +38,12 @@ class MDIonicStepper : public IonicStepper
{
private:
vector<vector<double> > taum_;
double th_temp_;
double th_time_;
double th_width_;
double eta_;
bool thermostat_;
void stoermer_start(const vector<vector<double> >& fion);
void stoermer_end(const vector<vector<double> >& fion);
vector<vector< double> > vhalf_; // vhalf_[nsp_][3*na_]: v(t+dt/2)
public:
......@@ -29,15 +52,21 @@ class MDIonicStepper : public IonicStepper
thermostat_ = ( s.ctrl.thermostat == "ON" );
th_temp_ = s.ctrl.th_temp;
th_time_ = s.ctrl.th_time;
th_width_ = s.ctrl.th_width;
eta_ = 0.0;
taum_.resize(tau0_.size());
for ( int is = 0; is < taum_.size(); is++ )
taum_[is].resize(tau0_[is].size());
vhalf_.resize(nsp_);
for ( int is = 0; is < nsp_; is++ )
{
const int nais = atoms_.na(is);
vhalf_[is].resize(3*nais);
}
atoms_.get_velocities(v0_);
}
void update(const vector<vector< double> >& fion);
void preprocess(const vector<vector<double> >& fion) { stoermer_start(fion);}
void postprocess(const vector<vector<double> >& fion) { stoermer_end(fion);}
void compute_rp(const vector<vector< double> >& f0);
void compute_v0(const vector<vector< double> >& f0);