Commit f3773072 by Francois Gygi

Implementation of the external potential interface


git-svn-id: http://qboxcode.org/svn/qb/branches/vext@1853 cba15fb0-1239-40c8-b417-11db7ca47a34
parent 5000b96f
......@@ -75,5 +75,7 @@ struct Control
std::string iter_cmd;
int iter_cmd_period;
std::string vext;
};
#endif
......@@ -30,6 +30,7 @@
#include "ConfinementPotential.h"
#include "D3vector.h"
#include "ElectricEnthalpy.h"
#include "ExternalPotential.h"
#include "Timer.h"
#include "blas.h"
......@@ -79,6 +80,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();
......@@ -336,6 +341,13 @@ void EnergyFunctional::update_vhxc(bool compute_stress)
// FT to tmpr_r
vft->backward(&vlocal_g[0],&tmp_r[0]);
// 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);
}
// add local potential in tmp_r to v_r[ispin][i]
// v_r contains the xc potential
const int size = tmp_r.size();
......
////////////////////////////////////////////////////////////////////////////////
//
// 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 "ExternalPotential.h"
#include "FourierTransform.h"
#include "Basis.h"
#include <iostream>
#include <fstream>
#include <cassert>
using namespace std;
void ExternalPotential::update(const ChargeDensity& cd)
{
// read vext on nrow=0 tasks
vector<double> vext_read;
FourierTransform *vft = cd.vft();
const int myrow = cd.context().myrow();
if ( myrow == 0 )
{
// read the external potential from file s_.ctrl.vext
if ( cd.context().onpe0() )
cout << "ExternalPotential::update: vext file: " << s_.ctrl.vext << endl;
ifstream vfile(s_.ctrl.vext.c_str());
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];
}
// 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);
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;
double ecut = min(min(ecut0,ecut1),ecut2);
if ( cd.context().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;
}
Basis basis(cd.vcomm(),D3vector(0,0,0));
basis.resize(cell,cell,ecut);
if ( cd.context().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;
}
// interpolate on grid compatible with the charge density cd
FourierTransform ft(basis,n_[0],n_[1],n_[2]);
vector<complex<double> > vext_read_g(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];
for ( int i = 0; i < tmp_r.size(); i++ )
tmp_r[i] = vext_read[istart+i];
// compute Fourier coefficients
ft.forward(&tmp_r[0],&vext_read_g[0]);
// vext_read_g now contains the Fourier coefficients of vext
// interpolate to vft grid
FourierTransform ft2(basis,vft->np0(),vft->np1(),vft->np2());
tmp_r.resize(vft->np012loc());
vext_r_.resize(tmp_r.size());
ft2.backward(&vext_read_g[0],&tmp_r[0]);
for ( int i = 0; i < vext_r_.size(); i++ )
vext_r_[i] = real(tmp_r[i]);
}
////////////////////////////////////////////////////////////////////////////////
//
// 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 "Sample.h"
#include "ChargeDensity.h"
#include <vector>
#include <string>
class ExternalPotential
{
private:
Sample& s_;
int n_[3];
std::vector<double> vext_r_;
std::vector<complex<double> > vext_g_;
public:
ExternalPotential(Sample& s): s_(s) {}
~ExternalPotential() {}
double v(size_t i) const { return vext_r_[i]; }
void update(const ChargeDensity& cd);
};
#endif
......@@ -50,6 +50,7 @@ OBJECTS=qb.o AtomSet.o Atom.o Species.o \
GlobalExtForce.o \
uuid_str.o sampling.o CGOptimizer.o LineMinimizer.o \
ElectricEnthalpy.o PartialChargeCmd.o \
ExternalPotential.o \
$(PLTOBJECTS)
CXXFLAGS += -DTARGET='"$(TARGET)"'
$(EXEC): $(OBJECTS)
......@@ -355,6 +356,7 @@ EnergyFunctional.o: UnitCell.h SlaterDet.h Basis.h Timer.h Sample.h AtomSet.h
EnergyFunctional.o: Atom.h D3tensor.h blas.h ConstraintSet.h ExtForceSet.h
EnergyFunctional.o: Control.h Species.h ChargeDensity.h FourierTransform.h
EnergyFunctional.o: XCOperator.h NonLocalPotential.h ConfinementPotential.h
EnergyFunctional.o: ExternalPotential.h
EnergyFunctional.o: StructureFactor.h ElectricEnthalpy.h Matrix.h Context.h
EnergyFunctional.o: blacs.h D3vector.h Wavefunction.h UnitCell.h SlaterDet.h
EnergyFunctional.o: Basis.h Timer.h Sample.h AtomSet.h Atom.h D3tensor.h
......@@ -379,6 +381,12 @@ 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: 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
FermiTemp.o: Sample.h AtomSet.h Context.h blacs.h Atom.h D3vector.h
FermiTemp.o: UnitCell.h D3tensor.h blas.h ConstraintSet.h ExtForceSet.h
FermiTemp.o: Wavefunction.h Control.h
......@@ -537,6 +545,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: 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
RseedCmd.o: ExtForceSet.h Wavefunction.h Control.h
......@@ -665,6 +680,9 @@ UnitCell.o: D3vector.h
UserInterface.o: UserInterface.h qbox_xmlns.h
VWNFunctional.o: VWNFunctional.h XCFunctional.h
VWNFunctional.o: XCFunctional.h
Vext.o: Sample.h AtomSet.h Context.h blacs.h Atom.h D3vector.h UnitCell.h
Vext.o: D3tensor.h blas.h ConstraintSet.h ExtForceSet.h Wavefunction.h
Vext.o: Control.h ExternalPotential.h
Wavefunction.o: Wavefunction.h D3vector.h UnitCell.h SlaterDet.h Context.h
Wavefunction.o: blacs.h Basis.h Matrix.h Timer.h jacobi.h SharedFilePtr.h
Wavefunction.o: D3vector.h UnitCell.h
......@@ -721,7 +739,7 @@ qb.o: ChargeMixNdim.h ChargeMixRcut.h Debug.h Dspin.h Ecut.h Ecutprec.h
qb.o: Ecuts.h Efield.h Polarization.h Emass.h ExtStress.h FermiTemp.h
qb.o: IterCmd.h IterCmdPeriod.h Dt.h Nempty.h NetCharge.h Nrowmax.h Nspin.h
qb.o: RefCell.h ScfTol.h Stress.h Thermostat.h ThTemp.h ThTime.h ThWidth.h
qb.o: WfDiag.h WfDyn.h Xc.h
qb.o: Vext.h ExternalPotential.h WfDiag.h WfDyn.h Xc.h
qbox_xmlns.o: qbox_xmlns.h
release.o: release.h
sinft.o: sinft.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/>.
//
////////////////////////////////////////////////////////////////////////////////
//
// ResponseCmd.C:
//
////////////////////////////////////////////////////////////////////////////////
#include "ResponseCmd.h"
#include<iostream>
using namespace std;
#include "BOSampleStepper.h"
#include<ctime>
#include<cassert>
int ResponseCmd::action(int argc, char **argv)
{
if ( argc < 3 || argc > 4)
{
if ( ui->onpe0() )
cout << " use: response vext_file nitscf [nite]" << endl;
return 1;
}
if ( s->wf.nst() == 0 )
{
if ( ui->onpe0() )
cout << " ResponseCmd: no states, cannot run" << endl;
return 1;
}
if ( s->wf.ecut() == 0.0 )
{
if ( ui->onpe0() )
cout << " ResponseCmd: ecut = 0.0, cannot run" << endl;
return 1;
}
string vext_filename(argv[1]);
int niter = atoi(argv[2]);
int nite = 0;
if ( argc == 4 )
{
// nitscf nite
nite = atoi(argv[3]);
}
SampleStepper* stepper = new BOSampleStepper(*s,nitscf,nite);
assert(stepper!=0);
// update external potential from file
if ( !s.vext )
{
if ( ui->onpe0() )
cout << " ResponseCmd: vext file not defined" << endl;
return 1;
}
ifstream vextfile(s.vext.c_str());
stepper->step(0);
// compute density
// change sign of potential
stepper->step(0);
// compute density
// compute density difference
// output density difference
delete stepper;
return 0;
}
////////////////////////////////////////////////////////////////////////////////
//
// 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/>.
//
////////////////////////////////////////////////////////////////////////////////
//
// ResponseCmd.h:
//
////////////////////////////////////////////////////////////////////////////////
#ifndef RESPONSECMD_H
#define RESPONSECMD_H
#include <iostream>
#include "UserInterface.h"
class Sample;
class ResponseCmd : public Cmd
{
private:
public:
Sample *s;
ResponseCmd(Sample *sample) : s(sample) {};
const char *name(void) const { return "response"; }
const char *help_msg(void) const
{
return
"\n response\n\n"
" syntax: response vext_file nitscf [nite]\n\n"
" The response command computes the response to the external potential\n"
" defined in the file vext_file.\n\n";
}
int action(int argc, char **argv);
};
#endif
......@@ -27,6 +27,7 @@
class Context;
class UserInterface;
class ExternalPotential;
class Sample
{
......@@ -39,6 +40,7 @@ class Sample
AtomSet atoms;
ConstraintSet constraints;
ExtForceSet extforces;
ExternalPotential* vext;
Wavefunction wf;
Wavefunction* wfv; // wavefunction velocity
Control ctrl;
......@@ -46,7 +48,7 @@ class Sample
Sample(const Context& ctxt, UserInterface *ui_ = 0) : ctxt_(ctxt), ui(ui_),
atoms(ctxt), constraints(ctxt),
extforces(ctxt), wf(ctxt), wfv(0) {}
extforces(ctxt), vext(0), wf(ctxt), wfv(0) {}
~Sample(void) { delete wfv; }
void reset(void)
{
......
////////////////////////////////////////////////////////////////////////////////
//
// 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/>.
//
////////////////////////////////////////////////////////////////////////////////
//
// Vext.h
//
////////////////////////////////////////////////////////////////////////////////
#ifndef VEXT_H
#define VEXT_H
#include<iostream>
#include<iomanip>
#include<sstream>
#include<stdlib.h>
#include "Sample.h"
#include "ExternalPotential.h"
class Vext : public Var
{
Sample *s;
public:
const char *name ( void ) const { return "vext"; };
int set ( int argc, char **argv )
{
if ( argc > 3 )
{
if ( ui->onpe0() )
cout << " vext takes only one value" << endl;
return 1;
}
if ( !strcmp(argv[1],"NULL") )
{
// set vext NULL
// reset file name to empty string
delete s->vext;
s->vext = 0;
s->ctrl.vext.clear();
}
else
{
s->ctrl.vext = argv[1];
s->vext = new ExternalPotential(*s);
}
return 0;
}
string print (void) const
{
ostringstream st;
st.setf(ios::left,ios::adjustfield);
st << setw(10) << name() << " = ";
st.setf(ios::right,ios::adjustfield);
st << setw(10) << s->ctrl.vext;
return st.str();
}
Vext(Sample *sample) : s(sample) { s->ctrl.vext = ""; };
};
#endif
......@@ -111,6 +111,7 @@ using namespace std;
#include "ThTemp.h"
#include "ThTime.h"
#include "ThWidth.h"
#include "Vext.h"
#include "WfDiag.h"
#include "WfDyn.h"
#include "Xc.h"
......@@ -326,6 +327,7 @@ int main(int argc, char **argv, char **envp)
ui.addVar(new ThTemp(s));
ui.addVar(new ThTime(s));
ui.addVar(new ThWidth(s));
ui.addVar(new Vext(s));
ui.addVar(new WfDiag(s));
ui.addVar(new WfDyn(s));
ui.addVar(new Xc(s));
......
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