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 std::string get_tag() const override {
37 return std::string("lrf");
38 }
39
40
41 void read_and_set_derived_values(World& world, const commandlineparser& parser, std::string tag) {
42 read_input_and_commandline_options(world,parser,tag);
43 }
44
45 double radius() const {return get<double>("radius");}
46 double gamma() const {return get<double>("gamma");}
47 double volume_element() const {return get<double>("volume_element");}
48 double tol() const {return get<double>("tol");}
49 int optimize() const {return get<int>("optimize");}
50 std::string gridtype() const {return get<std::string>("gridtype");}
51 std::string orthomethod() const {return get<std::string>("orthomethod");}
52 std::string rhsfunctiontype() const {return get<std::string>("rhsfunctiontype");}
53 std::string f12type() const {return get<std::string>("f12type");}
54 };
55
56
57 class gridbase {
58 public:
59 double get_volume_element() const {return volume_element;}
60 double get_radius() const {return radius;}
61
62 virtual ~gridbase() = default;
63 // visualize the grid in xyz format
64 template<std::size_t NDIM>
65 void visualize(const std::string filename, const std::vector<Vector<double,NDIM>>& grid) const {
66 print("visualizing grid to file",filename);
67 print("a total of",grid.size(),"grid points");
68 std::ofstream file(filename);
69 for (const auto& r : grid) {
70 // formatted output
71 file << std::fixed << std::setprecision(6);
72 for (int i=0; i<NDIM; ++i) file << r[i] << " ";
73 file << std::endl;
74 }
75 file.close();
76 }
77
78 protected:
79 double volume_element=0.1;
80 double radius=3;
81 bool do_print=false;
82 };
83
84 template<std::size_t NDIM>
85 class dftgrid : public gridbase {
86 public:
88 explicit dftgrid(const double volume_element, const double radius) {
89 // increase number of radial grid points until the volume element is below the threshold
90 double current_ve=1.0;
91 std::size_t nradial=10;
92 while (current_ve>volume_element) {
93 nradial+=10;
94 print("trying nradial",nradial);
95 GridBuilder tmp;
96 tmp.set_angular_order(7);
97 tmp.set_nradial(nradial);
98 tmp.make_grid();
99 double volume=4./3. *M_PI * std::pow(radius,3.0);
100 auto npoints=tmp.get_points().size();
101 current_ve=volume/npoints;
102 print("volume, npoints, volume element",volume,npoints,current_ve);
103 }
104 builder.set_nradial(nradial);
106 }
107
108
109 explicit dftgrid(const std::size_t nradial, const std::size_t angular_order, const coord_3d origin=coord_3d()) {
110 static_assert(NDIM==3,"DFT Grids only defined for NDIM=3");
111 builder.set_nradial(nradial);
112 builder.set_angular_order(angular_order);
113 builder.set_origin(origin);
115 }
116
117 std::vector<Vector<double,NDIM>> get_grid() const {
118 return builder.get_points();
119 }
120
121 };
122
123 /// grid with random points around the origin, with a Gaussian distribution
124 template<std::size_t NDIM>
125 class randomgrid : public gridbase {
126 public:
128 : gridbase(), origin(origin) {
129 this->volume_element=volume_element;
130 this->radius=radius;
131 }
132
133 std::vector<Vector<double,NDIM>> get_grid() const {
134 std::vector<Vector<double, NDIM>> grid;
135 long npoint_within_volume=volume()/volume_element;
136 if (do_print) print("npoint_within_volume",npoint_within_volume);
137
139 auto is_in_cell = [&cell](const Vector<double, NDIM>& r) {
140 for (int d = 0; d < NDIM; ++d) if (r[d] < cell(d, 0) or r[d] > cell(d, 1)) return false;
141 return true;
142 };
143 double rad=radius;
144 auto o=origin;
145 auto is_in_sphere = [&rad,&o](const Vector<double, NDIM>& r) {
146 return ((r-o).normf()<rad);
147 };
148
149 // set variance such that about 70% of all grid points sits within the radius
150 double variance=radius;
151 if (NDIM==2) variance=0.6*radius;
152 if (NDIM==3) variance=0.5*radius;
153 long maxrank=10*npoint_within_volume;
154 long rank=0;
155 for (int r=0; r<maxrank; ++r) {
156 auto tmp = gaussian_random_distribution(origin, variance);
157 if (not is_in_cell(tmp)) continue;
158 if (is_in_sphere(tmp)) ++rank;
159 grid.push_back(tmp);
160 if (rank==npoint_within_volume) break;
161 }
162 if (do_print) {
163 print("origin ",origin);
164 print("radius ",radius);
165 print("grid points in volume ",rank);
166 print("total grid points ",grid.size());
167 print("ratio ",rank/double(grid.size()));
168 print("volume element ",volume()/rank);
169 }
170 return grid;
171 }
172
174 return origin;
175 }
176
177 private:
178
179 double volume() const {
180 MADNESS_CHECK(NDIM>0 and NDIM<4);
181 if (NDIM==1) return 2.0*radius;
182 if (NDIM==2) return constants::pi*radius*radius;
183 if (NDIM==3) return 4.0 / 3.0 * constants::pi * std::pow(radius, 3.0);
184 }
185
187 std::random_device rd{};
188 std::mt19937 gen{rd()};
189 Vector<double,NDIM> result;
190 for (int i = 0; i < NDIM; ++i) {
191 std::normal_distribution<> d{origin[i], variance};
192 result[i]=d(gen);
193 }
194
195 return result;
196 }
197
199
200 };
201
202 template<std::size_t NDIM>
205 std::vector<long> stride;
206 long index=0;
210
211 cartesian_grid(const double volume_per_gridpoint, const double radius) {
212 double length_per_gridpoint=std::pow(volume_per_gridpoint,1.0/NDIM);
213 n_per_dim=ceil(2.0*radius/length_per_gridpoint);
214 print("length per gridpoint, n_per_dim",length_per_gridpoint,n_per_dim);
215 print("volume_per_gridpoint",std::pow(length_per_gridpoint,NDIM));
216 initialize(-radius,radius);
217 print("increment",increment);
218 }
219
220
221 cartesian_grid(const long n_per_dim, const double lo, const double hi)
223 initialize(lo,hi);
224 }
225
227 hivec(other.hivec), stride(other.stride), index(0), n_per_dim(other.n_per_dim),
228 total_n(other.total_n), increment(other.increment) {
229 }
230
232 cartesian_grid<NDIM> tmp(other);
233 std::swap(*this,other);
234 return *this;
235 }
236
237 void initialize(const double lo, const double hi) {
238 lovec.fill(lo);
239 hivec.fill(hi);
240 increment=(hivec-lovec)*(1.0/double(n_per_dim-1));
241 stride=std::vector<long>(NDIM,1l);
242 total_n=std::pow(n_per_dim,NDIM);
243 for (long i=NDIM-2; i>=0; --i) stride[i]=n_per_dim*stride[i+1];
244
245 }
246
247 double volume_per_gridpoint() const{
248 double volume=1.0;
249 for (int i=0; i<NDIM; ++i) volume*=(hivec[i]-lovec[i]);
250 return volume/total_n;
251 }
252
253 void operator++() {
254 index++;
255 }
256
257 bool operator()() const {
258 return index < total_n;
259 }
260
263 for (int idim=0; idim<NDIM; ++idim) {
264 tmp[idim]=(index/stride[idim])%n_per_dim;
265 }
266 return lovec+tmp*increment;
267 }
268
269 };
270
271 /// given a molecule, return a suitable grid
272 template<std::size_t NDIM>
273 class molecular_grid : public gridbase {
274
275 public:
276 /// ctor takes centers of the grid and the grid parameters
277 molecular_grid(const std::vector<Vector<double,NDIM>> origins, const LowRankFunctionParameters& params)
278 : centers(origins)
279 {
280 if (centers.size()==0) centers.push_back(Vector<double,NDIM>(0) );
281 if (params.gridtype()=="random") grid_builder=std::make_shared<randomgrid<NDIM>>(params.volume_element(),params.radius());
282 // else if (params.gridtype()=="cartesian") grid_builder=std::make_shared<cartesian_grid<NDIM>>(params.volume_element(),params.radius());
283 else if (params.gridtype()=="dftgrid") {
284 if constexpr (NDIM==3) {
285 grid_builder=std::make_shared<dftgrid<NDIM>>(params.volume_element(),params.radius());
286 } else {
287 MADNESS_EXCEPTION("no dft grid with NDIM != 3",1);
288 }
289 } else {
290 MADNESS_EXCEPTION("no such grid type",1);
291 }
292 }
293
294
295 /// ctor takes centers of the grid and the grid builder
296 molecular_grid(const std::vector<Vector<double,NDIM>> origins, std::shared_ptr<gridbase> grid)
297 : centers(origins), grid_builder(grid) {
298 if (centers.size()==0) centers.push_back({0,0,0});
299 }
300
301 /// ctor takes molecule and grid builder
302 molecular_grid(const Molecule& molecule, std::shared_ptr<gridbase> grid) : molecular_grid(molecule.get_all_coords_vec(),grid) {}
303
304 std::vector<Vector<double,NDIM>> get_grid() const {
305 MADNESS_CHECK_THROW(grid_builder,"no grid builder given in molecular_grid");
306 MADNESS_CHECK_THROW(centers.size()>0,"no centers given in molecular_grid");
307 std::vector<Vector<double,NDIM>> grid;
308 for (const auto& coords : centers) {
309 print("atom sites",coords);
310 if (auto builder=dynamic_cast<dftgrid<NDIM>*>(grid_builder.get())) {
311 if constexpr (NDIM==3) {
312 dftgrid<NDIM> b1(builder->builder.get_nradial(),builder->builder.get_angular_order(),coords);
313 auto atomgrid=b1.get_grid();
314 grid.insert(grid.end(),atomgrid.begin(),atomgrid.end());
315 } else {
316 MADNESS_EXCEPTION("no DFT grid for NDIM /= 3",1);
317 }
318 } else if (auto builder=dynamic_cast<randomgrid<NDIM>*>(grid_builder.get())) {
319 randomgrid<NDIM> b1(builder->get_volume_element(),builder->get_radius(),coords);
320 auto atomgrid=b1.get_grid();
321 grid.insert(grid.end(),atomgrid.begin(),atomgrid.end());
322 } else {
323 MADNESS_EXCEPTION("no such grid builder",1);
324 }
325 }
326 return grid;
327 }
328
329 private:
330 std::vector<Vector<double,NDIM>> centers;
331 std::shared_ptr<gridbase> grid_builder;
332
333 };
334
335template<std::size_t PDIM>
336struct particle {
337 std::array<int,PDIM> dims;
338
339 /// default constructor
340 particle() = default;
341
342 /// convenience for particle 1 (the left/first particle)
344 particle p;
345 for (int i=0; i<PDIM; ++i) p.dims[i]=i;
346 return p;
347 }
348 /// convenience for particle 2 (the right/second particle)
350 particle p;
351 for (int i=0; i<PDIM; ++i) p.dims[i]=i+PDIM;
352 return p;
353 }
354
355 /// return the other particle
358 if (is_first()) return particle2();
359 return particle1();
360 }
361
362 particle(const int p) : particle(std::vector<int>(1,p)) {}
363 particle(const int p1, const int p2) : particle(std::vector<int>({p1,p2})) {}
364 particle(const int p1, const int p2,const int p3) : particle(std::vector<int>({p1,p2,p3})) {}
365 particle(const std::vector<int> p) {
366 for (int i=0; i<PDIM; ++i) dims[i]=p[i];
367 }
368
369 std::string str() const {
370 std::stringstream ss;
371 ss << *this;
372 return ss.str();
373 }
374
375
376 /// type conversion to std::array
377 std::array<int,PDIM> get_array() const {
378 return dims;
379 }
380
381
382 /// assuming two particles only
383 bool is_first() const {return dims[0]==0;}
384 /// assuming two particles only
385 bool is_last() const {return dims[0]==(PDIM);}
386
387 template<std::size_t DUMMYDIM=PDIM>
388 typename std::enable_if_t<DUMMYDIM==1, std::tuple<int>>
389 get_tuple() const {return std::tuple<int>(dims[0]);}
390
391 template<std::size_t DUMMYDIM=PDIM>
392 typename std::enable_if_t<DUMMYDIM==2, std::tuple<int,int>>
393 get_tuple() const {return std::tuple<int,int>(dims[0],dims[1]);}
394
395 template<std::size_t DUMMYDIM=PDIM>
396 typename std::enable_if_t<DUMMYDIM==3, std::tuple<int,int,int>>
397 get_tuple() const {return std::tuple<int,int,int>(dims[0],dims[1],dims[2]);}
398};
399
400template<std::size_t PDIM>
401std::ostream& operator<<(std::ostream& os, const particle<PDIM>& p) {
402 os << "(";
403 for (auto i=0; i<PDIM-1; ++i) os << p.dims[i] << ";";
404 os << p.dims[PDIM-1] << ")";
405 return os;
406}
407
408/// the low-rank functor is what the LowRankFunction will represent
409
410/// derive from this class :
411/// must implement in inner product
412/// may implement an operator()(const coord_nd&)
413template<typename T, std::size_t NDIM, std::size_t LDIM=NDIM/2>
415
416 virtual ~LRFunctorBase() {};
417 virtual std::vector<Function<T,LDIM>> inner(const std::vector<Function<T,LDIM>>& rhs,
418 const particle<LDIM> p1, const particle<LDIM> p2) const =0;
419
420 virtual Function<T,LDIM> inner(const Function<T,LDIM>& rhs, const particle<LDIM> p1, const particle<LDIM> p2) const {
421 return inner(std::vector<Function<T,LDIM>>({rhs}),p1,p2)[0];
422 }
423
424 virtual T operator()(const Vector<T,NDIM>& r) const =0;
425 virtual typename Tensor<T>::scalar_type norm2() const {
426 MADNESS_EXCEPTION("L2 norm not implemented",1);
427 }
428
429 virtual World& world() const =0;
430 friend std::vector<Function<T,LDIM>> inner(const LRFunctorBase& functor, const std::vector<Function<T,LDIM>>& rhs,
431 const particle<LDIM> p1, const particle<LDIM> p2) {
432 return functor.inner(rhs,p1,p2);
433 }
434 friend Function<T,LDIM> inner(const LRFunctorBase& functor, const Function<T,LDIM>& rhs,
435 const particle<LDIM> p1, const particle<LDIM> p2) {
436 return functor.inner(rhs,p1,p2);
437 }
438
439};
440
441template<typename T, std::size_t NDIM, std::size_t LDIM=NDIM/2>
442struct LRFunctorF12 : public LRFunctorBase<T,NDIM> {
443// LRFunctorF12() = default;
444 LRFunctorF12(const std::shared_ptr<SeparatedConvolution<T,LDIM>> f12, const std::vector<Function<T,LDIM>>& a,
445 const std::vector<Function<T,LDIM>>& b) : f12(f12), a(a), b(b) {
446
447 // if a or b are missing, they are assumed to be 1
448 // you may not provide a or b, but if you do they have to have the same size because they are summed up
449 if (a.size()>0 and b.size()>0)
450 MADNESS_CHECK_THROW(a.size()==b.size(), "a and b must have the same size");
451 if (a.size()==0) this->a.resize(b.size());
452 if (b.size()==0) this->b.resize(a.size());
453 MADNESS_CHECK_THROW(this->a.size()==this->b.size(), "a and b must have the same size");
454 }
455
456 /// delegate to the other ctor with vector arguments
458 const Function<T,LDIM>& a, const Function<T,LDIM>& b)
459 : LRFunctorF12(f12,std::vector<Function<T,LDIM>>({a}), std::vector<Function<T,LDIM>>({b})) {}
460
461
462private:
463 std::shared_ptr<SeparatedConvolution<T,LDIM>> f12; ///< a two-particle function
464 std::vector<Function<T,LDIM>> a,b; ///< the lo-dim functions
465public:
466
467 World& world() const {return f12->get_world();}
468 std::vector<Function<T,LDIM>> inner(const std::vector<Function<T,LDIM>>& rhs,
469 const particle<LDIM> p1, const particle<LDIM> p2) const {
470
471 std::vector<Function<T,LDIM>> result;
472 // functor is now \sum_i a_i(1) b_i(2) f12
473 // result(1) = \sum_i \int a_i(1) f(1,2) b_i(2) rhs(2) d2
474 // = \sum_i a_i(1) \int f(1,2) b_i(2) rhs(2) d2
475 World& world=rhs.front().world();
476
477 const int nbatch=30;
478 for (int i=0; i<rhs.size(); i+=nbatch) {
479 std::vector<Function<T,LDIM>> rhs_batch;
480 auto begin= rhs.begin()+i;
481 auto end= (i+nbatch)<rhs.size() ? rhs.begin()+i+nbatch : rhs.end();
482 std::copy(begin,end, std::back_inserter(rhs_batch));
483 auto tmp2= zero_functions_compressed<T,LDIM>(world,rhs_batch.size());
484
485 if (a.size()==0) tmp2=apply(world,*(f12),rhs_batch);
486
487 for (int ia=0; ia<a.size(); ia++) {
488 auto premultiply= p1.is_first() ? a[ia] : b[ia];
489 auto postmultiply= p1.is_first() ? b[ia] : a[ia];
490
491 auto tmp=copy(world,rhs_batch);
492 if (premultiply.is_initialized()) tmp=rhs_batch*premultiply;
493 auto tmp1=apply(world,*(f12),tmp);
494 if (postmultiply.is_initialized()) tmp1=tmp1*postmultiply;
495 tmp2+=tmp1;
496 }
497
498 for (auto& t : tmp2) result.push_back(t);
499 }
500 return result;
501 }
502
503 typename Tensor<T>::scalar_type norm2() const {
505 [](const Vector<double, LDIM>& r) { return 1.0; });
506 std::vector<Function<T, LDIM>> pre, post;
507 std::size_t sz = a.size();
508 if (sz == 0) {
509 pre = std::vector<Function<T, LDIM>>(1, one);
510 post = std::vector<Function<T, LDIM>>(1, one);
511 } else {
512 pre = (a.front().is_initialized()) ? a : std::vector<Function<T, LDIM>>(sz, one);
513 post = (b.front().is_initialized()) ? b : std::vector<Function<T, LDIM>>(sz, one);
514 }
515
516 const SeparatedConvolution<T,LDIM>& f12a=*(f12);
518
519 // \int f(1,2)^2 d1d2 = \int f(1,2)^2 pre(1)^2 post(2)^2 d1 d2
520 // || \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
521 std::vector<Function<T,LDIM>> aa,bb;
522 for (std::size_t i=0; i<pre.size(); ++i) {
523 aa=append(aa,(pre[i]*pre));
524 bb=append(bb,(post[i]*post));
525 }
526// typename Tensor<T>::scalar_type term1 =madness::inner(mul(world(),post,post),f12sq(mul(world(),pre,pre)));
527 typename Tensor<T>::scalar_type term1 =madness::inner(bb,f12sq(aa));
528 return sqrt(term1);
529
530 }
531
533
534 if (a.size()==0) return 0.0;
535 auto split = [](const Vector<double,NDIM>& r) {
536 Vector<double,LDIM> first, second;
537 for (int i=0; i<LDIM; ++i) {
538 first[i]=r[i];
539 second[i]=r[i+LDIM];
540 }
541 return std::make_pair(first,second);
542 };
543
544 double gamma=f12->info.mu;
545 auto [first,second]=split(r);
546
547
548 double result=0.0;
549 for (std::size_t ia=0; ia<a.size(); ++ia) {
550 double result1=1.0;
551 if (a[ia].is_initialized()) result1*=a[ia](first);
552 if (b[ia].is_initialized()) result1*=b[ia](first);
553 if (f12->info.type==OT_SLATER) result1*=exp(-gamma*(first-second).normf());
554 else if (f12->info.type==OT_GAUSS) result1*=exp(-gamma* madness::inner(first-second,first-second));
555 else {
556 MADNESS_EXCEPTION("no such operator_type",1);
557 }
558 result+=result1;
559 }
560 return result;
561
562 }
563};
564
565template<typename T, std::size_t NDIM, std::size_t LDIM=NDIM/2>
566struct LRFunctorPure : public LRFunctorBase<T,NDIM> {
567 LRFunctorPure() = default;
569 World& world() const {return f.world();}
570
571 Function<T, NDIM> f; ///< a hi-dim function
572
573 std::vector<Function<T,LDIM>> inner(const std::vector<Function<T,LDIM>>& rhs,
574 const particle<LDIM> p1, const particle<LDIM> p2) const {
575 return madness::innerXX<LDIM>(f,rhs,p1.get_array(),p2.get_array());
576// std::vector<Function<T,LDIM>> result;
577// for (const auto& r : rhs) result.push_back(madness::inner(f,r,p1.get_tuple(),p2.get_tuple()));
578// return result;
579 }
580
582 return f(r);
583 }
584
585 typename Tensor<T>::scalar_type norm2() const {
586 return f.norm2();
587 }
588};
589
590
591 /// LowRankFunction represents a hi-dimensional (NDIM) function as a sum of products of low-dimensional (LDIM) functions
592
593 /// f(1,2) = \sum_i g_i(1) h_i(2)
594 /// a LowRankFunction can be created from a hi-dim function directly, or from a composite like f(1,2) phi(1) psi(2),
595 /// where f(1,2) is a two-particle function (e.g. a Slater function)
596 template<typename T, std::size_t NDIM, std::size_t LDIM=NDIM/2>
598 public:
599
601 double rank_revealing_tol=1.e-8; // rrcd tol
602 std::string orthomethod="canonical";
603 bool do_print=false;
604 std::vector<Function<T,LDIM>> g,h;
607
609
611 double tol, std::string orthomethod) : world(g.front().world()),
613
614 }
615
616 /// shallow copy ctor
619 g(other.g), h(other.h) {
620 }
621
622 /// deep copy
623 friend LowRankFunction copy(const LowRankFunction& other) {
625 }
626
627 LowRankFunction& operator=(const LowRankFunction& f) { // Assignment required for storage in vector
628 LowRankFunction ff(f);
629 std::swap(ff.g,g);
630 std::swap(ff.h,h);
631 return *this;
632 }
633
634 /// function evaluation
636 Vector<double,LDIM> first, second;
637 for (int i=0; i<LDIM; ++i) {
638 first[i]=r[i];
639 second[i]=r[i+LDIM];
640 }
641 double result=0.0;
642 for (int i=0; i<rank(); ++i) result+=g[i](first)*h[i](second);
643 return result;
644 }
645
646 /*
647 * arithmetic section
648 */
649
650 /// addition
652 LowRankFunction<T,NDIM> result=copy(*this);
653 result+=b;
654 return result;
655 }
656 /// subtraction
658 LowRankFunction<T,NDIM> result=copy(*this);
659 result-=b;
660 return result;
661 }
662
663 /// in-place addition
665
666 g=append(g,copy(b.g));
667 h=append(h,copy(b.h));
668 return *this;
669 }
670
671 /// in-place subtraction
673 g=append(g,-1.0*b.g); // operator* implies deep copy of b.g
674 h=append(h,copy(b.h));
675 return *this;
676 }
677
678 /// scale by a scalar
679 template<typename Q>
683
684 /// out-of-place scale by a scalar (no type conversion)
688
689 /// multiplication with a scalar
690 friend LowRankFunction operator*(const T a, const LowRankFunction& other) {
691 return other*a;
692 }
693
694 /// in-place scale by a scalar (no type conversion)
696 g=g*a;
697 return *this;
698 }
699
700 /// l2 norm
702 auto tmp1=matrix_inner(world,h,h);
703 auto tmp2=matrix_inner(world,g,g);
704 return sqrt(tmp1.trace(tmp2));
705 }
706
707 std::vector<Function<T,LDIM>> get_functions(const particle<LDIM>& p) const {
708 MADNESS_CHECK(p.is_first() or p.is_last());
709 if (p.is_first()) return g;
710 return h;
711 }
712
713 std::vector<Function<T,LDIM>> get_g() const {return g;}
714 std::vector<Function<T,LDIM>> get_h() const {return h;}
715
716 long rank() const {return g.size();}
717
718 /// return the size in GByte
719 double size() const {
720 double sz=get_size(world,g);
721 sz+=get_size(world,h);
722 return sz;
723 }
724
726 auto fapprox=hartree_product(g[0],h[0]);
727 for (int i=1; i<g.size(); ++i) fapprox+=hartree_product(g[i],h[i]);
728 return fapprox;
729 }
730
731 /// orthonormalize the argument vector
732 std::vector<Function<T,LDIM>> orthonormalize(const std::vector<Function<T,LDIM>>& g) const {
733
734 double tol=rank_revealing_tol;
735 std::vector<Function<T,LDIM>> g2;
736 auto ovlp=matrix_inner(world,g,g);
737 if (orthomethod=="canonical") {
738 tol*=0.01;
739 print("orthonormalizing with method/tol",orthomethod,tol);
740 g2=orthonormalize_canonical(g,ovlp,tol);
741 } else if (orthomethod=="cholesky") {
742 print("orthonormalizing with method/tol",orthomethod,tol);
743 g2=orthonormalize_rrcd(g,ovlp,tol);
744 }
745 else {
746 MADNESS_EXCEPTION("no such orthomethod",1);
747 }
748 double tight_thresh=FunctionDefaults<3>::get_thresh()*0.1;
749 return truncate(g2,tight_thresh);
750 }
751
752
753 /// optimize the lrf using the lrfunctor
754
755 /// @param[in] nopt number of iterations (wrt to Alg. 4.3 in Halko)
756 void optimize(const LRFunctorBase<T,NDIM>& lrfunctor1, const long nopt=1) {
757 timer t(world);
759 for (int i=0; i<nopt; ++i) {
760 // orthonormalize h
762 t.tag("ortho1");
763 g=truncate(inner(lrfunctor1,h,p2,p1));
764 t.tag("inner1");
766 t.tag("ortho2");
767 h=truncate(inner(lrfunctor1,g,p1,p1));
768 t.tag("inner2");
769 }
770 }
771
772 /// remove linear dependencies using rank-revealing cholesky decomposition
773 ///
774 /// @return this with g orthonormal and h orthogonal
775 void reorthonormalize_rrcd(double thresh=-1.0) {
776
777 // with g~ being the orthonormalized g and h~ the orthonormalized h (from rrcd)
778 // with
779 // ovlp = <g_i | g_j> = L L^T
780 // we get
781 // | g~_i(1) > = | g_j(1) > piv (L^T)^(-1) // with piv permuting/truncating, and LT^(-1) already rank-reduced
782 // | g_i(1) > = | g~_j(1) > L^T piv^T
783 // = | g_j(1) > piv LT^(-1) LT piv^T
784 // thus:
785 // f(1,2) = sum_i | g_i(1)>< h_i(2) |
786 // = | g_i(1)> piv_g LTg^(-1) LTg piv_g^(-1) <h_i(2)| (1)
787 // = | g_i(1)> piv_g LTg^(-1) ( LTg piv_g^(-1) <h_i(2)| ) (2)
788 // = | g_i(1)> piv_g LTg^(-1) <h'_i(2)| (3)
789 // = | g_i(1)> piv_g LTg^(-1) piv_h Lh Lh^(-1) piv_h^T <h'_i(2)| (4)
790 // = (| g_i(1)> piv_g LTg^(-1) piv_h Lh ) ( Lh^(-1) piv_h^T LTg piv_g^(-1) <h_i(2)| )
791 // Notes:
792 // piv^(-1) = piv^T since piv is a unitary permutation matrix
793 // Thus
794 // | h~_i(2) > = |h_i(2) > piv_g Lg piv_h LTh^(-1)
795
796 double tol=rank_revealing_tol;
797
798 // no pivot_inverse is necessary..
799 auto pivot_vec = [&] (const std::vector<Function<T,LDIM>>& v, const Tensor<integer>& ipiv) {
800 std::size_t rank=ipiv.size();
801 std::vector<Function<T,LDIM> > pv(rank);
802 for(int i=0;i<ipiv.size();++i) pv[i]=v[ipiv[i]]; // all elements in [rank,v.size()) = 0 due to lindep
803 return pv;
804 };
805
806 auto pivot_mat = [&] (const Tensor<T>& t, const Tensor<integer>& ipiv) {
807 // std::size_t rank=ipiv.size();
808 Tensor<T> pt(t.dim(0),t.dim(1));
809 // for(int i=0;i<ipiv.size();++i) pt(i,_)=t(ipiv[i],_); // all elements in [rank,v.size()) = 0 due to lindep
810 for(int i=0;i<ipiv.size();++i) pt(ipiv[i],_)=t(i,_); // all elements in [rank,v.size()) = 0 due to lindep
811 return pt;
812 };
813
814 auto pivot = [](const Tensor<integer>& ipiv) {
815 Tensor<T> piv(ipiv.size(),ipiv.size());
816 for (int i=0; i<ipiv.size(); ++i) piv(ipiv[i],i)=1.0;
817 return piv;
818 };
819
820
821 auto LT = [&](const Tensor<T>& ovlp1) {
822 Tensor<integer> piv;
823 auto ovlp=copy(ovlp1); // copy since ovlp will be destroyed
824 // auto ovlp=matrix_inner(world,v,v);
825 int rank;
826 rr_cholesky(ovlp,tol,piv,rank); // destroys ovlp and gives back Upper ∆ Matrix from CD
827 print("initial, final rank",ovlp.dim(0),rank);
828 // piv=piv(Slice(0,rank-1));
829 Tensor<T> p=pivot(piv);
830 Tensor<T> L=inner(p,transpose(ovlp));
831
832 double err1=(ovlp1-inner(L,L,1,1)).normf();
833 L=L(_,Slice(0,rank-1)); // (new, old)
834 double err2=(ovlp1-inner(L,L,1,1)).normf();
835 print("err1, err2", err1,err2);
836 return std::make_pair(L,piv);
837 };
838
839 // orthonormalize g: g_ortho= transform(g_piv,LT_g^(-1))
840 auto ovlp_g=matrix_inner(world,g,g);
841 // ovlp_g=0.5*(ovlp_g+transpose(ovlp_g));
842 // for (int i=0; i<ovlp_g.dim(0); ++i) ovlp_g(i,i)+=1.e-12;
843 auto [L_g,piv_g]=LT(ovlp_g); // (1)
844 L_g=pivot_mat(L_g,piv_g);
845 double zero=(ovlp_g-inner(L_g,L_g,1,1)).normf();
846 print("zero",zero);
847 // auto hpiv_g=pivot_vec(h,piv_g); // (2) : | h_i > piv_g
848 // auto ovlp_hsmall=matrix_inner(world,hpiv_g,hpiv_g); // (3) : < h_i | h_j >
849 // auto ovlp1=inner(LT_g,inner(ovlp_hsmall,LT_g,1,1)); // (3) : LTg piv_g <h_i | h_i> piv_g Lg
850 // auto [LTh,piv_h]=LT(ovlp1); // (3)
851
852 // auto gpiv_g=pivot_vec(g,piv_g); // | g_i > piv_g
853 // auto gtrans=inner(pivot_mat(inverse(LT_g),piv_h),LTh,1,1); // LTg^(-1) piv_g Lh
854 // auto htrans=inner(pivot_mat(transpose(LT_g),piv_h),inverse(LTh)); // LTg^(-1) piv_g Lh
855 auto gpiv_g=pivot_vec(g,piv_g); // | g_i > piv_g
856 auto gtrans=inverse(L_g);
857 auto htrans=L_g;
858 auto hpiv_g=pivot_vec(h,piv_g); // | h_i > piv_g
859
860 g=transform(world,gpiv_g,(gtrans));
861 h=transform(world,hpiv_g,transpose(htrans));
862 auto Gortho=matrix_inner(world,g,g);
863 for (int i=0; i<Gortho.dim(0); ++i) Gortho(i,i)-=1.0;
864 print("G_ortho.normf()",Gortho.normf(),Gortho.dim(0));
865 print("end of line");
866 // g=gpiv_g;
867 // h=hpiv_g;
868 }
869
870 void reorthonormalize(double thresh=-1.0) {
871 if (orthomethod=="canonical") {
873 } else if (orthomethod=="cholesky") {
875 } else {
876 MADNESS_EXCEPTION("no such orthomethod",1);
877 }
878 }
879 /// after external operations g might not be orthonormal and/or optimal -- reorthonormalize
880
881 /// orthonormalization similar to Bischoff, Harrison, Valeev, JCP 137 104103 (2012), Sec II C 3
882 /// f =\sum_i g_i h_i
883 /// = g X- (X+)^T (Y+)^T Y- h
884 /// = g X- U S V^T Y- h
885 /// = g (X- U) (S V^T Y-) h
886 /// requires 2 matrix_inner and 2 transforms. g and h are optimal
887 /// @param[in] thresh SVD threshold
890 Tensor<T> ovlp_g = matrix_inner(world, g, g);
891 Tensor<T> ovlp_h = matrix_inner(world, h, h);
892
893 ovlp_g=0.5*(ovlp_g+transpose(ovlp_g));
894 ovlp_h=0.5*(ovlp_h+transpose(ovlp_h));
895
896 auto [eval_g, evec_g] = syev(ovlp_g);
897 auto [eval_h, evec_h] = syev(ovlp_h);
898
899 // get relevant part of the eigenvalues and eigenvectors
900 // eigenvalues are sorted in ascending order
901 auto get_slice = [](auto eval, double thresh) {
902 // remove small/negative eigenvalues
903 eval.screen(thresh);
904 Slice s;
905 for (int i=0; i<eval.size(); ++i) {
906 MADNESS_CHECK_THROW(eval[i]>=0.0,"negative eigenvalues in reorthonormalize");
907 if (eval[i]>thresh) {
908 return s=Slice(i,-1); // from i to the end
909 break;
910 }
911 }
912 return s;
913 };
914
915 Slice gslice=get_slice(eval_g,1.e-12);
916 Slice hslice=get_slice(eval_h,1.e-12);
917
918 Tensor<T> Xplus=copy(evec_g(_,gslice));
919 Tensor<T> Xminus=copy(evec_g(_,gslice));
920 Tensor<T> Yplus=copy(evec_h(_,hslice));
921 Tensor<T> Yminus=copy(evec_h(_,hslice));
922 eval_g=copy(eval_g(gslice));
923 eval_h=copy(eval_h(hslice));
924
925 for (int i=0; i<eval_g.size(); ++i) Xplus(_,i)*=std::pow(eval_g(i),0.5);
926 for (int i=0; i<eval_g.size(); ++i) Xminus(_,i)*=std::pow(eval_g(i),-0.5);
927 for (int i=0; i<eval_h.size(); ++i) Yplus(_,i)*=std::pow(eval_h(i),0.5);
928 for (int i=0; i<eval_h.size(); ++i) Yminus(_,i)*=std::pow(eval_h(i),-0.5);
929
930 Tensor<T> M=madness::inner(Xplus,Yplus,0,0); // (X+)^T Y+
931 auto [U,s,VT]=svd(M);
932
933 // truncate
934 typename Tensor<T>::scalar_type s_accumulated=0.0;
935 int i=s.size()-1;
936 for (;i>=0; i--) {
937 s_accumulated+=s[i];
938 if (s_accumulated>thresh) {
939 i++;
940 break;
941 }
942 }
943 for (int j=0; j<s.size(); ++j) VT(j,_)*=s[j];
944 Tensor<T> XX=madness::inner(Xminus,U,1,0);
945 Tensor<T> YY=madness::inner(Yminus,VT,1,1);
946
949 }
950
951
952 double check_orthonormality(const std::vector<Function<T,LDIM>>& v) const {
954 return check_orthonormality(ovlp);
955 }
956
957 double check_orthonormality(const Tensor<T>& ovlp) const {
958 timer t(world);
960 Tensor<T> ovlp2=ovlp;
961 for (int i=0; i<ovlp2.dim(0); ++i) ovlp2(i,i)-=1.0;
962 if (world.rank()==0 and do_print) {
963 print("absmax",ovlp2.absmax());
964 print("l2",ovlp2.normf()/ovlp2.size());
965 }
966 t.tag("check_orthonoramality");
967 return ovlp.absmax();
968 }
969
970 /// compute the l2 error |functor - \sum_i g_ih_i|_2
971
972 /// \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
973 /// since we are subtracting large numbers the numerics are sensitive, and NaN may be returned..
974 double l2error(const LRFunctorBase<T,NDIM>& lrfunctor1) const {
975
976 timer t(world);
978
979 // \int f(1,2)^2 d1d2
980 double term1 = lrfunctor1.norm2();
981 term1=term1*term1;
982 t.tag("computing term1");
983
984 // \int f(1,2) pre(1) post(2) \sum_i g(1) h(2) d1d2
985// double term2=madness::inner(pre*g,f12(post*h));
986 double term2=madness::inner(g,inner(lrfunctor1,h,p2,p1));
987 t.tag("computing term2");
988
989 // g functions are orthonormal
990 // \int gh(1,2)^2 d1d2 = \int \sum_{ij} g_i(1) g_j(1) h_i(2) h_j(2) d1d2
991 // = \sum_{ij} \int g_i(1) g_j(1) d1 \int h_i(2) h_j(2) d2
992 // = \sum_{ij} delta_{ij} \int h_i(2) h_j(2) d2
993 // = \sum_{i} \int h_i(2) h_i(2) d2
994 double zero=check_orthonormality(g);
995 if (zero>1.e-10) print("g is not orthonormal",zero);
996 // double term3a=madness::inner(h,h);
997 auto tmp1=matrix_inner(world,h,h);
998 auto tmp2=matrix_inner(world,g,g);
999 double term3=tmp1.trace(tmp2);
1000// print("term3/a/diff",term3a,term3,term3-term3a);
1001 t.tag("computing term3");
1002
1003 double arg=term1-2.0*term2+term3;
1004 if (arg<0.0) {
1005 print("negative l2 error");
1006 arg*=-1.0;
1007// throw std::runtime_error("negative argument in l2error");
1008 }
1009 double error=sqrt(arg)/sqrt(term1);
1010 if (world.rank()==0 and do_print) {
1011 print("term1,2,3, error",term1, term2, term3, " --",error);
1012 }
1013
1014 return error;
1015 }
1016
1017 };
1018
1019 // This interface is necessary to compute inner products
1020 template<typename T, std::size_t NDIM>
1022 World& world=a.world;
1023 return (matrix_inner(world,a.g,b.g).emul(matrix_inner(world,a.h,b.h))).sum();
1024 }
1025
1026
1027
1028// template<typename T, std::size_t NDIM>
1029// LowRankFunction<T,NDIM> inner(const Function<T,NDIM>& lhs, const LowRankFunction<T,NDIM>& rhs,
1030// const std::tuple<int> v1, const std::tuple<int> v2) {
1031// World& world=rhs.world;
1032// // int lhs(1,2) rhs(2,3) d2 = \sum \int lhs(1,2) g_i(2) h_i(3) d2
1033// // = \sum \int lhs(1,2) g_i(2) d2 h_i(3)
1034// LowRankFunction<T, NDIM + NDIM - 2> result(world);
1035// result.h=rhs.h;
1036// decltype(rhs.g) g;
1037// for (int i=0; i<rhs.rank(); ++i) {
1038// g.push_back(inner(lhs,rhs.g[i],{v1},{0}));
1039// }
1040// result.g=g;
1041// return result;
1042// }
1043
1044 /**
1045 * inner product: LowRankFunction lrf; Function f, g; double d
1046 * lrf(1,3) = inner(lrf(1,2), lrf(2,3))
1047 * lrf(1,3) = inner(lrf(1,2), f(2,3))
1048 * g(1) = inner(lrf(1,2), f(2))
1049 * d = inner(lrf(1,2), f(1,2))
1050 * d = inner(lrf(1,2), lrf(1,2))
1051 */
1052
1053 /// lrf(1,3) = inner(full(1,2), lrf(2,3))
1054
1055 /// @param[in] f1 the first function
1056 /// @param[in] f2 the second function
1057 /// @param[in] p1 the integration variable of the first function
1058 /// @param[in] p2 the integration variable of the second function
1059 template<typename T, std::size_t NDIM, std::size_t PDIM>
1061 const particle<PDIM> p1, const particle<PDIM> p2) {
1062 auto result=inner(f2,f1,p2,p1);
1063 std::swap(result.g,result.h);
1064 return result;
1065 }
1066
1067 /// lrf(1,3) = inner(lrf(1,2), full(2,3))
1068
1069 /// @param[in] f1 the first function
1070 /// @param[in] f2 the second function
1071 /// @param[in] p1 the integration variable of the first function
1072 /// @param[in] p2 the integration variable of the second function
1073 template<typename T, std::size_t NDIM, std::size_t PDIM>
1075 const particle<PDIM> p1, const particle<PDIM> p2) {
1076 static_assert(TensorTypeData<T>::iscomplex==false, "complex inner in LowRankFunction not implemented");
1077 World& world=f1.world;
1078 static_assert(2*PDIM==NDIM);
1079 // int f(1,2) k(2,3) d2 = \sum \int g_i(1) h_i(2) k(2,3) d2
1080 // = \sum g_i(1) \int h_i(2) k(2,3) d2
1081 LowRankFunction<T, NDIM> result(world);
1082 if (p1.is_last()) { // integrate over 2: result(1,3) = lrf(1,2) f(2,3)
1083 result.g = f1.g;
1085 result.h=innerXX<PDIM>(f2,f1.h,p2.get_array(),particle<PDIM>::particle1().get_array());
1086 } else if (p1.is_first()) { // integrate over 1: result(2,3) = lrf(1,2) f(1,3)
1087 result.g = f1.h; // correct! second variable of f1 becomes first variable of result
1089 result.h=innerXX<PDIM>(f2,f1.g,p2.get_array(),particle<PDIM>::particle1().get_array());
1090 }
1091 return result;
1092 }
1093
1094 /// lrf(1,3) = inner(lrf(1,2), lrf(2,3))
1095
1096 /// @param[in] f1 the first function
1097 /// @param[in] f2 the second function
1098 /// @param[in] p1 the integration variable of the first function
1099 /// @param[in] p2 the integration variable of the second function
1100 template<typename T, std::size_t NDIM, std::size_t PDIM>
1102 const particle<PDIM> p1, const particle<PDIM> p2) {
1103 World& world=f1.world;
1104 static_assert(2*PDIM==NDIM);
1105
1106 // inner(lrf(1,2) ,lrf(2,3) ) = \sum_ij g1_i(1) <h1_i(2) g2_j(2)> h2_j(3)
1107 auto matrix=matrix_inner(world,f2.get_functions(p2),f1.get_functions(p1));
1108 auto htilde=transform(world,f2.get_functions(p2.complement()),matrix);
1109 auto gg=copy(world,f1.get_functions(p1.complement()));
1110 return LowRankFunction<T,NDIM>(gg,htilde,f1.rank_revealing_tol,f1.orthomethod);
1111 }
1112
1113 /// f(1) = inner(lrf(1,2), f(2))
1114
1115 /// @param[in] f1 the first function
1116 /// @param[in] vf vector of the second functions
1117 /// @param[in] p1 the integration variable of the first function
1118 /// @param[in] p2 the integration variable of the second function, dummy variable for consistent notation
1119 template<typename T, std::size_t NDIM, std::size_t PDIM>
1120 std::vector<Function<T,NDIM-PDIM>> inner(const LowRankFunction<T,NDIM>& f1, const std::vector<Function<T,PDIM>>& vf,
1122 World& world=f1.world;
1123 static_assert(2*PDIM==NDIM);
1124 MADNESS_CHECK(p2.is_first());
1125
1126 // inner(lrf(1,2), f_k(2) ) = \sum_i g1_i(1) <h1_i(2) f_k(2)>
1127 auto matrix=matrix_inner(world,f1.get_functions(p1),vf);
1128 return transform(world,f1.get_functions(p1.complement()),matrix);
1129 }
1130
1131 /// f(1) = inner(lrf(1,2), f(2))
1132
1133 /// @param[in] f1 the first function
1134 /// @param[in] vf the second function
1135 /// @param[in] p1 the integration variable of the first function
1136 /// @param[in] p2 the integration variable of the second function, dummy variable for consistent notation
1137 template<typename T, std::size_t NDIM, std::size_t PDIM>
1140 return inner(f1,std::vector<Function<T,PDIM>>({f2}),p1,p2)[0];
1141 }
1142
1143 template<typename T, std::size_t NDIM, std::size_t LDIM=NDIM/2>
1145 public:
1146
1149
1151 std::vector<Vector<double,LDIM>> origins; ///< origins of the molecular grid
1152
1156
1159
1161
1162 LowRankFunctionFactory& set_radius(const double radius) {
1163 parameters.set_user_defined_value("radius",radius);
1164 return *this;
1165 }
1166 LowRankFunctionFactory& set_volume_element(const double volume_element) {
1167 parameters.set_user_defined_value("volume_element",volume_element);
1168 return *this;
1169 }
1172 return *this;
1173 }
1174 LowRankFunctionFactory& set_orthomethod(const std::string orthomethod) {
1175 parameters.set_user_defined_value("orthomethod",orthomethod);
1176 return *this;
1177 }
1178
1180 World& world=lrfunctor.world();
1181 bool do_print=true;
1182 timer t1(world);
1183 t1.do_print=do_print;
1184 auto orthomethod=parameters.orthomethod();
1185 auto rank_revealing_tol=parameters.tol();
1186
1187 // get sampling grid
1189 auto grid=mgrid.get_grid();
1190 if (world.rank()==0) print("grid size",grid.size());
1191
1192 auto Y=Yformer(lrfunctor,grid,parameters.rhsfunctiontype());
1193 t1.tag("Yforming");
1194
1195 auto ovlp=matrix_inner(world,Y,Y); // error in symmetric matrix_inner, use non-symmetric form here!
1196 t1.tag("compute ovlp");
1197 auto g=truncate(orthonormalize_rrcd(Y,ovlp,rank_revealing_tol));
1198 t1.tag("rrcd/truncate/thresh");
1199 auto sz=get_size(world,g);
1200 if (world.rank()==0 and do_print) print("gsize",sz);
1201// check_orthonormality(g);
1202
1203 if (world.rank()==0 and do_print) {
1204 print("Y.size()",Y.size());
1205 print("g.size()",g.size());
1206 }
1207
1208 auto h=truncate(inner(lrfunctor,g,p1,p1));
1209 t1.tag("Y backprojection with truncation");
1211
1212 }
1213
1214 /// 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)
1215 std::vector<Function<T,LDIM>> Yformer(const LRFunctorBase<T,NDIM>& lrfunctor1, const std::vector<Vector<double,LDIM>>& grid,
1216 const std::string rhsfunctiontype, const double exponent=30.0) const {
1217
1218 World& world=lrfunctor1.world();
1219 std::vector<Function<double,LDIM>> Y;
1220 if (rhsfunctiontype=="exponential") {
1221 std::vector<Function<double,LDIM>> omega;
1222 double coeff=std::pow(2.0*exponent/constants::pi,0.25*LDIM);
1223 for (const auto& point : grid) {
1224 omega.push_back(FunctionFactory<double,LDIM>(world)
1225 .functor([&point,&exponent,&coeff](const Vector<double,LDIM>& r)
1226 {
1227 auto r_rel=r-point;
1228 return coeff*exp(-exponent*madness::inner(r_rel,r_rel));
1229 }));
1230 }
1231 Y=inner(lrfunctor1,omega,p2,p1);
1232 } else {
1233 MADNESS_EXCEPTION("confused rhsfunctiontype",1);
1234 }
1235 auto norms=norm2s(world,Y);
1236 std::vector<Function<double,LDIM>> Ynormalized;
1237
1238 for (int i=0; i<Y.size(); ++i) if (norms[i]>parameters.tol()) Ynormalized.push_back(Y[i]);
1239 normalize(world,Ynormalized);
1240 return Ynormalized;
1241 }
1242
1243 };
1244
1245
1246} // namespace madness
1247
1248#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
long dim(int i) const
Returns the size of dimension i.
Definition basetensor.h:147
long size() const
Returns the number of elements in the tensor.
Definition basetensor.h:138
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
FunctionFactory & f(T(*f)(const coordT &))
Definition function_factory.h:187
A multiresolution adaptive numerical function.
Definition mra.h:139
Definition lowrankfunction.h:1144
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:1215
LowRankFunctionFactory(const LowRankFunctionParameters param, const Molecule &molecule)
Definition lowrankfunction.h:1157
const particle< LDIM > p2
Definition lowrankfunction.h:1148
LowRankFunctionFactory & set_radius(const double radius)
Definition lowrankfunction.h:1162
const particle< LDIM > p1
Definition lowrankfunction.h:1147
LowRankFunctionFactory & set_rank_revealing_tol(const double rrtol)
Definition lowrankfunction.h:1170
LowRankFunctionFactory(const LowRankFunctionParameters param, const std::vector< Vector< double, LDIM > > origins={})
Definition lowrankfunction.h:1154
LowRankFunction< T, NDIM > project(const LRFunctorBase< T, NDIM > &lrfunctor) const
Definition lowrankfunction.h:1179
LowRankFunctionParameters parameters
Definition lowrankfunction.h:1150
LowRankFunctionFactory & set_volume_element(const double volume_element)
Definition lowrankfunction.h:1166
std::vector< Vector< double, LDIM > > origins
origins of the molecular grid
Definition lowrankfunction.h:1151
LowRankFunctionFactory & set_orthomethod(const std::string orthomethod)
Definition lowrankfunction.h:1174
LowRankFunctionFactory(const LowRankFunctionFactory &other)=default
LowRankFunction represents a hi-dimensional (NDIM) function as a sum of products of low-dimensional (...
Definition lowrankfunction.h:597
void optimize(const LRFunctorBase< T, NDIM > &lrfunctor1, const long nopt=1)
optimize the lrf using the lrfunctor
Definition lowrankfunction.h:756
friend LowRankFunction operator*(const T a, const LowRankFunction &other)
multiplication with a scalar
Definition lowrankfunction.h:690
double rank_revealing_tol
Definition lowrankfunction.h:601
std::vector< Function< T, LDIM > > get_h() const
Definition lowrankfunction.h:714
std::vector< Function< T, LDIM > > get_functions(const particle< LDIM > &p) const
Definition lowrankfunction.h:707
std::string orthomethod
Definition lowrankfunction.h:602
LowRankFunction operator-(const LowRankFunction &b) const
subtraction
Definition lowrankfunction.h:657
LowRankFunction operator*(const Q a) const
scale by a scalar
Definition lowrankfunction.h:680
std::vector< Function< T, LDIM > > orthonormalize(const std::vector< Function< T, LDIM > > &g) const
orthonormalize the argument vector
Definition lowrankfunction.h:732
double check_orthonormality(const std::vector< Function< T, LDIM > > &v) const
Definition lowrankfunction.h:952
LowRankFunction operator*(const T a) const
out-of-place scale by a scalar (no type conversion)
Definition lowrankfunction.h:685
void reorthonormalize_rrcd(double thresh=-1.0)
Definition lowrankfunction.h:775
TensorTypeData< T >::scalar_type norm2() const
l2 norm
Definition lowrankfunction.h:701
void reorthonormalize_canonical(double thresh=-1.0)
after external operations g might not be orthonormal and/or optimal – reorthonormalize
Definition lowrankfunction.h:888
long rank() const
Definition lowrankfunction.h:716
LowRankFunction(const LowRankFunction &other)
shallow copy ctor
Definition lowrankfunction.h:617
bool do_print
Definition lowrankfunction.h:603
void reorthonormalize(double thresh=-1.0)
Definition lowrankfunction.h:870
LowRankFunction & operator*=(const T a)
in-place scale by a scalar (no type conversion)
Definition lowrankfunction.h:695
LowRankFunction & operator=(const LowRankFunction &f)
Definition lowrankfunction.h:627
const particle< LDIM > p2
Definition lowrankfunction.h:606
Function< T, NDIM > reconstruct() const
Definition lowrankfunction.h:725
std::vector< Function< T, LDIM > > get_g() const
Definition lowrankfunction.h:713
const particle< LDIM > p1
Definition lowrankfunction.h:605
LowRankFunction & operator-=(const LowRankFunction &b)
in-place subtraction
Definition lowrankfunction.h:672
LowRankFunction(std::vector< Function< T, LDIM > > g, std::vector< Function< T, LDIM > > h, double tol, std::string orthomethod)
Definition lowrankfunction.h:610
std::vector< Function< T, LDIM > > h
Definition lowrankfunction.h:604
LowRankFunction & operator+=(const LowRankFunction &b)
in-place addition
Definition lowrankfunction.h:664
double l2error(const LRFunctorBase< T, NDIM > &lrfunctor1) const
compute the l2 error |functor - \sum_i g_ih_i|_2
Definition lowrankfunction.h:974
LowRankFunction(World &world)
Definition lowrankfunction.h:608
LowRankFunction operator+(const LowRankFunction &b) const
addition
Definition lowrankfunction.h:651
double check_orthonormality(const Tensor< T > &ovlp) const
Definition lowrankfunction.h:957
std::vector< Function< T, LDIM > > g
Definition lowrankfunction.h:604
World & world
Definition lowrankfunction.h:600
friend LowRankFunction copy(const LowRankFunction &other)
deep copy
Definition lowrankfunction.h:623
double size() const
return the size in GByte
Definition lowrankfunction.h:719
T operator()(const Vector< double, NDIM > &r) const
function evaluation
Definition lowrankfunction.h:635
Definition molecule.h:129
class for holding the parameters for calculation
Definition QCCalculationParametersBase.h:294
virtual void read_input_and_commandline_options(World &world, const commandlineparser &parser, const std::string tag)
Definition QCCalculationParametersBase.h:330
void set_user_defined_value(const std::string &key, const T &value)
Definition QCCalculationParametersBase.h:521
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
scalar_type absmax(long *ind=0) const
Return the absolute maximum value (and if ind is non-null, its index) in the Tensor.
Definition tensor.h:1754
float_scalar_type normf() const
Returns the Frobenius norm of the tensor.
Definition tensor.h:1726
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 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:85
std::vector< Vector< double, NDIM > > get_grid() const
Definition lowrankfunction.h:117
dftgrid(const double volume_element, const double radius)
Definition lowrankfunction.h:88
GridBuilder builder
Definition lowrankfunction.h:87
dftgrid(const std::size_t nradial, const std::size_t angular_order, const coord_3d origin=coord_3d())
Definition lowrankfunction.h:109
Definition lowrankfunction.h:57
bool do_print
Definition lowrankfunction.h:81
double get_radius() const
Definition lowrankfunction.h:60
double volume_element
Definition lowrankfunction.h:79
double radius
Definition lowrankfunction.h:80
double get_volume_element() const
Definition lowrankfunction.h:59
virtual ~gridbase()=default
void visualize(const std::string filename, const std::vector< Vector< double, NDIM > > &grid) const
Definition lowrankfunction.h:65
given a molecule, return a suitable grid
Definition lowrankfunction.h:273
std::vector< Vector< double, NDIM > > get_grid() const
Definition lowrankfunction.h:304
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:296
std::vector< Vector< double, NDIM > > centers
Definition lowrankfunction.h:330
std::shared_ptr< gridbase > grid_builder
Definition lowrankfunction.h:331
molecular_grid(const Molecule &molecule, std::shared_ptr< gridbase > grid)
ctor takes molecule and grid builder
Definition lowrankfunction.h:302
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:277
grid with random points around the origin, with a Gaussian distribution
Definition lowrankfunction.h:125
static Vector< double, NDIM > gaussian_random_distribution(const Vector< double, NDIM > origin, double variance)
Definition lowrankfunction.h:186
Vector< double, NDIM > get_origin() const
Definition lowrankfunction.h:173
std::vector< Vector< double, NDIM > > get_grid() const
Definition lowrankfunction.h:133
randomgrid(const double volume_element, const double radius, const Vector< double, NDIM > origin=Vector< double, NDIM >(0.0))
Definition lowrankfunction.h:127
Vector< double, NDIM > origin
Definition lowrankfunction.h:198
double volume() const
Definition lowrankfunction.h:179
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:28
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:2503
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:401
void rr_cholesky(Tensor< T > &A, typename Tensor< T >::scalar_type tol, Tensor< integer > &piv, int &rank)
Compute the rank-revealing Cholesky factorization.
Definition lapack.cc:1203
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:612
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})
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:845
void truncate(World &world, response_space &v, double tol, bool fence)
Definition basic_operators.cc:31
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:707
@ reconstructed
s coeffs at the leaves only
Definition funcdefaults.h:60
response_space transpose(response_space &f)
Definition basic_operators.cc:10
std::vector< Function< T, NDIM > > orthonormalize_canonical(const std::vector< Function< T, NDIM > > &v, const Tensor< T > &ovlp, double lindep)
Definition vmra.h:519
static const Slice _(0,-1, 1)
Tensor< T > inverse(const Tensor< T > &a_in)
invert general square matrix A
Definition lapack.cc:832
@ 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:43
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:667
NDIM & f
Definition mra.h:2481
void error(const char *msg)
Definition world.cc:139
const Function< T, NDIM > & change_tree_state(const Function< T, NDIM > &f, const TreeState finalstate, bool fence=true)
change tree state of a function
Definition mra.h:2807
NDIM const Function< R, NDIM > & g
Definition mra.h:2481
double inner(response_space &a, response_space &b)
Definition response_functions.h:640
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:1810
Vector< double, 3 > coord_3d
Definition funcplot.h:1047
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:1896
double get_size(World &world, const std::vector< Function< T, NDIM > > &v)
Definition vmra.h:1851
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:2066
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 L
Definition rk.cc:46
static const double thresh
Definition rk.cc:45
the low-rank functor is what the LowRankFunction will represent
Definition lowrankfunction.h:414
virtual ~LRFunctorBase()
Definition lowrankfunction.h:416
virtual Tensor< T >::scalar_type norm2() const
Definition lowrankfunction.h:425
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:430
virtual Function< T, LDIM > inner(const Function< T, LDIM > &rhs, const particle< LDIM > p1, const particle< LDIM > p2) const
Definition lowrankfunction.h:420
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:434
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:442
std::vector< Function< T, LDIM > > b
the lo-dim functions
Definition lowrankfunction.h:464
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:457
Tensor< T >::scalar_type norm2() const
Definition lowrankfunction.h:503
std::vector< Function< T, LDIM > > a
Definition lowrankfunction.h:464
T operator()(const Vector< double, NDIM > &r) const
Definition lowrankfunction.h:532
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:468
World & world() const
Definition lowrankfunction.h:467
std::shared_ptr< SeparatedConvolution< T, LDIM > > f12
a two-particle function
Definition lowrankfunction.h:463
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:444
Definition lowrankfunction.h:566
World & world() const
Definition lowrankfunction.h:569
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:573
LRFunctorPure(const Function< T, NDIM > &f)
Definition lowrankfunction.h:568
Function< T, NDIM > f
a hi-dim function
Definition lowrankfunction.h:571
Tensor< T >::scalar_type norm2() const
Definition lowrankfunction.h:585
T operator()(const Vector< double, NDIM > &r) const
Definition lowrankfunction.h:581
Definition lowrankfunction.h:20
std::string rhsfunctiontype() const
Definition lowrankfunction.h:52
double tol() const
Definition lowrankfunction.h:48
int optimize() const
Definition lowrankfunction.h:49
double gamma() const
Definition lowrankfunction.h:46
double volume_element() const
Definition lowrankfunction.h:47
std::string get_tag() const override
Definition lowrankfunction.h:36
std::string gridtype() const
Definition lowrankfunction.h:50
LowRankFunctionParameters()
Definition lowrankfunction.h:22
std::string f12type() const
Definition lowrankfunction.h:53
void read_and_set_derived_values(World &world, const commandlineparser &parser, std::string tag)
Definition lowrankfunction.h:41
std::string orthomethod() const
Definition lowrankfunction.h:51
double radius() const
Definition lowrankfunction.h:45
Definition lowrankfunction.h:203
long total_n
Definition lowrankfunction.h:208
long index
Definition lowrankfunction.h:206
void operator++()
Definition lowrankfunction.h:253
bool operator()() const
Definition lowrankfunction.h:257
Vector< double, NDIM > hivec
Definition lowrankfunction.h:204
long n_per_dim
Definition lowrankfunction.h:207
cartesian_grid(const double volume_per_gridpoint, const double radius)
Definition lowrankfunction.h:211
cartesian_grid & operator=(const cartesian_grid< NDIM > &other)
Definition lowrankfunction.h:231
double volume_per_gridpoint() const
Definition lowrankfunction.h:247
void initialize(const double lo, const double hi)
Definition lowrankfunction.h:237
Vector< double, NDIM > lovec
Definition lowrankfunction.h:204
Vector< double, NDIM > increment
Definition lowrankfunction.h:209
Vector< double, NDIM > get_coordinates() const
Definition lowrankfunction.h:261
std::vector< long > stride
Definition lowrankfunction.h:205
cartesian_grid(const long n_per_dim, const double lo, const double hi)
Definition lowrankfunction.h:221
cartesian_grid(const cartesian_grid< NDIM > &other)
Definition lowrankfunction.h:226
very simple command line parser
Definition commandlineparser.h:15
Definition lowrankfunction.h:336
std::array< int, PDIM > dims
Definition lowrankfunction.h:337
particle(const int p1, const int p2, const int p3)
Definition lowrankfunction.h:364
particle complement() const
return the other particle
Definition lowrankfunction.h:356
static particle particle1()
convenience for particle 1 (the left/first particle)
Definition lowrankfunction.h:343
bool is_first() const
assuming two particles only
Definition lowrankfunction.h:383
std::array< int, PDIM > get_array() const
type conversion to std::array
Definition lowrankfunction.h:377
bool is_last() const
assuming two particles only
Definition lowrankfunction.h:385
particle(const std::vector< int > p)
Definition lowrankfunction.h:365
static particle particle2()
convenience for particle 2 (the right/second particle)
Definition lowrankfunction.h:349
particle(const int p)
Definition lowrankfunction.h:362
std::string str() const
Definition lowrankfunction.h:369
particle()=default
default constructor
particle(const int p1, const int p2)
Definition lowrankfunction.h:363
std::enable_if_t< DUMMYDIM==2, std::tuple< int, int > > get_tuple() const
Definition lowrankfunction.h:393
std::enable_if_t< DUMMYDIM==3, std::tuple< int, int, int > > get_tuple() const
Definition lowrankfunction.h:397
std::enable_if_t< DUMMYDIM==1, std::tuple< int > > get_tuple() const
Definition lowrankfunction.h:389
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.