Commit 3d1dd7d1 by Francois Gygi

merged hfb branch

git-svn-id: http://qboxcode.org/svn/qb/trunk@1115 cba15fb0-1239-40c8-b417-11db7ca47a34
parent b54d90aa
......@@ -15,12 +15,12 @@
// XCPotential.C
//
////////////////////////////////////////////////////////////////////////////////
// $Id: XCPotential.C,v 1.9 2008-09-08 15:56:19 fgygi Exp $
#include "XCPotential.h"
#include "LDAFunctional.h"
#include "PBEFunctional.h"
#include "BLYPFunctional.h"
#include "B3LYPFunctional.h"
#include "Basis.h"
#include "FourierTransform.h"
#include "blas.h" // daxpy, dcopy
......@@ -43,6 +43,16 @@ cd_(cd), ctxt_(cd.vcontext()), vft_(*cd_.vft()), vbasis_(*cd_.vbasis())
{
xcf_ = new BLYPFunctional(cd_.rhor);
}
else if ( functional_name == "PBE0" )
{
const double x_coeff = 0.75;
const double c_coeff = 1.0;
xcf_ = new PBEFunctional(cd_.rhor,x_coeff,c_coeff);
}
else if ( functional_name == "B3LYP" )
{
xcf_ = new B3LYPFunctional(cd_.rhor);
}
else
{
throw XCPotentialException("unknown functional name");
......@@ -70,6 +80,12 @@ XCPotential::~XCPotential(void)
}
////////////////////////////////////////////////////////////////////////////////
bool XCPotential::isGGA(void)
{
return xcf_->isGGA();
}
////////////////////////////////////////////////////////////////////////////////
void XCPotential::update(vector<vector<double> >& vr)
{
// compute exchange-correlation energy and add vxc potential to vr[ispin][ir]
......@@ -276,7 +292,6 @@ void XCPotential::update(vector<vector<double> >& vr)
{
const double *const e = xcf_->exc;
const double *const v1 = xcf_->vxc1;
//const double *const v2 = xcf_->vxc2;
const double *const rh = xcf_->rho;
{
for ( int ir = 0; ir < np012loc_; ir++ )
......@@ -290,10 +305,6 @@ void XCPotential::update(vector<vector<double> >& vr)
{
const double *const v1_up = xcf_->vxc1_up;
const double *const v1_dn = xcf_->vxc1_dn;
//const double *const v2_upup = xcf_->vxc2_upup;
//const double *const v2_updn = xcf_->vxc2_updn;
//const double *const v2_dnup = xcf_->vxc2_dnup;
//const double *const v2_dndn = xcf_->vxc2_dndn;
const double *const eup = xcf_->exc_up;
const double *const edn = xcf_->exc_dn;
const double *const rh_up = xcf_->rho_up;
......@@ -323,13 +334,13 @@ void XCPotential::compute_stress(valarray<double>& sigma_exc)
dxc_ = 0.0;
const double *const e = xcf_->exc;
const int size = xcf_->np();
if ( nspin_ == 1 )
{
// unpolarized
const double *const rh = xcf_->rho;
const double *const v = xcf_->vxc1;
const int size = xcf_->np();
for ( int i = 0; i < size; i++ )
{
dxc_ += rh[i] * (e[i] - v[i]);
......@@ -342,7 +353,6 @@ void XCPotential::compute_stress(valarray<double>& sigma_exc)
const double *const rh_dn = xcf_->rho_dn;
const double *const v_up = xcf_->vxc1_up;
const double *const v_dn = xcf_->vxc1_dn;
const int size = xcf_->np();
for ( int i = 0; i < size; i++ )
{
const double rh = rh_up[i] + rh_dn[i];
......
......@@ -51,6 +51,7 @@ class XCPotential
public:
const XCFunctional* xcf() { return xcf_; }
bool isGGA(void);
XCPotential(const ChargeDensity& cd, const std::string functional_name);
~XCPotential();
void update(std::vector<std::vector<double> >& vr);
......
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