🎨 added a new function get_chargeless_potential_between_sidbs

include/fiction/technology/charge_distribution_surface.hpp

@@ -334,7 +334,7 @@ class charge_distribution_surface<Lyt, false> : public Lyt
     /**
      * Returns the charge state of a given cell.
      *
-     * @param c cell.
+     * @param c The cell.
      * @return The charge state of the given cell.
      */
     [[nodiscard]] sidb_charge_state get_charge_state(const typename Lyt::cell& c) const noexcept
@@ -363,10 +363,10 @@ class charge_distribution_surface<Lyt, false> : public Lyt
         return -1;
     }
     /**
-     *  Returns the distance between two cells
+     *  Returns the distance between two cells.
      *
-     *  @param c1 the first cell to compare
-     *  @param c2 the second cell to compare
+     *  @param c1 the first cell to compare.
+     *  @param c2 the second cell to compare.
      *  @return a constexpr double representing the distance between the two cells.
      */
     [[nodiscard]] double get_distance_between_cells(const typename Lyt::cell& c1,
@@ -391,27 +391,51 @@ class charge_distribution_surface<Lyt, false> : public Lyt
         return strg->dist_mat[index1][index2];
     }
     /**
-     * Calculates and returns the electrostatic potential between two cells.
+     * Returns the chargeless electrostatic potential between two cells.
      *
-     * @param c1 The first cell
-     * @param c2 The second cell
-     * @return The electrostatic potential between `c1` and `c2`.
+     * @note If the signed electrostatic potential \f$ V_{i,j} \f$ is required, use the `get_potential_between_sidbs`
+     * function.
+     *
+     * @param c1 The first cell.
+     * @param c2 The second cell.
+     * @return The chargeless electrostatic potential between `c1` and `c2`, i.e, \f$ \frac{V_{i,j}}{n_j} \f$.
+     */
+    [[nodiscard]] double get_chargeless_potential_between_sidbs(const typename Lyt::cell& c1,
+                                                                const typename Lyt::cell& c2) const noexcept
+    {
+        if (const auto index1 = cell_to_index(c1), index2 = cell_to_index(c2); (index1 != -1) && (index2 != -1))
+        {
+            return strg->pot_mat[static_cast<uint64_t>(index1)][static_cast<uint64_t>(index2)];
+        }
+
+        return 0;
+    }
+    /**
+     * Calculates and returns the electrostatic potential at one cell (`c1`) generated by another cell (`c2`).
+     *
+     * @note If the chargeless electrostatic potential \f$ \frac{V_{i,j}}{n_j} \f$ is required, use the
+     * `get_chargeless_potential_between_sidbs` function.
+     *
+     * @param c1 The first cell.
+     * @param c2 The second cell.
+     * @return The electrostatic potential between `c1` and `c2`, i.e., \f$ V_{i,j} \f$.
      */
     [[nodiscard]] double get_potential_between_sidbs(const typename Lyt::cell& c1,
                                                      const typename Lyt::cell& c2) const noexcept
     {
         if (const auto index1 = cell_to_index(c1), index2 = cell_to_index(c2); (index1 != -1) && (index2 != -1))
         {
-            return strg->pot_mat[static_cast<uint64_t>(index1)][static_cast<uint64_t>(index2)];
+            return strg->pot_mat[static_cast<uint64_t>(index1)][static_cast<uint64_t>(index2)] *
+                   charge_state_to_sign(get_charge_state(c2));
         }
 
         return 0;
     }
     /**
      * Calculates and returns the potential of two indices.
      *
-     * @param index1 The first index
-     * @param index2 The second index
+     * @param index1 The first index.
+     * @param index2 The second index.
      * @return The potential between `index1` and `index2`.
      */
     [[nodiscard]] double get_electrostatic_potential_by_indices(const uint64_t index1,
@@ -450,14 +474,7 @@ class charge_distribution_surface<Lyt, false> : public Lyt
         const auto index1 = static_cast<std::size_t>(cell_to_index(c1));
         const auto index2 = static_cast<std::size_t>(cell_to_index(c2));
 
-        if (strg->dist_mat[index1][index2] == 0)
-        {
-            return 0.0;
-        }
-
-        return (strg->phys_params.k / strg->dist_mat[index1][index2] *
-                std::exp(-strg->dist_mat[index1][index2] / strg->phys_params.lambda_tf) *
-                physical_constants::ELECTRIC_CHARGE);
+        return potential_between_sidbs_by_index(index1, index2);
     }
     /**
      * The function calculates the electrostatic potential for each SiDB position (local).

test/technology/charge_distribution_surface.cpp

@@ -3,6 +3,7 @@
 //
 
 #include <catch2/catch_template_test_macros.hpp>
+#include <catch2/matchers/catch_matchers_floating_point.hpp>
 
 #include <fiction/layouts/cartesian_layout.hpp>
 #include <fiction/layouts/cell_level_layout.hpp>
@@ -143,18 +144,18 @@ TEMPLATE_TEST_CASE(
         // all SiDBs' charge states are set to neutral
         charge_layout.set_all_charge_states(sidb_charge_state::NEUTRAL);
         //
-        //        // read SiDBs' charge states
+        // read SiDBs' charge states
         CHECK(charge_layout.get_charge_state({5, 4}) == sidb_charge_state::NEUTRAL);
         CHECK(charge_layout.get_charge_state({5, 5}) == sidb_charge_state::NEUTRAL);
         CHECK(charge_layout.get_charge_state({5, 6}) == sidb_charge_state::NEUTRAL);
         CHECK(charge_layout.get_charge_state({5, 1}) == sidb_charge_state::NONE);
         //
         charge_layout.set_all_charge_states(sidb_charge_state::NEUTRAL);
         //
-        //        // all SiDBs' charge states are set to negative
+        // all SiDBs' charge states are set to negative
         charge_layout.set_all_charge_states(sidb_charge_state::NEGATIVE);
         //
-        //        // read SiDBs' charge states
+        // read SiDBs' charge states
         CHECK(charge_layout.get_charge_state({5, 4}) == sidb_charge_state::NEGATIVE);
         CHECK(charge_layout.get_charge_state({5, 5}) == sidb_charge_state::NEGATIVE);
         CHECK(charge_layout.get_charge_state({5, 6}) == sidb_charge_state::NEGATIVE);
@@ -207,15 +208,15 @@ TEMPLATE_TEST_CASE(
         lyt.assign_cell_type({1, 10, 1}, TestType::cell_type::NORMAL);
         charge_distribution_surface charge_layout{lyt, sidb_simulation_parameters{}};
 
-        CHECK(charge_layout.get_potential_between_sidbs({0, 0, 0}, {0, 0, 0}) == 0.0);
-        CHECK(charge_layout.get_potential_between_sidbs({1, 8, 0}, {1, 8, 0}) == 0.0);
-        CHECK(charge_layout.get_potential_between_sidbs({1, 10, 1}, {1, 10, 1}) == 0.0);
-        CHECK((charge_layout.get_potential_between_sidbs({1, 8, 0}, {0, 0, 0}) - 0.0121934043) < 0.00000001);
-        CHECK(std::abs(charge_layout.get_potential_between_sidbs({0, 0, 0}, {1, 10, 1}) -
-                       charge_layout.get_potential_between_sidbs({1, 10, 1}, {0, 0, 0})) <
+        CHECK(charge_layout.get_chargeless_potential_between_sidbs({0, 0, 0}, {0, 0, 0}) == 0.0);
+        CHECK(charge_layout.get_chargeless_potential_between_sidbs({1, 8, 0}, {1, 8, 0}) == 0.0);
+        CHECK(charge_layout.get_chargeless_potential_between_sidbs({1, 10, 1}, {1, 10, 1}) == 0.0);
+        CHECK((charge_layout.get_chargeless_potential_between_sidbs({1, 8, 0}, {0, 0, 0}) - 0.0121934043) < 0.00000001);
+        CHECK(std::abs(charge_layout.get_chargeless_potential_between_sidbs({0, 0, 0}, {1, 10, 1}) -
+                       charge_layout.get_chargeless_potential_between_sidbs({1, 10, 1}, {0, 0, 0})) <
               physical_constants::POP_STABILITY_ERR);
-        CHECK(charge_layout.get_potential_between_sidbs({0, 0, 0}, {1, 8, 0}) >
-              charge_layout.get_potential_between_sidbs({1, 10, 1}, {0, 0, 0}));
+        CHECK(charge_layout.get_chargeless_potential_between_sidbs({0, 0, 0}, {1, 8, 0}) >
+              charge_layout.get_chargeless_potential_between_sidbs({1, 10, 1}, {0, 0, 0}));
     }
     //
     SECTION("Local Potential")
@@ -372,4 +373,41 @@ TEMPLATE_TEST_CASE(
         charge_layout_new.increase_charge_index_by_one();
         CHECK(charge_layout_new.get_charge_index().first == 15);
     }
+
+    SECTION("using chargeless and normal potential function")
+    {
+        TestType                         lyt_new{{11, 11}};
+        const sidb_simulation_parameters params{3, -0.32, 5.0 * 1E-9, 3.84 * 1E-10, 7.68 * 1E-10, 2.25 * 1E-10};
+
+        lyt_new.assign_cell_type({0, 0, 1}, TestType::cell_type::NORMAL);
+        lyt_new.assign_cell_type({1, 3, 0}, TestType::cell_type::NORMAL);
+        lyt_new.assign_cell_type({10, 5, 1}, TestType::cell_type::NORMAL);
+
+        charge_distribution_surface charge_layout_new{lyt_new, params};
+
+        charge_layout_new.assign_charge_state({0, 0, 1}, sidb_charge_state::NEGATIVE);
+        charge_layout_new.assign_charge_state({1, 3, 0}, sidb_charge_state::POSITIVE);
+        charge_layout_new.assign_charge_state({10, 5, 1}, sidb_charge_state::NEUTRAL);
+
+        CHECK(charge_layout_new.get_charge_state({0, 0, 1}) == sidb_charge_state::NEGATIVE);
+        CHECK(charge_layout_new.get_charge_state({1, 3, 0}) == sidb_charge_state::POSITIVE);
+        CHECK(charge_layout_new.get_charge_state({10, 5, 1}) == sidb_charge_state::NEUTRAL);
+
+        CHECK(charge_layout_new.get_chargeless_potential_between_sidbs({0, 0, 1}, {1, 3, 0}) > 0);
+        CHECK(charge_layout_new.get_potential_between_sidbs({0, 0, 1}, {1, 3, 0}) > 0);
+
+        CHECK(charge_layout_new.get_chargeless_potential_between_sidbs({0, 0, 1}, {10, 5, 1}) > 0);
+        CHECK(charge_layout_new.get_potential_between_sidbs({0, 0, 1}, {10, 5, 1}) == 0.0);
+
+        CHECK(charge_layout_new.get_chargeless_potential_between_sidbs({10, 5, 1}, {0, 0, 1}) > 0);
+        CHECK(charge_layout_new.get_potential_between_sidbs({10, 5, 1}, {0, 0, 1}) < 0);
+
+        CHECK_THAT(charge_layout_new.get_potential_between_sidbs({10, 5, 1}, {0, 0, 1}) +
+                       charge_layout_new.get_chargeless_potential_between_sidbs({10, 5, 1}, {0, 0, 1}),
+                   Catch::Matchers::WithinAbs(0.000000, 0.000001));
+
+        CHECK_THAT(charge_layout_new.get_potential_between_sidbs({0, 0, 1}, {1, 3, 0}) -
+                       charge_layout_new.get_chargeless_potential_between_sidbs({0, 0, 1}, {1, 3, 0}),
+                   Catch::Matchers::WithinAbs(0.000000, 0.000001));
+    }
 }