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
34 //// \file examples/vnucso.cc
35 /// \brief Solves the spin-orbit nuclear potential problem
36
37 /*!
38
39 \file vnucso.cc
40 \brief Solves the Hartree-Fock equation for the 2-cosh potential with spin-orbit in Nuclear
41 Density Functional Theory witough assumption on spatial symmetry.
42 \ingroup examples
43
44 Points of interest
45
46 - Forming a Hamiltonian and a Fock matrix
47 - Forming and solving a generalized eigensystem problem using LAPACK
48 - Refining the representation of real and complex functions
49 - Vectors of functions and operators, inner-product, gaxpy
50 - Application of the Helmholtz bound-state Green function as vector of operators and functions
51 - Projection and change of representation of functions from multiwavelets of degree k to k+1
52
53 The source is <a href=http://code.google.com/p/m-a-d-n-e-s-s/source/browse/local/trunk/src/apps/examples/vnucso.cc>here</a>.
54
55
56 This is a more involved example than the Hydrogen and Helium.
57 This example couples the traditional diagonalization approach with that of the integral equation approach.
58 The details are described in:
59 G. I. Fann , J. Pei, R. J. Harrison1, J. Jia, J. Hill1 , M. Ou, W. Nazarewicz, W. A. Shelton
60 and N. Schunck, "Fast multiresolution methods for density functional theory in nuclear physics,"
61 Journal of Physics, 180 (2009) 012080.
62
63
64 */
65
66
67 //#define WORLD_INSTANTIATE_STATIC_TEMPLATES
68 #include <madness/mra/mra.h>
69 #include <madness/mra/vmra.h>
70 #include <madness/mra/operator.h>
71 #include <madness/constants.h>
72
73 using namespace madness;
74 // using namespace std;
75
76 typedef Vector<double,3> coordT;
77 typedef std::shared_ptr< FunctionFunctorInterface<double,3> > real_functorT;
78 typedef std::shared_ptr< FunctionFunctorInterface<double_complex,3> > complex_functorT;
79 typedef Function<double,3> real_functionT;
80 typedef Function<double_complex,3> complex_functionT;
81 typedef FunctionFactory<double,3> real_factoryT;
82 typedef FunctionFactory<double_complex,3> complex_factoryT;
83 typedef SeparatedConvolution<double,3> operatorT;
84 typedef std::shared_ptr<operatorT> poperatorT;
85
86 const double L = 120.0; // box is [-L,L]
87 const double zeta = 7.5; // potential wells at +/-zeta
88 const double R1 = 2.0; // potential parameter
89 const double R2 = 2.0; // potential parameter
90 const double a1 = 1.0; // potential parameter
91 const double a2 = 1.0; // potential parameter
92 const double reduced = 0.04825964488415279478;
93 const double V1 = -50.0*reduced; // potential parameter
94 const double V2 = -50.0*reduced; // potential parameter
95 const double lambda_correct = 0.0026608048208104861/reduced; // SO potential parameter
96
97 //static const double lambda_fudge = lambda_correct*100.0;
98
99 static double lambda = lambda_correct;
100 static const double fac1=exp(-R1/a1);
101 static const double fac2=exp(-R2/a2);
102
103 struct Guess : FunctionFunctorInterface<double_complex,3> {
104 const double z; // z-coordinate center
105 const double exponent; // exponent
106 const int nx, ny, nz; // powers of x, y, z ... only 0 or 1 supported!
107
GuessGuess108 Guess(double z, double exponent, int nx, int ny, int nz)
109 : z(z), exponent(exponent), nx(nx), ny(ny), nz(nz) {}
110
operator ()Guess111 double_complex operator()(const coordT& r) const {
112 double rsq = r[0]*r[0] + r[1]*r[1] + (r[2]-z)*(r[2]-z);
113 double psi = exp(-exponent*rsq);
114 if (nx) psi *= r[0];
115 if (ny) psi *= r[1];
116 if (nz) psi *= (r[2]-z);
117 return double_complex(psi,0.0);
118 }
119 };
120
V(const coordT & r)121 static double V(const coordT& r)
122 {
123 const double x=r[0], y=r[1], z=r[2];
124 const double zp=(z+zeta), zm = (z-zeta);
125 const double rp = sqrt(x*x + y*y + zp*zp);
126 const double rm = sqrt(x*x + y*y + zm*zm);
127
128 return V1/(1.0 + fac1*cosh(rp/a1)) + V2/(1.0 + fac2*cosh(rm/a2));
129 }
130
moments(World & world,const std::vector<complex_functionT> & u,const std::vector<complex_functionT> & v)131 void moments(World& world, const std::vector<complex_functionT>& u, const std::vector<complex_functionT>& v) {
132 FunctionDefaults<3>::set_autorefine(true);
133 complex_functionT rho = complex_factoryT(world);
134 rho.compress();
135 reconstruct(world, u);
136 reconstruct(world, v);
137 for (unsigned int i=0; i<u.size(); i++) {
138 complex_functionT psisq = conj(u[i])*u[i] + conj(v[i])*v[i];
139 psisq.compress();
140 rho.gaxpy(1.0, psisq, 1.0);
141 }
142 double_complex moment1 = rho.trace();
143 complex_functionT rhosq = (rho*rho).scale(2.);
144 double_complex moment2a = rhosq.trace();
145 double_complex moment2b = inner(rho,rho);
146 double_complex moment3 = inner(rho,rhosq);
147 if (world.rank() == 0) {
148 print(" moment1", moment1);
149 print(" moment2", moment2a, moment2b);
150 print(" moment3", moment3);
151 }
152 FunctionDefaults<3>::set_autorefine(false);
153 }
154
gaxpy1(World & world,const double_complex alpha,std::vector<complex_functionT> & a,const double_complex beta,const std::vector<complex_functionT> & b,bool fence=true)155 void gaxpy1(World& world,
156 const double_complex alpha,
157 std::vector<complex_functionT>& a,
158 const double_complex beta,
159 const std::vector<complex_functionT>& b,
160 bool fence=true)
161 {
162 MADNESS_ASSERT(a.size() == b.size());
163
164 for (unsigned int i=0; i<a.size(); i++) {
165 a[i] = alpha*a[i] + beta*b[i];
166 }
167 if (fence) world.gop.fence();
168 }
169
make_bsh_operators(World & world,const Tensor<double> & evals,double tol)170 std::vector<poperatorT> make_bsh_operators(World& world, const Tensor<double>& evals, double tol)
171 {
172 int n = evals.dim(0);
173 std::vector<poperatorT> ops(n);
174 for (int i=0; i<n; i++) {
175 double eps = evals(i);
176 if (eps > 0) eps = -0.05;
177 double lo = 0.1*tol; // heuristic
178 ops[i] = poperatorT(BSHOperatorPtr3D(world, sqrt(-eps), lo, tol));
179 }
180 return ops;
181 }
182
hamiltonian_matrix(World & world,const std::vector<complex_functionT> & u,const std::vector<complex_functionT> & v,const std::vector<complex_functionT> & Vu,const std::vector<complex_functionT> & Vv,const std::vector<complex_functionT> du[3],const std::vector<complex_functionT> dv[3])183 Tensor<double_complex> hamiltonian_matrix(World& world,
184 const std::vector<complex_functionT>& u,
185 const std::vector<complex_functionT>& v,
186 const std::vector<complex_functionT>& Vu,
187 const std::vector<complex_functionT>& Vv,
188 const std::vector<complex_functionT> du[3],
189 const std::vector<complex_functionT> dv[3])
190 {
191 reconstruct(world, u);
192 reconstruct(world, v);
193 int n = u.size();
194 Tensor<double_complex> r(n,n);
195
196 for (int axis=0; axis<3; axis++) {
197 r += matrix_inner(world, du[axis], du[axis], true);
198 r += matrix_inner(world, dv[axis], dv[axis], true);
199 }
200 return r + matrix_inner(world, u, Vu, true) + matrix_inner(world, v, Vv, true);
201 }
202
apply_potential(World & world,const real_functionT & V0,const real_functionT & V0x,const real_functionT & V0y,const real_functionT & V0z,const std::vector<complex_functionT> & u,const std::vector<complex_functionT> & v,std::vector<complex_functionT> & Vu,std::vector<complex_functionT> & Vv,std::vector<complex_functionT> du[3],std::vector<complex_functionT> dv[3],bool doso)203 void apply_potential(World& world,
204 const real_functionT& V0,
205 const real_functionT& V0x,
206 const real_functionT& V0y,
207 const real_functionT& V0z,
208 const std::vector<complex_functionT>& u,
209 const std::vector<complex_functionT>& v,
210 std::vector<complex_functionT>& Vu,
211 std::vector<complex_functionT>& Vv,
212 std::vector<complex_functionT> du[3],
213 std::vector<complex_functionT> dv[3],
214 bool doso)
215 {
216 const double_complex lam(lambda,0.0);
217 const double_complex one(1.0,0.0);
218 const double_complex I(0.0,1.0);
219 const int x=0, y=1, z=2; // Attempt to make SO term more readable
220
221 reconstruct(world, u);
222 reconstruct(world, v);
223 V0.reconstruct();
224 V0x.reconstruct();
225 V0y.reconstruct();
226 V0z.reconstruct();
227
228 for (int axis=0; axis<3; axis++) {
229 complex_derivative_3d D = free_space_derivative<double_complex,3>(world, axis);
230 du[axis] = apply(world, D, u);
231 dv[axis] = apply(world, D, v);
232 }
233
234 Vu = mul(world, V0, u);
235 Vv = mul(world, V0, v );
236
237 if (!doso) return;
238
239 gaxpy(world, one, Vu, -I*lam, mul(world, V0y, du[x]));
240 gaxpy(world, one, Vu, I*lam, mul(world, V0x, du[y]));
241 gaxpy(world, one, Vu, lam, mul(world, V0z, dv[x]));
242 gaxpy(world, one, Vu, -I*lam, mul(world, V0z, dv[y]));
243 gaxpy(world, one, Vu, -lam, mul(world, V0x, dv[z]));
244 gaxpy(world, one, Vu, I*lam, mul(world, V0y, dv[z]));
245
246
247 gaxpy(world, one, Vv, I*lam, mul(world, V0y, dv[x]));
248 gaxpy(world, one, Vv, -I*lam, mul(world, V0x, dv[y]));
249 gaxpy(world, one, Vv, -lam, mul(world, V0z, du[x]));
250 gaxpy(world, one, Vv, -I*lam, mul(world, V0z, du[y]));
251 gaxpy(world, one, Vv, lam, mul(world, V0x, du[z]));
252 gaxpy(world, one, Vv, I*lam, mul(world, V0y, du[z]));
253
254 }
255
normalize2(World & world,std::vector<complex_functionT> & u,std::vector<complex_functionT> & v)256 void normalize2(World& world, std::vector<complex_functionT>& u, std::vector<complex_functionT>& v) {
257 std::vector<double> unorm = norm2s(world,u);
258 std::vector<double> vnorm = norm2s(world,v);
259 std::vector<double> normu(u.size());
260
261 for (unsigned int i=0; i<u.size(); i++) {
262 normu[i] = sqrt(unorm[i]*unorm[i] + vnorm[i]*vnorm[i]);
263 }
264 for (unsigned int i=0; i<u.size(); i++) {
265 u[i].scale(double_complex(1.0/normu[i],0.));
266 v[i].scale(double_complex(1.0/normu[i],0.));
267 }
268 }
269
270
271
doit(World & world)272 void doit(World& world) {
273 // We are not using time-reversal symmetry and hence are doing
274 // about 2x too much work. With a little bit more complexity we
275 // could readily exploit it. Time reversal symmetry implies that
276 // if (u,v) (two-component wave function) is an eigen function
277 // then (-v*, u*) is a degenerate solution. We could use this
278 // symmetry by ensuring that after updating (and before
279 // diagonalization) the higher energy states are orthogonal to
280 // lower states and their time-reversed form.
281
282 if (world.rank() == 0) print("entered solver at time", wall_time());
283 long k = 5; // wavelet order
284 double thresh = 1e-3; // precision for wave function
285 coordT zero;
286
287 zero[0]=0.; zero[1]=0.; zero[2]=0.;
288
289 // FunctionDefaults<3>::k = k;
290 // FunctionDefaults<3>::thresh = thresh;
291 // FunctionDefaults<3>::refine = true;
292 // FunctionDefaults<3>::autorefine = false;
293 // FunctionDefaults<3>::initial_level = 2;
294 // FunctionDefaults<3>::truncate_mode = 1;
295
296 FunctionDefaults<3>::set_k(k);
297 FunctionDefaults<3>::set_thresh(thresh);
298 FunctionDefaults<3>::set_refine(true);
299 FunctionDefaults<3>::set_autorefine(false);
300 FunctionDefaults<3>::set_initial_level(2);
301 FunctionDefaults<3>::set_truncate_mode(1);
302 FunctionDefaults<3>::set_cubic_cell(-L,L);
303
304 // for (int i=0; i<3; i++) {
305 // FunctionDefaults<3>::cell(i,0) = -L;
306 // FunctionDefaults<3>::cell(i,1) = L;
307 // }
308 if (world.rank() == 0) print("Making guesses");
309 std::vector<complex_functionT> u;
310 std::vector<complex_functionT> v;
311 std::vector<complex_functionT> w;
312
313 // Initial guess is as for SO-free but first with (u,0) and then (0,v)
314
315 // Could this get any more verbose? This needs to be fixed.
316 u.push_back(complex_factoryT(world).functor(complex_functorT(new Guess( zeta, 0.15, 0, 0, 0))).nofence()); // s
317 v.push_back(complex_functionT(complex_factoryT(world)));
318
319 u.push_back(complex_factoryT(world).functor(complex_functorT(new Guess(-zeta, 0.15, 0, 0, 0))).nofence()); // s
320 v.push_back(complex_functionT(complex_factoryT(world)));
321
322 u.push_back(complex_factoryT(world).functor(complex_functorT(new Guess( zeta, 0.15, 0, 0, 1))).nofence()); // pz
323 v.push_back(complex_functionT(complex_factoryT(world)));
324
325 u.push_back(complex_factoryT(world).functor(complex_functorT(new Guess(-zeta, 0.15, 0, 0, 1))).nofence()); // pz
326 v.push_back(complex_functionT(complex_factoryT(world)));
327 u.push_back(complex_factoryT(world).functor(complex_functorT(new Guess( zeta, 0.15, 0, 1, 0))).nofence()); // py
328 v.push_back(complex_functionT(complex_factoryT(world)));
329 u.push_back(complex_factoryT(world).functor(complex_functorT(new Guess(-zeta, 0.15, 0, 1, 0))).nofence()); // py
330 v.push_back(complex_functionT(complex_factoryT(world)));
331 u.push_back(complex_factoryT(world).functor(complex_functorT(new Guess( zeta, 0.15, 1, 0, 0))).nofence()); // px
332 v.push_back(complex_functionT(complex_factoryT(world)));
333
334 u.push_back(complex_factoryT(world).functor(complex_functorT(new Guess(-zeta, 0.15, 1, 0, 0))).nofence()); // px
335 v.push_back(complex_functionT(complex_factoryT(world)));
336 u.push_back(complex_factoryT(world).functor(complex_functorT(new Guess( zeta, 0.03, 0, 0, 0))).nofence()); // s
337 v.push_back(complex_functionT(complex_factoryT(world)));
338 u.push_back(complex_factoryT(world).functor(complex_functorT(new Guess(-zeta, 0.03, 0, 0, 0))).nofence()); // s
339 v.push_back(complex_functionT(complex_factoryT(world)));
340 //
341 v.push_back(complex_factoryT(world).functor(complex_functorT(new Guess( zeta, 0.15, 0, 0, 0))).nofence()); // s
342 u.push_back(complex_functionT(complex_factoryT(world)));
343 v.push_back(complex_factoryT(world).functor(complex_functorT(new Guess(-zeta, 0.15, 0, 0, 0))).nofence()); // s
344 u.push_back(complex_functionT(complex_factoryT(world)));
345 v.push_back(complex_factoryT(world).functor(complex_functorT(new Guess( zeta, 0.15, 0, 0, 1))).nofence()); // pz
346 u.push_back(complex_functionT(complex_factoryT(world)));
347 v.push_back(complex_factoryT(world).functor(complex_functorT(new Guess(-zeta, 0.15, 0, 0, 1))).nofence()); // pz
348 u.push_back(complex_functionT(complex_factoryT(world)));
349 v.push_back(complex_factoryT(world).functor(complex_functorT(new Guess( zeta, 0.15, 0, 1, 0))).nofence()); // py
350 u.push_back(complex_functionT(complex_factoryT(world)));
351 v.push_back(complex_factoryT(world).functor(complex_functorT(new Guess(-zeta, 0.15, 0, 1, 0))).nofence()); // py
352 u.push_back(complex_functionT(complex_factoryT(world)));
353 v.push_back(complex_factoryT(world).functor(complex_functorT(new Guess( zeta, 0.15, 1, 0, 0))).nofence()); // px
354 u.push_back(complex_functionT(complex_factoryT(world)));
355 v.push_back(complex_factoryT(world).functor(complex_functorT(new Guess(-zeta, 0.15, 1, 0, 0))).nofence()); // px
356 u.push_back(complex_functionT(complex_factoryT(world)));
357 v.push_back(complex_factoryT(world).functor(complex_functorT(new Guess( zeta, 0.03, 0, 0, 0))).nofence()); // s
358 u.push_back(complex_functionT(complex_factoryT(world)));
359 v.push_back(complex_factoryT(world).functor(complex_functorT(new Guess(-zeta, 0.03, 0, 0, 0))).nofence()); // s
360 u.push_back(complex_functionT(complex_factoryT(world)));
361 world.gop.fence();
362 normalize2(world, u, v);
363
364
365 int nvec=u.size();
366 Tensor<double_complex> H1(nvec,nvec);
367 Tensor<double_complex> S1(nvec,nvec);
368
369 Tensor<double_complex> c;
370 Tensor<double> e, e1(nvec);
371 double maxerr;
372 real_functionT V0, V0x, V0y, V0z;
373 //double shift=0.;
374 bool doso = false; // turned on once converged to 1e-4
375 if (world.rank() == 0) {
376 print(" u size ", u.size(), "v size ", v.size()); }
377
378 doitagain:
379
380 while (thresh > 0.9e-10) {
381 if (world.rank() == 0) {
382 print("\n Solving with thresh",thresh,"k",k,"at time", wall_time(), "\n");
383 printf(" iter root energy err time\n");
384 printf(" ---- ---- ------------------ ------- ------\n");
385 }
386
387 V0 = real_factoryT(world).f(V).thresh(thresh*.1);
388 V0x = real_derivative_3d(world,0)(V0);
389 V0y = real_derivative_3d(world,1)(V0);
390 V0z = real_derivative_3d(world,2)(V0);
391
392 // If it does not converge in 10 iters it probably needs more precision.
393 for (int iter=0; iter<10; iter++) {
394
395 // loadbalance(world, u, v, V0, V0x, V0y, V0z);
396
397 std::vector<complex_functionT> Vu, Vv, du[3], dv[3];
398
399 apply_potential(world, V0, V0x, V0y, V0z, u, v, Vu, Vv, du, dv, doso);
400 world.gop.fence(); // free memory
401
402 Tensor<double_complex> H = hamiltonian_matrix(world, u, v, Vu, Vv, du, dv);
403 Tensor<double_complex> S = matrix_inner(world, u, u) + matrix_inner(world, v, v);
404
405 if ( iter==0) {
406 for (int iii = 0; iii < nvec; iii++ )
407 for (int jjj = 0; jjj < nvec; jjj++ ){
408 H1(jjj,iii) = double_complex(0.5,0.0)*(H(iii,jjj)+H(jjj,iii));
409 }
410 }
411 else
412 for (int iii = 0; iii < nvec; iii++ )
413 for (int jjj = 0; jjj < nvec; jjj++ ){
414 H1(jjj,iii) = H(jjj,iii);
415 }
416
417 if ( iter==0) {
418 for (int iii = 0; iii < nvec; iii++ )
419 for (int jjj = 0; jjj < nvec; jjj++ ){
420 S1(jjj,iii) = double_complex(0.5,0.0)*(S(iii,jjj)+S(jjj,iii));
421 }
422 }
423 else
424 for (int iii = 0; iii < nvec; iii++ )
425 for (int jjj = 0; jjj < nvec; jjj++ ){
426 S1(jjj,iii) = S(jjj,iii);
427 }
428
429 world.gop.fence(); // free memory
430
431 sygv(H1, S1, 1, c, e);
432
433 for (int axis=0; axis<3; axis++) {
434 du[axis].clear();
435 dv[axis].clear();
436 }
437
438 world.gop.fence(); // free memory
439 e = e(Slice(0,nvec-1));
440 u = transform(world, u, c(_,Slice(0,nvec-1)));
441 v = transform(world, v, c(_,Slice(0,nvec-1)));
442
443 Vu = transform(world, Vu, c(_,Slice(0,nvec-1)));
444 Vv = transform(world, Vv, c(_,Slice(0,nvec-1)));
445 world.gop.fence(); // free memory
446
447 truncate(world, u);
448 truncate(world, v);
449 truncate(world, Vu);
450 truncate(world, Vv);
451 world.gop.fence();
452
453 normalize2(world, u, v);
454
455 world.gop.fence(); // free memory
456
457 for (int iii = 0; iii < nvec; iii++ )
458 e1[iii] = e[iii]; // +shift;
459
460
461 for (int i=0; i<nvec; i++) {
462 Vu[i].refine();
463 Vv[i].refine();
464 }
465
466 world.gop.fence();
467 std::vector<poperatorT> ops = make_bsh_operators(world, e1, thresh);
468
469 std::vector<complex_functionT> u_new = apply(world, ops, Vu);
470 std::vector<complex_functionT> v_new = apply(world, ops, Vv);
471
472 normalize2(world, u_new, v_new);
473
474 Vu.clear();
475 Vv.clear();
476 world.gop.fence();
477
478 std::vector<double> rnormu = norm2s(world,add(world, u, u_new));
479 std::vector<double> rnormv = norm2s(world,add(world, v, v_new));
480 std::vector<double> rnorm(nvec);
481
482 for (int i=0; i<nvec; i++)
483 rnorm[i] = sqrt(rnormu[i]*rnormu[i] + rnormv[i]*rnormv[i]);
484 world.gop.fence();
485
486 maxerr= 0.;
487 for (int i=0; i<nvec; i++) {
488 if (world.rank() == 0) printf(" %3d %3d %18.12f %.1e %7.1f\n",
489 iter, i, e[i]/reduced, rnorm[i], wall_time());
490 maxerr = std::max(maxerr,rnorm[i]);
491 }
492
493 u = u_new;
494 v = v_new;
495
496 u_new.clear();
497 v_new.clear();
498
499 if (maxerr < 1.e-3 && !doso) {
500 // We initially converged to a threshold of 1e-4 without SO
501 // and now we repeat the 1e-4 iteration with SO and continue
502 // with SO on.
503 if (world.rank() == 0) print("\n Turning on the SO interaction \n");
504 doso = true;
505 goto doitagain; // Gotta love it
506 break;
507 }
508 }
509
510 moments(world, u, v);
511
512 thresh *= 1e-1;
513 k += 1;
514 FunctionDefaults<3>::set_k(k);
515 FunctionDefaults<3>::set_thresh(thresh);
516 reconstruct(world,u);
517 reconstruct(world,v);
518 for (unsigned int i=0; i<u.size(); i++) {
519 u[i] = madness::project(u[i], k, thresh);
520 v[i] = madness::project(v[i], k, thresh);
521 }
522 V0 = real_factoryT(world).f(V).thresh(thresh*.1);
523 V0x = real_derivative_3d(world,0)(V0);
524 V0y = real_derivative_3d(world,1)(V0);
525 V0z = real_derivative_3d(world,2)(V0);
526 world.gop.fence();
527
528 }
529 }
530
main(int argc,char ** argv)531 int main(int argc, char** argv) {
532 initialize(argc, argv);
533 World world(SafeMPI::COMM_WORLD);
534 startup(world,argc,argv);
535
536 // cpu_set_t mask;
537 // CPU_ZERO(&mask);
538 // CPU_SET(MPI::COMM_WORLD.Get_rank(), &mask);
539 // if( sched_setaffinity( 0, sizeof(mask), &mask ) == -1 ) {
540 // printf("WARNING: Could not set CPU Affinity, continuing...\n");
541 // }
542
543 try {
544 doit(world);
545 } catch (const SafeMPI::Exception& e) {
546 print(e);
547 error("caught an MPI exception");
548 } catch (const madness::MadnessException& e) {
549 print(e);
550 error("caught a MADNESS exception");
551 } catch (const madness::TensorException& e) {
552 print(e);
553 error("caught a Tensor exception");
554 } catch (const char* s) {
555 print(s);
556 error("caught a string exception");
557 } catch (const std::string& s) {
558 print(s);
559 error("caught a string (class) exception");
560 } catch (const std::exception& e) {
561 print(e.what());
562 error("caught an STL exception");
563 } catch (...) {
564 error("caught unhandled exception");
565 }
566
567 finalize();
568 return 0;
569 }
570