diff --git a/.github/workflows/build-test.yml b/.github/workflows/build-test.yml
index 21a9034ee..4c5011800 100644
--- a/.github/workflows/build-test.yml
+++ b/.github/workflows/build-test.yml
@@ -27,10 +27,10 @@ jobs:
         os: [ubuntu-latest]  #, macos-12]
 
     steps:
-    - uses: actions/checkout@v2
+    - uses: actions/checkout@v4
 
     - name: Cache conda
-      uses: actions/cache@v1
+      uses: actions/cache@v4
       env:
         CACHE_NUMBER: 0  # Increase this value to reset cache if .github/mrcpp-gha.yml has not changed
       with:
@@ -45,7 +45,7 @@ jobs:
           activate-environment: mrcpp-gha
           environment-file: .github/mrcpp-gha.yml
           channel-priority: true
-          python-version: 3.6
+          python-version: 3.11
           use-only-tar-bz2: true  # IMPORTANT: This needs to be set for caching to work properly!
 
     - name: Configure
diff --git a/.github/workflows/code-coverage.yml b/.github/workflows/code-coverage.yml
index 26ea24f02..f25bb30d9 100644
--- a/.github/workflows/code-coverage.yml
+++ b/.github/workflows/code-coverage.yml
@@ -16,7 +16,7 @@ jobs:
     runs-on: ubuntu-latest
 
     steps:
-    - uses: actions/checkout@v2
+    - uses: actions/checkout@v4
 
     - name: Set up environment
       uses: conda-incubator/setup-miniconda@v3
@@ -26,7 +26,7 @@ jobs:
           activate-environment: mrcpp-codecov
           environment-file: .github/mrcpp-codecov.yml
           channel-priority: true
-          python-version: 3.6
+          python-version: 3.11
 
     - name: Configure
       shell: bash -l {0}
diff --git a/cmake/MRCPPConfig.cmake.in b/cmake/MRCPPConfig.cmake.in
index f7a06bd3d..adccf088b 100644
--- a/cmake/MRCPPConfig.cmake.in
+++ b/cmake/MRCPPConfig.cmake.in
@@ -94,7 +94,7 @@ if(NOT TARGET ${PN}::mrcpp)
   get_filename_component(_fext ${${PN}_LIBRARY} EXT)
   include("${CMAKE_CURRENT_LIST_DIR}/${PN}Targets.cmake")
 
-  find_package(Eigen3 3.3 CONFIG REQUIRED)
+  find_package(Eigen3 CONFIG REQUIRED)
 
   get_target_property(MRCPP_HAS_OMP MRCPP::mrcpp MRCPP_HAS_OMP)
   if(MRCPP_HAS_OMP)
diff --git a/docs/mrcpp_api/mwfunctions.rst b/docs/mrcpp_api/mwfunctions.rst
index 527497f5c..5b2b8c6fc 100644
--- a/docs/mrcpp_api/mwfunctions.rst
+++ b/docs/mrcpp_api/mwfunctions.rst
@@ -165,7 +165,7 @@ Constructing an MRA
 
 An MRA is defined in two steps, first the computational domain is given by a
 ``BoundingBox`` (D is the dimension), e.g. for a total domain of
