Commit 6702a3e0 by Francois Gygi

Prepare for metaGGA, SCAN functional implem.

parent 569b78be
......@@ -44,6 +44,7 @@ class B3LYPFunctional : public XCFunctional
B3LYPFunctional(const std::vector<std::vector<double> > &rhoe);
bool isGGA() const { return true; };
bool isMeta() const { return false; };
std::string name() const { return "B3LYP"; };
void setxc(void);
};
......
......@@ -49,6 +49,7 @@ class BLYPFunctional : public XCFunctional
double *vc2_upup, double *vc2_dndn, double *vc2_updn, double *vc2_dnup);
bool isGGA() const { return true; };
bool isMeta() const { return false; };
std::string name() const { return "BLYP"; };
void setxc(void);
};
......
......@@ -63,6 +63,7 @@ class LDAFunctional : public XCFunctional
};
bool isGGA() const { return false; };
bool isMeta() const { return false; };
std::string name() const { return "LDA"; };
void setxc(void);
};
......
......@@ -30,6 +30,7 @@ OBJECTS=qb.o AtomSet.o Atom.o Species.o \
XCOperator.o ExchangeOperator.o Bisection.o \
XCPotential.o LDAFunctional.o VWNFunctional.o \
PBEFunctional.o BLYPFunctional.o B3LYPFunctional.o \
SCANFunctional.o \
NonLocalPotential.o SampleReader.o StructuredDocumentHandler.o \
SampleHandler.o AtomSetHandler.o WavefunctionHandler.o \
XMLGFPreprocessor.o Base64Transcoder.o \
......@@ -560,6 +561,8 @@ RunCmd.o: Context.h blacs.h D3vector.h Wavefunction.h UnitCell.h SlaterDet.h
RunCmd.o: Basis.h Sample.h AtomSet.h Atom.h D3tensor.h blas.h ConstraintSet.h
RunCmd.o: ExtForceSet.h Control.h ChargeDensity.h CPSampleStepper.h
RunCmd.o: UserInterface.h
SCANFunctional.o: SCANFunctional.h XCFunctional.h
SCANFunctional.o: XCFunctional.h
SDAIonicStepper.o: SDAIonicStepper.h IonicStepper.h Sample.h AtomSet.h
SDAIonicStepper.o: Context.h blacs.h Atom.h D3vector.h UnitCell.h D3tensor.h
SDAIonicStepper.o: blas.h ConstraintSet.h ExtForceSet.h Wavefunction.h
......@@ -702,13 +705,12 @@ XCOperator.o: D3vector.h UnitCell.h D3tensor.h blas.h ConstraintSet.h
XCOperator.o: ExtForceSet.h Wavefunction.h Control.h ChargeDensity.h Timer.h
XCOperator.o: XCPotential.h ExchangeOperator.h SlaterDet.h Basis.h Matrix.h
XCOperator.o: FourierTransform.h
XCOperator.o: Sample.h AtomSet.h Context.h blacs.h Atom.h D3vector.h
XCOperator.o: UnitCell.h D3tensor.h blas.h ConstraintSet.h ExtForceSet.h
XCOperator.o: Wavefunction.h Control.h
XCPotential.o: XCPotential.h Control.h D3vector.h ChargeDensity.h Timer.h
XCPotential.o: Context.h blacs.h LDAFunctional.h XCFunctional.h
XCPotential.o: VWNFunctional.h PBEFunctional.h BLYPFunctional.h
XCPotential.o: B3LYPFunctional.h Basis.h UnitCell.h FourierTransform.h blas.h
XCPotential.o: B3LYPFunctional.h SCANFunctional.h Sample.h AtomSet.h Atom.h
XCPotential.o: UnitCell.h D3tensor.h blas.h ConstraintSet.h ExtForceSet.h
XCPotential.o: Wavefunction.h Basis.h FourierTransform.h
XCPotential.o: Control.h D3vector.h ChargeDensity.h Timer.h Context.h blacs.h
XMLGFPreprocessor.o: Timer.h Context.h blacs.h Base64Transcoder.h Matrix.h
XMLGFPreprocessor.o: XMLGFPreprocessor.h
......
......@@ -51,6 +51,7 @@ class PBEFunctional : public XCFunctional
double x_coeff=1.0, double c_coeff=1.0);
bool isGGA() const { return true; };
bool isMeta() const { return false; };
std::string name() const { return "PBE"; };
void setxc(void);
};
......
......@@ -75,6 +75,7 @@ class VWNFunctional : public XCFunctional
double &ec, double &vc_up, double &vc_dn);
bool isGGA() const { return false; };
bool isMeta() const { return false; };
std::string name() const { return "VWN"; };
void setxc(void);
};
......
......@@ -65,11 +65,14 @@ class XCFunctional
const double *rho, *rho_up, *rho_dn;
double *grad_rho[3], *grad_rho_up[3], *grad_rho_dn[3];
double *tau;
double *exc, *exc_up, *exc_dn;
double *vxc1, *vxc1_up, *vxc1_dn;
double *vxc2, *vxc2_upup, *vxc2_dndn, *vxc2_updn, *vxc2_dnup;
double *vxc3;
virtual bool isGGA(void) const = 0;
virtual bool isMeta(void) const = 0;
virtual std::string name(void) const = 0;
int np(void) const { return _np; };
int nspin(void) const { return _nspin; };
......@@ -80,9 +83,11 @@ class XCFunctional
grad_rho[0] = grad_rho[1] = grad_rho[2] = 0;
grad_rho_up[0] = grad_rho_up[1] = grad_rho_up[2] = 0;
grad_rho_dn[0] = grad_rho_dn[1] = grad_rho_dn[2] = 0;
tau = 0;
exc = exc_up = exc_dn = 0;
vxc1 = vxc1_up = vxc1_dn = 0;
vxc2 = vxc2_upup = vxc2_dndn = vxc2_updn = vxc2_dnup = 0;
vxc3 = 0;
}
// virtual destructor needed to ensure proper deallocation
......
......@@ -16,19 +16,22 @@
//
////////////////////////////////////////////////////////////////////////////////
#include "XCOperator.h"
#include "Sample.h"
#include "ChargeDensity.h"
#include "XCPotential.h"
#include "ExchangeOperator.h"
using namespace std;
////////////////////////////////////////////////////////////////////////////////
XCOperator::XCOperator(Sample& s, const ChargeDensity& cd) :cd_(cd)
XCOperator::XCOperator(Sample& s, const ChargeDensity& cd) :cd_(cd), s_(s)
{
// set initial values
xcp_ = 0;
xop_ = 0;
exc_ = 0.0 ;
dxc_ = 0.0 ;
hasHF_ = false;
hasMeta_ = false;
sigma_exc_.resize(6);
......@@ -41,10 +44,9 @@ XCOperator::XCOperator(Sample& s, const ChargeDensity& cd) :cd_(cd)
( functional_name == "BLYP" ) )
{
// create only an xc potential
xcp_ = new XCPotential(cd, functional_name, s.ctrl);
xcp_ = new XCPotential(cd, functional_name, s_);
hasPotential_ = true;
hasGGA_ = xcp_->isGGA();
hasHF_ = false;
HFmixCoeff_ = 0.0;
}
else if ( functional_name == "HF" )
......@@ -59,27 +61,35 @@ XCOperator::XCOperator(Sample& s, const ChargeDensity& cd) :cd_(cd)
else if ( functional_name == "PBE0" )
{
// create an exchange potential
xcp_ = new XCPotential(cd, functional_name, s.ctrl);
xcp_ = new XCPotential(cd, functional_name, s_);
// create the exchange operator with mixing coeff=0.25
xop_ = new ExchangeOperator(s, s.ctrl.alpha_PBE0);
xop_ = new ExchangeOperator(s_, s_.ctrl.alpha_PBE0);
hasPotential_ = true;
hasGGA_ = xcp_->isGGA();
hasHF_ = true;
HFmixCoeff_ = s.ctrl.alpha_PBE0;;
HFmixCoeff_ = s_.ctrl.alpha_PBE0;;
}
else if ( functional_name == "B3LYP" )
{
// create an exchange potential
xcp_ = new XCPotential(cd, functional_name, s.ctrl);
xcp_ = new XCPotential(cd, functional_name, s_);
// create the exchange operator with mixing coeff=0.20
xop_ = new ExchangeOperator(s, 0.20);
xop_ = new ExchangeOperator(s_, 0.20);
hasPotential_ = true;
hasGGA_ = xcp_->isGGA();
hasHF_ = true;
HFmixCoeff_ = 0.20;
}
else if ( functional_name == "SCAN" )
{
// create an exchange potential
xcp_ = new XCPotential(cd, functional_name, s_);
hasPotential_ = true;
hasGGA_ = xcp_->isGGA();
hasMeta_ = true;
}
else
{
throw XCOperatorException(
......@@ -101,7 +111,7 @@ void XCOperator::update(std::vector<std::vector<double> >& vr, bool compute_stre
// compute vxc potential and energy
if ( hasPotential_ )
{
// update LDA/GGA xc potential
// update LDA/GGA/MetaGGA xc potential
xcp_->update( vr );
// LDA/GGA exchange energy
......@@ -133,6 +143,8 @@ void XCOperator::apply_self_energy(Wavefunction &dwf)
{
if ( hasHF() )
xop_->apply_operator(dwf);
if ( hasMeta() )
xcp_->apply_meta_operator(dwf);
}
////////////////////////////////////////////////////////////////////////////////
......
......@@ -19,9 +19,12 @@
#ifndef XCOPERATOR_H
#define XCOPERATOR_H
#include "Sample.h"
#include <valarray>
#include <vector>
#include <string>
class Sample;
class Wavefunction;
class ChargeDensity;
class XCPotential;
class ExchangeOperator;
......@@ -32,6 +35,7 @@ class XCOperator
XCPotential* xcp_;
ExchangeOperator* xop_;
const Sample& s_;
const ChargeDensity& cd_;
double HFmixCoeff_ ;
double exc_; // XC energy: includes local and HF terms
......@@ -42,6 +46,7 @@ class XCOperator
bool hasPotential_;
bool hasGGA_;
bool hasHF_;
bool hasMeta_;
public:
......@@ -59,6 +64,7 @@ class XCOperator
bool hasGGA(void) { return hasGGA_; };
bool hasHF(void) { return hasHF_; };
bool hasMeta(void) { return hasMeta_; };
void update(std::vector<std::vector<double> >& vr, bool compute_stress);
void apply_self_energy(Wavefunction &dwf);
......
......@@ -22,6 +22,8 @@
#include "PBEFunctional.h"
#include "BLYPFunctional.h"
#include "B3LYPFunctional.h"
#include "SCANFunctional.h"
#include "Sample.h"
#include "Basis.h"
#include "FourierTransform.h"
#include "blas.h" // daxpy, dcopy
......@@ -30,7 +32,7 @@ using namespace std;
////////////////////////////////////////////////////////////////////////////////
XCPotential::XCPotential(const ChargeDensity& cd, const string functional_name,
const Control& ctrl): cd_(cd), vft_(*cd_.vft()), vbasis_(*cd_.vbasis())
const Sample& s): cd_(cd), vft_(*cd_.vft()), vbasis_(*cd_.vbasis()), s_(s)
{
if ( functional_name == "LDA" )
{
......@@ -50,7 +52,7 @@ XCPotential::XCPotential(const ChargeDensity& cd, const string functional_name,
}
else if ( functional_name == "PBE0" )
{
const double x_coeff = 1.0 - ctrl.alpha_PBE0;
const double x_coeff = 1.0 - s_.ctrl.alpha_PBE0;
const double c_coeff = 1.0;
xcf_ = new PBEFunctional(cd_.rhor,x_coeff,c_coeff);
}
......@@ -58,6 +60,10 @@ XCPotential::XCPotential(const ChargeDensity& cd, const string functional_name,
{
xcf_ = new B3LYPFunctional(cd_.rhor);
}
else if ( functional_name == "SCAN" )
{
xcf_ = new SCANFunctional(cd_.rhor);
}
else
{
throw XCPotentialException("unknown functional name");
......@@ -91,6 +97,12 @@ bool XCPotential::isGGA(void)
}
////////////////////////////////////////////////////////////////////////////////
bool XCPotential::isMeta(void)
{
return xcf_->isMeta();
}
////////////////////////////////////////////////////////////////////////////////
void XCPotential::update(vector<vector<double> >& vr)
{
// compute exchange-correlation energy and add vxc potential to vr[ispin][ir]
......@@ -216,6 +228,11 @@ void XCPotential::update(vector<vector<double> >& vr)
} // j
}
if ( xcf_->isMeta() )
{
// compute tau
}
xcf_->setxc();
// compute xc potential
......@@ -526,3 +543,9 @@ void XCPotential::compute_stress(valarray<double>& sigma_exc)
sigma_exc[5] = tsum[5];
}
}
////////////////////////////////////////////////////////////////////////////////
void XCPotential::apply_meta_operator(Wavefunction& dwf)
{
// apply meta operator to dwf using xcf_->vxc3
}
......@@ -28,11 +28,14 @@
class Basis;
class FourierTransform;
class XCFunctional;
class Sample;
class Wavefunction;
class XCPotential
{
private:
const Sample& s_;
const ChargeDensity& cd_;
XCFunctional* xcf_;
......@@ -52,11 +55,13 @@ class XCPotential
const XCFunctional* xcf() { return xcf_; }
bool isGGA(void);
bool isMeta(void);
XCPotential(const ChargeDensity& cd, const std::string functional_name,
const Control& ctrl);
const Sample& s);
~XCPotential();
void update(std::vector<std::vector<double> >& vr);
void compute_stress(std::valarray<double>& sigma_exc);
void apply_meta_operator(Wavefunction& dwf);
double exc(void) { return exc_; }
double dxc(void) { return dxc_; }
};
......
......@@ -50,10 +50,12 @@ class Xc : public Var
v == "BLYP" ||
v == "HF" ||
v == "PBE0" ||
v == "B3LYP" ) )
v == "B3LYP" ||
v == "SCAN" ) )
{
if ( ui->onpe0() )
cout << " xc must be LDA, VWN, PBE, BLYP, HF, PBE0 or B3LYP" << endl;
cout << " xc must be LDA, VWN, PBE, BLYP, HF, PBE0, B3LYP or SCAN"
<< endl;
return 1;
}
......
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