From df921d439f396a58ab750dfa05c105d80c738812 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?J=C3=A9r=C3=B4me=20H=C3=A9nin?= <heninj@gmail.com>
Date: Wed, 25 Dec 2024 11:24:17 +0100
Subject: [PATCH] Change weighted Poisson algorithm

---
 src/colvargrid.cpp | 216 ++++++++++++++++++++++-----------------------
 1 file changed, 107 insertions(+), 109 deletions(-)

diff --git a/src/colvargrid.cpp b/src/colvargrid.cpp
index e2c9fc3d4..d89e94248 100644
--- a/src/colvargrid.cpp
+++ b/src/colvargrid.cpp
@@ -630,17 +630,34 @@ int integrate_potential::integrate(const int itmax, const cvm::real &tol, cvm::r
       cvm::log("Integrated in " + cvm::to_str(iter) + " steps, error: " + cvm::to_str(err));
 
     // DEBUG ###########################
-    auto backup = data;
-    data = divergence;
-    std::ofstream os("div.dat");
-    write_multicol(os);
-    os.close();
-    data = weights;
-    os.open("weights.dat");
-    write_multicol(os);
-    os.close();
-    data = backup;
-    // DEBUG ###########################
+    // auto backup = data;
+    // data = divergence;
+    // std::ofstream os("div.dat");
+    // write_multicol(os);
+    // os.close();
+    // data = weights;
+    // os.open("weights.dat");
+    // write_multicol(os);
+    // os.close();
+    // data = backup;
+    // DEBUG 2 ###########################
+    // Compute terms of the Laplacian matrix
+    std::vector<cvm::real> lap_mat(nt, 0.);
+
+    std::vector<size_t> cols = { 0, 1, 2, 3, 4, 5, nt-6, nt-5, nt-4, nt-3, nt-2, nt-1 };
+
+    for (size_t i = 0; i < cols.size(); i++) {
+      this->reset();
+      data[cols[i]] = 1.;
+      laplacian_weighted(data, lap_mat);
+      printf("Col  %3li  | ", cols[i]);
+      for (size_t j = 0; j < cols.size(); j++) {
+        printf(" %6.1f", lap_mat[cols[j]]);
+      }
+      printf("\n");
+    }
+    // DEBUG 2 ###########################
+
 
   } else {
     cvm::error("Cannot integrate PMF in dimension > 3\n");
@@ -742,8 +759,10 @@ void integrate_potential::update_div_local(const std::vector<int> &ix0)
     count += get_grad(g10, ix);
 
     // cvm::real weight = count / 1000.L + min_weight;
-    // cvm::real const weight = 10;
-    cvm::real weight = cvm::logn(count + 1) + 1;
+    cvm::real const weight = (linear_index > nt/2 ? 1.1 : 1);
+    // cvm::real const weight = 1.;
+
+    // cvm::real weight = cvm::logn(count + 1) * 1e-8 + 1;
 
     divergence[linear_index] = ((g10[0]-g00[0] + g11[0]-g01[0]) / widths[0]
                               + (g01[1]-g00[1] + g11[1]-g10[1]) / widths[1]) * 0.5 * weight;
@@ -1241,16 +1260,13 @@ void integrate_potential::laplacian_weighted(const std::vector<cvm::real> &A, st
 
     for (i=1; i<w-1; i++) {
       // Full range of j, but factor may change on y edges (j == 0 and j == h-1)
-      LA[li] = fact * ffx * ((weights[li] + weights[li + xm]) * (A[li + xm] - A[li]) +
-                             (weights[li] + weights[li + xp]) * (A[li + xp] - A[li]));
+      LA[li] = fact * ffx * (weights[li + xm] * A[li + xm] + weights[li + xp] * A[li + xp] - 2.0 * weights[li] * A[li]);
       li++;
       for (j=1; j<h-1; j++) {
-        LA[li] = ffx * ((weights[li] + weights[li + xm]) * (A[li + xm] - A[li]) +
-                        (weights[li] + weights[li + xp]) * (A[li + xp] - A[li]));
+        LA[li] = ffx * (weights[li + xm] * A[li + xm] + weights[li + xp] * A[li + xp] - 2.0 * weights[li] * A[li]);
         li++;
       }
-      LA[li] = fact * ffx * ((weights[li] + weights[li + xm]) * (A[li + xm] - A[li]) +
-                             (weights[li] + weights[li + xp]) * (A[li + xp] - A[li]));
+      LA[li] = fact * ffx * (weights[li + xm] * A[li + xm] + weights[li + xp] * A[li + xp] - 2.0 * weights[li] * A[li]);
       li++;
     }
     // Edges along x (x components only)
@@ -1260,43 +1276,37 @@ void integrate_potential::laplacian_weighted(const std::vector<cvm::real> &A, st
       xm =  h * (w - 1);
       xp =  h;
       fact = periodic[1] ? 1.0 : 0.5;
-      LA[li] = fact * ffx * ((weights[li] + weights[li + xm]) * (A[li + xm] - A[li]) +
-                             (weights[li] + weights[li + xp]) * (A[li + xp] - A[li]));
-      LA[li2] = fact * ffx * ((weights[li2] + weights[li2 - xm]) * (A[li2 - xm] - A[li2]) +
-                              (weights[li2] + weights[li2 - xp]) * (A[li2 - xp] - A[li2]));
+      LA[li]  = fact * ffx * (weights[li + xm] * A[li + xm] + weights[li + xp] * A[li + xp] - 2.0 * weights[li] * A[li]);
+      LA[li2] = fact * ffx * (weights[li2 - xp] * A[li2 - xp] + weights[li2 - xm] * A[li2 - xm] - 2.0 * weights[li2] * A[li2]);
       li++;
       li2++;
       for (j=1; j<h-1; j++) {
-        LA[li] = ffx * ((weights[li] + weights[li + xm]) * (A[li + xm] - A[li]) +
-                        (weights[li] + weights[li + xp]) * (A[li + xp] - A[li]));
-        LA[li2] = ffx * ((weights[li2] + weights[li2 - xm]) * (A[li2 - xm] - A[li2]) +
-                         (weights[li2] + weights[li2 - xp]) * (A[li2 - xp] - A[li2]));
+        LA[li]  = ffx * (weights[li + xm] * A[li + xm] + weights[li + xp] * A[li + xp] - 2.0 * weights[li] * A[li]);
+        LA[li2] = ffx * (weights[li2 - xp] * A[li2 - xp] + weights[li2 - xm] * A[li2 - xm] - 2.0 * weights[li2] * A[li2]);
         li++;
         li2++;
       }
-      LA[li] = fact * ffx * ((weights[li] + weights[li + xm]) * (A[li + xm] - A[li]) +
-                             (weights[li] + weights[li + xp]) * (A[li + xp] - A[li]));
-      LA[li2] = fact * ffx * ((weights[li2] + weights[li2 - xm]) * (A[li2 - xm] - A[li2]) +
-                              (weights[li2] + weights[li2 - xp]) * (A[li2 - xp] - A[li2]));
-    } else { // Neumann BC along x
+      LA[li]  = fact * ffx * (weights[li + xm] * A[li + xm] + weights[li + xp] * A[li + xp] - 2.0 * weights[li] * A[li]);
+      LA[li2] = fact * ffx * (weights[li2 - xp] * A[li2 - xp] + weights[li2 - xm] * A[li2 - xm] - 2.0 * weights[li2] * A[li2]);
+    } else {
       xm = -h;
       xp =  h;
       fact = periodic[1] ? 1.0 : 0.5; // Halve in corners in full PBC only
       // lower corner, "j == 0"
-      LA[li]  = fact * ffx * (weights[li] + weights[li + xp]) * (A[li + xp] - A[li]);
-      LA[li2] = fact * ffx * (weights[li2] + weights[li2 + xm]) * (A[li2 + xm] - A[li2]);
+      LA[li]  = fact * ffx * (weights[li + xp] * A[li + xp] - weights[li] * A[li]);
+      LA[li2] = fact * ffx * (weights[li2 + xm] * A[li2 + xm] - weights[li2] * A[li2]);
       li++;
       li2++;
       for (j=1; j<h-1; j++) {
         // x gradient (+ y term of laplacian, calculated below)
-        LA[li]  = ffx * (weights[li] + weights[li + xp]) * (A[li + xp] - A[li]);
-        LA[li2] = ffx * (weights[li2] + weights[li2 + xm]) * (A[li2 + xm] - A[li2]);
+        LA[li]  = ffx * (weights[li + xp] * A[li + xp] - weights[li] * A[li]);
+        LA[li2] = ffx * (weights[li2 + xm] * A[li2 + xm] - weights[li2] * A[li2]);
         li++;
         li2++;
       }
       // upper corner, j == h-1
-      LA[li]  = fact * ffx * (weights[li] + weights[li + xp]) * (A[li + xp] - A[li]);
-      LA[li2] = fact * ffx * (weights[li2] + weights[li2 + xm]) * (A[li2 + xm] - A[li2]);
+      LA[li]  = fact * ffx * (weights[li + xp] * A[li + xp] - weights[li] * A[li]);
+      LA[li2] = fact * ffx * (weights[li2 + xm] * A[li2 + xm] - weights[li2] * A[li2]);
     }
 
     // Now adding all y components
@@ -1309,8 +1319,7 @@ void integrate_potential::laplacian_weighted(const std::vector<cvm::real> &A, st
       if (i == 1) fact = 1.0;
       if (i == w - 1) fact = periodic[0] ? 1.0 : 0.5;
       for (j=1; j<h-1; j++) {
-        LA[li] += fact * ffy * ((weights[li] + weights[li + ym]) * (A[li + ym] - A[li]) +
-                                (weights[li] + weights[li + yp]) * (A[li + yp] - A[li]));
+        LA[li] += fact * ffy * (weights[li + ym] * A[li + ym] + weights[li + yp] * A[li + yp] - 2.0 * weights[li] * A[li]);
         li++;
       }
       li += 2; // skip the edges and move to next column
@@ -1322,43 +1331,37 @@ void integrate_potential::laplacian_weighted(const std::vector<cvm::real> &A, st
       fact = periodic[0] ? 1.0 : 0.5;
       ym = h - 1;
       yp = 1;
-      LA[li] += fact * ffy * ((weights[li] + weights[li + ym]) * (A[li + ym] - A[li]) +
-                              (weights[li] + weights[li + yp]) * (A[li + yp] - A[li]));
-      LA[li2] += fact * ffy * ((weights[li2] + weights[li2 - ym]) * (A[li2 - ym] - A[li2]) +
-                               (weights[li2] + weights[li2 - yp]) * (A[li2 - yp] - A[li2]));
+      LA[li]  += fact * ffy * (weights[li + ym] * A[li + ym] + weights[li + yp] * A[li + yp] - 2.0 * weights[li] * A[li]);
+      LA[li2] += fact * ffy * (weights[li2 - yp] * A[li2 - yp] + weights[li2 - ym] * A[li2 - ym] - 2.0 * weights[li2] * A[li2]);
       li  += h;
       li2 += h;
       for (i=1; i<w-1; i++) {
-        LA[li] += ffy * ((weights[li] + weights[li + ym]) * (A[li + ym] - A[li]) +
-                         (weights[li] + weights[li + yp]) * (A[li + yp] - A[li]));
-        LA[li2] += ffy * ((weights[li2] + weights[li2 - ym]) * (A[li2 - ym] - A[li2]) +
-                          (weights[li2] + weights[li2 - yp]) * (A[li2 - yp] - A[li2]));
+        LA[li]  += ffy * (weights[li + ym] * A[li + ym] + weights[li + yp] * A[li + yp] - 2.0 * weights[li] * A[li]);
+        LA[li2] += ffy * (weights[li2 - yp] * A[li2 - yp] + weights[li2 - ym] * A[li2 - ym] - 2.0 * weights[li2] * A[li2]);
         li  += h;
         li2 += h;
       }
-      LA[li] += fact * ffy * ((weights[li] + weights[li + ym]) * (A[li + ym] - A[li]) +
-                              (weights[li] + weights[li + yp]) * (A[li + yp] - A[li]));
-      LA[li2] += fact * ffy * ((weights[li2] + weights[li2 - ym]) * (A[li2 - ym] - A[li2]) +
-                               (weights[li2] + weights[li2 - yp]) * (A[li2 - yp] - A[li2]));
-    } else { // Neumann BC along y
+      LA[li]  += fact * ffy * (weights[li + ym] * A[li + ym] + weights[li + yp] * A[li + yp] - 2.0 * weights[li] * A[li]);
+      LA[li2] += fact * ffy * (weights[li2 - yp] * A[li2 - yp] + weights[li2 - ym] * A[li2 - ym] - 2.0 * weights[li2] * A[li2]);
+    } else {
       ym = -1;
       yp = 1;
       fact = periodic[0] ? 1.0 : 0.5; // Halve in corners in full PBC only
       // Left corner
-      LA[li]  += fact * ffy * (weights[li] + weights[li + yp]) * (A[li + yp] - A[li]);
-      LA[li2] += fact * ffy * (weights[li2] + weights[li2 + ym]) * (A[li2 + ym] - A[li2]);
+      LA[li]  += fact * ffy * (weights[li + yp] * A[li + yp] - weights[li] * A[li]);
+      LA[li2] += fact * ffy * (weights[li2 + ym] * A[li2 + ym] - weights[li2] * A[li2]);
       li  += h;
       li2 += h;
       for (i=1; i<w-1; i++) {
         // y gradient (+ x term of laplacian, calculated above)
-        LA[li]  += ffy * (weights[li] + weights[li + yp]) * (A[li + yp] - A[li]);
-        LA[li2] += ffy * (weights[li2] + weights[li2 + ym]) * (A[li2 + ym] - A[li2]);
+        LA[li]  += ffy * (weights[li + yp] * A[li + yp] - weights[li] * A[li]);
+        LA[li2] += ffy * (weights[li2 + ym] * A[li2 + ym] - weights[li2] * A[li2]);
         li  += h;
         li2 += h;
       }
       // Right corner
-      LA[li]  += fact * ffy * (weights[li] + weights[li + yp]) * (A[li + yp] - A[li]);
-      LA[li2] += fact * ffy * (weights[li2] + weights[li2 + ym]) * (A[li2 + ym] - A[li2]);
+      LA[li]  += fact * ffy * (weights[li + yp] * A[li + yp] - weights[li] * A[li]);
+      LA[li2] += fact * ffy * (weights[li2 + ym] * A[li2 + ym] - weights[li2] * A[li2]);
     }
 
   } else if (nd == 3) {
@@ -1395,15 +1398,15 @@ void integrate_potential::laplacian_weighted(const std::vector<cvm::real> &A, st
       for (j=0; j<d; j++) { // full range of y
         if (j == 1) fact *= ifacty;
         if (j == d-1) fact *= facty;
-        LA[li] = fact * ffx * (A[li + xm] + A[li + xp] - 2.0 * A[li]);
+        LA[li] = fact * ffx * (weights[li + xm] * A[li + xm] + weights[li + xp] * A[li + xp] - 2.0 * weights[li] * A[li]);
         li++;
         fact *= ifactz;
         for (k=1; k<h-1; k++) { // full range of z
-          LA[li] = fact * ffx * (A[li + xm] + A[li + xp] - 2.0 * A[li]);
+          LA[li] = fact * ffx * (weights[li + xm] * A[li + xm] + weights[li + xp] * A[li + xp] - 2.0 * weights[li] * A[li]);
           li++;
         }
         fact *= factz;
-        LA[li] = fact * ffx * (A[li + xm] + A[li + xp] - 2.0 * A[li]);
+        LA[li] = fact * ffx * (weights[li + xm] * A[li + xm] + weights[li + xp] * A[li + xp] - 2.0 * weights[li] * A[li]);
         li++;
       }
     }
@@ -1417,20 +1420,20 @@ void integrate_potential::laplacian_weighted(const std::vector<cvm::real> &A, st
       for (j=0; j<d; j++) {
         if (j == 1) fact *= ifacty;
         if (j == d-1) fact *= facty;
-        LA[li]  = fact * ffx * (A[li + xm] + A[li + xp] - 2.0 * A[li]);
-        LA[li2] = fact * ffx * (A[li2 - xp] + A[li2 - xm] - 2.0 * A[li2]);
+        LA[li]  = fact * ffx * (weights[li + xm] * A[li + xm] + weights[li + xp] * A[li + xp] - 2.0 * weights[li] * A[li]);
+        LA[li2] = fact * ffx * (weights[li2 - xp] * A[li2 - xp] + weights[li2 - xm] * A[li2 - xm] - 2.0 * weights[li2] * A[li2]);
         li++;
         li2++;
         fact *= ifactz;
         for (k=1; k<h-1; k++) {
-          LA[li]  = fact * ffx * (A[li + xm] + A[li + xp] - 2.0 * A[li]);
-          LA[li2] = fact * ffx * (A[li2 - xp] + A[li2 - xm] - 2.0 * A[li2]);
+          LA[li]  = fact * ffx * (weights[li + xm] * A[li + xm] + weights[li + xp] * A[li + xp] - 2.0 * weights[li] * A[li]);
+          LA[li2] = fact * ffx * (weights[li2 - xp] * A[li2 - xp] + weights[li2 - xm] * A[li2 - xm] - 2.0 * weights[li2] * A[li2]);
           li++;
           li2++;
         }
         fact *= factz;
-        LA[li]  = fact * ffx * (A[li + xm] + A[li + xp] - 2.0 * A[li]);
-        LA[li2] = fact * ffx * (A[li2 - xp] + A[li2 - xm] - 2.0 * A[li2]);
+        LA[li]  = fact * ffx * (weights[li + xm] * A[li + xm] + weights[li + xp] * A[li + xp] - 2.0 * weights[li] * A[li]);
+        LA[li2] = fact * ffx * (weights[li2 - xp] * A[li2 - xp] + weights[li2 - xm] * A[li2 - xm] - 2.0 * weights[li2] * A[li2]);
         li++;
         li2++;
       }
@@ -1441,21 +1444,21 @@ void integrate_potential::laplacian_weighted(const std::vector<cvm::real> &A, st
       for (j=0; j<d; j++) {
         if (j == 1) fact *= ifacty;
         if (j == d-1) fact *= facty;
-        LA[li]  = fact * ffx * (A[li + xp] - A[li]);
-        LA[li2] = fact * ffx * (A[li2 + xm] - A[li2]);
+        LA[li]  = fact * ffx * (weights[li + xp] * A[li + xp] - weights[li] * A[li]);
+        LA[li2] = fact * ffx * (weights[li2 + xm] * A[li2 + xm] - weights[li2] * A[li2]);
         li++;
         li2++;
         fact *= ifactz;
         for (k=1; k<h-1; k++) {
           // x gradient (+ y, z terms of laplacian, calculated below)
-          LA[li]  = fact * ffx * (A[li + xp] - A[li]);
-          LA[li2] = fact * ffx * (A[li2 + xm] - A[li2]);
+          LA[li]  = fact * ffx * (weights[li + xp] * A[li + xp] - weights[li] * A[li]);
+          LA[li2] = fact * ffx * (weights[li2 + xm] * A[li2 + xm] - weights[li2] * A[li2]);
           li++;
           li2++;
         }
         fact *= factz;
-        LA[li]  = fact * ffx * (A[li + xp] - A[li]);
-        LA[li2] = fact * ffx * (A[li2 + xm] - A[li2]);
+        LA[li]  = fact * ffx * (weights[li + xp] * A[li + xp] - weights[li] * A[li]);
+        LA[li2] = fact * ffx * (weights[li2 + xm] * A[li2 + xm] - weights[li2] * A[li2]);
         li++;
         li2++;
       }
@@ -1469,15 +1472,15 @@ void integrate_potential::laplacian_weighted(const std::vector<cvm::real> &A, st
       if (i == 1) fact *= ifactx;
       if (i == w-1) fact *= factx;
       for (j=1; j<d-1; j++) {
-        LA[li] += fact * ffy * (A[li + ym] + A[li + yp] - 2.0 * A[li]);
+        LA[li] += fact * ffy * (weights[li + ym] * A[li + ym] + weights[li + yp] * A[li + yp] - 2.0 * weights[li] * A[li]);
         li++;
         fact *= ifactz;
         for (k=1; k<h-1; k++) {
-          LA[li] += fact * ffy * (A[li + ym] + A[li + yp] - 2.0 * A[li]);
+          LA[li] += fact * ffy * (weights[li + ym] * A[li + ym] + weights[li + yp] * A[li + yp] - 2.0 * weights[li] * A[li]);
           li++;
         }
         fact *= factz;
-        LA[li] += fact * ffy * (A[li + ym] + A[li + yp] - 2.0 * A[li]);
+        LA[li] += fact * ffy * (weights[li + ym] * A[li + ym] + weights[li + yp] * A[li + yp] - 2.0 * weights[li] * A[li]);
         li++;
       }
       li += 2 * h; // skip columns in front and back slabs
@@ -1492,20 +1495,20 @@ void integrate_potential::laplacian_weighted(const std::vector<cvm::real> &A, st
       for (i=0; i<w; i++) {
         if (i == 1) fact *= ifactx;
         if (i == w-1) fact *= factx;
-        LA[li]  += fact * ffy * (A[li + ym] + A[li + yp] - 2.0 * A[li]);
-        LA[li2] += fact * ffy * (A[li2 - yp] + A[li2 - ym] - 2.0 * A[li2]);
+        LA[li]  += fact * ffy * (weights[li + ym] * A[li + ym] + weights[li + yp] * A[li + yp] - 2.0 * weights[li] * A[li]);
+        LA[li2] += fact * ffy * (weights[li2 - yp] * A[li2 - yp] + weights[li2 - ym] * A[li2 - ym] - 2.0 * weights[li2] * A[li2]);
         li++;
         li2++;
         fact *= ifactz;
         for (k=1; k<h-1; k++) {
-          LA[li]  += fact * ffy * (A[li + ym] + A[li + yp] - 2.0 * A[li]);
-          LA[li2] += fact * ffy * (A[li2 - yp] + A[li2 - ym] - 2.0 * A[li2]);
+          LA[li]  += fact * ffy * (weights[li + ym] * A[li + ym] + weights[li + yp] * A[li + yp] - 2.0 * weights[li] * A[li]);
+          LA[li2] += fact * ffy * (weights[li2 - yp] * A[li2 - yp] + weights[li2 - ym] * A[li2 - ym] - 2.0 * weights[li2] * A[li2]);
           li++;
           li2++;
         }
         fact *= factz;
-        LA[li]  += fact * ffy * (A[li + ym] + A[li + yp] - 2.0 * A[li]);
-        LA[li2] += fact * ffy * (A[li2 - yp] + A[li2 - ym] - 2.0 * A[li2]);
+        LA[li]  += fact * ffy * (weights[li + ym] * A[li + ym] + weights[li + yp] * A[li + yp] - 2.0 * weights[li] * A[li]);
+        LA[li2] += fact * ffy * (weights[li2 - yp] * A[li2 - yp] + weights[li2 - ym] * A[li2 - ym] - 2.0 * weights[li2] * A[li2]);
         li++;
         li2++;
         li  += h * static_cast<size_t>(d - 1);
@@ -1518,21 +1521,21 @@ void integrate_potential::laplacian_weighted(const std::vector<cvm::real> &A, st
       for (i=0; i<w; i++) {
         if (i == 1) fact *= ifactx;
         if (i == w-1) fact *= factx;
-        LA[li]  += fact * ffy * (A[li + yp] - A[li]);
-        LA[li2] += fact * ffy * (A[li2 + ym] - A[li2]);
+        LA[li]  += fact * ffy * (weights[li + yp] * A[li + yp] - weights[li] * A[li]);
+        LA[li2] += fact * ffy * (weights[li2 + ym] * A[li2 + ym] - weights[li2] * A[li2]);
         li++;
         li2++;
         fact *= ifactz;
         for (k=1; k<h-1; k++) {
           // y gradient (+ x, z terms of laplacian, calculated above and below)
-          LA[li]  += fact * ffy * (A[li + yp] - A[li]);
-          LA[li2] += fact * ffy * (A[li2 + ym] - A[li2]);
+          LA[li]  += fact * ffy * (weights[li + yp] * A[li + yp] - weights[li] * A[li]);
+          LA[li2] += fact * ffy * (weights[li2 + ym] * A[li2 + ym] - weights[li2] * A[li2]);
           li++;
           li2++;
         }
         fact *= factz;
-        LA[li]  += fact * ffy * (A[li + yp] - A[li]);
-        LA[li2] += fact * ffy * (A[li2 + ym] - A[li2]);
+        LA[li]  += fact * ffy * (weights[li + yp] * A[li + yp] - weights[li] * A[li]);
+        LA[li2] += fact * ffy * (weights[li2 + ym] * A[li2 + ym] - weights[li2] * A[li2]);
         li++;
         li2++;
         li  += h * static_cast<size_t>(d - 1);
@@ -1548,21 +1551,21 @@ void integrate_potential::laplacian_weighted(const std::vector<cvm::real> &A, st
       if (i == 1) fact *= ifactx;
       if (i == w-1) fact *= factx;
       for (k=1; k<h-1; k++) {
-        LA[li] += fact * ffz * (A[li + zm] + A[li + zp] - 2.0 * A[li]);
+        LA[li] += fact * ffz * (weights[li + zm] * A[li + zm] + weights[li + zp] * A[li + zp] - 2.0 * weights[li] * A[li]);
         li++;
       }
       fact *= ifacty;
       li += 2; // skip edge slabs
       for (j=1; j<d-1; j++) { // full range of y
         for (k=1; k<h-1; k++) {
-          LA[li] += fact * ffz * (A[li + zm] + A[li + zp] - 2.0 * A[li]);
+          LA[li] += fact * ffz * (weights[li + zm] * A[li + zm] + weights[li + zp] * A[li + zp] - 2.0 * weights[li] * A[li]);
           li++;
         }
         li += 2; // skip edge slabs
       }
       fact *= facty;
       for (k=1; k<h-1; k++) {
-        LA[li] += fact * ffz * (A[li + zm] + A[li + zp] - 2.0 * A[li]);
+        LA[li] += fact * ffz * (weights[li + zm] * A[li + zm] + weights[li + zp] * A[li + zp] - 2.0 * weights[li] * A[li]);
         li++;
       }
       li += 2; // skip edge slabs
@@ -1577,20 +1580,20 @@ void integrate_potential::laplacian_weighted(const std::vector<cvm::real> &A, st
       for (i=0; i<w; i++) {
         if (i == 1) fact *= ifactx;
         if (i == w-1) fact *= factx;
-        LA[li]  += fact * ffz * (A[li + zm] + A[li + zp] - 2.0 * A[li]);
-        LA[li2] += fact * ffz * (A[li2 - zp] + A[li2 - zm] - 2.0 * A[li2]);
+        LA[li]  += fact * ffz * (weights[li + zm] * A[li + zm] + weights[li + zp] * A[li + zp] - 2.0 * weights[li] * A[li]);
+        LA[li2] += fact * ffz * (weights[li2 - zp] * A[li2 - zp] + weights[li2 - zm] * A[li2 - zm] - 2.0 * weights[li2] * A[li2]);
         li  += h;
         li2 += h;
         fact *= ifacty;
         for (j=1; j<d-1; j++) {
-          LA[li]  += fact * ffz * (A[li + zm] + A[li + zp] - 2.0 * A[li]);
-          LA[li2] += fact * ffz * (A[li2 - zp] + A[li2 - zm] - 2.0 * A[li2]);
+          LA[li]  += fact * ffz * (weights[li + zm] * A[li + zm] + weights[li + zp] * A[li + zp] - 2.0 * weights[li] * A[li]);
+          LA[li2] += fact * ffz * (weights[li2 - zp] * A[li2 - zp] + weights[li2 - zm] * A[li2 - zm] - 2.0 * weights[li2] * A[li2]);
           li  += h;
           li2 += h;
         }
         fact *= facty;
-        LA[li]  += fact * ffz * (A[li + zm] + A[li + zp] - 2.0 * A[li]);
-        LA[li2] += fact * ffz * (A[li2 - zp] + A[li2 - zm] - 2.0 * A[li2]);
+        LA[li]  += fact * ffz * (weights[li + zm] * A[li + zm] + weights[li + zp] * A[li + zp] - 2.0 * weights[li] * A[li]);
+        LA[li2] += fact * ffz * (weights[li2 - zp] * A[li2 - zp] + weights[li2 - zm] * A[li2 - zm] - 2.0 * weights[li2] * A[li2]);
         li  += h;
         li2 += h;
       }
@@ -1601,31 +1604,26 @@ void integrate_potential::laplacian_weighted(const std::vector<cvm::real> &A, st
       for (i=0; i<w; i++) {
         if (i == 1) fact *= ifactx;
         if (i == w-1) fact *= factx;
-        LA[li]  += fact * ffz * (A[li + zp] - A[li]);
-        LA[li2] += fact * ffz * (A[li2 + zm] - A[li2]);
+        LA[li]  += fact * ffz * (weights[li + zp] * A[li + zp] - weights[li] * A[li]);
+        LA[li2] += fact * ffz * (weights[li2 + zm] * A[li2 + zm] - weights[li2] * A[li2]);
         li  += h;
         li2 += h;
         fact *= ifacty;
         for (j=1; j<d-1; j++) {
           // z gradient (+ x, y terms of laplacian, calculated above)
-          LA[li]  += fact * ffz * (A[li + zp] - A[li]);
-          LA[li2] += fact * ffz * (A[li2 + zm] - A[li2]);
+          LA[li]  += fact * ffz * (weights[li + zp] * A[li + zp] - weights[li] * A[li]);
+          LA[li2] += fact * ffz * (weights[li2 + zm] * A[li2 + zm] - weights[li2] * A[li2]);
           li  += h;
           li2 += h;
         }
         fact *= facty;
-        LA[li]  += fact * ffz * (A[li + zp] - A[li]);
-        LA[li2] += fact * ffz * (A[li2 + zm] - A[li2]);
+        LA[li]  += fact * ffz * (weights[li + zp] * A[li + zp] - weights[li] * A[li]);
+        LA[li2] += fact * ffz * (weights[li2 + zm] * A[li2 + zm] - weights[li2] * A[li2]);
         li  += h;
         li2 += h;
       }
     }
   }
-
-  // Multiply all terms by weight - simple Laplacian scaling
-  // for (size_t li = 0; li < LA.size(); li++) {
-  //   LA[li] *= weights[li];
-  // }
 }
 
 /*