1 /*
2 This file is part of MADNESS.
3
4 Copyright (C) 2007,2010 Oak Ridge National Laboratory
5
6 This program is free software; you can redistribute it and/or modify
7 it under the terms of the GNU General Public License as published by
8 the Free Software Foundation; either version 2 of the License, or
9 (at your option) any later version.
10
11 This program is distributed in the hope that it will be useful,
12 but WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 GNU General Public License for more details.
15
16 You should have received a copy of the GNU General Public License
17 along with this program; if not, write to the Free Software
18 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
19
20 For more information please contact:
21
22 Robert J. Harrison
23 Oak Ridge National Laboratory
24 One Bethel Valley Road
25 P.O. Box 2008, MS-6367
26
27 email: harrisonrj@ornl.gov
28 tel: 865-241-3937
29 fax: 865-572-0680
30
31 $Id$
32 */
33 //#define WORLD_INSTANTIATE_STATIC_TEMPLATES
34 #include <madness/world/world.h>
35 #include <madness/mra/mra.h>
36 #include <madness/tensor/tensor.h>
37 #include <fstream>
38 #include "xcfunctional.h"
39
40 using namespace madness;
41
42 static std::string df_repo_functionals[] = {
43 "lda_x",
44 "lda_c_vwn",
45 "gga_x_pbe",
46 "gga_c_pbe",
47 "gga_x_pw91",
48 "gga_c_pw91",
49 "gga_x_b88",
50 "gga_c_lyp",
51 "lda_c_vwn_rpa",
52 "lda_c_vwn",
53 "lda_c_pz",
54 "lda_c_pw",
55 "hyb_gga_xc_b3lyp",
56 "gga_xc_hcth_93",
57 "gga_xc_hcth_407",
58 "gga_xc_hcth_147",
59 "gga_xc_hcth_120",
60 "gga_xc_edf1",
61 "gga_xc_b97_2",
62 "gga_xc_b97_1",
63 "gga_xc_b97",
64 "gga_x_pw91",
65 "gga_x_pbe",
66 "gga_x_ft97_b",
67 "gga_x_b88",
68 "gga_c_pw91",
69 "gga_c_pbe",
70 "gga_c_p86",
71 "gga_c_lyp"};
72
73
74 struct xcfunc_data_point
75 {
76 double rhoa, rhob;
77 double sigmaaa, sigmaab, sigmabb;
78 double zk;
79 double vrhoa, vrhob;
80 double vsigmaaa, vsigmaab, vsigmabb;
81 double v2rhoa2, v2rhoab, v2rhob2;
82 };
83
read_test_data(const std::string & dfname,bool spin_polarized)84 std::vector<xcfunc_data_point> read_test_data(const std::string& dfname,
85 bool spin_polarized)
86 {
87
88 std::vector<xcfunc_data_point> dps;
89
90 double rhoa=1e-25;
91 double siga;
92 while (rhoa<1e16)
93 {
94 siga=1e-25;
95
96 while (siga<1e16)
97 {
98 xcfunc_data_point dp;
99 dp.rhoa=rhoa;
100 dp.rhob=rhoa;
101 dp.sigmaaa=siga;
102 dp.sigmaab=siga;
103 dp.sigmabb=siga;
104
105 dp.zk=0.0;
106 dp.vrhoa=0.0;
107 dp.vrhob=0.0;
108 dp.vsigmaaa=0.0;
109 dp.vsigmaab=0.0;
110 dp.vsigmabb=0.0;
111 dp.v2rhoa2=0.0;
112 dp.v2rhoab=0.0;
113 dp.v2rhob2=0.0;
114
115 if (!spin_polarized)
116 {
117 if (std::abs(dp.rhoa-dp.rhob) <= 1e-10)
118 dps.push_back(dp);
119 }
120 else
121 {
122 dps.push_back(dp);
123 }
124
125 siga*=10;
126 }
127 rhoa*=10;
128 }
129
130 return dps;
131
132 }
133
test_xcfunctional(World & world)134 void test_xcfunctional(World& world)
135 {
136 bool spin_polarized ;
137 //spin_polarized = true;
138 spin_polarized = false;
139
140 int ispin = 0; //alpha=0 beta=1
141
142 //int istr = 0; // 0-7
143 for (int istr=2;istr<3;istr++){
144
145 //vama5 0 "lda_x",
146 //vama5 1 "lda_c_vwn_rpa",
147 //vama5 2 "lda_c_vwn",
148 //vama5 3 "lda_c_pz",
149 //vama5 4 "lda_c_pw",
150 //vama5 5 "hyb_gga_xc_b3lyp",
151 //vama5 6 "gga_xc_hcth_93",
152 //vama5 7 "gga_xc_hcth_407",
153 //vama5 8 "gga_xc_hcth_147",
154 //vama5 9 "gga_xc_hcth_120",
155 //vama5 10 "gga_xc_edf1",
156 //vama5 11 "gga_xc_b97_2",
157 //vama5 12 "gga_xc_b97_1",
158 //vama5 13 "gga_xc_b97",
159 //vama5 14 "gga_x_pw91",
160 //vama5 15 "gga_x_pbe",
161
162 //vama5 16 "gga_x_ft97_b",
163 //vama5 17 "gga_x_b88",
164 //vama5 18 "gga_c_pw91",
165 //vama5 19 "gga_c_pbe",
166
167 //vama5 20 "gga_c_p86",
168 //vama5 21 "gga_c_lyp"};
169
170 XCfunctional xcfunc;
171 std::string xcfuncstr = df_repo_functionals[istr];
172 std::cout << "Testing exchange-correlation functional: "<< xcfuncstr << std::endl;
173
174 xcfuncstr += " 1.0";
175 xcfuncstr += " gga_c_pbe 1.0"; //hack to do both x and c for PBE
176 xcfuncstr += " RHOTOL 0.0 RHOMIN 0.0 SIGTOL 0.0 SIGMIN 0.0";
177 xcfunc.initialize(xcfuncstr,spin_polarized,world);
178
179 std::string fpath("df_repo/");
180 fpath += df_repo_functionals[istr];
181 fpath += ".data";
182 std::vector<xcfunc_data_point> dps = read_test_data(fpath.c_str(),spin_polarized);
183
184 Tensor<double> rhoa_t((long)dps.size());
185 Tensor<double> rhob_t((long)dps.size());
186 Tensor<double> sigmaaa_t((long)dps.size());
187 Tensor<double> sigmaab_t((long)dps.size());
188 Tensor<double> sigmabb_t((long)dps.size());
189 Tensor<double> zk_t((long)dps.size());
190 Tensor<double> vrhoa_t((long)dps.size());
191 Tensor<double> vrhob_t((long)dps.size());
192 Tensor<double> vsigmaaa_t((long)dps.size());
193 Tensor<double> vsigmaab_t((long)dps.size());
194 Tensor<double> vsigmabb_t((long)dps.size());
195
196 Tensor<double> rhoa1_t((long)dps.size());
197 Tensor<double> sigmaaa1_t((long)dps.size());
198
199 std::vector<Tensor<double> > xc_args;
200 for (unsigned int idp = 0; idp < dps.size(); idp++)
201 {
202 rhoa_t(idp) = dps[idp].rhoa;
203 rhob_t(idp) = dps[idp].rhob;
204
205 sigmaaa_t(idp) = dps[idp].sigmaaa;
206 sigmaab_t(idp) = dps[idp].sigmaab;
207 sigmabb_t(idp) = dps[idp].sigmabb;
208
209 zk_t(idp) = dps[idp].zk;
210
211 vrhoa_t(idp) = dps[idp].vrhoa;
212 vrhob_t(idp) = dps[idp].vrhob;
213
214
215 vsigmaaa_t(idp) = dps[idp].vsigmaaa;
216 vsigmaab_t(idp) = dps[idp].vsigmaab;
217 vsigmabb_t(idp) = dps[idp].vsigmabb;
218 }
219 if (spin_polarized)
220 {
221 xc_args.push_back(rhoa_t);
222 xc_args.push_back(rhob_t);
223 if (xcfunc.is_gga()) {
224 xc_args.push_back(sigmaaa_t);
225 xc_args.push_back(sigmaab_t);
226 xc_args.push_back(sigmabb_t);
227 }
228 }
229 else
230 {
231 xc_args.push_back(rhoa_t);
232 if (xcfunc.is_gga()) {
233 xc_args.push_back(sigmaaa_t);
234 }
235 }
236
237 print("xc_args_size", xc_args.size());
238 print("ispin ", ispin);
239 print("spin polarized ", spin_polarized);
240 //vama1 std::cout << "Testing spin-polarized case: " << std::endl << std::endl;
241
242 // local vr[0] and semilocal vr[1-3] potential
243 std::vector<Tensor<double> > vr = xcfunc.vxc(xc_args,ispin);
244
245
246 #if 0
247 print("\n");
248 if(what == 0)
249 print("rhoa \t vr \t read\t comp \n");
250 if(what == 1)
251 print("rhoa \t vsigmaa \t read\t comp \n");
252 if(what == 2)
253 print("rhoa \t vsigmab \t read\t comp \n");
254 if(what == 3)
255 print("rhoa \t zk \t read\t comp \n");
256 for (unsigned int idp = 0; idp < dps.size(); idp++)
257 {
258 if (spin_polarized)
259 {
260 if(what == 0)
261 printf("%25.12e %25.12e %25.12e %25.12e \n",
262 rhoa_t[idp], vrhoa_t[idp], dps[idp].vrhoa, vr[idp] );
263 if(what == 1)
264 printf("%i %25.12e %25.12e %25.12e %25.12e \n",
265 idp, rhoa_t[idp], vsigmaaa_t[idp], dps[idp].vsigmaaa, vr[idp]);
266 if(what == 2)
267 printf("%i %25.12e %25.12e %25.12e %25.12e \n",
268 idp, rhoa_t[idp], vsigmaab_t[idp], dps[idp].vsigmaab, vr[idp]);
269 if(what == 3)
270 printf("%25.12e %25.12e %25.12e %25.12e \n",
271 rhoa_t[idp], zk_t[idp], dps[idp].zk, vr[idp]);
272 }
273 else
274 {
275 if(what == 0)
276 printf("%25.12e %25.12e %25.12e %25.12e \n",
277 rhoa_t[idp], vrhoa_t[idp], dps[idp].vrhoa, vr[idp] );
278 if(what == 1)
279 printf("vsga %i %25.12e %25.12e %25.12e %25.12e \n",
280 idp, rhob_t[idp], vsigmaaa_t[idp], dps[idp].vsigmaaa, vr[idp]);
281 // idp, rhob_t[idp], .5*(vsigmaaa_t[idp] + .5*vsigmaab_t[idp]), dps[idp].vsigmaaa, vr[idp]);
282 if(what == 2)
283 printf (" \n bad bad bad\n ") ;
284 if(what ==3)
285 printf("%25.12e %25.12e %25.12e %25.12e \n",
286 rhoa_t[idp], zk_t[idp], dps[idp].zk, vr[idp]);
287 }
288 }
289 print("\n\n");
290 #endif
291
292 FILE *f=0;
293 f = fopen((df_repo_functionals[istr]+".dat").c_str(), "w");
294
295 //if(xcfunc.is_gga()){
296 if (xcfunc.is_spin_polarized())
297 {
298 fprintf(f,"%25s %25s %25s %25s %25s %25s %25s\n","#rhoa","rhob","sigmaaa",
299 "sigmaab","sigmabb","vrhoa ","vsigmaaa * drho_x");
300 for (unsigned int idp = 0; idp < dps.size(); idp++) {
301 fprintf(f,"%25.12e %25.12e %25.12e %25.12e %25.12e %25.12e %25.12e\n",
302 rhoa_t[idp], rhob_t[idp], sigmaaa_t[idp], sigmaab_t[idp], sigmabb_t[idp],vr[0][idp], vr[1][idp]);
303 }
304 }
305 else
306 {
307 fprintf(f,"%25s %25s %25s %25s\n","#rho","sigma","vrho ","vsigma*drho_x");
308 for (unsigned int idp = 0; idp < dps.size(); idp++)
309 {
310 fprintf(f,"%25.12e %25.12e %25.12e %25.12e\n",
311 rhoa_t[idp], sigmaaa_t[idp], vr[0][idp], vr[1][idp]);
312 }
313 }
314 fprintf(f,"\n\n");
315 /*}
316 else
317 {
318 if (xcfunc.is_spin_polarized())
319 {
320 fprintf(f,"%25s %25s %25s\n","#rhoa","rhob","vrhoa");
321 for (unsigned int idp = 0; idp < dps.size(); idp++)
322 {
323 fprintf(f,"%25.12e %25.12e %25.12e\n",
324 rhoa_t[idp], rhob_t[idp], vr[idp]);
325 }
326 }
327 else
328 {
329 fprintf(f,"%25s %25s\n","#rho","vrho");
330 for (unsigned int idp = 0; idp < dps.size(); idp++)
331 {
332 fprintf(f,"%25.12e %25.12e\n",
333 rhoa_t[idp], vr[idp]);
334 }
335 }
336 fprintf(f,"\n\n");
337 }*/
338
339 fclose(f);
340
341 }
342
343 }
344
main(int argc,char ** argv)345 int main(int argc, char** argv) {
346 madness::initialize(argc, argv);
347
348 madness::World world(SafeMPI::COMM_WORLD);
349 world.gop.fence();
350
351 test_xcfunctional(world);
352
353 madness::finalize();
354 return 0;
355 }
356