Commit 3db3b8a0 by Francois Gygi

Restored correct preconditioning when using stress confinement


git-svn-id: http://qboxcode.org/svn/qb/trunk@1315 cba15fb0-1239-40c8-b417-11db7ca47a34
parent f7f2c0e6
......@@ -25,6 +25,7 @@
#include "PSDWavefunctionStepper.h"
#include "PSDAWavefunctionStepper.h"
#include "JDWavefunctionStepper.h"
#include "Preconditioner.h"
#include "SDIonicStepper.h"
#include "SDAIonicStepper.h"
#include "CGIonicStepper.h"
......@@ -210,6 +211,8 @@ void BOSampleStepper::step(int niter)
Timer tm_iter;
Preconditioner prec(wf,ef_,s_.ctrl.ecutprec);
WavefunctionStepper* wf_stepper = 0;
if ( wf_dyn == "SD" )
{
......@@ -225,11 +228,11 @@ void BOSampleStepper::step(int niter)
wf_stepper = new SDWavefunctionStepper(wf,dt2bye,tmap);
}
else if ( wf_dyn == "PSD" )
wf_stepper = new PSDWavefunctionStepper(wf,s_.ctrl.ecutprec,tmap);
wf_stepper = new PSDWavefunctionStepper(wf,prec,tmap);
else if ( wf_dyn == "PSDA" )
wf_stepper = new PSDAWavefunctionStepper(wf,s_.ctrl.ecutprec,tmap);
wf_stepper = new PSDAWavefunctionStepper(wf,prec,tmap);
else if ( wf_dyn == "JD" )
wf_stepper = new JDWavefunctionStepper(wf,s_.ctrl.ecutprec,ef_,tmap);
wf_stepper = new JDWavefunctionStepper(wf,prec,ef_,tmap);
// wf_stepper == 0 indicates that wf_dyn == LOCKED
......
......@@ -23,6 +23,7 @@
#include "SampleStepper.h"
#include "EnergyFunctional.h"
#include "Sample.h"
#include "ChargeDensity.h"
#include "Wavefunction.h"
class WavefunctionStepper;
......
......@@ -68,11 +68,19 @@ void ConfinementPotential::update(void)
// gfgp = G f'(G)
const double gfgp = gsq * fg * fp * sigmas_inv;
if ( ecuts_ > 0.0 )
{
// fstress[ig] = G^2 * f(G)
fstress_[ig] = gsq * fg;
// dfstress = 2 f(G) + G * f'(G)
dfstress_[ig] = 2.0 * fg + gfgp;
}
else
{
fstress_[ig] = 0.0;
dfstress_[ig] = 0.0;
}
// ekin = sum_G |c_G|^2 G^2
// econf = sum_G |c_G|^2 fstress[G]
......
......@@ -26,17 +26,13 @@ using namespace std;
////////////////////////////////////////////////////////////////////////////////
JDWavefunctionStepper::JDWavefunctionStepper(Wavefunction& wf,
double ecutprec, EnergyFunctional& ef, TimerMap& tmap) :
WavefunctionStepper(wf,tmap), prec_(0), wft_(wf), dwft_(wf), ef_(ef)
{
prec_ = new Preconditioner(wf,ecutprec);
}
Preconditioner& prec, EnergyFunctional& ef, TimerMap& tmap) :
WavefunctionStepper(wf,tmap), prec_(prec), wft_(wf), dwft_(wf), ef_(ef)
{}
////////////////////////////////////////////////////////////////////////////////
JDWavefunctionStepper::~JDWavefunctionStepper(void)
{
delete prec_;
}
{}
////////////////////////////////////////////////////////////////////////////////
void JDWavefunctionStepper::update(Wavefunction& dwf)
......@@ -85,7 +81,7 @@ void JDWavefunctionStepper::update(Wavefunction& dwf)
tmap_["jd_residual"].stop();
// dwf.sd->c() now contains the descent direction (HV-VA)
prec_->update(wf_);
prec_.update(wf_);
tmap_["jd_compute_z"].start();
for ( int ispin = 0; ispin < wf_.nspin(); ispin++ )
......@@ -114,7 +110,7 @@ void JDWavefunctionStepper::update(Wavefunction& dwf)
for ( int i = 0; i < ngwl; i++ )
{
const double fac = prec_->diag(ispin,ikp,n,i);
const double fac = prec_.diag(ispin,ikp,n,i);
const double f0 = -fac * dcn[2*i];
const double f1 = -fac * dcn[2*i+1];
cn[2*i] = f0;
......@@ -160,7 +156,7 @@ void JDWavefunctionStepper::update(Wavefunction& dwf)
for ( int i = 0; i < ngwl; i++ )
{
const double fac = prec_->diag(ispin,ikp,n,i);
const double fac = prec_.diag(ispin,ikp,n,i);
cn[i] = -fac * cpn[i];
}
}
......
......@@ -28,7 +28,7 @@ class JDWavefunctionStepper : public WavefunctionStepper
{
private:
Preconditioner *prec_;
Preconditioner& prec_;
EnergyFunctional& ef_;
Wavefunction wft_, dwft_;
......@@ -37,7 +37,7 @@ class JDWavefunctionStepper : public WavefunctionStepper
void update(Wavefunction& dwf);
virtual void preprocess(void) {}
JDWavefunctionStepper(Wavefunction& wf, double ecutprec,
JDWavefunctionStepper(Wavefunction& wf, Preconditioner& prec,
EnergyFunctional& ef, TimerMap& tmap);
~JDWavefunctionStepper();
};
......
......@@ -25,18 +25,14 @@ using namespace std;
////////////////////////////////////////////////////////////////////////////////
PSDAWavefunctionStepper::PSDAWavefunctionStepper(Wavefunction& wf,
double ecutprec, TimerMap& tmap) : prec_(0),
Preconditioner& prec, TimerMap& tmap) : prec_(prec),
WavefunctionStepper(wf,tmap), wf_last_(wf), dwf_last_(wf),
extrapolate_(false)
{
prec_ = new Preconditioner(wf,ecutprec);
}
{}
////////////////////////////////////////////////////////////////////////////////
PSDAWavefunctionStepper::~PSDAWavefunctionStepper(void)
{
delete prec_;
}
{}
////////////////////////////////////////////////////////////////////////////////
void PSDAWavefunctionStepper::update(Wavefunction& dwf)
......@@ -76,7 +72,7 @@ void PSDAWavefunctionStepper::update(Wavefunction& dwf)
// dwf.sd->c() now contains the descent direction (HV-VA) (residual)
// update the preconditioner using the residual
prec_->update(dwf);
prec_.update(dwf);
for ( int ispin = 0; ispin < wf_.nspin(); ispin++ )
{
......@@ -98,7 +94,7 @@ void PSDAWavefunctionStepper::update(Wavefunction& dwf)
double* dcn = &dc[2*mloc*n];
for ( int i = 0; i < ngwl; i++ )
{
const double fac = prec_->diag(ispin,ikp,n,i);
const double fac = prec_.diag(ispin,ikp,n,i);
const double f0 = -fac * dcn[2*i];
const double f1 = -fac * dcn[2*i+1];
dcn[2*i] = f0;
......
......@@ -27,7 +27,7 @@ class PSDAWavefunctionStepper : public WavefunctionStepper
{
private:
Preconditioner *prec_;
Preconditioner& prec_;
Wavefunction wf_last_, dwf_last_;
// Anderson acceleration flag
......@@ -38,7 +38,8 @@ class PSDAWavefunctionStepper : public WavefunctionStepper
void update(Wavefunction& dwf);
virtual void preprocess(void) { extrapolate_ = false; }
PSDAWavefunctionStepper(Wavefunction& wf, double ecutprec, TimerMap& tmap);
PSDAWavefunctionStepper(Wavefunction& wf, Preconditioner& prec,
TimerMap& tmap);
~PSDAWavefunctionStepper();
};
#endif
......@@ -25,16 +25,13 @@ using namespace std;
////////////////////////////////////////////////////////////////////////////////
PSDWavefunctionStepper::PSDWavefunctionStepper(Wavefunction& wf,
double ecutprec, TimerMap& tmap) : WavefunctionStepper(wf,tmap)
{
prec_ = new Preconditioner(wf,ecutprec);
}
Preconditioner& prec, TimerMap& tmap) : WavefunctionStepper(wf,tmap),
prec_(prec)
{}
////////////////////////////////////////////////////////////////////////////////
PSDWavefunctionStepper::~PSDWavefunctionStepper(void)
{
delete prec_;
}
{}
////////////////////////////////////////////////////////////////////////////////
void PSDWavefunctionStepper::update(Wavefunction& dwf)
......@@ -77,7 +74,7 @@ void PSDWavefunctionStepper::update(Wavefunction& dwf)
// dwf.sd->c() now contains the descent direction (HV-VA)
// update preconditioner using the residual
prec_->update(dwf);
prec_.update(dwf);
for ( int ispin = 0; ispin < wf_.nspin(); ispin++ )
{
......@@ -98,7 +95,7 @@ void PSDWavefunctionStepper::update(Wavefunction& dwf)
for ( int i = 0; i < ngwl; i++ )
{
const double fac = prec_->diag(ispin,ikp,n,i);
const double fac = prec_.diag(ispin,ikp,n,i);
const double delta_re = fac * dc[2*i];
const double delta_im = fac * dc[2*i+1];
c[2*i] -= delta_re;
......
......@@ -21,18 +21,20 @@
#include "WavefunctionStepper.h"
class Preconditioner;
class EnergyFunctional;
class PSDWavefunctionStepper : public WavefunctionStepper
{
private:
Preconditioner *prec_;
Preconditioner& prec_;
public:
void update(Wavefunction& dwf);
PSDWavefunctionStepper(Wavefunction& wf, double ecutprec, TimerMap& tmap);
PSDWavefunctionStepper(Wavefunction& wf, Preconditioner& prec,
TimerMap& tmap);
~PSDWavefunctionStepper();
};
#endif
......@@ -19,12 +19,14 @@
#include "Preconditioner.h"
#include "Basis.h"
#include "Wavefunction.h"
#include "EnergyFunctional.h"
#include "ConfinementPotential.h"
#include "SlaterDet.h"
using namespace std;
////////////////////////////////////////////////////////////////////////////////
Preconditioner::Preconditioner(const Wavefunction& wf, double ecutprec)
: ecutprec_(ecutprec)
Preconditioner::Preconditioner(const Wavefunction& wf, EnergyFunctional& ef,
double ecutprec) : ef_(ef), ecutprec_(ecutprec)
{
kpg2_.resize(wf.nspin());
ekin_.resize(wf.nspin());
......@@ -46,12 +48,13 @@ Preconditioner::Preconditioner(const Wavefunction& wf, double ecutprec)
////////////////////////////////////////////////////////////////////////////////
double Preconditioner::diag(int ispin, int ikp, int n, int ig) const
{
const valarray<double>& fstress = ef_.confpot(ikp)->fstress();
if ( ecutprec_ == 0.0 )
{
const double ekin_n = ekin_[ispin][ikp][n];
if ( ekin_n == 0.0 )
return 0.0;
const double q2 = kpg2_[ispin][ikp][ig];
const double q2 = kpg2_[ispin][ikp][ig] + fstress[ig];
#if 0
// Payne Teter Allan adaptive preconditioner
const double x = 0.5 * q2 / ekin_n;
......@@ -69,11 +72,7 @@ double Preconditioner::diag(int ispin, int ikp, int n, int ig) const
else
{
// next lines: inclusion of fstress if confinement potential is used
// const valarray<double>& fstress = ef_.confpot(ikp)->fstress();
// double e = 0.5 * ( kpg2_ptr[ig] + fstress[ig] );
// diag_[ispin][ikp][ig] = ( e < ecutpr ) ? 0.5 / ecutpr : 0.5 / e;
const double e = 0.5 * kpg2_[ispin][ikp][ig];
double e = 0.5 * ( kpg2_[ispin][ikp][ig] + fstress[ig] );
return ( e < ecutprec_ ) ? 0.5 / ecutprec_ : 0.5 / e;
}
}
......
......@@ -20,6 +20,7 @@
#define PRECONDITIONER_H
class Wavefunction;
class EnergyFunctional;
#include <vector>
#include <valarray>
......@@ -28,6 +29,8 @@ class Preconditioner
{
private:
EnergyFunctional& ef_;
// kpg2_[ispin][ikp][ig]
std::vector<std::vector<const double *> > kpg2_;
// ekin_[ispin][ikp][n]
......@@ -42,7 +45,7 @@ class Preconditioner
double diag(int ispin, int ikp, int n, int ig) const;
Preconditioner(const Wavefunction& wf, double ecutprec);
Preconditioner(const Wavefunction& wf, EnergyFunctional& ef, double ecutprec);
//~Preconditioner();
};
#endif
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