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 #include <madness/mra/mra.h>
34 #include <iostream>
35
36 #include "dft.h"
37 #include "hartreefock.h"
38
39 using std::cout;
40 using std::endl;
41
42 using namespace madness;
43
44 const double PI = 3.1415926535897932384;
45
46 typedef Vector<double,3> coordT;
47
48 /// Returns radius for smoothing nuclear potential with energy precision eprec
49 //*****************************************************************************
smoothing_parameter(double Z,double eprec)50 static double smoothing_parameter(double Z, double eprec) {
51 // The min is since asymptotic form not so good at low acc.
52 // The 2 is from two electrons in 1s closed shell.
53 if (Z == 0.0) return 1.0;
54 double Z5 = Z*Z*Z*Z*Z;
55 double c = pow(std::min(1e-3,eprec)/2.0/0.00435/Z5,1.0/3.0);
56 return c;
57 }
58 //*****************************************************************************
59
60
61 /// Regularized 1/r potential.
62
63 /// Invoke as \c u(r/c)/c where \c c is the radius of the
64 /// smoothed volume.
65 //*****************************************************************************
smoothed_potential(double r)66 static double smoothed_potential(double r) {
67 const double THREE_SQRTPI = 5.31736155271654808184;
68 double r2 = r*r, pot;
69 if (r > 6.5){
70 pot = 1.0/r;
71 } else if (r > 1e-8){
72 pot = erf(r)/r + (exp(-r2) + 16.0*exp(-4.0*r2))/(THREE_SQRTPI);
73 } else{
74 pot = (2.0 + 17.0/3.0)/sqrt(PI);
75 }
76
77 return pot;
78 }
79 //*****************************************************************************
80
81 //*****************************************************************************
psi_func_be1(const coordT & rr)82 static double psi_func_be1(const coordT& rr)
83 {
84 const double x=rr[0], y=rr[1], z=rr[2];
85 double r = sqrt(x*x+y*y+z*z);
86 return exp(-4.0*r+1e-4);
87 }
88 //*****************************************************************************
89
90 //*****************************************************************************
psi_func_be2(const coordT & rr)91 static double psi_func_be2(const coordT& rr)
92 {
93 const double x=rr[0], y=rr[1], z=rr[2];
94 double r = sqrt(x*x+y*y+z*z);
95 return (1.0 - 2.0*r*exp(-2.0*r));
96 }
97 //*****************************************************************************
98
99
100 //*****************************************************************************
V_func_be(const coordT & r)101 static double V_func_be(const coordT& r)
102 {
103 const double x=r[0], y=r[1], z=r[2];
104 double rr = sqrt(x*x + y*y + z*z);
105 double c = smoothing_parameter(4.0, 1e-7);
106 return -4.0 * smoothed_potential(rr/c) / c;
107 }
108 //*****************************************************************************
109
110 //*****************************************************************************
rho_func_be(const coordT & rr)111 static double rho_func_be(const coordT& rr)
112 {
113 const double x=rr[0], y=rr[1], z=rr[2];
114 // double e1 = 20.0;
115 // double coeff = pow(e1/PI, 1.5);
116 // return -4.0 * coeff * exp(-e1 * (x*x + y*y + z*z));
117 double c = 0.1;
118 double r = sqrt(x*x + y*y + z*z);
119 r = r / c;
120 const double RPITO1P5 = 0.1795871221251665617; // 1.0/Pi^1.5
121 return 4.0 * ((-3.0/2.0+(1.0/3.0)*r*r)*exp(-r*r)+(-32.0+(256.0/3.0)*r*r)*exp(-4.0*r*r))*RPITO1P5/c/c/c;
122 }
123 //*****************************************************************************
124
125 //*****************************************************************************
test_hf_be(World & world)126 void test_hf_be(World& world)
127 {
128 //if (world.rank() == 0) cout << "Running test application HartreeFock ..." << endl;
129
130 typedef Vector<double,3> coordT;
131 typedef std::shared_ptr< FunctionFunctorInterface<double,3> > functorT;
132
133 // Function defaults
134 double bsize = 42.4;
135 int funck = 6;
136 double thresh = 1e-4;
137 FunctionDefaults<3>::set_k(funck);
138 FunctionDefaults<3>::set_thresh(thresh);
139 FunctionDefaults<3>::set_autorefine(false);
140 FunctionDefaults<3>::set_cubic_cell(-bsize, bsize);
141
142 // Nuclear potential (Be)
143 const coordT origin(0.0);
144 if (world.rank() == 0) madness::print("Creating Function object for nuclear charge density ...");
145 Function<double,3> rhon = FunctionFactory<double,3>(world).f(rho_func_be).thresh(thresh).initial_level(4);
146 Function<double,3> vnuc = FunctionFactory<double,3>(world).f(V_func_be).thresh(thresh);
147
148 double rhontrace = rhon.trace();
149 if (world.rank() == 0)
150 printf("trace of rhon = %.8f\n\n", rhontrace);
151
152 if (world.rank() == 0)
153 cout << "Operating on nuclear charge density ..." << endl;
154 SeparatedConvolution<double, 3> op = CoulombOperator<double, 3> (world,
155 FunctionDefaults<3>::get_k(), 1e-8, thresh);
156 Function<double, 3> V_from_rho_nuc = apply(op, rhon);
157 if (world.rank() == 0) printf("\n");
158 double L = 2.0 * bsize;
159 double bstep = L / 100.0;
160 vnuc.reconstruct();
161 V_from_rho_nuc.reconstruct();
162 for (int i = 0; i < 101; i++)
163 {
164 coordT p(-L / 2 + i * bstep);
165 double error = fabs(vnuc(p) - V_from_rho_nuc(p));
166 if (world.rank() == 0)
167 printf("%.2f\t\t%.8f\t%.8f\t%.8f\t%.8f\n", p[0], vnuc(p), V_from_rho_nuc(p),
168 error, error / vnuc(p));
169 }
170 if (world.rank() == 0) printf("\n");
171
172
173 // Guess for the wavefunctions
174 if (world.rank() == 0) madness::print("Creating wavefunction's ...");
175 Function<double,3> psi1 = FunctionFactory<double,3>(world).f(psi_func_be1);
176 psi1.scale(1.0/psi1.norm2());
177 Function<double,3> psi2 = FunctionFactory<double,3>(world).f(psi_func_be2);
178 psi2.scale(1.0/psi2.norm2());
179
180 // Create list of wavefunctions
181 std::vector<Function<double,3> > phis;
182 phis.push_back(psi1);
183 phis.push_back(psi2);
184
185 // Creat list of eigenvalues
186 std::vector<double> eigs;
187 eigs.push_back(-5.0);
188 eigs.push_back(-0.5);
189
190 // Create DFT object
191 if (world.rank() == 0) madness::print("Creating DFT object...");
192 //HartreeFock hf(world, Vnuc, phis, eigs, true, true, thresh);
193 ElectronicStructureParams params;
194 params.periodic = false;
195 DFT<double,3> dftcalc(world, rhon, phis, eigs, params);
196 if (world.rank() == 0) madness::print("Running DFT object...");
197 dftcalc.solve(51);
198
199 }
200 //*****************************************************************************
201
202 #define TO_STRING(s) TO_STRING2(s)
203 #define TO_STRING2(s) #s
204
205 //*****************************************************************************
main(int argc,char ** argv)206 int main(int argc, char** argv)
207 {
208 SafeMPI::Init(argc, argv);
209 World world(SafeMPI::COMM_WORLD);
210 if (world.rank() == 0)
211 {
212 print("");
213 print("--------------------------------------------");
214 print(" MADNESS", " multiresolution testsuite");
215 print("--------------------------------------------");
216 print("");
217 print(" number of processors ...", world.size());
218 print(" processor frequency ...", cpu_frequency());
219 print(" host system ...", TO_STRING(HOST_SYSTEM));
220 print(" byte order ...", TO_STRING(MADNESS_BYTE_ORDER));
221 print(" configured by ...", MADNESS_CONFIGURATION_USER);
222 print(" configured on ...", MADNESS_CONFIGURATION_HOST);
223 print(" configured at ...", MADNESS_CONFIGURATION_DATE);
224 print(" CXX ...", MADNESS_CONFIGURATION_CXX);
225 print(" CXXFLAGS ...", MADNESS_CONFIGURATION_CXXFLAGS);
226 #ifdef OPTERON_TUNE
227 print(" tuning for ...", "opteron");
228 #elif defined(CORE_DUO_TUNE)
229 print(" tuning for ...", "core duo");
230 #else
231 print(" tuning for ...", "default");
232 #endif
233 #ifdef BOUNDS_CHECKING
234 print(" tensor bounds checking ...", "enabled");
235 #endif
236 #ifdef TENSOR_INSTANCE_COUNT
237 print(" tensor instance count ...", "enabled");
238 #endif
239 print(" ");
240 }
241
242 try
243 {
244 printf("WSTHORNTON: Starting up the world ... \n");
245
246 startup(world,argc,argv);
247 if (world.rank() == 0) print("Initial tensor instance count", BaseTensor::get_instance_count());
248 test_hf_be(world);
249 }
250 catch (const SafeMPI::Exception& e)
251 {
252 print(e);
253 error("caught an MPI exception");
254 }
255 catch (const madness::MadnessException& e)
256 {
257 print(e);
258 error("caught a MADNESS exception");
259 }
260 catch (const madness::TensorException& e)
261 {
262 print(e);
263 error("caught a Tensor exception");
264 }
265 catch (const char* s)
266 {
267 print(s);
268 error("caught a string exception");
269 }
270 catch (const std::string& s)
271 {
272 print(s);
273 error("caught a string (class) exception");
274 }
275 catch (const std::exception& e)
276 {
277 print(e.what());
278 error("caught an STL exception");
279 }
280 catch (...)
281 {
282 error("caught unhandled exception");
283 }
284
285 if (world.rank() == 0)
286 print("entering final fence");
287 world.gop.fence();
288 if (world.rank() == 0)
289 print("done with final fence");
290 if (world.rank() == 0)
291 print("Final tensor instance count", BaseTensor::get_instance_count());
292 SafeMPI::Finalize();
293
294 return 0;
295 }
296 //*****************************************************************************
297
298