Merge branch 'vext' into develop

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) {}
+  update_density_first_(true), update_vh_(true), update_vxc_(true) {}
 
 ////////////////////////////////////////////////////////////////////////////////
 BOSampleStepper::~BOSampleStepper()
@@ -158,7 +158,8 @@ void BOSampleStepper::initialize_density(void)
       cd_.rhog[1][i] = 0.5 * rhopst[i];
     }
   }
-  initial_atomic_density = true;
+  cd_.update_rhor();
+  update_density_first_ = false;
 }
 
 ////////////////////////////////////////////////////////////////////////////////
@@ -753,11 +754,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
+        // 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();
 
@@ -843,8 +845,16 @@ void BOSampleStepper::step(int niter)
           }
         } // if nite_ > 0
 
+        // update vhxc:
+        // at first scf step:
+        // - update both vh and vxc
+        // at later steps:
+        // - update depending of values of update_vh_ and update_vxc_
         tmap["update_vhxc"].start();
-        ef_.update_vhxc(compute_stress);
+        if ( itscf == 0 )
+          ef_.update_vhxc(compute_stress);
+        else
+          ef_.update_vhxc(compute_stress, update_vh_, update_vxc_);
         tmap["update_vhxc"].stop();
 
         // reset stepper only if multiple non-selfconsistent steps
@@ -1294,13 +1304,6 @@ void BOSampleStepper::step(int niter)
     ionic_stepper->compute_v(energy,fion);
     // positions r0 and velocities v0 are consistent
   }
-  else
-  {
-    // delete wavefunction velocities
-    if ( s_.wfv != 0 )
-      delete s_.wfv;
-    s_.wfv = 0;
-  }
 
   for ( int ispin = 0; ispin < nspin; ispin++ )
   {
@@ -1314,5 +1317,5 @@ void BOSampleStepper::step(int niter)
 
   if ( ntc_extrapolation || asp_extrapolation ) delete wfmm;
 
-  initial_atomic_density = false;
+  update_density_first_ = true;
 }

src/BOSampleStepper.h

@@ -42,7 +42,8 @@ class BOSampleStepper : public SampleStepper
   WavefunctionStepper* wf_stepper;
   IonicStepper* ionic_stepper;
 
-  bool initial_atomic_density;
+  bool update_density_first_;
+  bool update_vh_, update_vxc_;
 
   // Do not allow construction of BOSampleStepper unrelated to a Sample
   BOSampleStepper(void);
@@ -52,7 +53,15 @@ class BOSampleStepper : public SampleStepper
   mutable TimerMap tmap;
 
   void step(int niter);
+  // initialize density with sum of atomic densities
   void initialize_density(void);
+  void set_update_vh(bool update_vh) { update_vh_ = update_vh; }
+  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();

src/Basis.h

@@ -55,7 +55,7 @@ class Basis
 
   std::vector<int>    idx_;   // 3-d index of vectors idx[i*3+j]
   std::vector<double> g_;     // norm of g vectors g[localsize]
-  std::vector<double> kpg_;   // norm of g vectors g[localsize]
+  std::vector<double> kpg_;   // norm of k+g vectors kpg[localsize]
   std::vector<double> gi_;    // inverse norm of g vectors gi[localsize]
   std::vector<double> kpgi_;  // inverse norm of k+g vectors kpgi[localsize]
   std::vector<double> g2_;    // 2-norm of g vectors g2[localsize]

src/ChargeDensity.C

@@ -26,6 +26,7 @@
 #include <iomanip>
 #include <algorithm> // fill
 #include <functional>
+#include <fstream>
 using namespace std;
 
 ////////////////////////////////////////////////////////////////////////////////
@@ -250,3 +251,34 @@ ostream& operator<< ( ostream& os, const ChargeDensity& cd )
   cd.print(os);
   return os;
 }
