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 // BHandHLYPFunctional.cpp
16 //
17 ////////////////////////////////////////////////////////////////////////////////
18
19 #include <cmath>
20 #include <cassert>
21 #include "BHandHLYPFunctional.h"
22 #include "BLYPFunctional.h"
23 using namespace std;
24
25 ////////////////////////////////////////////////////////////////////////////////
BHandHLYPFunctional(const vector<vector<double>> & rhoe)26 BHandHLYPFunctional::BHandHLYPFunctional(const vector<vector<double> > &rhoe)
27 {
28 _nspin = rhoe.size();
29 if ( _nspin > 1 ) assert(rhoe[0].size() == rhoe[1].size());
30 _np = rhoe[0].size();
31
32 if ( _nspin == 1 )
33 {
34 _exc.resize(_np);
35 _vxc1.resize(_np);
36 _vxc2.resize(_np);
37 _grad_rho[0].resize(_np);
38 _grad_rho[1].resize(_np);
39 _grad_rho[2].resize(_np);
40 rho = &rhoe[0][0];
41 grad_rho[0] = &_grad_rho[0][0];
42 grad_rho[1] = &_grad_rho[1][0];
43 grad_rho[2] = &_grad_rho[2][0];
44 exc = &_exc[0];
45 vxc1 = &_vxc1[0];
46 vxc2 = &_vxc2[0];
47 }
48 else
49 {
50 _exc_up.resize(_np);
51 _exc_dn.resize(_np);
52 _vxc1_up.resize(_np);
53 _vxc1_dn.resize(_np);
54 _vxc2_upup.resize(_np);
55 _vxc2_updn.resize(_np);
56 _vxc2_dnup.resize(_np);
57 _vxc2_dndn.resize(_np);
58 _grad_rho_up[0].resize(_np);
59 _grad_rho_up[1].resize(_np);
60 _grad_rho_up[2].resize(_np);
61 _grad_rho_dn[0].resize(_np);
62 _grad_rho_dn[1].resize(_np);
63 _grad_rho_dn[2].resize(_np);
64
65 rho_up = &rhoe[0][0];
66 rho_dn = &rhoe[1][0];
67 grad_rho_up[0] = &_grad_rho_up[0][0];
68 grad_rho_up[1] = &_grad_rho_up[1][0];
69 grad_rho_up[2] = &_grad_rho_up[2][0];
70 grad_rho_dn[0] = &_grad_rho_dn[0][0];
71 grad_rho_dn[1] = &_grad_rho_dn[1][0];
72 grad_rho_dn[2] = &_grad_rho_dn[2][0];
73 exc_up = &_exc_up[0];
74 exc_dn = &_exc_dn[0];
75 vxc1_up = &_vxc1_up[0];
76 vxc1_dn = &_vxc1_dn[0];
77 vxc2_upup = &_vxc2_upup[0];
78 vxc2_updn = &_vxc2_updn[0];
79 vxc2_dnup = &_vxc2_dnup[0];
80 vxc2_dndn = &_vxc2_dndn[0];
81 }
82 }
83
84 ////////////////////////////////////////////////////////////////////////////////
setxc(void)85 void BHandHLYPFunctional::setxc(void)
86 {
87 if ( _np == 0 ) return;
88
89 if ( _nspin == 1 )
90 {
91 assert( rho != 0 );
92 assert( grad_rho[0] != 0 && grad_rho[1] != 0 && grad_rho[2] != 0 );
93 assert( exc != 0 );
94 assert( vxc1 != 0 );
95 assert( vxc2 != 0 );
96 for ( int i = 0; i < _np; i++ )
97 {
98 double grad = sqrt(grad_rho[0][i]*grad_rho[0][i] +
99 grad_rho[1][i]*grad_rho[1][i] +
100 grad_rho[2][i]*grad_rho[2][i] );
101 excbhandhlyp(rho[i],grad,&exc[i],&vxc1[i],&vxc2[i]);
102 }
103 }
104 else
105 {
106 assert( rho_up != 0 );
107 assert( rho_dn != 0 );
108 assert( grad_rho_up[0] != 0 && grad_rho_up[1] != 0 && grad_rho_up[2] != 0 );
109 assert( grad_rho_dn[0] != 0 && grad_rho_dn[1] != 0 && grad_rho_dn[2] != 0 );
110 assert( exc_up != 0 );
111 assert( exc_dn != 0 );
112 assert( vxc1_up != 0 );
113 assert( vxc1_dn != 0 );
114 assert( vxc2_upup != 0 );
115 assert( vxc2_updn != 0 );
116 assert( vxc2_dnup != 0 );
117 assert( vxc2_dndn != 0 );
118
119 for ( int i = 0; i < _np; i++ )
120 {
121 double grx_up = grad_rho_up[0][i];
122 double gry_up = grad_rho_up[1][i];
123 double grz_up = grad_rho_up[2][i];
124 double grx_dn = grad_rho_dn[0][i];
125 double gry_dn = grad_rho_dn[1][i];
126 double grz_dn = grad_rho_dn[2][i];
127 double grad_up2 = grx_up*grx_up + gry_up*gry_up + grz_up*grz_up;
128 double grad_dn2 = grx_dn*grx_dn + gry_dn*gry_dn + grz_dn*grz_dn;
129 double grad_up_grad_dn = grx_up*grx_dn + gry_up*gry_dn + grz_up*grz_dn;
130
131 excbhandhlyp(rho_up[i],rho_dn[i],grad_up2,grad_dn2,grad_up_grad_dn,
132 &exc_up[i],&exc_dn[i],&vxc1_up[i],&vxc1_dn[i],
133 &vxc2_upup[i],&vxc2_dndn[i],
134 &vxc2_updn[i],&vxc2_dnup[i]);
135 }
136 }
137 }
138
139 ////////////////////////////////////////////////////////////////////////////////
excbhandhlyp(double rho,double grad,double * exc,double * vxc1,double * vxc2)140 void BHandHLYPFunctional::excbhandhlyp(double rho, double grad,
141 double *exc, double *vxc1, double *vxc2)
142 {
143 // BHandHLYP unpolarized
144 // Coefficients of the BHandHLYP functional
145 // A. Becke, JCP 98, 5648 (1993)
146 // See also X.Xu and W. Goddard, J.Phys.Chem. A108, 2305 (2004)
147 // EcLSDA is the Vosko-Wilk-Nusair correlation energy
148 // dExBecke88 is the difference ExB88 - ExLDA
149 // B3LYP: 0.2 ExHF + 0.8 ExSlater + 0.19 EcLSDA + 0.81 EcLYP + 0.72 dExBecke88
150 // BHandHLYP: 0.5 ExHF + 0.5 ExB88 + 1.0 EcLYP
151 const double xb88_coeff = 0.5; // Becke88 exchange
152 const double clyp_coeff = 1.0 ;
153
154 double ex_b88,vx1_b88,vx2_b88;
155 double ec_lyp,vc1_lyp,vc2_lyp;
156
157 BLYPFunctional::exb88(rho,grad,&ex_b88,&vx1_b88,&vx2_b88);
158 BLYPFunctional::eclyp(rho,grad,&ec_lyp,&vc1_lyp,&vc2_lyp);
159
160 *exc = xb88_coeff * ex_b88 +
161 clyp_coeff * ec_lyp;
162
163 *vxc1 = xb88_coeff * vx1_b88 +
164 clyp_coeff * vc1_lyp;
165
166 *vxc2 = xb88_coeff * vx2_b88 +
167 clyp_coeff * vc2_lyp;
168 }
169
170 ////////////////////////////////////////////////////////////////////////////////
excbhandhlyp(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)171 void BHandHLYPFunctional::excbhandhlyp(double rho_up, double rho_dn,
172 double grad_up2, double grad_dn2, double grad_up_grad_dn,
173 double *exc_up, double *exc_dn, double *vxc1_up, double *vxc1_dn,
174 double *vxc2_upup, double *vxc2_dndn, double *vxc2_updn, double *vxc2_dnup)
175 {
176 // BHandHLYP spin-polarized
177 // Coefficients of the BHandHLYP functional
178 // A. Becke, JCP 98, 5648 (1993)
179 // See also X.Xu and W. Goddard, J.Phys.Chem. A108, 2305 (2004)
180 // EcLSDA is the Vosko-Wilk-Nusair correlation energy
181 // dExBecke88 is the difference ExB88 - ExLDA
182 // B3LYP: 0.2 ExHF + 0.8 ExSlater + 0.19 EcLSDA + 0.81 EcLYP + 0.72 dExBecke88
183 // BHandHLYP: 0.5 ExHF + 0.5 ExB88 + 1.0 EcLYP
184 const double xb88_coeff = 0.5; // Becke88 exchange
185 const double clyp_coeff = 1.0;
186
187 double ex_b88_up,ex_b88_dn,vx1_b88_up,vx1_b88_dn,
188 vx2_b88_upup,vx2_b88_dndn,vx2_b88_updn,vx2_b88_dnup;
189 double ec_lyp_up,ec_lyp_dn,vc1_lyp_up,vc1_lyp_dn,
190 vc2_lyp_upup,vc2_lyp_dndn,vc2_lyp_updn,vc2_lyp_dnup;
191
192 BLYPFunctional::exb88_sp(rho_up,rho_dn,grad_up2,grad_dn2,grad_up_grad_dn,
193 &ex_b88_up,&ex_b88_dn,&vx1_b88_up,&vx1_b88_dn,
194 &vx2_b88_upup,&vx2_b88_dndn,&vx2_b88_updn,&vx2_b88_dnup);
195 BLYPFunctional::eclyp_sp(rho_up,rho_dn,grad_up2,grad_dn2,grad_up_grad_dn,
196 &ec_lyp_up,&ec_lyp_dn,&vc1_lyp_up,&vc1_lyp_dn,
197 &vc2_lyp_upup,&vc2_lyp_dndn,&vc2_lyp_updn,&vc2_lyp_dnup);
198
199 *exc_up = xb88_coeff * ex_b88_up +
200 clyp_coeff * ec_lyp_up;
201
202 *exc_dn = xb88_coeff * ex_b88_dn +
203 clyp_coeff * ec_lyp_dn;
204
205 *vxc1_up = xb88_coeff * vx1_b88_up +
206 clyp_coeff * vc1_lyp_up;
207
208 *vxc1_dn = xb88_coeff * vx1_b88_dn +
209 clyp_coeff * vc1_lyp_dn;
210
211 *vxc2_upup = xb88_coeff * vx2_b88_upup + clyp_coeff * vc2_lyp_upup;
212 *vxc2_dndn = xb88_coeff * vx2_b88_dndn + clyp_coeff * vc2_lyp_dndn;
213 *vxc2_updn = xb88_coeff * vx2_b88_updn + clyp_coeff * vc2_lyp_updn;
214 *vxc2_dnup = xb88_coeff * vx2_b88_dnup + clyp_coeff * vc2_lyp_dnup;
215 }
216