Commit 2d6c6a91 by Francois Gygi

Add VWN functional


git-svn-id: http://qboxcode.org/svn/qb/trunk@1693 cba15fb0-1239-40c8-b417-11db7ca47a34
parent f1a28019
......@@ -28,7 +28,7 @@ OBJECTS=qb.o AtomSet.o Atom.o Species.o \
LoadCmd.o SaveCmd.o \
SpeciesCmd.o SpeciesReader.o SpeciesHandler.o \
XCOperator.o ExchangeOperator.o Bisection.o KPGridConnectivity.o \
XCPotential.o LDAFunctional.o \
XCPotential.o LDAFunctional.o VWNFunctional.o \
PBEFunctional.o BLYPFunctional.o B3LYPFunctional.o \
NonLocalPotential.o SampleReader.o StructuredDocumentHandler.o \
SampleHandler.o AtomSetHandler.o WavefunctionHandler.o \
......
////////////////////////////////////////////////////////////////////////////////
//
// Copyright (c) 2015 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/>.
//
////////////////////////////////////////////////////////////////////////////////
//
// VWNFunctional.C
//
// LDA Exchange-correlation energy and potential
// S.H.Vosko, L.Wilk, M.Nusair, Can. J. Phys. 58, 1200 (1980)
//
////////////////////////////////////////////////////////////////////////////////
#include <cmath>
#include <cassert>
#include <vector>
#include "VWNFunctional.h"
#include<iostream>
using namespace std;
void VWNFunctional::setxc(void)
{
if ( _np == 0 ) return;
if ( _nspin == 1 )
{
// unpolarized
assert(rho != 0);
assert(exc != 0);
assert(vxc1 != 0);
#pragma omp parallel for
for ( int ir = 0; ir < _np; ir++ )
{
double ex,vx,ec,vc;
xc_unpolarized(rho[ir],ex,vx,ec,vc);
exc[ir] = ex + ec;
vxc1[ir] = vx + vc;
}
}
else
{
// polarized
assert(rho_up != 0);
assert(rho_dn != 0);
assert(exc != 0);
assert(vxc1_up != 0);
assert(vxc1_dn != 0);
const double fz_prefac = 1.0 / ( cbrt(2.0)*2.0 - 2.0 );
const double dfz_prefac = (4.0/3.0) * fz_prefac;
#pragma omp parallel for
for ( int ir = 0; ir < _np; ir++ )
{
double ex = 0.0;
double ec = 0.0;
double vx_up = 0.0;
double vx_dn = 0.0;
double vc_up = 0.0;
double vc_dn = 0.0;
double roe_up = rho_up[ir];
double roe_dn = rho_dn[ir];
double roe = roe_up + roe_dn;
if ( roe > 0.0 )
{
const double zeta = ( roe_up - roe_dn ) / roe;
const double zp1 = 1.0 + zeta;
const double zm1 = 1.0 - zeta;
const double zp1_13 = cbrt(zp1);
const double zm1_13 = cbrt(zm1);
const double fz = fz_prefac * ( zp1_13 * zp1 + zm1_13 * zm1 - 2.0 );
const double dfz = dfz_prefac * ( zp1_13 - zm1_13 );
double ex_u,ec_u,vx_u,vc_u;
double ex_p,ec_p,vx_p,vc_p;
double a,da;
xc_unpolarized(roe,ex_u,vx_u,ec_u,vc_u);
xc_polarized(roe,ex_p,vx_p,ec_p,vc_p);
alpha_c(roe,a,da);
const double ex_pu = ex_p - ex_u;
ex = ex_u + fz * ex_pu;
double vx = vx_u + fz * ( vx_p - vx_u );
vx_up = vx + ex_pu * ( 1.0 - zeta ) * dfz;
vx_dn = vx - ex_pu * ( 1.0 + zeta ) * dfz;
const double zeta3 = zeta*zeta*zeta;
const double zeta4 = zeta3*zeta;
a *= (9.0/8.0)*fz_prefac;
da *= (9.0/8.0)*fz_prefac;
double ec_pu = ec_p - ec_u - a;
ec = ec_u + a * fz + ec_pu * fz * zeta4;
const double vc1 = vc_u + da * fz + ( vc_p - vc_u - da ) * fz * zeta4;
const double vc2 = a * dfz + ec_pu * zeta3 * ( 4.0*fz + zeta*dfz );
vc_up = vc1 + ( 1.0 - zeta ) * vc2;
vc_dn = vc1 - ( 1.0 + zeta ) * vc2;
}
vxc1_up[ir] = vx_up + vc_up;
vxc1_dn[ir] = vx_dn + vc_dn;
exc[ir] = ex + ec;
}
}
}
void VWNFunctional::xc_unpolarized(const double rh, double &ex, double &vx,
double &ec, double &vc)
{
// unpolarized xc energy and potential
// const double third=1.0/3.0;
// c1 is (3.D0/(4.D0*pi))**third
const double c1 = 0.6203504908994001;
// alpha = (4/(9*pi))**third = 0.521061761198
// const double alpha = 0.521061761198;
// c2 = -(3/(4*pi)) / alpha = -0.458165293283
// const double c2 = -0.458165293283;
// c3 = (4/3) * c2 = -0.610887057711
const double c3 = -0.610887057711;
ex = 0.0;
vx = 0.0;
ec = 0.0;
vc = 0.0;
if ( rh > 0.0 )
{
const double A = 0.0310907;
const double x0 = -0.10498;
const double b = 3.72744;
const double c = 12.9352;
const double Q = sqrt( 4.0 * c - b * b );
const double fac1 = 2.0 * b / Q;
const double fac2 = b * x0 / ( x0 * x0 + b * x0 + c );
const double fac3 = 2.0 * ( 2.0 * x0 + b ) / Q;
double ro13 = cbrt(rh);
double rs = c1 / ro13;
// Next line : exchange part in Hartree units
vx = c3 / rs;
ex = 0.75 * vx;
double sqrtrs = sqrt(rs);
double X = rs + b * sqrtrs + c;
double fatan = atan( Q / ( 2.0 * sqrtrs + b ) );
ec = A * ( log( rs / X ) + fac1 * fatan -
fac2 * ( log( (sqrtrs-x0)*(sqrtrs-x0) / X ) +
fac3 * fatan ));
double t = sqrtrs - x0;
vc = ec + ( A / 3.0 ) * ( b * sqrtrs * x0 - c * t ) / ( X * t );
}
}
void VWNFunctional::xc_polarized(const double rh, double &ex, double &vx,
double &ec, double &vc)
{
// polarized xc energy and potential
// const double third=1.0/3.0;
// c1 is (3.D0/(4.D0*pi))**third
const double c1 = 0.6203504908994001;
// alpha = (4/(9*pi))**third = 0.521061761198
// const double alpha = 0.521061761198;
// c2 = -(3/(4*pi)) / alpha = -0.458165293283
// const double c2 = -0.458165293283;
// c3 = (4/3) * c2 = -0.610887057711
// c4 = 2**third * c3
const double c4 = -0.769669463118;
ex = 0.0;
vx = 0.0;
ec = 0.0;
vc = 0.0;
if ( rh > 0.0 )
{
const double A = 0.01554535;
const double x0 = -0.32500;
const double b = 7.06042;
const double c = 18.0578;
const double Q = sqrt( 4.0 * c - b * b );
const double fac1 = 2.0 * b / Q;
const double fac2 = b * x0 / ( x0 * x0 + b * x0 + c );
const double fac3 = 2.0 * ( 2.0 * x0 + b ) / Q;
double ro13 = cbrt(rh);
double rs = c1 / ro13;
// Next line : exchange part in Hartree units
vx = c4 / rs;
ex = 0.75 * vx;
double sqrtrs = sqrt(rs);
double X = rs + b * sqrtrs + c;
double fatan = atan( Q / ( 2.0 * sqrtrs + b ) );
ec = A * ( log( rs / X ) + fac1 * fatan -
fac2 * ( log( (sqrtrs-x0)*(sqrtrs-x0) / X ) +
fac3 * fatan ));
double t = sqrtrs - x0;
vc = ec + ( A / 3.0 ) * ( b * sqrtrs * x0 - c * t ) / ( X * t );
}
}
void VWNFunctional::alpha_c(const double rh, double &a, double &da)
{
// VWN spin stiffness alpha_c(rh)
// a = spin stiffness
// da = d(rh * a)/drh
// const double third=1.0/3.0;
// c1 is (3.D0/(4.D0*pi))**third
const double c1 = 0.6203504908994001;
// alpha = (4/(9*pi))**third = 0.521061761198
// const double alpha = 0.521061761198;
// c2 = -(3/(4*pi)) / alpha = -0.458165293283
// const double c2 = -0.458165293283;
// c3 = (4/3) * c2 = -0.610887057711
a = 0.0;
da = 0.0;
if ( rh > 0.0 )
{
const double A = 1.0/(6.0*M_PI*M_PI);
const double x0 = -0.0047584;
const double b = 1.13107;
const double c = 13.0045;
const double Q = sqrt( 4.0 * c - b * b );
const double fac1 = 2.0 * b / Q;
const double fac2 = b * x0 / ( x0 * x0 + b * x0 + c );
const double fac3 = 2.0 * ( 2.0 * x0 + b ) / Q;
double ro13 = cbrt(rh);
double rs = c1 / ro13;
double sqrtrs = sqrt(rs);
double X = rs + b * sqrtrs + c;
double fatan = atan( Q / ( 2.0 * sqrtrs + b ) );
a = A * ( log( rs / X ) + fac1 * fatan -
fac2 * ( log( (sqrtrs-x0)*(sqrtrs-x0) / X ) +
fac3 * fatan ));
double t = sqrtrs - x0;
da = a + ( A / 3.0 ) * ( b * sqrtrs * x0 - c * t ) / ( X * t );
}
}
////////////////////////////////////////////////////////////////////////////////
//
// Copyright (c) 2015 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/>.
//
////////////////////////////////////////////////////////////////////////////////
//
// VWNFunctional.h
//
////////////////////////////////////////////////////////////////////////////////
#ifndef VWNFUNCTIONAL_H
#define VWNFUNCTIONAL_H
#include <vector>
#include <cassert>
#include "XCFunctional.h"
class VWNFunctional : public XCFunctional
{
void xc_unpolarized(const double rh, double &ex, double &vx,
double &ec, double &vc);
void xc_polarized(const double rh, double &ex, double &vx,
double &ec, double &vc);
void alpha_c(const double rh, double &a, double &da);
std::vector<double> _exc;
std::vector<std::vector<double> > _vxc;
VWNFunctional();
public:
VWNFunctional(const std::vector<std::vector<double> > &rhoe)
{
_nspin = rhoe.size();
if ( _nspin > 1 ) assert(rhoe[0].size() == rhoe[1].size());
_np = rhoe[0].size();
_exc.resize(_np);
_vxc.resize(_nspin);
for ( int i = 0; i < _nspin; i++ )
{
_vxc[i].resize(_np);
}
if ( _nspin == 1 )
{
rho = &rhoe[0][0];
exc = &_exc[0];
vxc1 = &_vxc[0][0];
}
else
{
rho_up = &rhoe[0][0];
rho_dn = &rhoe[1][0];
exc = &_exc[0];
vxc1_up = &_vxc[0][0];
vxc1_dn = &_vxc[1][0];
}
};
bool isGGA() const { return false; };
std::string name() const { return "VWN"; };
void setxc(void);
};
#endif
......@@ -35,6 +35,7 @@ XCOperator::XCOperator(Sample& s, const ChargeDensity& cd) :cd_(cd)
// check the name of the functional
if ( ( functional_name == "LDA" ) ||
( functional_name == "VWN" ) ||
( functional_name == "PBE" ) ||
( functional_name == "BLYP" ) )
{
......
......@@ -18,6 +18,7 @@
#include "XCPotential.h"
#include "LDAFunctional.h"
#include "VWNFunctional.h"
#include "PBEFunctional.h"
#include "BLYPFunctional.h"
#include "B3LYPFunctional.h"
......@@ -35,6 +36,10 @@ XCPotential::XCPotential(const ChargeDensity& cd, const string functional_name,
{
xcf_ = new LDAFunctional(cd_.rhor);
}
else if ( functional_name == "VWN" )
{
xcf_ = new VWNFunctional(cd_.rhor);
}
else if ( functional_name == "PBE" )
{
xcf_ = new PBEFunctional(cd_.rhor);
......
......@@ -45,6 +45,7 @@ class Xc : public Var
string v = argv[1];
if ( !( v == "LDA" ||
v == "VWN" ||
v == "PBE" ||
v == "BLYP" ||
v == "HF" ||
......@@ -52,7 +53,7 @@ class Xc : public Var
v == "B3LYP" ) )
{
if ( ui->onpe0() )
cout << " xc must be LDA, PBE, BLYP, HF, PBE0 or B3LYP" << endl;
cout << " xc must be LDA, VWN, PBE, BLYP, HF, PBE0 or B3LYP" << 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