+
+////////////////////////////////////////////////////////////////////////////////
+void ChargeDensity::update_rhog(void)
+{
+  // recalculate rhog from rhor
+  assert(rhor.size() == wf_.nspin());
+  const double omega = vbasis_->cell().volume();
+  assert(omega!=0.0);
+
+  for ( int ispin = 0; ispin < wf_.nspin(); ispin++ )
+  {
+    const int rhor_size = rhor[ispin].size();
+    double *const prhor = &rhor[ispin][0];
+    if ( rhocore_r )
+    {
+      #pragma omp parallel for
+      for ( int i = 0; i < rhor_size; i++ )
+        rhotmp[i] = complex<double> ( omega * (prhor[i] - rhocore_r[i]), 0);
+    }
+    else
+    {
+      #pragma omp parallel for
+      for ( int i = 0; i < rhor_size; i++ )
+        rhotmp[i] = complex<double> ( omega * prhor[i], 0);
+    }
+
+    assert(rhotmp.size() == vft_->np012loc() );
+
+    vft_->forward(&rhotmp[0],&rhog[ispin][0]);
+  }
+}

src/ChargeDensity.h

@@ -26,6 +26,7 @@
 #include <map>
 #include "Timer.h"
 #include "Context.h"
+#include "D3vector.h"
 
 class Wavefunction;
 class FourierTransform;
@@ -56,6 +57,7 @@ class ChargeDensity
   double* rhocore_r;
   void update_density(void);
   void update_rhor(void);
+  void update_rhog(void);
 
   const Context& context(void) const { return ctxt_; }
   MPI_Comm vcomm(void) const { return vcomm_; }

src/ConstraintSet.C

@@ -23,6 +23,7 @@
 #include "TorsionConstraint.h"
 #include "Atom.h"
 #include "AtomSet.h"
+#include <cstring>
 #include "Context.h"
 #include <iostream>
 #include <iomanip>

src/ElectricEnthalpy.C

@@ -731,6 +731,13 @@ void ElectricEnthalpy::print(ostream& os) const
 }
 
 ////////////////////////////////////////////////////////////////////////////////
+void ElectricEnthalpy::set_e_field(D3vector e_field_val)
+{
+  e_field_ = e_field_val;
+  finite_field_ = norm2(e_field_) != 0.0;
+}
+
+////////////////////////////////////////////////////////////////////////////////
 ostream& operator<< ( ostream& os, const ElectricEnthalpy& e )
 {
   e.print(os);

src/ElectricEnthalpy.h

@@ -90,6 +90,7 @@ class ElectricEnthalpy
 
   double enthalpy(Wavefunction& dwf, bool compute_hpsi);
 
+  void set_e_field(D3vector e_field_val);
   void update(void);
   void print(std::ostream& os) const;
 

src/EnergyFunctional.C

@@ -30,6 +30,7 @@
 #include "ConfinementPotential.h"
 #include "D3vector.h"
 #include "ElectricEnthalpy.h"
+#include "ExternalPotential.h"
 
 #include "Timer.h"
 #include "blas.h"
@@ -41,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;
@@ -65,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());
 
@@ -80,6 +83,10 @@ EnergyFunctional::EnergyFunctional( Sample& s, ChargeDensity& cd)
          << vft->np012() << endl;
   }
 
+  // external potential
+  if ( s_.vext )
+    s_.vext->update(cd_);
+
   const int ngloc = vbasis_->localsize();
 
   nsp_ = atoms.nsp();
@@ -220,7 +227,8 @@ EnergyFunctional::~EnergyFunctional(void)
 }
 
 ////////////////////////////////////////////////////////////////////////////////
-void EnergyFunctional::update_vhxc(bool compute_stress)
+void EnergyFunctional::update_vhxc(bool compute_stress, bool update_vh,
+  bool update_vxc)
 {
   // called when the charge density has changed
   // update Hartree and xc potentials using the charge density cd_
@@ -254,10 +262,10 @@ 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);
+  if ( update_vxc )
+    xco->update(vxc_r, compute_stress);
+
   exc_ = xco->exc();
   dxc_ = xco->dxc();
   if ( compute_stress )
@@ -268,7 +276,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]);
@@ -309,7 +317,7 @@ void EnergyFunctional::update_vhxc(bool compute_stress)
     ehesum += norm(r) * g2i[ig];
     ehepsum += 2.0*real(conj(r)*rp * g2i[ig]);
     ehpsum += norm(rp) * g2i[ig];
-    rhogt[ig] = r+rp;
+    if ( update_vh ) rhogt[ig] = r+rp;
   }
   // factor 1/2 from definition of Ehart cancels with half sum over G
   // Note: rhogt[ig] includes a factor 1/Omega
