fgygi-2016-09-09 updated implementation of response

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

src/BOSampleStepper.C

@@ -46,7 +46,7 @@ using namespace std;
 BOSampleStepper::BOSampleStepper(Sample& s, int nitscf, int nite) :
   SampleStepper(s), cd_(s.wf), ef_(s,cd_),
   dwf(s.wf), wfv(s.wfv), nitscf_(nitscf), nite_(nite),
-  initial_atomic_density(false), initial_density(false), first_step(true) {}
+  update_density_first_(true), update_vxc_(true) {}
 
 ////////////////////////////////////////////////////////////////////////////////
 BOSampleStepper::~BOSampleStepper()
@@ -155,21 +155,8 @@ void BOSampleStepper::initialize_density(void)
       cd_.rhog[1][i] = 0.5 * rhopst[i];
     }
   }
-  initial_atomic_density = true;
-}
-
-////////////////////////////////////////////////////////////////////////////////
-void BOSampleStepper:: initialize_density(const vector<vector<double> >& rhor)
-{
-  // initialize cd_ with a given density rhor
-  assert( rhor.size() == cd_.rhor.size() );
-  for ( int ispin = 0; ispin < rhor.size(); ispin++ )
-  {
-    assert( rhor[ispin].size() == cd_.rhor[ispin].size() );
-    for ( int ir = 0; ir < rhor[ispin].size(); ir++ )
-      cd_.rhor[ispin][ir] = rhor[ispin][ir];
-  }
-  initial_density = true;
+  cd_.update_rhor();
+  update_density_first_ = false;
 }
 
 ////////////////////////////////////////////////////////////////////////////////
@@ -762,13 +749,12 @@ void BOSampleStepper::step(int niter)
         if ( nite_ > 0 && onpe0 )
           cout << "  BOSampleStepper: start scf iteration" << endl;
 
-        // compute new density in cd_.rhog
+        // update charge density
         tmap["charge"].start();
-        if ( itscf==0 && initial_atomic_density )
-          cd_.update_rhor();
-        else if ( itscf==0 && initial_density )
-          cd_.update_rhog();
-        else
+        // The density is updated at the first step if update_density_first_
+        // is true.
+        // It is always updated after the first step
+        if ( ( update_density_first_ || itscf>0 ) )
           cd_.update_density();
         tmap["charge"].stop();
 
@@ -851,21 +837,24 @@ void BOSampleStepper::step(int niter)
           }
         } // if nite_ > 0
 
+        // update vhxc:
+        // at first scf step:
+        // - update both vh and vxc
+        // at later steps:
+        // - update vh
+        // - update vxc only if update_vxc_ is true
         tmap["update_vhxc"].start();
-        bool freeze_vh;
-        bool freeze_vxc;
-        if ( itscf == 0 && first_step )
-        {
-          // at first SCF iteration, vhxc must be updated
-          freeze_vh = false;
-          freeze_vxc = false;
-        }
+        if ( itscf == 0 )
+          ef_.update_vhxc(compute_stress);
         else
         {
-          freeze_vh = s_.ctrl.freeze_vh;
-          freeze_vxc = s_.ctrl.freeze_vxc;
+          if ( update_vxc_ )
+            ef_.update_vhxc(compute_stress);
+          else
+            // update vh only
+            ef_.update_vh(compute_stress);
         }
-        ef_.update_vhxc(compute_stress, freeze_vh, freeze_vxc);
+
         tmap["update_vhxc"].stop();
 
         // reset stepper only if multiple non-selfconsistent steps
@@ -1128,7 +1117,7 @@ void BOSampleStepper::step(int niter)
         tmap["charge"].stop();
 
         tmap["update_vhxc"].start();
-        ef_.update_vhxc(compute_stress,s_.ctrl.freeze_vh,s_.ctrl.freeze_vxc);
+        ef_.update_vhxc(compute_stress);
         tmap["update_vhxc"].stop();
         const bool compute_forces = true;
         tmap["energy"].start();
