Skip to content

Commit

Permalink
repeated orthonormalization of the initial Y projection
Browse files Browse the repository at this point in the history
  • Loading branch information
@fbischoff
fbischoff committed Sep 8, 2023
1 parent a0593f0 commit c33a65a
Show file tree
Hide file tree
Showing 2 changed files with 96 additions and 43 deletions.
105 changes: 69 additions & 36 deletions src/madness/chem/lowrankfunction.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -106,7 +106,7 @@ namespace madness {
/// a LowRankFunction can be created from a hi-dim function directly, or from a composite like f(1,2) phi(1) psi(2),
/// where f(1,2) is a two-particle function (e.g. a Slater function)
template<typename T, std::size_t NDIM, std::size_t LDIM=NDIM/2>
class LowRank {
class LowRankFunction {
public:

/// what the LowRankFunction will represent
Expand All @@ -123,14 +123,20 @@ namespace madness {
return (f12.get());
}
T operator()(const Vector<double,NDIM>& r) const {
Vector<double,LDIM> first, second;
for (int i=0; i<LDIM; ++i) {
first[i]=r[i];
second[i]=r[i+LDIM];

if (f12->info.type==OT_SLATER) {
double gamma=f12->info.mu;
Vector<double,LDIM> first, second;
for (int i=0; i<LDIM; ++i) {
first[i]=r[i];
second[i]=r[i+LDIM];
}
MADNESS_CHECK(has_f12());
double result=a(first)*b(second)*exp(-gamma*(first-second).normf());
return result;
} else {
return 1.0;
}
MADNESS_CHECK(has_f12());
double result=a(first)*b(second)*exp(-(first-second).normf());
return result;
}
};

Expand Down Expand Up @@ -242,6 +248,7 @@ namespace madness {
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)
Expand Down Expand Up @@ -280,41 +287,42 @@ namespace madness {


World& world;
double rank_revealing_tol; // rrcd tol
double rank_revealing_tol=1.e-12; // rrcd tol
bool do_print=true;
bool stable_power_iteration=true;
std::vector<Function<T,LDIM>> g,h;
LRFunctor lrfunctor;
particle<LDIM> p1=particle<LDIM>::particle1();
particle<LDIM> p2=particle<LDIM>::particle2();

LowRank(World& world) : world(world) {}
LowRankFunction(World& world) : world(world) {}

/// construct from the hi-dim function f
LowRank(const Function<T,NDIM>& f) : LowRank(f.world()) {
LowRankFunction(const Function<T,NDIM>& f) : LowRankFunction(f.world()) {
lrfunctor.f=f;
}

/// construct from the hi-dim function f12*a(1)(b(2)
LowRank(const std::shared_ptr<SeparatedConvolution<T,LDIM>> f12, const Function<T,LDIM>& a,
const Function<T,LDIM>& b) : LowRank(a.world()) {
LowRankFunction(const std::shared_ptr<SeparatedConvolution<T,LDIM>> f12, const Function<T,LDIM>& a,
const Function<T,LDIM>& b) : LowRankFunction(a.world()) {
lrfunctor.a=a;
lrfunctor.b=b;
lrfunctor.f12=f12;
}

LowRank(std::vector<Function<T,LDIM>> g, std::vector<Function<T,LDIM>> h)
LowRankFunction(std::vector<Function<T,LDIM>> g, std::vector<Function<T,LDIM>> h)
: world(g.front().world()), g(g), h(h) {}

LowRank(const LowRank& a) : world(a.world), g(copy(world,a.g)), h(copy(world,a.h)) {} // Copy constructor necessary
LowRankFunction(const LowRankFunction& a) : world(a.world), g(copy(world, a.g)), h(copy(world, a.h)) {} // Copy constructor necessary

LowRank& operator=(const LowRank& f) { // Assignment required for storage in vector
LowRank ff(f);
LowRankFunction& operator=(const LowRankFunction& f) { // Assignment required for storage in vector
LowRankFunction ff(f);
std::swap(ff.g,g);
std::swap(ff.h,h);
return *this;
}

/// function evaluation
T operator()(const Vector<double,NDIM>& r) const {
Vector<double,LDIM> first, second;
for (int i=0; i<LDIM; ++i) {
Expand All @@ -326,25 +334,42 @@ namespace madness {
return result;
}

LowRank operator-(const LowRank& b) const { // Operator- necessary
return LowRank(g-b.g,h-b.h);
/*
* arithmetic section
*/

/// addition
LowRankFunction operator+(const LowRankFunction& b) const {
return LowRankFunction(g + b.g, h + b.h);
}
/// subtraction
LowRankFunction operator-(const LowRankFunction& b) const {
return LowRankFunction(g - b.g, h - b.h);
}

LowRank& operator+=(const LowRank& b) { // Operator+= necessary
/// in-place addition
LowRankFunction& operator+=(const LowRankFunction& b) {
g+=b.g;
h+=b.h;
return *this;
}

LowRank operator*(double a) const { // Scale by a constant necessary
return LowRank(g*a,h);
/// in-place subtraction
LowRankFunction& operator-=(const LowRankFunction& b) {
g-=b.g;
h-=b.h;
return *this;
}

LowRank multiply(const std::vector<Function<T,LDIM>>& vec, const long particle) {
auto result=*this; // deep copy
if (particle==0) result.g=g*vec;
if (particle==1) result.h=h*vec;
return *this;
/// scale by a scalar
template<typename Q>
LowRankFunction operator*(const Q a) const {
return LowRankFunction<TensorResultType<T,Q>,NDIM>(g * a, Q(h));
}

/// in-place scale by a scalar (no type conversion)
LowRankFunction operator*(const T a) const {
return LowRankFunction(g * a, h);
}

void project(const double volume_per_point, const double radius, const std::string gridtype, std::string rhsfunctiontype,
Expand Down Expand Up @@ -411,8 +436,16 @@ namespace madness {
std::ostringstream oss;
oss << std::scientific << std::setprecision(1) << rank_revealing_tol;
std::string scientificString = oss.str();
g=orthonormalize_rrcd(Y,rank_revealing_tol);
t1.tag("Y orthonormalizing with rank_revealing_tol "+scientificString);
double err=1.0;
double tight_thresh=FunctionDefaults<3>::get_thresh()*0.1;
g=Y;
while (err>1.e-10) {
g=truncate(orthonormalize_rrcd(g,rank_revealing_tol),std::min(1.e-3,tight_thresh*100.0));
err=check_orthonormality(g);
print("error of non-orthonormality",err);
t1.tag("Y orthonormalizing with rank_revealing_tol loose_thresh "+scientificString);
}
g=truncate(g);

if (world.rank()==0 and do_print) {
print("Y.size()",Y.size());
Expand Down Expand Up @@ -482,7 +515,7 @@ namespace madness {
double check_orthonormality(const std::vector<Function<T,LDIM>>& v) const {
Tensor<T> ovlp=matrix_inner(world,v,v);
for (int i=0; i<ovlp.dim(0); ++i) ovlp(i,i)-=1.0;
return ovlp.normf()/ovlp.dim(0);
return ovlp.normf()/ovlp.size();
}

double explicit_error() const {
Expand Down Expand Up @@ -531,7 +564,7 @@ namespace madness {
double term3=madness::inner(h,h);
t.tag("computing term3");

double error=sqrt(term1-2.0*term2+term3)/term1;
double error=sqrt(term1-2.0*term2+term3)/sqrt(term1);
if (world.rank()==0 and do_print) {
print("term1,2,3, error",term1, term2, term3, " --",error);
}
Expand All @@ -543,19 +576,19 @@ namespace madness {

// This interface is necessary to compute inner products
template<typename T, std::size_t NDIM>
double inner(const LowRank<T,NDIM>& a, const LowRank<T,NDIM>& b) {
double inner(const LowRankFunction<T,NDIM>& a, const LowRankFunction<T,NDIM>& b) {
World& world=a.world;
return (matrix_inner(world,a.g,b.g).emul(matrix_inner(world,a.h,b.h))).sum();
}



template<typename T, std::size_t NDIM>
LowRank<T,NDIM> inner(const Function<T,NDIM>& lhs, const LowRank<T,NDIM>& rhs, const std::tuple<int> v1, const std::tuple<int> v2) {
LowRankFunction<T,NDIM> inner(const Function<T,NDIM>& lhs, const LowRankFunction<T,NDIM>& rhs, const std::tuple<int> v1, const std::tuple<int> v2) {
World& world=rhs.world;
// int lhs(1,2) rhs(2,3) d2 = \sum \int lhs(1,2) g_i(2) h_i(3) d2
// = \sum \int lhs(1,2) g_i(2) d2 h_i(3)
LowRank<T,NDIM+NDIM-2> result(world);
LowRankFunction<T, NDIM + NDIM - 2> result(world);
result.h=rhs.h;
decltype(rhs.g) g;
for (int i=0; i<rhs.rank(); ++i) {
Expand All @@ -566,11 +599,11 @@ namespace madness {
}

template<typename T, std::size_t NDIM>
LowRank<T,NDIM> inner(const LowRank<T,NDIM>& f, const Function<T,NDIM>& g, const std::tuple<int> v1, const std::tuple<int> v2) {
LowRankFunction<T,NDIM> inner(const LowRankFunction<T,NDIM>& f, const Function<T,NDIM>& g, const std::tuple<int> v1, const std::tuple<int> v2) {
World& world=f.world;
// int f(1,2) k(2,3) d2 = \sum \int g_i(1) h_i(2) k(2,3) d2
// = \sum g_i(1) \int h_i(2) k(2,3) d2
LowRank<T,NDIM+NDIM-2> result(world);
LowRankFunction<T, NDIM + NDIM - 2> result(world);
result.g=f.g;
decltype(f.h) h;
for (int i=0; i<f.rank(); ++i) {
Expand Down
34 changes: 27 additions & 7 deletions src/madness/chem/test_low_rank_function.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,7 +5,6 @@

#include<madness.h>
#include<madness/chem/lowrankfunction.h>
#include<madness/chem/electronic_correlation_factor.h>


#include<madness/world/test_utilities.h>
Expand All @@ -25,7 +24,9 @@ struct LowRankFunctionParameters : QCCalculationParametersBase {
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<bool>("stable_power_iteration",true,"use stable power iteration algorithm (orthonormalize)");
initialize<double>("tol",1.e-12,"rank-reduced choleski tolerance");
initialize<double>("tol",1.e-12,"rank-reduced cholesky tolerance");
initialize<std::string>("f12type","Slater","correlation factor",{"Slater","SlaterF12"});
initialize<std::string>("transpose","slater2","transpose of the matrix",{"slater1","slater2"});
initialize<std::string>("gridtype","random","the grid type",{"random","cartesian"});
initialize<std::string>("rhsfunctiontype","exponential","the type of function",{"delta","exponential"});
initialize<int>("optimize",1,"number of optimization iterations");
Expand Down Expand Up @@ -60,15 +61,18 @@ int test_lowrank_function(World& world, LowRankFunctionParameters& parameters) {
print("eps, k, NDIM, id",FunctionDefaults<NDIM>::get_thresh(),FunctionDefaults<NDIM>::get_k(),NDIM,id);

parameters.print("grid","end");
std::string f12type=parameters.get<std::string>("f12type");
std::string transpose=parameters.get<std::string>("transpose");

json j;
std::string jsonfilename="test_low_rank_function."+id+".json";
j["radius"]=parameters.radius();
j["f12type"]="SlaterF12";
j["f12type"]=parameters.get<std::string>("f12type");
j["gamma"]=parameters.gamma();
j["volume_element"]=parameters.volume_element();
j["tol"]=parameters.tol();
j["rank"]=parameters.rank();
j["transpose"]=transpose;
j["stable_power_iteration"]=parameters.stable_power_iteration();
j["gridtype"]=parameters.gridtype();
j["rhsfunctiontype"]=parameters.rhsfunctiontype();
Expand All @@ -79,17 +83,33 @@ int test_lowrank_function(World& world, LowRankFunctionParameters& parameters) {

Vector<double,LDIM> offset;
offset.fill(0.0);
// Function<double,LDIM> phi1=FunctionFactory<double,LDIM>(world).functor([](const Vector<double,LDIM>& r)
// { return exp(-r.normf());});
Function<double,LDIM> phi1=FunctionFactory<double,LDIM>(world).functor([](const Vector<double,LDIM>& r)
{ return 1.0;});
Function<double,LDIM> phi2=FunctionFactory<double,LDIM>(world).functor([&offset](const Vector<double,LDIM>& r)
{ return exp(-1.0*(r-offset).normf());});
if (transpose=="slater1") std::swap(phi1,phi2);
{
double n1 = phi1.norm2();
double n2 = phi2.norm2();
bool first_one = (fabs(phi1({1, 1, 1}) - 1.0) < 1.e-6);
if (world.rank() == 0) {
if (first_one) print("1(1) phi(2)");
else print("phi(1) 1(2)");
print("norms", n1, n2);
}
}

Function<double,LDIM> one=FunctionFactory<double,LDIM>(world)
.functor([](const Vector<double,LDIM>& r) { return 1.0;});


std::shared_ptr<real_convolution_3d> f12(SlaterF12OperatorPtr(world,parameters.gamma(),1.e-6,FunctionDefaults<LDIM>::get_thresh()));
std::shared_ptr<real_convolution_3d> f12;
if (f12type=="slaterf12") f12.reset(SlaterF12OperatorPtr(world,parameters.gamma(),1.e-6,FunctionDefaults<LDIM>::get_thresh()));
else if (f12type=="slater") f12.reset(SlaterOperatorPtr(world,parameters.gamma(),1.e-6,FunctionDefaults<LDIM>::get_thresh()));
else {
MADNESS_EXCEPTION(std::string("unknown f12type"+f12type).c_str(),1);
}


auto compute_result = [&world, &one](const auto& lrf) {
real_function_3d result=real_factory_3d(world);
Expand Down Expand Up @@ -119,7 +139,7 @@ int test_lowrank_function(World& world, LowRankFunctionParameters& parameters) {
print("reference.norm2() = int f12 phi2 d2",n2);
output(0.0,0.0,0.0,0.0,0.0,0.0);

LowRank<double,6> lrf(f12,copy(phi1),copy(phi2));
LowRankFunction<double,6> lrf(f12, copy(phi1), copy(phi2));
lrf.stable_power_iteration=parameters.stable_power_iteration();
plot_plane<6>(world,lrf.lrfunctor,"plot_original."+id,PlotParameters(world).set_plane({"x1","x4"}));
double cpu0=cpu_time();
Expand Down

0 comments on commit c33a65a

Please sign in to comment.