@@ -337,14 +345,23 @@ void EnergyFunctional::update_vhxc(bool compute_stress)
   // FT to tmpr_r
   vft->backward(&vlocal_g[0],&tmp_r[0]);
 
-  // add local potential in tmp_r to v_r[ispin][i]
-  // v_r contains the xc potential
+  // add external potential vext to tmp_r
+  if ( s_.vext )
+  {
+    for ( int i = 0; i < tmp_r.size(); i++ )
+      tmp_r[i] += s_.vext->v(i);
+    // update eext_
+    eext_ = s_.vext->compute_eext(cd_);
+  }
+
+  // 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
@@ -352,8 +369,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;
     }
   }
 
@@ -677,6 +694,8 @@ double EnergyFunctional::energy(bool compute_hpsi, Wavefunction& dwf,
     ets_ = - wf_entropy * s_.ctrl.fermi_temp * boltz;
   }
   etotal_ = ekin_ + econf_ + eps_ + enl_ + ecoul_ + exc_ + ets_ + eexf_;
+  if ( s_.vext )
+    etotal_ += eext_;
   enthalpy_ = etotal_;
 
   // Electric enthalpy
@@ -1187,6 +1206,8 @@ void EnergyFunctional::print(ostream& os) const
      << "  <epv>     " << setw(15) << epv() << " </epv>\n"
      << "  <eefield> " << setw(15) << eefield() << " </eefield>\n"
      << "  <enthalpy>" << setw(15) << enthalpy() << " </enthalpy>" << endl;
+  if ( s_.vext )
+    os << "  <eext>    " << setw(15) << eext() << " </eext>" << endl;
 }
 
 ////////////////////////////////////////////////////////////////////////////////

src/EnergyFunctional.h

@@ -68,7 +68,7 @@ class EnergyFunctional
   std::vector<int> na_;
   int nsp_;
   double ekin_, econf_, eps_, enl_, ehart_, ehart_e_, ehart_ep_, ehart_p_,
-         ecoul_, exc_, esr_, eself_, ets_, eexf_, etotal_;
+         ecoul_, exc_, esr_, eself_, ets_, eexf_, eext_, etotal_;
   double dxc_;
   double epv_, eefield_, enthalpy_;
   std::valarray<double> sigma_ekin,sigma_econf,sigma_eps,sigma_ehart,sigma_exc,
@@ -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,
@@ -104,12 +105,18 @@ class EnergyFunctional
   double eefield(void) const { return eefield_; }
   double epv(void) const { return epv_; }
   double enthalpy(void) const { return enthalpy_; }
+  double eext(void) const { return eext_; }
 
   ElectricEnthalpy* el_enth() { return el_enth_; }
 
   const ConfinementPotential *confpot(int ikp) const { return cfp[ikp]; }
 
-  void update_vhxc(bool compute_stress);
+  void update_vhxc(bool compute_stress, bool update_vh, bool update_vxc);
+  // update both vh and vxc
+  void update_vhxc(bool compute_stress)
+  {
+    update_vhxc(compute_stress, true, true);
+  }
 
   void atoms_moved(void);
   void cell_moved(void);

src/ExternalPotential.C