@@ -1327,7 +1316,5 @@ void BOSampleStepper::step(int niter)
 
   if ( ntc_extrapolation || asp_extrapolation ) delete wfmm;
 
-  initial_atomic_density = false;
-  initial_density = false;
-  first_step = false;
+  update_density_first_ = false;
 }

src/BOSampleStepper.h

@@ -42,9 +42,8 @@ class BOSampleStepper : public SampleStepper
   WavefunctionStepper* wf_stepper;
   IonicStepper* ionic_stepper;
 
-  bool initial_atomic_density;
-  bool initial_density;
-  bool first_step;  // true if step() function has never been called before
+  bool update_density_first_;
+  bool update_vxc_;
 
   // Do not allow construction of BOSampleStepper unrelated to a Sample
   BOSampleStepper(void);
@@ -54,15 +53,16 @@ class BOSampleStepper : public SampleStepper
   mutable TimerMap tmap;
 
   void step(int niter);
-  void initialize_density(void);         // initialize density by atomic density
-  void initialize_density(const vector<vector<double> >& rhor); // initialize density by given density
+  // initialize density with sum of atomic densities
+  void initialize_density(void);
+  void set_update_vxc(bool update_vxc) { update_vxc_ = update_vxc; }
+  void set_update_density_first(bool update_density_first)
+    { update_density_first_ = update_density_first; }
 
   EnergyFunctional& ef(void) { return ef_; }
   ChargeDensity& cd(void) { return cd_; }
 
   BOSampleStepper(Sample& s, int nitscf, int nite);
-  BOSampleStepper(Sample& s, const vector<vector<double> > * rhor_initial,
-                  int nitscf, int nite);
   ~BOSampleStepper();
 };
 #endif

src/Control.h

@@ -75,8 +75,5 @@ struct Control
 
   std::string iter_cmd;
   int iter_cmd_period;
-
-  bool freeze_vh;  // only valid at GS-only calculation
-  bool freeze_vxc;  // only valid at GS-only calculation
 };
 #endif

src/EnergyFunctional.C

@@ -42,7 +42,7 @@
 using namespace std;
 
 ////////////////////////////////////////////////////////////////////////////////
-EnergyFunctional::EnergyFunctional( Sample& s, ChargeDensity& cd)
+EnergyFunctional::EnergyFunctional(Sample& s, ChargeDensity& cd)
  : s_(s), cd_(cd)
 {
   const AtomSet& atoms = s_.atoms;
@@ -66,9 +66,11 @@ EnergyFunctional::EnergyFunctional( Sample& s, ChargeDensity& cd)
   int np2v = vft->np2();
 
   v_r.resize(wf.nspin());
+  vxc_r.resize(wf.nspin());
   for ( int ispin = 0; ispin < wf.nspin(); ispin++ )
   {
     v_r[ispin].resize(vft->np012loc());
+    vxc_r[ispin].resize(vft->np012loc());
   }
   tmp_r.resize(vft->np012loc());
 
@@ -224,8 +226,7 @@ EnergyFunctional::~EnergyFunctional(void)
 }
 
 ////////////////////////////////////////////////////////////////////////////////
