MADNESS 0.10.1
DKops.h
Go to the documentation of this file.
1#ifndef MADNESS_APPS_MOLDFT_DKOPS_H_INCLUDED
2#define MADNESS_APPS_MOLDFT_DKOPS_H_INCLUDED
3
4
5#include <madness/mra/mra.h>
6#include <madness/constants.h>
7#include <madness/mra/operator.h>
8#include <iostream>
9#include <fstream>
10#include <exception>
11
12//#include "../../../../mpfrc++-3/mpreal.h"
13
14using namespace madness;
15
16static const double opthresh = 1e-16;
17
18double q(double t){
19 return exp(-t);
20}
21
22double w(double t, double eps){
23 double c = 137.0359895;
24 double PI = constants::pi;
25 return 1.0/(c*sqrt(PI))*exp(-1.0/2.0*t-c*c*exp(-t))-(1.0+eps/(c*c))*exp((2.0*eps+(eps*eps)/(c*c))*exp(-t)-t)*(erfc((c+eps/c)*(exp(-t/2.0)))-2.0);
26}
27
28double Ebark(double t, double eps){
29 double R = 1e-8;
30 return pow(1.0/(2.0*q(t)), 1.5)*w(t, eps)*exp(-1.0/(4.0*q(t))*R*R);
31}
32
33real_convolution_3d Ebar(World& world, double eps){
34 Vector<double,3> args = vec(0.0,0.0,0.0);
35
36
37 std::vector<double> cvec;
38 std::vector<double> tvec;
39
40 double fac = pow(2.0*constants::pi,-1.5);
41 double dt = 1.0/8.0;
42 double tval = -10.0;
43 while(tval < 100.0){
44
45 if(Ebark(tval, eps) > opthresh){
46 //if(world.rank()==0) print(" Ebark(",tval,",",eps," = ",Ebark(tval,eps));
47 cvec.push_back(dt/pow(2.0*q(tval),1.5)*w(tval,eps)*fac);
48 tvec.push_back(1.0/(4.0*q(tval)));
49 }
50 //else{
51 //if(world.rank()==0) print("Ebark(",tval,",",eps," = ",Ebark(tval,eps));
52
53 //}
54 tval = tval + dt;
55 }
56
57 Tensor<double> ctensor(cvec.size());
58 Tensor<double> ttensor(tvec.size());
59 for(int i = 0; i < cvec.size(); i++){
60 ctensor[i] = cvec[i];
61 ttensor[i] = tvec[i];
62 }
63
64 int n = cvec.size();
65 //if(world.rank()==0) print("Made an Ebar! n = ",n);
66 //if(world.rank()==0) print("Made an Ebar! Here's what's inside:\n\nc:\n",ctensor,"\nt:\n",ttensor);
67
68 //test ebar?
69 /*
70 if(world.rank()==0) print("Testing Ebar. Value should be ~22.680321967196");
71 double testvalue = 0.0;
72 for(int i = 0; i < cvec.size(); i++){
73 testvalue += ctensor[i]*exp(-1e-2*ttensor[i]);
74 }
75 if(world.rank()==0) print("And the result is: ", testvalue);
76 */
77
78 //return real_convolution_3d(world, args, ctensor, ttensor);
79 return real_convolution_3d(world, ctensor, ttensor);
80}
81
82//real_convolution_3d Ebar_fixed(World& world){
83// //Tensor<double> c(305l), t(305l);
84// std::vector<double> c(365), t(365);
85// Vector<double,3> args = vec(0.0,0.0,0.0);
86// #include "RelCoeffs/Ebar_coeffs.dat"
87// int n = 305;
88//
89// for(int i=0; i < n; i++){
90// if(c[i]*exp(-t[i]*1e-16) < opthresh){
91// c.erase(c.begin()+i);
92// t.erase(t.begin()+i);
93// i--;
94// n--;
95// }
96// }
97//
98// //print("n = ", n);
99//
100// Tensor<double> ctens(n), ttens(n);
101// for(int i = 0; i < n; i++){
102// ctens[i] = c[i];
103// ttens[i] = t[i];
104// }
105// //if(world.rank()==0) print("Made a Ebar (fixed)! n = ", n);
106// //if(world.rank()==0) print("Made a Pbar! Here's what's inside:\n\nc:\n",ctens,"\nt:\n",ttens);
107// //return real_convolution_3d(world, args, ctens, ttens);
108// return real_convolution_3d(world, ctens, ttens);
109//}
110
111real_convolution_3d Pbar(World& world){
112 //Tensor<double> c(305l), t(305l);
113 std::vector<double> c;
114 std::vector<double> t;
115 Vector<double,3> args = vec(0.0,0.0,0.0);
116
117 std::ifstream inf("/gpfs/home/jscanderson/DKproject/Pbar_t.csv");
118 if(!inf){
119 if(world.rank() == 0) std::cerr << "Unable to open Pbar_t.csv" << std::endl;
120 exit(1);
121 }
122 std::string strInput;
123 getline(inf, strInput);
124 while(inf){
125 t.push_back(std::stod(strInput));
126 getline(inf, strInput);
127 }
128 inf.close();
129 //if(world.rank()==0) print("Done Reading Pbar_t");
130
131 inf.open("/gpfs/home/jscanderson/DKproject/Pbar_c.csv");
132 if(!inf){
133 if(world.rank() == 0) std::cerr << "Unable to open Pbar_c.csv" << std::endl;
134 exit(1);
135 }
136 getline(inf, strInput);
137 while(inf){
138 //if(world.rank()==0) print("a", strInput, "b");
139 c.push_back(std::stod(strInput));
140 getline(inf, strInput);
141 }
142 inf.close();
143 //if(world.rank()==0) print("Done Reading Pbar_c");
144
145
146 int n = c.size();
147 /*
148 for(int i=0; i < n; i++){
149 if(c[i]*exp(-t[i]*1e-30) < opthresh){
150 c.erase(c.begin()+i);
151 t.erase(t.begin()+i);
152 i--;
153 n--;
154 }
155 }
156 */
157 //print("n = ", n);
158
159 Tensor<double> ctens(n), ttens(n);
160 for(int i = 0; i < n; i++){
161 ctens[i] = c[i];
162 ttens[i] = t[i];
163 }
164 //if(world.rank()==0) print("Made a Pbar! n = ", n);
165 //if(world.rank()==0) print("Made a Pbar! Here's what's inside:\n\nc:\n",ctens,"\nt:\n",ttens);
166 //return real_convolution_3d(world, args, ctens, ttens);
167 return real_convolution_3d(world, ctens, ttens);
168}
169
170real_convolution_3d Tbar(World& world){
171 //Tensor<double> c(305l), t(305l);
172 std::vector<double> c;
173 std::vector<double> t;
174 Vector<double,3> args = vec(0.0,0.0,0.0);
175
176 std::ifstream inf("/gpfs/home/jscanderson/DKproject/Tbar_t.csv");
177 if(!inf){
178 if(world.rank() == 0) std::cerr << "Unable to open Tbar_t.csv" << std::endl;
179 exit(1);
180 }
181 std::string strInput;
182 getline(inf, strInput);
183 while(inf){
184 t.push_back(std::stod(strInput));
185 getline(inf, strInput);
186 }
187 inf.close();
188 //if(world.rank()==0) print("Done Reading Tbar_t");
189
190 inf.open("/gpfs/home/jscanderson/DKproject/Tbar_c.csv");
191 if(!inf){
192 if(world.rank() == 0) std::cerr << "Unable to open Tbar_c.csv" << std::endl;
193 exit(1);
194 }
195 getline(inf, strInput);
196 while(inf){
197 //if(world.rank()==0) print("a", strInput, "b");
198 c.push_back(std::stod(strInput));
199 getline(inf, strInput);
200 }
201 inf.close();
202 //if(world.rank()==0) print("Done Reading Tbar_c");
203
204
205 int n = c.size();
206 /*
207 for(int i=0; i < n; i++){
208 if(c[i]*exp(-t[i]*1e-30) < opthresh){
209 c.erase(c.begin()+i);
210 t.erase(t.begin()+i);
211 i--;
212 n--;
213 }
214 }
215 */
216 //print("n = ", n);
217
218 Tensor<double> ctens(n), ttens(n);
219 for(int i = 0; i < n; i++){
220 ctens[i] = c[i];
221 ttens[i] = t[i];
222 }
223 //if(world.rank()==0) print("Made a Tbar! n = ", n);
224 //if(world.rank()==0) print("Made a Pbar! Here's what's inside:\n\nc:\n",ctens,"\nt:\n",ttens);
225 //return real_convolution_3d(world, args, ctens, ttens);
226 return real_convolution_3d(world, ctens, ttens);
227}
228
229
230real_convolution_3d A(World& world){
231 //Tensor<double> c(301l), t(301l);
232 std::vector<double> c, t;
233 Vector<double,3> args = vec(0.0,0.0,0.0);
234
235 std::ifstream inf("/gpfs/home/jscanderson/DKproject/A_t.csv");
236 if(!inf){
237 if(world.rank() == 0) std::cerr << "Unable to open A_t.csv" << std::endl;
238 exit(1);
239 }
240 std::string strInput;
241 getline(inf, strInput);
242 while(inf){
243 t.push_back(std::stod(strInput));
244 getline(inf, strInput);
245 }
246 inf.close();
247 //if(world.rank()==0) print("Done Reading A_t");
248
249 inf.open("/gpfs/home/jscanderson/DKproject/A_c.csv");
250 if(!inf){
251 if(world.rank() == 0) std::cerr << "Unable to open A_c.csv" << std::endl;
252 exit(1);
253 }
254 getline(inf, strInput);
255 while(inf){
256 c.push_back(std::stod(strInput));
257 getline(inf, strInput);
258 }
259 inf.close();
260 //if(world.rank()==0) print("Done Reading A_c");
261
262 int n = c.size();
263 /*
264 for(int i=0; i < n; i++){
265 if(c[i]*exp(-t[i]*1e-30) < opthresh){
266 c.erase(c.begin()+i);
267 t.erase(t.begin()+i);
268 i--;
269 n--;
270 }
271 }
272 */
273 double fac = pow(2.0*constants::pi,1.5);
274 Tensor<double> ctens(n), ttens(n);
275 for(int i = 0; i < n; i++){
276 ctens[i] = c[i]*fac;
277 ttens[i] = t[i];
278 }
279 //if(world.rank()==0) print("Made an A!, n = ", n);
280 //if(world.rank()==0) print("Made an A! Here's what's inside:\n\nc:\n",ctens,"\nt:\n",ttens);
281 //return real_convolution_3d(world, args, ctens, ttens);
282 return real_convolution_3d(world, ctens, ttens);
283}
284
285//Same as A(World& world) above but adds 1/sqrt(2) to the largest coefficient
286real_convolution_3d A2(World& world){
287 //Tensor<double> c(301l), t(301l);
288 std::vector<double> c, t;
289 Vector<double,3> args = vec(0.0,0.0,0.0);
290
291 std::ifstream inf("/gpfs/home/jscanderson/DKproject/A_t.csv");
292 if(!inf){
293 if(world.rank() == 0) std::cerr << "Unable to open A_t.csv" << std::endl;
294 exit(1);
295 }
296 std::string strInput;
297 getline(inf, strInput);
298 while(inf){
299 t.push_back(std::stod(strInput));
300 getline(inf, strInput);
301 }
302 inf.close();
303 //if(world.rank()==0) print("Done Reading A_t");
304
305 inf.open("/gpfs/home/jscanderson/DKproject/A_c.csv");
306 if(!inf){
307 if(world.rank() == 0) std::cerr << "Unable to open A_c.csv" << std::endl;
308 exit(1);
309 }
310 getline(inf, strInput);
311 while(inf){
312 c.push_back(std::stod(strInput));
313 getline(inf, strInput);
314 }
315 inf.close();
316 //if(world.rank()==0) print("Done Reading A_c");
317
318 int n = c.size();
319 /*
320 for(int i=0; i < n; i++){
321 if(c[i]*exp(-t[i]*1e-30) < opthresh){
322 c.erase(c.begin()+i);
323 t.erase(t.begin()+i);
324 i--;
325 n--;
326 }
327 }
328 */
329 double fac = pow(2.0*constants::pi,1.5);
330 Tensor<double> ctens(n), ttens(n);
331 for(int i = 0; i < n; i++){
332 ctens[i] = c[i]*fac;
333 ttens[i] = t[i];
334 }
335 c[n-1] += 1/sqrt(2)*fac*std::pow(t[n-1]*constants::pi,1.5); //largest exponent is at the end, so add 1/sqrt(2) to the coefficient there
336 //if(world.rank()==0) print("Made an A!, n = ", n);
337 //if(world.rank()==0) print("Made an A! Here's what's inside:\n\nc:\n",ctens,"\nt:\n",ttens);
338 //return real_convolution_3d(world, args, ctens, ttens);
339 return real_convolution_3d(world, ctens, ttens);
340}
341
342real_convolution_3d PbarA(World& world){
343 //Tensor<double> c(441l), t(441l);
344 std::vector<double> c, t;
345 Vector<double,3> args = vec(0.0,0.0,0.0);
346
347 std::ifstream inf("/gpfs/home/jscanderson/DKproject/PbarA_t.csv");
348 if(!inf){
349 if(world.rank() == 0) std::cerr << "Unable to open PbarA_t.csv" << std::endl;
350 exit(1);
351 }
352 std::string strInput;
353 getline(inf, strInput);
354 while(inf){
355 t.push_back(std::stod(strInput));
356 getline(inf, strInput);
357 }
358 inf.close();
359 //if(world.rank()==0) print("Done Reading PbarA_t");
360
361 inf.open("/gpfs/home/jscanderson/DKproject/PbarA_c.csv");
362 if(!inf){
363 if(world.rank() == 0) std::cerr << "Unable to open PbarA_c.csv" << std::endl;
364 exit(1);
365 }
366 getline(inf, strInput);
367 while(inf){
368 c.push_back(std::stod(strInput));
369 getline(inf, strInput);
370 }
371 inf.close();
372 //if(world.rank()==0) print("Done Reading PbarA_c");
373
374 int n = c.size();
375 /*
376 for(int i=0; i < n; i++){
377 if(c[i]*exp(-t[i]*1e-30) < opthresh){
378 c.erase(c.begin()+i);
379 t.erase(t.begin()+i);
380 i--;
381 n--;
382 }
383 }
384 */
385 double fac = pow(2.0*constants::pi,1.5);
386 Tensor<double> ctens(n), ttens(n);
387 for(int i = 0; i < n; i++){
388 ctens[i] = c[i]*fac;
389 ttens[i] = t[i];
390 }
391 //if(world.rank()==0) print("Made a PbarA!, n = ", n);
392 //if(world.rank()==0) print("Made a PbarA! Here's what's inside:\n\nc:\n",ctens,"\nt:\n",ttens);
393 //return real_convolution_3d(world, args, ctens, ttens);
394 return real_convolution_3d(world, ctens, ttens);
395}
396
397std::vector<real_convolution_3d_ptr> gradPbarA(World& world){
398 const double pi = constants::pi;
399 std::vector<double> c, t;
400
401 std::ifstream inf("/gpfs/home/jscanderson/DKproject/PbarA_t.csv");
402 if(!inf){
403 if(world.rank() == 0) std::cerr << "Unable to open PbarA_t.csv" << std::endl;
404 exit(1);
405 }
406 std::string strInput;
407 getline(inf, strInput);
408 while(inf){
409 t.push_back(std::stod(strInput));
410 getline(inf, strInput);
411 }
412 inf.close();
413 //if(world.rank()==0) print("Done Reading PbarA_t");
414
415 inf.open("/gpfs/home/jscanderson/DKproject/PbarA_c.csv");
416 if(!inf){
417 if(world.rank() == 0) std::cerr << "Unable to open PbarA_c.csv" << std::endl;
418 exit(1);
419 }
420 getline(inf, strInput);
421 while(inf){
422 c.push_back(std::stod(strInput));
423 getline(inf, strInput);
424 }
425 inf.close();
426 //if(world.rank()==0) print("Done Reading PbarA_c");
427
428 int n = c.size();
429 /*
430 for(int i=0; i < n; i++){
431 if(c[i]*exp(-t[i]*1e-30) < opthresh){
432 c.erase(c.begin()+i);
433 t.erase(t.begin()+i);
434 i--;
435 n--;
436 }
437 }
438 */
439 double fac = pow(2.0*pi,1.5);
440 for(int i = 0; i < n; i++){
441 c[i] = c[i]*fac;
442 }
443
444 const Tensor<double> width = FunctionDefaults<3>::get_cell_width();
445 const int k = FunctionDefaults<3>::get_k();
446 double hi = width.normf(); // Diagonal width of cell
447
448 std::vector<real_convolution_3d_ptr> gradG(3);
449
450 for (int dir=0; dir<3; dir++) {
451 std::vector< ConvolutionND<double,3> > ops(n);
452 for (int mu=0; mu<n; mu++) {
453 // We cache the normalized operator so the factor is the value we must multiply
454 // by to recover the coeff we want.
455 double cc = std::pow(t[mu]/pi,1.5); // Normalization coeff
456 ops[mu].setfac(c[mu]/cc/width[dir]);
457
458 for (int d=0; d<3; d++) {
459 if (d != dir) ops[mu].setop(d,GaussianConvolution1DCache<double>::get(k, t[mu]*width[d]*width[d], 0, false));
460 }
461 ops[mu].setop(dir,GaussianConvolution1DCache<double>::get(k, t[mu]*width[dir]*width[dir], 1, false));
462 }
463 gradG[dir] = real_convolution_3d_ptr(new SeparatedConvolution<double,3>(world, ops));
464 }
465
466 return gradG;
467
468}
469
470#endif
Tbar
real_convolution_3d Tbar(World &world)
Definition DKops.h:170
A2
real_convolution_3d A2(World &world)
Definition DKops.h:286
Ebar
real_convolution_3d Ebar(World &world, double eps)
Definition DKops.h:33
w
double w(double t, double eps)
Definition DKops.h:22
Pbar
real_convolution_3d Pbar(World &world)
Definition DKops.h:111
PbarA
real_convolution_3d PbarA(World &world)
Definition DKops.h:342
opthresh
static const double opthresh
Definition DKops.h:16
q
double q(double t)
Definition DKops.h:18
gradPbarA
std::vector< real_convolution_3d_ptr > gradPbarA(World &world)
Definition DKops.h:397
Ebark
double Ebark(double t, double eps)
Definition DKops.h:28
A
Definition test_ar.cc:118
madness::BaseTensor::size
long size() const
Returns the number of elements in the tensor.
Definition basetensor.h:138
madness::FunctionDefaults::get_k
static int get_k()
Returns the default wavelet order.
Definition funcdefaults.h:266
madness::FunctionDefaults::get_cell_width
static const Tensor< double > & get_cell_width()
Returns the width of each user cell dimension.
Definition funcdefaults.h:468
madness::Tensor
A tensor is a multidimension array.
Definition tensor.h:317
madness::Tensor::normf
float_scalar_type normf() const
Returns the Frobenius norm of the tensor.
Definition tensor.h:1726
madness::Vector
A simple, fixed dimension vector.
Definition vector.h:64
madness::World
A parallel world class.
Definition world.h:132
madness::World::rank
ProcessID rank() const
Returns the process rank in this World (same as MPI_Comm_rank()).
Definition world.h:318
constants.h
Defines common mathematical and physical constants.
R
static const double R
Definition csqrt.cc:46
pow
static double pow(const double *a, const double *b)
Definition lda.h:74
mra.h
Main include file for MADNESS and defines Function interface.
madness::constants::pi
const double pi
Mathematical constant .
Definition constants.h:48
madness
Namespace for all elements and tools of MADNESS.
Definition DFParameters.h:10
madness::real_convolution_3d_ptr
std::shared_ptr< real_convolution_3d > real_convolution_3d_ptr
Definition functypedefs.h:135
madness::real_convolution_3d
SeparatedConvolution< double, 3 > real_convolution_3d
Definition functypedefs.h:121
madness::dir
static const string dir
Definition corepotential.cc:249
madness::vec
Vector< T, sizeof...(Ts)+1 > vec(T t, Ts... ts)
Factory function for creating a madness::Vector.
Definition vector.h:711
cc
const double cc
Definition navstokes_cosines.cc:107
mu
const double mu
Definition navstokes_cosines.cc:95
d
static const double d
Definition nonlinschro.cc:121
operator.h
Implements most functionality of separated operators.
c
static const double c
Definition relops.cc:10
PI
static const double PI
Definition relops.cc:12
k
static const long k
Definition rk.cc:44
madness::GaussianConvolution1DCache
Definition convolution1d.h:849
e
void e()
Definition test_sig.cc:75
pi
static const double pi
Definition testcosine.cc:6
Generated on Thu Dec 19 2024 21:22:52 for MADNESS by doxygen 1.9.8