Skip to content

Commit

Permalink
some cleanup
Browse files Browse the repository at this point in the history
  • Loading branch information
@fbischoff
fbischoff committed Sep 19, 2023
1 parent 6000124 commit 336c966
Show file tree
Hide file tree
Showing 2 changed files with 113 additions and 214 deletions.
187 changes: 52 additions & 135 deletions src/madness/chem/lowrankfunction.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -174,11 +174,8 @@ namespace madness {
double exponent;
double radius=2;
randomgaussian(double exponent, double radius) : exponent(exponent), radius(radius) {
// Vector<double,NDIM> ran; // [0,1]
// RandomVector(NDIM,ran.data());
Vector<double,NDIM> ran= this->gaussian_random_distribution(0,radius);
random_origin=2.0*radius*ran-Vector<double,NDIM>(radius);
// print("origin at ",random_origin, ", exponent",exponent);
}
double operator()(const Vector<double,NDIM>& r) const {
// return exp(-exponent*inner(r-random_origin,r-random_origin));
Expand Down Expand Up @@ -328,64 +325,11 @@ 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()) {
MADNESS_CHECK(functor.f12->info.type==OT_SLATER);
// 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-8));
};
result.push_back(FunctionFactory<T,LDIM>(world).functor(sl));
}
return result;
}


World& world;
double rank_revealing_tol=1.e-8; // rrcd tol
std::string orthomethod="symmetric";
std::string orthomethod="canonical";
bool do_print=true;
std::vector<Function<T,LDIM>> g,h;
LRFunctor lrfunctor;
Expand Down Expand Up @@ -481,36 +425,8 @@ namespace madness {
orthomethod=params.orthomethod();
rank_revealing_tol=params.tol();


// make grid
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);
}


// get sampling grid
std::vector<Vector<double,LDIM>> grid=make_grid(params);
auto Y=Yformer(grid,params.rhsfunctiontype());
t1.tag("Yforming");

Expand All @@ -519,7 +435,7 @@ namespace madness {
g=truncate(orthonormalize_rrcd(Y,ovlp,rank_revealing_tol));
t1.tag("rrcd/truncate/thresh");
auto sz=get_size(world,g);
print("gsize",sz);
if (world.rank()==0 and do_print) print("gsize",sz);
check_orthonormality(g);

if (world.rank()==0 and do_print) {
Expand All @@ -537,10 +453,7 @@ namespace madness {
const double exponent=30.0) {

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

} else if (rhsfunctiontype=="exponential") {
if (rhsfunctiontype=="exponential") {
std::vector<Function<double,LDIM>> omega;
double coeff=std::pow(2.0*exponent/constants::pi,0.25*LDIM);
for (const auto& point : grid) {
Expand All @@ -565,13 +478,21 @@ namespace madness {

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

/// return the size in GByte
double size() const {
double sz=get_size(world,g);
sz+=get_size(world,h);
return sz;
}

Function<T,NDIM> reconstruct() const {
auto fapprox=hartree_product(g[0],h[0]);
for (int i=1; i<g.size(); ++i) fapprox+=hartree_product(g[i],h[i]);
return fapprox;
}


/// orthonormalize the argument vector
std::vector<Function<T,LDIM>> orthonormalize(const std::vector<Function<T,LDIM>>& g) const {

double tol=rank_revealing_tol;
Expand All @@ -596,7 +517,7 @@ namespace madness {
}


/// optimize using Cholesky decomposition
/// optimize the lrf using the lrfunctor

/// @param[in] nopt number of iterations (wrt to Alg. 4.3 in Halko)
void optimize(const long nopt=1) {
Expand All @@ -615,7 +536,7 @@ namespace madness {
}
}

/// after external operations g might not be orthonormal and/or optimal
/// after external operations g might not be orthonormal and/or optimal -- reorthonormalize

/// orthonormalization similar to Bischoff, Harrison, Valeev, JCP 137 104103 (2012), Sec II C 3
/// f =\sum_i g_i h_i
Expand All @@ -634,34 +555,13 @@ namespace madness {
Tensor<T> Xminus=copy(evec_g);
Tensor<T> Yplus=copy(evec_h);
Tensor<T> Yminus=copy(evec_h);
// print("eval_g",eval_g);
// print("eval_h",eval_h);
for (int i=0; i<eval_g.size(); ++i) Xplus(_,i)*=std::pow(eval_g(i),0.5);
for (int i=0; i<eval_g.size(); ++i) Xminus(_,i)*=std::pow(eval_g(i),-0.5);
for (int i=0; i<eval_h.size(); ++i) Yplus(_,i)*=std::pow(eval_h(i),0.5);
for (int i=0; i<eval_h.size(); ++i) Yminus(_,i)*=std::pow(eval_h(i),-0.5);

// test
//{
// auto gortho = transform(world, h, Yminus);
// auto one=matrix_inner(world,gortho,gortho);
// check_orthonormality(one);
//}
//{
// auto gortho = transform(world, g, Xminus);
// auto one=matrix_inner(world,gortho,gortho);
// check_orthonormality(one);
//}
//{
// auto one=madness::inner(Xminus,Xplus,1,1);
// check_orthonormality(one);
// one=madness::inner(Yminus,Yplus,1,1);
// check_orthonormality(one);
//}

Tensor<T> M=madness::inner(Xplus,Yplus,0,0); // (X+)^T Y+
auto [U,s,VT]=svd(M);
// print("s",s);

// truncate
typename Tensor<T>::scalar_type s_accumulated=0.0;
Expand All @@ -673,44 +573,61 @@ namespace madness {
break;
}
}
print("accumulated s",s_accumulated,thresh,s.size(),i);
for (int j=0; j<s.size(); ++j) VT(j,_)*=s[j];
// test
// auto mm=madness::inner(U,VT,1,0);
// double none=(mm-M).normf();
// print("U S V^T - M",none);


Tensor<T> XX=madness::inner(Xminus,U,1,0);
Tensor<T> YY=madness::inner(Yminus,VT,1,1);

g=truncate(transform(world,g,XX));
h=truncate(transform(world,h,YY));
// {
// auto one=matrix_inner(world,g,g);
// check_orthonormality(one);
// }
// {
// auto one=matrix_inner(world,h,h);
// check_orthonormality(one);
// }
}

/// 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;
}


double check_orthonormality(const std::vector<Function<T,LDIM>>& v) const {
Tensor<T> ovlp=matrix_inner(world,v,v);
return check_orthonormality(ovlp);
}

double check_orthonormality(const Tensor<T>& ovlp) const {
timer t(world);
t.do_print=do_print;
Tensor<T> ovlp2=ovlp;
for (int i=0; i<ovlp2.dim(0); ++i) ovlp2(i,i)-=1.0;
print("absmax",ovlp2.absmax());
print("l2",ovlp2.normf()/ovlp2.size());
if (world.rank()==0 and do_print) {
print("absmax",ovlp2.absmax());
print("l2",ovlp2.normf()/ovlp2.size());
}
t.tag("check_orthonoramality");
return ovlp.absmax();
}
Expand Down Expand Up @@ -756,8 +673,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
Loading

0 comments on commit 336c966

Please sign in to comment.