MADNESS 0.10.1
lowrankfunction.h
Go to the documentation of this file.
1//
2// Created by Florian Bischoff on 8/10/23.
3//
4
5#ifndef MADNESS_LOWRANKFUNCTION_H
6#define MADNESS_LOWRANKFUNCTION_H
7
8
9#include<madness/mra/mra.h>
10#include<madness/mra/vmra.h>
14#include <random>
15
16
17
18namespace madness {
19
21
23
24 // initialize with: key, value, comment (optional), allowed values (optional)
25 initialize<double>("radius",2.0,"the radius");
26 initialize<double>("gamma",1.0,"the exponent of the correlation factor");
27 initialize<double>("volume_element",0.1,"volume covered by each grid point");
28 initialize<double>("tol",1.e-8,"rank-reduced cholesky tolerance");
29 initialize<std::string>("f12type","Slater","correlation factor",{"Slater","SlaterF12"});
30 initialize<std::string>("orthomethod","cholesky","orthonormalization",{"cholesky","canonical","symmetric"});
31 initialize<std::string>("transpose","slater2","transpose of the matrix",{"slater1","slater2"});
32 initialize<std::string>("gridtype","random","the grid type",{"random","cartesian","dftgrid"});
33 initialize<std::string>("rhsfunctiontype","exponential","the type of function",{"exponential"});
34 initialize<int>("optimize",1,"number of optimization iterations");
35 }
36
37 void read_and_set_derived_values(World& world, const commandlineparser& parser, std::string tag) {
38 read_input_and_commandline_options(world,parser,tag);
39 }
40
41 double radius() const {return get<double>("radius");}
42 double gamma() const {return get<double>("gamma");}
43 double volume_element() const {return get<double>("volume_element");}
44 double tol() const {return get<double>("tol");}
45 int optimize() const {return get<int>("optimize");}
46 std::string gridtype() const {return get<std::string>("gridtype");}
47 std::string orthomethod() const {return get<std::string>("orthomethod");}
48 std::string rhsfunctiontype() const {return get<std::string>("rhsfunctiontype");}
49 std::string f12type() const {return get<std::string>("f12type");}
50 };
51
52
53 class gridbase {
54 public:
55 double get_volume_element() const {return volume_element;}
56 double get_radius() const {return radius;}
57
58 virtual ~gridbase() = default;
59 // visualize the grid in xyz format
60 template<std::size_t NDIM>
61 void visualize(const std::string filename, const std::vector<Vector<double,NDIM>>& grid) const {
62 print("visualizing grid to file",filename);
63 print("a total of",grid.size(),"grid points");
64 std::ofstream file(filename);
65 for (const auto& r : grid) {
66 // formatted output
67 file << std::fixed << std::setprecision(6);
68 for (int i=0; i<NDIM; ++i) file << r[i] << " ";
69 file << std::endl;
70 }
71 file.close();
72 }
73
74 protected:
75 double volume_element=0.1;
76 double radius=3;
77 bool do_print=false;
78 };
79
80 template<std::size_t NDIM>
81 class dftgrid : public gridbase {
82 public:
84 explicit dftgrid(const double volume_element, const double radius) {
85 // increase number of radial grid points until the volume element is below the threshold
86 double current_ve=1.0;
87 std::size_t nradial=10;
89 nradial+=10;
90 print("trying nradial",nradial);
93 tmp.set_nradial(nradial);
94 tmp.make_grid();
95 double volume=4./3. *M_PI * std::pow(radius,3.0);
96 auto npoints=tmp.get_points().size();
97 current_ve=volume/npoints;
98 print("volume, npoints, volume element",volume,npoints,current_ve);
99 }
100 builder.set_nradial(nradial);
102 }
103
104
105 explicit dftgrid(const std::size_t nradial, const std::size_t angular_order, const coord_3d origin=coord_3d()) {
106 static_assert(NDIM==3,"DFT Grids only defined for NDIM=3");
107 builder.set_nradial(nradial);
108 builder.set_angular_order(angular_order);
109 builder.set_origin(origin);
111 }
112
113 std::vector<Vector<double,NDIM>> get_grid() const {
114 return builder.get_points();
115 }
116
117 };
118
119 /// grid with random points around the origin, with a Gaussian distribution
120 template<std::size_t NDIM>
121 class randomgrid : public gridbase {
122 public:
124 : gridbase(), origin(origin) {
125 this->volume_element=volume_element;
126 this->radius=radius;
127 }
128
129 std::vector<Vector<double,NDIM>> get_grid() const {
130 std::vector<Vector<double, NDIM>> grid;
132 if (do_print) print("npoint_within_volume",npoint_within_volume);
133
135 auto is_in_cell = [&cell](const Vector<double, NDIM>& r) {
136 for (int d = 0; d < NDIM; ++d) if (r[d] < cell(d, 0) or r[d] > cell(d, 1)) return false;
137 return true;
138 };
139 double rad=radius;
140 auto o=origin;
141 auto is_in_sphere = [&rad,&o](const Vector<double, NDIM>& r) {
142 return ((r-o).normf()<rad);
143 };
144
145 // set variance such that about 70% of all grid points sits within the radius
146 double variance=radius;
147 if (NDIM==2) variance=0.6*radius;
148 if (NDIM==3) variance=0.5*radius;
149 long maxrank=10*npoint_within_volume;
150 long rank=0;
151 for (int r=0; r<maxrank; ++r) {
153 if (not is_in_cell(tmp)) continue;
154 if (is_in_sphere(tmp)) ++rank;
155 grid.push_back(tmp);
156 if (rank==npoint_within_volume) break;
157 }
158 if (do_print) {
159 print("origin ",origin);
160 print("radius ",radius);
161 print("grid points in volume ",rank);
162 print("total grid points ",grid.size());
163 print("ratio ",rank/double(grid.size()));
164 print("volume element ",volume()/rank);
165 }
166 return grid;
167 }
168
170 return origin;
171 }
172
173 private:
174
175 double volume() const {
177 if (NDIM==1) return 2.0*radius;
178 if (NDIM==2) return constants::pi*radius*radius;
179 if (NDIM==3) return 4.0 / 3.0 * constants::pi * std::pow(radius, 3.0);
180 }
181
183 std::random_device rd{};
184 std::mt19937 gen{rd()};
185 Vector<double,NDIM> result;
186 for (int i = 0; i < NDIM; ++i) {
187 std::normal_distribution<> d{origin[i], variance};
188 result[i]=d(gen);
189 }
190
191 return result;
192 }
193
195
196 };
197
198 template<std::size_t NDIM>
201 std::vector<long> stride;
202 long index=0;
206
207 cartesian_grid(const double volume_per_gridpoint, const double radius) {
210 print("length per gridpoint, n_per_dim",length_per_gridpoint,n_per_dim);
211 print("volume_per_gridpoint",std::pow(length_per_gridpoint,NDIM));
212 initialize(-radius,radius);
213 print("increment",increment);
214 }
215
216
217 cartesian_grid(const long n_per_dim, const double lo, const double hi)
219 initialize(lo,hi);
220 }
221
223 hivec(other.hivec), stride(other.stride), index(0), n_per_dim(other.n_per_dim),
224 total_n(other.total_n), increment(other.increment) {
225 }
226
229 std::swap(*this,other);
230 return *this;
231 }
232
233 void initialize(const double lo, const double hi) {
234 lovec.fill(lo);
235 hivec.fill(hi);
236 increment=(hivec-lovec)*(1.0/double(n_per_dim-1));
237 stride=std::vector<long>(NDIM,1l);
238 total_n=std::pow(n_per_dim,NDIM);
239 for (long i=NDIM-2; i>=0; --i) stride[i]=n_per_dim*stride[i+1];
240
241 }
242
243 double volume_per_gridpoint() const{
244 double volume=1.0;
245 for (int i=0; i<NDIM; ++i) volume*=(hivec[i]-lovec[i]);
246 return volume/total_n;
247 }
248
249 void operator++() {
250 index++;
251 }
252
253 bool operator()() const {
254 return index < total_n;
255 }
256
259 for (int idim=0; idim<NDIM; ++idim) {
260 tmp[idim]=(index/stride[idim])%n_per_dim;
261 }
262 return lovec+tmp*increment;
263 }
264
265 };
266
267 /// given a molecule, return a suitable grid
268 template<std::size_t NDIM>
269 class molecular_grid : public gridbase {
270
271 public:
272 /// ctor takes centers of the grid and the grid parameters
274 : centers(origins)
275 {
276 if (centers.size()==0) centers.push_back(Vector<double,NDIM>(0) );
277 if (params.gridtype()=="random") grid_builder=std::make_shared<randomgrid<NDIM>>(params.volume_element(),params.radius());
278 // else if (params.gridtype()=="cartesian") grid_builder=std::make_shared<cartesian_grid<NDIM>>(params.volume_element(),params.radius());
279 else if (params.gridtype()=="dftgrid") {
280 if constexpr (NDIM==3) {
281 grid_builder=std::make_shared<dftgrid<NDIM>>(params.volume_element(),params.radius());
282 } else {
283 MADNESS_EXCEPTION("no dft grid with NDIM != 3",1);
284 }
285 } else {
286 MADNESS_EXCEPTION("no such grid type",1);
287 }
288 }
289
290
291 /// ctor takes centers of the grid and the grid builder
292 molecular_grid(const std::vector<Vector<double,NDIM>> origins, std::shared_ptr<gridbase> grid)
293 : centers(origins), grid_builder(grid) {
294 if (centers.size()==0) centers.push_back({0,0,0});
295 }
296
297 /// ctor takes molecule and grid builder
298 molecular_grid(const Molecule& molecule, std::shared_ptr<gridbase> grid) : molecular_grid(molecule.get_all_coords_vec(),grid) {}
299
300 std::vector<Vector<double,NDIM>> get_grid() const {
301 MADNESS_CHECK_THROW(grid_builder,"no grid builder given in molecular_grid");
302 MADNESS_CHECK_THROW(centers.size()>0,"no centers given in molecular_grid");
303 std::vector<Vector<double,NDIM>> grid;
304 for (const auto& coords : centers) {
305 print("atom sites",coords);
306 if (auto builder=dynamic_cast<dftgrid<NDIM>*>(grid_builder.get())) {
307 if constexpr (NDIM==3) {
308 dftgrid<NDIM> b1(builder->builder.get_nradial(),builder->builder.get_angular_order(),coords);
309 auto atomgrid=b1.get_grid();
310 grid.insert(grid.end(),atomgrid.begin(),atomgrid.end());
311 } else {
312 MADNESS_EXCEPTION("no DFT grid for NDIM /= 3",1);
313 }
314 } else if (auto builder=dynamic_cast<randomgrid<NDIM>*>(grid_builder.get())) {
315 randomgrid<NDIM> b1(builder->get_volume_element(),builder->get_radius(),coords);
316 auto atomgrid=b1.get_grid();
317 grid.insert(grid.end(),atomgrid.begin(),atomgrid.end());
318 } else {
319 MADNESS_EXCEPTION("no such grid builder",1);
320 }
321 }
322 return grid;
323 }
324
325 private:
326 std::vector<Vector<double,NDIM>> centers;
327 std::shared_ptr<gridbase> grid_builder;
328
329 };
330
331template<std::size_t PDIM>
332struct particle {
333 std::array<int,PDIM> dims;
334
335 /// default constructor
336 particle() = default;
337
338 /// convenience for particle 1 (the left/first particle)
340 particle p;
341 for (int i=0; i<PDIM; ++i) p.dims[i]=i;
342 return p;
343 }
344 /// convenience for particle 2 (the right/second particle)
346 particle p;
347 for (int i=0; i<PDIM; ++i) p.dims[i]=i+PDIM;
348 return p;
349 }
350
351 /// return the other particle
354 if (is_first()) return particle2();
355 return particle1();
356 }
357
358 particle(const int p) : particle(std::vector<int>(1,p)) {}
359 particle(const int p1, const int p2) : particle(std::vector<int>({p1,p2})) {}
360 particle(const int p1, const int p2,const int p3) : particle(std::vector<int>({p1,p2,p3})) {}
361 particle(const std::vector<int> p) {
362 for (int i=0; i<PDIM; ++i) dims[i]=p[i];
363 }
364
365 std::string str() const {
366 std::stringstream ss;
367 ss << *this;
368 return ss.str();
369 }
370
371
372 /// type conversion to std::array
373 std::array<int,PDIM> get_array() const {
374 return dims;
375 }
376
377
378 /// assuming two particles only
379 bool is_first() const {return dims[0]==0;}
380 /// assuming two particles only
381 bool is_last() const {return dims[0]==(PDIM);}
382
383 template<std::size_t DUMMYDIM=PDIM>
384 typename std::enable_if_t<DUMMYDIM==1, std::tuple<int>>
385 get_tuple() const {return std::tuple<int>(dims[0]);}
386
387 template<std::size_t DUMMYDIM=PDIM>
388 typename std::enable_if_t<DUMMYDIM==2, std::tuple<int,int>>
389 get_tuple() const {return std::tuple<int,int>(dims[0],dims[1]);}
390
391 template<std::size_t DUMMYDIM=PDIM>
392 typename std::enable_if_t<DUMMYDIM==3, std::tuple<int,int,int>>
393 get_tuple() const {return std::tuple<int,int,int>(dims[0],dims[1],dims[2]);}
394};
395
396template<std::size_t PDIM>
397std::ostream& operator<<(std::ostream& os, const particle<PDIM>& p) {
398 os << "(";
399 for (auto i=0; i<PDIM-1; ++i) os << p.dims[i] << ";";
400 os << p.dims[PDIM-1] << ")";
401 return os;
402}
403
404/// the low-rank functor is what the LowRankFunction will represent
405
406/// derive from this class :
407/// must implement in inner product
408/// may implement an operator()(const coord_nd&)
409template<typename T, std::size_t NDIM, std::size_t LDIM=NDIM/2>
411
412 virtual ~LRFunctorBase() {};
413 virtual std::vector<Function<T,LDIM>> inner(const std::vector<Function<T,LDIM>>& rhs,
414 const particle<LDIM> p1, const particle<LDIM> p2) const =0;
415
416 virtual Function<T,LDIM> inner(const Function<T,LDIM>& rhs, const particle<LDIM> p1, const particle<LDIM> p2) const {
417 return inner(std::vector<Function<T,LDIM>>({rhs}),p1,p2)[0];
418 }
419
420 virtual T operator()(const Vector<T,NDIM>& r) const =0;
421 virtual typename Tensor<T>::scalar_type norm2() const {
422 MADNESS_EXCEPTION("L2 norm not implemented",1);
423 }
424
425 virtual World& world() const =0;
426 friend std::vector<Function<T,LDIM>> inner(const LRFunctorBase& functor, const std::vector<Function<T,LDIM>>& rhs,
427 const particle<LDIM> p1, const particle<LDIM> p2) {
428 return functor.inner(rhs,p1,p2);
429 }
431 const particle<LDIM> p1, const particle<LDIM> p2) {
432 return functor.inner(rhs,p1,p2);
433 }
434
435};
436
437template<typename T, std::size_t NDIM, std::size_t LDIM=NDIM/2>
438struct LRFunctorF12 : public LRFunctorBase<T,NDIM> {
439// LRFunctorF12() = default;
440 LRFunctorF12(const std::shared_ptr<SeparatedConvolution<T,LDIM>> f12, const std::vector<Function<T,LDIM>>& a,
441 const std::vector<Function<T,LDIM>>& b) : f12(f12), a(a), b(b) {
442
443 // if a or b are missing, they are assumed to be 1
444 // you may not provide a or b, but if you do they have to have the same size because they are summed up
445 if (a.size()>0 and b.size()>0)
446 MADNESS_CHECK_THROW(a.size()==b.size(), "a and b must have the same size");
447 if (a.size()==0) this->a.resize(b.size());
448 if (b.size()==0) this->b.resize(a.size());
449 MADNESS_CHECK_THROW(this->a.size()==this->b.size(), "a and b must have the same size");
450 }
451
452 /// delegate to the other ctor with vector arguments
454 const Function<T,LDIM>& a, const Function<T,LDIM>& b)
455 : LRFunctorF12(f12,std::vector<Function<T,LDIM>>({a}), std::vector<Function<T,LDIM>>({b})) {}
456
457
458private:
459 std::shared_ptr<SeparatedConvolution<T,LDIM>> f12; ///< a two-particle function
460 std::vector<Function<T,LDIM>> a,b; ///< the lo-dim functions
461public:
462
463 World& world() const {return f12->get_world();}
464 std::vector<Function<T,LDIM>> inner(const std::vector<Function<T,LDIM>>& rhs,
465 const particle<LDIM> p1, const particle<LDIM> p2) const {
466
467 std::vector<Function<T,LDIM>> result;
468 // functor is now \sum_i a_i(1) b_i(2) f12
469 // result(1) = \sum_i \int a_i(1) f(1,2) b_i(2) rhs(2) d2
470 // = \sum_i a_i(1) \int f(1,2) b_i(2) rhs(2) d2
471 World& world=rhs.front().world();
472
473 const int nbatch=30;
474 for (int i=0; i<rhs.size(); i+=nbatch) {
475 std::vector<Function<T,LDIM>> rhs_batch;
476 auto begin= rhs.begin()+i;
477 auto end= (i+nbatch)<rhs.size() ? rhs.begin()+i+nbatch : rhs.end();
478 std::copy(begin,end, std::back_inserter(rhs_batch));
480
481 if (a.size()==0) tmp2=apply(world,*(f12),rhs_batch);
482
483 for (int ia=0; ia<a.size(); ia++) {
484 auto premultiply= p1.is_first() ? a[ia] : b[ia];
485 auto postmultiply= p1.is_first() ? b[ia] : a[ia];
486
487 auto tmp=copy(world,rhs_batch);
488 if (premultiply.is_initialized()) tmp=rhs_batch*premultiply;
489 auto tmp1=apply(world,*(f12),tmp);
490 if (postmultiply.is_initialized()) tmp1=tmp1*postmultiply;
491 tmp2+=tmp1;
492 }
493
494 for (auto& t : tmp2) result.push_back(t);
495 }
496 return result;
497 }
498
499 typename Tensor<T>::scalar_type norm2() const {
501 [](const Vector<double, LDIM>& r) { return 1.0; });
502 std::vector<Function<T, LDIM>> pre, post;
503 std::size_t sz = a.size();
504 if (sz == 0) {
505 pre = std::vector<Function<T, LDIM>>(1, one);
506 post = std::vector<Function<T, LDIM>>(1, one);
507 } else {
508 pre = (a.front().is_initialized()) ? a : std::vector<Function<T, LDIM>>(sz, one);
509 post = (b.front().is_initialized()) ? b : std::vector<Function<T, LDIM>>(sz, one);
510 }
511
514
515 // \int f(1,2)^2 d1d2 = \int f(1,2)^2 pre(1)^2 post(2)^2 d1 d2
516 // || \sum_i f(1,2) a_i(1) b_i(2) || = \int ( \sum_{ij} a_i(1) a_j(1) f(1,2)^2 b_i(1) b_j(2) ) d1d2
517 std::vector<Function<T,LDIM>> aa,bb;
518 for (std::size_t i=0; i<pre.size(); ++i) {
519 aa=append(aa,(pre[i]*pre));
520 bb=append(bb,(post[i]*post));
521 }
522// typename Tensor<T>::scalar_type term1 =madness::inner(mul(world(),post,post),f12sq(mul(world(),pre,pre)));
524 return sqrt(term1);
525
526 }
527
529
530 if (a.size()==0) return 0.0;
531 auto split = [](const Vector<double,NDIM>& r) {
533 for (int i=0; i<LDIM; ++i) {
534 first[i]=r[i];
535 second[i]=r[i+LDIM];
536 }
537 return std::make_pair(first,second);
538 };
539
540 double gamma=f12->info.mu;
541 auto [first,second]=split(r);
542
543
544 double result=0.0;
545 for (std::size_t ia=0; ia<a.size(); ++ia) {
546 double result1=1.0;
547 if (a[ia].is_initialized()) result1*=a[ia](first);
548 if (b[ia].is_initialized()) result1*=b[ia](first);
549 if (f12->info.type==OT_SLATER) result1*=exp(-gamma*(first-second).normf());
550 else if (f12->info.type==OT_GAUSS) result1*=exp(-gamma* madness::inner(first-second,first-second));
551 else {
552 MADNESS_EXCEPTION("no such operator_type",1);
553 }
554 result+=result1;
555 }
556 return result;
557
558 }
559};
560
561template<typename T, std::size_t NDIM, std::size_t LDIM=NDIM/2>
562struct LRFunctorPure : public LRFunctorBase<T,NDIM> {
563 LRFunctorPure() = default;
565 World& world() const {return f.world();}
566
567 Function<T, NDIM> f; ///< a hi-dim function
568
569 std::vector<Function<T,LDIM>> inner(const std::vector<Function<T,LDIM>>& rhs,
570 const particle<LDIM> p1, const particle<LDIM> p2) const {
572// std::vector<Function<T,LDIM>> result;
573// for (const auto& r : rhs) result.push_back(madness::inner(f,r,p1.get_tuple(),p2.get_tuple()));
574// return result;
575 }
576
578 return f(r);
579 }
580
581 typename Tensor<T>::scalar_type norm2() const {
582 return f.norm2();
583 }
584};
585
586
587 /// LowRankFunction represents a hi-dimensional (NDIM) function as a sum of products of low-dimensional (LDIM) functions
588
589 /// f(1,2) = \sum_i g_i(1) h_i(2)
590 /// a LowRankFunction can be created from a hi-dim function directly, or from a composite like f(1,2) phi(1) psi(2),
591 /// where f(1,2) is a two-particle function (e.g. a Slater function)
592 template<typename T, std::size_t NDIM, std::size_t LDIM=NDIM/2>
594 public:
595
597 double rank_revealing_tol=1.e-8; // rrcd tol
598 std::string orthomethod="canonical";
599 bool do_print=false;
600 std::vector<Function<T,LDIM>> g,h;
603
605
607 double tol, std::string orthomethod) : world(g.front().world()),
609
610 }
611
612 /// shallow copy ctor
615 g(other.g), h(other.h) {
616 }
617
618 /// deep copy
619 friend LowRankFunction copy(const LowRankFunction& other) {
621 }
622
623 LowRankFunction& operator=(const LowRankFunction& f) { // Assignment required for storage in vector
625 std::swap(ff.g,g);
626 std::swap(ff.h,h);
627 return *this;
628 }
629
630 /// function evaluation
633 for (int i=0; i<LDIM; ++i) {
634 first[i]=r[i];
635 second[i]=r[i+LDIM];
636 }
637 double result=0.0;
638 for (int i=0; i<rank(); ++i) result+=g[i](first)*h[i](second);
639 return result;
640 }
641
642 /*
643 * arithmetic section
644 */
645
646 /// addition
648 LowRankFunction<T,NDIM> result=copy(*this);
649 result+=b;
650 return result;
651 }
652 /// subtraction
654 LowRankFunction<T,NDIM> result=copy(*this);
655 result-=b;
656 return result;
657 }
658
659 /// in-place addition
661
662 g=append(g,copy(b.g));
663 h=append(h,copy(b.h));
664 return *this;
665 }
666
667 /// in-place subtraction
669 g=append(g,-1.0*b.g); // operator* implies deep copy of b.g
670 h=append(h,copy(b.h));
671 return *this;
672 }
673
674 /// scale by a scalar
675 template<typename Q>
679
680 /// out-of-place scale by a scalar (no type conversion)
684
685 /// multiplication with a scalar
686 friend LowRankFunction operator*(const T a, const LowRankFunction& other) {
687 return other*a;
688 }
689
690 /// in-place scale by a scalar (no type conversion)
692 g=g*a;
693 return *this;
694 }
695
696 /// l2 norm
698 auto tmp1=matrix_inner(world,h,h);
699 auto tmp2=matrix_inner(world,g,g);
700 return sqrt(tmp1.trace(tmp2));
701 }
702
703 std::vector<Function<T,LDIM>> get_functions(const particle<LDIM>& p) const {
704 MADNESS_CHECK(p.is_first() or p.is_last());
705 if (p.is_first()) return g;
706 return h;
707 }
708
709 std::vector<Function<T,LDIM>> get_g() const {return g;}
710 std::vector<Function<T,LDIM>> get_h() const {return h;}
711
712 long rank() const {return g.size();}
713
714 /// return the size in GByte
715 double size() const {
716 double sz=get_size(world,g);
717 sz+=get_size(world,h);
718 return sz;
719 }
720
722 auto fapprox=hartree_product(g[0],h[0]);
723 for (int i=1; i<g.size(); ++i) fapprox+=hartree_product(g[i],h[i]);
724 return fapprox;
725 }
726
727 /// orthonormalize the argument vector
728 std::vector<Function<T,LDIM>> orthonormalize(const std::vector<Function<T,LDIM>>& g) const {
729
730 double tol=rank_revealing_tol;
731 std::vector<Function<T,LDIM>> g2;
732 auto ovlp=matrix_inner(world,g,g);
733 if (orthomethod=="canonical") {
734 tol*=0.01;
735 print("orthonormalizing with method/tol",orthomethod,tol);
737 } else if (orthomethod=="cholesky") {
738 print("orthonormalizing with method/tol",orthomethod,tol);
739 g2=orthonormalize_rrcd(g,ovlp,tol);
740 }
741 else {
742 MADNESS_EXCEPTION("no such orthomethod",1);
743 }
745 return truncate(g2,tight_thresh);
746 }
747
748
749 /// optimize the lrf using the lrfunctor
750
751 /// @param[in] nopt number of iterations (wrt to Alg. 4.3 in Halko)
752 void optimize(const LRFunctorBase<T,NDIM>& lrfunctor1, const long nopt=1) {
753 timer t(world);
755 for (int i=0; i<nopt; ++i) {
756 // orthonormalize h
758 t.tag("ortho1");
760 t.tag("inner1");
762 t.tag("ortho2");
764 t.tag("inner2");
765 }
766 }
767
768 /// remove linear dependencies without orthonormalization
770
771 // use rank-revealing cholesky decomposition to remove linear dependencies
772
773
774 }
775
776 /// after external operations g might not be orthonormal and/or optimal -- reorthonormalize
777
778 /// orthonormalization similar to Bischoff, Harrison, Valeev, JCP 137 104103 (2012), Sec II C 3
779 /// f =\sum_i g_i h_i
780 /// = g X- (X+)^T (Y+)^T Y- h
781 /// = g X- U S V^T Y- h
782 /// = g (X- U) (S V^T Y-) h
783 /// requires 2 matrix_inner and 2 transforms. g and h are optimal, but contain all cusps etc..
784 /// @param[in] thresh SVD threshold
785 void reorthonormalize(double thresh=-1.0) {
789 auto [eval_g, evec_g] = syev(ovlp_g);
790 auto [eval_h, evec_h] = syev(ovlp_h);
791
792 // get relevant part of the eigenvalues and eigenvectors
793 // eigenvalues are sorted in ascending order
794 auto get_slice = [](auto eval, double thresh) {
795 // remove small/negative eigenvalues
796 eval.screen(thresh);
797 Slice s;
798 for (int i=0; i<eval.size(); ++i) {
799 MADNESS_CHECK_THROW(eval[i]>=0.0,"negative eigenvalues in reorthonormalize");
800 if (eval[i]>thresh) {
801 return s=Slice(i,-1); // from i to the end
802 break;
803 }
804 }
805 return s;
806 };
807
810
817
818 for (int i=0; i<eval_g.size(); ++i) Xplus(_,i)*=std::pow(eval_g(i),0.5);
819 for (int i=0; i<eval_g.size(); ++i) Xminus(_,i)*=std::pow(eval_g(i),-0.5);
820 for (int i=0; i<eval_h.size(); ++i) Yplus(_,i)*=std::pow(eval_h(i),0.5);
821 for (int i=0; i<eval_h.size(); ++i) Yminus(_,i)*=std::pow(eval_h(i),-0.5);
822
823 Tensor<T> M=madness::inner(Xplus,Yplus,0,0); // (X+)^T Y+
824 auto [U,s,VT]=svd(M);
825
826 // truncate
828 int i=s.size()-1;
829 for (;i>=0; i--) {
830 s_accumulated+=s[i];
831 if (s_accumulated>thresh) {
832 i++;
833 break;
834 }
835 }
836 for (int j=0; j<s.size(); ++j) VT(j,_)*=s[j];
839
842 }
843
844
845 double check_orthonormality(const std::vector<Function<T,LDIM>>& v) const {
848 }
849
850 double check_orthonormality(const Tensor<T>& ovlp) const {
851 timer t(world);
854 for (int i=0; i<ovlp2.dim(0); ++i) ovlp2(i,i)-=1.0;
855 if (world.rank()==0 and do_print) {
856 print("absmax",ovlp2.absmax());
857 print("l2",ovlp2.normf()/ovlp2.size());
858 }
859 t.tag("check_orthonoramality");
860 return ovlp.absmax();
861 }
862
863 /// compute the l2 error |functor - \sum_i g_ih_i|_2
864
865 /// \int (f(1,2) - gh(1,2))^2 = \int f(1,2)^2 - 2\int f(1,2) gh(1,2) + \int gh(1,2)^2
866 /// since we are subtracting large numbers the numerics are sensitive, and NaN may be returned..
868
869 timer t(world);
871
872 // \int f(1,2)^2 d1d2
873 double term1 = lrfunctor1.norm2();
875 t.tag("computing term1");
876
877 // \int f(1,2) pre(1) post(2) \sum_i g(1) h(2) d1d2
878// double term2=madness::inner(pre*g,f12(post*h));
880 t.tag("computing term2");
881
882 // g functions are orthonormal
883 // \int gh(1,2)^2 d1d2 = \int \sum_{ij} g_i(1) g_j(1) h_i(2) h_j(2) d1d2
884 // = \sum_{ij} \int g_i(1) g_j(1) d1 \int h_i(2) h_j(2) d2
885 // = \sum_{ij} delta_{ij} \int h_i(2) h_j(2) d2
886 // = \sum_{i} \int h_i(2) h_i(2) d2
887 double zero=check_orthonormality(g);
888 if (zero>1.e-10) print("g is not orthonormal",zero);
889 // double term3a=madness::inner(h,h);
890 auto tmp1=matrix_inner(world,h,h);
891 auto tmp2=matrix_inner(world,g,g);
892 double term3=tmp1.trace(tmp2);
893// print("term3/a/diff",term3a,term3,term3-term3a);
894 t.tag("computing term3");
895
896 double arg=term1-2.0*term2+term3;
897 if (arg<0.0) {
898 print("negative l2 error");
899 arg*=-1.0;
900// throw std::runtime_error("negative argument in l2error");
901 }
902 double error=sqrt(arg)/sqrt(term1);
903 if (world.rank()==0 and do_print) {
904 print("term1,2,3, error",term1, term2, term3, " --",error);
905 }
906
907 return error;
908 }
909
910 };
911
912 // This interface is necessary to compute inner products
913 template<typename T, std::size_t NDIM>
915 World& world=a.world;
916 return (matrix_inner(world,a.g,b.g).emul(matrix_inner(world,a.h,b.h))).sum();
917 }
918
919
920
921// template<typename T, std::size_t NDIM>
922// LowRankFunction<T,NDIM> inner(const Function<T,NDIM>& lhs, const LowRankFunction<T,NDIM>& rhs,
923// const std::tuple<int> v1, const std::tuple<int> v2) {
924// World& world=rhs.world;
925// // int lhs(1,2) rhs(2,3) d2 = \sum \int lhs(1,2) g_i(2) h_i(3) d2
926// // = \sum \int lhs(1,2) g_i(2) d2 h_i(3)
927// LowRankFunction<T, NDIM + NDIM - 2> result(world);
928// result.h=rhs.h;
929// decltype(rhs.g) g;
930// for (int i=0; i<rhs.rank(); ++i) {
931// g.push_back(inner(lhs,rhs.g[i],{v1},{0}));
932// }
933// result.g=g;
934// return result;
935// }
936
937 /**
938 * inner product: LowRankFunction lrf; Function f, g; double d
939 * lrf(1,3) = inner(lrf(1,2), lrf(2,3))
940 * lrf(1,3) = inner(lrf(1,2), f(2,3))
941 * g(1) = inner(lrf(1,2), f(2))
942 * d = inner(lrf(1,2), f(1,2))
943 * d = inner(lrf(1,2), lrf(1,2))
944 */
945
946 /// lrf(1,3) = inner(full(1,2), lrf(2,3))
947
948 /// @param[in] f1 the first function
949 /// @param[in] f2 the second function
950 /// @param[in] p1 the integration variable of the first function
951 /// @param[in] p2 the integration variable of the second function
952 template<typename T, std::size_t NDIM, std::size_t PDIM>
954 const particle<PDIM> p1, const particle<PDIM> p2) {
955 auto result=inner(f2,f1,p2,p1);
956 std::swap(result.g,result.h);
957 return result;
958 }
959
960 /// lrf(1,3) = inner(lrf(1,2), full(2,3))
961
962 /// @param[in] f1 the first function
963 /// @param[in] f2 the second function
964 /// @param[in] p1 the integration variable of the first function
965 /// @param[in] p2 the integration variable of the second function
966 template<typename T, std::size_t NDIM, std::size_t PDIM>
968 const particle<PDIM> p1, const particle<PDIM> p2) {
969 static_assert(TensorTypeData<T>::iscomplex==false, "complex inner in LowRankFunction not implemented");
970 World& world=f1.world;
971 static_assert(2*PDIM==NDIM);
972 // int f(1,2) k(2,3) d2 = \sum \int g_i(1) h_i(2) k(2,3) d2
973 // = \sum g_i(1) \int h_i(2) k(2,3) d2
974 LowRankFunction<T, NDIM> result(world);
975 if (p1.is_last()) { // integrate over 2: result(1,3) = lrf(1,2) f(2,3)
976 result.g = f1.g;
979 } else if (p1.is_first()) { // integrate over 1: result(2,3) = lrf(1,2) f(1,3)
980 result.g = f1.h; // correct! second variable of f1 becomes first variable of result
983 }
984 return result;
985 }
986
987 /// lrf(1,3) = inner(lrf(1,2), lrf(2,3))
988
989 /// @param[in] f1 the first function
990 /// @param[in] f2 the second function
991 /// @param[in] p1 the integration variable of the first function
992 /// @param[in] p2 the integration variable of the second function
993 template<typename T, std::size_t NDIM, std::size_t PDIM>
995 const particle<PDIM> p1, const particle<PDIM> p2) {
996 World& world=f1.world;
997 static_assert(2*PDIM==NDIM);
998
999 // inner(lrf(1,2) ,lrf(2,3) ) = \sum_ij g1_i(1) <h1_i(2) g2_j(2)> h2_j(3)
1000 auto matrix=matrix_inner(world,f2.get_functions(p2),f1.get_functions(p1));
1001 auto htilde=transform(world,f2.get_functions(p2.complement()),matrix);
1002 auto gg=copy(world,f1.get_functions(p1.complement()));
1003 return LowRankFunction<T,NDIM>(gg,htilde,f1.rank_revealing_tol,f1.orthomethod);
1004 }
1005
1006 /// f(1) = inner(lrf(1,2), f(2))
1007
1008 /// @param[in] f1 the first function
1009 /// @param[in] vf vector of the second functions
1010 /// @param[in] p1 the integration variable of the first function
1011 /// @param[in] p2 the integration variable of the second function, dummy variable for consistent notation
1012 template<typename T, std::size_t NDIM, std::size_t PDIM>
1013 std::vector<Function<T,NDIM-PDIM>> inner(const LowRankFunction<T,NDIM>& f1, const std::vector<Function<T,PDIM>>& vf,
1015 World& world=f1.world;
1016 static_assert(2*PDIM==NDIM);
1017 MADNESS_CHECK(p2.is_first());
1018
1019 // inner(lrf(1,2), f_k(2) ) = \sum_i g1_i(1) <h1_i(2) f_k(2)>
1020 auto matrix=matrix_inner(world,f1.get_functions(p1),vf);
1021 return transform(world,f1.get_functions(p1.complement()),matrix);
1022 }
1023
1024 /// f(1) = inner(lrf(1,2), f(2))
1025
1026 /// @param[in] f1 the first function
1027 /// @param[in] vf the second function
1028 /// @param[in] p1 the integration variable of the first function
1029 /// @param[in] p2 the integration variable of the second function, dummy variable for consistent notation
1030 template<typename T, std::size_t NDIM, std::size_t PDIM>
1033 return inner(f1,std::vector<Function<T,PDIM>>({f2}),p1,p2)[0];
1034 }
1035
1036 template<typename T, std::size_t NDIM, std::size_t LDIM=NDIM/2>
1038 public:
1039
1042
1044 std::vector<Vector<double,LDIM>> origins; ///< origins of the molecular grid
1045
1049
1052
1054
1055 LowRankFunctionFactory& set_radius(const double radius) {
1056 parameters.set_user_defined_value("radius",radius);
1057 return *this;
1058 }
1059 LowRankFunctionFactory& set_volume_element(const double volume_element) {
1060 parameters.set_user_defined_value("volume_element",volume_element);
1061 return *this;
1062 }
1065 return *this;
1066 }
1067 LowRankFunctionFactory& set_orthomethod(const std::string orthomethod) {
1068 parameters.set_user_defined_value("orthomethod",orthomethod);
1069 return *this;
1070 }
1071
1073 World& world=lrfunctor.world();
1074 bool do_print=true;
1075 timer t1(world);
1076 t1.do_print=do_print;
1077 auto orthomethod=parameters.orthomethod();
1078 auto rank_revealing_tol=parameters.tol();
1079
1080 // get sampling grid
1082 auto grid=mgrid.get_grid();
1083 if (world.rank()==0) print("grid size",grid.size());
1084
1086 t1.tag("Yforming");
1087
1088 auto ovlp=matrix_inner(world,Y,Y); // error in symmetric matrix_inner, use non-symmetric form here!
1089 t1.tag("compute ovlp");
1090 auto g=truncate(orthonormalize_rrcd(Y,ovlp,rank_revealing_tol));
1091 t1.tag("rrcd/truncate/thresh");
1092 auto sz=get_size(world,g);
1093 if (world.rank()==0 and do_print) print("gsize",sz);
1094// check_orthonormality(g);
1095
1096 if (world.rank()==0 and do_print) {
1097 print("Y.size()",Y.size());
1098 print("g.size()",g.size());
1099 }
1100
1101 auto h=truncate(inner(lrfunctor,g,p1,p1));
1102 t1.tag("Y backprojection with truncation");
1104
1105 }
1106
1107 /// apply a rhs (delta or exponential) on grid points to the hi-dim function and form Y = A_ij w_j (in Halko's language)
1108 std::vector<Function<T,LDIM>> Yformer(const LRFunctorBase<T,NDIM>& lrfunctor1, const std::vector<Vector<double,LDIM>>& grid,
1109 const std::string rhsfunctiontype, const double exponent=30.0) const {
1110
1111 World& world=lrfunctor1.world();
1112 std::vector<Function<double,LDIM>> Y;
1113 if (rhsfunctiontype=="exponential") {
1114 std::vector<Function<double,LDIM>> omega;
1115 double coeff=std::pow(2.0*exponent/constants::pi,0.25*LDIM);
1116 for (const auto& point : grid) {
1117 omega.push_back(FunctionFactory<double,LDIM>(world)
1118 .functor([&point,&exponent,&coeff](const Vector<double,LDIM>& r)
1119 {
1120 auto r_rel=r-point;
1121 return coeff*exp(-exponent*madness::inner(r_rel,r_rel));
1122 }));
1123 }
1125 } else {
1126 MADNESS_EXCEPTION("confused rhsfunctiontype",1);
1127 }
1128 auto norms=norm2s(world,Y);
1129 std::vector<Function<double,LDIM>> Ynormalized;
1130
1131 for (int i=0; i<Y.size(); ++i) if (norms[i]>parameters.tol()) Ynormalized.push_back(Y[i]);
1132 normalize(world,Ynormalized);
1133 return Ynormalized;
1134 }
1135
1136 };
1137
1138
1139} // namespace madness
1140
1141#endif //MADNESS_LOWRANKFUNCTION_H
Definition IntegratorXX.h:29
std::vector< madness::Vector< double, 3 > > get_points() const
get the grid points
Definition IntegratorXX.h:105
void set_nradial(const std::size_t nradial1)
number of radial gridpoints on the interval [0,\inf)
Definition IntegratorXX.h:69
void set_angular_order(const std::size_t order)
a quadrature of a give order will integrate a spherical harmonic of that order exactly
Definition IntegratorXX.h:59
std::size_t get_angular_order() const
a quadrature of a give order will integrate a spherical harmonic of that order exactly
Definition IntegratorXX.h:64
void set_origin(const madness::Vector< double, 3 > &o)
the origin/center of the grid
Definition IntegratorXX.h:79
std::size_t get_nradial() const
number of radial gridpoints on the interval [0,\inf)
Definition IntegratorXX.h:74
void make_grid()
Definition IntegratorXX.h:83
static const double & get_thresh()
Returns the default threshold.
Definition funcdefaults.h:176
static const Tensor< double > & get_cell()
Gets the user cell for the simulation.
Definition funcdefaults.h:347
FunctionFactory implements the named-parameter idiom for Function.
Definition function_factory.h:86
FunctionFactory & functor(const std::shared_ptr< FunctionFunctorInterface< T, NDIM > > &f)
Definition function_factory.h:141
A multiresolution adaptive numerical function.
Definition mra.h:139
Definition lowrankfunction.h:1037
std::vector< Function< T, LDIM > > Yformer(const LRFunctorBase< T, NDIM > &lrfunctor1, const std::vector< Vector< double, LDIM > > &grid, const std::string rhsfunctiontype, const double exponent=30.0) const
apply a rhs (delta or exponential) on grid points to the hi-dim function and form Y = A_ij w_j (in Ha...
Definition lowrankfunction.h:1108
LowRankFunctionFactory(const LowRankFunctionParameters param, const Molecule &molecule)
Definition lowrankfunction.h:1050
const particle< LDIM > p2
Definition lowrankfunction.h:1041
LowRankFunctionFactory & set_radius(const double radius)
Definition lowrankfunction.h:1055
const particle< LDIM > p1
Definition lowrankfunction.h:1040
LowRankFunctionFactory & set_rank_revealing_tol(const double rrtol)
Definition lowrankfunction.h:1063
LowRankFunctionFactory(const LowRankFunctionParameters param, const std::vector< Vector< double, LDIM > > origins={})
Definition lowrankfunction.h:1047
LowRankFunction< T, NDIM > project(const LRFunctorBase< T, NDIM > &lrfunctor) const
Definition lowrankfunction.h:1072
LowRankFunctionParameters parameters
Definition lowrankfunction.h:1043
LowRankFunctionFactory & set_volume_element(const double volume_element)
Definition lowrankfunction.h:1059
std::vector< Vector< double, LDIM > > origins
origins of the molecular grid
Definition lowrankfunction.h:1044
LowRankFunctionFactory & set_orthomethod(const std::string orthomethod)
Definition lowrankfunction.h:1067
LowRankFunctionFactory(const LowRankFunctionFactory &other)=default
LowRankFunction represents a hi-dimensional (NDIM) function as a sum of products of low-dimensional (...
Definition lowrankfunction.h:593
void optimize(const LRFunctorBase< T, NDIM > &lrfunctor1, const long nopt=1)
optimize the lrf using the lrfunctor
Definition lowrankfunction.h:752
friend LowRankFunction operator*(const T a, const LowRankFunction &other)
multiplication with a scalar
Definition lowrankfunction.h:686
double rank_revealing_tol
Definition lowrankfunction.h:597
std::vector< Function< T, LDIM > > get_h() const
Definition lowrankfunction.h:710
std::vector< Function< T, LDIM > > get_functions(const particle< LDIM > &p) const
Definition lowrankfunction.h:703
std::string orthomethod
Definition lowrankfunction.h:598
LowRankFunction operator-(const LowRankFunction &b) const
subtraction
Definition lowrankfunction.h:653
LowRankFunction operator*(const Q a) const
scale by a scalar
Definition lowrankfunction.h:676
std::vector< Function< T, LDIM > > orthonormalize(const std::vector< Function< T, LDIM > > &g) const
orthonormalize the argument vector
Definition lowrankfunction.h:728
double check_orthonormality(const std::vector< Function< T, LDIM > > &v) const
Definition lowrankfunction.h:845
LowRankFunction operator*(const T a) const
out-of-place scale by a scalar (no type conversion)
Definition lowrankfunction.h:681
TensorTypeData< T >::scalar_type norm2() const
l2 norm
Definition lowrankfunction.h:697
long rank() const
Definition lowrankfunction.h:712
LowRankFunction(const LowRankFunction &other)
shallow copy ctor
Definition lowrankfunction.h:613
bool do_print
Definition lowrankfunction.h:599
void reorthonormalize(double thresh=-1.0)
after external operations g might not be orthonormal and/or optimal – reorthonormalize
Definition lowrankfunction.h:785
LowRankFunction & operator*=(const T a)
in-place scale by a scalar (no type conversion)
Definition lowrankfunction.h:691
LowRankFunction & operator=(const LowRankFunction &f)
Definition lowrankfunction.h:623
const particle< LDIM > p2
Definition lowrankfunction.h:602
Function< T, NDIM > reconstruct() const
Definition lowrankfunction.h:721
std::vector< Function< T, LDIM > > get_g() const
Definition lowrankfunction.h:709
const particle< LDIM > p1
Definition lowrankfunction.h:601
LowRankFunction & operator-=(const LowRankFunction &b)
in-place subtraction
Definition lowrankfunction.h:668
LowRankFunction(std::vector< Function< T, LDIM > > g, std::vector< Function< T, LDIM > > h, double tol, std::string orthomethod)
Definition lowrankfunction.h:606
std::vector< Function< T, LDIM > > h
Definition lowrankfunction.h:600
LowRankFunction & operator+=(const LowRankFunction &b)
in-place addition
Definition lowrankfunction.h:660
double l2error(const LRFunctorBase< T, NDIM > &lrfunctor1) const
compute the l2 error |functor - \sum_i g_ih_i|_2
Definition lowrankfunction.h:867
LowRankFunction(World &world)
Definition lowrankfunction.h:604
LowRankFunction operator+(const LowRankFunction &b) const
addition
Definition lowrankfunction.h:647
double check_orthonormality(const Tensor< T > &ovlp) const
Definition lowrankfunction.h:850
std::vector< Function< T, LDIM > > g
Definition lowrankfunction.h:600
World & world
Definition lowrankfunction.h:596
friend LowRankFunction copy(const LowRankFunction &other)
deep copy
Definition lowrankfunction.h:619
void remove_linear_depdencies(double thresh=-1.0)
remove linear dependencies without orthonormalization
Definition lowrankfunction.h:769
double size() const
return the size in GByte
Definition lowrankfunction.h:715
T operator()(const Vector< double, NDIM > &r) const
function evaluation
Definition lowrankfunction.h:631
Definition molecule.h:124
class for holding the parameters for calculation
Definition QCCalculationParametersBase.h:290
virtual void read_input_and_commandline_options(World &world, const commandlineparser &parser, const std::string tag)
Definition QCCalculationParametersBase.h:325
void set_user_defined_value(const std::string &key, const T &value)
Definition QCCalculationParametersBase.h:533
Convolutions in separated form (including Gaussian)
Definition operator.h:139
friend SeparatedConvolution< Q, NDIM > combine(const std::shared_ptr< SeparatedConvolution< Q, NDIM > > left, const std::shared_ptr< SeparatedConvolution< Q, NDIM > > right)
combine 2 convolution operators to one
Definition operator.h:1754
A slice defines a sub-range or patch of a dimension.
Definition slice.h:103
Traits class to specify support of numeric types.
Definition type_data.h:56
A tensor is a multidimensional array.
Definition tensor.h:317
TensorTypeData< T >::scalar_type scalar_type
C++ typename of the real type associated with a complex type.
Definition tensor.h:409
A simple, fixed dimension vector.
Definition vector.h:64
constexpr size_type size() const
Accessor for the number of elements in the Vector.
Definition vector.h:292
constexpr void fill(const T &t)
Fill the Vector with the specified value.
Definition vector.h:371
A parallel world class.
Definition world.h:132
ProcessID rank() const
Returns the process rank in this World (same as MPI_Comm_rank()).
Definition world.h:320
ProcessID size() const
Returns the number of processes in this World (same as MPI_Comm_size()).
Definition world.h:330
Definition lowrankfunction.h:81
std::vector< Vector< double, NDIM > > get_grid() const
Definition lowrankfunction.h:113
dftgrid(const double volume_element, const double radius)
Definition lowrankfunction.h:84
GridBuilder builder
Definition lowrankfunction.h:83
dftgrid(const std::size_t nradial, const std::size_t angular_order, const coord_3d origin=coord_3d())
Definition lowrankfunction.h:105
Definition lowrankfunction.h:53
bool do_print
Definition lowrankfunction.h:77
double get_radius() const
Definition lowrankfunction.h:56
double volume_element
Definition lowrankfunction.h:75
double radius
Definition lowrankfunction.h:76
double get_volume_element() const
Definition lowrankfunction.h:55
virtual ~gridbase()=default
void visualize(const std::string filename, const std::vector< Vector< double, NDIM > > &grid) const
Definition lowrankfunction.h:61
given a molecule, return a suitable grid
Definition lowrankfunction.h:269
std::vector< Vector< double, NDIM > > get_grid() const
Definition lowrankfunction.h:300
molecular_grid(const std::vector< Vector< double, NDIM > > origins, std::shared_ptr< gridbase > grid)
ctor takes centers of the grid and the grid builder
Definition lowrankfunction.h:292
std::vector< Vector< double, NDIM > > centers
Definition lowrankfunction.h:326
std::shared_ptr< gridbase > grid_builder
Definition lowrankfunction.h:327
molecular_grid(const Molecule &molecule, std::shared_ptr< gridbase > grid)
ctor takes molecule and grid builder
Definition lowrankfunction.h:298
molecular_grid(const std::vector< Vector< double, NDIM > > origins, const LowRankFunctionParameters &params)
ctor takes centers of the grid and the grid parameters
Definition lowrankfunction.h:273
grid with random points around the origin, with a Gaussian distribution
Definition lowrankfunction.h:121
static Vector< double, NDIM > gaussian_random_distribution(const Vector< double, NDIM > origin, double variance)
Definition lowrankfunction.h:182
Vector< double, NDIM > get_origin() const
Definition lowrankfunction.h:169
std::vector< Vector< double, NDIM > > get_grid() const
Definition lowrankfunction.h:129
randomgrid(const double volume_element, const double radius, const Vector< double, NDIM > origin=Vector< double, NDIM >(0.0))
Definition lowrankfunction.h:123
Vector< double, NDIM > origin
Definition lowrankfunction.h:194
double volume() const
Definition lowrankfunction.h:175
Definition y.cc:25
double(* f1)(const coord_3d &)
Definition derivatives.cc:55
char * p(char *buf, const char *name, int k, int initial_level, double thresh, int order)
Definition derivatives.cc:72
double(* f2)(const coord_3d &)
Definition derivatives.cc:56
static double lo
Definition dirac-hatom.cc:23
auto T(World &world, response_space &f) -> response_space
Definition global_functions.cc:34
Tensor< typename Tensor< T >::scalar_type > arg(const Tensor< T > &t)
Return a new tensor holding the argument of each element of t (complex types only)
Definition tensor.h:2502
static const double v
Definition hatom_sf_dirac.cc:20
#define MADNESS_CHECK(condition)
Check a condition — even in a release build the condition is always evaluated so it can have side eff...
Definition madness_exception.h:182
#define MADNESS_EXCEPTION(msg, value)
Macro for throwing a MADNESS exception.
Definition madness_exception.h:119
#define MADNESS_CHECK_THROW(condition, msg)
Check a condition — even in a release build the condition is always evaluated so it can have side eff...
Definition madness_exception.h:207
Main include file for MADNESS and defines Function interface.
constexpr double pi
Mathematical constant .
Definition constants.h:48
Namespace for all elements and tools of MADNESS.
Definition DFParameters.h:10
std::ostream & operator<<(std::ostream &os, const particle< PDIM > &p)
Definition lowrankfunction.h:397
static const char * filename
Definition legendre.cc:96
std::vector< Function< T, NDIM > > orthonormalize_rrcd(const std::vector< Function< T, NDIM > > &v, Tensor< T > &ovlp, const double tol, Tensor< integer > &piv, int &rank)
Definition vmra.h:594
std::vector< double > grid({0.00000000000000e+00, 7.74564534039568e-10, 2.47256327461087e-08, 1.86997862087340e-07, 7.83531011839395e-07, 2.37369448442132e-06, 5.85385578447464e-06, 1.25192693640128e-05, 2.41120079015646e-05, 4.28530598081574e-05, 7.14571735608571e-05, 1.13129528705579e-04, 1.71543841230856e-04, 2.50802052874775e-04, 3.55376294767219e-04, 4.90034340610786e-04, 6.59750258771429e-04, 8.69602422403842e-04, 1.12466142949409e-03, 1.42987080942386e-03, 1.78992364014940e-03, 2.20913836324745e-03, 2.69133715810767e-03, 3.23973021906902e-03, 3.85680917042725e-03, 4.54425265659531e-03, 5.30284686331574e-03, 6.13242336803442e-03, 7.03181629276860e-03, 7.99884025229928e-03, 9.03029006708050e-03, 1.01219626578657e-02, 1.12687009724886e-02, 1.24644592296948e-02, 1.37023882155871e-02, 1.49749388498821e-02, 1.62739817657269e-02, 1.75909402310371e-02, 1.89169333923840e-02, 2.02429265537309e-02, 2.15598850190411e-02, 2.28589279348859e-02, 2.41314785691809e-02, 2.53694075550732e-02, 2.65651658122794e-02, 2.77119041269022e-02, 2.88035767176875e-02, 2.98350265324687e-02, 3.08020504919994e-02, 3.17014434167336e-02, 3.25310199214522e-02, 3.32896141281727e-02, 3.39770576143407e-02, 3.45941365656990e-02, 3.51425296266603e-02, 3.56247284215205e-02, 3.60439431446186e-02, 3.64039959753441e-02, 3.67092053552739e-02, 3.69642643623641e-02, 3.71741165259965e-02, 3.73438324441572e-02, 3.74784904900007e-02, 3.75830647318932e-02, 3.76623229435371e-02, 3.77207372560624e-02, 3.77624096112071e-02, 3.77910137249598e-02, 3.78097547768664e-02, 3.78213475154040e-02, 3.78280129289835e-02, 3.78314930902835e-02, 3.78330832537561e-02, 3.78336797869059e-02, 3.78338420591352e-02, 3.78338660102034e-02, 3.78338692887670e-02, 3.78339467173722e-02, 3.78344713082705e-02, 3.78363997704733e-02, 3.78415404238543e-02, 3.78527910300104e-02, 3.78743388654655e-02, 3.79118156730933e-02, 3.79724014263881e-02, 3.80648723179728e-02, 3.81995899976738e-02, 3.83884307919462e-02, 3.86446553881524e-02, 3.89827212156132e-02, 3.94180414519212e-02, 3.99666961805422e-02, 4.06451026796855e-02, 4.14696530921280e-02, 4.24563287754274e-02, 4.36203014320401e-02, 4.49755316463183e-02, 4.65343756946706e-02, 4.83072114386611e-02, 5.03020937588608e-02, 5.25244493483274e-02, 5.49768197749344e-02, 5.76586605651183e-02, 6.05662026884007e-02, 6.36923812686800e-02, 6.70268346561297e-02, 7.05559752077735e-02, 7.42631312932333e-02, 7.81287582138778e-02, 8.21307139477169e-02, 8.62445939566306e-02, 9.04441177619647e-02, 9.47015586503668e-02, 9.89882067501835e-02, 1.03274854850000e-01, 1.07532295738402e-01, 1.11731819543736e-01, 1.15845699552650e-01, 1.19847655286489e-01, 1.23713282207134e-01, 1.27420438292594e-01, 1.30949578844237e-01, 1.34284032231687e-01, 1.37410210811966e-01, 1.40317752935249e-01, 1.42999593725433e-01, 1.45451964152040e-01, 1.47674319741506e-01, 1.49669202061706e-01, 1.51442037805696e-01, 1.53000881854049e-01, 1.54356112068327e-01, 1.55520084724940e-01, 1.56506760408239e-01, 1.57331310820682e-01, 1.58009717319825e-01, 1.58558372048446e-01, 1.58993692284754e-01, 1.59331758112215e-01, 1.59587982708421e-01, 1.59776823502693e-01, 1.59911541182394e-01, 1.60004012073979e-01, 1.60064597827274e-01, 1.60102074634902e-01, 1.60123622470357e-01, 1.60134873076513e-01, 1.60140013729894e-01, 1.60141942192096e-01, 1.60142466782995e-01, 1.60142544211600e-01, 1.60142551561710e-01, 1.60142701413062e-01, 1.60143716679379e-01, 1.60147448926215e-01, 1.60157397883103e-01, 1.60179171729209e-01, 1.60220874304751e-01, 1.60293404996033e-01, 1.60410659555492e-01, 1.60589622976194e-01, 1.60850348663762e-01, 1.61215821451274e-01, 1.61711705392906e-01, 1.62365980655838e-01, 1.63208477113433e-01, 1.64270315334478e-01, 1.65583268477170e-01, 1.67179061053879e-01, 1.69088622564332e-01, 1.71341315543301e-01, 1.73964158589697e-01, 1.76981065407121e-01, 1.80412120776561e-01, 1.84272913700736e-01, 1.88573946723051e-01, 1.93320138663612e-01, 1.98510435776210e-01, 2.04137543672572e-01, 2.10187789353851e-01, 2.16641119414426e-01, 2.23471237026998e-01, 2.30645876773227e-01, 2.38127212845837e-01, 2.45872392711143e-01, 2.53834185077787e-01, 2.61961728055316e-01, 2.70201360784824e-01, 2.78497519653254e-01, 2.86793678521684e-01, 2.95033311251191e-01, 3.03160854228721e-01, 3.11122646595364e-01, 3.18867826460670e-01, 3.26349162533280e-01, 3.33523802279509e-01, 3.40353919892081e-01, 3.46807249952656e-01, 3.52857495633936e-01, 3.58484603530297e-01, 3.63674900642895e-01, 3.68421092583456e-01, 3.72722125605771e-01, 3.76582918529946e-01, 3.80013973899386e-01, 3.83030880716811e-01, 3.85653723763206e-01, 3.87906416742175e-01, 3.89815978252628e-01, 3.91411770829337e-01, 3.92724723972029e-01, 3.93786562193074e-01, 3.94629058650669e-01, 3.95283333913601e-01, 3.95779217855233e-01, 3.96144690642745e-01, 3.96405416330313e-01, 3.96584379751016e-01, 3.96701634310474e-01, 3.96774165001756e-01, 3.96815867577298e-01, 3.96837641423404e-01, 3.96847590380292e-01, 3.96851322627128e-01, 3.96852337893446e-01, 3.96852487744797e-01, 3.96852501080501e-01, 3.96852763595511e-01, 3.96854542175721e-01, 3.96861080460548e-01, 3.96878509402610e-01, 3.96916653612757e-01, 3.96989709691015e-01, 3.97116771576032e-01, 3.97322182348883e-01, 3.97635696933973e-01, 3.98092445611658e-01, 3.98732694042811e-01, 3.99601401444597e-01, 4.00747584484612e-01, 4.02223500212562e-01, 4.04083666765013e-01, 4.06383745508240e-01, 4.09179312589000e-01, 4.12524551422180e-01, 4.16470900356858e-01, 4.21065691550599e-01, 4.26350817893211e-01, 4.32361464629548e-01, 4.39124941137697e-01, 4.46659646152579e-01, 4.54974196640938e-01, 4.64066746612100e-01, 4.73924517493193e-01, 4.84523556430928e-01, 4.95828733144949e-01, 5.07793979903267e-01, 5.20362772980513e-01, 5.33468847761119e-01, 5.47037133628609e-01, 5.60984889100643e-01, 5.75223012480244e-01, 5.89657498736383e-01, 6.04191009524582e-01, 6.18724520312781e-01, 6.33159006568921e-01, 6.47397129948521e-01, 6.61344885420555e-01, 6.74913171288045e-01, 6.88019246068651e-01, 7.00588039145897e-01, 7.12553285904216e-01, 7.23858462618236e-01, 7.34457501555971e-01, 7.44315272437064e-01, 7.53407822408226e-01, 7.61722372896585e-01, 7.69257077911467e-01, 7.76020554419616e-01, 7.82031201155953e-01, 7.87316327498565e-01, 7.91911118692306e-01, 7.95857467626984e-01, 7.99202706460164e-01, 8.01998273540924e-01, 8.04298352284151e-01, 8.06158518836602e-01, 8.07634434564552e-01, 8.08780617604567e-01, 8.09649325006353e-01, 8.10289573437506e-01, 8.10746322115191e-01, 8.11059836700281e-01, 8.11265247473132e-01, 8.11392309358149e-01, 8.11465365436407e-01, 8.11503509646554e-01, 8.11520938588616e-01, 8.11527476873443e-01, 8.11529255453653e-01, 8.11529517968663e-01, 8.11529541176345e-01, 8.11529996288764e-01, 8.11533079746512e-01, 8.11544414925490e-01, 8.11574630828117e-01, 8.11640760013437e-01, 8.11767414594025e-01, 8.11987697014585e-01, 8.12343809961018e-01, 8.12887338428862e-01, 8.13679186466493e-01, 8.14789161138720e-01, 8.16295206552722e-01, 8.18282301065180e-01, 8.20841040761609e-01, 8.24065941688939e-01, 8.28053501868532e-01, 8.32900071579870e-01, 8.38699586575474e-01, 8.45541223590310e-01, 8.53507040609252e-01, 8.62669665762836e-01, 8.73090098389325e-01, 8.84815683732401e-01, 8.97878318988171e-01, 9.12292943068421e-01, 9.28056355648415e-01, 9.45146402996031e-01, 9.63521558948635e-01, 9.83120919457849e-01, 1.00386461862596e+00, 1.02565466339202e+00, 1.04837617327936e+00, 1.07189900117812e+00, 1.09607970128622e+00, 1.12076380133619e+00, 1.14578832833438e+00, 1.17098453044675e+00, 1.19618073255911e+00, 1.22120525955731e+00, 1.24588935960728e+00, 1.27007005971537e+00, 1.29359288761414e+00, 1.31631439750148e+00, 1.33810444226754e+00, 1.35884814143565e+00, 1.37844750194486e+00, 1.39682265789747e+00, 1.41391270524508e+00, 1.42967611782508e+00, 1.44409074190533e+00, 1.45715337716110e+00, 1.46887896250417e+00, 1.47929939513066e+00, 1.48846202028425e+00, 1.49642783730319e+00, 1.50326947431802e+00, 1.50906898931363e+00, 1.51391555902497e+00, 1.51790311920456e+00, 1.52112802013189e+00, 1.52368675982832e+00, 1.52567385434078e+00, 1.52717989975478e+00, 1.52828987442701e+00, 1.52908172246464e+00, 1.52962525093248e+00, 1.52998136387891e+00, 1.53020164629947e+00, 1.53032830088006e+00, 1.53039443006538e+00, 1.53042464596801e+00, 1.53043598114699e+00, 1.53043906460474e+00, 1.53043951971715e+00, 1.53043956057345e+00, 1.53044036881866e+00, 1.53044584480621e+00, 1.53046597522602e+00, 1.53051963638178e+00, 1.53063707680903e+00, 1.53086200575658e+00, 1.53125321066044e+00, 1.53188564029607e+00, 1.53285090571065e+00, 1.53425716788373e+00, 1.53622839887814e+00, 1.53890302152796e+00, 1.54243195096192e+00, 1.54697607896990e+00, 1.55270325889653e+00, 1.55978486392298e+00, 1.56839200485204e+00, 1.57869150446938e+00, 1.59084173390576e+00, 1.60498842193063e+00, 1.62126055060600e+00, 1.63976645013905e+00, 1.66059020209859e+00, 1.68378845349035e+00, 1.70938773469079e+00, 1.73738236216545e+00, 1.76773299256309e+00, 1.80036587856233e+00, 1.83517285918360e+00, 1.87201209863831e+00, 1.91070956866812e+00, 1.95106125024279e+00, 1.99283601194713e+00, 2.03577910489527e+00, 2.07961619803354e+00, 2.12405786366203e+00, 2.16880441129773e+00, 2.21355095893343e+00, 2.25799262456193e+00, 2.30182971770019e+00, 2.34477281064834e+00, 2.38654757235268e+00, 2.42689925392735e+00, 2.46559672395716e+00, 2.50243596341186e+00, 2.53724294403314e+00, 2.56987583003238e+00, 2.60022646043001e+00, 2.62822108790467e+00, 2.65382036910512e+00, 2.67701862049687e+00, 2.69784237245642e+00, 2.71634827198947e+00, 2.73262040066483e+00, 2.74676708868971e+00, 2.75891731812608e+00, 2.76921681774343e+00, 2.77782395867248e+00, 2.78490556369893e+00, 2.79063274362556e+00, 2.79517687163354e+00, 2.79870580106750e+00, 2.80138042371732e+00, 2.80335165471173e+00, 2.80475791688482e+00, 2.80572318229939e+00, 2.80635561193503e+00, 2.80674681683888e+00, 2.80697174578644e+00, 2.80708918621368e+00, 2.80714284736945e+00, 2.80716297778926e+00, 2.80716845377681e+00, 2.80716926202201e+00, 2.80716933652781e+00, 2.80717083214942e+00, 2.80718096521951e+00, 2.80721821566121e+00, 2.80731751323041e+00, 2.80753483148951e+00, 2.80795105244894e+00, 2.80867495963022e+00, 2.80984524239094e+00, 2.81163142287735e+00, 2.81423364814557e+00, 2.81788132295423e+00, 2.82283059256804e+00, 2.82936071868432e+00, 2.83776942436573e+00, 2.84836731472069e+00, 2.86147150815776e+00, 2.87739863756603e+00, 2.89645740105067e+00, 2.91894085730732e+00, 2.94511867090741e+00, 2.97522951738903e+00, 3.00947385695645e+00, 3.04800727879293e+00, 3.09093460564947e+00, 3.13830493080142e+00, 3.19010773712238e+00, 3.24627022150018e+00, 3.30665591781433e+00, 3.37106467900872e+00, 3.43923404429894e+00, 3.51084198217487e+00, 3.58551096454407e+00, 3.66281329305810e+00, 3.74227756629515e+00, 3.82339614690749e+00, 3.90563346187864e+00, 3.98843494737088e+00, 4.07123643286313e+00, 4.15347374783428e+00, 4.23459232844661e+00, 4.31405660168366e+00, 4.39135893019769e+00, 4.46602791256689e+00, 4.53763585044283e+00, 4.60580521573305e+00, 4.67021397692744e+00, 4.73059967324159e+00, 4.78676215761938e+00, 4.83856496394034e+00, 4.88593528909229e+00, 4.92886261594884e+00, 4.96739603778531e+00, 5.00164037735273e+00, 5.03175122383435e+00, 5.05792903743445e+00, 5.08041249369110e+00, 5.09947125717573e+00, 5.11539838658400e+00, 5.12850258002107e+00, 5.13910047037603e+00, 5.14750917605745e+00, 5.15403930217373e+00, 5.15898857178754e+00, 5.16263624659620e+00, 5.16523847186442e+00, 5.16702465235083e+00, 5.16819493511155e+00, 5.16891884229283e+00, 5.16933506325225e+00, 5.16955238151135e+00, 5.16965167908056e+00, 5.16968892952225e+00, 5.16969906259234e+00, 5.16970055821395e+00, 5.16970070074519e+00, 5.16970361247163e+00, 5.16972333987315e+00, 5.16979586028538e+00, 5.16998917613145e+00, 5.17041225862441e+00, 5.17122257158139e+00, 5.17263189841936e+00, 5.17491024424232e+00, 5.17838764046389e+00, 5.18345374010175e+00, 5.19055515605250e+00, 5.20019056052898e+00, 5.21290362959218e+00, 5.22927398050981e+00, 5.24990630974935e+00, 5.27541799408998e+00, 5.30642546508522e+00, 5.34352970658452e+00, 5.38730125510908e+00, 5.43826510271255e+00, 5.49688591095665e+00, 5.56355394250644e+00, 5.63857210361471e+00, 5.72214446673726e+00, 5.81436660831299e+00, 5.91521805324665e+00, 6.02455706599185e+00, 6.14211796971732e+00, 6.26751111140535e+00, 6.40022552357701e+00, 6.53963426446189e+00, 6.68500234967725e+00, 6.83549712169932e+00, 6.99020084039182e+00, 7.14812522029963e+00, 7.30822758986854e+00, 7.46942830557509e+00, 7.63062902128164e+00, 7.79073139085055e+00, 7.94865577075836e+00, 8.10335948945085e+00, 8.25385426147293e+00, 8.39922234668829e+00, 8.53863108757317e+00, 8.67134549974483e+00, 8.79673864143286e+00, 8.91429954515833e+00, 9.02363855790353e+00, 9.12449000283719e+00, 9.21671214441292e+00, 9.30028450753547e+00, 9.37530266864374e+00, 9.44197070019353e+00, 9.50059150843763e+00, 9.55155535604110e+00, 9.59532690456566e+00, 9.63243114606496e+00, 9.66343861706020e+00, 9.68895030140082e+00, 9.70958263064037e+00, 9.72595298155800e+00, 9.73866605062120e+00, 9.74830145509768e+00, 9.75540287104843e+00, 9.76046897068629e+00, 9.76394636690786e+00, 9.76622471273082e+00, 9.76763403956879e+00, 9.76844435252577e+00, 9.76886743501873e+00, 9.76906075086480e+00, 9.76913327127703e+00, 9.76915299867855e+00, 9.76915591040499e+00, 9.76915619860098e+00, 9.76916219846504e+00, 9.76920284848150e+00, 9.76935228306099e+00, 9.76975062707579e+00, 9.77062242513200e+00, 9.77229214503906e+00, 9.77519618483832e+00, 9.77989091339106e+00, 9.78705638895934e+00, 9.79749552527532e+00, 9.81212860682633e+00, 9.83198319082184e+00, 9.85817956879146e+00, 9.89191209222870e+00, 9.93442679048652e+00, 9.98699582179647e+00, 1.00508893966709e+01, 1.01273458942923e+01, 1.02175409544910e+01, 1.03225563687840e+01, 1.04433496124935e+01, 1.05807248555829e+01, 1.07353062625751e+01, 1.09075143424092e+01, 1.10975460386012e+01, 1.13053591604470e+01, 1.15306616495989e+01, 1.17729060559758e+01, 1.20312894658455e+01, 1.23047589865396e+01, 1.25920227503353e+01, 1.28915662583656e+01, 1.32016737478116e+01, 1.35204541357752e+01, 1.38458709746296e+01, 1.41757757494884e+01, 1.45079437615262e+01, 1.48401117735640e+01, 1.51700165484228e+01, 1.54954333872772e+01, 1.58142137752408e+01, 1.61243212646868e+01, 1.64238647727170e+01, 1.67111285365127e+01, 1.69845980572069e+01, 1.72429814670766e+01, 1.74852258734535e+01, 1.77105283626053e+01, 1.79183414844512e+01, 1.81083731806431e+01, 1.82805812604773e+01, 1.84351626674695e+01, 1.85725379105589e+01, 1.86933311542684e+01, 1.87983465685614e+01, 1.88885416287600e+01, 1.89649981263815e+01, 1.90288917012559e+01, 1.90814607325659e+01, 1.91239754308237e+01, 1.91577079542609e+01, 1.91839043322306e+01, 1.92037589162261e+01, 1.92183919977771e+01, 1.92288311340931e+01, 1.92359966096613e+01, 1.92406913382141e+01, 1.92435953780133e+01, 1.92452650979204e+01, 1.92461368959766e+01, 1.92465352399914e+01, 1.92466846745709e+01, 1.92467253245873e+01, 1.92467313244514e+01, 1.92467319433609e+01, 1.92467450814034e+01, 1.92468340936939e+01, 1.92471613140832e+01, 1.92480335772770e+01, 1.92499425738442e+01, 1.92535987985263e+01, 1.92599578423078e+01, 1.92702379990799e+01, 1.92859284083041e+01, 1.93087872289342e+01, 1.93408296294082e+01, 1.93843056757497e+01, 1.94416684964889e+01, 1.95155333909890e+01, 1.96086288188291e+01, 1.97237404546042e+01, 1.98636497079429e+01, 2.00310682866674e+01, 2.02285705167748e+01, 2.04585252224186e+01, 2.07230290096769e+01, 2.10238427883022e+01, 2.13623333059298e+01, 2.17394213608060e+01, 2.21555382047478e+01, 2.26105914517851e+01, 2.31039415749319e+01, 2.36343898099582e+01, 2.42001779979100e+01, 2.47990005951181e+01, 2.54280287686536e+01, 2.60839461849717e+01, 2.67629957981511e+01, 2.74610366597973e+01, 2.81736095129742e+01, 2.88960097044501e+01, 2.96233657592419e+01, 3.03507218140337e+01, 3.10731220055095e+01, 3.17856948586865e+01, 3.24837357203327e+01, 3.31627853335120e+01, 3.38187027498302e+01, 3.44477309233657e+01, 3.50465535205737e+01, 3.56123417085256e+01, 3.61427899435518e+01, 3.66361400666986e+01, 3.70911933137359e+01, 3.75073101576777e+01, 3.78843982125540e+01, 3.82228887301816e+01, 3.85237025088068e+01, 3.87882062960652e+01, 3.90181610017090e+01, 3.92156632318163e+01, 3.93830818105409e+01, 3.95229910638796e+01, 3.96381026996546e+01, 3.97311981274948e+01, 3.98050630219949e+01, 3.98624258427341e+01, 3.99059018890755e+01, 3.99379442895496e+01, 3.99608031101797e+01, 3.99764935194039e+01, 3.99867736761760e+01, 3.99931327199574e+01, 3.99967889446395e+01, 3.99986979412068e+01, 3.99995702044006e+01, 3.99998974247898e+01, 3.99999864370804e+01, 3.99999995751228e+01})
const std::vector< Function< T, NDIM > > & change_tree_state(const std::vector< Function< T, NDIM > > &v, const TreeState finalstate, const bool fence=true)
change tree state of the functions
Definition vmra.h:277
std::vector< double > norm2s(World &world, const std::vector< Function< T, NDIM > > &v)
Computes the 2-norms of a vector of functions.
Definition vmra.h:827
void truncate(World &world, response_space &v, double tol, bool fence)
Definition basic_operators.cc:30
std::vector< Function< TENSOR_RESULT_TYPE(T, R), NDIM > > transform(World &world, const std::vector< Function< T, NDIM > > &v, const Tensor< R > &c, bool fence=true)
Transforms a vector of functions according to new[i] = sum[j] old[j]*c[j,i].
Definition vmra.h:689
@ reconstructed
s coeffs at the leaves only
Definition funcdefaults.h:60
std::vector< Function< T, NDIM > > orthonormalize_canonical(const std::vector< Function< T, NDIM > > &v, const Tensor< T > &ovlp, double lindep)
Definition vmra.h:501
static const Slice _(0,-1, 1)
@ OT_SLATER
1/r
Definition operatorinfo.h:15
@ OT_GAUSS
exp(-r)
Definition operatorinfo.h:16
void svd(const Tensor< T > &a, Tensor< T > &U, Tensor< typename Tensor< T >::scalar_type > &s, Tensor< T > &VT)
Compute the singluar value decomposition of an n-by-m matrix using *gesvd.
Definition lapack.cc:739
void print(const T &t, const Ts &... ts)
Print items to std::cout (items separated by spaces) and terminate with a new line.
Definition print.h:225
response_space apply(World &world, std::vector< std::vector< std::shared_ptr< real_convolution_3d > > > &op, response_space &f)
Definition basic_operators.cc:39
std::vector< Function< T, NDIM > > append(const std::vector< Function< T, NDIM > > &lhs, const std::vector< Function< T, NDIM > > &rhs)
combine two vectors
Definition vmra.h:649
NDIM & f
Definition mra.h:2448
void error(const char *msg)
Definition world.cc:139
NDIM const Function< R, NDIM > & g
Definition mra.h:2448
double inner(response_space &a, response_space &b)
Definition response_functions.h:442
void normalize(World &world, std::vector< Function< T, NDIM > > &v, bool fence=true)
Normalizes a vector of functions — v[i] = v[i].scale(1.0/v[i].norm2())
Definition vmra.h:1680
Vector< double, 3 > coord_3d
Definition funcplot.h:1042
static XNonlinearSolver< std::vector< Function< T, NDIM > >, T, vector_function_allocator< T, NDIM > > nonlinear_vector_solver(World &world, const long nvec)
Definition nonlinsol.h:284
void matrix_inner(DistributedMatrix< T > &A, const std::vector< Function< T, NDIM > > &f, const std::vector< Function< T, NDIM > > &g, bool sym=false)
Definition distpm.cc:46
Function< T, KDIM+LDIM > hartree_product(const std::vector< Function< T, KDIM > > &left, const std::vector< Function< T, LDIM > > &right)
Performs a Hartree/outer product on the two given low-dimensional function vectors.
Definition mra.h:1864
double get_size(World &world, const std::vector< Function< T, NDIM > > &v)
Definition vmra.h:1721
Function< T, NDIM > copy(const Function< T, NDIM > &f, const std::shared_ptr< WorldDCPmapInterface< Key< NDIM > > > &pmap, bool fence=true)
Create a new copy of the function with different distribution and optional fence.
Definition mra.h:2034
void syev(const Tensor< T > &A, Tensor< T > &V, Tensor< typename Tensor< T >::scalar_type > &e)
Real-symmetric or complex-Hermitian eigenproblem.
Definition lapack.cc:969
Definition mraimpl.h:50
static const double b
Definition nonlinschro.cc:119
static const double d
Definition nonlinschro.cc:121
static const double a
Definition nonlinschro.cc:118
double Q(double a)
Definition relops.cc:20
static const double thresh
Definition rk.cc:45
the low-rank functor is what the LowRankFunction will represent
Definition lowrankfunction.h:410
virtual ~LRFunctorBase()
Definition lowrankfunction.h:412
virtual Tensor< T >::scalar_type norm2() const
Definition lowrankfunction.h:421
friend std::vector< Function< T, LDIM > > inner(const LRFunctorBase &functor, const std::vector< Function< T, LDIM > > &rhs, const particle< LDIM > p1, const particle< LDIM > p2)
Definition lowrankfunction.h:426
virtual Function< T, LDIM > inner(const Function< T, LDIM > &rhs, const particle< LDIM > p1, const particle< LDIM > p2) const
Definition lowrankfunction.h:416
virtual World & world() const =0
friend Function< T, LDIM > inner(const LRFunctorBase &functor, const Function< T, LDIM > &rhs, const particle< LDIM > p1, const particle< LDIM > p2)
Definition lowrankfunction.h:430
virtual std::vector< Function< T, LDIM > > inner(const std::vector< Function< T, LDIM > > &rhs, const particle< LDIM > p1, const particle< LDIM > p2) const =0
virtual T operator()(const Vector< T, NDIM > &r) const =0
Definition lowrankfunction.h:438
std::vector< Function< T, LDIM > > b
the lo-dim functions
Definition lowrankfunction.h:460
LRFunctorF12(const std::shared_ptr< SeparatedConvolution< T, LDIM > > f12, const Function< T, LDIM > &a, const Function< T, LDIM > &b)
delegate to the other ctor with vector arguments
Definition lowrankfunction.h:453
Tensor< T >::scalar_type norm2() const
Definition lowrankfunction.h:499
std::vector< Function< T, LDIM > > a
Definition lowrankfunction.h:460
T operator()(const Vector< double, NDIM > &r) const
Definition lowrankfunction.h:528
std::vector< Function< T, LDIM > > inner(const std::vector< Function< T, LDIM > > &rhs, const particle< LDIM > p1, const particle< LDIM > p2) const
Definition lowrankfunction.h:464
World & world() const
Definition lowrankfunction.h:463
std::shared_ptr< SeparatedConvolution< T, LDIM > > f12
a two-particle function
Definition lowrankfunction.h:459
LRFunctorF12(const std::shared_ptr< SeparatedConvolution< T, LDIM > > f12, const std::vector< Function< T, LDIM > > &a, const std::vector< Function< T, LDIM > > &b)
Definition lowrankfunction.h:440
Definition lowrankfunction.h:562
World & world() const
Definition lowrankfunction.h:565
std::vector< Function< T, LDIM > > inner(const std::vector< Function< T, LDIM > > &rhs, const particle< LDIM > p1, const particle< LDIM > p2) const
Definition lowrankfunction.h:569
LRFunctorPure(const Function< T, NDIM > &f)
Definition lowrankfunction.h:564
Function< T, NDIM > f
a hi-dim function
Definition lowrankfunction.h:567
Tensor< T >::scalar_type norm2() const
Definition lowrankfunction.h:581
T operator()(const Vector< double, NDIM > &r) const
Definition lowrankfunction.h:577
Definition lowrankfunction.h:20
std::string rhsfunctiontype() const
Definition lowrankfunction.h:48
double tol() const
Definition lowrankfunction.h:44
int optimize() const
Definition lowrankfunction.h:45
double gamma() const
Definition lowrankfunction.h:42
double volume_element() const
Definition lowrankfunction.h:43
std::string gridtype() const
Definition lowrankfunction.h:46
LowRankFunctionParameters()
Definition lowrankfunction.h:22
std::string f12type() const
Definition lowrankfunction.h:49
void read_and_set_derived_values(World &world, const commandlineparser &parser, std::string tag)
Definition lowrankfunction.h:37
std::string orthomethod() const
Definition lowrankfunction.h:47
double radius() const
Definition lowrankfunction.h:41
Definition lowrankfunction.h:199
long total_n
Definition lowrankfunction.h:204
long index
Definition lowrankfunction.h:202
void operator++()
Definition lowrankfunction.h:249
bool operator()() const
Definition lowrankfunction.h:253
Vector< double, NDIM > hivec
Definition lowrankfunction.h:200
long n_per_dim
Definition lowrankfunction.h:203
cartesian_grid(const double volume_per_gridpoint, const double radius)
Definition lowrankfunction.h:207
cartesian_grid & operator=(const cartesian_grid< NDIM > &other)
Definition lowrankfunction.h:227
double volume_per_gridpoint() const
Definition lowrankfunction.h:243
void initialize(const double lo, const double hi)
Definition lowrankfunction.h:233
Vector< double, NDIM > lovec
Definition lowrankfunction.h:200
Vector< double, NDIM > increment
Definition lowrankfunction.h:205
Vector< double, NDIM > get_coordinates() const
Definition lowrankfunction.h:257
std::vector< long > stride
Definition lowrankfunction.h:201
cartesian_grid(const long n_per_dim, const double lo, const double hi)
Definition lowrankfunction.h:217
cartesian_grid(const cartesian_grid< NDIM > &other)
Definition lowrankfunction.h:222
very simple command line parser
Definition commandlineparser.h:15
Definition lowrankfunction.h:332
std::array< int, PDIM > dims
Definition lowrankfunction.h:333
particle(const int p1, const int p2, const int p3)
Definition lowrankfunction.h:360
particle complement() const
return the other particle
Definition lowrankfunction.h:352
static particle particle1()
convenience for particle 1 (the left/first particle)
Definition lowrankfunction.h:339
bool is_first() const
assuming two particles only
Definition lowrankfunction.h:379
std::array< int, PDIM > get_array() const
type conversion to std::array
Definition lowrankfunction.h:373
bool is_last() const
assuming two particles only
Definition lowrankfunction.h:381
particle(const std::vector< int > p)
Definition lowrankfunction.h:361
static particle particle2()
convenience for particle 2 (the right/second particle)
Definition lowrankfunction.h:345
particle(const int p)
Definition lowrankfunction.h:358
std::string str() const
Definition lowrankfunction.h:365
particle()=default
default constructor
particle(const int p1, const int p2)
Definition lowrankfunction.h:359
std::enable_if_t< DUMMYDIM==2, std::tuple< int, int > > get_tuple() const
Definition lowrankfunction.h:389
std::enable_if_t< DUMMYDIM==3, std::tuple< int, int, int > > get_tuple() const
Definition lowrankfunction.h:393
std::enable_if_t< DUMMYDIM==1, std::tuple< int > > get_tuple() const
Definition lowrankfunction.h:385
Definition timing_utilities.h:9
double tag(const std::string msg)
Definition timing_utilities.h:45
bool do_print
Definition timing_utilities.h:12
InputParameters param
Definition tdse.cc:203
static const double omega
Definition tdse_example.cc:51
void split(const Range< ConcurrentHashMap< int, int >::iterator > &range)
Definition test_hashthreaded.cc:63
double aa
Definition testbsh.cc:68
constexpr coord_t one(1.0)
constexpr std::size_t NDIM
Definition testgconv.cc:54
double h(const coord_1d &r)
Definition testgconv.cc:175
static Molecule molecule
Definition testperiodicdft.cc:39
Defines operations on vectors of Functions.