Commit beb58eb3 by Francois Gygi

Modified vext implementation.


git-svn-id: http://qboxcode.org/svn/qb/branches/vext@1879 cba15fb0-1239-40c8-b417-11db7ca47a34
parent b4b20231
......@@ -75,7 +75,5 @@ struct Control
std::string iter_cmd;
int iter_cmd_period;
//std::string vext; // filename for external potential
};
#endif
......@@ -81,7 +81,7 @@ EnergyFunctional::EnergyFunctional( Sample& s, ChargeDensity& cd)
}
// external potential
if ( ! s_.vext->filename.empty() )
if ( s_.vext )
s_.vext->update(cd_);
const int ngloc = vbasis_->localsize();
......@@ -342,7 +342,7 @@ void EnergyFunctional::update_vhxc(bool compute_stress)
vft->backward(&vlocal_g[0],&tmp_r[0]);
// add external potential vext to tmp_r
if ( ! s_.vext->filename.empty() )
if ( s_.vext )
{
for ( int i = 0; i < tmp_r.size(); i++ )
tmp_r[i] += s_.vext->v(i);
......
......@@ -36,16 +36,16 @@ void ExternalPotential::update(const ChargeDensity& cd)
if ( myrow == 0 )
{
// read the external potential from file s_.ctrl.vext
ifstream vfile(filename.c_str());
if ( cd.context().onpe0() )
ifstream vfile(filename_.c_str());
if ( ctxt.onpe0() )
{
if ( vfile )
cout << "ExternalPotential::update: read external "
<< "potential from file: " << filename << endl;
<< "potential from file: " << filename_ << endl;
else
{
cout << "ExternalPotential::update: file not found: "
<< filename << endl;
<< filename_ << endl;
ctxt.abort(1);
}
}
......@@ -104,7 +104,6 @@ void ExternalPotential::update(const ChargeDensity& cd)
// interpolate on grid compatible with the charge density cd
FourierTransform ft(basis,n_[0],n_[1],n_[2]);
vext_g_.resize(basis.localsize());
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];
......@@ -112,6 +111,7 @@ void ExternalPotential::update(const ChargeDensity& cd)
tmp_r[i] = vext_read[istart+i];
// 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
......@@ -123,16 +123,3 @@ void ExternalPotential::update(const ChargeDensity& cd)
for ( int i = 0; i < vext_r_.size(); i++ )
vext_r_[i] = real(tmp_r[i]);
}
void ExternalPotential::reverse()
{
for ( int ir = 0; ir < vext_r_.size(); ir++ )
vext_r_[ir] *= -1;
for ( int ig = 0; ig < vext_g_.size(); ig++ )
vext_g_[ig] *= -1;
//transform(vext_r_.begin(), vext_r_.end(), vext_r_.begin(),
// bind1st(multiplies<double>(), -1));
//transform(vext_g_.begin(), vext_g_.end(), vext_g_.begin(),
// bind1st(multiplies<complex<double> >(), -1));
}
......@@ -20,8 +20,6 @@
#include <vector>
#include <string>
using namespace std;
#include "Sample.h"
#include "ChargeDensity.h"
......@@ -32,22 +30,24 @@ class ExternalPotential
private:
Sample& s_;
int n_[3]; // real space grid size in 3 dimensions
int n_[3]; // real space grid size in 3 dimensions
double ecut_;
double amplitude_;
vector<double> vext_r_; // vext in real space
vector<complex<double> > vext_g_; // vext in g space
std::string filename_; // filename for external potential
public:
ExternalPotential(Sample& s): s_(s), ecut_(0.0) {}
ExternalPotential(Sample& s,std::string name): s_(s),
filename_(name), ecut_(0.0), amplitude_(1.0) {}
~ExternalPotential() {}
string filename; // filename for external potential
int n(int i) const { return n_[i]; }
double ecut(void) const { return ecut_; }
double v(size_t i) const { return vext_r_[i]; }
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 reverse();
void set_amplitude(double a) { amplitude_ = a; }
};
#endif
......@@ -381,10 +381,10 @@ ExtForceSet.o: blacs.h UnitCell.h D3tensor.h blas.h
ExtStress.o: Sample.h AtomSet.h Context.h blacs.h Atom.h D3vector.h
ExtStress.o: UnitCell.h D3tensor.h blas.h ConstraintSet.h ExtForceSet.h
ExtStress.o: Wavefunction.h Control.h
ExternalPotential.o: ExternalPotential.h Sample.h AtomSet.h Context.h blacs.h
ExternalPotential.o: Atom.h D3vector.h UnitCell.h D3tensor.h blas.h
ExternalPotential.o: ConstraintSet.h ExtForceSet.h Wavefunction.h Control.h
ExternalPotential.o: ChargeDensity.h Timer.h FourierTransform.h Basis.h
ExternalPotential.o: Basis.h D3vector.h UnitCell.h ExternalPotential.h
ExternalPotential.o: Sample.h AtomSet.h Context.h blacs.h Atom.h D3tensor.h
ExternalPotential.o: blas.h ConstraintSet.h ExtForceSet.h Wavefunction.h
ExternalPotential.o: Control.h ChargeDensity.h Timer.h FourierTransform.h
ExternalPotential.o: Sample.h AtomSet.h Context.h blacs.h Atom.h D3vector.h
ExternalPotential.o: UnitCell.h D3tensor.h blas.h ConstraintSet.h
ExternalPotential.o: ExtForceSet.h Wavefunction.h Control.h ChargeDensity.h
......@@ -547,12 +547,13 @@ RescaleVCmd.o: ExtForceSet.h Wavefunction.h Control.h
ResetVcmCmd.o: UserInterface.h Sample.h AtomSet.h Context.h blacs.h Atom.h
ResetVcmCmd.o: D3vector.h UnitCell.h D3tensor.h blas.h ConstraintSet.h
ResetVcmCmd.o: ExtForceSet.h Wavefunction.h Control.h
ResponseCmd.o: ResponseCmd.h UserInterface.h BOSampleStepper.h
ResponseCmd.o: SampleStepper.h Timer.h EnergyFunctional.h StructureFactor.h
ResponseCmd.o: ElectricEnthalpy.h Matrix.h Context.h blacs.h D3vector.h
ResponseCmd.o: Wavefunction.h UnitCell.h SlaterDet.h Basis.h Sample.h
ResponseCmd.o: AtomSet.h Atom.h D3tensor.h blas.h ConstraintSet.h
ResponseCmd.o: ExtForceSet.h Control.h ChargeDensity.h
ResponseCmd.o: BOSampleStepper.h SampleStepper.h Timer.h EnergyFunctional.h
ResponseCmd.o: StructureFactor.h ElectricEnthalpy.h Matrix.h Context.h
ResponseCmd.o: blacs.h D3vector.h Wavefunction.h UnitCell.h SlaterDet.h
ResponseCmd.o: Basis.h Sample.h AtomSet.h Atom.h D3tensor.h blas.h
ResponseCmd.o: ConstraintSet.h ExtForceSet.h Control.h ChargeDensity.h
ResponseCmd.o: ExternalPotential.h FourierTransform.h ResponseCmd.h
ResponseCmd.o: UserInterface.h
ResponseCmd.o: UserInterface.h
RseedCmd.o: UserInterface.h Sample.h AtomSet.h Context.h blacs.h Atom.h
RseedCmd.o: D3vector.h UnitCell.h D3tensor.h blas.h ConstraintSet.h
......
......@@ -28,6 +28,7 @@ using namespace std;
#include "ResponseCmd.h"
#include "mpi.h"
////////////////////////////////////////////////////////////////////////////////
int ResponseCmd::action(int argc, char **argv)
{
// " syntax: response amplitude nitscf [nite]\n\n"
......@@ -61,13 +62,15 @@ int ResponseCmd::action(int argc, char **argv)
if ( ui->onpe0() )
cout << " ResponseCmd: start computing charge density response under "
<< " external potential from " << argv[2] << endl;
s->vext->filename = string(argv[2]);
if ( s->vext )
delete s->vext;
s->vext = new ExternalPotential(*s,argv[2]);
nitscf = atoi(argv[3]);
if ( argc == 4 )
if ( argc == 5 )
nite = atoi(argv[4]);
responseVext(nitscf,nite);
//s->vext->filename.clear();
}
else
{
......@@ -86,6 +89,7 @@ int ResponseCmd::action(int argc, char **argv)
return 0;
}
////////////////////////////////////////////////////////////////////////////////
void ResponseCmd::responseEfield(double amplitude, int nitscf, int nite)
{
// compute dipole change
......@@ -153,30 +157,28 @@ void ResponseCmd::responseEfield(double amplitude, int nitscf, int nite)
delete stepper;
}
////////////////////////////////////////////////////////////////////////////////
void ResponseCmd::responseVext(int nitscf, int nite)
{
SampleStepper* stepper = new BOSampleStepper(*s,nitscf,nite);
ChargeDensity &cd = stepper->cd();
stepper->step(0);
ChargeDensity rho1(s->wf);
rho1.update_density();
const vector<vector<double> > &rhor1 = rho1.rhor;
const vector<vector<double> > rhor1 = cd.rhor;
s->vext->reverse();
s->vext->set_amplitude(-1.0);
if ( ui->onpe0() )
cout << " ResponseCmd: external potential is reversed, "
<< " starting another scf iteration" << endl;
stepper->step(0);
ChargeDensity rho2(s->wf);
rho2.update_density();
const vector<vector<double> > &rhor2 = rho2.rhor;
const vector<vector<double> > rhor2 = cd.rhor;
// compute drho_r as rhor1 - rhor2
const int nspin = s->wf.nspin();
const int np012loc = rho1.vft()->np012loc();
const double omega = rho1.vbasis()->cell().volume();
const int np012loc = cd.vft()->np012loc();
const double omega = cd.vbasis()->cell().volume();
const double omega_inv = 1.0 / omega;
vector<vector<double> > drho_r;
vector<vector<complex<double> > > drho_r_tmp; // to be used for FT
......@@ -197,12 +199,12 @@ void ResponseCmd::responseVext(int nitscf, int nite)
// Fourier (forward) transform drho_r to the basis that is compatible
// with the grid size of the external potential.
const UnitCell& cell = rho1.vbasis()->cell();
const UnitCell& cell = cd.vbasis()->cell();
const double ecut = s->vext->ecut();
const int np0 = rho1.vft()->np0();
const int np1 = rho1.vft()->np1();
const int np2 = rho1.vft()->np2();
Basis basis(rho1.vcomm(),D3vector(0,0,0));
const int np0 = cd.vft()->np0();
const int np1 = cd.vft()->np1();
const int np2 = cd.vft()->np2();
Basis basis(cd.vcomm(),D3vector(0,0,0));
basis.resize(cell,cell,ecut);
FourierTransform ft1 (basis,np0,np1,np2);
......@@ -210,9 +212,8 @@ void ResponseCmd::responseVext(int nitscf, int nite)
drho_g.resize(nspin);
for ( int ispin = 0; ispin < nspin; ispin++ )
{
drho_g[ispin].resize(basis.localsize());
ft1.forward(&drho_r_tmp[ispin][0],&drho_g[ispin][0]);
drho_g[ispin].resize(basis.localsize());
ft1.forward(&drho_r_tmp[ispin][0],&drho_g[ispin][0]);
}
// Fourier (backward) transform drho_g to drho_r the grid
......@@ -225,25 +226,25 @@ void ResponseCmd::responseVext(int nitscf, int nite)
const int n012loc = ft2.np012loc();
for ( int ispin = 0; ispin < nspin; ispin++ )
{
drho_r[ispin].resize(n012loc);
drho_r_tmp[ispin].resize(n012loc);
ft2.backward(&drho_g[ispin][0],&drho_r_tmp[ispin][0]);
for ( int ir = 0; ir < drho_r[ispin].size(); ir++ )
{
drho_r[ispin].resize(n012loc);
drho_r_tmp[ispin].resize(n012loc);
ft2.backward(&drho_g[ispin][0],&drho_r_tmp[ispin][0]);
for ( int ir = 0; ir < drho_r[ispin].size(); ir++ )
{
drho_r[ispin][ir] = real( drho_r_tmp[ispin][ir] * omega_inv );
}
drho_r[ispin][ir] = real( drho_r_tmp[ispin][ir] * omega_inv );
}
}
MPI_Barrier(MPI_COMM_WORLD);
const Context& ctxt = cd.context();
ctxt.barrier();
if ( ui->onpe0() )
cout << " ResponseCmd: density response has been"
<< " interpolated from grid size \n"
<< " interpolated from grid size \n"
<< " (" << np0 << ", " << np1 << ", " << np2 << ")"
<< " to (" << n0 << ", " << n1 << ", " << n2 << ") \n";
<< " to (" << n0 << ", " << n1 << ", " << n2 << ") \n";
// process 0 collects the drho_r from all processors on column 0
const Context& ctxt = *(s->wf.spincontext());
if ( ctxt.onpe0() )
{
for ( int ispin = 0; ispin < nspin; ispin++ )
......@@ -306,9 +307,9 @@ void ResponseCmd::responseVext(int nitscf, int nite)
// process 0 output density difference
if ( ctxt.onpe0() )
{
ofstream os;
string filename = s->vext->filename + ".response";
os.open(filename.c_str());
ofstream os;
string filename = s->vext->filename() + ".response";
os.open(filename.c_str());
os << n0 << " " << n1 << " " << n2 << " " << endl;
for ( int i = 0; i < drho_r[0].size(); i++ )
if ( nspin == 2 )
......@@ -316,7 +317,7 @@ void ResponseCmd::responseVext(int nitscf, int nite)
else
os << drho_r[0][i] * omega / ft2.np012() << endl;
cout << " ResponseCmd: charge density response has been written in: "
<< filename << endl;
<< filename << endl;
}
delete stepper;
......
......@@ -47,12 +47,14 @@ class Vext : public Var
if ( !strcmp(argv[1],"NULL") )
{
// set vext NULL
// reset file name to empty string
s->vext->filename.clear();
delete s->vext;
s->vext = 0;
}
else
{
s->vext->filename = argv[1];
if ( s->vext )
delete s->vext;
s->vext = new ExternalPotential(*s,argv[1]);
}
return 0;
......@@ -63,12 +65,14 @@ class Vext : public Var
ostringstream st;
st.setf(ios::left,ios::adjustfield);
st << setw(10) << name() << " = ";
st.setf(ios::right,ios::adjustfield);
st << setw(10) << s->vext->filename;
if ( s->vext )
{
st.setf(ios::right,ios::adjustfield);
st << setw(10) << s->vext->filename();
}
return st.str();
}
Vext(Sample *sample) : s(sample) { s->vext = new ExternalPotential(*s); }
~Vext(void) { delete s->vext; }
Vext(Sample *sample) : s(sample) {}
};
#endif
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