+////////////////////////////////////////////////////////////////////////////////
+//
+// Copyright (c) 2016 The Regents of the University of California
+//
+// This file is part of Qbox
+//
+// Qbox is distributed under the terms of the GNU General Public License
+// as published by the Free Software Foundation, either version 2 of
+// the License, or (at your option) any later version.
+// See the file COPYING in the root directory of this distribution
+// or <http://www.gnu.org/licenses/>.
+//
+////////////////////////////////////////////////////////////////////////////////
+//
+// ExternalPotential.C
+//
+////////////////////////////////////////////////////////////////////////////////
+
+#include <iostream>
+#include <fstream>
+#include <cassert>
+using namespace std;
+
+#include "Basis.h"
+#include "ExternalPotential.h"
+#include "FourierTransform.h"
+#include "Function3d.h"
+#include "Base64Transcoder.h"
+
+////////////////////////////////////////////////////////////////////////////////
+bool abs_compare(const double &a, const double &b)
+{
+  return (abs(a) < abs(b));
+}
+
+////////////////////////////////////////////////////////////////////////////////
+void ExternalPotential::update(const ChargeDensity& cd)
+{
+  const Context* ctxt = s_.wf.spincontext();
+  bool onpe0 = ctxt->onpe0();
+  int nprow = ctxt->nprow();
+  int myrow = ctxt->myrow();
+  int mycol = ctxt->mycol();
+  MPI_Comm vcomm = cd.vcomm();
+
+  Timer tm_read_vext;
+  double time, tmin, tmax;
+
+  // In cube mode and xml mode, the external potential is
+  // read by processors on row 0 and stored in vext_read
+  vector<double> vext_read, vext_read_loc;
+
+  tm_read_vext.start();
+  if ( fmt_ == "cube" )
+  {
+    // read cube file, n_'s are determined by cube file
+    if ( myrow == 0 )
+    {
+      ifstream vfile(filename_.c_str());
+      if (!vfile)
+      {
+        if (mycol == 0)
+          cout << "  ExternalPotential::update: file not found: "
+               << filename_ << endl;
+        ctxt->abort(1);
+      }
+      string tmpstr;
+      for (int i = 0; i < 2; i++)
+        getline(vfile, tmpstr);  // skip comments
+      int natom;
+      vfile >> natom;
+      getline(vfile, tmpstr);
+      for (int i = 0; i < 3; i++)
+      {
+        vfile >> n_[i];
+        getline(vfile, tmpstr);
+      }
+      for (int i = 0; i < natom; i++)
+        getline(vfile, tmpstr);  // skip atom coordinates
+      vext_read.resize(n_[0] * n_[1] * n_[2]);
+      for (int nx = 0; nx < n_[0]; nx++)
+        for (int ny = 0; ny < n_[1]; ny++)
+          for (int nz = 0; nz < n_[2]; nz++)
+          {
+            const int ir = nx + ny * n_[0] + nz * n_[0] * n_[1];
+            vfile >> vext_read[ir];
+          }
+      vfile.close();
+    }
+    MPI_Bcast(&n_[0],3,MPI_INT,0,vcomm);
+  }
+  else if ( fmt_ == "xml" )
+  {
+    if (myrow == 0)
+    {
+      Function3d f;
+      f.read(filename_);
+      vext_read = f.val;
+      n_[0] = f.nx;
+      n_[1] = f.ny;
+      n_[2] = f.nz;
+    }
+    MPI_Bcast(&n_[0],3,MPI_INT,0,vcomm);
+  }
+  tm_read_vext.stop();
+  // at this point, all processes have correct n_ regardless of io mode
+
+  // now construct 2 Fourier transforms ft1 and ft2
+  // ft1 is used to transform vext_read_loc to vext_g (G space)
+  // ft2 is used to transform vext_g to vext_r_ (R space)
+  // the whole process is a Fourier interpolation/extrapolation
+  // of the external potential on the charge density grid
+
+  // create a Basis with largest possible ecut that is compatible with
+  // the external potential grid from file
+  const UnitCell& cell = cd.vbasis()->cell();
+  int n0max = (n_[0]-2)/2;
+  int n1max = (n_[1]-2)/2;
+  int n2max = (n_[2]-2)/2;
+  double ecut0 = 0.5 * norm2(cell.b(0)) * n0max*n0max;
+  double ecut1 = 0.5 * norm2(cell.b(1)) * n1max*n1max;
+  double ecut2 = 0.5 * norm2(cell.b(2)) * n2max*n2max;
+  ecut_ = min(min(ecut0,ecut1),ecut2);
+  Basis basis(vcomm,D3vector(0,0,0));
+  basis.resize(cell,cell,ecut_);
+
+  FourierTransform *vft = cd.vft();
+
+  FourierTransform ft1(basis,n_[0],n_[1],n_[2]);
+  vext_read_loc.resize(ft1.np012loc());
+  vector<complex<double> > vext_g(basis.localsize());
+
+  // check that the basis fits in the vft grid
+  //assert(basis.fits_in_grid(vft->np0(),vft->np1(),vft->np2()));
+
+  FourierTransform ft2(basis,vft->np0(),vft->np1(),vft->np2());
+  vext_r_.resize(ft2.np012loc());
+
+  // xml mode or cube mode: processors on row 0 scatter
+  // vext to other rows
+  vector<int> scounts(nprow,0);
+  vector<int> sdispls(nprow,0);
+  int displ = 0;
+  for ( int iproc = 0; iproc < nprow; iproc++ )
+  {
+    sdispls[iproc] = displ;
+    scounts[iproc] = ft1.np012loc(iproc);
+    displ += ft1.np012loc(iproc);
+  }
+  MPI_Scatterv(&vext_read[0],&scounts[0],&sdispls[0],MPI_DOUBLE,
+               &vext_read_loc[0],ft1.np012loc(),MPI_DOUBLE,0,vcomm);
+
+  // now vext_read_loc on all processors contains the correct portion of vext
+  // Fourier forward transform vext_read_loc to vext_g
+  vector<complex<double> > tmp_r(ft1.np012loc());
+  for ( int ir = 0; ir < tmp_r.size(); ir++ )
+    tmp_r[ir] = complex<double>(vext_read_loc[ir],0);
+  ft1.forward(&tmp_r[0],&vext_g[0]);
+  // Fourier backward transform vext_g to vext_r_
+  tmp_r.resize(ft2.np012loc());
+  ft2.backward(&vext_g[0],&tmp_r[0]);
+  for ( int i = 0; i < vext_r_.size(); i++ )
+    vext_r_[i] = real(tmp_r[i]);
+  if ( onpe0 )
+  {
+    cout << "  ExternalPotential::update: read external potential from file:  "
+         << filename_ << endl;
+    cout << "  ExternalPotential::update: grid size n0 = "
+         << n_[0] << ", n1 = " << n_[1] << ", n2 = " << n_[2] << endl;
+    cout << "  ExternalPotential::update: ecut:  " << ecut_ << endl;
+    if ( amplitude_ != 0 )
+      cout << "  ExternalPotential::update: amplitude:  " << amplitude_ << endl;
+  }
+  if ( amplitude_ == 0.0 )
+  {
+    // If amplitude_ = 0.0, use following scheme to get an amplitude.
+    // Empirically, an absolute magnitude of 1.0E-3 ~ 1.0E-5 hartree for Vext
+    // would be suitable. Here the amplitude is set to scale the
+    // magnitude of vext to be 1.0E-3 hartree
+    if (vext_read_loc.size() > 0)
+      magnitude_ = abs(*max_element(vext_read_loc.begin(),
+                       vext_read_loc.end(), abs_compare));
+    ctxt->dmax('C',1,1,&magnitude_,1);
+    MPI_Bcast(&magnitude_,1,MPI_DOUBLE,0,vcomm);
+    amplitude_ = 1.0E-3 / magnitude_;
+    if ( onpe0 )
+      cout << "  ExternalPotential::update: amplitude automatically"
+           << " determined to be " << amplitude_
+           << " (max abs value of vext = " << magnitude_ << ")" << endl;
+  }
+
+  time = tm_read_vext.real();
+  tmin = time;
+  tmax = time;
+  ctxt->dmin(1,1,&tmin,1);
+  ctxt->dmax(1,1,&tmax,1);
+  if ( onpe0 )
+  {
+    cout << "  ExternalPotential::update: vext read time "
+         << "min: " << tmin << " max: " << tmax << endl;
+  }
+}
+
+////////////////////////////////////////////////////////////////////////////////
+double ExternalPotential::compute_eext(const ChargeDensity& cd)
+{
+  // Eext =  integral ( rhor * vext_r )
+  double eext = 0.0;
+  double omega = s_.wf.cell().volume();
+  const int np012loc = cd.vft()->np012loc();
+  const int np012 = cd.vft()->np012();
+  const std::vector<std::vector<double> > &rhor = cd.rhor;
+  for ( int ispin = 0; ispin < s_.wf.nspin(); ispin++ )
+  {
+    for ( int ir = 0; ir < np012loc; ir++ )
+    {
+      assert( rhor[ispin].size() == np012loc );
+      assert( vext_r_.size() == np012loc );
+      eext += rhor[ispin][ir] * v(ir);
+    }
+  }
+  double eext_sum = 0.0;
+  MPI_Allreduce(&eext,&eext_sum,1,MPI_DOUBLE,MPI_SUM,cd.vbasis()->comm());
+  eext_sum =  eext_sum * omega / np012;
+  return eext_sum;
+}

