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 while(!fstr.eof())
88 {
89 xcfunc_data_point dp;
90 fstr >> tmpstr; fstr >> dp.rhoa;
91 fstr >> tmpstr; fstr >> dp.rhob;
92 fstr >> tmpstr; fstr >> dp.sigmaaa;
93 fstr >> tmpstr; fstr >> dp.sigmaab;
94 fstr >> tmpstr; fstr >> dp.sigmabb;
95
96 fstr >> tmpstr; fstr >> tmpstr; fstr >> dp.zk;
97 fstr >> tmpstr; fstr >> tmpstr; fstr >> dp.vrhoa;
98 fstr >> tmpstr; fstr >> tmpstr; fstr >> dp.vrhob;
99 fstr >> tmpstr; fstr >> tmpstr; fstr >> dp.vsigmaaa;
100 fstr >> tmpstr; fstr >> tmpstr; fstr >> dp.vsigmaab;
101 fstr >> tmpstr; fstr >> tmpstr; fstr >> dp.vsigmabb;
102 fstr >> tmpstr; fstr >> tmpstr; fstr >> dp.v2rhoa2;
103 fstr >> tmpstr; fstr >> tmpstr; fstr >> dp.v2rhoab;
104 fstr >> tmpstr; fstr >> tmpstr; fstr >> dp.v2rhob2;
105
106 // skip unwanted lines for now
107 for (int iskip = 0; iskip < 15; iskip++)
108 fstr.getline(buffer,120);
109
110 if (!spin_polarized)
111 {
112 if (std::abs(dp.vrhoa-dp.vrhob) <= 1e-10)
113 dps.push_back(dp);
114 }
115 else
116 {
117 dps.push_back(dp);
118 }
119 }
120
121 return dps;
122
123 }
124
test_xcfunctional()125 void test_xcfunctional()
126 {
127 bool spin_polarized = false;
128
129 for (int istr = 6; istr < 22; istr++)
130 {
131 XCfunctional xcfunc;
132 std::string xcfuncstr = df_repo_functionals[istr];
133 std::cout << "Testing exchange-correlation functional: "<< xcfuncstr << std::endl;
134
135 xcfuncstr += " 1.0";
136 xcfunc.initialize(xcfuncstr,spin_polarized);
137
138 std::string fpath("df_repo/");
139 fpath += df_repo_functionals[istr];
140 fpath += ".data";
141 std::vector<xcfunc_data_point> dps = read_test_data(fpath.c_str(),spin_polarized);
142
143 Tensor<double> rhoa_t((long)dps.size());
144 Tensor<double> rhob_t((long)dps.size());
145 Tensor<double> sigmaaa_t((long)dps.size());
146 Tensor<double> sigmaab_t((long)dps.size());
147 Tensor<double> sigmabb_t((long)dps.size());
148 std::vector<Tensor<double> > xc_args;
149 for (unsigned int idp = 0; idp < dps.size(); idp++)
150 {
151 rhoa_t(idp) = dps[idp].rhoa;
152 rhob_t(idp) = dps[idp].rhob;
153 sigmaaa_t(idp) = dps[idp].sigmaaa;
154 sigmaab_t(idp) = dps[idp].sigmaab;
155 sigmabb_t(idp) = dps[idp].sigmabb;
156 }
157 if (spin_polarized)
158 {
159 xc_args.push_back(rhoa_t);
160 xc_args.push_back(rhob_t);
161 xc_args.push_back(sigmaaa_t);
162 xc_args.push_back(sigmaab_t);
163 xc_args.push_back(sigmabb_t);
164 }
165 else
166 {
167 xc_args.push_back(rhoa_t);
168 xc_args.push_back(sigmaaa_t);
169 }
170 Tensor<double> vr = xcfunc.vxc(xc_args, 0, 0);
171
172 for (unsigned int idp = 0; idp < dps.size(); idp++)
173 {
174 printf("%25.12e %25.12e %25.12e %25.12e %25.12e\n",
175 rhoa_t[idp], sigmaaa_t[idp], dps[idp].vrhoa, vr[idp],
176 std::abs(dps[idp].vrhoa - vr[idp]));
177 }
178 print("\n\n");
179 }
180
181 spin_polarized = true;
182
183 std::cout << "Testing spin-polarized case: " << std::endl << std::endl;
184
185 for (int istr = 6; istr < 22; istr++)
186 {
187 XCfunctional xcfunc;
188 std::string xcfuncstr = df_repo_functionals[istr];
189 std::cout << "Testing exchange-correlation functional: "<< xcfuncstr << std::endl;
190
191 xcfuncstr += " 1.0";
192 xcfunc.initialize(xcfuncstr,spin_polarized);
193
194 std::string fpath("df_repo/");
195 fpath += df_repo_functionals[istr];
196 fpath += ".data";
197 std::vector<xcfunc_data_point> dps = read_test_data(fpath.c_str(),spin_polarized);
198
199 Tensor<double> rhoa_t((long)dps.size());
200 Tensor<double> rhob_t((long)dps.size());
201 Tensor<double> sigmaaa_t((long)dps.size());
202 Tensor<double> sigmaab_t((long)dps.size());
203 Tensor<double> sigmabb_t((long)dps.size());
204 std::vector<Tensor<double> > xc_args;
205 for (unsigned int idp = 0; idp < dps.size(); idp++)
206 {
207 rhoa_t(idp) = dps[idp].rhoa;
208 rhob_t(idp) = dps[idp].rhob;
209 sigmaaa_t(idp) = dps[idp].sigmaaa;
210 sigmaab_t(idp) = dps[idp].sigmaab;
211 sigmabb_t(idp) = dps[idp].sigmabb;
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 {
219 xc_args.push_back(sigmaaa_t);
220 xc_args.push_back(sigmaab_t);
221 xc_args.push_back(sigmabb_t);
222 }
223 }
224 else
225 {
226 xc_args.push_back(rhoa_t);
227 if (xcfunc.is_gga())
228 {
229 xc_args.push_back(sigmaaa_t);
230 }
231 }
232 Tensor<double> vr = xcfunc.vxc(xc_args,0, 0);
233
234 if (xcfunc.is_spin_polarized())
235 {
236 for (unsigned int idp = 0; idp < dps.size(); idp++)
237 {
238 printf("%25.12e %25.12e %25.12e %25.12e %25.12e %25.12e %25.12e %25.12e\n",
239 rhoa_t[idp], rhob_t[idp], sigmaaa_t[idp], sigmaab_t[idp], sigmabb_t[idp],
240 dps[idp].vrhoa, vr[idp],
241 std::abs(dps[idp].vrhoa - vr[idp]));
242 }
243 }
244 else
245 {
246 for (unsigned int idp = 0; idp < dps.size(); idp++)
247 {
248 printf("%25.12e %25.12e %25.12e %25.12e %25.12e\n",
249 rhoa_t[idp], sigmaaa_t[idp], dps[idp].vrhoa, vr[idp],
250 std::abs(dps[idp].vrhoa - vr[idp]));
251 }
252 }
253 print("\n\n");
254 }
255
256 }
257
main(int argc,char ** argv)258 int main(int argc, char** argv) {
259 madness::initialize(argc, argv);
260
261 madness::World world(SafeMPI::COMM_WORLD);
262 world.gop.fence();
263
264 #ifdef MADNESS_HAS_LIBXC
265 test_xcfunctional();
266 #else
267 std::cout << "WARNING: To run this program you need libXC. All tests will be skipped." << std::endl;
268 #endif
269
270 madness::finalize();
271 return 0;
272 }
273