From 299cd3db3e49280926961371cfbef3b98b4b53df Mon Sep 17 00:00:00 2001 From: fbischoff Date: Wed, 20 Dec 2023 10:49:12 +0100 Subject: [PATCH] use CompositeFactory with vectors of functions --- src/madness/chem/CMakeLists.txt | 2 +- src/madness/chem/projector.h | 114 +++++++++++-------------- src/madness/chem/test_projector.cc | 131 +++++++++++++++++++++++++++++ src/madness/mra/funcimpl.h | 10 ++- src/madness/mra/mra.h | 42 ++++----- src/madness/mra/test6.cc | 2 +- 6 files changed, 209 insertions(+), 92 deletions(-) create mode 100644 src/madness/chem/test_projector.cc diff --git a/src/madness/chem/CMakeLists.txt b/src/madness/chem/CMakeLists.txt index e6f9d51fc6e..14b375a22f1 100644 --- a/src/madness/chem/CMakeLists.txt +++ b/src/madness/chem/CMakeLists.txt @@ -141,7 +141,7 @@ if(BUILD_TESTING) SET(CMAKE_INSTALL_RPATH_USE_LINK_PATH TRUE) # The list of unit test source files set(CHEM_TEST_SOURCES_SHORT test_pointgroupsymmetry.cc test_masks_and_boxes.cc - test_qc.cc test_MolecularOrbitals.cc test_BSHApply.cc) + test_qc.cc test_MolecularOrbitals.cc test_BSHApply.cc test_projector.cc) set(CHEM_TEST_SOURCES_LONG test_localizer.cc test_ccpairfunction.cc) if (LIBXC_FOUND) list(APPEND CHEM_TEST_SOURCES_SHORT test_dft.cc ) diff --git a/src/madness/chem/projector.h b/src/madness/chem/projector.h index fb263c9ba73..d7664cc4fe6 100644 --- a/src/madness/chem/projector.h +++ b/src/madness/chem/projector.h @@ -108,18 +108,21 @@ namespace madness { /// @param[in] f the 6D function to be projected /// @param[in] the particle that is projected (1 or 2) /// @return the projected function - real_function_6d operator()(const real_function_6d& f, const size_t particle) const { - real_function_6d result = real_factory_6d(f.world()); - MADNESS_ASSERT(particle == 1 or particle == 2); + template + typename std::enable_if >::type + operator()(const Function& f, size_t particle1=-1) const { + Function result = FunctionFactory(f.world()); + if (particle1==-1) particle1=particle; + MADNESS_ASSERT(particle1 == 1 or particle1 == 2); for (size_t i = 0; i < mo_ket_.size(); i++) { - real_function_3d tmp1 = mo_ket_[i]; - real_function_3d tmp2 = f.project_out(mo_bra_[i], particle - 1); - real_function_6d tmp12; - if (particle == 1) { - tmp12 = CompositeFactory(f.world()).particle1(copy(tmp1)).particle2(copy(tmp2)); + Function tmp1 = mo_ket_[i]; + Function tmp2 = f.project_out(mo_bra_[i], particle - 1); + Function tmp12; + if (particle1 == 1) { + tmp12 = CompositeFactory(f.world()).particle1(copy(tmp1)).particle2(copy(tmp2)); tmp12.fill_tree(); } else { - tmp12 = CompositeFactory(f.world()).particle1(copy(tmp2)).particle2(copy(tmp1)); + tmp12 = CompositeFactory(f.world()).particle1(copy(tmp2)).particle2(copy(tmp1)); tmp12.fill_tree(); } result += tmp12; @@ -128,7 +131,8 @@ namespace madness { } template - argT operator()(const argT& argument) const { + typename std::enable_if >::value, argT>::type + operator()(const argT& argument) const { return madness::apply(*this,argument); } @@ -275,78 +279,58 @@ namespace madness { // imprecise for lower accuracies // return (f-O1(f)-O2(f)+O1(O2(f))).truncate().reduce_rank(); - const double thresh=FunctionDefaults<2*NDIM>::get_thresh(); - const double tight_thresh=FunctionDefaults<2*NDIM>::get_thresh()*0.1; - // Eq. (A9): g_kl = < k(1) l(2) | f(1,2) > // note no (kl) symmetry here! + reconstruct(world,bra1_,false); + reconstruct(world,bra2_,true); Tensor g_kl(bra1_.size(),bra2_.size()); for (size_t k=0; k kl=CompositeFactory(world) - .particle1(copy(bra1_[k])).particle2(copy(bra2_[l])); + .particle1(bra1_[k]).particle2(bra2_[l]); g_kl(k,l)=inner(f,kl); } } -// if (world.rank()==0) {print(g_kl);}; // Eq. (A12) // project out the mainly first particle: O1 (1 - 1/2 O2) - Function r1=FunctionFactory(world).thresh(tight_thresh); + std::vector> h2(bra1_.size()); + std::vector> h1(ket2_.size()); + reconstruct(world,bra1_,false); + reconstruct(world,bra2_,true); for (size_t k=0; k h2=f.project_out(bra1_[k],0); - - // Eq. (A12) - for (size_t l=0; l::set_thresh(tight_thresh); - r1=(r1+hartree_product(ket1_[k],h2)); - FunctionDefaults<2*NDIM>::set_thresh(thresh); - r1.set_thresh(thresh); - r1.print_size("r1"+stringify(k)); - } + h2[k]=f.project_out(bra1_[k],0); + + // project out the mainly second particle: O2 (1 - 1/2 O1): first term + // Eq. (A11) + std::size_t l=k; + h1[l]=f.project_out(bra2_[l],1); - // project out the mainly second particle: O2 (1 - 1/2 O1) - Function r2=FunctionFactory(world).thresh(tight_thresh); - for (size_t l=0; l h1=f.project_out(bra2_[l],1); - - // Eq. (A13) - for (size_t k=0; k::set_thresh(tight_thresh); - r2=(r2+hartree_product(h1,ket2_[l])); - r2.set_thresh(thresh); - FunctionDefaults<2*NDIM>::set_thresh(thresh); - r2.print_size("r2"+stringify(l)); } - FunctionDefaults<2*NDIM>::set_thresh(tight_thresh); - Function result=(f-r1-r2).truncate().reduce_rank(); - FunctionDefaults<2*NDIM>::set_thresh(thresh); - -// // for debugging purposes only: check orthogonality -// for (size_t k=0; knocc(); ++k) { -// for (size_t l=0; lnocc(); ++l) { -// real_function_6d kl=CompositeFactory(world) -// .particle1(copy(O1_mos[k])).particle2(copy(O2_mos[l])); -// g_kl(k,l)=inner(result,kl); -// } -// } -// if (world.rank()==0) {print(g_kl);}; - - return result; + + // Transforms a vector of functions according to new[i] = sum[j] old[j]*c[j,i] + // project out the mainly first particle: O1 (1 - 1/2 O2): second term + // Eq. (A12), (A13) + h2-=transform(world,ket2_,0.5*transpose(g_kl),false); // ordering g(k,l) is correct + h1-=transform(world,ket1_,0.5*g_kl,true); // ordering g(k,l) is correct + // aka + // for (size_t l=0; l result=copy(f); + result-=hartree_product(ket1_,h2); + result-=hartree_product(h1,ket2_); + result.truncate(); + return result; } private: diff --git a/src/madness/chem/test_projector.cc b/src/madness/chem/test_projector.cc new file mode 100644 index 00000000000..03a57052b42 --- /dev/null +++ b/src/madness/chem/test_projector.cc @@ -0,0 +1,131 @@ +// +// Created by Florian Bischoff on 12/19/23. +// + +#include +#include +#include + + +using namespace madness; + +template +int test_Q12_projector(World& world) { + test_output t1("testing Q12 projector for dimension "+std::to_string(NDIM)); + t1.set_cout_to_terminal(); + constexpr std::size_t LDIM=NDIM/2; + static_assert(NDIM==LDIM*2); + double thresh=FunctionDefaults::get_thresh(); + FunctionDefaults::set_tensor_type(TT_2D); + FunctionDefaults::set_tensor_type(TT_FULL); + + auto g1=[](const Vector& r){return exp(-inner(r,r));}; + auto g2=[](const Vector& r){return exp(-2.0*inner(r,r));}; + auto g_hidim=[](const Vector& r){return exp(-3.0*inner(r,r));}; + + // set up an orthonormal basis for projection + std::vector> vphi; + vphi.push_back(FunctionFactory(world).functor(g1)); + vphi.push_back(FunctionFactory(world).functor(g2)); + vphi=orthonormalize_symmetric(vphi); +// auto S=matrix_inner(world,vphi,vphi); +// print("overlap"); +// print(S); + + StrongOrthogonalityProjector SO(world); + SO.set_spaces(vphi); + + Function f=FunctionFactory(world).functor(g_hidim); + double fnorm=f.norm2(); + print("fnorm",fnorm); + + // check that the projector is indeed a projector + Function f1=SO(f); + double sonorm=f1.norm2(); + print("Q12(f) norm",sonorm); + double refnorm=0.0028346312885398958; // according to mathematica + print("err0",fabs(sonorm-refnorm)); + t1.checkpoint(fabs(sonorm-refnorm) f2=SO(f1); + double err=(f1-f2).norm2()/f.norm2(); + print("err",err); + t1.checkpoint(fabs(err) O1(vphi); + Projector O2(vphi); + O1.set_particle(1); + O2.set_particle(2); + Function f3=f-O1(f)-O2(f)+O1(O2(f)); + double err1=(f1-f3).norm2()/f.norm2(); + print("err1",err1); + t1.checkpoint(fabs(err1)::set_thresh(thresh); + FunctionDefaults<2>::set_thresh(thresh); + FunctionDefaults<3>::set_thresh(thresh); + FunctionDefaults<4>::set_thresh(thresh); + FunctionDefaults<5>::set_thresh(thresh); + FunctionDefaults<6>::set_thresh(thresh); + + FunctionDefaults<1>::set_cubic_cell(-L/2,L/2); + FunctionDefaults<2>::set_cubic_cell(-L/2,L/2); + FunctionDefaults<3>::set_cubic_cell(-L/2,L/2); + FunctionDefaults<4>::set_cubic_cell(-L/2,L/2); + FunctionDefaults<5>::set_cubic_cell(-L/2,L/2); + FunctionDefaults<6>::set_cubic_cell(-L/2,L/2); + + FunctionDefaults<1>::set_k(k); + FunctionDefaults<2>::set_k(k); + FunctionDefaults<3>::set_k(k); + FunctionDefaults<4>::set_k(k); + FunctionDefaults<5>::set_k(k); + FunctionDefaults<6>::set_k(k); + + print("entering testsuite for 6-dimensional functions\n"); + print("k ",k); + print("thresh ",thresh); + print("boxsize ",L); + print("tensor type: ", FunctionDefaults<6>::get_tensor_type()); + print(""); + + + int error=0; + { + error+=test_Q12_projector(world); + error+=test_Q12_projector(world); + error+=test_Q12_projector(world); + + } + + + world.gop.fence(); + finalize(); + + return error; +} + diff --git a/src/madness/mra/funcimpl.h b/src/madness/mra/funcimpl.h index 317703b7f24..149427df009 100644 --- a/src/madness/mra/funcimpl.h +++ b/src/madness/mra/funcimpl.h @@ -720,7 +720,8 @@ namespace madness { accumulate_op(FunctionImpl* f) : impl(f) {} accumulate_op(const accumulate_op& other) = default; void operator()(const Key& key, const coeffT& coeff, const bool& is_leaf) const { - impl->get_coeffs().task(key, &nodeT::accumulate, coeff, impl->get_coeffs(), key, impl->get_tensor_args()); + if (coeff.has_data()) + impl->get_coeffs().task(key, &nodeT::accumulate, coeff, impl->get_coeffs(), key, impl->get_tensor_args()); } template void serialize (Archive& ar) { ar & impl; @@ -5003,7 +5004,7 @@ template /// forces fence double finalize_apply(); - /// after summing upwe need to do some cleanup; + /// after summing up we need to do some cleanup; /// forces fence void finalize_sum(); @@ -6484,11 +6485,12 @@ template std::vector other_s(fcoeff.get_svdtensor().dim_per_vector(otherdim)+1,_); // do the actual contraction + std::vector shape(LDIM,k); for (int r=0; r leaf_op(impl.get(),left,op); impl->hartree_product(left,right,leaf_op,true); impl->finalize_sum(); - this->truncate(0.0,false); + this->truncate(); } @@ -1271,7 +1271,7 @@ namespace madness { hartree_leaf_op leaf_op(impl.get(),k()); impl->hartree_product(left,right,leaf_op,true); impl->finalize_sum(); - this->truncate(0.0,false); + this->truncate(); } @@ -1861,41 +1861,41 @@ namespace madness { return mul(left,right,true); } - /// Performs a Hartree product on the two given low-dimensional functions + /// Performs a Hartree/outer product on the two given low-dimensional function vectors + + /// @return result(x,y) = \sum_i f_i(x) g_i(y) template Function - hartree_product(const Function& left2, const Function& right2) { + hartree_product(const std::vector>& left, const std::vector>& right) { - // we need both sum and difference coeffs for error estimation - Function& left = const_cast< Function& >(left2); - Function& right = const_cast< Function& >(right2); + MADNESS_CHECK_THROW(left.size()==right.size(), "hartree_product: left and right must have same size"); + if (left.size()==0) return Function(); const double thresh=FunctionDefaults::get_thresh(); - FunctionFactory factory=FunctionFactory(left.world()) - .k(left.k()).thresh(thresh); + FunctionFactory factory=FunctionFactory(left.front().world()) + .k(left.front().k()).thresh(thresh); Function result=factory.empty(); - bool same=(left2.get_impl()==right2.get_impl()); - // some prep work - left.make_nonstandard(true,true); - right.make_nonstandard(true,true); + change_tree_state(left,nonstandard_with_leaves); + change_tree_state(right,nonstandard_with_leaves); + std::vector>> vleft=get_impl(left); + std::vector>> vright=get_impl(right); - std::vector>> vleft; - std::vector>> vright; - vleft.push_back(left.get_impl()); - vright.push_back(right.get_impl()); result.do_hartree_product(vleft,vright); - left.standard(false); - if (not same) right.standard(false); - left2.world().gop.fence(); - return result; } + /// Performs a Hartree product on the two given low-dimensional functions + template + Function + hartree_product(const Function& left2, const Function& right2) { + typedef std::vector> vector; + return hartree_product(vector({left2}),vector({right2})); + } /// Performs a Hartree product on the two given low-dimensional functions template diff --git a/src/madness/mra/test6.cc b/src/madness/mra/test6.cc index 3895410f090..c9024a387d3 100644 --- a/src/madness/mra/test6.cc +++ b/src/madness/mra/test6.cc @@ -807,7 +807,7 @@ int main(int argc, char**argv) { error+=test_vector_composite(world,k,thresh); // test(world,k,thresh); error+=test_hartree_product(world,k,thresh); - error+=test_hartree_product(world,k,thresh); +// error+=test_hartree_product(world,k,thresh); error+=test_convolution(world,k,thresh); error+=test_multiply(world,k,thresh); error+=test_add(world,k,thresh);