src/ExternalPotential.h

+////////////////////////////////////////////////////////////////////////////////
+//
+// Copyright (c) 2016 The Regents of the University of California
+//
+// This file is part of Qbox
+//
+// Qbox is distributed under the terms of the GNU General Public License
+// as published by the Free Software Foundation, either version 2 of
+// the License, or (at your option) any later version.
+// See the file COPYING in the root directory of this distribution
+// or <http://www.gnu.org/licenses/>.
+//
+////////////////////////////////////////////////////////////////////////////////
+//
+// ExternalPotential.h
+//
+////////////////////////////////////////////////////////////////////////////////
+#ifndef EXTERNALPOTENTIAL_H
+#define EXTERNALPOTENTIAL_H
+
+#include <vector>
+#include <string>
+#include "Sample.h"
+#include "ChargeDensity.h"
+
+class Sample;
+
+class ExternalPotential
+{
+  private:
+
+  Sample& s_;
+  int n_[3];               // real space grid size in 3 dimensions
+                           // read from cube file in cube file mode,
+                           // otherwise must be given in constructor
+  double ecut_;
+  double magnitude_;       // the magnitude of external potential, defined as
+                           // the average of its largest 0.1% (absolute) values
+  double amplitude_;       // overall scaling factor of external potential
+  vector<double> vext_r_;  // vext in real space
+  std::string filename_;   // file name for external potential
+  std::string fmt_;        // file format: "cube" or "xml"
+
+  public:
+
+  ExternalPotential(Sample& s, std::string name, std::string fmt="xml"):
+    s_(s), filename_(name), ecut_(0.0), amplitude_(1.0), magnitude_(0.0),
+    fmt_(fmt)
+  {
+    assert( fmt_ == "cube" || fmt_ == "xml" );
+  }
+  ~ExternalPotential() {}
+
+  int n(int i) const { return n_[i]; }
+  double ecut(void) const { return ecut_; }
+  double magnitude(void) const { return magnitude_; }
+  double amplitude(void) const { return amplitude_; }
+  std::string filename(void) const { return filename_; }
+  double v(size_t i) const { return amplitude_ * vext_r_[i]; }
+  void update(const ChargeDensity& cd);
+  void set_amplitude(double a) { amplitude_ = a; }
+  void reverse(void) {amplitude_ *= -1; }
+  double compute_eext(const ChargeDensity& cd);
+};
+#endif

