From 31706ad74ed2fa30a080dfd994bc0eed46f495da Mon Sep 17 00:00:00 2001 From: fbischoff Date: Tue, 19 Sep 2023 15:00:20 +0200 Subject: [PATCH] working on the grid --- src/madness/chem/lowrankfunction.h | 111 +++++++++++++++------ src/madness/chem/test_low_rank_function.cc | 42 +++++++- src/madness/mra/funcplot.h | 4 +- 3 files changed, 119 insertions(+), 38 deletions(-) diff --git a/src/madness/chem/lowrankfunction.h b/src/madness/chem/lowrankfunction.h index d2eb6cc601f..a56388d333f 100644 --- a/src/madness/chem/lowrankfunction.h +++ b/src/madness/chem/lowrankfunction.h @@ -50,15 +50,88 @@ namespace madness { }; + class gridbase { + protected: + double volume_element=0.1; + double radius=3; + }; + + template + class randomgrid : public gridbase { + public: + randomgrid(const double volume_element, const double radius) : gridbase() { + this->volume_element=volume_element; + this->radius=radius; + } + + std::vector> get_grid() const { + std::vector> grid; + long npoint_within_volume=volume()/volume_element; + print("npoint_within_volume",npoint_within_volume); + + auto cell = FunctionDefaults::get_cell(); + auto is_in_cell = [&cell](const Vector& r) { + for (int d = 0; d < NDIM; ++d) if (r[d] < cell(d, 0) or r[d] > cell(d, 1)) return false; + return true; + }; + double rad=radius; + auto is_in_sphere = [&rad](const Vector& r) { + return (r.normf()0 and NDIM<4); + if (NDIM==1) return 2.0*radius; + if (NDIM==2) return constants::pi*radius*radius; + if (NDIM==3) return 4.0 / 3.0 * constants::pi * std::pow(radius, 3.0); + } + + static Vector gaussian_random_distribution(double mean, double variance) { + std::random_device rd{}; + std::mt19937 gen{rd()}; + std::normal_distribution<> d{mean, variance}; + Vector result; + for (int i = 0; i < NDIM; ++i) result[i]=d(gen); + return result; + } + +// double computed_volume_element() const { +// double volume = 4.0 / 3.0 * constants::pi * std::pow(radius, 3.0); +// print("volume element in random grid", volume / (0.67 * rank)); +// } + + }; template - class cgrid { + class sphericalgrid : public gridbase { public: double volume_element=0.1; double radius=3; bool hard_shell=true; - cgrid(const double volume_element, const double radius, bool hard_shell) + sphericalgrid(const double volume_element, const double radius, bool hard_shell) :volume_element(volume_element), radius(radius), hard_shell(hard_shell) { }; @@ -584,34 +657,8 @@ namespace madness { /// make a sampling grid Omega_i(r) for the Halko algorithm std::vector> make_grid(const LowRankFunctionParameters& params) const { - std::vector> grid; - if (params.gridtype()=="random") { - double volume=4.0/3.0*constants::pi *std::pow(params.radius(),3.0); - long rank = lround(volume/(0.67*params.volume_element() )); - for (int i=0; i::gaussian_random_distribution(0,params.radius()); - auto cell=FunctionDefaults::get_cell(); - auto is_in_cell = [&cell](const Vector& r) { - for (int d=0; dcell(d,1)) return false; - return true; - }; - if (not is_in_cell(tmp)) continue; - grid.push_back(tmp); - } - if (world.rank()==0 and do_print) { - double volume = 4.0/3.0*constants::pi * std::pow(params.radius(),3.0); - print("volume element in random grid",volume/(0.67*rank)); - } - - - } else if (params.gridtype()=="spherical") { - cgrid cg(params.volume_element(),params.radius(),params.hard_shell()); - grid=cg.get_coordinates(); - } else { - MADNESS_EXCEPTION("unknown grid type in project",1); - } - - return grid; + randomgrid grid(params.volume_element(),params.radius()); + return grid.get_grid(); } double check_orthonormality(const std::vector>& v) const { @@ -673,8 +720,8 @@ namespace madness { // = \sum_{ij} \int g_i(1) g_j(1) d1 \int h_i(2) h_j(2) d2 // = \sum_{ij} delta_{ij} \int h_i(2) h_j(2) d2 // = \sum_{i} \int h_i(2) h_i(2) d2 -// double zero=check_orthonormality(g); -// if (zero>1.e-10) print("g is not orthonormal",zero); + double zero=check_orthonormality(g); + if (zero>1.e-10) print("g is not orthonormal",zero); double term3=madness::inner(h,h); t.tag("computing term3"); diff --git a/src/madness/chem/test_low_rank_function.cc b/src/madness/chem/test_low_rank_function.cc index dbf25701f49..8c1593e2730 100644 --- a/src/madness/chem/test_low_rank_function.cc +++ b/src/madness/chem/test_low_rank_function.cc @@ -303,6 +303,35 @@ int test_Kcommutator(World& world, LowRankFunctionParameters& parameters) { } +template +int test_grids(World& world, LowRankFunctionParameters& parameters) { + randomgrid g(parameters.volume_element(),parameters.radius()); + g.get_grid(); + + return 0; +} + +template +int test_construction_optimization(World& world, LowRankFunctionParameters parameters) { + parameters.set_user_defined_value("volume_element",0.05); + constexpr std::size_t NDIM=2*LDIM; + test_output t1("LowRankFunction::construction/optimization in dimension "+std::to_string(NDIM)); + t1.set_cout_to_terminal(); + OperatorInfo info(1.0,1.e-6,FunctionDefaults::get_thresh(),OT_SLATER); + auto slater=std::shared_ptr >(new SeparatedConvolution(world,info)); + Function one=FunctionFactory(world).functor([](const Vector& r){return exp(-0.2*inner(r,r));}); + + LowRankFunction lrf(slater,one,one); + lrf.project(parameters); + double error=lrf.l2error(); + t1.checkpoint(error<2.e-2,"l2 error in projection "+std::to_string(error)); + print("l2 error project ",error); + lrf.optimize(); + error=lrf.l2error(); + print("l2 error optimize",error); + t1.checkpoint(error<1.e-2,"l2 error in optimization "+std::to_string(error)); + return t1.end(); +} int main(int argc, char **argv) { @@ -310,7 +339,7 @@ int main(int argc, char **argv) { startup(world, argc, argv); commandlineparser parser(argc, argv); int k = parser.key_exists("k") ? std::atoi(parser.value("k").c_str()) : 6; - double thresh = parser.key_exists("thresh") ? std::stod(parser.value("thresh")) : 1.e-4; + double thresh = parser.key_exists("thresh") ? std::stod(parser.value("thresh")) : 1.e-5; FunctionDefaults<6>::set_tensor_type(TT_2D); FunctionDefaults<1>::set_thresh(thresh); @@ -346,12 +375,17 @@ int main(int argc, char **argv) { try { - parser.set_keyval("geometry", "he"); - parser.print_map(); + isuccess+=test_grids<1>(world,parameters); + isuccess+=test_grids<2>(world,parameters); + isuccess+=test_grids<3>(world,parameters); + isuccess+=test_construction_optimization<1>(world,parameters); + isuccess+=test_construction_optimization<2>(world,parameters); +// isuccess+=test_arithmetic<1>(world,parameters); +// isuccess+=test_arithmetic<2>(world,parameters); // isuccess+=test_lowrank_function(world,parameters); - isuccess+=test_Kcommutator(world,parameters); + isuccess+=test_Kcommutator(world,parameters); } catch (std::exception& e) { madness::print("an error occured"); madness::print(e.what()); diff --git a/src/madness/mra/funcplot.h b/src/madness/mra/funcplot.h index 4a826b2e82a..c0bc02c9e40 100644 --- a/src/madness/mra/funcplot.h +++ b/src/madness/mra/funcplot.h @@ -55,7 +55,7 @@ namespace madness { // initialize with: key, value, comment (optional), allowed values (optional) initialize("zoom",2,"zoom into the simulation cell"); initialize("npoints",151,"number of plot points per dimension"); - initialize>("origin",{0.0,0.0,0.0},"origin of the plot"); + initialize>("origin",{},"origin of the plot"); initialize>("plane",{"x1","x2"},"plot plane: x1, x2, .., x6"); } @@ -84,7 +84,7 @@ namespace madness { Vector origin() const { auto origin_vec=get>("origin"); // fill in zeros if the default origin has fewer dimensions than the actual origin - std::size_t missing=NDIM-origin_vec.size(); + int missing=NDIM-origin_vec.size(); for (auto i=0; i o(origin_vec); return o;