Commit 44449954 by Francois Gygi

Implementation of the response command for constant field


git-svn-id: http://qboxcode.org/svn/qb/branches/vext@1854 cba15fb0-1239-40c8-b417-11db7ca47a34
parent f3773072
......@@ -15,7 +15,6 @@
// BOSampleStepper.h
//
////////////////////////////////////////////////////////////////////////////////
// $Id: BOSampleStepper.h,v 1.9 2009-09-08 05:35:39 fgygi Exp $
#ifndef BOSAMPLESTEPPER_H
#define BOSAMPLESTEPPER_H
......@@ -55,6 +54,9 @@ class BOSampleStepper : public SampleStepper
void step(int niter);
void initialize_density(void);
EnergyFunctional& ef(void) { return ef_; }
ChargeDensity& cd(void) { return cd_; }
BOSampleStepper(Sample& s, int nitscf, int nite);
~BOSampleStepper();
};
......
......@@ -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);
......
......@@ -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;
......
......@@ -50,7 +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 \
ExternalPotential.o ResponseCmd.o \
$(PLTOBJECTS)
CXXFLAGS += -DTARGET='"$(TARGET)"'
$(EXEC): $(OBJECTS)
......@@ -384,9 +384,11 @@ 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: 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
ExternalPotential.o: ExtForceSet.h Wavefunction.h Control.h ChargeDensity.h
ExternalPotential.o: Timer.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
......@@ -601,11 +603,16 @@ SampleReader.o: UnitCell.h D3tensor.h blas.h ConstraintSet.h ExtForceSet.h
SampleReader.o: Wavefunction.h Control.h SampleReader.h XMLGFPreprocessor.h
SampleReader.o: Matrix.h Timer.h SampleHandler.h StructureHandler.h
SampleReader.o: StructuredDocumentHandler.h StrX.h
SampleStepper.o: SampleStepper.h Timer.h Sample.h AtomSet.h Context.h blacs.h
SampleStepper.o: Atom.h D3vector.h UnitCell.h D3tensor.h blas.h
SampleStepper.o: ConstraintSet.h ExtForceSet.h Wavefunction.h Control.h
SampleStepper.o: Species.h
SampleStepper.o: Timer.h
SampleStepper.o: SampleStepper.h Timer.h EnergyFunctional.h StructureFactor.h
SampleStepper.o: ElectricEnthalpy.h Matrix.h Context.h blacs.h D3vector.h
SampleStepper.o: Wavefunction.h UnitCell.h SlaterDet.h Basis.h Sample.h
SampleStepper.o: AtomSet.h Atom.h D3tensor.h blas.h ConstraintSet.h
SampleStepper.o: ExtForceSet.h Control.h ChargeDensity.h Species.h
SampleStepper.o: Timer.h EnergyFunctional.h StructureFactor.h
SampleStepper.o: ElectricEnthalpy.h Matrix.h Context.h blacs.h D3vector.h
SampleStepper.o: Wavefunction.h UnitCell.h SlaterDet.h Basis.h Sample.h
SampleStepper.o: AtomSet.h Atom.h D3tensor.h blas.h ConstraintSet.h
SampleStepper.o: ExtForceSet.h Control.h ChargeDensity.h
SampleWriter.o: SampleWriter.h Context.h blacs.h Sample.h AtomSet.h Atom.h
SampleWriter.o: D3vector.h UnitCell.h D3tensor.h blas.h ConstraintSet.h
SampleWriter.o: ExtForceSet.h Wavefunction.h Control.h qbox_xmlns.h Timer.h
......@@ -682,7 +689,7 @@ 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
Vext.o: Control.h ExternalPotential.h ChargeDensity.h Timer.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
......@@ -731,15 +738,16 @@ qb.o: BasisMapping.h ConstraintCmd.h DistanceCmd.h ExtForceCmd.h
qb.o: FoldInWsCmd.h HelpCmd.h KpointCmd.h ListAtomsCmd.h ListSpeciesCmd.h
qb.o: LoadCmd.h MoveCmd.h PartialChargeCmd.h PlotCmd.h PrintCmd.h QuitCmd.h
qb.o: RandomizeRCmd.h RandomizeVCmd.h RandomizeWfCmd.h ResetVcmCmd.h
qb.o: RescaleVCmd.h RseedCmd.h RunCmd.h SaveCmd.h SetCmd.h SetVelocityCmd.h
qb.o: SpeciesCmd.h StatusCmd.h StrainCmd.h TorsionCmd.h BisectionCmd.h
qb.o: Bisection.h SlaterDet.h Basis.h Matrix.h AlphaPBE0.h AtomsDyn.h BlHF.h
qb.o: BtHF.h Cell.h CellDyn.h CellLock.h CellMass.h ChargeMixCoeff.h
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: Vext.h ExternalPotential.h WfDiag.h WfDyn.h Xc.h
qb.o: RescaleVCmd.h ResponseCmd.h RseedCmd.h RunCmd.h SaveCmd.h SetCmd.h
qb.o: SetVelocityCmd.h SpeciesCmd.h StatusCmd.h StrainCmd.h TorsionCmd.h
qb.o: BisectionCmd.h Bisection.h SlaterDet.h Basis.h Matrix.h AlphaPBE0.h
qb.o: AtomsDyn.h BlHF.h BtHF.h Cell.h CellDyn.h CellLock.h CellMass.h
qb.o: ChargeMixCoeff.h ChargeMixNdim.h ChargeMixRcut.h Debug.h Dspin.h Ecut.h
qb.o: Ecutprec.h Ecuts.h Efield.h Polarization.h Emass.h ExtStress.h
qb.o: FermiTemp.h IterCmd.h IterCmdPeriod.h Dt.h Nempty.h NetCharge.h
qb.o: Nrowmax.h Nspin.h RefCell.h ScfTol.h Stress.h Thermostat.h ThTemp.h
qb.o: ThTime.h ThWidth.h Vext.h ExternalPotential.h ChargeDensity.h WfDiag.h
qb.o: WfDyn.h Xc.h
qbox_xmlns.o: qbox_xmlns.h
release.o: release.h
sinft.o: sinft.h
......
......@@ -18,6 +18,7 @@
#include "ResponseCmd.h"
#include<iostream>
#include<sstream>
using namespace std;
#include "BOSampleStepper.h"
......@@ -26,10 +27,13 @@ using namespace std;
int ResponseCmd::action(int argc, char **argv)
{
if ( argc < 3 || argc > 4)
// " syntax: response amplitude nitscf [nite]\n\n"
// " syntax: response -vext vext_file nitscf [nite]\n\n"
if ( argc < 2 || argc > 5)
{
if ( ui->onpe0() )
cout << " use: response vext_file nitscf [nite]" << endl;
cout << " use: response amplitude nitscf [nite]" << endl
<< " response [-vext vext_file] nitscf [nite]" << endl;
return 1;
}
......@@ -46,39 +50,86 @@ int ResponseCmd::action(int argc, char **argv)
return 1;
}
string vext_filename(argv[1]);
int niter = atoi(argv[2]);
double amplitude = 0.0;
bool vext = false;
string vext_filename;
int nitscf = 0;
int nite = 0;
if ( argc == 4 )
if ( !strcmp(argv[1],"-vext") )
{
assert(argc==4 || argc==5);
vext = true;
vext_filename = atoi(argv[2]);
nitscf = atoi(argv[3]);
if ( argc == 5 )
nite = atoi(argv[4]);
if ( ui->onpe0() )
cout << " ResponseCmd: -vext option not implemented" << endl;
return 1;
}
else
{
// nitscf nite
nite = atoi(argv[3]);
assert(argc==3 || argc==4);
amplitude = atof(argv[1]);
if ( amplitude == 0.0 )
{
if ( ui->onpe0() )
cout << " ResponseCmd: amplitude is 0.0" << endl;
return 1;
}
nitscf = atoi(argv[2]);
if ( argc == 4 )
nite = atoi(argv[3]);
}
// compute dipole change
string cmd;
istringstream cmdstream(cmd);
cmdstream.clear();
cmdstream.str("print polarization\n");
ui->processCmds(cmdstream,"[response_cmd]",true);
SampleStepper* stepper = new BOSampleStepper(*s,nitscf,nite);
assert(stepper!=0);
ElectricEnthalpy* el_enth = stepper->ef().el_enth();
// 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());
D3vector e_field[3] = { D3vector(amplitude,0.0,0.0),
D3vector(0.0,amplitude,0.0),
D3vector(0.0,0.0,amplitude) };
D3vector dipole_p[3], dipole_m[3];
stepper->step(0);
// compute density
// compute change in dipole in 3 directions by finite difference
for ( int idir = 0; idir < 3; idir++ )
{
el_enth->set_e_field(e_field[idir]);
stepper->step(0);
dipole_p[idir] = el_enth->dipole_total();
// change sign of potential
el_enth->set_e_field(-e_field[idir]);
stepper->step(0);
dipole_m[idir] = el_enth->dipole_total();
}
stepper->step(0);
// compute density
D3vector ddx = 0.5 * (dipole_p[0]-dipole_m[0])/amplitude;
D3vector ddy = 0.5 * (dipole_p[1]-dipole_m[1])/amplitude;
D3vector ddz = 0.5 * (dipole_p[2]-dipole_m[2])/amplitude;
// compute density difference
// output density difference
if ( ui->onpe0() )
{
cout << "<polarizability>" << endl;
cout << setprecision(6);
cout << " <a_xx> " << setw(10) << ddx.x << " </a_xx>"
" <a_yx> " << setw(10) << ddx.y << " </a_yx>"
" <a_zx> " << setw(10) << ddx.z << " </a_zx>" << endl;
cout << " <a_xy> " << setw(10) << ddy.x << " </a_xy>"
" <a_yy> " << setw(10) << ddy.y << " </a_yy>"
" <a_zy> " << setw(10) << ddy.z << " </a_zy>" << endl;
cout << " <a_xz> " << setw(10) << ddz.x << " </a_xz>"
" <a_yz> " << setw(10) << ddz.y << " </a_yz>"
" <a_zz> " << setw(10) << ddz.z << " </a_zz>" << endl;
cout << "</polarizability>" << endl;
}
delete stepper;
......
......@@ -38,9 +38,12 @@ class ResponseCmd : public Cmd
{
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";
" syntax: response amplitude nitscf [nite]\n\n"
" syntax: response -vext vext_file nitscf [nite]\n\n"
" The response command computes the first order change in dipole caused\n"
" by an external electric field defined by the amplitude argument.\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\n";
}
int action(int argc, char **argv);
......
......@@ -15,12 +15,13 @@
// SampleStepper.h
//
////////////////////////////////////////////////////////////////////////////////
// $Id: SampleStepper.h,v 1.14 2009-09-08 05:38:31 fgygi Exp $
#ifndef SAMPLESTEPPER_H
#define SAMPLESTEPPER_H
#include "Timer.h"
#include "EnergyFunctional.h"
#include "ChargeDensity.h"
#include <map>
#include <string>
#include <vector>
......@@ -49,6 +50,8 @@ class SampleStepper
void print_stress(void);
void compute_sigma(void); // compute kinetic contribution to stress
virtual void initialize_density() {}
virtual EnergyFunctional& ef(void) {}
virtual ChargeDensity& cd(void) {}
SampleStepper(Sample& s);
virtual ~SampleStepper(void);
......
......@@ -65,6 +65,7 @@ using namespace std;
#include "RandomizeWfCmd.h"
#include "ResetVcmCmd.h"
#include "RescaleVCmd.h"
#include "ResponseCmd.h"
#include "RseedCmd.h"
#include "RunCmd.h"
#include "SaveCmd.h"
......@@ -281,6 +282,7 @@ int main(int argc, char **argv, char **envp)
ui.addCmd(new RandomizeVCmd(s));
ui.addCmd(new RandomizeWfCmd(s));
ui.addCmd(new RescaleVCmd(s));
ui.addCmd(new ResponseCmd(s));
ui.addCmd(new ResetVcmCmd(s));
ui.addCmd(new RseedCmd(s));
ui.addCmd(new RunCmd(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