Merge pull request #98 from stigrj/gauge

CartesianConvolution

cmake/custom/fetch_mrcpp.cmake

@@ -2,7 +2,7 @@
 set(_build_type ${CMAKE_BUILD_TYPE})
 set(CMAKE_BUILD_TYPE Release)
 
-find_package(MRCPP 1.4 CONFIG QUIET)
+find_package(MRCPP 1.5 CONFIG QUIET)
 
 # whether MRCPP was fetched and built locally
 set(MRCPP_FETCHED FALSE)
@@ -18,7 +18,7 @@ else()
     GIT_REPOSITORY
       https://github.com/MRChemSoft/mrcpp.git
     GIT_TAG
-      1a6d450ac356feffba5c5cb2d02b92a875dbb9aa
+      v1.5.0
     )
 
   set(CMAKE_CXX_COMPILER ${CMAKE_CXX_COMPILER})

src/vampyr/operators/convolutions.h

@@ -5,10 +5,12 @@
 #include <MRCPP/operators/HelmholtzOperator.h>
 #include <MRCPP/operators/IdentityConvolution.h>
 #include <MRCPP/operators/PoissonOperator.h>
+#include <MRCPP/operators/CartesianConvolution.h>
 #include <MRCPP/treebuilders/apply.h>
 
 namespace vampyr {
 
+void cartesian_convolution(pybind11::module &);
 void helmholtz_operator(pybind11::module &);
 void poisson_operator(pybind11::module &);
 
@@ -17,7 +19,19 @@ template <int D> void convolutions(pybind11::module &m) {
     using namespace mrcpp;
     using namespace pybind11::literals;
 
-    py::class_<ConvolutionOperator<D>>(m, "ConvolutionOperator");
+    py::class_<ConvolutionOperator<D>>(m, "ConvolutionOperator")
+        .def(py::init<const MultiResolutionAnalysis<D> &, GaussExp<1> &, double>(),
+             "mra"_a,
+             "kernel"_a,
+             "prec"_a)
+        .def(
+            "__call__",
+            [](ConvolutionOperator<D> &O, FunctionTree<D> *inp) {
+                auto out = std::make_unique<FunctionTree<D>>(inp->getMRA());
+                apply<D>(O.getBuildPrec(), *out, O, *inp);
+                return out;
+            },
+            "inp"_a);
 
     py::class_<IdentityConvolution<D>, ConvolutionOperator<D>>(m, "IdentityConvolution")
         .def(py::init<const MultiResolutionAnalysis<D> &, double>(), "mra"_a, "prec"_a)
@@ -35,10 +49,32 @@ template <int D> void convolutions(pybind11::module &m) {
             },
             "inp"_a);
 
+    if constexpr (D == 3) cartesian_convolution(m);
     if constexpr (D == 3) helmholtz_operator(m);
     if constexpr (D == 3) poisson_operator(m);
 }
 
