Skip to content

Commit

Permalink
converging the error of the RI approximation by increasing the volume…
Browse files Browse the repository at this point in the history 
… (best value 4) and decreasing the volume per gridpoint (best value 0.02)
  • Loading branch information
@fbischoff
fbischoff committed Aug 18, 2023
1 parent 742ae93 commit d5403f5
Show file tree
Hide file tree
Showing 3 changed files with 141 additions and 39 deletions.
142 changes: 117 additions & 25 deletions src/madness/chem/lowrankfunction.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -87,7 +87,7 @@ namespace madness {
return exp(-exponent*r2);
}

Vector<double,NDIM> gaussian_random_distribution(double mean, double variance) {
static Vector<double,NDIM> gaussian_random_distribution(double mean, double variance) {
std::random_device rd{};
std::mt19937 gen{rd()};
std::normal_distribution<> d{mean, variance};
Expand Down Expand Up @@ -209,6 +209,60 @@ namespace madness {

}

/// inner product: result(1) = \int f(1,2) delta(2) d2

/// @param[in] functor the hidim function
/// @param[in] grid grid points with delta functions
/// @param[in] p1 the variable in f(1,2) to be integrated over
/// @param[in] p2 the variable in rhs to be integrated over (usually all of them)
static std::vector<Function<T,LDIM>> inner(LRFunctor& functor, const std::vector<Vector<double,LDIM>>& grid,
const particle<LDIM> p1, const particle<LDIM> p2) {
std::vector<Function<T,LDIM>> result;
MADNESS_CHECK(functor.has_f() xor functor.has_f12());
MADNESS_CHECK(p1.is_first() xor p1.is_last());
MADNESS_CHECK(p2.is_first());

if (functor.has_f()) {
MADNESS_EXCEPTION("no grid points with an explicit hi-dim function",1);

} else if (functor.has_f12()) {
// functor is now a(1) b(2) f12
// result(1) = \int a(1) f(1,2) b(2) delta(R-2) d2
// = a(1) f(1,R) b(R)
World& world=functor.f12->get_world();
auto premultiply= p1.is_first() ? functor.a : functor.b;
auto postmultiply= p1.is_first() ? functor.b : functor.a;

std::vector<Function<T,LDIM>> f1R= slater_functions_on_grid(world,grid);
print("bla1");
if (premultiply.is_initialized()) {
for (int i=0; i<grid.size(); ++i) f1R[i] = f1R[i]*premultiply(grid[i]);
}
print("bla2");
if (postmultiply.is_initialized()) f1R=f1R*postmultiply;
print("bla3");
result=f1R;

} else {
MADNESS_EXCEPTION("confused functor in LowRankFunction",1);
}
return result;

}


static std::vector<Function<T,LDIM>> slater_functions_on_grid(World& world, const std::vector<Vector<double,LDIM>>& grid) {
std::vector<Function<T,LDIM>> result;
for (const auto& point : grid) {
auto sl=[&point](const Vector<double,LDIM>& r) {
return exp(-sqrt(madness::inner(r-point,r-point)+1.e-12));
};
result.push_back(FunctionFactory<T,LDIM>(world).functor(sl));
}
return result;
}


World& world;
std::vector<Function<T,LDIM>> g,h;
LRFunctor lrfunctor;
Expand Down Expand Up @@ -274,45 +328,58 @@ namespace madness {
return *this;
}

void project(const double volume_per_point, const double radius, const std::string gridtype, std::string rhsfunctiontype) {
long rank=0;
if (gridtype=="random") {
// number of points within radius = variance: 0.67 * #total points = 0.67*rank
// volume of sphere 4/3 pi *r^3
// volume per point=volume/#points in radius = volume / (0.67 * rank)
// rank= volume/(0.67/volume_per_point)
double volume=4.0/3.0*constants::pi *std::pow(radius,3.0);
rank = lround(volume/(0.67*volume_per_point ));
}
project(rank,radius,gridtype,rhsfunctiontype);
}

/// following Halko

/// ||A - Q Q^T A||< epsilon
/// Y = A Omega && Q = QR(Y)
/// || f(1,2) - \sum_i g_i(1) h_i(2) || < epsilon
/// Y_i(1) = \int f(1,2) Omega_i(2) d2 && g_i(1) = QR(Y_i(1)) && h_i(2) = \int g_i^*(1) f(1,2) d1
void project(const long rank, const double radius, const std::string gridtype) {
void project(const long rank, const double radius, const std::string gridtype, std::string rhsfunctiontype) {
timer t1(world);
std::vector<Function<double,LDIM>> omega2;

// make grid
std::vector<Vector<double,LDIM>> grid;
if (gridtype=="random") {
for (long i=0; i<rank; ++i) {
// omega2.push_back(FunctionFactory<double,LDIM>(world).functor(randomgaussian<LDIM>(RandomValue<double>()+10.0,radius)));
omega2.push_back(FunctionFactory<double,LDIM>(world).functor(randomgaussian<LDIM>(50.0,radius)));
for (int i=0; i<rank; ++i) {
auto tmp=randomgaussian<LDIM>::gaussian_random_distribution(0,radius);
auto cell=FunctionDefaults<LDIM>::get_cell();
auto is_in_cell = [&cell](const Vector<double,LDIM>& r) {
for (int d=0; d<LDIM; ++d) if (r[d]<cell(d,0) or r[d]>cell(d,1)) return false;
return true;
};
if (not is_in_cell(tmp)) continue;
grid.push_back(tmp);
}
print("using random gaussian distribution");
} else if (gridtype=="cartesian") {
double volume = 4.0/3.0*constants::pi * std::pow(radius,3.0);
print("volume element in random grid",volume/(0.67*rank));

cartesian_grid<LDIM> cg(9,-radius,radius);
auto c=cg;
c.index=0;
for (; c(); ++c) {
omega2.push_back(FunctionFactory<double,LDIM>(world)
.functor([&c](const Vector<double,LDIM>& r)
{
auto r_rel=r-c.get_coordinates();
return exp(-10*madness::inner(r_rel,r_rel));
}));
}

} else if (gridtype=="cartesian") {
long nperdim=std::lround(std::pow(double(rank),1.0/3.0));
cartesian_grid<LDIM> cg(nperdim,-radius,radius);
for (cg.index=0; cg(); ++cg) grid.push_back(cg.get_coordinates());
print("volume element in cartesian grid",cg.volume_per_gridpoint());
} else {
MADNESS_EXCEPTION("confused gridtype",1);
MADNESS_EXCEPTION("unknown grid type in project",1);
}

print("grid is",gridtype,"with radius",radius,"and",grid.size(),"gridpoints");
print("rhs functions are",rhsfunctiontype);

t1.tag("projection lowdim functions");
print("radius",radius);
print("initial rank",omega2.size());

auto Y=inner(lrfunctor,omega2,p2,p1);
auto Y=Yformer(grid,rhsfunctiontype);
t1.tag("Yforming");

double tol=1.e-12;
Expand All @@ -327,6 +394,31 @@ namespace madness {

}

/// 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)
std::vector<Function<T,LDIM>> Yformer(const std::vector<Vector<double,LDIM>>& grid, const std::string rhsfunctiontype,
const double exponent=30.0) {

std::vector<Function<double,LDIM>> Y;
if (rhsfunctiontype=="delta") {
Y=inner(lrfunctor,grid,p2,p1);

} else if (rhsfunctiontype=="exponential") {
std::vector<Function<double,LDIM>> omega;
for (const auto& point : grid) {
omega.push_back(FunctionFactory<double,LDIM>(world)
.functor([&point,&exponent](const Vector<double,LDIM>& r)
{
auto r_rel=r-point;
return exp(-exponent*madness::inner(r_rel,r_rel));
}));
}
Y=inner(lrfunctor,omega,p2,p1);
} else {
MADNESS_EXCEPTION("confused rhsfunctiontype",1);
}
return Y;
}

long rank() const {return g.size();}

Function<T,NDIM> reconstruct() const {
Expand Down
35 changes: 21 additions & 14 deletions src/madness/chem/test_ccpairfunction.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -21,18 +21,24 @@ struct XParameters : QCCalculationParametersBase {
XParameters() : QCCalculationParametersBase() {

// initialize with: key, value, comment (optional), allowed values (optional)
initialize<double>("radius",5.0,"the radius");
initialize<int>("rank",500,"the number of grid points in random grids");
initialize<double>("radius",2.0,"the radius");
initialize<double>("volume_element",0.1,"volume covered by each grid point");
initialize<long>("rank",500,"the number of grid points in random grids");
initialize<std::string>("gridtype","random","the grid type",{"random","cartesian"});
initialize<std::string>("rhsfunctiontype","delta","the type of function",{"delta","exponential"});
initialize<int>("optimize",0,"number of optimization iterations");
}

void read_and_set_derived_values(World& world, const commandlineparser& parser, std::string tag) {
read_input_and_commandline_options(world,parser,tag);
}

double radius() const {return get<double>("radius");}
int rank() const {return get<int>("rank");}
double volume_element() const {return get<double>("volume_element");}
long rank() const {return get<long>("rank");}
int optimize() const {return get<int>("optimize");}
std::string gridtype() const {return get<std::string>("gridtype");}
std::string rhsfunctiontype() const {return get<std::string>("rhsfunctiontype");}
};


Expand Down Expand Up @@ -161,6 +167,8 @@ int test_lowrank_function3(World& world, XParameters& parameters) {
constexpr std::size_t NDIM=2*LDIM;
print("eps, k, NDIM",FunctionDefaults<NDIM>::get_thresh(),FunctionDefaults<NDIM>::get_k(),NDIM);

parameters.print("grid","end");

Vector<double,LDIM> offset;
offset.fill(0.0);
// Function<double,LDIM> phi1=FunctionFactory<double,LDIM>(world).functor([](const Vector<double,LDIM>& r)
Expand All @@ -181,12 +189,14 @@ int test_lowrank_function3(World& world, XParameters& parameters) {
// };

LowRank<double,6> lrf(f12,copy(phi1),copy(phi2));
lrf.project(parameters.rank(),parameters.radius(),parameters.gridtype());
lrf.optimize(1);
// plot_plane<6>(world,lrf.lrfunctor,"plot_f12_r2",PlotParameters(world).set_plane({"x1","x4"}));
// lrf.project(parameters.rank(),parameters.radius(),parameters.gridtype(),parameters.rhsfunctiontype());
lrf.project(parameters.volume_element(),parameters.radius(),parameters.gridtype(),parameters.rhsfunctiontype());
lrf.optimize(parameters.optimize());
print("lrf.rank()",lrf.rank());

// compare
// \phi(1) \bar \phi(1) = \int phi(1) \phi(2) f(1,2) d2
// \phi(1) \bar \phi(1) = \intn phi(1) \phi(2) f(1,2) d2
// = \int \sum_r g_r(1) h_r(2) d2
// = \sum_r g_r(1) <\phi|h_r>
auto reference = phi1* (*f12)(phi2);
Expand All @@ -198,12 +208,10 @@ int test_lowrank_function3(World& world, XParameters& parameters) {
double resultnorm=result.norm2();
double error=diff.norm2();
print("refnorm, resultnorm, abs. error, rel. error",refnorm, resultnorm, error, error/refnorm);
print("radius, initial/final rank, rel. error",parameters.radius(),parameters.rank(),lrf.rank(), error/refnorm);

plot<LDIM>({reference, result, diff}, "f_and_approx", std::vector<std::string>({"adsf", "asdf", "diff"}));
print("radius, initial/final rank, vol. el, rel. error",parameters.radius(),parameters.rank(),lrf.rank(), parameters.volume_element(), error/refnorm);

plot_plane<6>(world,lrf.lrfunctor,"plot_f12_r2");
plot_plane<6>(world,lrf,"lrf_6d");
// plot_plane<LDIM>(world,{reference, result, diff}, "f_and_approx", PlotParameters(world).set_plane({"x1","x2"}));
// plot_plane<6>(world,lrf,"lrf_6d",PlotParameters(world).set_plane({"x1","x4"}));



Expand Down Expand Up @@ -263,7 +271,7 @@ int test_lowrank_function(World& world) {
f12.reconstruct();
Function<double,NDIM> reference=inner(lhs,rhs,{0},{0});
LowRank<double,NDIM> lrf(rhs);
lrf.project(50,3.0,"random");
lrf.project(50l,3.0,"random","delta");
lrf.optimize();

double ferror=lrf.error();
Expand Down Expand Up @@ -297,7 +305,7 @@ int test_lowrank_function(World& world) {
double radius = 5.0;

LowRank<double, NDIM> lr(copy(f));
lr.project(rank, radius,"random");
lr.project(rank, radius,"random","delta");
double err = lr.error();
print("error in f_approx", err);
lr.orthonormalize(true);
Expand Down Expand Up @@ -1237,7 +1245,6 @@ int main(int argc, char **argv) {
FunctionDefaults<6>::get_thresh());
XParameters parameters;
parameters.read_and_set_derived_values(world,parser,"grid");
parameters.print("grid parameters");
int isuccess=0;
#ifdef USE_GENTENSOR

Expand Down
3 changes: 3 additions & 0 deletions src/madness/mra/funcplot.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -83,6 +83,9 @@ namespace madness {
template<std::size_t NDIM>
Vector<double,NDIM> origin() const {
auto origin_vec=get<std::vector<double>>("origin");
// fill in zeros if the default origin has fewer dimensions than the actual origin
std::size_t missing=NDIM-origin_vec.size();
for (auto i=0; i<missing; ++i) origin_vec.push_back(0.0);
Vector<double,NDIM> o(origin_vec);
return o;
}
Expand Down

0 comments on commit d5403f5

Please sign in to comment.