Commit edf03449 by mahe

Added IPA to response command. Made slight modification to vext variable to…

Added IPA to response command. Made slight modification to vext variable to facilitate cDFT implementaion.
parent e39dbb70
...@@ -46,7 +46,7 @@ using namespace std; ...@@ -46,7 +46,7 @@ using namespace std;
BOSampleStepper::BOSampleStepper(Sample& s, int nitscf, int nite) : BOSampleStepper::BOSampleStepper(Sample& s, int nitscf, int nite) :
SampleStepper(s), cd_(s.wf), ef_(s,cd_), SampleStepper(s), cd_(s.wf), ef_(s,cd_),
dwf(s.wf), wfv(s.wfv), nitscf_(nitscf), nite_(nite), dwf(s.wf), wfv(s.wfv), nitscf_(nitscf), nite_(nite),
update_density_first_(true), update_vxc_(true) {} update_density_first_(true), update_vh_(true), update_vxc_(true) {}
//////////////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////////////
BOSampleStepper::~BOSampleStepper() BOSampleStepper::~BOSampleStepper()
...@@ -848,20 +848,12 @@ void BOSampleStepper::step(int niter) ...@@ -848,20 +848,12 @@ void BOSampleStepper::step(int niter)
// at first scf step: // at first scf step:
// - update both vh and vxc // - update both vh and vxc
// at later steps: // at later steps:
// - update vh // - update depending of values of update_vh_ and update_vxc_
// - update vxc only if update_vxc_ is true
tmap["update_vhxc"].start(); tmap["update_vhxc"].start();
if ( itscf == 0 ) if ( itscf == 0 )
ef_.update_vhxc(compute_stress); ef_.update_vhxc(compute_stress);
else else
{ ef_.update_vhxc(compute_stress, update_vh_, update_vxc_);
if ( update_vxc_ )
ef_.update_vhxc(compute_stress);
else
// update vh only
ef_.update_vh(compute_stress);
}
tmap["update_vhxc"].stop(); tmap["update_vhxc"].stop();
// reset stepper only if multiple non-selfconsistent steps // reset stepper only if multiple non-selfconsistent steps
......
...@@ -43,7 +43,7 @@ class BOSampleStepper : public SampleStepper ...@@ -43,7 +43,7 @@ class BOSampleStepper : public SampleStepper
IonicStepper* ionic_stepper; IonicStepper* ionic_stepper;
bool update_density_first_; bool update_density_first_;
bool update_vxc_; bool update_vh_, update_vxc_;
// Do not allow construction of BOSampleStepper unrelated to a Sample // Do not allow construction of BOSampleStepper unrelated to a Sample
BOSampleStepper(void); BOSampleStepper(void);
...@@ -55,6 +55,7 @@ class BOSampleStepper : public SampleStepper ...@@ -55,6 +55,7 @@ class BOSampleStepper : public SampleStepper
void step(int niter); void step(int niter);
// initialize density with sum of atomic densities // initialize density with sum of atomic densities
void initialize_density(void); 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_vxc(bool update_vxc) { update_vxc_ = update_vxc; }
void set_update_density_first(bool update_density_first) void set_update_density_first(bool update_density_first)
{ update_density_first_ = update_density_first; } { update_density_first_ = update_density_first; }
......
...@@ -25,6 +25,7 @@ ...@@ -25,6 +25,7 @@
#include <iomanip> #include <iomanip>
#include <algorithm> // fill #include <algorithm> // fill
#include <functional> #include <functional>
#include <fstream>
using namespace std; using namespace std;
//////////////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////////////
...@@ -267,4 +268,4 @@ void ChargeDensity::update_rhog(void) ...@@ -267,4 +268,4 @@ void ChargeDensity::update_rhog(void)
vft_->forward(&rhotmp[0],&rhog[ispin][0]); vft_->forward(&rhotmp[0],&rhog[ispin][0]);
} }
} }
\ No newline at end of file
...@@ -26,6 +26,7 @@ ...@@ -26,6 +26,7 @@
#include <map> #include <map>
#include "Timer.h" #include "Timer.h"
#include "Context.h" #include "Context.h"
#include "D3vector.h"
class Wavefunction; class Wavefunction;
class FourierTransform; class FourierTransform;
......
...@@ -226,7 +226,7 @@ EnergyFunctional::~EnergyFunctional(void) ...@@ -226,7 +226,7 @@ EnergyFunctional::~EnergyFunctional(void)
} }
//////////////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////////////
void EnergyFunctional::update_vhxc(bool compute_stress, bool update_vxc) void EnergyFunctional::update_vhxc(bool compute_stress, bool update_vh, bool update_vxc)
{ {
// called when the charge density has changed // called when the charge density has changed
// update Hartree and xc potentials using the charge density cd_ // update Hartree and xc potentials using the charge density cd_
...@@ -315,7 +315,7 @@ void EnergyFunctional::update_vhxc(bool compute_stress, bool update_vxc) ...@@ -315,7 +315,7 @@ void EnergyFunctional::update_vhxc(bool compute_stress, bool update_vxc)
ehesum += norm(r) * g2i[ig]; ehesum += norm(r) * g2i[ig];
ehepsum += 2.0*real(conj(r)*rp * g2i[ig]); ehepsum += 2.0*real(conj(r)*rp * g2i[ig]);
ehpsum += norm(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 // factor 1/2 from definition of Ehart cancels with half sum over G
// Note: rhogt[ig] includes a factor 1/Omega // Note: rhogt[ig] includes a factor 1/Omega
......
...@@ -111,11 +111,12 @@ class EnergyFunctional ...@@ -111,11 +111,12 @@ class EnergyFunctional
const ConfinementPotential *confpot(int ikp) const { return cfp[ikp]; } const ConfinementPotential *confpot(int ikp) const { return cfp[ikp]; }
void update_vhxc(bool compute_stress, bool update_vxc); void update_vhxc(bool compute_stress, bool update_vh, bool update_vxc);
// update both vh and vxc // update both vh and vxc
void update_vhxc(bool compute_stress) { update_vhxc(compute_stress, true); } void update_vhxc(bool compute_stress)
// update only vh {
void update_vh(bool compute_stress) { update_vhxc(compute_stress, false); } update_vhxc(compute_stress, true, true);
}
void atoms_moved(void); void atoms_moved(void);
void cell_moved(void); void cell_moved(void);
......
...@@ -46,7 +46,7 @@ class ExternalPotential ...@@ -46,7 +46,7 @@ class ExternalPotential
ExternalPotential(Sample& s, std::string name, std::string io="cube", ExternalPotential(Sample& s, std::string name, std::string io="cube",
int nx=0, int ny=0, int nz=0): int nx=0, int ny=0, int nz=0):
s_(s), filename_(name), ecut_(0.0), amplitude_(0.0), magnitude_(0.0), io_(io){ s_(s), filename_(name), ecut_(0.0), amplitude_(1.0), magnitude_(0.0), io_(io){
assert( io_ == "cube" || io == "base64_serial" || io == "base64_parallel" ); assert( io_ == "cube" || io == "base64_serial" || io == "base64_parallel" );
if (io != "cube") if (io != "cube")
{ {
......
...@@ -38,8 +38,8 @@ using namespace std; ...@@ -38,8 +38,8 @@ using namespace std;
int ResponseCmd::action(int argc, char **argv) int ResponseCmd::action(int argc, char **argv)
{ {
// " syntax: response amplitude nitscf [nite]\n\n" // " syntax: response amplitude nitscf [nite]\n\n"
// " syntax: response -vext vext_file [-RPA] [-amplitude a] // " syntax: response -vext vext_file [-RPA|-IPA] [-amplitude a]
// [-io iomode -nx nx -ny ny -nz nz] nitscf [nite]\n\n" // [-io iomode -nx nx -ny ny -nz nz] [-q qx qy qz] nitscf [nite]\n\n"
if ( s->wf.nst() == 0 ) if ( s->wf.nst() == 0 )
{ {
...@@ -60,23 +60,42 @@ int ResponseCmd::action(int argc, char **argv) ...@@ -60,23 +60,42 @@ int ResponseCmd::action(int argc, char **argv)
int iarg = 1; int iarg = 1;
if ( !strcmp(argv[iarg],"-vext") ) if ( !strcmp(argv[iarg],"-vext") )
{ {
// response to vext
if ( s->vext )
{
if ( ui->onpe0() )
cout << " ResponseCmd: cannot run when vext is already set" << endl;
return 1;
}
string filename; string filename;
bool rpa = false; bool rpa = false;
double amplitude = atof(argv[iarg]); bool ipa = false;
double amplitude = 0.0;
string io = "cube"; string io = "cube";
int nx, ny, nz; int nx, ny, nz;
nx = ny = nz = 0; nx = ny = nz = 0;
iarg++; iarg++;
// response to vext
if ( s->vext )
delete s->vext;
filename = argv[iarg]; filename = argv[iarg];
iarg++; iarg++;
if ( !strcmp(argv[iarg],"-RPA") ) if ( !strcmp(argv[iarg],"-RPA") )
{ {
rpa = true; rpa = true;
iarg++; iarg++;
} }
if ( !strcmp(argv[iarg],"-IPA") )
{
ipa = true;
iarg++;
}
if ( rpa && ipa )
{
if ( ui->onpe0() )
cout << " Only one of -RPA or -IPA can be specified" << endl;
return 1;
}
if ( !strcmp(argv[iarg],"-amplitude") ) if ( !strcmp(argv[iarg],"-amplitude") )
{ {
iarg++; iarg++;
...@@ -106,6 +125,7 @@ int ResponseCmd::action(int argc, char **argv) ...@@ -106,6 +125,7 @@ int ResponseCmd::action(int argc, char **argv)
iarg++; iarg++;
} }
} }
nitscf = atoi(argv[iarg]); nitscf = atoi(argv[iarg]);
iarg++; iarg++;
...@@ -121,7 +141,9 @@ int ResponseCmd::action(int argc, char **argv) ...@@ -121,7 +141,9 @@ int ResponseCmd::action(int argc, char **argv)
s->vext = new ExternalPotential(*s, filename, io, nx, ny, nz); s->vext = new ExternalPotential(*s, filename, io, nx, ny, nz);
s->vext->set_amplitude(amplitude); s->vext->set_amplitude(amplitude);
responseVext(rpa, nitscf, nite, io, nx, ny, nz); responseVext(rpa, ipa, nitscf, nite, io);
delete s->vext;
s->vext = 0;
} }
else else
{ {
...@@ -222,9 +244,15 @@ void ResponseCmd::responseEfield(double amplitude, int nitscf, int nite) ...@@ -222,9 +244,15 @@ void ResponseCmd::responseEfield(double amplitude, int nitscf, int nite)
} }
//////////////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////////////
void ResponseCmd::responseVext(bool rpa, int nitscf, int nite, void ResponseCmd::responseVext(bool rpa, bool ipa, int nitscf, int nite, string io)
string io, int nxx, int nyy, int nzz)
{ {
// if (ui->onpe0())
// {
// cout << "ResponseCmd:responseVext:\n"
// << "RPA = " << rpa << ", IPA = " << ipa << ", io = " << io
// << ", q = " << q << ", nitscf = " << nitscf << ", nite = " << nite << "\n";
// }
s->wf.info(cout, "wavefunction"); s->wf.info(cout, "wavefunction");
const int nspin = s->wf.nspin(); const int nspin = s->wf.nspin();
...@@ -233,6 +261,12 @@ void ResponseCmd::responseVext(bool rpa, int nitscf, int nite, ...@@ -233,6 +261,12 @@ void ResponseCmd::responseVext(bool rpa, int nitscf, int nite,
if (rpa) if (rpa)
stepper->set_update_vxc(false); stepper->set_update_vxc(false);
if (ipa)
{
stepper->set_update_vh(false);
stepper->set_update_vxc(false);
}
// save a copy of initial wave functions // save a copy of initial wave functions
Wavefunction wf0(s->wf); Wavefunction wf0(s->wf);
wf0 = s->wf; wf0 = s->wf;
...@@ -241,13 +275,15 @@ void ResponseCmd::responseVext(bool rpa, int nitscf, int nite, ...@@ -241,13 +275,15 @@ void ResponseCmd::responseVext(bool rpa, int nitscf, int nite,
MPI_Comm vcomm = cd.vcomm(); MPI_Comm vcomm = cd.vcomm();
stepper->step(0); stepper->step(0);
const vector<vector<double> > rhor1 = cd.rhor; // density with +Vext vector<vector<double> > rhor1; // density with +Vext
rhor1 = cd.rhor;
s->wf = wf0; s->wf = wf0;
s->vext->reverse(); s->vext->reverse();
stepper->step(0); stepper->step(0);
const vector<vector<double> > rhor2 = cd.rhor; // density with -Vext vector<vector<double> > rhor2; // density with -Vext
rhor2 = cd.rhor;
s->vext->reverse(); s->vext->reverse();
// restore initial wave functions // restore initial wave functions
......
...@@ -27,9 +27,8 @@ class ResponseCmd : public Cmd ...@@ -27,9 +27,8 @@ class ResponseCmd : public Cmd
{ {
private: private:
void responseVext(bool rpa, int nitscf, int nite,
string io, int nx, int ny, int nz);
void responseEfield(double amplitude, int nitscf, int nite); void responseEfield(double amplitude, int nitscf, int nite);
void responseVext(bool rpa, bool ipa, int nitscf, int nite, string io);
public: public:
...@@ -43,22 +42,25 @@ class ResponseCmd : public Cmd ...@@ -43,22 +42,25 @@ class ResponseCmd : public Cmd
return return
"\n response\n\n" "\n response\n\n"
" syntax: response amplitude nitscf [nite]\n" " syntax: response amplitude nitscf [nite]\n"
" response -vext vext_file [-RPA] [-amplitude a] \n" " response -vext vext_file [-RPA|-IPA] [-amplitude a] \n"
" [-io iomode -nx nx -ny ny -nz nz] nitscf [nite]\n\n" " [-io iomode -nx nx -ny ny -nz nz] nitscf [nite]\n\n"
" The response command computes the polarizability tensor by\n" " The response command computes the polarizability tensor by\n"
" finite differences using external electric fields in the x,y,z\n" " finite differences using external electric fields in the x,y,z\n"
" directions with magnitude defined by the amplitude argument.\n" " directions with magnitude defined by the amplitude argument.\n"
" If the -vext option is used, the response command computes the\n" " If the -vext option is used, the response command computes the\n"
" response to the external potential defined in the file vext_file.\n" " response to the external potential defined in the file vext_file.\n"
" Acceptable control flags are:\n" " Control flags:\n"
" 1. -RPA. Compute response within Random Phase Approximation, Vxc is frozen.\n" " -RPA Compute response within the Random Phase Approximation,\n"
" 2. -amplitude a. Scale the Vext by a before any calculations, \n" " Vxc is frozen.\n"
" then scale the charge density response by 1/a before output.\n" " -IPA Compute response within the Independent Particle Approximation,\n"
" 3. -io iomode. How the vext/response file shall be read/write. Possible choice:\n" " VHartree and Vxc are frozen.\n"
" cube: Gaussian cube format.\n" " -amplitude a\n"
" base64_serial or base64_parallel: base64-encoded binary grid function.\n" " Scale the external potential by a before any calculation, \n"
" grid size need to be specified by nx, ny, nz.\n" " then scale the charge density response by 1/a before output.\n"
" use cube for best compatibility and base64_parallel for best performance\n\n"; " -io iomode\n"
" Valid choices of iomode: cube, base64_serial, base64_parallel\n"
" -nx nx, -ny ny, -nz nz\n"
" grid size (for base64_serial and base64_parallel only)\n\n";
} }
int action(int argc, char **argv); int action(int argc, char **argv);
......
...@@ -37,26 +37,37 @@ class Vext : public Var ...@@ -37,26 +37,37 @@ class Vext : public Var
int set ( int argc, char **argv ) int set ( int argc, char **argv )
{ {
if ( argc > 3 ) if ( argc == 2 )
{ {
if ( ui->onpe0() ) if ( !strcmp(argv[1],"NULL") )
cout << " vext takes only one value" << endl; {
return 1; // set vext NULL
delete s->vext;
s->vext = 0;
}
else
{
if ( s->vext )
delete s->vext;
s->vext = new ExternalPotential(*s,argv[1]);
}
} }
else if ( argc == 3 and !strcmp(argv[1],"amplitude") )
if ( !strcmp(argv[1],"NULL") )
{ {
// set vext NULL if ( s->vext )
delete s->vext; s->vext->set_amplitude(atof(argv[2]));
s->vext = 0; else
{
cout << "vext not set yet" << endl;
return 1;
}
} }
else else
{ {
if ( s->vext ) if ( ui->onpe0() )
delete s->vext; cout << "unknown option for vext variable" << endl;
s->vext = new ExternalPotential(*s,argv[1]); return 1;
} }
return 0; return 0;
} }
...@@ -68,7 +79,8 @@ class Vext : public Var ...@@ -68,7 +79,8 @@ class Vext : public Var
if ( s->vext ) if ( s->vext )
{ {
st.setf(ios::right,ios::adjustfield); st.setf(ios::right,ios::adjustfield);
st << setw(10) << s->vext->filename(); st << setw(10) << s->vext->filename()
<< " (amplitude" << " = " << s->vext->amplitude() << ")";
} }
return st.str(); return st.str();
} }
......
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