Skip to content

Commit

Permalink
working on the grid
Browse files Browse the repository at this point in the history
  • Loading branch information
@fbischoff
fbischoff committed Sep 20, 2023
1 parent 336c966 commit 31706ad
Show file tree
Hide file tree
Showing 3 changed files with 119 additions and 38 deletions.
111 changes: 79 additions & 32 deletions src/madness/chem/lowrankfunction.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -50,15 +50,88 @@ namespace madness {
};


class gridbase {
protected:
double volume_element=0.1;
double radius=3;
};

template<std::size_t NDIM>
class randomgrid : public gridbase {
public:
randomgrid(const double volume_element, const double radius) : gridbase() {
this->volume_element=volume_element;
this->radius=radius;
}

std::vector<Vector<double,NDIM>> get_grid() const {
std::vector<Vector<double, NDIM>> grid;
long npoint_within_volume=volume()/volume_element;
print("npoint_within_volume",npoint_within_volume);

auto cell = FunctionDefaults<NDIM>::get_cell();
auto is_in_cell = [&cell](const Vector<double, NDIM>& 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<double, NDIM>& r) {
return (r.normf()<rad);
};

// set variance such that about 70% of all grid points sits within the radius
double variance=radius;
if (NDIM==2) variance=0.6*radius;
if (NDIM==3) variance=0.5*radius;
long maxrank=10*npoint_within_volume;
long rank=0;
for (int r=0; r<maxrank; ++r) {
auto tmp = gaussian_random_distribution(0, variance);
if (not is_in_cell(tmp)) continue;
if (is_in_sphere(tmp)) ++rank;
grid.push_back(tmp);
if (rank==npoint_within_volume) break;
}
print("grid points in volume ",rank);
print("total grid points ",grid.size());
print("ratio ",rank/double(grid.size()));
print("volume element ",volume()/rank);
return grid;
}

private:

double volume() const {
MADNESS_CHECK(NDIM>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<double,NDIM> gaussian_random_distribution(double mean, double variance) {
std::random_device rd{};
std::mt19937 gen{rd()};
std::normal_distribution<> d{mean, variance};
Vector<double,NDIM> 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<std::size_t NDIM>
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) {
};

Expand Down Expand Up @@ -584,34 +657,8 @@ namespace madness {
/// make a sampling grid Omega_i(r) for the Halko algorithm
std::vector<Vector<double,LDIM>> make_grid(const LowRankFunctionParameters& params) const {

std::vector<Vector<double,LDIM>> 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<rank; ++i) {
auto tmp=randomgaussian<LDIM>::gaussian_random_distribution(0,params.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);
}
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<LDIM> 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<LDIM> grid(params.volume_element(),params.radius());
return grid.get_grid();
}

double check_orthonormality(const std::vector<Function<T,LDIM>>& v) const {
Expand Down Expand Up @@ -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");

Expand Down
42 changes: 38 additions & 4 deletions src/madness/chem/test_low_rank_function.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -303,14 +303,43 @@ int test_Kcommutator(World& world, LowRankFunctionParameters& parameters) {

}

template<std::size_t LDIM>
int test_grids(World& world, LowRankFunctionParameters& parameters) {
randomgrid<LDIM> g(parameters.volume_element(),parameters.radius());
g.get_grid();

return 0;
}

template<std::size_t LDIM>
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<LDIM>::get_thresh(),OT_SLATER);
auto slater=std::shared_ptr<SeparatedConvolution<double,LDIM> >(new SeparatedConvolution<double,LDIM>(world,info));
Function<double,LDIM> one=FunctionFactory<double,LDIM>(world).functor([](const Vector<double,LDIM>& r){return exp(-0.2*inner(r,r));});

LowRankFunction<double,NDIM> 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) {

madness::World& world = madness::initialize(argc, 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);
Expand Down Expand Up @@ -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());
Expand Down
4 changes: 2 additions & 2 deletions src/madness/mra/funcplot.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -55,7 +55,7 @@ namespace madness {
// initialize with: key, value, comment (optional), allowed values (optional)
initialize<double>("zoom",2,"zoom into the simulation cell");
initialize<long>("npoints",151,"number of plot points per dimension");
initialize<std::vector<double>>("origin",{0.0,0.0,0.0},"origin of the plot");
initialize<std::vector<double>>("origin",{},"origin of the plot");
initialize<std::vector<std::string>>("plane",{"x1","x2"},"plot plane: x1, x2, .., x6");
}

Expand Down Expand Up @@ -84,7 +84,7 @@ namespace madness {
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();
int 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 31706ad

Please sign in to comment.