Commit d583b5d3 by Francois Gygi

New classes for the SampleStepper hierarchy


git-svn-id: http://qboxcode.org/svn/qb/trunk@123 cba15fb0-1239-40c8-b417-11db7ca47a34
parent 02d4895e
////////////////////////////////////////////////////////////////////////////////
//
// BOSampleStepper.C
//
////////////////////////////////////////////////////////////////////////////////
// $Id: BOSampleStepper.C,v 1.1 2003-11-21 20:01:06 fgygi Exp $
#include "BOSampleStepper.h"
#include "EnergyFunctional.h"
#include "SlaterDet.h"
#include "Basis.h"
#include "WavefunctionStepper.h"
#include "SDWavefunctionStepper.h"
#include "PSDWavefunctionStepper.h"
#include "SDIonicStepper.h"
#include "MDIonicStepper.h"
#include <iostream>
#include <iomanip>
using namespace std;
////////////////////////////////////////////////////////////////////////////////
BOSampleStepper::BOSampleStepper(Sample& s, int nite) : SampleStepper(s),
dwf(s.wf), wfv(s.wfv), nite_(nite)
{
fion.resize(s_.atoms.nsp());
for ( int is = 0; is < fion.size(); is++ )
fion[is].resize(3*s_.atoms.na(is));
}
////////////////////////////////////////////////////////////////////////////////
void BOSampleStepper::step(EnergyFunctional& e, int niter)
{
AtomSet& atoms = s_.atoms;
Wavefunction& wf = s_.wf;
UnitCell dcell;
atoms.get_positions(tau0);
atoms.get_velocities(vel);
const double dt = s_.ctrl.dt;
double ekin_ion=0.0, temp_ion=0.0, eta;
const string wf_dyn = s_.ctrl.wf_dyn;
const string atoms_dyn = s_.ctrl.atoms_dyn;
const bool compute_hpsi = ( wf_dyn != "LOCKED" );
const bool compute_forces = ( atoms_dyn != "LOCKED" );
const bool compute_stress = ( s_.ctrl.stress == "ON" );
const bool ttherm = ( s_.ctrl.thermostat == "ON" );
const int ndofs = 3 * s_.atoms.size();
const double boltz = 1.0 / ( 11605.0 * 2.0 * 13.6058 );
const double th_temp = s_.ctrl.th_temp;
const double th_time = s_.ctrl.th_time;
Timer tm_iter;
WavefunctionStepper* wf_stepper = 0;
if ( wf_dyn == "SD" )
wf_stepper = new SDWavefunctionStepper(s_,tmap);
else if ( wf_dyn == "PSD" )
wf_stepper = new PSDWavefunctionStepper(s_,tmap);
assert(wf_stepper!=0);
IonicStepper* ionic_stepper = 0;
if ( atoms_dyn == "SD" )
ionic_stepper = new SDIonicStepper(s_);
else if ( atoms_dyn == "MD" )
ionic_stepper = new MDIonicStepper(s_);
for ( int iter = 0; iter < niter; iter++ )
{
// ionic iteration
tm_iter.start();
if ( s_.ctxt_.onpe0() )
cout << " <iteration count=\"" << iter+1 << "\">\n";
// wavefunction extrapolation
if ( s_.wfv != 0 && atoms_dyn != "LOCKED" )
{
//wf_stepper->extrapolate(wf,wfv);
}
double energy;
// do nite-1 electronic steps without forces
if ( wf_stepper != 0 )
{
for ( int ite = 0; ite < nite_-1; ite++ )
{
energy = e.energy(true,dwf,false,fion,false,dcell);
wf_stepper->update(dwf);
if ( s_.ctxt_.onpe0() )
{
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"
<< flush;
}
}
}
// last electronic iteration
energy =
e.energy(compute_hpsi,dwf,compute_forces,fion,compute_stress,dcell);
if ( s_.ctxt_.onpe0() )
{
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"
<< flush;
}
// make ionic step
if ( compute_forces )
{
if ( s_.wf.context().onpe0() )
{
for ( int is = 0; is < atoms.atom_list.size(); is++ )
{
int i = 0;
for ( int ia = 0; ia < atoms.atom_list[is].size(); ia++ )
{
Atom* pa = atoms.atom_list[is][ia];
cout << " <atom name=\"" << pa->name() << "\""
<< " species=\"" << pa->species()
<< "\">\n"
<< " <position> "
<< ionic_stepper->tau0(is,i) << " "
<< ionic_stepper->tau0(is,i+1) << " "
<< ionic_stepper->tau0(is,i+2) << " </position>\n"
<< " <force> " << fion[is][i] << " "
<< fion[is][i+1] << " " << fion[is][i+2]
<< " </force>\n </atom>" << endl;
i += 3;
}
}
}
if ( iter == 0 )
ionic_stepper->preprocess(fion);
ionic_stepper->update(fion);
ekin_ion = ionic_stepper->ekin();
if ( iter == niter-1 )
ionic_stepper->postprocess(fion);
e.atoms_moved();
}
if ( s_.ctxt_.onpe0() )
{
cout << " <econst> " << energy+ekin_ion << " </econst>\n";
}
// print iteration time
double time = tm_iter.real();
double tmin = time;
double tmax = time;
s_.ctxt_.dmin(1,1,&tmin,1);
s_.ctxt_.dmax(1,1,&tmax,1);
if ( s_.ctxt_.onpe0() )
{
cout << " <!-- timing "
<< setw(15) << "iteration"
<< " : " << setprecision(3) << setw(9) << tmin
<< " " << setprecision(3) << setw(9) << tmax << " -->" << endl;
cout << " </iteration>" << endl;
}
} // for iter
}
////////////////////////////////////////////////////////////////////////////////
//
// BOSampleStepper.h
//
////////////////////////////////////////////////////////////////////////////////
// $Id: BOSampleStepper.h,v 1.1 2003-11-21 20:01:06 fgygi Exp $
#ifndef BOSAMPLESTEPPER_H
#define BOSAMPLESTEPPER_H
#include "SampleStepper.h"
#include "Sample.h"
#include "Wavefunction.h"
class EnergyFunctional;
class WavefunctionStepper;
class IonicStepper;
using namespace std;
class BOSampleStepper : public SampleStepper
{
private:
Wavefunction dwf;
Wavefunction* wfv;
vector<vector<double> > tau0,taum,vel,fion;
int nite_;
WavefunctionStepper* wf_stepper;
IonicStepper* ionic_stepper;
// Do not allow construction of BOSampleStepper unrelated to a Sample
BOSampleStepper(void);
public:
mutable TimerMap tmap;
void step(EnergyFunctional& e, int niter);
BOSampleStepper(Sample& s, int nite);
//~BOSampleStepper();
};
#endif
////////////////////////////////////////////////////////////////////////////////
//
// SampleStepper.C
//
////////////////////////////////////////////////////////////////////////////////
// $Id: CPSampleStepper.C,v 1.1 2003-11-21 20:01:06 fgygi Exp $
#include "CPSampleStepper.h"
#include "EnergyFunctional.h"
#include "SlaterDet.h"
#include "MDWavefunctionStepper.h"
#include "MDIonicStepper.h"
#include "Basis.h"
#include <iostream>
#include <iomanip>
using namespace std;
////////////////////////////////////////////////////////////////////////////////
CPSampleStepper::CPSampleStepper(Sample& s) : SampleStepper(s), 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)
{
AtomSet& atoms = s_.atoms;
Wavefunction& wf = s_.wf;
UnitCell dcell;
const double dt = s_.ctrl.dt;
double ekin_ion=0.0,ekin_e, temp_ion, eta;
const string wf_dyn = s_.ctrl.wf_dyn;
assert(wf_dyn=="MD");
const string atoms_dyn = s_.ctrl.atoms_dyn;
const bool compute_hpsi = true;
const bool compute_forces = ( atoms_dyn != "LOCKED" );
const bool compute_stress = false;
Timer tm_iter;
double energy =
e.energy(compute_hpsi,dwf,compute_forces,fion,compute_stress,dcell);
for ( int iter = 0; iter < niter; iter++ )
{
tm_iter.reset();
tm_iter.start();
if ( s_.ctxt_.mype() == 0 )
cout << " <iteration count=\"" << iter+1 << "\">\n";
if ( iter == 0 )
mdwf_stepper->stoermer_start(dwf);
else
mdwf_stepper->update(dwf);
ekin_e = mdwf_stepper->ekin();
if ( s_.ctxt_.onpe0() )
{
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"
<< flush;
}
if ( compute_forces )
{
if ( s_.wf.context().onpe0() )
{
for ( int is = 0; is < atoms.atom_list.size(); is++ )
{
int i = 0;
for ( int ia = 0; ia < atoms.atom_list[is].size(); ia++ )
{
Atom* pa = atoms.atom_list[is][ia];
cout << " <atom name=\"" << pa->name() << "\""
<< " 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]
<< " </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() )
{
cout << " <ekin_e> " << ekin_e << " </ekin_e>\n";
if ( compute_forces )
{
cout << " <ekin_ion> " << ekin_ion << " </ekin_ion>\n";
cout << " <temp_ion> " << mdionic_stepper->temp() << " </temp_ion>\n";
cout << " <eta_ion> " << mdionic_stepper->eta() << " </eta_ion>\n";
}
cout << " <econst> " << energy+ekin_ion+ekin_e << " </econst>\n";
cout << " <ekin_ec> " << energy+ekin_ion+2*ekin_e << " </ekin_ec>\n";
}
energy =
e.energy(compute_hpsi,dwf,compute_forces,fion,compute_stress,dcell);
if ( s_.ctxt_.mype() == 0 )
cout << " </iteration>" << endl;
// print iteration time
tm_iter.stop();
double time = tm_iter.real();
double tmin = time;
double tmax = time;
s_.ctxt_.dmin(1,1,&tmin,1);
s_.ctxt_.dmax(1,1,&tmax,1);
if ( s_.ctxt_.myproc()==0 )
{
cout << " <!-- timing "
<< setw(15) << "iteration"
<< " : " << setprecision(3) << setw(9) << tmin
<< " " << setprecision(3) << setw(9) << tmax << " -->" << endl;
}
} // iter
// dwf and fion now contain the forces on wavefunctions and ions at the
// endpoint
mdwf_stepper->stoermer_end(dwf); // replace wfm by wfv
if ( compute_forces )
{
mdionic_stepper->postprocess(fion);
}
}
////////////////////////////////////////////////////////////////////////////////
//
// CPSampleStepper.h
//
////////////////////////////////////////////////////////////////////////////////
// $Id: CPSampleStepper.h,v 1.1 2003-11-21 20:01:06 fgygi Exp $
#ifndef CPSAMPLESTEPPER_H
#define CPSAMPLESTEPPER_H
#include "SampleStepper.h"
#include "Sample.h"
#include "Wavefunction.h"
class EnergyFunctional;
class MDWavefunctionStepper;
class MDIonicStepper;
using namespace std;
class CPSampleStepper : public SampleStepper
{
private:
Wavefunction dwf;
Wavefunction* wfv;
vector<vector<double> > fion;
MDWavefunctionStepper* mdwf_stepper;
MDIonicStepper* mdionic_stepper;
// Do not allow construction of CPSampleStepper unrelated to a Sample
CPSampleStepper(void);
public:
mutable TimerMap tmap;
void step(EnergyFunctional& e, int niter);
CPSampleStepper(Sample& s);
//~CPSampleStepper();
};
#endif
////////////////////////////////////////////////////////////////////////////////
//
// IonicStepper.h:
//
////////////////////////////////////////////////////////////////////////////////
// $Id: IonicStepper.h,v 1.1 2003-11-21 20:01:06 fgygi Exp $
#ifndef IONICSTEPPER_H
#define IONICSTEPPER_H
#include "Sample.h"
#include <vector>
using namespace std;
class IonicStepper
{
protected:
AtomSet& atoms_;
double dt_;
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<double> pmass_; // pmass_[nsp_]
public:
IonicStepper (Sample& s) : atoms_(s.atoms), dt_(s.ctrl.dt)
{
ndofs_ = 3 * s.atoms.size();
nsp_ = atoms_.nsp();
na_.resize(nsp_);
tau0_.resize(nsp_);
vel_.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);
pmass_[is] = atoms_.species_list[is]->mass() * 1822.89;
}
atoms_.get_positions(tau0_);
atoms_.get_velocities(vel_);
}
double tau0(int is, int i) const { return tau0_[is][i]; }
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;
virtual double ekin(void) const = 0;
virtual double temp(void) const = 0;
virtual ~IonicStepper() {}
};
#endif
////////////////////////////////////////////////////////////////////////////////
//
// MDIonicStepper.C
//
////////////////////////////////////////////////////////////////////////////////
// $Id: MDIonicStepper.C,v 1.1 2003-11-21 20:01:06 fgygi Exp $
#include "MDIonicStepper.h"
////////////////////////////////////////////////////////////////////////////////
double MDIonicStepper::temp(void) const
{
const double boltz = 1.0 / ( 11605.0 * 2.0 * 13.6058 );
if ( ndofs_ > 0.0 )
return 2.0 * ( ekin_ / boltz ) / ndofs_;
else
return 0.0;
}
////////////////////////////////////////////////////////////////////////////////
void MDIonicStepper::update(const vector<vector< double> >& fion)
{
// Verlet update
if ( thermostat_ )
{
eta_ = 0.0;
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];
if ( dt_ != 0.0 )
{
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;
}
}
}
// Next line: linear thermostat, width of tanh = 100.0
eta_ = tanh ( ( temp() - th_temp_ ) / 100. ) / 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++ )
{
const double dt2bym = dt_ * dt_ / pmass_[is];
for ( int i = 0; i < tau0_[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;
}
}
atoms_.set_positions(tau0_);
atoms_.set_velocities(vel_);
}
////////////////////////////////////////////////////////////////////////////////
void MDIonicStepper::stoermer_start(const vector<vector< double> >& fion)
{
// 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
// x-1 = x0 - h * v0 + 0.5 * dt2/m * f0
atoms_.get_velocities(vel_);
for ( int is = 0; is < tau0_.size(); is++ )
{
const double dt2bym = dt_ * dt_ / pmass_[is];
for ( int i = 0; i < tau0_[is].size(); i++ )
{
const double taum = tau0_[is][i] - dt_ * vel_[is][i] +
0.5 * dt2bym * fion[is][i];
taum_[is][i] = taum;
}
}
}
////////////////////////////////////////////////////////////////////////////////
void MDIonicStepper::stoermer_end(const vector<vector< double> >& fion)
{
// 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
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];
}
}
}
atoms_.set_velocities(vel_);
}
////////////////////////////////////////////////////////////////////////////////
//
// MDIonicStepper.h:
//
////////////////////////////////////////////////////////////////////////////////
// $Id: MDIonicStepper.h,v 1.1 2003-11-21 20:01:06 fgygi Exp $
#ifndef MDIONICSTEPPER_H
#define MDIONICSTEPPER_H
#include "IonicStepper.h"
class MDIonicStepper : public IonicStepper
{
private:
vector<vector<double> > taum_;
double ekin_;
double th_temp_;
double th_time_;
double eta_;
bool thermostat_;
void stoermer_start(const vector<vector<double> >& fion);
void stoermer_end(const vector<vector<double> >& fion);
public:
MDIonicStepper(Sample& s) : IonicStepper(s)
{
thermostat_ = ( s.ctrl.thermostat == "ON" );
th_temp_ = s.ctrl.th_temp;
th_time_ = s.ctrl.th_time;
eta_ = s.ctrl.pdamp;
taum_.resize(tau0_.size());
for ( int is = 0; is < taum_.size(); is++ )
taum_[is].resize(tau0_[is].size());
}
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);}
double ekin(void) const { return ekin_; }
double temp(void) const;
double eta(void) const { return eta_; }
void set_eta(double x) { eta_ = x; }
};
#endif
////////////////////////////////////////////////////////////////////////////////
//
// MDWavefunctionStepper.h
//
////////////////////////////////////////////////////////////////////////////////
// $Id: MDWavefunctionStepper.h,v 1.1 2003-11-21 20:01:06 fgygi Exp $
#ifndef MDWAVEFUNCTIONSTEPPER_H
#define MDWAVEFUNCTIONSTEPPER_H
class Sample;
class Wavefunction;
#include "WavefunctionStepper.h"
#include "Timer.h"
#include <map>
#include <string>
using namespace std;
typedef map<string,Timer> TimerMap;
class MDWavefunctionStepper : public WavefunctionStepper
{
private:
Sample& s_;
Wavefunction& wf_;
double dt_, dt2bye_;
TimerMap& tmap_;
double ekin_ep_, ekin_em_;
double ekin_eh(void);
public:
void update(Wavefunction& dwf);
void stoermer_start(Wavefunction& dwf);
void stoermer_end(Wavefunction& dwf);
double ekin(void) const { return 0.5*(ekin_ep_ + ekin_em_); }
MDWavefunctionStepper(Sample& s, TimerMap& tmap);
~MDWavefunctionStepper() {};
};
#endif