+void cartesian_convolution(pybind11::module &m) {
+    namespace py = pybind11;
+    using namespace mrcpp;
+    using namespace pybind11::literals;
+
+    py::class_<CartesianConvolution, ConvolutionOperator<3>>(m, "CartesianConvolution")
+        .def(py::init<const MultiResolutionAnalysis<3> &, GaussExp<1> &, double>(),
+             "mra"_a,
+             "kernel"_a,
+             "prec"_a)
+        .def(
+            "__call__",
+            [](CartesianConvolution &O, FunctionTree<3> *inp) {
+                auto out = std::make_unique<FunctionTree<3>>(inp->getMRA());
+                apply<3>(O.getBuildPrec(), *out, O, *inp);
+                return out;
+            },
+            "inp"_a)
+        .def("setCartesianComponents", &CartesianConvolution::setCartesianComponents);
+}
+
 void poisson_operator(pybind11::module &m) {
     namespace py = pybind11;
     using namespace mrcpp;

src/vampyr/operators/derivatives.h

@@ -20,7 +20,6 @@ template <int D> void derivatives(pybind11::module &m) {
         An abstract base class for derivative operators
     )mydelimiter")
         // clang-format on
-        .def(py::init<const MultiResolutionAnalysis<D> &>(), "mra"_a)
         .def("getOrder", &DerivativeOperator<D>::getOrder)
         .def(
             "__call__",

src/vampyr/tests/test_convolutions_1d.py

@@ -17,6 +17,26 @@
 ffunc = vp.GaussFunc(alpha=alpha, beta=beta, position=r0)
 
 
+def test_GaussKernel():
+    beta = 1.0e5
+    alpha = (beta / np.pi) ** (1.0 / 2.0)
+    ifunc = vp.GaussFunc(alpha=alpha, beta=beta)
+    iexp = vp.GaussExp()
+    iexp.append(ifunc)
+    I = vp.ConvolutionOperator(mra, iexp, prec=epsilon)
+
+    ftree = vp.FunctionTree(mra)
+    vp.advanced.build_grid(out=ftree, inp=ffunc)
+    vp.advanced.project(prec=epsilon, out=ftree, inp=ffunc)
+
+    gtree = vp.FunctionTree(mra)
+    vp.advanced.apply(prec=epsilon, out=gtree, oper=I, inp=ftree)
+    assert gtree.integrate() == pytest.approx(ftree.integrate(), rel=epsilon)
+
+    gtree2 = I(ftree)
+    assert gtree2.integrate() == pytest.approx(ftree.integrate(), rel=epsilon)
+
+
 def test_Identity():
     I = vp.IdentityConvolution(mra, prec=epsilon)
 

src/vampyr/tests/test_convolutions_3d.py

@@ -2,6 +2,7 @@
 import pytest
 
 from vampyr import vampyr3d as vp
+from vampyr import vampyr1d as vp1
 
 epsilon = 1.0e-3
 mu = epsilon / 10
@@ -22,6 +23,46 @@
 vp.advanced.build_grid(out=ftree, inp=ffunc)
 vp.advanced.project(prec=epsilon, out=ftree, inp=ffunc)
 
+def test_GaussKernel():
+    b = 1.0e4
+    a = (b / np.pi) ** (D / 2.0)
+    ifunc = vp1.GaussFunc(alpha=a, beta=b)
+    iexp = vp1.GaussExp()
+    iexp.append(ifunc)
+    I = vp.ConvolutionOperator(mra, iexp, prec=epsilon)
+
+    ftree = vp.FunctionTree(mra)
+    vp.advanced.build_grid(out=ftree, inp=ffunc)
+    vp.advanced.project(prec=epsilon, out=ftree, inp=ffunc)
+
+    gtree = vp.FunctionTree(mra)
+    vp.advanced.apply(prec=epsilon, out=gtree, oper=I, inp=ftree)
+    assert gtree.integrate() == pytest.approx(ftree.integrate(), rel=epsilon)
+
+    gtree2 = I(ftree)
+    assert gtree2.integrate() == pytest.approx(ftree.integrate(), rel=epsilon)
+
+
+def test_CartesianConvolution():
+    b = 1.0e4
+    a = (b / np.pi) ** (D / 2.0)
+    ifunc = vp1.GaussFunc(alpha=a, beta=b)
+    iexp = vp1.GaussExp()
+    iexp.append(ifunc)
+    O = vp.CartesianConvolution(mra, iexp, prec=epsilon)
+    O.setCartesianComponents(0, 0, 0)
+
+    ftree = vp.FunctionTree(mra)
+    vp.advanced.build_grid(out=ftree, inp=ffunc)
+    vp.advanced.project(prec=epsilon, out=ftree, inp=ffunc)
+
+    gtree = vp.FunctionTree(mra)
+    vp.advanced.apply(prec=epsilon, out=gtree, oper=O, inp=ftree)
+    assert gtree.integrate() == pytest.approx(ftree.integrate(), rel=epsilon)
+
+    gtree2 = O(ftree)
+    assert gtree2.integrate() == pytest.approx(ftree.integrate(), rel=epsilon)
+
 
 def test_Identity():
     I = vp.IdentityConvolution(mra, prec=epsilon)