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
32 //#define WORLD_INSTANTIATE_STATIC_TEMPLATES
33 #include <madness/world/MADworld.h>
34 #include <madness/mra/mra.h>
35 #include <madness/tensor/tensor.h>
36 #include <fstream>
37 #include "xcfunctional.h"
38
39 using namespace madness;
40
41 static std::string df_repo_functionals[] = {
42 "lda_x",
43 "lda_c_vwn_rpa",
44 "lda_c_vwn",
45 "lda_c_pz",
46 "lda_c_pw",
47 "hyb_gga_xc_b3lyp",
48 "gga_xc_hcth_93",
49 "gga_xc_hcth_407",
50 "gga_xc_hcth_147",
51 "gga_xc_hcth_120",
52 "gga_xc_edf1",
53 "gga_xc_b97_2",
54 "gga_xc_b97_1",
55 "gga_xc_b97",
56 "gga_x_pw91",
57 "gga_x_pbe",
58 "gga_x_ft97_b",
59 "gga_x_b88",
60 "gga_c_pw91",
61 "gga_c_pbe",
62 "gga_c_p86",
63 "gga_c_lyp"};
64
65
66 struct xcfunc_data_point
67 {
68 double rhoa, rhob;
69 double sigmaaa, sigmaab, sigmabb;
70 double zk;
71 double vrhoa, vrhob;
72 double vsigmaaa, vsigmaab, vsigmabb;
73 double v2rhoa2, v2rhoab, v2rhob2;
74 };
75
read_test_data(const std::string & dfname,bool spin_polarized)76 std::vector<xcfunc_data_point> read_test_data(const std::string& dfname,
77 bool spin_polarized)
78 {
79 std::ifstream fstr(dfname.c_str());
80 char buffer[120];
81 fstr.getline(buffer, 120);
82 fstr.getline(buffer, 120);
83
84 std::string tmpstr;
85
86 std::vector<xcfunc_data_point> dps;
87 print(dfname.c_str());
88 int i=1;
89 while(!fstr.eof())
90 {
91 xcfunc_data_point dp;
92
93 fstr >> tmpstr; fstr >> dp.rhoa;
94 fstr >> tmpstr; fstr >> dp.rhob;
95 fstr >> tmpstr; fstr >> dp.sigmaaa;
96 fstr >> tmpstr; fstr >> dp.sigmaab;
97 fstr >> tmpstr; fstr >> dp.sigmabb;
98
99
100
101 fstr >> tmpstr; fstr >> tmpstr; fstr >> dp.zk;
102 fstr >> tmpstr; fstr >> tmpstr; fstr >> dp.vrhoa;
103 fstr >> tmpstr; fstr >> tmpstr; fstr >> dp.vrhob;
104 fstr >> tmpstr; fstr >> tmpstr; fstr >> dp.vsigmaaa;
105 fstr >> tmpstr; fstr >> tmpstr; fstr >> dp.vsigmaab;
106 fstr >> tmpstr; fstr >> tmpstr; fstr >> dp.vsigmabb;
107 fstr >> tmpstr; fstr >> tmpstr; fstr >> dp.v2rhoa2;
108 fstr >> tmpstr; fstr >> tmpstr; fstr >> dp.v2rhoab;
109 fstr >> tmpstr; fstr >> tmpstr; fstr >> dp.v2rhob2;
110
111
112 for (int iskip = 0; iskip < 13; iskip++)
113 fstr.getline(buffer,120);
114 if (!spin_polarized)
115 {
116 if (std::abs(dp.rhoa-dp.rhob) <= 1e-10)
117 dps.push_back(dp);
118 }
119 else
120 {
121 dps.push_back(dp);
122 }
123 i=i+1;
124 if(i == 26) break;
125 }
126
127 return dps;
128
129 }
130
test_xcfunctional(World & world)131 void test_xcfunctional(World& world)
132 {
133 bool spin_polarized ;
134 //spin_polarized = true;
135 spin_polarized = false;
136
137 int ispin = 0; //alpha=0 beta=1
138 int istr = 0;
139
140 //vama5 0 "lda_x",
141 //vama5 1 "lda_c_vwn_rpa",
142 //vama5 2 "lda_c_vwn",
143 //vama5 3 "lda_c_pz",
144 //vama5 4 "lda_c_pw",
145 //vama5 5 "hyb_gga_xc_b3lyp",
146 //vama5 6 "gga_xc_hcth_93",
147 //vama5 7 "gga_xc_hcth_407",
148 //vama5 8 "gga_xc_hcth_147",
149 //vama5 9 "gga_xc_hcth_120",
150 //vama5 10 "gga_xc_edf1",
151 //vama5 11 "gga_xc_b97_2",
152 //vama5 12 "gga_xc_b97_1",
153 //vama5 13 "gga_xc_b97",
154 //vama5 14 "gga_x_pw91",
155 //vama5 15 "gga_x_pbe",
156
157 //vama5 16 "gga_x_ft97_b",
158 //vama5 17 "gga_x_b88",
159 //vama5 18 "gga_c_pw91",
160 //vama5 19 "gga_c_pbe",
161
162 //vama5 20 "gga_c_p86",
163 //vama5 21 "gga_c_lyp"};
164
165 XCfunctional xcfunc;
166 std::string xcfuncstr = df_repo_functionals[istr];
167 std::cout << "Testing exchange-correlation functional: "<< xcfuncstr << std::endl;
168
169 xcfuncstr += " 1.0";
170 xcfunc.initialize(xcfuncstr,spin_polarized,world);
171
172 std::string fpath("df_repo/");
173 fpath += df_repo_functionals[istr];
174 fpath += ".data";
175 std::vector<xcfunc_data_point> dps = read_test_data(fpath.c_str(),spin_polarized);
176
177 print("hola");
178 Tensor<double> rhoa_t((long)dps.size());
179 Tensor<double> rhob_t((long)dps.size());
180 Tensor<double> sigmaaa_t((long)dps.size());
181 Tensor<double> sigmaab_t((long)dps.size());
182 Tensor<double> sigmabb_t((long)dps.size());
183 Tensor<double> zk_t((long)dps.size());
184 Tensor<double> vrhoa_t((long)dps.size());
185 Tensor<double> vrhob_t((long)dps.size());
186 Tensor<double> vsigmaaa_t((long)dps.size());
187 Tensor<double> vsigmaab_t((long)dps.size());
188 Tensor<double> vsigmabb_t((long)dps.size());
189
190 Tensor<double> rhoa1_t((long)dps.size());
191 Tensor<double> sigmaaa1_t((long)dps.size());
192
193 std::vector<Tensor<double> > xc_args;
194 for (unsigned int idp = 0; idp < dps.size(); idp++)
195 {
196 rhoa_t(idp) = dps[idp].rhoa;
197 rhob_t(idp) = dps[idp].rhob;
198
199 sigmaaa_t(idp) = dps[idp].sigmaaa;
200 sigmaab_t(idp) = dps[idp].sigmaab;
201 sigmabb_t(idp) = dps[idp].sigmabb;
202
203 zk_t(idp) = dps[idp].zk;
204
205 vrhoa_t(idp) = dps[idp].vrhoa;
206 vrhob_t(idp) = dps[idp].vrhob;
207
208
209 vsigmaaa_t(idp) = dps[idp].vsigmaaa;
210 vsigmaab_t(idp) = dps[idp].vsigmaab;
211 vsigmabb_t(idp) = dps[idp].vsigmabb;
212 }
213 if (spin_polarized)
214 {
215 xc_args.push_back(rhoa_t);
216 xc_args.push_back(rhob_t);
217 if (xcfunc.is_gga()) {
218 xc_args.push_back(sigmaaa_t);
219 xc_args.push_back(sigmaab_t);
220 xc_args.push_back(sigmabb_t);
221 }
222 }
223 else
224 {
225 xc_args.push_back(rhoa_t);
226 if (xcfunc.is_gga()) {
227 xc_args.push_back(sigmaaa_t);
228 }
229 }
230
231 print("xc_args_size", xc_args.size());
232 print("ispin ", ispin);
233 print("spin polarized ", spin_polarized);
234 //vama1 std::cout << "Testing spin-polarized case: " << std::endl << std::endl;
235
236 // xc local vr[0] and semilocal vr[1-3] potential
237 std::vector<Tensor<double> > vr = xcfunc.vxc(xc_args,ispin);
238
239 #if 0
240 print("\n");
241 if(what == 0)
242 print("rhoa \t vr \t read\t comp \n");
243 if(what == 1)
244 print("rhoa \t vsigmaa \t read\t comp \n");
245 if(what == 2)
246 print("rhoa \t vsigmab \t read\t comp \n");
247 if(what == 3)
248 print("rhoa \t zk \t read\t comp \n");
249 for (unsigned int idp = 0; idp < dps.size(); idp++)
250 {
251 if (spin_polarized)
252 {
253 if(what == 0)
254 printf("%25.12e %25.12e %25.12e %25.12e \n",
255 rhoa_t[idp], vrhoa_t[idp], dps[idp].vrhoa, vr[idp] );
256 if(what == 1)
257 printf("%i %25.12e %25.12e %25.12e %25.12e \n",
258 idp, rhoa_t[idp], vsigmaaa_t[idp], dps[idp].vsigmaaa, vr[idp]);
259 if(what == 2)
260 printf("%i %25.12e %25.12e %25.12e %25.12e \n",
261 idp, rhoa_t[idp], vsigmaab_t[idp], dps[idp].vsigmaab, vr[idp]);
262 if(what == 3)
263 printf("%25.12e %25.12e %25.12e %25.12e \n",
264 rhoa_t[idp], zk_t[idp], dps[idp].zk, vr[idp]);
265 }
266 else
267 {
268 if(what == 0)
269 printf("%25.12e %25.12e %25.12e %25.12e \n",
270 rhoa_t[idp], vrhoa_t[idp], dps[idp].vrhoa, vr[idp] );
271 if(what == 1)
272 printf("vsga %i %25.12e %25.12e %25.12e %25.12e \n",
273 idp, rhob_t[idp], vsigmaaa_t[idp], dps[idp].vsigmaaa, vr[idp]);
274 // idp, rhob_t[idp], .5*(vsigmaaa_t[idp] + .5*vsigmaab_t[idp]), dps[idp].vsigmaaa, vr[idp]);
275 if(what == 2)
276 printf (" \n bad bad bad\n ") ;
277 if(what ==3)
278 printf("%25.12e %25.12e %25.12e %25.12e \n",
279 rhoa_t[idp], zk_t[idp], dps[idp].zk, vr[idp]);
280 }
281 }
282 print("\n\n");
283 #endif
284
285 if (xcfunc.is_spin_polarized())
286 {
287 printf("%25s %25s %25s %25s %25s %25s %25s %25s\n","#rhoa","rhob","sigmaaa","sigmaab","sigmabb","vrhoa (input)","vr (output)","vrhoa-vr");
288 for (unsigned int idp = 0; idp < dps.size(); idp++)
289 {
290 printf("%25.12e %25.12e %25.12e %25.12e %25.12e %25.12e %25.12e %25.12e\n",
291 rhoa_t[idp], rhob_t[idp], sigmaaa_t[idp], sigmaab_t[idp], sigmabb_t[idp],
292 dps[idp].vrhoa, vr[0][idp],
293 std::abs(dps[idp].vrhoa - vr[0][idp]));
294 }
295 }
296 else
297 {
298 printf("%25s %25s %25s %25s %25s\n","#rhoa","sigmaaa","vrhoa (input)","vr (output)","vrhoa-vr");
299 for (unsigned int idp = 0; idp < dps.size(); idp++)
300 {
301 printf("%25.12e %25.12e %25.12e %25.12e %25.12e\n",
302 rhoa_t[idp], sigmaaa_t[idp], dps[idp].vrhoa, vr[0][idp],
303 std::abs(dps[idp].vrhoa - vr[0][idp]));
304 }
305 }
306 print("\n\n");
307
308 }
309
main(int argc,char ** argv)310 int main(int argc, char** argv) {
311 madness::initialize(argc, argv);
312
313 madness::World world(SafeMPI::COMM_WORLD);
314 world.gop.fence();
315
316 test_xcfunctional(world);
317
318 madness::finalize();
319 return 0;
320 }
321