Commit 2ca7ee41 by He Ma

1. New functionality:

  1.1. Multiple options added to response command. Now the way to call response command is:
    response -vext vext_file [-RPA or -IPA] [-amplitude amplitude] [-parallel_write] nitscf [nite]
    a) -RPA option will freeze the update of Vxc (only work for semilocal functional now); IPA will freeze the update of Vxc and Vhartree (not implemented yet since it will involve major change of EnergyFunctional class). In order to implment RPA, some changes are made to EnergyFunctional class and XCPotential class. Additionally, two boolean flags are added in Control class to indicate whether Vxc or Vhartree is freezed.
    b) -amplitude option can set the amplitude of Vext. Right now it seems that this option is very convenient since it facilitates the test on whether the magnetitude for Vext is in linear-response regime. So I suggest keeping this option.    c) -parallel_write: not implemented yet. By this option I intend to realize parallel I/O, like in SampleReader or SlaterDet::write. We may need to change the format of file into binary (e.g. base64) to implement this.


2. Improvement of implementations:
  2.1. Now response command will use the GS charge density and charge density under +Vext to generate an initial guess (by 2*rho0 - rho(+Vext)) for the calculation under -Vext. By doing this improvement it is observed that for silane the convergence of SCF iterations under -Vext becomes 3 times faster. In order to implement this there are some minor changes made to BOSamplerStepper.
  2.2. Now in ExternalPotential::update process 0 will first read the file and MPI_Scatterv to other processes in column 0; then processes in column 0 will perform in-row MPI_Bcast to send Vext to all other processes. Likewise, now in ResponseCmd::responseVext process 0 will MPI_Gatherv drho_r from all processes on column 0. These functionalities are encapsulated in two functions in ExternalPotential.C and ResponseCmd.C respectively, to make the code more modular. As mentioned before, the IO performance might be further improved by resorting to binary format and MPI file operations.

3. Other comments:
  3.1.  Regarding whether Vext should be updated when "set vext" command is called: right now our construction of Vext strongly depend on the ChargeDensity instance passed to update() function, so it seems that we can only update vext during "run" or "response" command. It would be interesting to consider if we can have a ChargeDensity member in Sample class, as you mentioned earlier. But it involves major changes. 
  As a compromise, would it be possible if we keep a ChargeDensity member in Sample class (and we still have another ChargeDensity instance created in a BO sampler) and only update it manually by explicit command, or at the destruction of BO sampler? In this way, the awkwardness that Vext can only be updated inside a BO sampler where a charge density exists can be removed.
  3.2. In ResponseCmd::responseVext function, a previous modification that changes s->wf.spincontext() to cd.context() seems to cause some tricky problem... Right now I rearlly don't understand why this will happen, but it appears that things work fine after I reverted this change. There might be some hidden bugs here that may be significant when we test against large calculations where the context is complicated. I will return to this problem later.
  3.3. I understand that the option list of response command is too long compared to other commands in Qbox and it may be unpleasant. Later we may simplify it.


