/*
DFT++ is a density functional package developed by the research group
of Professor Tomas Arias
Copyright 1996-2003 Sohrab Ismail-Beigi
This file is part of DFT++.
DFT++ is free software; you can redistribute it and/or modify
it 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.
DFT++ is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with DFT++; if not, write to the Free Software
Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
Please see the file CREDITS for a list of authors.
For academic users, we request that publications using results obtained with
this software reference
"New algebraic formulation of density functional calculation," by Sohrab Ismail-Beigi
and T.A. Arias, Computer Physics Communications 128:1-2, 1-45 (June 2000).
and, if using the wavelet basis, further reference
"Multiresolution analysis of electronic structure: semicardinal and wavelet bases,"
T.A. Arias, Reviews of Modern Physics 71:1, 267-311 (January 1999).
and
"Robust ab initio calculation of condensed matter: transparent convergence through
semicardinal multiresolution analysis,'' I.P. Daykov, T.A. Arias, and
Torkel D. Engeness, Physical Review Letters, 90:21, 216402 (May 2003).
For your convenience, preprints of the above articles may be obtained from
http://arXiv.org/abs/cond-mat/9909130, 9805262, and 0204411, respectively.
*/
/*
* Sohrab Ismail-Beigi Original code, Sept. 9, 1996
* Revised, Jan. 1997
*
* Calculate exchange-correlation energy and its derivative versus
* the electron (number) density.
*
*/
/* $Id: exc.cpp,v 1.9.2.15 2003/05/29 18:54:23 ivan Exp $ */
#include "header.h"
/*
real exc_PZ(const real n);
real excprime_PZ(const real n);
real exc_VWN(const real n);
real excprime_VWN(const real n);
*/
/*
* Return the exchange-correlation energy density for density n.
*
* Perdew-Zunger Approximation
*
* The formula for Exc for rs < 1 are the standard expansions in
* rs; that for rs > 1 was copied from a code by Fatih Yanik.
*
*/
real exc_PZ(const real &n)
{
/* constants */
static const real pi = M_PI,
X1 = 0.75*pow(3.0/(2.0*pi),2.0/3.0), /* Exchange energy coeff */
AA = 0.0622/2.0, /* Correl. energy rs < 1 expansion coeffs */
BB = -0.0960/2.0,
CC = 0.0040/2.0,
DD = -0.0232/2.0,
GA = -0.2846/2.0, /* Coerrel. energy rs > 1 expansion coeffs */
B1 = 1.0529,
B2 = 0.3334;
/* Temporary variables */
real rs;
real excPZ; //output
/* Calculate rs */
if (n <= 0.0)
// rs = 1.0e15;
return excPZ = 0.0;
else
rs = pow(4.0*pi/3.0*n,-1.0/3.0);
/* Calculate exc */
/* Exchange part */
excPZ = -X1/rs;
/* Correlation part */
if (rs > 1.0)
excPZ += GA/(1.0 + B1*sqrt(rs) + B2*rs );
else
excPZ += AA*log(rs) + BB + CC*rs*log(rs) + DD*rs;
return excPZ;
}
/*
* Derivative of above versus n
*/
real
excprime_PZ(const real &n)
{
/* constants */
static const real pi = M_PI,
X1 = 0.75*pow(3.0/(2.0*pi),2.0/3.0), /* Exchange energy coeff */
AA = 0.0622/2.0, /* Correl. energy rs < 1 expansion coeffs */
// BB = -0.0960/2.0,
CC = 0.0040/2.0,
DD = -0.0232/2.0,
GA = -0.2846/2.0, /* Coerrel. energy rs > 1 expansion coeffs */
B1 = 1.0529,
B2 = 0.3334;
/* temp. vars */
real rs;
real excprimePZ;
/* Calculate rs */
if (n <= 0.0)
//rs = 1.0e15;
return excprimePZ = 0.0;
else
rs = pow(4.0*pi/3.0*n,-1.0/3.0);
/* Calculate dexc/dn: drs/dn = -rs/(3*n) is handy */
/* Exchange part */
excprimePZ = -X1/(3.0*rs*n);
/* Correlation part */
if (rs > 1.0)
excprimePZ += GA*(0.5*B1*sqrt(rs)+B2*rs)/
(3.0*n*pow(1.0+B1*sqrt(rs)+B2*rs,2.0));
else
excprimePZ +=
-(AA + CC*rs*(1.0+log(rs)) + DD*rs)/(3.0*n);
return excprimePZ;
}
real
exc_VWN(const real &n)
{
/* constants */
static const real pi = M_PI,
X1 = 0.75*pow(3.0/(2.0*pi),2.0/3.0), /* Exchange energy coeff */
A = 0.0310907,
x0 = -0.10498,
b = 3.72744,
c = 12.9352,
Q = sqrt(4*c-b*b),
X0 = x0*x0+b*x0+c;
real rs,x,X;
if (n <= 0.0)
return 0.;
else{
rs = pow(4.0*pi/3.0*(n),-1.0/3.0);
x=sqrt(rs); X=x*x+b*x+c;
return -X1/rs /* Exchange part */
+ A*( /* Correlation part */
+log(x*x/X)+2*b/Q*atan(Q/(2*x+b))
-(b*x0)/X0*(
log((x-x0)*(x-x0)/X)+2*(2*x0+b)/Q*atan(Q/(2*x+b))
)
);
}
}
/* Correlation part ONLY */
real
ec_VWN(const real &n)
{
/* constants */
static const real pi = M_PI,
// X1 = 0.75*pow(3.0/(2.0*pi),2.0/3.0), /* Exchange energy coeff */
A = 0.0310907,
x0 = -0.10498,
b = 3.72744,
c = 12.9352,
Q = sqrt(4*c-b*b),
X0 = x0*x0+b*x0+c;
real rs,x,X;
if (n <= 0.0)
return 0.;
else{
rs = pow(4.0*pi/3.0*(n),-1.0/3.0);
x=sqrt(rs); X=x*x+b*x+c;
return A*( /* Correlation part */
+log(x*x/X)+2*b/Q*atan(Q/(2*x+b))
-(b*x0)/X0*(
log((x-x0)*(x-x0)/X)+2*(2*x0+b)/Q*atan(Q/(2*x+b))
)
);
}
}
real
excprime_VWN(const real &n)
{
/* constants */
static const real pi = M_PI,
X1 = 0.75*pow(3.0/(2.0*pi),2.0/3.0), /* Exchange energy coeff */
A = 0.0310907,
x0 = -0.10498,
b = 3.72744,
c = 12.9352,
Q = sqrt(4*c-b*b),
X0 = x0*x0+b*x0+c;
real rs,x,X,dx;
if (n <= 0.0)
return 0.;
else {
rs = pow(4.0*pi/3.0*(n),-1.0/3.0);
x=sqrt(rs); X=x*x+b*x+c; dx=0.5/x;
return -rs/(3*n)*dx* /* Factors from chain rule */
(
2*X1/(rs*x)+ /* Exchange part */
A*( /* Correlation part */
2./x-(2*x+b)/X-4*b/(Q*Q+(2*x+b)*(2*x+b))
-(b*x0)/X0*(2/(x-x0)-(2*x+b)/X-4*(2*x0+b)/
(Q*Q+(2*x+b)*(2*x+b)) )
)
);
}
}
/* Correleation part ONLY */
real
ecprime_VWN(const real &n)
{
/* constants */
static const real pi = M_PI,
// X1 = 0.75*pow(3.0/(2.0*pi),2.0/3.0), /* Exchange energy coeff */
A = 0.0310907,
x0 = -0.10498,
b = 3.72744,
c = 12.9352,
Q = sqrt(4*c-b*b),
X0 = x0*x0+b*x0+c;
real rs,x,X,dx;
if (n <= 0.0)
return 0.;
else {
rs = pow(4.0*pi/3.0*(n),-1.0/3.0);
x=sqrt(rs); X=x*x+b*x+c; dx=0.5/x;
return -rs/(3*n)*dx* /* Factors from chain rule */
(
A*( /* Correlation part */
2./x-(2*x+b)/X-4*b/(Q*Q+(2*x+b)*(2*x+b))
-(b*x0)/X0*(2/(x-x0)-(2*x+b)/X-4*(2*x0+b)/
(Q*Q+(2*x+b)*(2*x+b)) )
)
);
}
}
////////// The ones below loop over the data
/* TAA: handy 'square' function */
#define SQ(a) ((a)*(a))
/* TAA: parameters and constants */
#define DENS_CUT 1e-12 /* low-end density to avoid numerical instabilities */
double beta=0.066725, gam=(1-log(2))/SQ(M_PI); /* Corr constants */
double mu=beta*SQ(M_PI)/3.,kappa=0.804; /* Exchange constants */
double bog=beta/gam; /* Common combinations of constants */
/* TAA: exchange function */
double ex(double kF) {
return -3./(4.*M_PI)*kF;
}
/* TAA: non-local exhange-enhancement function */
double Fx(double s) {
return (1+mu*SQ(s)*(1+1./kappa))/(1+mu*SQ(s)/kappa);
}
/* TAA: derivative of non-local echange-enhancement function */
double Fxp(double s) {
return
2*s*mu/SQ(1+mu*SQ(s)/kappa);
}
/* TAA: PBE's magic non-local correlation contrib "H" */
double H_PBE(double A,double t)
{
double u=A*SQ(t);
return gam*log(1+bog*SQ(t)*(1+u)/(1+u+SQ(u)));
}
/* TAA: Partial derivs of H and internal variables */
double HA_PBE(double A,double t)
{
/* TAA: simplify factors later */
double u=A*SQ(t);
double denom=1./(1+u+SQ(u));
return beta/(1+bog*SQ(t)*(1+u)*denom)*
SQ(SQ(t))*(-u*(2+u))*SQ(denom);
}
double Ht_PBE(double A,double t)
{
/* TAA: simplify factors later */
double u=A*SQ(t);
double denom=1./(1+u+SQ(u));
return beta/(1+bog*SQ(t)*(1+u)*denom)*
2*t*(1+2*u)*SQ(denom);
}
void
exc_GGA(const RealSpaceScalarFieldColumn &n,
const RealSpaceScalarFieldColumn &Dx,
const RealSpaceScalarFieldColumn &Dy,
const RealSpaceScalarFieldColumn &Dz,
RealSpaceScalarFieldColumn &excn)
{
if(n.representation != REALSPACE)
die("Charge density is not stored in the realspace representation!\n");
if(n.length!=excn.length)
die("Incompatible lengths in exc_LDA()\n");
int i;
double dens,kF,s;
double ks,t,ec,A;
/* Calculate exc */
for (i=0; i < n.length; i++)
{
/* Preliminaries */
dens=n.data.d[i].x+DENS_CUT;
kF=pow(3*SQ(M_PI)*dens,1./3.);
s=DENS_CUT+sqrt(SQ(Dx.data.d[i].x)+SQ(Dy.data.d[i].x)+SQ(Dz.data.d[i].x))/
(2*kF*dens);
ks=sqrt(4.*kF/M_PI);
t=s*(kF/ks);
ec=ec_VWN(dens);
A=bog/(exp(-ec/gam)-1); /* TAA:
reexpress as exp/sinh for num'cs */
/* Exchange contribution */
excn.data.d[i]=ex(kF)*Fx(s);
/* Correlation contrib */
excn.data.d[i]+=ec+H_PBE(A,t);
}
excn.representation=REALSPACE;
}
void
exc_LDA(const RealSpaceScalarFieldColumn &n, RealSpaceScalarFieldColumn &excn)
{
if(n.representation != REALSPACE)
die("Charge density is not stored in the realspace representation!\n");
if(n.length!=excn.length)
die("Incompatible lengths in exc_LDA()\n");
int i;
/* Calculate exc */
for (i=0; i < n.length; i++) {
excn.data.d[i].x = exc_PZ(n.data.d[i].x);
excn.data.d[i].y = 0.0;
}
excn.representation=REALSPACE;
}
RealSpaceScalarFieldColumn
exc_GGA(const RealSpaceScalarFieldColumn &n,
const RealSpaceScalarFieldColumn &Dx,
const RealSpaceScalarFieldColumn &Dy,
const RealSpaceScalarFieldColumn &Dz
)
{
RealSpaceScalarFieldColumn excn(n);
exc_GGA(n,Dx,Dy,Dz,excn);
return excn;
}
RealSpaceScalarFieldColumn
exc_LDA(const RealSpaceScalarFieldColumn &n)
{
RealSpaceScalarFieldColumn excn(n);
exc_LDA(n,excn);
return excn;
}
void
Dnexc_GGA(const RealSpaceScalarFieldColumn &n,
const RealSpaceScalarFieldColumn &Dx,
const RealSpaceScalarFieldColumn &Dy,
const RealSpaceScalarFieldColumn &Dz,
RealSpaceScalarFieldColumn &excpn)
{
if(n.representation != REALSPACE)
die("Charge density is not stored in the realspace representation!\n");
if(n.length!=excpn.length)
die("Incompatible lengths in exc_LDA()\n");
int i;
/* Calculate exc */
double dens,kF,s;
double ks,t,ec,A;
for (i=0; i < n.length; i++)
{
/* Preliminaries */
dens=n.data.d[i].x+DENS_CUT;
kF=pow(3*SQ(M_PI)*dens,1./3.);
s=DENS_CUT+sqrt(SQ(Dx.data.d[i].x)+SQ(Dy.data.d[i].x)+SQ(Dz.data.d[i].x))/
(2*kF*dens);
ks=sqrt(4.*kF/M_PI);
t=s*(kF/ks);
ec=ec_VWN(dens);
A=bog/(exp(-ec/gam)-1); /* TAA:
reexpress as exp/sinh for num'cs */
/* Exchange contrib */
excpn.data.d[i] =
ex(kF)*( Fx(s)-4*Fxp(s)*s )/(3*dens);
/* Correlation contrib */
excpn.data.d[i]+=ecprime_VWN(dens)
/* dH/dA*(dA/dn) */
+HA_PBE(A,t)*SQ(A)*exp(-ec/gam)*ecprime_VWN(dens)/beta
/* dH/dt*(dt/dn) */
+Ht_PBE(A,t)*t*(-7./6.)/dens;
}
excpn.representation=REALSPACE;
}
void
DDexc_GGA(int dir,
const RealSpaceScalarFieldColumn &n,
const RealSpaceScalarFieldColumn &Dx,
const RealSpaceScalarFieldColumn &Dy,
const RealSpaceScalarFieldColumn &Dz,
RealSpaceScalarFieldColumn &excpn)
{
if(n.representation != REALSPACE)
die("Charge density is not stored in the realspace representation!\n");
if(n.length!=excpn.length)
die("Incompatible lengths in exc_LDA()\n");
/* Pick out component which is varying */
RealSpaceScalarFieldColumn D(Dx);
switch (dir)
{
case 0:
D=Dx;
break;
case 1:
D=Dy;
break;
case 2:
D=Dz;
break;
}
int i;
double dens,kF,s;
double ks,t,ec,A;
/* Calculate exc */
for (i=0; i < n.length; i++)
{
/* Preliminaries */
dens=n.data.d[i].x+DENS_CUT;
kF=pow(3*SQ(M_PI)*dens,1./3.);
s=DENS_CUT+sqrt(SQ(Dx.data.d[i].x)+SQ(Dy.data.d[i].x)+SQ(Dz.data.d[i].x))/
(2*kF*dens);
ks=sqrt(4.*kF/M_PI);
t=s*(kF/ks);
ec=ec_VWN(dens);
A=bog/(exp(-ec/gam)-1); /* TAA:
reexpress as exp/sinh for num'cs */
/* Exchange contrib */
excpn.data.d[i] = ex(kF)*Fxp(s)/s/SQ(2*kF*dens)*D.data.d[i];
/* Correlation contrib */
excpn.data.d[i] += Ht_PBE(A,t)/t/SQ(2*ks*dens)*D.data.d[i];
}
excpn.representation=REALSPACE;
}
void
excprime_LDA(const RealSpaceScalarFieldColumn &n,
RealSpaceScalarFieldColumn &excpn)
{
if(n.representation != REALSPACE)
die("Charge density is not stored in the realspace representation!\n");
if(n.length!=excpn.length)
die("Incompatible lengths in exc_LDA()\n");
int i;
/* Calculate exc */
for (i=0; i < n.length; i++) {
excpn.data.d[i].x = excprime_PZ(n.data.d[i].x);
excpn.data.d[i].y = 0.0;
}
excpn.representation=REALSPACE;
}
RealSpaceScalarFieldColumn
Dnexc_GGA(const RealSpaceScalarFieldColumn &n,
const RealSpaceScalarFieldColumn &Dx,
const RealSpaceScalarFieldColumn &Dy,
const RealSpaceScalarFieldColumn &Dz)
{
RealSpaceScalarFieldColumn excpn(n);
Dnexc_GGA(n,Dx,Dy,Dz,excpn);
return excpn;
}
RealSpaceScalarFieldColumn
DDexc_GGA(int dir,const RealSpaceScalarFieldColumn &n,
const RealSpaceScalarFieldColumn &Dx,
const RealSpaceScalarFieldColumn &Dy,
const RealSpaceScalarFieldColumn &Dz)
{
RealSpaceScalarFieldColumn excpn(n);
DDexc_GGA(dir,n,Dx,Dy,Dz,excpn);
return excpn;
}
RealSpaceScalarFieldColumn
excprime_LDA(const RealSpaceScalarFieldColumn &n)
{
RealSpaceScalarFieldColumn excpn(n);
excprime_LDA(n,excpn);
return excpn;
}
// Everything below is commented out for now! -- mhe, 2/11/02
/**********************************************************
* *
* New code for Generalized Gradient Correction *
* *
*********************************************************/
/* The following Generalized gradient correction
* is based on PW91
*
*/
/*
#include "parallel.h"
#define IMAG(z) ((z).y)
// cutoff charge density: to avoid true vacuum.
#define CHD_CUTOFF 1e-12
//
// input: d - density
// s - abs(grad d)/(2*kf*d)
// output: gradient correction to exchange energy per electron.
//
static real
exe(real d, real s)
{
// constants block
const real a1 = 0.19645;
const real a2 = 0.27430;
const real a3 = 0.15084;
const real a4 = 100.0;
const real ax = -0.7385588;
const real a = 7.7956;
const real b1 = 0.004;
const real thrd = 1.0/3.0;
// local variable block
real fac;
real s2, s3, s4;
real p0, p1, p2, p3, p4;
fac = ax*pow(d,thrd);
s2 = s*s;
s3 = s2*s;
s4 = s2*s2;
p0 = 1.0/sqrt(1.0+a*a*s2);
p1 = log(a*s+1.0/p0);
p2 = exp(-a4*s2);
p3 = 1.0/(1.0+a1*s*p1+b1*s4);
p4 = 1.0+a1*s*p1+(a2-a3*p2)*s2;
return (fac*(p3*p4-1.0));
}
*/
//
// input: d - density
// s - abs(grad d)/(2*kf*d)
// u - (grad d)*grad(abs(grad d))/(d^2*(2*kf)^3)
// v - lap(d)/(d*(2*kf)^2)
// output: gradient correction to exchange potential per electron.
//
/*
static real
exvx(real d, real s, real u, real v)
{
// constants block
const real a1 = 0.19645;
const real a2 = 0.27430;
const real a3 = 0.15084;
const real a4 = 100.0;
const real ax = -0.7385588;
const real a = 7.7956;
const real b1 = 0.004;
const real thrd = 1.0/3.0;
const real thrd4 = 4.0/3.0;
// local variable block
real fac;
real s2, s3, s4;
real p0, p1, p2, p3, p4, p5, p6, p7, p8, p9, p10, p11;
real f;
real fss;
real fs;
fac = ax*pow(d,thrd);
s2 = s*s;
s3 = s2*s;
s4 = s2*s2;
p0 = 1.0/sqrt(1.0+a*a*s2);
p1 = log(a*s+1.0/p0);
p2 = exp(-a4*s2);
p3 = 1.0/(1.0+a1*s*p1+b1*s4);
p4 = 1.0+a1*s*p1+(a2-a3*p2)*s2;
f = p3*p4;
p5 = b1*s2-(a2-a3*p2);
p6 = a1*s*(p1+a*s*p0);
p7 = 2.0*(a2-a3*p2)+2.0*a3*a4*s2*p2-4.0*b1*s2*f;
fs = p3*(p3*p5*p6+p7);
p8 = 2.0*s*(b1-a3*a4*p2);
p9 = a1*p1+a*a1*s*p0*(3.0-a*a*s2*p0*p0);
p10 = 4.0*a3*a4*s*p2*(2.0-a4*s2)-8.0*b1*s*f-4.0*b1*s3*fs;
p11 = -p3*p3*(a1*p1+a*a1*s*p0+4.0*b1*s3);
fss = p3*p3*(p5*p9+p6*p8)+2.0*p3*p5*p6*p11+p3*p10+p7*p11;
return fac*(thrd4*(f-1.0)-(u-thrd4*s3)*fss-v*fs);
}
*/
//
// input : a, a1, b1, b2, b3, b4, p, rs
// rs - seitz radius
// output: gg, ggrs (implicit!!)
//
/*
static void
gcor(real a, real a1, real b1, real b2, real b3, real b4, real p, real rs,
real &gg, real &ggrs)
{
// local variables
real p1;
real rs12, rs32, rsp;
real q0, q1, q2, q3;
p1 = p + 1.0;
q0 = -2.0*a*(1.0+a1*rs);
rs12 = sqrt(rs);
rs32 = rs12*rs12*rs12;
rsp = pow(rs,p);
q1 = 2.0*a*(b1*rs12+b2*rs+b3*rs32+b4*rs*rsp);
q2 = log(1.0+1.0/q1);
gg = q0*q2;
q3 = a*(b1/rs12+2.0*b2+3.0*b3*rs12+2.0*b4*p1*rsp);
ggrs = -2.0*a*a1*q2-q0*q3/(q1*q1+q1);
return;
}
*/
//
// input : rs - seitz radius
// zet - relative polarization (here assume zet = 0.0)
// output : ec - local correlation energy per electron
// vcup, vcdn - local up & down correlational potential
// ecrs - derivative of ec w.r.t rs
// eczet - derivative of ec w.r.t zet
// alfc - correlation contribution to the spin stiffness
//
/*
static void
corlsd(real rs, real zet, real &ec, real &vcup, real &vcdn,
real &ecrs, real &eczet, real &alfc)
{
// constants block
const real gam = 0.5198421;
const real fzz = 1.709921;
const real thrd = 1.0/3.0;
const real thrd4 = 4.0/3.0;
// local variables
real f;
real eu, eurs;
real ep, eprs;
real alfm, alfrsm;
real z4;
real fz;
real comm;
f = (pow(1.0+zet,thrd4)+pow(1.0-zet,thrd4)-2.0)/gam;
gcor(0.03109070,0.213700,7.59570,3.58760,1.63820,
0.492940,1.000,rs,eu,eurs);
gcor(0.015545350,0.205480,14.11890,6.19770,3.36620,
0.625170,1.000,rs,ep,eprs);
gcor(0.01688690,0.111250,10.3570,3.62310,0.880260,
0.496710,1.000,rs,alfm,alfrsm);
// alfm is minus the spin stiffness alfc
alfc = -alfm;
z4 = (zet*zet)*(zet*zet);
ec = eu*(1.0-f*z4)+ep*f*z4-alfm*f*(1.0-z4)/fzz;
// energy done. now the potential:
ecrs = eurs*(1.0-f*z4)+eprs*f*z4-alfrsm*f*(1.0-z4)/fzz;
fz = thrd4*(pow(1.0+zet,thrd)-pow(1.0-zet,thrd))/gam;
eczet = 4.0*(zet*zet*zet)*f*(ep-eu+alfm/fzz)+fz*(z4*ep-z4*eu
-(1.0-z4)*alfm/fzz);
comm = ec -rs*ecrs/3.0-zet*eczet;
vcup = comm + eczet;
vcdn = comm - eczet;
return;
}
*/
//
// input : rs - seitz radius
// zet - relative polarization
// t - abs(grad d)/(d*2*ks*g) (g=ga=1.0 here)
// g (=1.0)
// ec - local correlation energy
// sk - ks
// fk - kf
// output: non-local correlation energy per electron
//
/*
static real
core(real rs, real zet, real t, real g, real ec, real sk, real fk)
{
// constants block
const real xnu = 15.75592;
const real cc0 = 0.004235;
const real cx = -0.001667212;
const real alf = 0.09;
const real c1 = 0.002568;
const real c2 = 0.023266;
const real c3 = 7.389e-6;
const real c4 = 8.723;
const real c5 = 0.472;
const real c6 = 7.389e-2;
const real a4 = 100.0;
// const real thrm = -1.0/3.0;
// const real thrd2 = 2.0/3.0;
// local variables
real bet;
real delt;
real g3, g4;
real pon;
real b, b2;
real t2, t4, t6;
real rs2, rs3;
real q4, q5, q6, q7;
real cc;
real r0, r1, r2, r3;
real coeff;
real h0, h1;
bet = xnu*cc0;
delt = 2.0*alf/bet;
g3 = g*g*g;
g4 = g3*g;
pon = -delt*ec/(g3*bet);
b = delt/(exp(pon)-1.0);
b2 = b*b;
t2 = t*t;
t4 = t2*t2;
t6 = t4*t2;
rs2 = rs*rs;
rs3 = rs2*rs;
q4 = 1.0+b*t2;
q5 = 1.0+b*t2+b2*t4;
q6 = c1+c2*rs+c3*rs2;
q7 = 1.0+c4*rs+c5*rs2+c6*rs3;
cc = -cx + q6/q7;
r0 = (sk/fk);
r0 *= r0;
r1 = a4*r0*g4;
coeff = cc-cc0-3.0*cx/7.0;
r2 = xnu*coeff*g3;
r3 = exp(-r1*t2);
h0 = g3*(bet/delt)*log(1.0+delt*q4*t2/q5);
h1 = r3*r2*t2;
return (h0+h1);
}
*/
//
// input : rs - seitz radius
// zet - relative polarization
// t - abs(grad d)/(d*2*ks*g) (g=ga=1.0 here)
// uu - (grad d)*grad(abs(grad d))/(d^2*(2*ks*g)^3)
// vv - lap(d)/(d*(2*ks*g)^2)
// ww - (grad d)*(grad zet)/(d*(2*ks*g)^2) (=0.0 here)
// g (=1.0)
// ec - local correlation energy
// ecrs - derivative of ec w.r.t rs
// sk - ks
// fk - kf
// output: non-local correlation energy per electron
//
/*
static real
corvx(real rs, real zet, real t, real uu, real vv, real ww,
real g, real ec, real eczet, real ecrs, real sk, real fk)
{
// constants block
const real xnu = 15.75592;
const real cc0 = 0.004235;
const real cx = -0.001667212;
const real alf = 0.09;
const real c1 = 0.002568;
const real c2 = 0.023266;
const real c3 = 7.389e-6;
const real c4 = 8.723;
const real c5 = 0.472;
const real c6 = 7.389e-2;
const real a4 = 100.0;
const real thrdm = -1.0/3.0;
const real thrd2 = 2.0/3.0;
// local variables
real bet;
real delt;
real g3, g4;
real pon;
real b, b2;
real t2, t4, t6;
real rs2, rs3;
real q4, q5, q6, q7, q8, q9;
real cc;
real r0, r1, r2, r3, r4;
real coeff;
real h0, h1, h, h0b, h0rs, h0bt, h0rst, h0z, h0t, h0zt, h0tt, h1rs, h1rst;
// Some version of Unix predefines "hz" (!?!#&@) as some number
// of Hertz... this caused compilation problems on the SP2.
#ifdef hz
#undef hz
#endif
real h1z, h1t, h1zt, h1tt, hrs, hrst, ht, htt, hz, hzt;
real ccrs;
real rsthrd;
real gz, fac;
real bg, bec;
real fact0, fact1, fact2, fact3, fact4, fact5;
real comm, pref;
bet = xnu*cc0;
delt = 2.0*alf/bet;
g3 = g*g*g;
g4 = g3*g;
pon = -delt*ec/(g3*bet);
b = delt/(exp(pon)-1.0);
b2 = b*b;
t2 = t*t;
t4 = t2*t2;
t6 = t4*t2;
rs2 = rs*rs;
rs3 = rs2*rs;
q4 = 1.0+b*t2;
q5 = 1.0+b*t2+b2*t4;
q6 = c1+c2*rs+c3*rs2;
q7 = 1.0+c4*rs+c5*rs2+c6*rs3;
cc = -cx + q6/q7;
r0 = (sk/fk);
r0 *= r0;
r1 = a4*r0*g4;
coeff = cc-cc0-3.0*cx/7.0;
r2 = xnu*coeff*g3;
r3 = exp(-r1*t2);
h0 = g3*(bet/delt)*log(1.0+delt*q4*t2/q5);
h1 = r3*r2*t2;
h = h0 + h1;
ccrs = (c2+2.0*c3*rs)/q7 - q6*(c4+2.0*c5*rs+3.0*c6*rs2)/(q7*q7);
rsthrd = rs/3.0;
r4 = rsthrd*ccrs/coeff;
// if(abs(zet) .ge. 1.d0) then
// if(zet .lt. 0.d0) zet=-1.d0+1.d-15
// if(zet .gt. 0.d0) zet=1.d0-1.d-15
// if(flag .eqv. .true.) write(*,50)
// flag=.false.
//50 format(/'warning: corrga - zet substituted')
// endif
gz = (pow(1.0+zet,thrdm) - pow(1.0-zet,thrdm))/3.0;
fac = delt/b+1.0;
bg = -3.0*b2*ec*fac/(bet*g4);
bec = b2*fac/(bet*g3);
q8 = q5*q5+delt*q4*q5*t2;
q9 = 1.0+2.0*b*t2;
h0b = -bet*g3*b*t6*(2.0+b*t2)/q8;
h0rs = -rsthrd*h0b*bec*ecrs;
fact0 = 2.0*delt-6.0*b;
fact1 = q5*q9+q4*q9*q9;
h0bt = 2.0*bet*g3*t4*((q4*q5*fact0-delt*fact1)/q8)/q8;
h0rst = rsthrd*t2*h0bt*bec*ecrs;
h0z = 3.0*gz*h0/g + h0b*(bg*gz+bec*eczet);
h0t = 2.0*bet*g3*q9/q8;
h0zt = 3.0*gz*h0t/g+h0bt*(bg*gz+bec*eczet);
fact2 = q4*q5+b*t2*(q4*q9+q5);
fact3 = 2.0*b*q5*q9+delt*fact2;
h0tt = 4.0*bet*g3*t*(2.0*b/q8-(q9*fact3/q8)/q8);
h1rs = r3*r2*t2*(-r4+r1*t2/3.0);
fact4 = 2.0-r1*t2;
h1rst = r3*r2*t2*(2.0*r4*(1.0-r1*t2)-thrd2*r1*t2*fact4);
h1z = gz*r3*r2*t2*(3.0-4.0*r1*t2)/g;
h1t = 2.0*r3*r2*(1.0-r1*t2);
h1zt = 2.0*gz*r3*r2*(3.0-11.0*r1*t2+4.0*r1*r1*t4)/g;
h1tt = 4.0*r3*r2*r1*t*(-2.0+r1*t2);
hrs = h0rs+h1rs;
hrst = h0rst+h1rst;
ht = h0t+h1t;
htt = h0tt+h1tt;
hz = h0z+h1z;
hzt = h0zt+h1zt;
comm = h+hrs+hrst+t2*ht/6.0+7.0*t2*t*htt/6.0;
pref = hz-gz*t2*ht/g;
fact5 = gz*(2.0*ht+t*htt)/g;
comm = comm-pref*zet-uu*htt-vv*ht-ww*(hzt-fact5);
return (comm+pref);
}
*/
//IPD: This works for Plane-Waves, where n is complex
// Keep it basis dependent
// If needed for wavelets, generalize it later
/*
RealSpaceColumn
ex_GC(RealSpaceColumn &n)
{
int Nx,Ny,Nz,Nx2,Ny2,Nz2,NyNz;
int i,j,k,index;
Basis *basis = n.basis;
Lattice *lattice = basis->lattice;
Nx = basis->Nx;
Ny = basis->Ny;
Nz = basis->Nz;
Nx2 = Nx/2;
Ny2 = Ny/2;
Nz2 = Nz/2;
NyNz = Ny * Nz;
int NxNyNz = basis->NxNyNz;
real *d;
d = (real *) mymalloc(sizeof(real)*n.length, "d", "exGC");
RealSpaceColumn exGC(n.length,n.basis);
basis = n.basis;
CoeffSpaceColumn Jn(basis->NxNyNz, basis);
// regulating charge density
//
// one of the two ways in doing the cutoff.
// The relevant value in GGA is grad(n)/n
// Big error can occur if n is too small.
// 1. One way of cuting off is what's done here,
// by providing a lower bound on n, so that
// we deal with grad(n_new)/n_new.
// 2. The second way of doing this is by dealing
// with grad(n)/n_new.
// or you can provide your own versions here.
for (i = 0; i < NxNyNz; i++) {
if (REAL(n.data.d[i]) < CHD_CUTOFF)
exGC.data.d[i] = CHD_CUTOFF;
else
exGC.data.d[i] = n.data.d[i];
}
apply_J(exGC, Jn);
real g0, g1, g2;
int l;
real dtemp;
for (i = 0; i < NxNyNz; i++)
d[i] = 0.0;
for (l = 0; l < 3; l++) {
g0 = lattice->G.m[0][l];
g1 = lattice->G.m[1][l];
g2 = lattice->G.m[2][l];
for (i=-Nx2; i < Nx2; i++)
for (j=-Ny2; j < Ny2; j++)
for (k=-Nz2; k < Nz2; k++)
{
index = 0;
if (k < 0) index += k+Nz; else index += k;
if (j < 0) index += Nz*(j+Ny); else index += Nz*j;
if (i < 0) index += NyNz*(i+Nx); else index += NyNz*i;
exGC.data.d[index] = (i * g0 + j * g1 + k * g2) * Jn.data.d[index];
}
exGC = I(exGC);
for (i = 0; i < NxNyNz; i++) {
dtemp = IMAG(exGC.data.d[i]);
d[i] += dtemp * dtemp;
}
}
// Preparing Variables
real rho, rs, kf, ks, t, sa;
real fk, sk;
real zet, ga, ww;
real ec, vcup, vcdn, ecrs, eczet, alfc;
const real cc0 = 1.0/3.0;
const real cc1 = 1.919158293;
const real cc2 = 1.563185284;
const real cc3 = 0.6203504909;
for (i = 0; i < NxNyNz; i++) {
rho = REAL(n.data.d[i]);
if (rho < CHD_CUTOFF) rho = CHD_CUTOFF;
rs = cc3 * pow(rho, -cc0 );
kf = cc1 / rs;
ks = cc2 / sqrt(rs);
// consistency
fk = kf;
sk = ks;
// spin polarization not implemented
zet = 0.0;
ga = 1.0;
ww = 0.0;
t = sqrt(d[i]);
sa = t / (2.0 * kf * rho);
corlsd(rs, zet, ec, vcup, vcdn, ecrs, eczet, alfc);
t = t/(2.0*rho*ks*ga);
exGC.data.d[i] = exe(rho, sa) + core(rs, zet, t, ga, ec, sk, fk);
}
myfree(d);
return exGC;
}
*/
// IPD: There is conversion between Real/Coeff Coluumns,
// which may be separate classes. Also uses INDEX, which may nned to be
// handled as **map()** in our new concept
/*
RealSpaceColumn
exprime_GC(RealSpaceColumn& n)
{
int Nx,Ny,Nz,Nx2,Ny2,Nz2,NyNz;
int i,j,k,index;
Basis *basis = n.basis;
Lattice *lattice = basis->lattice;
Nx = basis->Nx;
Ny = basis->Ny;
Nz = basis->Nz;
Nx2 = Nx/2;
Ny2 = Ny/2;
Nz2 = Nz/2;
NyNz = Ny * Nz;
int NxNyNz = basis->NxNyNz;
real *(d[3]);
real *uu;
real *vv;
for (i = 0; i < 3; i++)
d[i] = (real *) mymalloc(sizeof(real)*n.length, "d[]", "exGCprime");
uu = (real *) mymalloc(sizeof(real)*n.length, "uu", "exGC");
vv = (real *) mymalloc(sizeof(real)*n.length, "vv", "exGC");
RealSpaceColumn exGCprime(n.length,n.basis); // used temporary sometimes...
basis = n.basis;
CoeffSpaceColumn Jn(basis->NxNyNz, basis);
// regulating charge density
for (i = 0; i < NxNyNz; i++) {
if (REAL(n.data.d[i]) < CHD_CUTOFF)
exGCprime.data.d[i] = CHD_CUTOFF;
else
exGCprime.data.d[i] = n.data.d[i];
}
//IPD: The second apply_J() call overrides the first:
dft_log("***** exprime_GC() unnecessary apply_J!!!");
apply_J(exGCprime, Jn);
apply_J(n, Jn);
real g0, g1, g2, q0, q1, q2, gl, gm;
int l, m;
for (l = 0; l < 3; l++) {
g0 = lattice->G.m[0][l];
g1 = lattice->G.m[1][l];
g2 = lattice->G.m[2][l];
for (i=-Nx2; i < Nx2; i++)
for (j=-Ny2; j < Ny2; j++)
for (k=-Nz2; k < Nz2; k++)
{
index = 0;
if (k < 0) index += k+Nz; else index += k;
if (j < 0) index += Nz*(j+Ny); else index += Nz*j;
if (i < 0) index += NyNz*(i+Nx); else index += NyNz*i;
exGCprime.data.d[index] = (i * g0 + j * g1 + k * g2) * Jn.data.d[index];
}
exGCprime = I(exGCprime);
for (i = 0; i < NxNyNz; i++)
d[l][i] = IMAG(exGCprime.data.d[i]);
}
for (i = 0; i < NxNyNz; i++)
vv[i] = uu[i] = 0.0;
// diagonal part
real ddtemp;
for (l = 0; l < 3; l++) {
g0 = lattice->G.m[0][l];
g1 = lattice->G.m[1][l];
g2 = lattice->G.m[2][l];
for (i=-Nx2; i < Nx2; i++)
for (j=-Ny2; j < Ny2; j++)
for (k=-Nz2; k < Nz2; k++)
{
index = 0;
if (k < 0) index += k+Nz; else index += k;
if (j < 0) index += Nz*(j+Ny); else index += Nz*j;
if (i < 0) index += NyNz*(i+Nx); else index += NyNz*i;
gl = (i * g0 + j * g1 + k * g2);
exGCprime.data.d[index] = - gl * gl * Jn.data.d[index];
}
exGCprime = I(exGCprime);
for (i = 0; i < NxNyNz; i++) {
ddtemp = REAL(exGCprime.data.d[i]);
vv[i] += ddtemp;
uu[i] += d[l][i] * d[l][i] * ddtemp;
}
for (m = l + 1; m < 3; m++) {
q0 = lattice->G.m[0][m];
q1 = lattice->G.m[1][m];
q2 = lattice->G.m[2][m];
for (i=-Nx2; i < Nx2; i++)
for (j=-Ny2; j < Ny2; j++)
for (k=-Nz2; k < Nz2; k++)
{
index = 0;
if (k < 0) index += k+Nz; else index += k;
if (j < 0) index += Nz*(j+Ny); else index += Nz*j;
if (i < 0) index += NyNz*(i+Nx); else index += NyNz*i;
gl = (i * g0 + j * g1 + k * g2);
gm = (i * q0 + j * q1 + k * q2);
exGCprime.data.d[index] = - gl * gm * Jn.data.d[index];
}
exGCprime = I(exGCprime);
for (i = 0; i < NxNyNz; i++) {
ddtemp = REAL(exGCprime.data.d[i]);
uu[i] += 2.0 * d[l][i] * d[m][i] * ddtemp;
}
}
}
// Preparing Variables
real rho, sa, rs, kf, ks, t;
real fk, sk;
real zet, ga, ww;
real ec, vcup, vcdn, ecrs, eczet, alfc;
real u, ub, delta, v;
const real cc0 = 1.0/3.0;
const real cc1 = 1.919158293;
const real cc2 = 1.563185284;
const real cc3 = 0.6203504909;
for (i = 0; i < NxNyNz; i++) {
rho = REAL(n.data.d[i]);
if (rho < CHD_CUTOFF) rho = CHD_CUTOFF;
rs = cc3 * pow(rho, -cc0);
kf = cc1 / rs;
ks = cc2/ sqrt(rs);
// consistency
fk = kf;
sk = ks;
// spin polarization not implemented
zet = 0.0;
ga = 1.0;
ww = 0.0;
t = 0.0;
for (j = 0; j < 3; j++)
t += d[j][i] * d[j][i];
t = sqrt(t);
if (t > 0.0) {
u = uu[i] / (t * rho * rho * 8.0 * kf * kf * kf);
ub = uu[i] / t;
}
else {
u = 0.0;
uu[i] = 0.0;
ub = 0.0;
}
delta = vv[i];
v = delta / (rho * 4.0 * kf * kf);
vv[i] = vv[i] / (rho*pow(2.0*ks*ga,2));
sa = t / (2.0 * kf * rho);
corlsd(rs, zet, ec, vcup, vcdn, ecrs, eczet, alfc);
ub = ub/(rho*rho*pow(2.0*ks*ga,3));
t = t/(2.0*rho*ks*ga);
exGCprime.data.d[i] = exvx(rho, sa, u, v) +
corvx(rs, zet, t, ub, vv[i], ww, ga, ec, eczet, ecrs, sk, fk);
}
myfree(d[0]);
myfree(d[1]);
myfree(d[2]);
myfree(uu);
myfree(vv);
return exGCprime;
}
*/
// #############################################################
// Jason Cline's Teter LSD. I only changed a few cosmetic things.
// #############################################################
/* this is a sub routine that calculates either the exchange correlation
* energy density, or its derivative, depending on the flag for the LSD_TETER
* functional
*/
#define EXC 0
#define EXCPRIME 1
static void
exc1( RealSpaceScalarFieldColumn &n_up, RealSpaceScalarFieldColumn &n_dn,
RealSpaceScalarFieldColumn &out_up, RealSpaceScalarFieldColumn &out_dn,
real vol, int prime_flag)
{
// Mike Teter's parameters of 8 April 1993.
//(a0=(3/4)(3/(2 Pi))^(2/3))
//New parameters which accomodate spin polarization (fit to P-W)
//Paramagmetic limit:
const real a0p=.4581652932831429e0,
a1p=2.217058676663745e0,
a2p=0.7405551735357053e0,
a3p=0.01968227878617998e0,
b1p=1.0e0,
b2p=4.504130959426697e0,
b3p=1.110667363742916e0,
b4p=0.02359291751427506e0;
// Differences, ferromagnetic - paramagnetic (delta params):
const real da0=.119086804055547e0,
da1=0.6157402568883345e0,
da2=0.1574201515892867e0,
da3=0.003532336663397157e0;
// (note that b1=1 is fixed)
const real db1=0.0e0,
db2=0.2673612973836267e0,
db3=0.2052004607777787e0,
db4=0.004200005045691381e0;
const real ft=4.0/3.0,
ot=1.0/3.0;
const real PI = M_PI;
int i;
for(i = 0; i < n_up.length; i++){
real rs, zeta, rsprime, ntotsqinv, ntotinv;
scalar n_tot = n_up.data.d[i] + n_dn.data.d[i];
ntotinv = 1. / REAL(n_tot);
if (REAL(n_tot) <= 0.0)
rs = 1.0e15;
else
rs = pow(ft*PI*REAL(n_tot),(-ot));
zeta = REAL((n_up.data.d[i]-n_dn.data.d[i])/n_tot); // difference / sum
// Exchange energy spin interpolation function f(zeta)
real fxc;
real a0,a1,a2,a3,b1,b2,b3,b4;
fxc=(pow((1.0+zeta),ft)+pow((1.0-zeta),ft)-2.0)/(pow(2.0,ft)-2.0);
a0 = a0p+fxc*da0;
a1 = a1p+fxc*da1;
a2 = a2p+fxc*da2;
a3 = a3p+fxc*da3;
b1 = b1p+fxc*db1;
b2 = b2p+fxc*db2;
b3 = b3p+fxc*db3;
b4 = b4p+fxc*db4;
real n1, d1, exc;
n1 = -(a0+rs*(a1+rs*(a2+rs*a3)));
d1 = rs*(b1+rs*(b2+rs*(b3+rs*b4)));
// Exchange-correlation energy
exc = n1/d1;
// Exchange-correlation potential
real n2a, n2b, vxcp, fxcp, n3a, n3b, dexcdz;
n2a = -(a1+rs*(2.0*a2+rs*(3.0*a3)));
n2b = -(b1+rs*(2.0*b2+rs*(3.0*b3+rs*(4.0*b4))));
// only vxcp contributes when paramagnetic
vxcp = exc - ot*rs*(n2a+exc*n2b)/d1;
// dexcdz is d(exc)/d(zeta)
// proportional to d(fxc)/d(zeta) (which is 0 at zeta=0)
fxcp = ft*(pow((1.0+zeta),ot)-pow((1.0-zeta),ot))/
(pow(2.0,ft)-2.0);
n3a = -(da0+rs*(da1+rs*(da2+rs*da3)));
n3b = -rs*(db1+rs*(db2+rs*(db3+rs*db4)));
dexcdz = (n3a+exc*n3b)*fxcp/d1;
// if we are computing exc,
// then we only need an energy density per electron,
// which doesn't depend on electron spin
switch(prime_flag){
case EXC:
out_up.data.d[i] = exc;
out_dn.data.d[i] = exc;
break;
case EXCPRIME:
rsprime = -ot*ft*PI*pow(ft*PI*REAL(n_tot),-ft);
ntotsqinv = 1./REAL(n_tot*n_tot);
out_up.data.d[i] = (scalar)((dexcdz)*(ntotinv-REAL(n_up.data.d[i]-n_dn.data.d[i])*ntotsqinv)
+ ((n2a+exc*n2b)/d1)*(rsprime));
out_dn.data.d[i] = (scalar)((dexcdz)*(-ntotinv-REAL(n_up.data.d[i]-n_dn.data.d[i])*ntotsqinv)
+ ((n2a+exc*n2b)/d1)*(rsprime));
break;
default:
die("you idiot: either EXC or EXCPRIME into exc1\n");
}
}
}
// Cline's Teter LSD
RealSpaceScalarFieldColumn
exc_LSD_TETER(RealSpaceScalarFieldColumn &n_up,
RealSpaceScalarFieldColumn &n_dn)
{
real vol = n_up.basis_spec->lattice->unit_cell_volume;
RealSpaceScalarFieldColumn temp(n_up);
exc1(n_up,n_dn,temp, temp,vol,EXC);
return(temp);
}
// Cline's Teter LSD
void
excprime_LSD_TETER(RealSpaceScalarFieldColumn &n_up,
RealSpaceScalarFieldColumn &n_dn,
RealSpaceScalarFieldColumn &excprime_up,
RealSpaceScalarFieldColumn &excprime_dn)
{
real vol = n_up.basis_spec->lattice->unit_cell_volume;
exc1(n_up,n_dn,excprime_up, excprime_dn, vol,EXCPRIME);
}
syntax highlighted by Code2HTML, v. 0.9.1