src/Function3d.C

+//////////////////////////////////////////////////////////////////////////////// //
+// Copyright (c) 2018 The Regents of the University of California
+//
+// This file is part of Qbox
+//
+// Qbox is distributed under the terms of the GNU General Public License
+// as published by the Free Software Foundation, either version 2 of
+// the License, or (at your option) any later version.
+// See the file COPYING in the root directory of this distribution
+// or <http://www.gnu.org/licenses/>.
+//
+////////////////////////////////////////////////////////////////////////////////
+//
+// Function3d.C
+//
+////////////////////////////////////////////////////////////////////////////////
+
+#include "Function3d.h"
+#include "Function3dHandler.h"
+#include "Base64Transcoder.h"
+#include "qbox_xmlns.h"
+#include "Timer.h"
+#include <string>
+#include <iostream>
+#include <fstream>
+#include <sys/stat.h>
+#include <cstdio>
+using namespace std;
+#include <xercesc/util/XMLUniDefs.hpp>
+#include <xercesc/sax2/Attributes.hpp>
+#include <xercesc/util/PlatformUtils.hpp>
+#include <xercesc/sax2/SAX2XMLReader.hpp>
+#include <xercesc/sax2/XMLReaderFactory.hpp>
+#include <xercesc/framework/MemBufInputSource.hpp>
+using namespace xercesc;
+
+////////////////////////////////////////////////////////////////////////////////
+void Function3d::read(string filename)
+{
+  Timer tm;
+  struct stat statbuf;
+  bool filestat = !stat(filename.c_str(),&statbuf);
+  if ( !filestat )
+  {
+    cout << "Function3d::read: could not stat file " << filename << endl;
+    return;
+  }
+
+  FILE* infile;
+  infile = fopen(filename.c_str(),"r");
+  if ( !infile )
+  {
+    cout << "Function3d::read: could not open file " << filename