1 ////////////////////////////////////////////////////////////////////////////////
2 //
3 // Copyright (c) 2008 The Regents of the University of California
4 //
5 // This file is part of Qbox
6 //
7 // Qbox is distributed under the terms of the GNU General Public License
8 // as published by the Free Software Foundation, either version 2 of
9 // the License, or (at your option) any later version.
10 // See the file COPYING in the root directory of this distribution
11 // or <http://www.gnu.org/licenses/>.
12 //
13 ////////////////////////////////////////////////////////////////////////////////
14 //
15 // LDAFunctional.cpp
16 //
17 // LDA Exchange-correlation energy and potential
18 // Ceperley & Alder, parametrized by Perdew and Zunger
19 //
20 ////////////////////////////////////////////////////////////////////////////////
21
22 #include <cmath>
23 #include <cassert>
24 #include <vector>
25 #include "LDAFunctional.h"
26
setxc(void)27 void LDAFunctional::setxc(void)
28 {
29 if ( _np == 0 ) return;
30 if ( _nspin == 1 )
31 {
32 assert(rho != 0);
33 assert(exc != 0);
34 assert(vxc1 != 0);
35 #pragma omp parallel for
36 for ( int ir = 0; ir < _np; ir++ )
37 {
38 xc_unpolarized(rho[ir],exc[ir],vxc1[ir]);
39 }
40 }
41 else
42 {
43 // spin polarized
44 assert(rho_up != 0);
45 assert(rho_dn != 0);
46 assert(exc != 0);
47 assert(vxc1_up != 0);
48 assert(vxc1_dn != 0);
49 const double fz_prefac = 1.0 / ( cbrt(2.0)*2.0 - 2.0 );
50 const double dfz_prefac = (4.0/3.0) * fz_prefac;
51 #pragma omp parallel for
52 for ( int ir = 0; ir < _np; ir++ )
53 {
54 double excir = 0.0;
55 double v_up = 0.0;
56 double v_dn = 0.0;
57
58 double roe_up = rho_up[ir];
59 double roe_dn = rho_dn[ir];
60 double roe = roe_up + roe_dn;
61
62 if ( roe > 0.0 )
63 {
64 double zeta = ( roe_up - roe_dn ) / roe;
65
66 double zp1 = 1.0 + zeta;
67 double zm1 = 1.0 - zeta;
68 double zp1_13 = cbrt(zp1);
69 double zm1_13 = cbrt(zm1);
70 double fz = fz_prefac * ( zp1_13 * zp1 + zm1_13 * zm1 - 2.0 );
71 double dfz = dfz_prefac * ( zp1_13 - zm1_13 );
72
73 double xc_u, xc_p, v_u, v_p;
74 xc_unpolarized(roe,xc_u,v_u);
75 xc_polarized(roe,xc_p,v_p);
76
77 double xc_pu = xc_p - xc_u;
78 excir = xc_u + fz * xc_pu;
79
80 double v = v_u + fz * ( v_p - v_u );
81 v_up = v + xc_pu * ( 1.0 - zeta ) * dfz;
82 v_dn = v + xc_pu * ( -1.0 - zeta ) * dfz;
83 }
84
85 vxc1_up[ir] = v_up;
86 vxc1_dn[ir] = v_dn;
87 exc[ir] = excir;
88 }
89 }
90 }
91
xc_unpolarized(const double rh,double & ee,double & vv)92 void LDAFunctional::xc_unpolarized(const double rh, double &ee, double &vv)
93 {
94 // compute LDA xc energy and potential, unpolarized
95 // const double third=1.0/3.0;
96 // c1 is (3.D0/(4.D0*pi))**third
97 const double c1 = 0.6203504908994001;
98 // alpha = (4/(9*pi))**third = 0.521061761198
99 // const double alpha = 0.521061761198;
100 // c2 = -(3/(4*pi)) / alpha = -0.458165293283
101 // const double c2 = -0.458165293283;
102 // c3 = (4/3) * c2 = -0.610887057711
103 const double c3 = -0.610887057711;
104
105 const double A = 0.0311;
106 const double B = -0.048;
107 const double b1 = 1.0529;
108 const double b2 = 0.3334;
109 const double G = -0.1423;
110
111 // C from the PZ paper: const double C = 0.0020;
112 // D from the PZ paper: const double D = -0.0116;
113 // C and D by matching Ec and Vc at rs=1
114 const double D = G / ( 1.0 + b1 + b2 ) - B;
115 const double C = -A - D - G * ( (b1/2.0 + b2) / ((1.0+b1+b2)*(1.0+b1+b2)));
116
117 ee = 0.0;
118 vv = 0.0;
119
120 if ( rh > 0.0 )
121 {
122 double ro13 = cbrt(rh);
123 double rs = c1 / ro13;
124
125 double ex=0.0,vx=0.0,ec=0.0,vc=0.0;
126
127 // Next line : exchange part in Hartree units
128 vx = c3 / rs;
129 ex = 0.75 * vx;
130
131 // Next lines : Ceperley & Alder correlation (Zunger & Perdew)
132 if ( rs < 1.0 )
133 {
134 double logrs = log(rs);
135 ec = A * logrs + B + C * rs * logrs + D * rs;
136 vc = A * logrs + ( B - A / 3.0 ) +
137 (2.0/3.0) * C * rs * logrs +
138 ( ( 2.0 * D - C ) / 3.0 ) * rs;
139 }
140 else
141 {
142 double sqrtrs = sqrt(rs);
143 double den = 1.0 + b1 * sqrtrs + b2 * rs;
144 ec = G / den;
145 vc = ec * ( 1.0 + (7.0/6.0) * b1 * sqrtrs +
146 (4.0/3.0) * b2 * rs ) / den;
147 }
148 ee = ex + ec;
149 vv = vx + vc;
150 }
151 }
152
xc_polarized(const double rh,double & ee,double & vv)153 void LDAFunctional::xc_polarized(const double rh, double &ee, double &vv)
154 {
155 // compute LDA polarized XC energy and potential
156
157 // const double third=1.0/3.0;
158 // c1 is (3.D0/(4.D0*pi))**third
159 const double c1 = 0.6203504908994001;
160 // alpha = (4/(9*pi))**third = 0.521061761198
161 // const double alpha = 0.521061761198;
162 // c2 = -(3/(4*pi)) / alpha = -0.458165293283
163 // const double c2 = -0.458165293283;
164 // c3 = (4/3) * c2 = -0.610887057711
165 // const double c3 = -0.610887057711;
166 // c4 = 2**third * c3
167 const double c4 = -0.769669463118;
168
169 const double A = 0.01555;
170 const double B = -0.0269;
171 const double b1 = 1.3981;
172 const double b2 = 0.2611;
173 const double G = -0.0843;
174 // C from PZ paper: const double C = 0.0007;
175 // D from PZ paper: const double D = -0.0048;
176 // C and D by matching Ec and Vc at rs=1
177 const double D = G / ( 1.0 + b1 + b2 ) - B;
178 const double C = -A - D - G * ( (b1/2.0 + b2) / ((1.0+b1+b2)*(1.0+b1+b2)));
179
180 ee = 0.0;
181 vv = 0.0;
182
183 if ( rh > 0.0 )
184 {
185 double ro13 = cbrt(rh);
186 double rs = c1 / ro13;
187
188 double ex=0.0,vx=0.0,ec=0.0,vc=0.0;
189
190 // Next line : exchange part in Hartree units
191 vx = c4 / rs;
192 ex = 0.75 * vx;
193
194 // Next lines : Ceperley & Alder correlation (Zunger & Perdew)
195 if ( rs < 1.0 )
196 {
197 double logrs = log(rs);
198 ec = A * logrs + B + C * rs * logrs + D * rs;
199 vc = A * logrs + ( B - A / 3.0 ) +
200 (2.0/3.0) * C * rs * logrs +
201 ( ( 2.0 * D - C ) / 3.0 ) * rs;
202 }
203 else
204 {
205 double sqrtrs = sqrt(rs);
206 double den = 1.0 + b1 * sqrtrs + b2 * rs;
207 ec = G / den;
208 vc = ec * ( 1.0 + (7.0/6.0) * b1 * sqrtrs +
209 (4.0/3.0) * b2 * rs ) / den;
210 }
211 ee = ex + ec;
212 vv = vx + vc;
213 }
214 }
215