Commit 05c64241 authored by Francois Gygi's avatar Francois Gygi
Browse files

Implementation of spin-polarized BLYP functional


git-svn-id: http://qboxcode.org/svn/qb/trunk@1835 cba15fb0-1239-40c8-b417-11db7ca47a34
parent 3abe34cb
......@@ -15,7 +15,6 @@
// BLYPFunctional.C
//
////////////////////////////////////////////////////////////////////////////////
// $Id: BLYPFunctional.C,v 1.6 2008-09-08 15:56:18 fgygi Exp $
#include <cmath>
#include <cassert>
......@@ -47,8 +46,37 @@ BLYPFunctional::BLYPFunctional(const vector<vector<double> > &rhoe)
}
else
{
// not implemented
assert(false);
_exc_up.resize(_np);
_exc_dn.resize(_np);
_vxc1_up.resize(_np);
_vxc1_dn.resize(_np);
_vxc2_upup.resize(_np);
_vxc2_updn.resize(_np);
_vxc2_dnup.resize(_np);
_vxc2_dndn.resize(_np);
_grad_rho_up[0].resize(_np);
_grad_rho_up[1].resize(_np);
_grad_rho_up[2].resize(_np);
_grad_rho_dn[0].resize(_np);
_grad_rho_dn[1].resize(_np);
_grad_rho_dn[2].resize(_np);
rho_up = &rhoe[0][0];
rho_dn = &rhoe[1][0];
grad_rho_up[0] = &_grad_rho_up[0][0];
grad_rho_up[1] = &_grad_rho_up[1][0];
grad_rho_up[2] = &_grad_rho_up[2][0];
grad_rho_dn[0] = &_grad_rho_dn[0][0];
grad_rho_dn[1] = &_grad_rho_dn[1][0];
grad_rho_dn[2] = &_grad_rho_dn[2][0];
exc_up = &_exc_up[0];
exc_dn = &_exc_dn[0];
vxc1_up = &_vxc1_up[0];
vxc1_dn = &_vxc1_dn[0];
vxc2_upup = &_vxc2_upup[0];
vxc2_updn = &_vxc2_updn[0];
vxc2_dnup = &_vxc2_dnup[0];
vxc2_dndn = &_vxc2_dndn[0];
}
}
......@@ -73,7 +101,6 @@ void BLYPFunctional::setxc(void)
}
else
{
#if 0 // not implemented
assert( rho_up != 0 );
assert( rho_dn != 0 );
assert( grad_rho_up[0] != 0 && grad_rho_up[1] != 0 && grad_rho_up[2] != 0 );
......@@ -98,29 +125,29 @@ void BLYPFunctional::setxc(void)
double grx = grx_up + grx_dn;
double gry = gry_up + gry_dn;
double grz = grz_up + grz_dn;
double grad_up = sqrt(grx_up*grx_up + gry_up*gry_up + grz_up*grz_up);
double grad_dn = sqrt(grx_dn*grx_dn + gry_dn*gry_dn + grz_dn*grz_dn);
double grad = sqrt(grx*grx + gry*gry + grz*grz);
excpbe_sp(rho_up[i],rho_dn[i],grad_up,grad_dn,grad,&exc_up[i],&exc_dn[i],
&vxc1_up[i],&vxc1_dn[i],&vxc2_upup[i],&vxc2_dndn[i],
&vxc2_updn[i], &vxc2_dnup[i]);
double grad_up2 = grx_up*grx_up + gry_up*gry_up + grz_up*grz_up;
double grad_dn2 = grx_dn*grx_dn + gry_dn*gry_dn + grz_dn*grz_dn;
double grad_up_grad_dn = grx_up*grx_dn + gry_up*gry_dn + grz_up*grz_dn;
excblyp_sp(rho_up[i],rho_dn[i],grad_up2,grad_dn2,grad_up_grad_dn,
&exc_up[i],&exc_dn[i],&vxc1_up[i],&vxc1_dn[i],
&vxc2_upup[i],&vxc2_dndn[i],&vxc2_updn[i], &vxc2_dnup[i]);
}
#endif
}
}
////////////////////////////////////////////////////////////////////////////////
void BLYPFunctional::excblyp(double rho, double grad,
double *exc, double *vxc1, double *vxc2)
{
/* Becke exchange constants */
// Becke exchange constants
const double third = 1.0 / 3.0;
const double fourthirds = 4.0 / 3.0;
const double fivethirds = 5.0 / 3.0;
const double beta=0.0042;
const double ax = -0.7385587663820224058; /* -0.75*pow(3.0/pi,third) */
const double axa = -0.9305257363490999; /* -1.5*pow(3.0/(4*pi),third) */
const double ax = -0.7385587663820224058; // -0.75*pow(3.0/pi,third)
const double axa = -0.9305257363490999; // -1.5*pow(3.0/(4*pi),third)
/* LYP constants */
// LYP constants
const double a = 0.04918;
const double b = 0.132;
const double ab36 = a * b / 36.0;
......@@ -128,22 +155,17 @@ void BLYPFunctional::excblyp(double rho, double grad,
const double c_third = c / 3.0;
const double d = 0.349;
const double d_third = d / 3.0;
const double cf = 2.87123400018819; /* (3/10)*pow(3*pi*pi,2/3) */
const double cf = 2.87123400018819; // (3/10)*pow(3*pi*pi,2/3)
const double cfb = cf * b;
*exc = 0.0;
*vxc1 = 0.0;
*vxc2 = 0.0;
if ( rho < 1.e-18 )
{
return;
}
if ( rho < 1.e-18 ) return;
/*
* Becke's exchange
* A.D.Becke, Phys.Rev. B38, 3098 (1988)
*/
// Becke's exchange
// A.D.Becke, Phys.Rev. B38, 3098 (1988)
const double rha = 0.5 * rho;
const double grada = 0.5 * grad;
......@@ -155,19 +177,18 @@ void BLYPFunctional::excblyp(double rho, double grad,
const double asinhxa = asinh(xa);
const double frac = 1.0 / ( 1.0 + 6.0 * beta * xa * asinhxa );
const double ga = axa - beta * xa2 * frac;
/* N.B. in next line, ex is the energy density, hence rh13 */
// in next line, ex is the energy density, hence rh13
const double ex = rha13 * ga;
/* energy done, now the potential */
// potential
const double gpa = ( 6.0*beta*beta*xa2 * ( xa/sqrt(xa2+1.0) - asinhxa )
- 2.0*beta*xa ) * frac*frac;
const double vx1 = rha13 * fourthirds * ( ga - xa * gpa );
const double vx2 = - 0.5 * gpa / grada;
/*------------------------------------------------------------*/
/* LYP correlation */
/* Phys. Rev. B 37, 785 (1988). */
/* next lines specialized to the unpolarized case */
// LYP correlation
// Phys. Rev. B 37, 785 (1988).
// next lines specialized to the unpolarized case
const double rh13 = pow ( rho, third );
const double rhm13 = 1.0 / rh13;
const double rhm43 = rhm13 / rho;
......@@ -178,6 +199,8 @@ void BLYPFunctional::excblyp(double rho, double grad,
const double cfrac = num * deninv;
const double delta = rhm13 * ( c + d * deninv );
const double ddelta = - third * ( c * rhm43
+ d * rhm13 * rhm13 / ((d+rh13)*(d+rh13)) );
const double rhm53 = rhm43 * rhm13;
const double t1 = e * deninv;
const double t2 = rhm53;
......@@ -185,17 +208,15 @@ void BLYPFunctional::excblyp(double rho, double grad,
const double g = ab36 * t1 * t2 * t3;
/* next line, ec is the energy density, hence divide the energy by rho */
// ec is the energy density, hence divide the energy by rho
const double ec = - a * cfrac + 0.25 * g * grad * grad / rho;
/* energy done, now the potential */
// potential
const double de = c_third * rhm43 * e;
const double dnum = cfb * de;
const double dden = - d_third * rhm43;
const double dfrac = ( dnum * den - dden * num ) * deninv * deninv;
const double ddelta = - third * rhm43 * ( c + d * deninv ) -
rhm13 * d * dden * deninv * deninv;
const double dt1 = de * deninv - e * dden * deninv * deninv;
const double dt2 = - fivethirds * rhm53/rho;
const double dt3 = 14.0 * ddelta;
......@@ -209,3 +230,191 @@ void BLYPFunctional::excblyp(double rho, double grad,
*vxc1 = vx1 + vc1;
*vxc2 = vx2 + vc2;
}
////////////////////////////////////////////////////////////////////////////////
void BLYPFunctional::excblyp_sp(double rho_up, double rho_dn,
double grad_up2, double grad_dn2, double grad_up_grad_dn,
double *exc_up, double *exc_dn, double *vxc1_up, double *vxc1_dn,
double *vxc2_upup, double *vxc2_dndn, double *vxc2_updn, double *vxc2_dnup)
{
*exc_up = 0.0;
*exc_dn = 0.0;
*vxc1_up = 0.0;
*vxc1_dn = 0.0;
*vxc2_upup = 0.0;
*vxc2_updn = 0.0;
*vxc2_dnup = 0.0;
*vxc2_dndn = 0.0;
if ( rho_up < 1.e-18 && rho_dn < 1.e-18 ) return;
// Becke exchange constants
const double third = 1.0 / 3.0;
const double fourthirds = 4.0 / 3.0;
const double fivethirds = 5.0 / 3.0;
const double beta = 0.0042;
const double ax = -0.7385587663820224058; // -0.75*pow(3.0/pi,third)
const double cx = -0.9305257363490999; // -1.5*pow(3.0/(4*pi),third)
// LYP constants
const double a = 0.04918;
const double b = 0.132;
const double c = 0.2533;
const double d = 0.349;
const double cf = 2.87123400018819; // (3/10)*pow(3*pi*pi,2/3)
const double ninth = 1.0 / 9.0;
const double two113 = pow(2.0,11.0/3.0);
// Becke's exchange
// A.D.Becke, Phys.Rev. B38, 3098 (1988)
const double& rha = rho_up;
const double& rhb = rho_dn;
const double grada = sqrt(grad_up2);
const double gradb = sqrt(grad_dn2);
const double rha13 = pow ( rha, third );
const double rhb13 = pow ( rhb, third );
const double rha43 = rha * rha13;
const double rhb43 = rhb * rhb13;
const double rha83 = rha43 * rha43;
const double rhb83 = rhb43 * rhb43;
const double rha113 = rha83 * rha;
const double rhb113 = rhb83 * rhb;
double ex1a = 0.0;
double ex1b = 0.0;
double vx1a = 0.0;
double vx1b = 0.0;
double vx2a = 0.0;
double vx2b = 0.0;
if ( rho_up > 1.e-18 )
{
const double xa = grada / rha43;
const double xa2 = xa*xa;
const double asinhxa = asinh(xa);
const double fraca = 1.0 / ( 1.0 + 6.0 * beta * xa * asinhxa );
const double ga = cx - beta * xa2 * fraca;
// next line, ex is the energy density, hence rh13
ex1a = rha13 * ga;
const double gpa = ( 6.0*beta*beta*xa2 * ( xa/sqrt(xa2+1.0) - asinhxa )
- 2.0*beta*xa ) * fraca * fraca;
vx1a = rha13 * fourthirds * ( ga - xa * gpa );
vx2a = - gpa / grada;
}
if ( rho_dn > 1.e-18 )
{
const double xb = gradb / rhb43;
const double xb2 = xb*xb;
const double asinhxb = asinh(xb);
const double fracb = 1.0 / ( 1.0 + 6.0 * beta * xb * asinhxb );
const double gb = cx - beta * xb2 * fracb;
// next line, ex is the energy density, hence rh13
ex1b = rhb13 * gb;
const double gpb = ( 6.0*beta*beta*xb2 * ( xb/sqrt(xb2+1.0) - asinhxb )
- 2.0*beta*xb ) * fracb * fracb;
vx1b = rhb13 * fourthirds * ( gb - xb * gpb );
vx2b = - gpb / gradb;
}
// LYP correlation
// Phys. Rev. B 37, 785 (1988).
const double rho = rha + rhb;
const double rhoinv = 1.0/rho;
const double rhab = rha * rhb;
const double rhabrhm = rhab * rhoinv;
const double rh13 = pow ( rho, third );
const double rhm13 = 1.0 / rh13;
const double rhm43 = rhm13 / rho;
const double rhm113 = rhm43 * rhm43 * rhm43 * rh13;
const double rhm143 = rhm113 * rhoinv;
const double e = exp ( - c * rhm13 );
const double den = 1.0 + d * rhm13;
const double deninv = 1.0 / den;
const double w = e * deninv * rhm113;
const double dw = third*(c*d+(c-10.0*d)*rh13-11.0*rh13*rh13)*e*rhm143/
((d+rh13)*(d+rh13));
const double abw = a * b * w;
const double f1 = -4.0 * a * deninv * rhabrhm;
const double f2 = - two113 * cf * abw * rha * rhb * ( rha83 + rhb83 );
const double f = f1 + f2;
const double delta = rhm13 * ( c + d * deninv );
const double ddelta = - third * ( c * rhm43
+ d * rhm13 * rhm13 / ((d+rh13)*(d+rh13)) );
const double taa = 1.0 - 3.0 * delta + ( 11.0 - delta ) * rha * rhoinv;
const double tbb = 1.0 - 3.0 * delta + ( 11.0 - delta ) * rhb * rhoinv;
const double dtaa_drha = -3.0*ddelta - ddelta*rha*rhoinv +
(11.0-delta)*rhoinv*(1.0-rha*rhoinv);
const double dtaa_drhb = -3.0*ddelta - ddelta*rha*rhoinv -
(11.0-delta)*rhoinv*rha*rhoinv;
const double dtbb_drha = -3.0*ddelta - ddelta*rhb*rhoinv -
(11.0-delta)*rhoinv*rhb*rhoinv;
const double dtbb_drhb = -3.0*ddelta - ddelta*rhb*rhoinv +
(11.0-delta)*rhoinv*(1.0-rhb*rhoinv);
const double gaa = - abw * ( rhab * ninth * taa - rhb * rhb );
const double gbb = - abw * ( rhab * ninth * tbb - rha * rha );
const double gab = - abw * ( rhab * ninth * ( 47.0 - 7.0 * delta )
- fourthirds * rho * rho );
// next line, ec is the energy density, hence divide the energy by rho
const double ec1a = ( f1 + f2 + gaa * grad_up2
+ gab * grad_up_grad_dn
+ gbb * grad_dn2 ) / rho;
const double ec1b = ( f1 + f2 + gaa * grad_up2
+ gab * grad_up_grad_dn
+ gbb * grad_dn2 ) / rho;
const double A = -two113*cf*a*b;
const double df1_drha = -fourthirds*a*d*deninv*deninv*rhm43*rhabrhm
- 4.0*a*deninv*rhoinv*(rhb-rhabrhm);
const double df2_drha = A*dw*(rha113*rhb+rha*rhb113)
+ A*w*((11.0/3.0)*rha83*rhb+rhb113);
const double df1_drhb = -fourthirds*a*d*deninv*deninv*rhm43*rhabrhm
- 4.0*a*deninv*rhoinv*(rha-rhabrhm);
const double df2_drhb = A*dw*(rha113*rhb+rha*rhb113)
+ A*w*(rha113+rha*(11.0/3.0)*rhb83);
const double dgaa_drha = -a*b*dw*(rhab*ninth*taa-rhb*rhb)
- abw*(rhb*ninth*taa+rhab*ninth*dtaa_drha);
const double dgaa_drhb = -a*b*dw*(rhab*ninth*taa-rhb*rhb)
- abw*(rha*ninth*taa+rhab*ninth*dtaa_drhb-2.0*rhb);
const double dgbb_drha = -a*b*dw*(rhab*ninth*tbb-rha*rha)
- abw*(rhb*ninth*tbb+rhab*ninth*dtbb_drha-2.0*rha);
const double dgbb_drhb = -a*b*dw*(rhab*ninth*tbb-rha*rha)
- abw*(rha*ninth*tbb+rhab*ninth*dtbb_drhb);
const double dgab_drha = -a*b*dw*(rhab*ninth*(47.0-7.0*delta)
-fourthirds*rho*rho)
-abw*(rhb*ninth*(47.0-7.0*delta)
-rhab*ninth*7.0*ddelta-2.0*fourthirds*rho);
const double dgab_drhb = -a*b*dw*(rhab*ninth*(47.0-7.0*delta)
-fourthirds*rho*rho)
-abw*(rha*ninth*(47.0-7.0*delta)
-rhab*ninth*7.0*ddelta-2.0*fourthirds*rho);
const double vc1a = df1_drha + df2_drha
+ dgaa_drha * grad_up2
+ dgab_drha * grad_up_grad_dn
+ dgbb_drha * grad_dn2;
const double vc1b = df1_drhb + df2_drhb
+ dgaa_drhb * grad_up2
+ dgab_drhb * grad_up_grad_dn
+ dgbb_drhb * grad_dn2;
*exc_up = ex1a + ec1a;
*exc_dn = ex1b + ec1b;
*vxc1_up = vx1a + vc1a;
*vxc1_dn = vx1b + vc1b;
*vxc2_upup = vx2a - 2.0 * gaa;
*vxc2_updn = - gab;
*vxc2_dnup = - gab;
*vxc2_dndn = vx2b - 2.0 * gbb;
}
......@@ -15,8 +15,6 @@
// BLYPFunctional.h
//
////////////////////////////////////////////////////////////////////////////////
// $Id: BLYPFunctional.h,v 1.7 2009-06-29 09:57:57 fgygi Exp $
#ifndef BLYPFUNCTIONAL_H
#define BLYPFUNCTIONAL_H
......@@ -36,10 +34,9 @@ class BLYPFunctional : public XCFunctional
double *exc, double *vxc1, double *vxc2);
void excblyp_sp(double rho_up, double rho_dn,
double grad_up, double grad_dn, double grad,
double *exc_up, double *exc_dn,
double *vxc1_up, double *vxc1_dn, double *vxc2_upup, double *vxc2_dndn,
double *vxc2_updn, double *vxc2_dnup);
double grad_up2, double grad_dn2, double grad_up_grad_dn,
double *exc_up, double *exc_dn, double *vxc1_up, double *vxc1_dn,
double *vxc2_upup, double *vxc2_dndn, double *vxc2_updn, double *vxc2_dnup);
public:
......
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