From 1b3b3e418a12496823d81226213bdb9578b56565 Mon Sep 17 00:00:00 2001 From: Jonathon Misiewicz Date: Mon, 14 Oct 2024 09:34:27 -0400 Subject: [PATCH] testfilter --- src/examples/periodic/CMakeLists.txt | 3 +- src/examples/periodic/testfilter.cc | 122 +++++++++++++++++++++++++++ 2 files changed, 124 insertions(+), 1 deletion(-) create mode 100644 src/examples/periodic/testfilter.cc diff --git a/src/examples/periodic/CMakeLists.txt b/src/examples/periodic/CMakeLists.txt index 92c9ea76e9c..0361c169bea 100644 --- a/src/examples/periodic/CMakeLists.txt +++ b/src/examples/periodic/CMakeLists.txt @@ -1,7 +1,8 @@ # src/examples/periodic set(EXAMPLE_SOURCES - testpercoul) + testpercoul + testfilter) # Create executables for example applications foreach(example ${EXAMPLE_SOURCES}) diff --git a/src/examples/periodic/testfilter.cc b/src/examples/periodic/testfilter.cc new file mode 100644 index 00000000000..0ddd58098c6 --- /dev/null +++ b/src/examples/periodic/testfilter.cc @@ -0,0 +1,122 @@ +#include +#include + +const double L = 10; + +static double gaussian_3d(const madness::coord_3d& r, const double expnt) { + const double x=r[0], y=r[1], z=r[2]; + const double coeff = pow(expnt/M_PI, 1.5); + return coeff * exp(-expnt * (x*x + y*y + z*z)); +} + +static double rho_electronic_3d(const madness::coord_3d& r) { return gaussian_3d(r, 1); } +static double rho_nuclear_3d(const madness::coord_3d& r) { return -gaussian_3d(r, 10000); } + +static double rho_gaussian_func_3d(const madness::coord_3d& r) { + return rho_electronic_3d(r) + rho_nuclear_3d(r); +} + +template +void filter_moments_inplace(madness::Function& f, const int k, const bool fence=true) { + auto filter_op = [&](madness::Tensor &coeff) { + if constexpr (NDIM == 3) { + if (k >= 0) { + coeff(0, 0, 0) = 0.; + } + if (k >= 1) { + coeff(1, 0, 0) = 0.; + coeff(0, 1, 0) = 0.; + coeff(0, 0, 1) = 0.; + } + if (k >= 2) { + coeff(2, 0, 0) = 0.; + coeff(0, 2, 0) = 0.; + coeff(0, 0, 2) = 0.; + coeff(1, 1, 0) = 0.; + coeff(1, 0, 1) = 0.; + coeff(0, 1, 1) = 0.; + } + if (k >= 3) + abort();// TODO implement support for higher moments + } else + static_assert("filter_moments_inplace at the moment is implemented for NDIM=3"); + }; + + // on [-L,L] normalized scaling function k is Sqrt[(2 k + 1)/(2 L)] LegendreP[k, x/L] + // the coefficient of x^k in LegendreP[k, x] is 2^(-k) Binomial[2 k, k], hence + // coefficient of x^k in scaling function k is Sqrt[(2 k + 1)/(2 L)] Binomial[2 k, k]/(2 L)^k + // thus if all moments of up to k-1 vanish, k-th moment (=expectation value of x^k) of + // scaling function k is its coefficient times Sqrt[(2 L)/ (2 k + 1)] (2 L)^k / Binomial[2 k, k] + f.unaryop_coeff([&](const madness::Key &key, madness::Tensor &coeff) { + if (f.is_reconstructed()) { + filter_op(coeff); + } else if (f.is_compressed()) { + // for compressed form only need 1 node, but public interface does not allow mutable access to the data + if (key == f.get_impl()->key0()) { + filter_op(coeff); + } + } else { + MADNESS_EXCEPTION("filter_moments_inplace(f): f must be either compressed or reconstructed", 1); + } + }, + fence); +} + +int main(int argc, char**argv) { + using namespace madness; + { + auto &world = initialize(argc, argv); + startup(world, argc, argv, true); + + // Function defaults + int k = 11; + double thresh = 1e-9; + double eps = 1e-9; + FunctionDefaults<3>::set_k(k); + FunctionDefaults<3>::set_cubic_cell(-L / 2, L / 2); + FunctionDefaults<3>::set_thresh(thresh); + FunctionDefaults<3>::set_refine(true); + FunctionDefaults<3>::set_initial_level(2); + FunctionDefaults<3>::set_truncate_mode(1); + + Function rho, rhoE, rhoN, V_periodic, V_periodicE, V_periodicN; + + { + constexpr double filter_level = 0; + printf("building gaussian diff charge distribution ...\n\n"); + std::vector special_pt{coord_3d(std::array{0.0, 0.0, 0.0})}; + rho = FunctionFactory(world).special_points(special_pt).initial_level(4).f(rho_gaussian_func_3d); + filter_moments_inplace(rho, filter_level); + rho.truncate(); + rhoE = FunctionFactory(world).special_points(special_pt).initial_level(4).f(rho_electronic_3d); + filter_moments_inplace(rhoE, filter_level); + rhoE.truncate(); + rhoN = FunctionFactory(world).special_points(special_pt).initial_level(4).f(rho_nuclear_3d); + filter_moments_inplace(rhoN, filter_level); + rhoN.truncate(); + + BoundaryConditions<3> bc_periodic(BC_PERIODIC); + SeparatedConvolution pop = CoulombOperator(world, 1e-4, eps, bc_periodic); + printf("applying periodic operator ...\n\n"); + V_periodic = apply(pop, rho); + V_periodic.truncate(); + V_periodicE = apply(pop, rhoE); + V_periodicE.truncate(); + V_periodicN = apply(pop, rhoN); + V_periodicN.truncate(); + + V_periodic.reconstruct(); + V_periodicE.reconstruct(); + V_periodicN.reconstruct(); + } + + double bstep = L / 1000.0; + printf(" z\t\tV_(E+N)[z]\tV_E[z]\t\tV_N[z]\t\tV_E[z]+V_n[z]\trho\t\terror\n"); + for (int i = 0; i < 1001; i++) { + double x = -L / 2 + i * bstep; + coord_3d p(std::array{0,0,x}); + printf("%.3f\t\t%.8f\t%.8f\t%.8f\t%.8f\t%.8f\t%.8f\n", p[2], V_periodic(p), V_periodicE(p), V_periodicN(p), V_periodicE(p) + V_periodicN(p), rho(p), V_periodic(p) - V_periodicN(p) - V_periodicE(p)); + } + } + madness::finalize(); +}