rel1_57_8

git-svn-id: http://qboxcode.org/svn/qb/branches/rel1_57@1319 cba15fb0-1239-40c8-b417-11db7ca47a34

src/AngleCmd.h

@@ -33,8 +33,8 @@ class AngleCmd : public Cmd
 
   AngleCmd(Sample *sample) : s(sample) {};
 
-  char *name(void) const { return "angle"; }
-  char *help_msg(void) const
+  const char *name(void) const { return "angle"; }
+  const char *help_msg(void) const
   {
     return
     "\n angle\n\n"

src/AtomCmd.h

@@ -34,8 +34,8 @@ class AtomCmd : public Cmd
 
   AtomCmd(Sample *sample) : s(sample) {};
 
-  char *name(void) const { return "atom"; }
-  char *help_msg(void) const
+  const char *name(void) const { return "atom"; }
+  const char *help_msg(void) const
   {
     return
     "\n atom\n\n"

src/AtomsDyn.h

@@ -33,7 +33,7 @@ class AtomsDyn : public Var
 
   public:
 
-  char *name ( void ) const { return "atoms_dyn"; };
+  const char *name ( void ) const { return "atoms_dyn"; };
 
   int set ( int argc, char **argv )
   {

src/BOSampleStepper.C

@@ -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
 

src/BOSampleStepper.h

@@ -23,6 +23,7 @@
 #include "SampleStepper.h"
 #include "EnergyFunctional.h"
 #include "Sample.h"
+#include "ChargeDensity.h"
 #include "Wavefunction.h"
 
 class WavefunctionStepper;

src/BisectionCmd.h

@@ -35,8 +35,8 @@ class BisectionCmd : public Cmd
 
   BisectionCmd(Sample *sample) : s(sample) {};
 
-  char *name(void) const { return "bisection"; }
-  char *help_msg(void) const
+  const char *name(void) const { return "bisection"; }
+  const char *help_msg(void) const
   {
     return
     "\n bisection\n\n"

src/BlHF.h

@@ -32,7 +32,7 @@ class BlHF : public Var
 
   public:
 
-  char *name ( void ) const { return "blHF"; };
+  const char *name ( void ) const { return "blHF"; };
 
   int set ( int argc, char **argv )
   {

src/BtHF.h

@@ -32,7 +32,7 @@ class BtHF : public Var
 
   public:
 
-  char *name ( void ) const { return "btHF"; };
+  const char *name ( void ) const { return "btHF"; };
 
   int set ( int argc, char **argv )
   {

src/Cell.h

@@ -33,7 +33,7 @@ class Cell : public Var
 
   public:
 
-  char *name ( void ) const { return "cell"; };
+  const char *name ( void ) const { return "cell"; };
 
   int set ( int argc, char **argv )
   {

src/CellDyn.h

@@ -34,7 +34,7 @@ class CellDyn : public Var
 
   public:
 
-  char *name ( void ) const { return "cell_dyn"; };
+  const char *name ( void ) const { return "cell_dyn"; };
 
   int set ( int argc, char **argv )
   {

src/CellLock.h

@@ -33,7 +33,7 @@ class CellLock : public Var
 
   public:
 
-  char *name ( void ) const { return "cell_lock"; };
+  const char *name ( void ) const { return "cell_lock"; };
 
   int set ( int argc, char **argv )
   {

src/CellMass.h

@@ -33,7 +33,7 @@ class CellMass : public Var
 
   public:
 
-  char *name ( void ) const { return "cell_mass"; };
+  const char *name ( void ) const { return "cell_mass"; };
 
   int set ( int argc, char **argv )
   {

src/ChargeMixCoeff.h

@@ -33,7 +33,7 @@ class ChargeMixCoeff : public Var
 
   public:
 
-  char *name ( void ) const { return "charge_mix_coeff"; };
+  const char *name ( void ) const { return "charge_mix_coeff"; };
 
   int set ( int argc, char **argv )
   {

src/ChargeMixNdim.h

@@ -33,7 +33,7 @@ class ChargeMixNdim : public Var
 
   public:
 
-  char *name ( void ) const { return "charge_mix_ndim"; };
+  const char *name ( void ) const { return "charge_mix_ndim"; };
 
   int set ( int argc, char **argv )
   {

src/ChargeMixRcut.h

@@ -33,7 +33,7 @@ class ChargeMixRcut : public Var
 
   public:
 
-  char *name ( void ) const { return "charge_mix_rcut"; };
+  const char *name ( void ) const { return "charge_mix_rcut"; };
 
   int set ( int argc, char **argv )
   {

src/ComputeMLWFCmd.h

@@ -37,9 +37,8 @@ class ComputeMLWFCmd : public Cmd
 
   ComputeMLWFCmd(Sample *sample) : s(sample) {};
 
-  char *name(void) const { return "compute_mlwf"; }
-
-  char *help_msg(void) const
+  const char *name(void) const { return "compute_mlwf"; }
+  const char *help_msg(void) const
   {
     return
     "\n compute_mlwf\n\n"

src/ConfinementPotential.C

@@ -68,11 +68,19 @@ void ConfinementPotential::update(void)
     // gfgp = G f'(G)
     const double gfgp = gsq * fg * fp * sigmas_inv;
 
-    // fstress[ig] = G^2 * f(G)
-    fstress_[ig] = gsq * fg;
+    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;
+      // 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]

src/ConstraintCmd.h

@@ -31,8 +31,8 @@ class ConstraintCmd : public Cmd
 
   ConstraintCmd(Sample *sample) : s(sample) {};
 
-  char *name(void) const { return "constraint"; }
-  char *help_msg(void) const
+  const char *name(void) const { return "constraint"; }
+  const char *help_msg(void) const
   {
     return
     "\n constraint\n\n"

src/Debug.h

@@ -33,7 +33,7 @@ class Debug : public Var
 
   public:
 
-  char *name ( void ) const { return "debug"; };
+  const char *name ( void ) const { return "debug"; };
 
   int set ( int argc, char **argv )
   {

src/DistanceCmd.h

@@ -33,8 +33,8 @@ class DistanceCmd : public Cmd
 
   DistanceCmd(Sample *sample) : s(sample) {};
 
-  char *name(void) const { return "distance"; }
-  char *help_msg(void) const
+  const char *name(void) const { return "distance"; }
+  const char *help_msg(void) const
   {
     return
     "\n distance\n\n"

src/Dspin.h

@@ -32,7 +32,7 @@ class Dspin : public Var
 
   public:
 
-  char *name ( void ) const { return "delta_spin"; };
+  const char *name ( void ) const { return "delta_spin"; };
 
   int set ( int argc, char **argv )
   {

src/Dt.h

@@ -33,7 +33,7 @@ class Dt : public Var
 
   public:
 
-  char *name ( void ) const { return "dt"; };
+  const char *name ( void ) const { return "dt"; };
 
   int set ( int argc, char **argv )
   {

src/Ecut.h

@@ -33,7 +33,7 @@ class Ecut : public Var
 
   public:
 
-  char *name ( void ) const { return "ecut"; };
+  const char *name ( void ) const { return "ecut"; };
 
   int set ( int argc, char **argv )
   {

src/Ecutprec.h

@@ -33,7 +33,7 @@ class Ecutprec : public Var
 
   public:
 
-  char *name ( void ) const { return "ecutprec"; };
+  const char *name ( void ) const { return "ecutprec"; };
 
   int set ( int argc, char **argv )
   {

src/Ecuts.h

@@ -33,7 +33,7 @@ class Ecuts : public Var
 
   public:
 
-  char *name ( void ) const { return "ecuts"; };
+  const char *name ( void ) const { return "ecuts"; };
 
   int set ( int argc, char **argv )
   {

src/Emass.h

@@ -33,7 +33,7 @@ class Emass : public Var
 
   public:
 
-  char *name ( void ) const { return "emass"; };
+  const char *name ( void ) const { return "emass"; };
 
   int set ( int argc, char **argv )
   {

src/EnergyFunctional.C

@@ -20,6 +20,7 @@
 #include "Sample.h"
 #include "Species.h"
 #include "Wavefunction.h"
+#include "ChargeDensity.h"
 #include "SlaterDet.h"
 #include "Basis.h"
 #include "FourierTransform.h"

src/EnergyFunctional.h

@@ -24,7 +24,6 @@
 #include <valarray>
 #include <map>
 #include <string>
-#include "ChargeDensity.h"
 #include "StructureFactor.h"
 #include "Timer.h"
 
@@ -32,6 +31,7 @@ class Sample;
 class Basis;
 class AtomSet;
 class Wavefunction;
+class ChargeDensity;
 class UnitCell;
 class FourierTransform;
 class XCOperator;

src/ExtForceCmd.h

@@ -31,8 +31,8 @@ class ExtForceCmd : public Cmd
 
   ExtForceCmd(Sample *sample) : s(sample) {};
 
-  char *name(void) const { return "extforce"; }
-  char *help_msg(void) const
+  const char *name(void) const { return "extforce"; }
+  const char *help_msg(void) const
   {
     return
     "\n extforce\n\n"

src/ExtStress.h

@@ -33,7 +33,7 @@ class ExtStress : public Var
 
   public:
 
-  char *name ( void ) const { return "ext_stress"; };
+  const char *name ( void ) const { return "ext_stress"; };
 
   int set ( int argc, char **argv )
   {

src/FermiTemp.h

@@ -33,7 +33,7 @@ class FermiTemp : public Var
 
   public:
 
-  char *name ( void ) const { return "fermi_temp"; };
+  const char *name ( void ) const { return "fermi_temp"; };
 
   int set ( int argc, char **argv )
   {

src/FoldInWsCmd.h

@@ -36,8 +36,8 @@ class FoldInWsCmd : public Cmd
 
   FoldInWsCmd(Sample *sample) : s(sample) {};
 
-  char *name(void) const { return "fold_in_ws"; }
-  char *help_msg(void) const
+  const char *name(void) const { return "fold_in_ws"; }
+  const char *help_msg(void) const
   {
     return
     "\n fold_in_ws\n\n"

src/HelpCmd.h

@@ -36,9 +36,8 @@ class HelpCmd : public Cmd
 
   HelpCmd(Sample *sample) : s(sample) {};
 
-  char *name(void) const { return "help"; }
-
-  char *help_msg(void) const
+  const char *name(void) const { return "help"; }
+  const char *help_msg(void) const
   {
     return
     "\n help\n\n"

src/JDWavefunctionStepper.C

@@ -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];
           }
         }

src/JDWavefunctionStepper.h

@@ -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();
 };

src/KpointCmd.h

@@ -32,8 +32,8 @@ class KpointCmd : public Cmd
 
   KpointCmd(Sample *sample) : s(sample) {};
 
-  char *name(void) const { return "kpoint"; }
-  char *help_msg(void) const
+  const char *name(void) const { return "kpoint"; }
+  const char *help_msg(void) const
   {
     return
     "\n kpoint\n\n"

src/ListAtomsCmd.h

@@ -33,8 +33,8 @@ class ListAtomsCmd : public Cmd
 
   ListAtomsCmd(Sample *sample) : s(sample) {};
 
-  char *name(void) const { return "list_atoms"; }
-  char *help_msg(void) const
+  const char *name(void) const { return "list_atoms"; }
+  const char *help_msg(void) const
   {
     return
     "\n list_atoms\n\n"

src/ListSpeciesCmd.h

@@ -33,8 +33,8 @@ class ListSpeciesCmd : public Cmd
 
   ListSpeciesCmd(Sample *sample) : s(sample) {};
 
-  char *name(void) const { return "list_species"; }
-  char *help_msg(void) const
+  const char *name(void) const { return "list_species"; }
+  const char *help_msg(void) const
   {
     return
     "\n list_species\n\n"

src/LoadCmd.h

@@ -36,9 +36,8 @@ class LoadCmd : public Cmd
 
   LoadCmd(Sample *sample) : s(sample) {};
 
-  char *name(void) const { return "load"; }
-
-  char *help_msg(void) const
+  const char *name(void) const { return "load"; }
+  const char *help_msg(void) const
   {
     return
     "\n load\n\n"

src/MDWavefunctionStepper.C

@@ -108,14 +108,10 @@ void MDWavefunctionStepper::compute_wfm(Wavefunction& dwf)
       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();
       const int mloc = sd->c().mloc();
       const int nloc = sd->c().nloc();
-      const bool onrow0 = ( wf_.context().myrow() == 0 );
       for ( int n = 0; n < nloc; n++ )
       {
-        const int nglobal = sd->c().j(0,n);
-        const double occn = occ[nglobal];
         // note: double mloc length for complex<double> indices
         double* c = &cptr[2*mloc*n];
         double* cv = &cptrv[2*mloc*n];

src/MoveCmd.h

@@ -36,8 +36,8 @@ class MoveCmd : public Cmd
 
   MoveCmd(Sample *sample) : s(sample) {};
 
-  char *name(void) const { return "move"; }
-  char *help_msg(void) const
+  const char *name(void) const { return "move"; }
+  const char *help_msg(void) const
   {
     return
     "\n move\n\n"

src/Nempty.h

@@ -33,7 +33,7 @@ class Nempty : public Var
 
   public:
 
-  char *name ( void ) const { return "nempty"; };
+  const char *name ( void ) const { return "nempty"; };
 
   int set ( int argc, char **argv )
   {

src/NetCharge.h

@@ -33,7 +33,7 @@ class NetCharge : public Var
 
   public:
 
-  char *name ( void ) const { return "net_charge"; };
+  const char *name ( void ) const { return "net_charge"; };
 
   int set ( int argc, char **argv )
   {

src/NonLocalPotential.C

@@ -252,7 +252,6 @@ void NonLocalPotential::update_twnl(void)
   const double s3 = sqrt(3.0);
 
   const double *kpg   = basis_.kpg_ptr();
-  const double *kpg2  = basis_.kpg2_ptr();
   const double *kpgi  = basis_.kpgi_ptr();
   const double *kpg_x = basis_.kpgx_ptr(0);
   const double *kpg_y = basis_.kpgx_ptr(1);
@@ -299,9 +298,6 @@ void NonLocalPotential::update_twnl(void)
               const double tgx = kpg_x[ig];
               const double tgy = kpg_y[ig];
               const double tgz = kpg_z[ig];
-              const double tgx2 = tgx * tgx;
-              const double tgy2 = tgy * tgy;
-              const double tgz2 = tgz * tgz;
 
               const double tmp = kpgi[ig] * s14pi * dv;
               dt0_xx[ig] = tmp * tgx * tgx;
@@ -690,7 +686,6 @@ void NonLocalPotential::update_twnl(void)
               const double tgz2 = tgz * tgz;
 
               const double tgi = kpgi[ig];
-              const double tg2 = tg * tg;
               const double tgi2 = tgi * tgi;
 
               const double tgxx = tgx2 * tgi2;
@@ -868,7 +863,6 @@ void NonLocalPotential::update_twnl(void)
                 const double tgz2 = tgz * tgz;
 
                 const double tgi = kpgi[ig];
-                const double tg2 = tg * tg;
                 const double tgi2 = tgi * tgi;
 
                 const double tgxx = tgx2 * tgi2;

src/Nrowmax.h

@@ -33,7 +33,7 @@ class Nrowmax : public Var
 
   public:
 
-  char *name ( void ) const { return "nrowmax"; };
+  const char *name ( void ) const { return "nrowmax"; };
 
   int set ( int argc, char **argv )
   {

src/Nspin.h

@@ -32,7 +32,7 @@ class Nspin : public Var
 
   public:
 
-  char *name ( void ) const { return "nspin"; };
+  const char *name ( void ) const { return "nspin"; };
 
   int set ( int argc, char **argv )
   {

src/PSDAWavefunctionStepper.C

@@ -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;

src/PSDAWavefunctionStepper.h

@@ -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

src/PSDWavefunctionStepper.C

@@ -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;

src/PSDWavefunctionStepper.h

@@ -21,18 +21,20 @@
 
 #include "WavefunctionStepper.h"
 class Preconditioner;
+class EnergyFunctional;