Commit ae688b96 by Francois Gygi

removed Sample.h dependency


git-svn-id: http://qboxcode.org/svn/qb/trunk@482 cba15fb0-1239-40c8-b417-11db7ca47a34
parent 56cc7d07
......@@ -3,7 +3,7 @@
// BOSampleStepper.C
//
////////////////////////////////////////////////////////////////////////////////
// $Id: BOSampleStepper.C,v 1.29 2006-11-04 20:19:41 fgygi Exp $
// $Id: BOSampleStepper.C,v 1.30 2007-01-27 23:46:30 fgygi Exp $
#include "BOSampleStepper.h"
#include "EnergyFunctional.h"
......@@ -105,11 +105,22 @@ void BOSampleStepper::step(int niter)
WavefunctionStepper* wf_stepper = 0;
if ( wf_dyn == "SD" )
wf_stepper = new SDWavefunctionStepper(s_,tmap);
{
const double emass = s_.ctrl.emass;
double dt2bye = (emass == 0.0) ? 0.5 / wf.ecut() : dt*dt/emass;
// divide dt2bye by facs coefficient if stress == ON
const double facs = 2.0;
if ( s_.ctrl.stress == "ON" )
{
dt2bye /= facs;
}
wf_stepper = new SDWavefunctionStepper(wf,dt2bye,tmap);
}
else if ( wf_dyn == "PSD" )
wf_stepper = new PSDWavefunctionStepper(s_,*preconditioner,tmap);
wf_stepper = new PSDWavefunctionStepper(wf,*preconditioner,tmap);
else if ( wf_dyn == "PSDA" )
wf_stepper = new PSDAWavefunctionStepper(s_,*preconditioner,tmap);
wf_stepper = new PSDAWavefunctionStepper(wf,*preconditioner,tmap);
// wf_stepper == 0 indicates that wf_dyn == LOCKED
......@@ -639,13 +650,42 @@ void BOSampleStepper::step(int niter)
{
energy = ef_.energy(true,dwf,false,fion,false,sigma_eks);
s_.wf.diag(dwf,compute_eigvec);
if ( onpe0 )
{
// print eigenvalues
for ( int ispin = 0; ispin < wf.nspin(); ispin++ )
{
for ( int ikp = 0; ikp < wf.nkp(); ikp++ )
{
if ( wf.sd(ispin,ikp) != 0 )
{
if ( wf.sdcontext(ispin,ikp)->active() )
{
const int nst = wf.sd(ispin,ikp)->nst();
const double eVolt = 2.0 * 13.6058;
cout << " <eigenvalues spin=\"" << ispin
<< "\" kpoint=\"" << wf.sd(ispin,ikp)->kpoint()
<< "\" n=\"" << nst << "\">" << endl;
for ( int i = 0; i < nst; i++ )
{
cout << setw(12) << setprecision(5)
<< wf.sd(ispin,ikp)->eig(i)*eVolt;
if ( i%5 == 4 ) cout << endl;
}
if ( nst%5 != 0 ) cout << endl;
cout << " </eigenvalues>" << endl;
}
}
}
}
}
}
// update occupation numbers if fractionally occupied states
if ( fractional_occ )
{
s_.wf.update_occ(s_.ctrl.fermi_temp);
const double wf_entropy = s_.wf.entropy();
wf.update_occ(s_.ctrl.fermi_temp);
const double wf_entropy = wf.entropy();
if ( onpe0 )
{
cout << " <!-- Wavefunction entropy: " << wf_entropy
......
......@@ -3,7 +3,7 @@
// CPSampleStepper.C
//
////////////////////////////////////////////////////////////////////////////////
// $Id: CPSampleStepper.C,v 1.11 2005-09-16 23:08:11 fgygi Exp $
// $Id: CPSampleStepper.C,v 1.12 2007-01-27 23:46:31 fgygi Exp $
#include "CPSampleStepper.h"
#include "SlaterDet.h"
......@@ -21,7 +21,17 @@ using namespace std;
CPSampleStepper::CPSampleStepper(Sample& s) :
SampleStepper(s), cd_(s.wf), ef_(s,cd_), dwf(s.wf), wfv(s.wfv)
{
mdwf_stepper = new MDWavefunctionStepper(s,tmap);
const double emass = s.ctrl.emass;
const double dt = s.ctrl.dt;
double dt2bye = (emass == 0.0) ? 0.5 / s.wf.ecut() : dt*dt/emass;
// divide dt2bye by facs coefficient if stress == ON
const double facs = 2.0;
if ( s.ctrl.stress == "ON" )
{
dt2bye /= facs;
}
mdwf_stepper = new MDWavefunctionStepper(s.wf,s.wfv,dt,dt2bye,tmap);
assert(mdwf_stepper!=0);
mdionic_stepper = 0;
if ( s.ctrl.atoms_dyn != "LOCKED" )
......
......@@ -3,7 +3,7 @@
// EnergyFunctional.C
//
////////////////////////////////////////////////////////////////////////////////
// $Id: EnergyFunctional.C,v 1.22 2005-06-27 22:25:04 fgygi Exp $
// $Id: EnergyFunctional.C,v 1.23 2007-01-27 23:46:31 fgygi Exp $
#include "EnergyFunctional.h"
#include "Sample.h"
......@@ -138,8 +138,10 @@ EnergyFunctional::EnergyFunctional(const Sample& s, const ChargeDensity& cd)
create_cfp |= wf.sd(ispin,ikp) != 0 && wf.sdcontext(ispin,ikp)->active();
if ( create_cfp )
{
const double facs = 2.0;
const double sigmas = 0.5;
cfp[ikp] =
new ConfinementPotential(s_.ctrl.ecuts,s_.ctrl.facs,s_.ctrl.sigmas,
new ConfinementPotential(s_.ctrl.ecuts,facs,sigmas,
wf.sd(0,ikp)->basis());
}
}
......
......@@ -3,7 +3,7 @@
// MDWavefunctionStepper.C
//
////////////////////////////////////////////////////////////////////////////////
// $Id: MDWavefunctionStepper.C,v 1.6 2005-04-29 18:13:48 fgygi Exp $
// $Id: MDWavefunctionStepper.C,v 1.7 2007-01-27 23:46:31 fgygi Exp $
#include "MDWavefunctionStepper.h"
#include "Wavefunction.h"
......@@ -13,24 +13,17 @@
using namespace std;
////////////////////////////////////////////////////////////////////////////////
MDWavefunctionStepper::MDWavefunctionStepper(Sample& s, TimerMap& tmap) :
WavefunctionStepper(s,tmap)
MDWavefunctionStepper::MDWavefunctionStepper(Wavefunction& wf,
Wavefunction *wfv, double dt, double dt2bye, TimerMap& tmap) :
wfv_(wfv), dt_(dt), dt2bye_(dt2bye), WavefunctionStepper(wf,tmap)
{
dt_ = s_.ctrl.dt;
const double emass = s_.ctrl.emass;
dt2bye_ = (emass == 0.0) ? 0.5 / wf_.ecut() : dt_*dt_/emass;
// divide dt2bye by facs coefficient if confinement is on (ecuts>0)
if ( s_.ctrl.ecuts > 0.0 )
dt2bye_ /= s_.ctrl.facs;
assert(wfv!=0);
}
////////////////////////////////////////////////////////////////////////////////
void MDWavefunctionStepper::update(Wavefunction& dwf)
{
// Verlet update of wf using force dwf and wfm stored in *wfv
Wavefunction* const wfv = s_.wfv;
assert(wfv!=0);
for ( int ispin = 0; ispin < wf_.nspin(); ispin++ )
{
for ( int ikp = 0; ikp < wf_.nkp(); ikp++ )
......@@ -48,7 +41,7 @@ void MDWavefunctionStepper::update(Wavefunction& dwf)
// c = cp
SlaterDet* sd = wf_.sd(ispin,ikp);
double* cptr = (double*) sd->c().valptr();
double* cptrm = (double*) wfv->sd(ispin,ikp)->c().valptr();
double* cptrm = (double*) wfv_->sd(ispin,ikp)->c().valptr();
const double* dcptr =
(const double*) dwf.sd(ispin,ikp)->c().cvalptr();
const int mloc = sd->c().mloc();
......@@ -79,7 +72,7 @@ void MDWavefunctionStepper::update(Wavefunction& dwf)
tmap_["md_update_wf"].stop();
tmap_["riccati"].start();
wf_.sd(ispin,ikp)->riccati(*(wfv->sd(ispin,ikp)));
wf_.sd(ispin,ikp)->riccati(*(wfv_->sd(ispin,ikp)));
tmap_["riccati"].stop();
}
}
......@@ -95,8 +88,6 @@ void MDWavefunctionStepper::compute_wfm(Wavefunction& dwf)
// Compute wfm for first MD step using wf, wfv and dwf (= Hpsi)
// Replace then wfv by wfm
Wavefunction* const wfv = s_.wfv;
assert(wfv!=0);
// Compute cm using c and wavefunction velocity
// cm = c - dt * v - 0.5 * dt2/m * hpsi
// replace wfv by wfm
......@@ -112,7 +103,7 @@ void MDWavefunctionStepper::compute_wfm(Wavefunction& dwf)
SlaterDet* sd = wf_.sd(ispin,ikp);
double* cptr = (double*) sd->c().valptr();
double* cptrv = (double*) wfv->sd(ispin,ikp)->c().valptr();
double* cptrv = (double*) wfv_->sd(ispin,ikp)->c().valptr();
const double* dcptr =
(const double*) dwf.sd(ispin,ikp)->c().cvalptr();
const vector<double>& occ = sd->occ();
......@@ -140,7 +131,7 @@ void MDWavefunctionStepper::compute_wfm(Wavefunction& dwf)
}
}
tmap_["riccati"].start();
wfv->sd(ispin,ikp)->riccati(*wf_.sd(ispin,ikp));
wfv_->sd(ispin,ikp)->riccati(*wf_.sd(ispin,ikp));
tmap_["riccati"].stop();
}
}
......@@ -150,7 +141,7 @@ void MDWavefunctionStepper::compute_wfm(Wavefunction& dwf)
ekin_ep_ = ekin_eh();
// Note: ekin_ep is a first-order approximation of ekin_e using wf and wfm
// only. This makes ekin_e consistent with the following update() call
// Note: *wfv now contains wf(t-dt)
// Note: *wfv_ now contains wf(t-dt)
}
////////////////////////////////////////////////////////////////////////////////
......@@ -158,8 +149,6 @@ void MDWavefunctionStepper::compute_wfv(Wavefunction& dwf)
{
// Compute wfv = (wf - wfm)/dt - 0.5*dtbye*dwf
Wavefunction * const wfv = s_.wfv;
assert(wfv!=0);
assert(dt_!=0.0);
for ( int ispin = 0; ispin < wf_.nspin(); ispin++ )
{
......@@ -172,7 +161,7 @@ void MDWavefunctionStepper::compute_wfv(Wavefunction& dwf)
// compute final velocity wfv
// v = ( c - cm ) / dt - 0.5 * dt/m * hpsi
// Note: At this point, *wfv contains wf(t-dt)
// Note: At this point, *wfv_ contains wf(t-dt)
// hpsi must be orthogonal to the subspace spanned by c
// compute descent direction H psi - psi (psi^T H psi)
......@@ -203,7 +192,7 @@ void MDWavefunctionStepper::compute_wfv(Wavefunction& dwf)
const double dt_inv = 1.0/dt_;
const double half_dtbye = 0.5 * dt2bye_ / dt_;
double* cptr = (double*) sd->c().valptr();
double* cptrv = (double*) wfv->sd(ispin,ikp)->c().valptr();
double* cptrv = (double*) wfv_->sd(ispin,ikp)->c().valptr();
const double* dcptr =
(const double*) dwf.sd(ispin,ikp)->c().cvalptr();
const int mloc = sd->c().mloc();
......@@ -229,7 +218,7 @@ void MDWavefunctionStepper::compute_wfv(Wavefunction& dwf)
- half_dtbye * dctmp1;
}
}
// Note: *wfv now contains the wavefunction velocity
// Note: *wfv_ now contains the wavefunction velocity
}
}
}
......@@ -242,7 +231,7 @@ double MDWavefunctionStepper::ekin_eh(void)
// compute ekin at time t - 0.5*dt using wf and wfm
tmap_["ekin_e"].start();
double ekin_e = 0.0;
Wavefunction* const wfv = s_.wfv; // assume that wfv contains wfm
// assume that wfv contains wfm
for ( int ispin = 0; ispin < wf_.nspin(); ispin++ )
{
for ( int ikp = 0; ikp < wf_.nkp(); ikp++ )
......@@ -253,7 +242,7 @@ double MDWavefunctionStepper::ekin_eh(void)
{
SlaterDet* sd = wf_.sd(ispin,ikp);
double* cptr = (double*) sd->c().valptr();
double* cptrm = (double*) wfv->sd(ispin,ikp)->c().valptr();
double* cptrm = (double*) wfv_->sd(ispin,ikp)->c().valptr();
const int mloc = sd->c().mloc();
const int nloc = sd->c().nloc();
// compute electronic kinetic energy at time t-1/2
......
......@@ -3,7 +3,7 @@
// MDWavefunctionStepper.h
//
////////////////////////////////////////////////////////////////////////////////
// $Id: MDWavefunctionStepper.h,v 1.3 2004-03-11 21:52:31 fgygi Exp $
// $Id: MDWavefunctionStepper.h,v 1.4 2007-01-27 23:46:31 fgygi Exp $
#ifndef MDWAVEFUNCTIONSTEPPER_H
#define MDWAVEFUNCTIONSTEPPER_H
......@@ -14,7 +14,9 @@ class MDWavefunctionStepper : public WavefunctionStepper
{
private:
double dt_, dt2bye_;
double dt_;
double dt2bye_;
Wavefunction *wfv_;
double ekin_ep_, ekin_em_;
double ekin_eh(void);
......@@ -26,7 +28,8 @@ class MDWavefunctionStepper : public WavefunctionStepper
void compute_wfv(Wavefunction& dwf);
double ekin(void) const { return 0.5*(ekin_ep_ + ekin_em_); }
MDWavefunctionStepper(Sample& s, TimerMap& tmap);
MDWavefunctionStepper(Wavefunction& wf, Wavefunction* wfv,
double dt, double dt2bye, TimerMap& tmap);
~MDWavefunctionStepper() {};
};
#endif
......@@ -3,20 +3,19 @@
// PSDAWavefunctionStepper.C
//
////////////////////////////////////////////////////////////////////////////////
// $Id: PSDAWavefunctionStepper.C,v 1.10 2004-12-18 23:20:17 fgygi Exp $
// $Id: PSDAWavefunctionStepper.C,v 1.11 2007-01-27 23:46:31 fgygi Exp $
#include "PSDAWavefunctionStepper.h"
#include "Wavefunction.h"
#include "SlaterDet.h"
#include "Sample.h"
#include "Preconditioner.h"
#include <iostream>
using namespace std;
////////////////////////////////////////////////////////////////////////////////
PSDAWavefunctionStepper::PSDAWavefunctionStepper(Sample& s,
PSDAWavefunctionStepper::PSDAWavefunctionStepper(Wavefunction& wf,
Preconditioner& p, TimerMap& tmap) :
WavefunctionStepper(s,tmap), prec_(p), wf_last_(s.wf), dwf_last_(s.wf),
WavefunctionStepper(wf,tmap), prec_(p), wf_last_(wf), dwf_last_(wf),
extrapolate_(false)
{}
......
......@@ -3,12 +3,13 @@
// PSDAWavefunctionStepper.h
//
////////////////////////////////////////////////////////////////////////////////
// $Id: PSDAWavefunctionStepper.h,v 1.3 2004-03-11 21:52:31 fgygi Exp $
// $Id: PSDAWavefunctionStepper.h,v 1.4 2007-01-27 23:46:31 fgygi Exp $
#ifndef PSDAWAVEFUNCTIONSTEPPER_H
#define PSDAWAVEFUNCTIONSTEPPER_H
#include "WavefunctionStepper.h"
#include "Wavefunction.h"
class Preconditioner;
class PSDAWavefunctionStepper : public WavefunctionStepper
......@@ -26,7 +27,7 @@ class PSDAWavefunctionStepper : public WavefunctionStepper
void update(Wavefunction& dwf);
virtual void preprocess(void) { extrapolate_ = false; }
PSDAWavefunctionStepper(Sample& s, Preconditioner& p, TimerMap& tmap);
PSDAWavefunctionStepper(Wavefunction& wf, Preconditioner& p, TimerMap& tmap);
~PSDAWavefunctionStepper() {};
};
#endif
......@@ -3,20 +3,19 @@
// PSDWavefunctionStepper.C
//
////////////////////////////////////////////////////////////////////////////////
// $Id: PSDWavefunctionStepper.C,v 1.6 2004-11-10 22:35:23 fgygi Exp $
// $Id: PSDWavefunctionStepper.C,v 1.7 2007-01-27 23:46:31 fgygi Exp $
#include "PSDWavefunctionStepper.h"
#include "Wavefunction.h"
#include "SlaterDet.h"
#include "Sample.h"
#include "Preconditioner.h"
#include <iostream>
using namespace std;
////////////////////////////////////////////////////////////////////////////////
PSDWavefunctionStepper::PSDWavefunctionStepper(Sample& s,
PSDWavefunctionStepper::PSDWavefunctionStepper(Wavefunction& wf,
Preconditioner& p, TimerMap& tmap) :
WavefunctionStepper(s,tmap), prec_(p)
WavefunctionStepper(wf,tmap), prec_(p)
{}
////////////////////////////////////////////////////////////////////////////////
......
......@@ -3,7 +3,7 @@
// PSDWavefunctionStepper.h
//
////////////////////////////////////////////////////////////////////////////////
// $Id: PSDWavefunctionStepper.h,v 1.3 2004-03-11 21:52:32 fgygi Exp $
// $Id: PSDWavefunctionStepper.h,v 1.4 2007-01-27 23:46:31 fgygi Exp $
#ifndef PSDWAVEFUNCTIONSTEPPER_H
#define PSDWAVEFUNCTIONSTEPPER_H
......@@ -21,7 +21,7 @@ class PSDWavefunctionStepper : public WavefunctionStepper
void update(Wavefunction& dwf);
PSDWavefunctionStepper(Sample& s, Preconditioner& p, TimerMap& tmap);
PSDWavefunctionStepper(Wavefunction& wf, Preconditioner& p, TimerMap& tmap);
~PSDWavefunctionStepper() {};
};
#endif
......@@ -3,7 +3,7 @@
// SDWavefunctionStepper.C
//
////////////////////////////////////////////////////////////////////////////////
// $Id: SDWavefunctionStepper.C,v 1.3 2004-11-10 22:35:23 fgygi Exp $
// $Id: SDWavefunctionStepper.C,v 1.4 2007-01-27 23:46:31 fgygi Exp $
#include "SDWavefunctionStepper.h"
#include "Wavefunction.h"
......@@ -13,17 +13,10 @@
using namespace std;
////////////////////////////////////////////////////////////////////////////////
SDWavefunctionStepper::SDWavefunctionStepper(Sample& s, TimerMap& tmap) :
WavefunctionStepper(s,tmap)
{
dt_ = s_.ctrl.dt;
const double emass = s_.ctrl.emass;
dt2bye_ = (emass == 0.0) ? 0.5 / wf_.ecut() : dt_*dt_/emass;
// divide dt2bye by facs coefficient if stress == ON
if ( s_.ctrl.stress == "ON" )
dt2bye_ /= s_.ctrl.facs;
}
SDWavefunctionStepper::SDWavefunctionStepper(Wavefunction& wf, double alpha,
TimerMap& tmap) :
alpha_(alpha), WavefunctionStepper(wf,tmap)
{}
////////////////////////////////////////////////////////////////////////////////
void SDWavefunctionStepper::update(Wavefunction& dwf)
......@@ -38,7 +31,7 @@ void SDWavefunctionStepper::update(Wavefunction& dwf)
{
// c = c - dt2bye * hpsi
tmap_["sd_update_wf"].start();
wf_.sd(ispin,ikp)->c().axpy(-dt2bye_,dwf.sd(ispin,ikp)->c());
wf_.sd(ispin,ikp)->c().axpy(-alpha_,dwf.sd(ispin,ikp)->c());
tmap_["sd_update_wf"].stop();
tmap_["gram"].start();
wf_.sd(ispin,ikp)->gram();
......
......@@ -3,7 +3,7 @@
// SDWavefunctionStepper.h
//
////////////////////////////////////////////////////////////////////////////////
// $Id: SDWavefunctionStepper.h,v 1.2 2004-02-04 19:55:16 fgygi Exp $
// $Id: SDWavefunctionStepper.h,v 1.3 2007-01-27 23:46:31 fgygi Exp $
#ifndef SDWAVEFUNCTIONSTEPPER_H
#define SDWAVEFUNCTIONSTEPPER_H
......@@ -14,13 +14,13 @@ class SDWavefunctionStepper : public WavefunctionStepper
{
private:
double dt_, dt2bye_;
double alpha_;
public:
void update(Wavefunction& dwf);
SDWavefunctionStepper(Sample& s, TimerMap& tmap);
SDWavefunctionStepper(Wavefunction& wf, double alpha, TimerMap& tmap);
~SDWavefunctionStepper() {};
};
#endif
......@@ -3,11 +3,10 @@
// WavefunctionStepper.h
//
////////////////////////////////////////////////////////////////////////////////
// $Id: WavefunctionStepper.h,v 1.5 2004-03-11 21:52:31 fgygi Exp $
// $Id: WavefunctionStepper.h,v 1.6 2007-01-27 23:46:31 fgygi Exp $
#ifndef WAVEFUNCTIONSTEPPER_H
#define WAVEFUNCTIONSTEPPER_H
#include "Sample.h"
#include "Timer.h"
#include <map>
#include <string>
......@@ -21,7 +20,6 @@ class WavefunctionStepper
private:
protected:
Sample& s_;
Wavefunction& wf_;
TimerMap& tmap_;
......@@ -31,8 +29,8 @@ class WavefunctionStepper
virtual void preprocess(void) {}
virtual void postprocess(void) {}
WavefunctionStepper(Sample& s, TimerMap& tmap) :
s_(s), wf_(s.wf), tmap_(tmap)
WavefunctionStepper(Wavefunction& wf, TimerMap& tmap) :
wf_(wf), tmap_(tmap)
{}
virtual ~WavefunctionStepper() {}
};
......
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