-void EnergyFunctional::update_vhxc(bool compute_stress,
-                                   bool freeze_vh, bool freeze_vxc)
+void EnergyFunctional::update_vhxc(bool compute_stress, bool update_vxc)
 {
   // called when the charge density has changed
   // update Hartree and xc potentials using the charge density cd_
@@ -259,10 +260,9 @@ void EnergyFunctional::update_vhxc(bool compute_stress,
   // update XC operator
   // compute xc energy, update self-energy operator and potential
   tmap["exc"].start();
-  for ( int ispin = 0; ispin < wf.nspin(); ispin++ )
-    memset((void*)&v_r[ispin][0], 0, vft->np012loc()*sizeof(double));
 
-  xco->update(v_r, compute_stress, freeze_vxc);
+  if ( update_vxc )
+    xco->update(vxc_r, compute_stress);
 
   exc_ = xco->exc();
   dxc_ = xco->dxc();
@@ -274,7 +274,7 @@ void EnergyFunctional::update_vhxc(bool compute_stress,
     // compute Fourier coefficients of Vxc
     for ( int ispin = 0; ispin < wf.nspin(); ispin++ )
     {
-      copy(v_r[ispin].begin(),v_r[ispin].end(),tmp_r.begin());
+      copy(vxc_r[ispin].begin(),vxc_r[ispin].end(),tmp_r.begin());
       if ( ispin == 0 )
       {
         vft->forward(&tmp_r[0],&vxc_g[0]);
@@ -346,20 +346,18 @@ void EnergyFunctional::update_vhxc(bool compute_stress,
   // add external potential vext to tmp_r
   if ( s_.vext )
   {
-    //if( s_.ctxt_.onpe0() )
-    //  cout << "  EnergyFunctional::update_vhxc: vext is applied " << s_.vext->filename() << endl;
     for ( int i = 0; i < tmp_r.size(); i++ )
       tmp_r[i] += s_.vext->v(i);
   }
 
-  // add local potential in tmp_r to v_r[ispin][i]
-  // v_r contains the xc potential
+  // compute local potential v_r[ispin][i]
+  // vxc_r contains the xc potential
   const int size = tmp_r.size();
   if ( wf.nspin() == 1 )
   {
     for ( int i = 0; i < size; i++ )
     {
-      v_r[0][i] += real(tmp_r[i]);
+      v_r[0][i] = vxc_r[0][i] + real(tmp_r[i]);
     }
   }
   else
@@ -367,8 +365,8 @@ void EnergyFunctional::update_vhxc(bool compute_stress,
     for ( int i = 0; i < size; i++ )
     {
       const double vloc = real(tmp_r[i]);
-      v_r[0][i] += vloc;
-      v_r[1][i] += vloc;
+      v_r[0][i] = vxc_r[0][i] + vloc;
+      v_r[1][i] = vxc_r[1][i] + vloc;
     }
   }
 

src/EnergyFunctional.h

@@ -79,6 +79,7 @@ class EnergyFunctional
   public:
 
   std::vector<std::vector<double> > v_r;
+  std::vector<std::vector<double> > vxc_r;
   mutable TimerMap tmap;
 
   double energy(bool compute_hpsi, Wavefunction& dwf,
@@ -109,8 +110,11 @@ class EnergyFunctional
 
   const ConfinementPotential *confpot(int ikp) const { return cfp[ikp]; }
 
-  void update_vhxc(bool compute_stress,
-                   bool freeze_vh=false, bool freeze_vxc=false);
+  void update_vhxc(bool compute_stress, bool update_vxc);
+  // update both vh and vxc
+  void update_vhxc(bool compute_stress) { update_vhxc(compute_stress, true); }
+  // update only vh
+  void update_vh(bool compute_stress) { update_vhxc(compute_stress, false); }
 
   void atoms_moved(void);
   void cell_moved(void);

src/ExternalPotential.C

-
 ////////////////////////////////////////////////////////////////////////////////
 //
 // Copyright (c) 2016 The Regents of the University of California
@@ -26,33 +25,47 @@ using namespace std;
 #include "ExternalPotential.h"
 #include "FourierTransform.h"
 
-void readRealSpaceFunction(const Context& ctxt, vector<double>& fr,
-                           const FourierTransform& ft,
-                           const string filename, const bool parallel);
-
 void ExternalPotential::update(const ChargeDensity& cd)
 {
-  const Context& ctxt = *(s_.wf.spincontext());
+  const Context &cd_ctxt = cd.context();
+  MPI_Comm vcomm = cd.vcomm();
+  int mype, vcomm_size;
+  MPI_Comm_size(vcomm,&vcomm_size);
+  MPI_Comm_rank(vcomm,&mype);
+  FourierTransform *vft = cd.vft();
+  vector<double> vext_read;
 
-  // process 0 reads the grid size and bcast to all processes
-  if ( ctxt.onpe0() )
+  // All tasks with myrow==0 read the potential
+  if ( mype == 0 )
   {
     ifstream vfile(filename_.c_str());
     if ( vfile )
-      cout << "  ExternalPotential::update: read external "
-           << "potential from file: " << filename_ << endl;
+    {
+      vfile >> n_[0] >> n_[1] >> n_[2];
+      int n012 = n_[0] * n_[1] * n_[2];
+      vext_read.resize(n012);
+      for ( int i = 0; i < n012; i++ )
+         vfile >> vext_read[i];
+      vfile.close();
+    }
     else
     {
-      cout << "  ExternalPotential::update: file not found: "
-           << filename_ << endl;
-      ctxt.abort(1);
+      if ( cd_ctxt.onpe0() )
+        cout << "  ExternalPotential::update: file not found: "
+             << filename_ << endl;
+      cd.context().abort(1);
     }
-    vfile >> n_[0] >> n_[1] >> n_[2];
+  }
+  if ( cd_ctxt.onpe0() )
+  {
+    cout << "  ExternalPotential::update: read external "
+         << "potential from file: " << filename_ << endl;
     cout << "  ExternalPotential::update: grid size "
          << n_[0] << " " << n_[1] << " " << n_[2] << endl;
-    vfile.close();
   }
-  MPI_Bcast(&n_[0],3,MPI_INT,0,cd.vcomm());
+
+  // broadcast sizes to lower rows
+  MPI_Bcast(&n_[0],3,MPI_INT,0,vcomm);
 
   // create a Basis compatible with the vext grid read from file
   const Basis* vbasis = cd.vbasis();
@@ -70,7 +83,7 @@ void ExternalPotential::update(const ChargeDensity& cd)
   double ecut2 = 0.5 * norm2(b2) * n2max*n2max;
   ecut_ = min(min(ecut0,ecut1),ecut2);
 
-  if ( ctxt.onpe0() )
+  if ( cd_ctxt.onpe0() )
   {
     cout << "  ExternalPotential::update: ecut0: " << ecut0 << endl;
     cout << "  ExternalPotential::update: ecut1: " << ecut1 << endl;
@@ -78,35 +91,47 @@ void ExternalPotential::update(const ChargeDensity& cd)
     cout << "  ExternalPotential::update: ecut:  " << ecut_ << endl;
   }
 
-  Basis basis(cd.vcomm(),D3vector(0,0,0));
+  Basis basis(vcomm,D3vector(0,0,0));
   basis.resize(cell,cell,ecut_);
-  if ( ctxt.onpe0() )
+  if ( cd_ctxt.onpe0() )
   {
     cout << "  ExternalPotential::update: np0: " << basis.np(0) << endl;
     cout << "  ExternalPotential::update: np1: " << basis.np(1) << endl;
     cout << "  ExternalPotential::update: np2: " << basis.np(2) << endl;
   }
 
-  // build FT to transform vext to g space
+  assert(basis.np(0)<=vft->np0());
+  assert(basis.np(1)<=vft->np1());
+  assert(basis.np(2)<=vft->np2());
+
+  // FourierTransform for vext grid
   FourierTransform ft(basis,n_[0],n_[1],n_[2]);
 
-  // read vext from file
+  // scatter parts of vext_read to lower rows
+  vector<double> vext_read_loc(ft.np012loc());
+  vector<int> scounts(vcomm_size,0);
+  vector<int> sdispls(vcomm_size,0);
+  int displ = 0;
+  for ( int iproc = 0; iproc < vcomm_size; iproc++ )
+  {
+    sdispls[iproc] = displ;
+    scounts[iproc] = ft.np012loc(iproc);
+    displ += ft.np012loc(iproc);
+  }
 
-  vector<double> vext_r_read;
-  bool parallel_read = false; // later parallel_read should be determined by vext file
-  readRealSpaceFunction(ctxt,vext_r_read,ft,filename_, parallel_read);
+  MPI_Scatterv(&vext_read[0],&scounts[0],&sdispls[0],MPI_DOUBLE,
+               &vext_read_loc[0],ft.np012loc(),MPI_DOUBLE,0,vcomm);
 
   vector<complex<double> > tmp_r(ft.np012loc());
   for ( int ir = 0; ir < tmp_r.size(); ir++ )
-    tmp_r[ir] = complex<double>(vext_r_read[ir],0);
+    tmp_r[ir] = complex<double>(vext_read_loc[ir],0);
 
   // compute Fourier coefficients
   vector<complex<double> > vext_g_(basis.localsize());
   ft.forward(&tmp_r[0],&vext_g_[0]);
   // vext_g_ now contains the Fourier coefficients of vext
 
-  // interpolate to vft grid of charge density
-  FourierTransform *vft = cd.vft();
+  // interpolate to vft grid
   FourierTransform ft2(basis,vft->np0(),vft->np1(),vft->np2());
   tmp_r.resize(vft->np012loc());
   vext_r_.resize(tmp_r.size());
@@ -114,97 +139,3 @@ void ExternalPotential::update(const ChargeDensity& cd)
   for ( int i = 0; i < vext_r_.size(); i++ )
     vext_r_[i] = real(tmp_r[i]);
 }
-
-void readRealSpaceFunction(const Context& ctxt, vector<double>& fr,
-                           const FourierTransform& ft,
-                           const string filename, const bool parallel)
-{
-  // This function read a real space function from a file and distribute
-  // it to fr variable among all process of the context
-  // ft is used to determine how fr is stored at each process
-  // If parallel is true, file will be read in parallel assuming base64 encoding,
-  // otherwise file will be read by process 0 in text format,
-  // and then distribute to other processes
-
-  // generate column and row communicators
-  MPI_Comm comm = ctxt.comm();
-  MPI_Comm comm_col;
-  MPI_Comm comm_row;
-  MPI_Comm_split(comm,ctxt.mycol(),ctxt.myrow(),&comm_col);
-  MPI_Comm_split(comm,ctxt.myrow(),ctxt.mycol(),&comm_row);
-
-  if ( ctxt.onpe0())
-  {
-    assert( ctxt.mycol() == 0 );
-    fr.resize(ft.np012());
-  }
-  else
-    fr.resize(ft.np012loc());
-
-  if ( parallel )
-  {
-    if ( ctxt.mycol() == 0 )
-    {
-      MPI_Info info;
-      MPI_File vfile;
-      int err;
-      err = MPI_File_open(comm_col,filename.c_str(),MPI_MODE_RDONLY,info,&vfile);
-      if ( err != 0 )
-      {
-        cout << " Error opening file " << filename << endl;
-        ctxt.abort(1);
-      }
-
-      // todo: implement MPI read
-
-      MPI_File_close(&vfile);
-    }
-  } // if parallel
-  else
-  {
-    if ( ctxt.onpe0() )
-    {
-      ifstream vfile(filename.c_str());
-      int n1, n2, n3;
-      vfile >> n1 >> n2 >> n3;
-      for ( int ir = 0; ir < fr.size(); ir++ )
-        vfile >> fr[ir];
-      vfile.close();
-    }
-
-    void* sbuf;
-    void* rbuf;
-    if ( ctxt.onpe0() )
-    {
-      sbuf = &fr[0];
-      rbuf = MPI_IN_PLACE;
-    }
-    else
-      rbuf = &fr[0];
-
-    int displ = 0;
-    vector<int> scounts(ctxt.nprow(),0);
-    vector<int> displs(ctxt.nprow(),0);
-    for ( int iproc = 0; iproc < ctxt.nprow(); iproc++ )
-    {
-      displs[iproc] = displ;
-      displ += ft.np012loc(iproc);
-      scounts[iproc] = ft.np012loc(iproc);
-    }
-
-    // process 0 scatter fr to other processes on first column
-    MPI_Barrier(comm);
-    if ( ctxt.mycol() == 0 )
-    {
-      MPI_Scatterv(sbuf,&scounts[0],&displs[0],MPI_DOUBLE,
-                   rbuf,ft.np012loc(),MPI_DOUBLE,0,comm_col);
-    }
-    if ( ctxt.onpe0() )
-    {
-      fr.resize(ft.np012loc());
-    }
-
-    // first column bcast fr to other columns
-    MPI_Bcast(&fr[0],ft.np012loc(),MPI_DOUBLE,0,comm_row);
-  } // if parallel
-}

src/ResponseCmd.h

@@ -27,18 +27,14 @@ class ResponseCmd : public Cmd
 {
   private:
 
-  void responseVext(int nitscf, int nite, bool parallel_output);
+  void responseVext(bool rpa, int nitscf, int nite);
   void responseEfield(double amplitude, int nitscf, int nite);
 
-
   public:
 
   Sample *s;
 
-  ResponseCmd(Sample *sample) : s(sample) {
-    s->ctrl.freeze_vh = false;
-    s->ctrl.freeze_vxc = false;
-  };
+  ResponseCmd(Sample *sample) : s(sample) {};
 
   const char *name(void) const { return "response"; }
   const char *help_msg(void) const
@@ -46,7 +42,7 @@ class ResponseCmd : public Cmd
     return
     "\n response\n\n"
     " syntax: response amplitude nitscf [nite]\n"
-    "         response -vext vext_file [-RPA or -IPA] [-amplitude amplitude] [-parallel_write] nitscf [nite]\n\n"
+    "         response -vext vext_file [-RPA] [-amplitude a] nitscf [nite]\n\n"
     "   The response command computes the polarizability tensor by\n"
     "   finite differences using external electric fields in the x,y,z\n"
     "   directions with magnitude defined by the amplitude argument.\n"

src/Vext.h

@@ -49,20 +49,12 @@ class Vext : public Var
       // set vext NULL
       delete s->vext;
       s->vext = 0;
-#ifdef DEBUG
-      if ( ui->onpe0() )
-        cout << " vext set to zero" << endl;
-#endif
     }
     else
     {
       if ( s->vext )
         delete s->vext;
       s->vext = new ExternalPotential(*s,argv[1]);
-#ifdef DEBUG
-      if ( ui->onpe0() )
-        cout << " vext set to" << argv[1] << endl;
-#endif
     }
 
     return 0;

src/XCOperator.C

@@ -94,8 +94,7 @@ XCOperator::~XCOperator()
 }
 
 ////////////////////////////////////////////////////////////////////////////////
-void XCOperator::update(std::vector<std::vector<double> >& vr,
-                        bool compute_stress, bool freeze_vxc)
+void XCOperator::update(std::vector<std::vector<double> >& vr, bool compute_stress)
 {
   // update xc potential and self-energy
   // used whenever the charge density and/or wave functions have changed
@@ -103,7 +102,7 @@ void XCOperator::update(std::vector<std::vector<double> >& vr,
   if ( hasPotential_ )
   {
     // update LDA/GGA xc potential
-    xcp_->update(vr, freeze_vxc);
+    xcp_->update( vr );
 
     // LDA/GGA exchange energy
     exc_ = xcp_->exc();

src/XCOperator.h

@@ -60,8 +60,7 @@ class XCOperator
   bool hasGGA(void) { return hasGGA_; };
   bool hasHF(void) { return hasHF_; };
 
-  void update(std::vector<std::vector<double> >& vr,
-              bool compute_stress, bool freeze_vxc);
+  void update(std::vector<std::vector<double> >& vr, bool compute_stress);
   void apply_self_energy(Wavefunction &dwf);
   void compute_stress(std::valarray<double>& sigma);
   void cell_moved(void);

src/XCPotential.C

@@ -91,9 +91,9 @@ bool XCPotential::isGGA(void)
 }
 
 ////////////////////////////////////////////////////////////////////////////////
-void XCPotential::update(vector<vector<double> >& vr, bool freeze_vxc)
+void XCPotential::update(vector<vector<double> >& vr)
 {
-  // compute exchange-correlation energy and add vxc potential to vr[ispin][ir]
+  // compute exchange-correlation energy and vxc potential vr[ispin][ir]
 
   // Input: total electronic density in:
   //   vector<vector<double> >           cd_.rhor[ispin][ir] (real space)
@@ -122,17 +122,7 @@ void XCPotential::update(vector<vector<double> >& vr, bool freeze_vxc)
   {
     // LDA functional
 
-    if ( ! freeze_vxc )
-    {
-      xcf_->setxc();
-      //if ( cd_.context().onpe0() )
-      //  cout << "  XCPotential::update(LDA): vxc is updated" << endl;
-    }
-    else
-    {
-      if ( cd_.context().onpe0() )
-        cout << "  XCPotential::update(LDA): vxc is freezed" << endl;
-    }
+    xcf_->setxc();
 
     exc_ = 0.0;
     dxc_ = 0.0;
@@ -151,7 +141,7 @@ void XCPotential::update(vector<vector<double> >& vr, bool freeze_vxc)
         const double rh_i = rh[i];
         exc_ += rh_i * e_i;
         dxc_ += rh_i * ( e_i - v_i );
-        vr[0][i] += v_i;
+        vr[0][i] = v_i;
       }
     }
     else
@@ -166,8 +156,8 @@ void XCPotential::update(vector<vector<double> >& vr, bool freeze_vxc)
         const double r_i = rh_up[i] + rh_dn[i];
         exc_ += r_i * e[i];
         dxc_ += r_i * e[i] - rh_up[i] * v_up[i] - rh_dn[i] * v_dn[i];
-        vr[0][i] += v_up[i];
-        vr[1][i] += v_dn[i];
+        vr[0][i] = v_up[i];
+        vr[1][i] = v_dn[i];
       }
     }
     double sum[2],tsum[2];
@@ -180,8 +170,6 @@ void XCPotential::update(vector<vector<double> >& vr, bool freeze_vxc)
   else
   {
     // GGA functional
-   if ( ! freeze_vxc )
-   {
     exc_ = 0.0;
 
     // compute grad_rho
@@ -229,14 +217,6 @@ void XCPotential::update(vector<vector<double> >& vr, bool freeze_vxc)
     }
 
     xcf_->setxc();
-    //if ( cd_.context().onpe0() )
-    //  cout << "  XCPotential::update(GGA): vxc is updated" << endl;
-    }
-    else
-    {
-      if ( cd_.context().onpe0() )
-         cout << "  XCPotential::update(GGA): vxc is freezed" << endl;
-    }
 
     // compute xc potential
     // take divergence of grad(rho)*vxc2
@@ -319,7 +299,7 @@ void XCPotential::update(vector<vector<double> >& vr, bool freeze_vxc)
       } // j
     }
 
-    // add xc potential to local potential in vr[i]
+    // xc potential vr[i]
     // div(vxc2*grad_rho) is stored in vxctmp[ispin][ir]
 
     double esum=0.0;
@@ -337,7 +317,7 @@ void XCPotential::update(vector<vector<double> >& vr, bool freeze_vxc)
           const double v_i = v1[ir] + vxctmp[0][ir];
           esum += rh_i * e_i;
           dsum += rh_i * ( e_i - v_i );
-          vr[0][ir] += v_i;
+          vr[0][ir] = v_i;
         }
       }
     }
@@ -357,8 +337,8 @@ void XCPotential::update(vector<vector<double> >& vr, bool freeze_vxc)
         const double v_up = v1_up[ir] + vxctmp[0][ir];
         const double v_dn = v1_dn[ir] + vxctmp[1][ir];
         dsum += r_up_i * ( eup[ir] - v_up ) + r_dn_i * ( edn[ir] - v_dn );
-        vr[0][ir] += v_up;
-        vr[1][ir] += v_dn;
+        vr[0][ir] = v_up;
+        vr[1][ir] = v_dn;
       }
     }
     double sum[2], tsum[2];

src/XCPotential.h

@@ -55,7 +55,7 @@ class XCPotential
   XCPotential(const ChargeDensity& cd, const std::string functional_name,
     const Control& ctrl);
   ~XCPotential();
-  void update(std::vector<std::vector<double> >& vr, bool freeze_vxc);
+  void update(std::vector<std::vector<double> >& vr);
   void compute_stress(std::valarray<double>& sigma_exc);
   double exc(void) { return exc_; }
   double dxc(void) { return dxc_; }