git-svn-id: http://qboxcode.org/svn/qb/branches/vext@1884 cba15fb0-1239-40c8-b417-11db7ca47a34
parent 5d982e59
......@@ -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_atomic_density(false), initial_density(false) {}
////////////////////////////////////////////////////////////////////////////////
BOSampleStepper::~BOSampleStepper()
......@@ -159,6 +159,20 @@ void BOSampleStepper::initialize_density(void)
}
////////////////////////////////////////////////////////////////////////////////
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;
}
////////////////////////////////////////////////////////////////////////////////
void BOSampleStepper::step(int niter)
{
const Context& ctxt = s_.ctxt_;
......@@ -752,6 +766,8 @@ void BOSampleStepper::step(int niter)
tmap["charge"].start();
if ( itscf==0 && initial_atomic_density )
cd_.update_rhor();
else if ( itscf==0 && initial_density )
cd_.update_rhog();
else
cd_.update_density();
tmap["charge"].stop();
......@@ -836,7 +852,20 @@ void BOSampleStepper::step(int niter)
} // if nite_ > 0
tmap["update_vhxc"].start();
ef_.update_vhxc(compute_stress);
bool freeze_vh;
bool freeze_vxc;
if ( itscf == 0 )
{
// at first SCF iteration, vhxc must be updated
freeze_vh = false;
freeze_vxc = false;
}
else
{
freeze_vh = s_.ctrl.freeze_vh;
freeze_vxc = s_.ctrl.freeze_vxc;
}
ef_.update_vhxc(compute_stress, freeze_vh, freeze_vxc);
tmap["update_vhxc"].stop();
// reset stepper only if multiple non-selfconsistent steps
......@@ -1299,4 +1328,5 @@ void BOSampleStepper::step(int niter)
if ( ntc_extrapolation || asp_extrapolation ) delete wfmm;
initial_atomic_density = false;
initial_density = false;
}
......@@ -43,6 +43,7 @@ class BOSampleStepper : public SampleStepper
IonicStepper* ionic_stepper;
bool initial_atomic_density;
bool initial_density;
// Do not allow construction of BOSampleStepper unrelated to a Sample
BOSampleStepper(void);
......@@ -52,12 +53,15 @@ class BOSampleStepper : public SampleStepper
mutable TimerMap tmap;
void step(int niter);
void initialize_density(void);
void initialize_density(void); // initialize density by atomic density
void initialize_density(const vector<vector<double> >& rhor); // initialize density by given density
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
......@@ -214,3 +214,34 @@ void ChargeDensity::update_rhor(void)
}
}
}
////////////////////////////////////////////////////////////////////////////////
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]);
}
}
......@@ -55,6 +55,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_; }
......
......@@ -75,5 +75,8 @@ 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
......@@ -224,7 +224,8 @@ EnergyFunctional::~EnergyFunctional(void)
}
////////////////////////////////////////////////////////////////////////////////
void EnergyFunctional::update_vhxc(bool compute_stress)
void EnergyFunctional::update_vhxc(bool compute_stress,
bool freeze_vh, bool freeze_vxc)
{
// called when the charge density has changed
// update Hartree and xc potentials using the charge density cd_
......@@ -261,7 +262,8 @@ void EnergyFunctional::update_vhxc(bool compute_stress)
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);
xco->update(v_r, compute_stress, freeze_vxc);
exc_ = xco->exc();
dxc_ = xco->dxc();
if ( compute_stress )
......@@ -344,6 +346,8 @@ 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);
}
......
......@@ -109,7 +109,8 @@ class EnergyFunctional
const ConfinementPotential *confpot(int ikp) const { return cfp[ikp]; }
void update_vhxc(bool compute_stress);
void update_vhxc(bool compute_stress,
bool freeze_vh=false, bool freeze_vxc=false);
void atoms_moved(void);
void cell_moved(void);
......
////////////////////////////////////////////////////////////////////////////////
//
// Copyright (c) 2016 The Regents of the University of California
......@@ -25,97 +26,87 @@ 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)
{
// read vext on nrow=0 tasks
vector<double> vext_read;
FourierTransform *vft = cd.vft();
const Context& ctxt = cd.context();
const int myrow = ctxt.myrow();
if ( myrow == 0 )
// process 0 reads the grid size and bcast to all processes
if ( ctxt.onpe0() )
{
// read the external potential from file s_.ctrl.vext
ifstream vfile(filename_.c_str());
if ( ctxt.onpe0() )
if ( vfile )
cout << " ExternalPotential::update: read external "
<< "potential from file: " << filename_ << endl;
else
{
if ( vfile )
cout << "ExternalPotential::update: read external "
<< "potential from file: " << filename_ << endl;
else
{
cout << "ExternalPotential::update: file not found: "
<< filename_ << endl;
ctxt.abort(1);
}
cout << " ExternalPotential::update: file not found: "
<< filename_ << endl;
ctxt.abort(1);
}
vfile >> n_[0] >> n_[1] >> n_[2];
if ( cd.context().onpe0() )
cout << "ExternalPotential::update: grid size "
<< n_[0] << " " << n_[1] << " " << n_[2] << endl;
int n012 = n_[0] * n_[1] * n_[2];
vext_read.resize(n012);
for ( int i = 0; i < n012; i++ )
vfile >> vext_read[i];
cout << " ExternalPotential::update: grid size "
<< n_[0] << " " << n_[1] << " " << n_[2] << endl;
vfile.close();
}
// broadcast sizes to lower rows
MPI_Bcast(&n_[0],3,MPI_INT,0,cd.vcomm());
int n012 = n_[0] * n_[1] * n_[2];
vext_read.resize(n012);
MPI_Bcast(&vext_read[0],n012,MPI_DOUBLE,0,cd.vcomm());
// vext_read now contains vext data on all tasks
// todo: replace Bcast with Scatter
// create a Basis compatible with the vext grid read from file
// determine largest ecut compatible with grid and current unit cell
const Basis* vbasis = cd.vbasis();
const UnitCell& cell = vbasis->cell();
const D3vector b0 = cell.b(0);
const D3vector b1 = cell.b(1);
const D3vector b2 = cell.b(2);
// determine largest ecut compatible with grid and current unit cell
int n0max = (n_[0]-2)/2;
int n1max = (n_[1]-2)/2;
int n2max = (n_[2]-2)/2;
double ecut0 = 0.5 * norm2(b0) * n0max*n0max;
double ecut1 = 0.5 * norm2(b1) * n1max*n1max;
double ecut2 = 0.5 * norm2(b2) * n2max*n2max;
ecut_ = min(min(ecut0,ecut1),ecut2);
if ( cd.context().onpe0() )
if ( ctxt.onpe0() )
{
cout << "ExternalPotential::update: ecut0: " << ecut0 << endl;
cout << "ExternalPotential::update: ecut1: " << ecut1 << endl;
cout << "ExternalPotential::update: ecut2: " << ecut2 << endl;
cout << "ExternalPotential::update: ecut: " << ecut_ << endl;
cout << " ExternalPotential::update: ecut0: " << ecut0 << endl;
cout << " ExternalPotential::update: ecut1: " << ecut1 << endl;
cout << " ExternalPotential::update: ecut2: " << ecut2 << endl;
cout << " ExternalPotential::update: ecut: " << ecut_ << endl;
}
Basis basis(cd.vcomm(),D3vector(0,0,0));
basis.resize(cell,cell,ecut_);
if ( cd.context().onpe0() )
if ( 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;
cout << " ExternalPotential::update: np0: " << basis.np(0) << endl;
cout << " ExternalPotential::update: np1: " << basis.np(1) << endl;
cout << " ExternalPotential::update: np2: " << basis.np(2) << endl;
}
// interpolate on grid compatible with the charge density cd
// build FT to transform vext to g space
FourierTransform ft(basis,n_[0],n_[1],n_[2]);
// read vext from file
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);
vector<complex<double> > tmp_r(ft.np012loc());
// index of local vext slice in global vext array
int istart = ft.np2_first() * n_[0] * n_[1];
for ( int i = 0; i < tmp_r.size(); i++ )
tmp_r[i] = vext_read[istart+i];
for ( int ir = 0; ir < tmp_r.size(); ir++ )
tmp_r[ir] = complex<double>(vext_r_read[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
// interpolate to vft grid of charge density
FourierTransform *vft = cd.vft();
FourierTransform ft2(basis,vft->np0(),vft->np1(),vft->np2());
tmp_r.resize(vft->np012loc());
vext_r_.resize(tmp_r.size());
......@@ -123,3 +114,99 @@ 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())
{
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++ )
{
displ += ft.np012loc(iproc);
displs[iproc] = displ;
scounts[iproc] = ft.np012loc(iproc);
}
// process 0 scatter fr to other processes on first column
MPI_Barrier(comm);
if ( ctxt.mycol() == 0 )
{
//MPI_Barrier(comm_col);
MPI_Scatterv(sbuf,&scounts[0],&displs[0],MPI_DOUBLE,
rbuf,ft.np012loc(),MPI_DOUBLE,0,comm_col);
//MPI_Barrier(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
}
......@@ -49,5 +49,6 @@ class ExternalPotential
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; }
};
#endif
......@@ -27,14 +27,18 @@ class ResponseCmd : public Cmd
{
private:
void responseVext(int nitscf, int nite);
void responseVext(int nitscf, int nite, bool parallel_output);
void responseEfield(double amplitude, int nitscf, int nite);
public:
Sample *s;
ResponseCmd(Sample *sample) : s(sample) {};
ResponseCmd(Sample *sample) : s(sample) {
s->ctrl.freeze_vh = false;
s->ctrl.freeze_vxc = false;
};
const char *name(void) const { return "response"; }
const char *help_msg(void) const
......@@ -42,7 +46,7 @@ class ResponseCmd : public Cmd
return
"\n response\n\n"
" syntax: response amplitude nitscf [nite]\n"
" response -vext vext_file nitscf [nite]\n\n"
" response -vext vext_file [-RPA or -IPA] [-amplitude amplitude] [-parallel_write] 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"
......
......@@ -49,12 +49,16 @@ class Vext : public Var
// set vext NULL
delete s->vext;
s->vext = 0;
if ( ui->onpe0() )
cout << " vext set to zero" << endl;
}
else
{
if ( s->vext )
delete s->vext;
s->vext = new ExternalPotential(*s,argv[1]);
if ( ui->onpe0() )
cout << " vext set to" << argv[1] << endl;
}
return 0;
......
......@@ -94,7 +94,8 @@ XCOperator::~XCOperator()
}
////////////////////////////////////////////////////////////////////////////////
void XCOperator::update(std::vector<std::vector<double> >& vr, bool compute_stress)
void XCOperator::update(std::vector<std::vector<double> >& vr,
bool compute_stress, bool freeze_vxc)
{
// update xc potential and self-energy
// used whenever the charge density and/or wave functions have changed
......@@ -102,7 +103,7 @@ void XCOperator::update(std::vector<std::vector<double> >& vr, bool compute_stre
if ( hasPotential_ )
{
// update LDA/GGA xc potential
xcp_->update( vr );
xcp_->update(vr, freeze_vxc);
// LDA/GGA exchange energy
exc_ = xcp_->exc();
......
......@@ -60,7 +60,8 @@ class XCOperator
bool hasGGA(void) { return hasGGA_; };
bool hasHF(void) { return hasHF_; };
void update(std::vector<std::vector<double> >& vr, bool compute_stress);
void update(std::vector<std::vector<double> >& vr,
bool compute_stress, bool freeze_vxc);
void apply_self_energy(Wavefunction &dwf);
void compute_stress(std::valarray<double>& sigma);
void cell_moved(void);
......
......@@ -91,7 +91,7 @@ bool XCPotential::isGGA(void)
}
////////////////////////////////////////////////////////////////////////////////
void XCPotential::update(vector<vector<double> >& vr)
void XCPotential::update(vector<vector<double> >& vr, bool freeze_vxc)
{
// compute exchange-correlation energy and add vxc potential to vr[ispin][ir]
......@@ -122,7 +122,17 @@ void XCPotential::update(vector<vector<double> >& vr)
{
// LDA functional
xcf_->setxc();
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;
}
exc_ = 0.0;
dxc_ = 0.0;
......@@ -170,6 +180,8 @@ void XCPotential::update(vector<vector<double> >& vr)
else
{
// GGA functional
if ( ! freeze_vxc )
{
exc_ = 0.0;
// compute grad_rho
......@@ -217,6 +229,14 @@ void XCPotential::update(vector<vector<double> >& vr)
}
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
......
......@@ -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);
void update(std::vector<std::vector<double> >& vr, bool freeze_vxc);
void compute_stress(std::valarray<double>& sigma_exc);
double exc(void) { return exc_; }
double dxc(void) { return dxc_; }
......
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