/* 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); }