-:math:`[-32,32]^3` in three dimensions (eight root boxes of size :math:`[32]^3`
+:math:`[-16,16]^3` in three dimensions (eight root boxes of size :math:`[16]^3`
 each):
 
 .. code-block:: cpp
diff --git a/external/upstream/fetch_eigen3.cmake b/external/upstream/fetch_eigen3.cmake
index be8f9e0bc..68a6683c8 100644
--- a/external/upstream/fetch_eigen3.cmake
+++ b/external/upstream/fetch_eigen3.cmake
@@ -1,4 +1,4 @@
-find_package(Eigen3 3.4 CONFIG QUIET
+find_package(Eigen3 CONFIG QUIET
   NO_CMAKE_PATH
   NO_CMAKE_PACKAGE_REGISTRY
   )
diff --git a/src/core/CrossCorrelation.cpp b/src/core/CrossCorrelation.cpp
index 8b04cf18d..788491efd 100644
--- a/src/core/CrossCorrelation.cpp
+++ b/src/core/CrossCorrelation.cpp
@@ -106,11 +106,11 @@ void CrossCorrelation::readCCCBin() {
     int K = this->order + 1;
     this->Left = MatrixXd::Zero(K * K, 2 * K);
     this->Right = MatrixXd::Zero(K * K, 2 * K);
-    double dL[2 * K];
-    double dR[2 * K];
+    std::vector<double> dL(2 * K);
+    std::vector<double> dR(2 * K);
     for (int i = 0; i < K * K; i++) {
-        L_fis.read((char *)dL, sizeof(double) * 2 * K);
-        R_fis.read((char *)dR, sizeof(double) * 2 * K);
+        L_fis.read(reinterpret_cast<char *>(dL.data()), sizeof(double) * 2 * K);
+        R_fis.read(reinterpret_cast<char *>(dR.data()), sizeof(double) * 2 * K);
         for (int j = 0; j < 2 * K; j++) {
             if (std::abs(dL[j]) < MachinePrec) dL[j] = 0.0;
             if (std::abs(dR[j]) < MachinePrec) dR[j] = 0.0;
diff --git a/src/core/MWFilter.cpp b/src/core/MWFilter.cpp
index 8ada4f2ca..19ed535cd 100644
--- a/src/core/MWFilter.cpp
+++ b/src/core/MWFilter.cpp
@@ -206,14 +206,14 @@ void MWFilter::generateBlocks() {
 
     int K = this->order + 1;
 
-    double dH[K];
-    double dG[K];
+    std::vector<double> dH(K);
+    std::vector<double> dG(K);
     /* read H0 and G0 from disk */
     this->G0 = Eigen::MatrixXd::Zero(K, K);
     this->H0 = Eigen::MatrixXd::Zero(K, K);
     for (int i = 0; i < K; i++) {
-        H_fis.read((char *)dH, sizeof(double) * K);
-        G_fis.read((char *)dG, sizeof(double) * K);
+        H_fis.read(reinterpret_cast<char *>(dH.data()), sizeof(double) * K);
+        G_fis.read(reinterpret_cast<char *>(dG.data()), sizeof(double) * K);
         for (int j = 0; j < K; j++) {
             this->G0(i, j) = dG[j]; // G0
             this->H0(i, j) = dH[j]; // H0
diff --git a/src/functions/GaussExp.cpp b/src/functions/GaussExp.cpp
index a57fe6708..21432686d 100644
--- a/src/functions/GaussExp.cpp
+++ b/src/functions/GaussExp.cpp
@@ -41,7 +41,7 @@ namespace mrcpp {
 
 template <int D> double GaussExp<D>::defaultScreening = 10.0;
 
-template <int D> GaussExp<D>::GaussExp(int nTerms, double prec) {
+template <int D> GaussExp<D>::GaussExp(int nTerms) {
     for (int i = 0; i < nTerms; i++) { this->funcs.push_back(nullptr); }
 }
 
diff --git a/src/functions/GaussExp.h b/src/functions/GaussExp.h
index a4315e381..382213214 100644
--- a/src/functions/GaussExp.h
+++ b/src/functions/GaussExp.h
@@ -35,8 +35,6 @@
 
 namespace mrcpp {
 
-#define GAUSS_EXP_PREC 1.e-10
-
 /** @class GaussExp
  *
  * @brief Gaussian expansion in D dimensions
@@ -53,7 +51,7 @@ namespace mrcpp {
 
 template <int D> class GaussExp : public RepresentableFunction<D, double> {
 public:
-    GaussExp(int nTerms = 0, double prec = GAUSS_EXP_PREC);
+    GaussExp(int nTerms = 0);
     GaussExp(const GaussExp<D> &gExp);
     GaussExp &operator=(const GaussExp<D> &gExp);
     ~GaussExp() override;
diff --git a/src/functions/Gaussian.cpp b/src/functions/Gaussian.cpp
index 6dbfa7c5b..0b8e3214b 100644
--- a/src/functions/Gaussian.cpp
+++ b/src/functions/Gaussian.cpp
@@ -194,7 +194,6 @@ template <int D> GaussExp<D> Gaussian<D>::periodify(const std::array<double, D>
         auto needed_cells_vec = std::vector<int>();
         for (auto i = 0; i < D; i++) {
             auto upper_bound = pos[i] + x_std;
-            auto lower_bound = pos[i] - x_std;
             // number of cells upp and down relative to the center of the Gaussian
             needed_cells_vec.push_back(std::ceil(upper_bound / period[i]));
         }
diff --git a/src/treebuilders/ABGVCalculator.cpp b/src/treebuilders/ABGVCalculator.cpp
index 10252f46a..376068b5d 100644
--- a/src/treebuilders/ABGVCalculator.cpp
+++ b/src/treebuilders/ABGVCalculator.cpp
@@ -48,7 +48,7 @@ ABGVCalculator::ABGVCalculator(const ScalingBasis &basis, double a, double b)
 
 void ABGVCalculator::calcValueVectors(const ScalingBasis &basis) {
     int kp1 = basis.getQuadratureOrder();
-    double sqrtCoef[kp1];
+    std::vector<double> sqrtCoef(kp1);
     for (int i = 0; i < kp1; i++) { sqrtCoef[i] = std::sqrt(2.0 * i + 1.0); }
 
     switch (basis.getScalingType()) {
@@ -74,7 +74,7 @@ void ABGVCalculator::calcValueVectors(const ScalingBasis &basis) {
 
 void ABGVCalculator::calcKMatrix(const ScalingBasis &basis) {
     int kp1 = basis.getQuadratureOrder();
-    double sqrtCoef[kp1];
+    std::vector<double> sqrtCoef(kp1);
     for (int i = 0; i < kp1; i++) { sqrtCoef[i] = std::sqrt(2.0 * i + 1.0); }
     getQuadratureCache(qCache);
     const VectorXd &roots = qCache.getRoots(kp1);
@@ -105,7 +105,6 @@ void ABGVCalculator::calcNode(MWNode<2> &node) {
     node.zeroCoefs();
 
     const auto &idx = node.getNodeIndex();
-    int l = idx[1] - idx[0];
     int np1 = idx.getScale() + 1;
     int kp1 = node.getKp1();
     int kp1_d = node.getKp1_d();
diff --git a/src/treebuilders/ConvolutionCalculator.cpp b/src/treebuilders/ConvolutionCalculator.cpp
index 497fe0dd8..8cfcac1e6 100644
--- a/src/treebuilders/ConvolutionCalculator.cpp
+++ b/src/treebuilders/ConvolutionCalculator.cpp
@@ -143,7 +143,6 @@ template <int D, typename T> MWNodeVector<D, T> *ConvolutionCalculator<D, T>::ma
     auto *band = new MWNodeVector<D, T>;
 
     int o_depth = gNode.getScale() - this->oper->getOperatorRoot();
-    int g_depth = gNode.getDepth();
     int width = this->oper->getMaxBandWidth(o_depth);
 
     bool periodic = gNode.getMWTree().isPeriodic();
@@ -228,8 +227,8 @@ template <int D, typename T> void ConvolutionCalculator<D, T>::calcNode(MWNode<D
 
     int o_depth = gNode.getScale() - this->oper->getOperatorRoot();
     if (manipulateOperator and this->oper->getOperatorRoot() < 0) o_depth = gNode.getDepth();
-    T tmpCoefs[gNode.getNCoefs()];
-    OperatorState<D, T> os(gNode, tmpCoefs);
+    std::vector<T> tmpCoefs(gNode.getNCoefs());
+    OperatorState<D, T> os(gNode, tmpCoefs.data());
     this->operStat.incrementGNodeCounters(gNode);
 
     // Get all nodes in f within the bandwith of O in g
diff --git a/src/treebuilders/DerivativeCalculator.cpp b/src/treebuilders/DerivativeCalculator.cpp
index b298d1b6e..f4936c01d 100644
--- a/src/treebuilders/DerivativeCalculator.cpp
+++ b/src/treebuilders/DerivativeCalculator.cpp
@@ -90,8 +90,8 @@ template <int D, typename T> void DerivativeCalculator<D, T>::calcNode(MWNode<D,
     // if (this->oper->getMaxBandWidth() > 1) MSG_ABORT("Only implemented for zero bw");
     outNode.zeroCoefs();
     int nComp = (1 << D);
-    T tmpCoefs[outNode.getNCoefs()];
-    OperatorState<D, T> os(outNode, tmpCoefs);
+    std::vector<T> tmpCoefs(outNode.getNCoefs());
+    OperatorState<D, T> os(outNode, tmpCoefs.data());
 
     os.setFNode(inpNode);
     os.setFIndex(inpNode.nodeIndex);
@@ -116,8 +116,8 @@ template <int D, typename T> void DerivativeCalculator<D, T>::calcNode(MWNode<D,
     gNode.zeroCoefs();
 
     int nComp = (1 << D);
-    T tmpCoefs[gNode.getNCoefs()];
-    OperatorState<D, T> os(gNode, tmpCoefs);
+    std::vector<T> tmpCoefs(gNode.getNCoefs());
+    OperatorState<D, T> os(gNode, tmpCoefs.data());
     this->operStat.incrementGNodeCounters(gNode);
 
     // Get all nodes in f within the bandwith of O in g
@@ -182,11 +182,8 @@ template <int D, typename T> void DerivativeCalculator<D, T>::applyOperator_bw0(
     // cout<<" applyOperator "<<endl;
     MWNode<D, T> &gNode = *os.gNode;
     MWNode<D, T> &fNode = *os.fNode;
-    const NodeIndex<D> &fIdx = *os.fIdx;
-    const NodeIndex<D> &gIdx = gNode.getNodeIndex();
     int depth = gNode.getDepth();
 
-    double oNorm = 1.0;
     double **oData = os.getOperData();
 
     for (int d = 0; d < D; d++) {
@@ -218,7 +215,6 @@ template <int D, typename T> void DerivativeCalculator<D, T>::applyOperator(Oper
     const NodeIndex<D> &gIdx = gNode.getNodeIndex();
     int depth = gNode.getDepth();
 
-    double oNorm = 1.0;
     double **oData = os.getOperData();
 
     for (int d = 0; d < D; d++) {
@@ -236,8 +232,6 @@ template <int D, typename T> void DerivativeCalculator<D, T>::applyOperator(Oper
 
         const OperatorNode &oNode = oTree.getNode(depth, oTransl);
         int oIdx = os.getOperIndex(d);
-        double ocn = oNode.getComponentNorm(oIdx);
-        oNorm *= ocn;
         if (this->applyDir == d) {
             oData[d] = const_cast<double *>(oNode.getCoefs()) + oIdx * os.kp1_2;
         } else {
diff --git a/src/treebuilders/ProjectionCalculator.cpp b/src/treebuilders/ProjectionCalculator.cpp
index 931733232..15ba014da 100644
--- a/src/treebuilders/ProjectionCalculator.cpp
+++ b/src/treebuilders/ProjectionCalculator.cpp
@@ -25,6 +25,7 @@
 
 #include "ProjectionCalculator.h"
 #include "trees/MWNode.h"
+#include <cassert>
 
 using Eigen::MatrixXd;
 
diff --git a/src/treebuilders/apply.cpp b/src/treebuilders/apply.cpp
index 77d8674be..bfcbf164e 100644
--- a/src/treebuilders/apply.cpp
+++ b/src/treebuilders/apply.cpp
@@ -118,6 +118,7 @@ template <int D, typename T> void apply(double prec, FunctionTree<D, T> &out, Co
  */
 template <int D> void apply(double prec, CompFunction<D> &out, ConvolutionOperator<D> &oper, const CompFunction<D> &inp, const ComplexDouble (*metric)[4], int maxIter, bool absPrec) {
 
+    out = inp.paramCopy(true);
     for (int icomp = 0; icomp < inp.Ncomp(); icomp++) {
         for (int ocomp = 0; ocomp < 4; ocomp++) {
             if (std::norm(metric[icomp][ocomp]) > MachinePrec) {
@@ -252,6 +253,7 @@ template <int D, typename T> void apply(double prec, FunctionTree<D, T> &out, Co
 template <int D, typename T>
 void apply(double prec, CompFunction<D> &out, ConvolutionOperator<D> &oper, CompFunction<D> &inp, FunctionTreeVector<D, T> *precTrees, const ComplexDouble (*metric)[4], int maxIter, bool absPrec) {
 
+    out = inp.paramCopy(true);
     for (int icomp = 0; icomp < inp.Ncomp(); icomp++) {
         for (int ocomp = 0; ocomp < 4; ocomp++) {
             if (std::norm(metric[icomp][ocomp]) > MachinePrec) {
@@ -295,6 +297,7 @@ template <int D, typename T> void apply_far_field(double prec, FunctionTree<D, T
 
 template <int D> void apply_far_field(double prec, CompFunction<D> &out, ConvolutionOperator<D> &oper, CompFunction<D> &inp, const ComplexDouble (*metric)[4], int maxIter, bool absPrec) {
 
+    out = inp.paramCopy(true);
     for (int icomp = 0; icomp < 4; icomp++) {
         if (inp.Comp[icomp] != nullptr) {
             for (int ocomp = 0; ocomp < 4; ocomp++) {
@@ -411,6 +414,7 @@ template <int D, typename T> void apply(FunctionTree<D, T> &out, DerivativeOpera
 template <int D> void apply(CompFunction<D> &out, DerivativeOperator<D> &oper, CompFunction<D> &inp, int dir, const ComplexDouble (*metric)[4]) {
     // TODO: sums and not only each components independently, when concrete examples with non diagonal metric are tested
 
+    out = inp.paramCopy(true); // note that this will copy the factor of inp (inp.func_ptr->data.c1)
     for (int icomp = 0; icomp < inp.Ncomp(); icomp++) {
         for (int ocomp = 0; ocomp < 4; ocomp++) {
             if (std::norm(metric[icomp][ocomp]) > MachinePrec) {
@@ -469,7 +473,7 @@ std::vector<CompFunction<3> *> gradient(DerivativeOperator<3> &oper, CompFunctio
         for (int icomp = 0; icomp < inp.Ncomp(); icomp++) {
             for (int ocomp = 0; ocomp < 4; ocomp++) {
                 if (std::norm(metric[icomp][ocomp]) > MachinePrec) {
-                    grad_d->func_ptr->Ncomp = ocomp;
+                    grad_d->func_ptr->Ncomp = ocomp + 1;
                     if (inp.isreal()) {
                         grad_d->func_ptr->isreal = 1;
                         grad_d->func_ptr->iscomplex = 0;
@@ -609,9 +613,6 @@ template void divergence<3, ComplexDouble>(FunctionTree<3, ComplexDouble> &out,
 template void divergence<1, ComplexDouble>(FunctionTree<1, ComplexDouble> &out, DerivativeOperator<1> &oper, std::vector<FunctionTree<1, ComplexDouble> *> &inp);
 template void divergence<2, ComplexDouble>(FunctionTree<2, ComplexDouble> &out, DerivativeOperator<2> &oper, std::vector<FunctionTree<2, ComplexDouble> *> &inp);
 template void divergence<3, ComplexDouble>(FunctionTree<3, ComplexDouble> &out, DerivativeOperator<3> &oper, std::vector<FunctionTree<3, ComplexDouble> *> &inp);
-template FunctionTreeVector<1, ComplexDouble> gradient<1>(DerivativeOperator<1> &oper, FunctionTree<1, ComplexDouble> &inp);
-template FunctionTreeVector<2, ComplexDouble> gradient<2>(DerivativeOperator<2> &oper, FunctionTree<2, ComplexDouble> &inp);
-template FunctionTreeVector<3, ComplexDouble> gradient<3>(DerivativeOperator<3> &oper, FunctionTree<3, ComplexDouble> &inp);
 
 template void apply(CompFunction<3> &out, DerivativeOperator<3> &oper, CompFunction<3> &inp, int dir = -1, const ComplexDouble (*metric)[4]);
 
diff --git a/src/treebuilders/apply.h b/src/treebuilders/apply.h
index e50b3b870..fa5a43661 100644
--- a/src/treebuilders/apply.h
+++ b/src/treebuilders/apply.h
@@ -53,8 +53,7 @@ template <int D, typename T> void divergence(CompFunction<D> &out, DerivativeOpe
 template <int D, typename T> void divergence(FunctionTree<D, T> &out, DerivativeOperator<D> &oper, std::vector<FunctionTree<D, T> *> &inp);
 template <int D, typename T> void divergence(CompFunction<D> &out, DerivativeOperator<D> &oper, std::vector<FunctionTree<D, T> *> *inp, const ComplexDouble (*metric)[4]  = defaultMetric);
 template <int D, typename T> FunctionTreeVector<D, T> gradient(DerivativeOperator<D> &oper, FunctionTree<D, T> &inp);
-// template <int D>
-std::vector<CompFunction<3>*> gradient(DerivativeOperator<3> &oper, CompFunction<3> &inp, ComplexDouble  (*metric)[4] = nullptr);
+std::vector<CompFunction<3>*> gradient(DerivativeOperator<3> &oper, CompFunction<3> &inp, const ComplexDouble  (*metric)[4] = defaultMetric);
 // clang-format on
 
 } // namespace mrcpp
diff --git a/src/treebuilders/grid.cpp b/src/treebuilders/grid.cpp
index 0e7fb968b..eaaf3696b 100644
--- a/src/treebuilders/grid.cpp
+++ b/src/treebuilders/grid.cpp
@@ -113,7 +113,6 @@ template <int D> void build_grid(FunctionTree<D> &out, const GaussExp<D> &inp, i
             builder.build(out, calculator, adaptor, maxIter);
         }
     } else {
-        auto period = out.getMRA().getWorldBox().getScalingFactors();
         for (auto i = 0; i < inp.size(); i++) {
             auto *gauss = inp.getFunc(i).copy();
             build_grid(out, *gauss, maxIter);
diff --git a/src/treebuilders/multiply.cpp b/src/treebuilders/multiply.cpp
index 4e046126e..46587ffaa 100644
--- a/src/treebuilders/multiply.cpp
+++ b/src/treebuilders/multiply.cpp
@@ -334,8 +334,8 @@ template <int D, typename T> double node_norm_dot(FunctionTree<D, T> &bra, Funct
 
     double result = 0.0;
     int ncoef = bra.getKp1_d() * bra.getTDim();
-    T valA[ncoef];
-    T valB[ncoef];
+    std::vector<T> valA(ncoef);
+    std::vector<T> valB(ncoef);
     int nNodes = bra.getNEndNodes();
 
     for (int n = 0; n < nNodes; n++) {
@@ -345,8 +345,8 @@ template <int D, typename T> double node_norm_dot(FunctionTree<D, T> &bra, Funct
             // convert to interpolating coef, take abs, convert back
             FunctionNode<D, T> *mwNode = static_cast<FunctionNode<D, T> *>(ket.findNode(idx));
             if (mwNode == nullptr) MSG_ABORT("Trees must have same grid");
-            node.getAbsCoefs(valA);
-            mwNode->getAbsCoefs(valB);
+            node.getAbsCoefs(valA.data());
+            mwNode->getAbsCoefs(valB.data());
             for (int i = 0; i < ncoef; i++) result += std::norm(valA[i] * valB[i]);
         } else {
             // approximate by product of node norms
diff --git a/src/trees/FunctionNode.cpp b/src/trees/FunctionNode.cpp
index ff23fb394..102e5bd26 100644
--- a/src/trees/FunctionNode.cpp
+++ b/src/trees/FunctionNode.cpp
@@ -130,7 +130,7 @@ template <int D, typename T> T FunctionNode<D, T>::integrateInterpolating() cons
     int qOrder = this->getKp1();
     getQuadratureCache(qc);
     const VectorXd &weights = qc.getWeights(qOrder);
-    double sqWeights[qOrder];
+    std::vector<double> sqWeights(qOrder);
     for (int i = 0; i < qOrder; i++) sqWeights[i] = std::sqrt(weights[i]);
 
     int kp1_p[D];
@@ -170,7 +170,7 @@ template <int D, typename T> T FunctionNode<D, T>::integrateValues() const {
     this->getCoefs(coefs);
     int ncoefs = coefs.size();
     int ncoefChild = ncoefs / (1 << D);
-    T cc[ncoefChild];
+    std::vector<T> cc(ncoefChild);
     // factorize out the children
     for (int i = 0; i < ncoefChild; i++) cc[i] = coefs[i];
     for (int j = 1; j < (1 << D); j++)
diff --git a/src/trees/FunctionTree.cpp b/src/trees/FunctionTree.cpp
index a41581692..995578436 100644
--- a/src/trees/FunctionTree.cpp
+++ b/src/trees/FunctionTree.cpp
@@ -290,9 +290,6 @@ template <int D, typename T> void FunctionTree<D, T>::loadTreeTXT(const std::str
  * @param[in] file: File name
  */
 template <int D, typename T> void FunctionTree<D, T>::saveTreeTXT(const std::string &fname) {
-    int nRoots = this->getRootBox().size();
-    MWNode<D, T> **roots = this->getRootBox().getNodes();
-
     std::ofstream out(fname);
     out << std::setprecision(14);
     out << D << std::endl;
@@ -306,6 +303,7 @@ template <int D, typename T> void FunctionTree<D, T>::saveTreeTXT(const std::str
     int ncoefs = 1;
     for (int d = 0; d < D; d++) ncoefs *= kp1;
     int Tdim = std::pow(2, D);
+    std::vector<T> values(ncoefs * Tdim);
 
     int nout = this->endNodeTable.size();
     out << Tdim * nout << std::endl; // could output only scaling coeff?
@@ -331,7 +329,9 @@ template <int D, typename T> void FunctionTree<D, T>::saveTreeTXT(const std::str
         std::array<int, D> l;
         NodeIndex<D> idx = this->endNodeTable[count]->getNodeIndex();
         MWNode<D, T> *node = &(this->getNode(idx, false));
-        T *values = node->getCoefs();
+	T *coefs = node->getCoefs();
+	for (int i = 0; i < ncoefs * Tdim; i++) values[i] = coefs[i];
+	    node->attachCoefs(values.data());
         int n = idx.getScale();
         node->mwTransform(Reconstruction);
         node->cvTransform(Forward);
@@ -350,7 +350,8 @@ template <int D, typename T> void FunctionTree<D, T>::saveTreeTXT(const std::str
             for (int i = 0; i < ncoefs; i++) out << values[cix * ncoefs + mapMRC[i]] << " ";
             out << std::endl;
         }
-    }
+	node->attachCoefs(coefs); // put back original coeff
+   }
     out.close();
 }
 /** @brief Write the tree structure to disk, for later use
@@ -358,9 +359,9 @@ template <int D, typename T> void FunctionTree<D, T>::saveTreeTXT(const std::str
  */
 template <int D, typename T> void FunctionTree<D, T>::saveTree(const std::string &file) {
     Timer t1;
+
     this->deleteGenerated();
     auto &allocator = this->getNodeAllocator();
-
     std::stringstream fname;
     fname << file << ".tree";
     std::fstream f;
@@ -370,14 +371,12 @@ template <int D, typename T> void FunctionTree<D, T>::saveTree(const std::string
     // Write size of tree
     int nChunks = allocator.getNChunksUsed();
     f.write((char *)&nChunks, sizeof(int));
-
     // Write tree data, chunk by chunk
     for (int iChunk = 0; iChunk < nChunks; iChunk++) {
-        f.write((char *)allocator.getNodeChunk(iChunk), allocator.getNodeChunkSize());
-        f.write((char *)allocator.getCoefChunk(iChunk), allocator.getCoefChunkSize());
+       f.write((char *)allocator.getNodeChunk(iChunk), allocator.getNodeChunkSize());
+       f.write((char *)allocator.getCoefChunk(iChunk), allocator.getCoefChunkSize());
     }
     f.close();
-    this->saveTreeTXT("MRC.dat");
     print::time(10, "Time write", t1);
 }
 
@@ -441,7 +440,6 @@ template <> double FunctionTree<3, double>::integrateEndNodes(RepresentableFunct
     // traverse tree, and treat end nodes only
     std::vector<FunctionNode<3> *> stack; // node from this
     for (int i = 0; i < this->getRootBox().size(); i++) stack.push_back(&(this->getRootFuncNode(i)));
-    int basis = getMRA().getScalingBasis().getScalingType();
     double result = 0.0;
     while (stack.size() > 0) {
         FunctionNode<3> *Node = stack.back();
@@ -454,9 +452,9 @@ template <> double FunctionTree<3, double>::integrateEndNodes(RepresentableFunct
             double *coefs = Node->getCoefs(); // save position of coeff, but do not use them!
             // The data in fmat is not organized so that two consecutive points are stored after each other in memory, so needs to copy before mwtransform, cannot use memory adress directly.
             int nc = fmat.cols();
-            double cc[nc];
+            std::vector<double> cc(nc);
             for (int i = 0; i < nc; i++) cc[i] = fmat(0, i);
-            Node->attachCoefs(cc);
+            Node->attachCoefs(cc.data());
             result += Node->integrateValues();
             Node->attachCoefs(coefs); // put back original coeff
         }
@@ -878,8 +876,6 @@ void FunctionTree<D, T>::makeCoeffVector(std::vector<T *> &coefs,
     indices.clear();
     parent_indices.clear();
     max_index = 0;
-    int sizecoeff = (1 << refTree.getDim()) * refTree.getKp1_d();
-    int sizecoeffW = ((1 << refTree.getDim()) - 1) * refTree.getKp1_d();
     std::vector<MWNode<D, double> *> refstack; // nodes from refTree
     std::vector<MWNode<D, T> *> thisstack;     // nodes from this Tree
     for (int rIdx = 0; rIdx < this->getRootBox().size(); rIdx++) {
@@ -1094,7 +1090,6 @@ template <> int FunctionTree<3, double>::saveNodesAndRmCoeff() {
         for (int rIdx = 0; rIdx < this->getRootBox().size(); rIdx++) { stack.push_back(this->getRootBox().getNodes()[rIdx]); }
         while (stack.size() > stack_p) {
             MWNode<3, double> *Node = stack[stack_p++];
-            int id = 0;
             NodesCoeff->put_data(Node->getNodeIndex(), sizecoeff, Node->getCoefs());
             for (int i = 0; i < Node->getNChildren(); i++) { stack.push_back(Node->children[i]); }
         }
@@ -1117,7 +1112,6 @@ template <> int FunctionTree<3, ComplexDouble>::saveNodesAndRmCoeff() {
         for (int rIdx = 0; rIdx < this->getRootBox().size(); rIdx++) { stack.push_back(this->getRootBox().getNodes()[rIdx]); }
         while (stack.size() > stack_p) {
             MWNode<3, ComplexDouble> *Node = stack[stack_p++];
-            int id = 0;
             NodesCoeff->put_data(Node->getNodeIndex(), sizecoeff, Node->getCoefs());
             for (int i = 0; i < Node->getNChildren(); i++) { stack.push_back(Node->children[i]); }
         }
@@ -1142,19 +1136,22 @@ template <int D, typename T> void FunctionTree<D, T>::deep_copy(FunctionTree<D,
  */
 template <int D, typename T> FunctionTree<D, double> *FunctionTree<D, T>::Real() {
     FunctionTree<D, double> *out = new FunctionTree<D, double>(this->getMRA(), this->getName());
-#pragma omp parallel num_threads(mrcpp_get_num_threads())
+    out->setZero();
+    // TODO: why does the omp parallelism not always work here?
+    //#pragma omp parallel num_threads(mrcpp_get_num_threads())
     {
         int nNodes = this->getNEndNodes();
-#pragma omp for schedule(guided)
+	//#pragma omp for schedule(guided)
         for (int n = 0; n < nNodes; n++) {
             MWNode<D, T> &inp_node = *this->endNodeTable[n];
-            MWNode<D, double> out_node = out->getNode(out_node.getNodeIndex()); // Full copy
+            MWNode<D, double> &out_node = out->getNode(inp_node.getNodeIndex(), true);
             double *out_coefs = out_node.getCoefs();
             const T *inp_coefs = inp_node.getCoefs();
             for (int i = 0; i < inp_node.getNCoefs(); i++) { out_coefs[i] = std::real(inp_coefs[i]); }
             out_node.calcNorms();
         }
     }
+    out->resetEndNodeTable();
     out->mwTransform(BottomUp);
     out->calcSquareNorm();
     return out;
@@ -1164,13 +1161,14 @@ template <int D, typename T> FunctionTree<D, double> *FunctionTree<D, T>::Real()
  */
 template <int D, typename T> FunctionTree<D, double> *FunctionTree<D, T>::Imag() {
     FunctionTree<D, double> *out = new FunctionTree<D, double>(this->getMRA(), this->getName());
-#pragma omp parallel num_threads(mrcpp_get_num_threads())
+    out->setZero();
+    //#pragma omp parallel num_threads(mrcpp_get_num_threads())
     {
         int nNodes = this->getNEndNodes();
-#pragma omp for schedule(guided)
+	//#pragma omp for schedule(guided)
         for (int n = 0; n < nNodes; n++) {
             MWNode<D, T> &inp_node = *this->endNodeTable[n];
-            MWNode<D, double> out_node = out->getNode(out_node.getNodeIndex()); // Full copy
+            MWNode<D, double> &out_node = out->getNode(inp_node.getNodeIndex(), true);
             double *out_coefs = out_node.getCoefs();
             const T *inp_coefs = inp_node.getCoefs();
             for (int i = 0; i < inp_node.getNCoefs(); i++) { out_coefs[i] = std::imag(inp_coefs[i]); }
@@ -1189,7 +1187,6 @@ template <int D, typename T> FunctionTree<D, double> *FunctionTree<D, T>::Imag()
 
 template <> void FunctionTree<3, double>::CopyTreeToComplex(FunctionTree<3, ComplexDouble> *&outTree) {
     delete outTree;
-    double ref = 0.0;
     outTree = new FunctionTree<3, ComplexDouble>(this->getMRA());
     std::vector<MWNode<3, double> *> instack;         // node from this
     std::vector<MWNode<3, ComplexDouble> *> outstack; // node from outTree
@@ -1198,7 +1195,6 @@ template <> void FunctionTree<3, double>::CopyTreeToComplex(FunctionTree<3, Comp
         instack.push_back(this->getRootBox().getNodes()[rIdx]);
         outstack.push_back(outTree->getRootBox().getNodes()[rIdx]);
     }
-    int nNodes = std::min(this->getNNodes(), this->getNodeAllocator().getMaxNodesPerChunk());
     int ncoefs = this->getNodeAllocator().getNCoefs();
     while (instack.size() > 0) {
         // inNode and outNode are the same node in space, but on different trees
@@ -1229,7 +1225,6 @@ template <> void FunctionTree<3, double>::CopyTreeToComplex(FunctionTree<3, Comp
 
 template <> void FunctionTree<2, double>::CopyTreeToComplex(FunctionTree<2, ComplexDouble> *&outTree) {
     delete outTree;
-    double ref = 0.0;
     outTree = new FunctionTree<2, ComplexDouble>(this->getMRA());
     std::vector<MWNode<2, double> *> instack;         // node from this
     std::vector<MWNode<2, ComplexDouble> *> outstack; // node from outTree
@@ -1238,7 +1233,6 @@ template <> void FunctionTree<2, double>::CopyTreeToComplex(FunctionTree<2, Comp
         instack.push_back(this->getRootBox().getNodes()[rIdx]);
         outstack.push_back(outTree->getRootBox().getNodes()[rIdx]);
     }
-    int nNodes = std::min(this->getNNodes(), this->getNodeAllocator().getMaxNodesPerChunk());
     int ncoefs = this->getNodeAllocator().getNCoefs();
     while (instack.size() > 0) {
         // inNode and outNode are the same node in space, but on different trees
@@ -1269,7 +1263,6 @@ template <> void FunctionTree<2, double>::CopyTreeToComplex(FunctionTree<2, Comp
 
 template <> void FunctionTree<1, double>::CopyTreeToComplex(FunctionTree<1, ComplexDouble> *&outTree) {
     delete outTree;
-    double ref = 0.0;
     outTree = new FunctionTree<1, ComplexDouble>(this->getMRA());
     std::vector<MWNode<1, double> *> instack;         // node from this
     std::vector<MWNode<1, ComplexDouble> *> outstack; // node from outTree
@@ -1278,7 +1271,6 @@ template <> void FunctionTree<1, double>::CopyTreeToComplex(FunctionTree<1, Comp
         instack.push_back(this->getRootBox().getNodes()[rIdx]);
         outstack.push_back(outTree->getRootBox().getNodes()[rIdx]);
     }
-    int nNodes = std::min(this->getNNodes(), this->getNodeAllocator().getMaxNodesPerChunk());
     int ncoefs = this->getNodeAllocator().getNCoefs();
     while (instack.size() > 0) {
         // inNode and outNode are the same node in space, but on different trees
@@ -1310,7 +1302,6 @@ template <> void FunctionTree<1, double>::CopyTreeToComplex(FunctionTree<1, Comp
 // for testing
 template <> void FunctionTree<3, double>::CopyTreeToReal(FunctionTree<3, double> *&outTree) {
     delete outTree;
-    double ref = 0.0;
     // FunctionTree<3, double>* inTree = this;
     outTree = new FunctionTree<3, double>(this->getMRA());
     std::vector<MWNode<3, double> *> instack;  // node from this
@@ -1320,7 +1311,6 @@ template <> void FunctionTree<3, double>::CopyTreeToReal(FunctionTree<3, double>
         instack.push_back(this->getRootBox().getNodes()[rIdx]);
         outstack.push_back(outTree->getRootBox().getNodes()[rIdx]);
     }
-    int nNodes = std::min(this->getNNodes(), this->getNodeAllocator().getMaxNodesPerChunk());
     int ncoefs = this->getNodeAllocator().getNCoefs();
     while (instack.size() > 0) {
         // inNode and outNode are the same node in space, but on different trees
diff --git a/src/trees/HilbertPath.cpp b/src/trees/HilbertPath.cpp
index b052064a8..ac2529852 100644
--- a/src/trees/HilbertPath.cpp
+++ b/src/trees/HilbertPath.cpp
@@ -114,6 +114,21 @@ const int HilbertPath<3>::hTable[12][8] = {
     {4,5,7,6,3,2,0,1}
 };
 
+template <int D>
+short int HilbertPath<D>::getChildPath(int hIdx) const {
+    return pTable[this->path][hIdx];
+}
+
+template <int D>
+int HilbertPath<D>::getZIndex(int hIdx) const {
+    return zTable[this->path][hIdx];
+}
+
+template <int D>
+int HilbertPath<D>::getHIndex(int zIdx) const {
+    return hTable[this->path][zIdx];
+}
+
 // clang-format on
 
 template class HilbertPath<1>;
diff --git a/src/trees/HilbertPath.h b/src/trees/HilbertPath.h
index 519e7ac73..981c5354f 100644
--- a/src/trees/HilbertPath.h
+++ b/src/trees/HilbertPath.h
@@ -42,10 +42,10 @@ template <int D> class HilbertPath final {
     }
 
     short int getPath() const { return this->path; }
-    short int getChildPath(int hIdx) const { return this->pTable[this->path][hIdx]; }
+    short int getChildPath(int hIdx) const;
 
-    int getZIndex(int hIdx) const { return this->zTable[this->path][hIdx]; }
-    int getHIndex(int zIdx) const { return this->hTable[this->path][zIdx]; }
+    int getZIndex(int hIdx) const;
+    int getHIndex(int zIdx) const;
 
 private:
     short int path{0};
@@ -54,4 +54,4 @@ template <int D> class HilbertPath final {
     static const int hTable[][8];
 };
 
-} // namespace mrcpp
+} // namespace mrcpp
\ No newline at end of file
diff --git a/src/trees/MWNode.cpp b/src/trees/MWNode.cpp
index 2d521b468..9f40bed23 100644
--- a/src/trees/MWNode.cpp
+++ b/src/trees/MWNode.cpp
@@ -346,13 +346,9 @@ template <int D, typename T> void MWNode<D, T>::giveChildrenCoefs(bool overwrite
  * copied/summed in the correct child node.
  */
 template <int D, typename T> void MWNode<D, T>::giveChildCoefs(int cIdx, bool overwrite) {
-
     MWNode<D, T> node_i = *this;
-
     node_i.mwTransform(Reconstruction);
-
     int kp1_d = this->getKp1_d();
-    int nChildren = this->getTDim();
 
     if (this->children[cIdx] == nullptr) MSG_ABORT("Child does not exist!");
     MWNode<D, T> &child = getMWChild(cIdx);
@@ -457,8 +453,8 @@ template <int D, typename T> void MWNode<D, T>::cvTransform(int operation, bool
 
     auto sb = this->getMWTree().getMRA().getScalingBasis();
     const MatrixXd &S = sb.getCVMap(operation);
-    T o_vec[nCoefs];
-    T *out_vec = o_vec;
+    std::vector<T> o_vec(nCoefs);
+    T *out_vec = o_vec.data();
     T *in_vec = this->coefs;
 
     int nChildren = this->getTDim();
@@ -566,8 +562,8 @@ template <int D, typename T> void MWNode<D, T>::mwTransform(int operation) {
     const MWFilter &filter = getMWTree().getMRA().getFilter();
     double overwrite = 0.0;
 
-    T o_vec[nCoefs];
-    T *out_vec = o_vec;
+    std::vector<T> o_vec(nCoefs);
+    T *out_vec = o_vec.data();
     T *in_vec = this->coefs;
 
     for (int i = 0; i < D; i++) {
@@ -1259,7 +1255,6 @@ template <int D, typename T> std::ostream &MWNode<D, T>::print(std::ostream &o)
  * i.e. *not* same normalization as a squareNorm
  */
 template <int D, typename T> void MWNode<D, T>::setMaxSquareNorm() {
-    auto n = this->getScale();
     this->maxWSquareNorm = calcScaledWSquareNorm();
     this->maxSquareNorm = calcScaledSquareNorm();
 
@@ -1275,7 +1270,6 @@ template <int D, typename T> void MWNode<D, T>::setMaxSquareNorm() {
 /** @brief recursively reset maxSquaredNorm and maxWSquareNorm of parent and descendants to value -1
  */
 template <int D, typename T> void MWNode<D, T>::resetMaxSquareNorm() {
-    auto n = this->getScale();
     this->maxSquareNorm = -1.0;
     this->maxWSquareNorm = -1.0;
     if (not this->isEndNode()) {
diff --git a/src/trees/NodeAllocator.cpp b/src/trees/NodeAllocator.cpp
index 5079459aa..cd3cf9fcc 100644
--- a/src/trees/NodeAllocator.cpp
+++ b/src/trees/NodeAllocator.cpp
@@ -128,10 +128,6 @@ template <int D, typename T> int NodeAllocator<D, T>::alloc(int nNodes, bool coe
     // return value is index of first new node
     auto sIdx = this->topStack;
 
-    // we require that the index for first child is a multiple of 2**D
-    // so that we can find the sibling rank using rank=sIdx%(2**D)
-    if (sIdx % nNodes != 0) MSG_ERROR(" node allocate error");
-
     // fill stack status
     auto &status = this->stackStatus;
     for (int i = sIdx; i < sIdx + nNodes; i++) {
diff --git a/src/trees/OperatorNode.cpp b/src/trees/OperatorNode.cpp
index 37f576eac..6253ba51b 100644
--- a/src/trees/OperatorNode.cpp
+++ b/src/trees/OperatorNode.cpp
@@ -94,10 +94,6 @@ double OperatorNode::calcComponentNorm(int i) const {
  *
  */
 MatrixXd OperatorNode::getComponent(int i) {
-    int depth = getDepth();
-    double prec = getOperTree().getNormPrecision();
-    double thrs = std::max(MachinePrec, prec / (8.0 * (1 << depth)));
-
     VectorXd coef_vec;
     this->getCoefs(coef_vec);
 
diff --git a/src/utils/Bank.cpp b/src/utils/Bank.cpp
index f8c111a53..94152c06e 100644
--- a/src/utils/Bank.cpp
+++ b/src/utils/Bank.cpp
@@ -8,8 +8,10 @@ namespace mrcpp {
 using namespace Eigen;
 using namespace std;
 
-int metadata_block[3]; // can add more metadata in future
-int const size_metadata = 3;
+#ifdef MRCPP_HAS_MPI
+    int metadata_block[3]; // can add more metadata in future
+    int const size_metadata = 3;
+#endif
 
 Bank::~Bank() {
     // delete all data and accounts
@@ -57,7 +59,9 @@ std::map<int, std::map<int, OrbBlock_struct> *> get_orbid2block;   // to get blo
 
 std::map<int, mem_struct *> mem;
 
+#ifdef MRCPP_HAS_MPI
 int const MIN_SCALE = -999; // Smaller than smallest scale
+#endif
 int naccounts = 0;
 
 void Bank::open() {
diff --git a/src/utils/Bank.h b/src/utils/Bank.h
index 69719c530..d0e032c3f 100644
--- a/src/utils/Bank.h
+++ b/src/utils/Bank.h
@@ -77,7 +77,10 @@ class Bank {
     void clear_bank();
     void remove_account(int account); // remove the content and the account
 
+    #ifdef MRCPP_HAS_MPI
     long long totcurrentsize = 0ll;                     // number of kB used by all accounts
+    long long maxsize = 0;                              // max total deposited data size (without containers)
+    #endif
     std::vector<int> accounts;                          // open bank accounts
     std::map<int, std::vector<deposit> *> get_deposits; // gives deposits of an account
     std::map<int, std::map<int, int> *> get_id2ix;
@@ -85,7 +88,6 @@ class Bank {
     std::map<int, std::vector<queue_struct> *> get_queue;            // gives deposits of an account
     std::map<int, std::map<int, std::vector<int>> *> get_readytasks; // used by task manager
     std::map<int, long long> currentsize;                            // total deposited data size (without containers)
-    long long maxsize = 0;                                           // max total deposited data size (without containers)
 };
 
 class BankAccount {
diff --git a/src/utils/CompFunction.cpp b/src/utils/CompFunction.cpp
index 43306a346..250796c25 100644
--- a/src/utils/CompFunction.cpp
+++ b/src/utils/CompFunction.cpp
@@ -122,7 +122,10 @@ template <int D>
  * Returns a copy without defined trees.
  */
 CompFunction<D> CompFunction<D>::paramCopy(bool alloc) const {
-    return CompFunction<D>(func_ptr->data, alloc);
+    CompFunction<D> out(func_ptr->data, alloc);
+    // we do not copy tree sizes:
+    for (int i = 0; i < 4; i++) out.func_ptr->data.Nchunks[i] = 0;
+    return out;
 }
 
 template <int D> void CompFunction<D>::flushMRAData() {
@@ -333,7 +336,6 @@ template <int D> const FunctionTree<D, ComplexDouble> &CompFunction<D>::complex(
 /* for backwards compatibility */
 template <int D> void CompFunction<D>::setReal(FunctionTree<D, double> *tree, int i) {
     func_ptr->isreal = 1;
-    // if (CompD[i] != nullptr) delete CompD[i];
     CompD[i] = tree;
     if (tree != nullptr) {
         func_ptr->Ncomp = std::max(Ncomp(), i + 1);
@@ -344,7 +346,6 @@ template <int D> void CompFunction<D>::setReal(FunctionTree<D, double> *tree, in
 
 template <int D> void CompFunction<D>::setCplx(FunctionTree<D, ComplexDouble> *tree, int i) {
     func_ptr->iscomplex = 1;
-    // if (CompC[i] != nullptr) delete CompC[i];
     CompC[i] = tree;
     if (tree != nullptr) {
         func_ptr->Ncomp = std::max(Ncomp(), i + 1);
@@ -359,36 +360,91 @@ template <int D> void CompFunction<D>::setCplx(FunctionTree<D, ComplexDouble> *t
  *
  */
 template <int D> void CompFunction<D>::add(ComplexDouble c, CompFunction<D> inp) {
-
-    if (Ncomp() < inp.Ncomp()) {
+    if (inp.getSquareNorm() < MachineZero) {
+        // nothing to add
+    } else if(this->getSquareNorm() < MachineZero) {
+        //copy
         func_ptr->data = inp.func_ptr->data;
         alloc(inp.Ncomp(), true);
-    }
-
-    for (int i = 0; i < inp.Ncomp(); i++) {
-        if (inp.isreal() and c.imag() < MachineZero) {
-            CompD[i]->add_inplace(c.real(), *inp.CompD[i]);
-        } else {
-            if (this->isreal()) {
+        for (int i = 0; i < inp.Ncomp(); i++) {
+            if(inp.isreal()) CompD[i]->add_inplace(c.real(), *inp.CompD[i]);
+            if(inp.iscomplex()) CompC[i]->add_inplace(c, *inp.CompC[i]);
+        }
+    } else {
+        // we must check if the result is real or complex
+        if (inp.isreal() and this->isreal() and c.imag() < MachineZero) {
+            if ((std::norm(func_ptr->data.c1[0]-inp.func_ptr->data.c1[0]) < MachineZero) and
+                c.imag() < MachineZero){
+                // real (possibly with complex c1)
+                for (int i = 0; i < inp.Ncomp(); i++) {
+                    CompD[i]->add_inplace(c.real(), *inp.CompD[i]);
+                }
+            }else if (func_ptr->data.c1[0].imag() > MachineZero or inp.func_ptr->data.c1[0].imag()  > MachineZero){
+                MSG_ABORT("Not implemented");
+            } else {
+                if(std::abs(1.0 - func_ptr->data.c1[0].real()) > MachineZero){
+                    rescale(func_ptr->data.c1[0].real());
+                    for (int i = 1; i < Ncomp(); i++) {
+                        if (std::norm(func_ptr->data.c1[0]-func_ptr->data.c1[i]) > MachineZero)
+                            MSG_ABORT("different scaling not implemented");
+                        func_ptr->data.c1[i] = {1.0, 0.0};
+                    }
+                    func_ptr->data.c1[0] = {1.0, 0.0};
+                }
+                for (int i = 0; i < inp.Ncomp(); i++) {
+                    CompD[i]->add_inplace(c.real()*inp.func_ptr->data.c1[i].real(), *inp.CompD[i]);
+                }
+            }
+        } else if (inp.isreal() and this->isreal() and c.imag() > MachineZero) {
+            MSG_ABORT("Not implemented");
+        } else if( this->iscomplex() and inp.iscomplex()) {
+            if (std::norm(1.0 - func_ptr->data.c1[0]) > MachineZero) {
+                rescale(func_ptr->data.c1[0]);
+                for (int i = 1; i < Ncomp(); i++) {
+                    if (std::norm(func_ptr->data.c1[0]-func_ptr->data.c1[i]) > MachineZero)
+                        MSG_ABORT("different scaling not implemented");
+                    func_ptr->data.c1[i] = {1.0, 0.0};
+                }
+                func_ptr->data.c1[0] = {1.0, 0.0};
+            }
+            for (int i = 0; i < inp.Ncomp(); i++) {
+                CompC[i]->add_inplace(c*inp.func_ptr->data.c1[i], *inp.CompC[i]);
+            }
+        } else if( this->isreal()) {
+            // we set as complex
+            for (int i = 0; i < Ncomp(); i++) {
                 CompD[i]->CopyTreeToComplex(CompC[i]);
                 delete CompD[i];
                 CompD[i] = nullptr;
-                func_ptr->iscomplex = true;
-                func_ptr->isreal = false;
             }
-            CompC[i]->add_inplace(c, *inp.CompC[i]);
+            func_ptr->iscomplex = 1;
+            func_ptr->isreal = 0;
+            if (std::norm(1.0-func_ptr->data.c1[0]) > MachineZero) {
+                rescale(func_ptr->data.c1[0]);
+                for (int i = 1; i < Ncomp(); i++) {
+                    if (std::norm(func_ptr->data.c1[0]-func_ptr->data.c1[i]) > MachineZero)
+                        MSG_ABORT("different scaling not implemented");
+                    func_ptr->data.c1[i] = {1.0, 0.0};
+                }
+                func_ptr->data.c1[0] = {1.0, 0.0};
+            }
+            for (int i = 0; i < inp.Ncomp(); i++) {
+                CompC[i]->add_inplace(c*inp.func_ptr->data.c1[i], *inp.CompC[i]);
+            }
+        } else {
+            MSG_ABORT("Not implemented");
         }
     }
 }
 
-template <int D> int CompFunction<D>::crop(double prec) {
+template <int D> int CompFunction<D>::crop(double prec, bool absPrec) {
     if (prec < 0.0) return 0;
     int nChunksremoved = 0;
     for (int i = 0; i < Ncomp(); i++) {
         if (isreal()) {
-            nChunksremoved += CompD[i]->crop(prec, 1.0, false);
+            nChunksremoved += CompD[i]->crop(prec, 1.0, absPrec);
         } else {
-            nChunksremoved += CompC[i]->crop(prec, 1.0, false);
+            nChunksremoved += CompC[i]->crop(prec, 1.0, absPrec);
         }
     }
     return nChunksremoved;
@@ -425,6 +481,26 @@ template class CompFunction<1>;
 template class CompFunction<2>;
 template class CompFunction<3>;
 
+/** @brief Deep copy that changes type from real to complex
+ *
+ * Deep copy: makes an exact copy with type complex from a real input
+ */
+template <int D> void CopyToComplex(CompFunction<D> &out, const CompFunction<D> &inp) {
+    out.func_ptr->data = inp.func_ptr->data;
+    out.defcomplex();
+    out.func_ptr->data.isreal = 0;
+    out.alloc(inp.Ncomp());
+    if (inp.getNNodes() == 0) return;
+    for (int i = 0; i < inp.Ncomp(); i++) {
+        if (inp.isreal()) {
+            inp.CompD[i]->CopyTreeToComplex(out.CompC[i]);
+        } else {
+            inp.CompC[i]->deep_copy(out.CompC[i]);
+        }
+    }
+}
+
+
 /** @brief Deep copy
  *
  * Deep copy: meta data is copied along with the content of each component.
@@ -545,6 +621,26 @@ template <int D> void linear_combination(CompFunction<D> &out, const std::vector
     }
 }
 
+/** @brief out = conj(inp) * inp
+ *
+ *  Note that output is always real
+ *
+ */
+template <int D> void make_density(CompFunction<D> &out, CompFunction<D> inp, double prec) {
+    multiply(prec, out, 1.0, inp, inp, -1, false, false, true);
+    if (out.iscomplex()) {
+        // copy onto real components
+        for (int i = 0; i < out.Ncomp(); i++) {
+            out.CompD[i] = out.CompC[i]->Real();
+            delete out.CompC[i];
+            out.CompC[i] = nullptr;
+        }
+        out.func_ptr->isreal = 1;
+        out.func_ptr->iscomplex = 0;
+    }
+}
+
+
 /** @brief out = inp_a * inp_b
  *
  */
@@ -607,6 +703,8 @@ template <int D> void multiply(double prec, CompFunction<D> &out, double coef, C
                     out.func_ptr->isreal = 0;
                     delete out.CompD[comp];
                     delete out.CompC[comp];
+                    out.CompD[comp] = nullptr;
+                    out.CompC[comp] = nullptr;
                     if (!out_allocated) out.alloc(out.Ncomp());
                     build_grid(*out.CompC[comp], *inp_a.CompC[comp]);
                     build_grid(*out.CompC[comp], *inp_b.CompC[comp]);
@@ -707,7 +805,7 @@ template <int D> void multiply(CompFunction<D> &out, FunctionTree<D, ComplexDoub
 /** @brief Compute <bra|ket> = int bra^\dag(r) * ket(r) dr.
  *
  *  Sum of component dots.
- *  Notice that the <bra| position is already complex conjugated.
+ *  Notice that the <bra| position is complex conjugated in the tree multiplication.
  *
  */
 template <int D> ComplexDouble dot(CompFunction<D> bra, CompFunction<D> ket) {
@@ -725,14 +823,10 @@ template <int D> ComplexDouble dot(CompFunction<D> bra, CompFunction<D> ket) {
         } else {
             dotprod += mrcpp::dot(*bra.CompC[comp], *ket.CompC[comp]);
         }
-        dotprod *= bra.func_ptr->data.c1[comp] * ket.func_ptr->data.c1[comp];
+        dotprod *= std::conj(bra.func_ptr->data.c1[comp]) * ket.func_ptr->data.c1[comp];
         dotprodtot += dotprod;
     }
-    if (bra.isreal() and ket.isreal()) {
-        return dotprodtot.real();
-    } else {
-        return dotprodtot;
-    }
+    return dotprodtot;
 }
 
 /** @brief Compute  <bra|ket> = int |bra^\dag(r)| * |ket(r)| dr.
@@ -767,7 +861,15 @@ void project(CompFunction<3> &out, std::function<double(const Coord<3> &r)> f, d
     mpi::share_function(out, 0, 123123, mpi::comm_share);
 }
 
-// template <int D, typename T>
+void project_real(CompFunction<3> &out, std::function<double(const Coord<3> &r)> f, double prec) {
+    bool need_to_project = not(out.isShared()) or mpi::share_master();
+    out.func_ptr->isreal = 1;
+    out.func_ptr->iscomplex = 0;
+    if (out.Ncomp() < 1) out.alloc(1);
+    if (need_to_project) mrcpp::project<3>(prec, *out.CompD[0], f);
+    mpi::share_function(out, 0, 123123, mpi::comm_share);
+}
+
 void project(CompFunction<3> &out, std::function<ComplexDouble(const Coord<3> &r)> f, double prec) {
     bool need_to_project = not(out.isShared()) or mpi::share_master();
     out.func_ptr->isreal = 0;
@@ -777,11 +879,21 @@ void project(CompFunction<3> &out, std::function<ComplexDouble(const Coord<3> &r
     mpi::share_function(out, 0, 123123, mpi::comm_share);
 }
 
+void project_cplx(CompFunction<3> &out, std::function<ComplexDouble(const Coord<3> &r)> f, double prec) {
+    bool need_to_project = not(out.isShared()) or mpi::share_master();
+    out.func_ptr->isreal = 0;
+    out.func_ptr->iscomplex = 1;
+    if (out.Ncomp() < 1) out.alloc(1);
+    if (need_to_project) mrcpp::project<3>(prec, *out.CompC[0], f);
+    mpi::share_function(out, 0, 123123, mpi::comm_share);
+}
+
 template <int D> void project(CompFunction<D> &out, RepresentableFunction<D, double> &f, double prec) {
     bool need_to_project = not(out.isShared()) or mpi::share_master();
     out.func_ptr->isreal = 1;
     out.func_ptr->iscomplex = 0;
     if (out.Ncomp() < 1) out.alloc(1);
+    if (need_to_project) build_grid(*out.CompD[0], f);
     if (need_to_project) mrcpp::project<D, double>(prec, *out.CompD[0], f);
     mpi::share_function(out, 0, 132231, mpi::comm_share);
 }
@@ -790,6 +902,7 @@ template <int D> void project(CompFunction<D> &out, RepresentableFunction<D, Com
     out.func_ptr->isreal = 0;
     out.func_ptr->iscomplex = 1;
     if (out.Ncomp() < 1) out.alloc(1);
+    if (need_to_project) build_grid(*out.CompC[0], f);
     if (need_to_project) mrcpp::project<D, ComplexDouble>(prec, *out.CompC[0], f);
     mpi::share_function(out, 0, 132231, mpi::comm_share);
 }
@@ -823,7 +936,10 @@ void rotate_cplx(CompFunctionVector &Phi, const ComplexMatrix &U, CompFunctionVe
     int N = Phi.size();
     int M = Psi.size();
     for (int i = 0; i < M; i++) {
-        for (int j; j < 4; j++) delete Psi[i].CompD[j];
+        for (int j = 0; j < 4; j++) {
+            delete Psi[i].CompD[j];
+            Psi[i].CompD[j] = nullptr;
+        }
         Psi[i].func_ptr->isreal = 0;
         Psi[i].func_ptr->iscomplex = 1;
     }
@@ -915,6 +1031,7 @@ void rotate_cplx(CompFunctionVector &Phi, const ComplexMatrix &U, CompFunctionVe
             int node_ix = indexVec_ref[n]; // SerialIx for this node in the reference tree
             // 4a) make a dense contiguous matrix with the coefficient from all the orbitals using node n
             std::vector<int> orbjVec; // to remember which orbital correspond to each orbVec.size();
+            if (node2orbVec.count(node_ix) == 0) continue;
             if (node2orbVec[node_ix].size() <= 0) continue;
             csize = sizecoeffW;
             if (parindexVec_ref[n] < 0) csize = sizecoeff; // for root nodes we include scaling coeff
@@ -922,6 +1039,7 @@ void rotate_cplx(CompFunctionVector &Phi, const ComplexMatrix &U, CompFunctionVe
             int shift = sizecoeff - sizecoeffW; // to copy only wavelet part
             if (parindexVec_ref[n] < 0) shift = 0;
             ComplexMatrix coeffBlock(csize, node2orbVec[node_ix].size());
+
             for (int j : node2orbVec[node_ix]) { // loop over indices of the orbitals using this node
                 int orb_node_ix = orb2node[j][node_ix];
                 for (int k = 0; k < csize; k++) coeffBlock(k, orbjVec.size()) = coeffVec[j][orb_node_ix][k + shift];
@@ -1203,6 +1321,7 @@ void rotate(CompFunctionVector &Phi, const ComplexMatrix &U, CompFunctionVector
             int node_ix = indexVec_ref[n]; // SerialIx for this node in the reference tree
             // 4a) make a dense contiguous matrix with the coefficient from all the orbitals using node n
             std::vector<int> orbjVec; // to remember which orbital correspond to each orbVec.size();
+            if (node2orbVec.count(node_ix) == 0) continue;
             if (node2orbVec[node_ix].size() <= 0) continue;
             csize = sizecoeffW;
             if (parindexVec_ref[n] < 0) csize = sizecoeff; // for root nodes we include scaling coeff
@@ -1398,7 +1517,6 @@ void rotate(CompFunctionVector &Phi, const ComplexMatrix &U, double prec) {
 void save_nodes(CompFunctionVector &Phi, FunctionTree<3> &refTree, BankAccount &account, int sizes) {
     int sizecoeff = (1 << refTree.getDim()) * refTree.getKp1_d();
     int sizecoeffW = ((1 << refTree.getDim()) - 1) * refTree.getKp1_d();
-    int max_nNodes = refTree.getNNodes();
     std::vector<double *> coeffVec;
     std::vector<ComplexDouble *> coeffVec_cplx;
     std::vector<double> scalefac;
@@ -1475,7 +1593,6 @@ CompFunctionVector multiply(CompFunctionVector &Phi, RepresentableFunction<3> &f
     // refine_grid(refTree, f); //to test
     mpi::allreduce_Tree_noCoeff(refTree, Phi, mpi::comm_wrk);
 
-    int kp1 = refTree.getKp1();
     int kp1_d = refTree.getKp1_d();
     int nCoefs = refTree.getTDim() * kp1_d;
 
@@ -1564,11 +1681,12 @@ CompFunctionVector multiply(CompFunctionVector &Phi, RepresentableFunction<3> &f
         for (int n = 0; n < max_n; n++) {
             MWNode<D> node(*(refNodes[n]), false);
             int node_ix = indexVec_ref[n]; // SerialIx for this node in the reference tree
+            if (node2orbVec.count(node_ix) == 0) continue;
 
             // 3a) make values for f at this node
             // 3a1) get coordinates of quadrature points for this node
             Eigen::MatrixXd pts; // Eigen::Zero(D, nCoefs);
-            double fval[nCoefs];
+            std::vector<double> fval(nCoefs);
             Coord<D> r;
             double *originalCoef = nullptr;
             MWNode<3> *Fnode = nullptr;
@@ -1589,7 +1707,7 @@ CompFunctionVector multiply(CompFunctionVector &Phi, RepresentableFunction<3> &f
                 } else {
                     originalCoef = Fnode->getCoefs();
                     for (int j = 0; j < nCoefs; j++) fval[j] = originalCoef[j];
-                    Fnode->attachCoefs(fval); // note that each thread has its own copy
+                    Fnode->attachCoefs(fval.data()); // note that each thread has its own copy
                     Fnode->mwTransform(Reconstruction);
                     Fnode->cvTransform(Forward);
                 }
@@ -1635,8 +1753,6 @@ CompFunctionVector multiply(CompFunctionVector &Phi, RepresentableFunction<3> &f
         }
     } else {
         // MPI
-        int count1 = 0;
-        int count2 = 0;
         TaskManager tasks(max_n);
         for (int nn = 0; nn < max_n; nn++) {
             int n = tasks.next_task();
@@ -1646,7 +1762,7 @@ CompFunctionVector multiply(CompFunctionVector &Phi, RepresentableFunction<3> &f
             // 3a1) get coordinates of quadrature points for this node
             Eigen::MatrixXd pts;           // Eigen::Zero(D, nCoefs);
             node.getExpandedChildPts(pts); // TODO: use getPrimitiveChildPts (less cache).
-            double fval[nCoefs];
+            std::vector<double> fval(nCoefs);
             Coord<D> r;
             MWNode<D> Fnode(*(refNodes[n]), false);
             if (Func == nullptr) {
@@ -1656,17 +1772,15 @@ CompFunctionVector multiply(CompFunctionVector &Phi, RepresentableFunction<3> &f
                 }
             } else {
                 int nIdx = Func->real().getIx(node.getNodeIndex());
-                count1++;
                 if (nIdx < 0) {
                     // use the function f instead of Func
-                    count2++;
                     for (int j = 0; j < nCoefs; j++) {
                         for (int d = 0; d < D; d++) r[d] = pts(d, j);
                         fval[j] = f.evalf(r);
                     }
                 } else {
-                    Func->real().getNodeCoeff(nIdx, fval); // fetch coef from Bank
-                    Fnode.attachCoefs(fval);
+                    Func->real().getNodeCoeff(nIdx, fval.data()); // fetch coef from Bank
+                    Fnode.attachCoefs(fval.data());
                     Fnode.mwTransform(Reconstruction);
                     Fnode.cvTransform(Forward);
                 }
@@ -1874,7 +1988,6 @@ ComplexMatrix calc_overlap_matrix_cplx(CompFunctionVector &BraKet) {
     }
 
     // 3) make dot product for all the nodes and accumulate into S
-    int ibank = 0;
 #pragma omp parallel if (serial)
     {
         ComplexMatrix S_omp = ComplexMatrix::Zero(N, N); // copy for each thread
@@ -1885,6 +1998,7 @@ ComplexMatrix calc_overlap_matrix_cplx(CompFunctionVector &BraKet) {
             int csize;
             int node_ix = indexVec_ref[n]; // SerialIx for this node in the reference tree
             std::vector<int> orbVec;       // identifies which orbitals use this node
+            if (serial and node2orbVec.count(node_ix) == 0) continue;
             if (serial and node2orbVec[node_ix].size() <= 0) continue;
             if (parindexVec_ref[n] < 0)
                 csize = sizecoeff;
@@ -1926,7 +2040,7 @@ ComplexMatrix calc_overlap_matrix_cplx(CompFunctionVector &BraKet) {
                             if (BraKet[orbVec[i]].func_ptr->data.n1[0] != BraKet[orbVec[j]].func_ptr->data.n1[0] and BraKet[orbVec[i]].func_ptr->data.n1[0] != 0 and
                                 BraKet[orbVec[j]].func_ptr->data.n1[0] != 0)
                                 continue;
-                            S_omp(orbVec[i], orbVec[j]) += S_temp(i, j);
+                            S(orbVec[i], orbVec[j]) += S_temp(i, j);
                         }
                     }
                 }
@@ -2022,7 +2136,6 @@ ComplexMatrix calc_overlap_matrix(CompFunctionVector &BraKet) {
     }
 
     // 3) make dot product for all the nodes and accumulate into S
-    int ibank = 0;
 #pragma omp parallel if (serial)
     {
         ComplexMatrix S_omp = ComplexMatrix::Zero(N, N); // copy for each thread
@@ -2033,6 +2146,7 @@ ComplexMatrix calc_overlap_matrix(CompFunctionVector &BraKet) {
             int csize;
             int node_ix = indexVec_ref[n]; // SerialIx for this node in the reference tree
             std::vector<int> orbVec;       // identifies which orbitals use this node
+            if (serial and node2orbVec.count(node_ix) == 0) continue;
             if (serial and node2orbVec[node_ix].size() <= 0) continue;
             if (parindexVec_ref[n] < 0)
                 csize = sizecoeff;
@@ -2227,7 +2341,6 @@ ComplexMatrix calc_overlap_matrix_cplx(CompFunctionVector &Bra, CompFunctionVect
     int totsiz = 0;
     int totget = 0;
     int mxtotsiz = 0;
-    int ibank = 0;
     // the omp crashes sometime for unknown reasons?
 #pragma omp parallel if (serial)
     {
@@ -2245,6 +2358,8 @@ ComplexMatrix calc_overlap_matrix_cplx(CompFunctionVector &Bra, CompFunctionVect
                 csize = sizecoeffW;
             if (serial) {
                 int node_ix = indexVec_ref[n];      // SerialIx for this node in the reference tree
+                if (node2orbVecBra.count(node_ix) == 0) continue;
+                if (node2orbVecKet.count(node_ix) == 0) continue;
                 int shift = sizecoeff - sizecoeffW; // to copy only wavelet part
                 ComplexMatrix coeffBlockBra(csize, node2orbVecBra[node_ix].size());
                 ComplexMatrix coeffBlockKet(csize, node2orbVecKet[node_ix].size());
@@ -2454,7 +2569,6 @@ ComplexMatrix calc_overlap_matrix(CompFunctionVector &Bra, CompFunctionVector &K
     int totsiz = 0;
     int totget = 0;
     int mxtotsiz = 0;
-    int ibank = 0;
 #pragma omp parallel if (serial)
     {
         DoubleMatrix S_omp = DoubleMatrix::Zero(N, M); // copy for each thread
@@ -2471,6 +2585,8 @@ ComplexMatrix calc_overlap_matrix(CompFunctionVector &Bra, CompFunctionVector &K
                 csize = sizecoeffW;
             if (serial) {
                 int node_ix = indexVec_ref[n];      // SerialIx for this node in the reference tree
+                if (node2orbVecBra.count(node_ix) == 0) continue;
+                if (node2orbVecKet.count(node_ix) == 0) continue;
                 int shift = sizecoeff - sizecoeffW; // to copy only wavelet part
                 DoubleMatrix coeffBlockBra(csize, node2orbVecBra[node_ix].size());
                 DoubleMatrix coeffBlockKet(csize, node2orbVecKet[node_ix].size());
@@ -2603,15 +2719,18 @@ template ComplexDouble dot(CompFunction<3> bra, CompFunction<3> ket);
 template void project(CompFunction<3> &out, RepresentableFunction<3, double> &f, double prec);
 template void project(CompFunction<3> &out, RepresentableFunction<3, ComplexDouble> &f, double prec);
 template void multiply(CompFunction<3> &out, CompFunction<3> inp_a, CompFunction<3> inp_b, double prec, bool absPrec, bool useMaxNorms, bool conjugate);
+template void multiply(double prec, CompFunction<3> &out, double coef, CompFunction<3> inp_a, CompFunction<3> inp_b, int maxIter = -1, bool absPrec = false, bool useMaxNorms = false, bool conjugate = false);
 template void multiply(CompFunction<3> &out, FunctionTree<3, double> &inp_a, RepresentableFunction<3, double> &f, double prec, int nrefine = 0, bool conjugate);
 template void multiply(CompFunction<3> &out, FunctionTree<3, ComplexDouble> &inp_a, RepresentableFunction<3, ComplexDouble> &f, double prec, int nrefine = 0, bool conjugate);
 template void multiply(CompFunction<3> &out, CompFunction<3> &inp_a, RepresentableFunction<3, double> &f, double prec, int nrefine = 0, bool conjugate);
 template void multiply(CompFunction<3> &out, CompFunction<3> &inp_a, RepresentableFunction<3, ComplexDouble> &f, double prec, int nrefine = 0, bool conjugate);
+template void CopyToComplex(CompFunction<3> &out, const CompFunction<3> &inp);
 template void deep_copy(CompFunction<3> *out, const CompFunction<3> &inp);
 template void deep_copy(CompFunction<3> &out, const CompFunction<3> &inp);
 template void add(CompFunction<3> &out, ComplexDouble a, CompFunction<3> inp_a, ComplexDouble b, CompFunction<3> inp_b, double prec, bool conjugate);
 template void linear_combination(CompFunction<3> &out, const std::vector<ComplexDouble> &c, std::vector<CompFunction<3>> &inp, double prec, bool conjugate);
 template double node_norm_dot(CompFunction<3> bra, CompFunction<3> ket);
 template void orthogonalize(double prec, CompFunction<3> &Bra, CompFunction<3> &Ket);
+template void make_density(CompFunction<3> &out, CompFunction<3> inp, double prec);
 
 } // namespace mrcpp
diff --git a/src/utils/CompFunction.h b/src/utils/CompFunction.h
index 82209819f..3b63a63bc 100644
--- a/src/utils/CompFunction.h
+++ b/src/utils/CompFunction.h
@@ -115,10 +115,14 @@ template <int D> class CompFunction {
     int conj() const { return func_ptr->data.conj; }           // soft conjugate
     int isreal() const { return func_ptr->data.isreal; }       // T=double
     int iscomplex() const { return func_ptr->data.iscomplex; } // T=DoubleComplex
-    void defreal() { func_ptr->data.isreal = 1; }              // define as real
-    void defcomplex() { func_ptr->data.iscomplex = 1; }        // define as complex
+    void defreal() { func_ptr->data.isreal = 1;                // define as real
+                     func_ptr->data.iscomplex = 0;}
+    void defcomplex() { func_ptr->data.isreal = 0;             // define as complex
+                        func_ptr->data.iscomplex = 1;}
     int share() const { return func_ptr->data.shared; }
     int *Nchunks() const { return func_ptr->data.Nchunks; } // number of chunks of each component tree
+    ComplexDouble getFac() const { return func_ptr->data.c1[0]; } // returns the overall multiplicative factor
+    void setFac(ComplexDouble fac) { func_ptr->data.c1[0] = fac; } // sets the overall multiplicative factor
 
     CompFunction paramCopy(bool alloc = false) const;
     ComplexDouble integrate() const;
@@ -132,7 +136,7 @@ template <int D> class CompFunction {
     const int getRank() const { return func_ptr->rank; };
     void add(ComplexDouble c, CompFunction<D> inp);
 
-    int crop(double prec);
+    int crop(double prec, bool absPrec = true);
     void rescale(ComplexDouble c);
     void free();
     int getSizeNodes() const;
@@ -158,6 +162,7 @@ template <int D> class CompFunction {
     std::shared_ptr<mrcpp::TreePtr<D>> func_ptr;
 };
 
+template <int D> void CopyToComplex(CompFunction<D> &out, const CompFunction<D> &inp);
 template <int D> void deep_copy(CompFunction<D> *out, const CompFunction<D> &inp);
 template <int D> void deep_copy(CompFunction<D> &out, const CompFunction<D> &inp);
 template <int D> void add(CompFunction<D> &out, ComplexDouble a, CompFunction<D> inp_a, ComplexDouble b, CompFunction<D> inp_b, double prec, bool conjugate = false);
@@ -165,6 +170,7 @@ template <int D> void linear_combination(CompFunction<D> &out, const std::vector
 template <int D> void multiply(CompFunction<D> &out, CompFunction<D> inp_a, CompFunction<D> inp_b, double prec, bool absPrec = false, bool useMaxNorms = false, bool conjugate = false);
 template <int D>
 void multiply(double prec, CompFunction<D> &out, double coef, CompFunction<D> inp_a, CompFunction<D> inp_b, int maxIter = -1, bool absPrec = false, bool useMaxNorms = false, bool conjugate = false);
+template <int D> void make_density(CompFunction<D> &out, CompFunction<D> inp, double prec);
 template <int D> void multiply(CompFunction<D> &out, CompFunction<D> inp_a, CompFunction<D> inp_b, bool absPrec = false, bool useMaxNorms = false, bool conjugate = false);
 template <int D> void multiply(CompFunction<D> &out, CompFunction<D> &inp_a, RepresentableFunction<D, double> &f, double prec, int nrefine = 0, bool conjugate = false);
 template <int D> void multiply(CompFunction<D> &out, CompFunction<D> &inp_a, RepresentableFunction<D, ComplexDouble> &f, double prec, int nrefine = 0, bool conjugate = false);
@@ -173,7 +179,9 @@ template <int D> void multiply(CompFunction<D> &out, FunctionTree<D, ComplexDoub
 template <int D> ComplexDouble dot(CompFunction<D> bra, CompFunction<D> ket);
 template <int D> double node_norm_dot(CompFunction<D> bra, CompFunction<D> ket);
 void project(CompFunction<3> &out, std::function<double(const Coord<3> &r)> f, double prec);
+void project_real(CompFunction<3> &out, std::function<double(const Coord<3> &r)> f, double prec); //overload of project is not always recognized by the compiler
 void project(CompFunction<3> &out, std::function<ComplexDouble(const Coord<3> &r)> f, double prec);
+void project_cplx(CompFunction<3> &out, std::function<ComplexDouble(const Coord<3> &r)> f, double prec); //overload of project is not always recognized by the compiler
 template <int D> void project(CompFunction<D> &out, RepresentableFunction<D, double> &f, double prec);
 template <int D> void project(CompFunction<D> &out, RepresentableFunction<D, ComplexDouble> &f, double prec);
 template <int D> void orthogonalize(double prec, CompFunction<D> &Bra, CompFunction<D> &Ket);
diff --git a/src/utils/tree_utils.cpp b/src/utils/tree_utils.cpp
index 333544f6e..b01c03347 100644
--- a/src/utils/tree_utils.cpp
+++ b/src/utils/tree_utils.cpp
@@ -118,8 +118,8 @@ template <int D, typename T> void tree_utils::mw_transform(const MWTree<D, T> &t
     int kp1_dm1 = math_utils::ipow(kp1, D - 1);
     const MWFilter &filter = tree.getMRA().getFilter();
     double overwrite = 0.0;
-    T tmpcoeff[kp1_d * tDim];
-    T tmpcoeff2[kp1_d * tDim];
+    std::vector<T> tmpcoeff(kp1_d * tDim);
+    std::vector<T> tmpcoeff2(kp1_d * tDim);
     int ftlim = tDim;
     int ftlim2 = tDim;
     int ftlim3 = tDim;
@@ -135,7 +135,7 @@ template <int D, typename T> void tree_utils::mw_transform(const MWTree<D, T> &t
     int i = 0;
     int mask = 1;
     for (int gt = 0; gt < tDim; gt++) {
-        T *out = tmpcoeff + gt * kp1_d;
+        T *out = tmpcoeff.data() + gt * kp1_d;
         for (int ft = 0; ft < ftlim; ft++) {
             // Operate in direction i only if the bits along other
             // directions are identical. The bit of the direction we
@@ -155,13 +155,13 @@ template <int D, typename T> void tree_utils::mw_transform(const MWTree<D, T> &t
         i++;
         mask = 2; // 1 << i;
         for (int gt = 0; gt < tDim; gt++) {
-            T *out = tmpcoeff2 + gt * kp1_d;
+            T *out = tmpcoeff2.data() + gt * kp1_d;
             for (int ft = 0; ft < ftlim2; ft++) {
                 // Operate in direction i only if the bits along other
                 // directions are identical. The bit of the direction we
                 // operate on determines the appropriate filter/operator
                 if ((gt | mask) == (ft | mask)) {
-                    T *in = tmpcoeff + ft * kp1_d;
+                    T *in = tmpcoeff.data() + ft * kp1_d;
                     int filter_index = 2 * ((gt >> i) & 1) + ((ft >> i) & 1);
                     const Eigen::MatrixXd &oper = filter.getSubFilter(filter_index, operation);
 
@@ -184,7 +184,7 @@ template <int D, typename T> void tree_utils::mw_transform(const MWTree<D, T> &t
                 // directions are identical. The bit of the direction we
                 // operate on determines the appropriate filter/operator
                 if ((gt | mask) == (ft | mask)) {
-                    T *in = tmpcoeff2 + ft * kp1_d;
+                    T *in = tmpcoeff2.data() + ft * kp1_d;
                     int filter_index = 2 * ((gt >> i) & 1) + ((ft >> i) & 1);
                     const Eigen::MatrixXd &oper = filter.getSubFilter(filter_index, operation);
 
@@ -201,8 +201,8 @@ template <int D, typename T> void tree_utils::mw_transform(const MWTree<D, T> &t
 
     if (D < 3) {
         T *out;
-        if (D == 1) out = tmpcoeff;
-        if (D == 2) out = tmpcoeff2;
+        if (D == 1) out = tmpcoeff.data();
+        if (D == 2) out = tmpcoeff2.data();
         if (b_overwrite) {
             for (int j = 0; j < tDim; j++) {
                 for (int i = 0; i < kp1_d; i++) { coeff_out[i + j * stride] = out[i + j * kp1_d]; }
@@ -234,7 +234,7 @@ template <typename T> void tree_utils::mw_transform_back(MWTree<3, T> &tree, T *
     int kp1_dm1 = math_utils::ipow(kp1, 2);
     const MWFilter &filter = tree.getMRA().getFilter();
     double overwrite = 0.0;
-    T tmpcoeff[kp1_d * tDim];
+    std::vector<T> tmpcoeff(kp1_d * tDim);
 
     int ftlim = tDim;
     int ftlim2 = tDim;
@@ -262,7 +262,7 @@ template <typename T> void tree_utils::mw_transform_back(MWTree<3, T> &tree, T *
     i++;
     mask = 2; // 1 << i;
     for (int gt = 0; gt < tDim; gt++) {
-        T *out = tmpcoeff + gt * kp1_d;
+        T *out = tmpcoeff.data() + gt * kp1_d;
         for (int ft = 0; ft < ftlim2; ft++) {
             // Operate in direction i only if the bits along other
             // directions are identical. The bit of the direction we
@@ -288,7 +288,7 @@ template <typename T> void tree_utils::mw_transform_back(MWTree<3, T> &tree, T *
             // directions are identical. The bit of the direction we
             // operate on determines the appropriate filter/operator
             if ((gt | mask) == (ft | mask)) {
-                T *in = tmpcoeff + ft * kp1_d;
+                T *in = tmpcoeff.data() + ft * kp1_d;
                 int filter_index = 2 * ((gt >> i) & 1) + ((ft >> i) & 1);
                 const Eigen::MatrixXd &oper = filter.getSubFilter(filter_index, operation);
 
diff --git a/tests/factory_functions.h b/tests/factory_functions.h
index d1c14a98f..571123263 100644
--- a/tests/factory_functions.h
+++ b/tests/factory_functions.h
@@ -129,7 +129,6 @@ void initialize(mrcpp::GaussFunc<D> **func) {
     double beta = 1.0e4;
     double alpha = std::pow(beta / mrcpp::pi, D / 2.0);
     double pos_data[3] = {-0.2, 0.5, 1.0};
-    const auto pow = std::array<int, D>{};
 
     auto pos = mrcpp::details::convert_to_std_array<double, D>(pos_data);
 
diff --git a/tests/functions/polynomial.cpp b/tests/functions/polynomial.cpp
index 48b09cc51..e643d1639 100644
--- a/tests/functions/polynomial.cpp
+++ b/tests/functions/polynomial.cpp
@@ -223,7 +223,6 @@ SCENARIO("Polynomials can be added and multiplied", "[poly_arithmetics], [polyno
             Polynomial R;
             R = P * Q;
             THEN("The coefficients of R are known") {
-                VectorXd &cr = R.getCoefs();
                 REQUIRE(R.getCoefs()[0] == Catch::Approx(0.0));
                 REQUIRE(R.getCoefs()[1] == Catch::Approx(0.0));
                 REQUIRE(R.getCoefs()[2] == Catch::Approx(1.0));
diff --git a/tests/operators/helmholtz_operator.cpp b/tests/operators/helmholtz_operator.cpp
index 8f570691d..ce669e27c 100644
--- a/tests/operators/helmholtz_operator.cpp
+++ b/tests/operators/helmholtz_operator.cpp
@@ -52,7 +52,6 @@ TEST_CASE("Helmholtz' kernel", "[init_helmholtz], [helmholtz_operator], [mw_oper
     const double exp_prec = 1.0e-4;
     const double proj_prec = 1.0e-3;
     const double ccc_prec = 1.0e-3;
-    const double band_prec = 1.0e-3;
 
     const int n = -3;
     const int k = 5;
diff --git a/tests/operators/poisson_operator.cpp b/tests/operators/poisson_operator.cpp
index 23bb22a06..e6ad04f69 100644
--- a/tests/operators/poisson_operator.cpp
+++ b/tests/operators/poisson_operator.cpp
@@ -50,7 +50,6 @@ TEST_CASE("Initialize Poisson operator", "[init_poisson], [poisson_operator], [m
     const double exp_prec = 1.0e-4;
     const double proj_prec = 1.0e-3;
     const double ccc_prec = 1.0e-3;
-    const double band_prec = 1.0e-3;
 
     const int n = -3;
     const int k = 5;
@@ -161,7 +160,6 @@ TEST_CASE("Apply Periodic Poisson' operator", "[apply_periodic_Poisson], [poisso
 
     // 2.0*pi periodic in all dirs // UPDATE ME
     auto scaling_factor = std::array<double, 3>{pi, pi, pi};
-    auto periodic = true;
 
     auto corner = std::array<int, 3>{-1, -1, -1};
     auto boxes = std::array<int, 3>{2, 2, 2};
diff --git a/tests/treebuilders/map.cpp b/tests/treebuilders/map.cpp
index 6be2153f0..64ced1e9c 100644
--- a/tests/treebuilders/map.cpp
+++ b/tests/treebuilders/map.cpp
@@ -80,7 +80,6 @@ template <int D> void testMapping() {
 
     const double ref_int = ref_tree.integrate();
     const double ref_norm = ref_tree.getSquareNorm();
-    const double inp_int = inp_tree.integrate();
     const double inp_norm = inp_tree.getSquareNorm();
 
     FMap<double, double> fmap = [](double val) { return val * val; };
diff --git a/tests/trees/function_tree.cpp b/tests/trees/function_tree.cpp
index d9492342e..99e5e709a 100644
--- a/tests/trees/function_tree.cpp
+++ b/tests/trees/function_tree.cpp
@@ -69,8 +69,7 @@ SCENARIO("Generating FunctionTree nodes", "[function_tree_generating], [function
 }
 
 template <int D> void testGeneratedNodes() {
-    const int depth = 3;
-
+    unsigned int depth = 3;
     Coord<D> r;
     if (r.size() >= 1) r[0] = -0.3;
     if (r.size() >= 2) r[1] = 0.6;
@@ -86,7 +85,6 @@ template <int D> void testGeneratedNodes() {
 
     WHEN("a non-existing node is fetched") {
         MWNode<D> &node = tree.getNode(r, depth);
-
         THEN("there will be allocated GenNodes") {
             REQUIRE(tree.getNGenNodes() > 0);
 
@@ -95,6 +93,7 @@ template <int D> void testGeneratedNodes() {
                 THEN("there will be no GenNodes") { REQUIRE(tree.getNGenNodes() == 0); }
             }
         }
+        (void)&node; // Clean up the fetched node
     }
     finalize(&mra);
 }
diff --git a/tests/trees/node_index.cpp b/tests/trees/node_index.cpp
index 791c94396..78d684a6c 100644
--- a/tests/trees/node_index.cpp
+++ b/tests/trees/node_index.cpp
@@ -72,7 +72,6 @@ template <int D> void testConstructors() {
     }
 
     SECTION("Parent constructor") {
-        int i = D;
         auto cIdx = nIdx->parent();
         REQUIRE(cIdx.getScale() == (nIdx->getScale() - 1));
     }