charge-charge interaction

src/scf/coulombs_law.cpp

@@ -0,0 +1,81 @@
+/*
+ * Copyright 2024 NWChemEx-Project
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+#include "scf_modules.hpp"
+#include <cmath>
+#include <simde/simde.hpp>
+
+namespace scf {
+
+using pt = simde::charge_charge_interaction;
+
+const auto desc = R"(
+Charge-Charge Interaction via Coulomb's Law
+-------------------------------------------
+
+Let :math:`Q` be a set of :math:`N_Q` point charges such that the :math:`i`-th
+point charge has charge :math:`q_i` and position :math:`\mathbf{r_i}`. The
+electrostatic potential at the point :math:`\mathbf{r}` generated by :math:`Q`,
+:math:`E_Q(\mathbf{r})` is given by:
+
+.. math::
+
+   E_Q(\mathbf{r}) = \sum_{i=1}^{N_Q} 
+                     \frac{q_i}{\left|\mathbf{r} - \mathbf{r_i}\right|}
+
+Let :math:`S` be a set of :math:`N_S` point charges that is disjoint with 
+:math:`Q`, then the interaction of :math:`S` with :math:`Q`, :math:`V_{SQ}` is 
+given by:
+
+.. math::
+
+   V_{SQ} =& \sum_{i=1}^{N_S} q_i E_Q(\mathbf{r_i})\\
+          =& \sum_{i=1}^{N_S}\sum_{j=1}^{N_Q} 
+             \frac{q_i q_j}{\left|\mathbf{r_i} - \mathbf{r_j}\right|}
+
+This module will compute :math:`V_{SQ}`, including in the common scenario where
+:math:`S` equals :math:`Q` (in which case terms for which the denominator is
+zero are skipped).
+)";
+
+MODULE_CTOR(CoulombsLaw) {
+    description(desc);
+    satisfies_property_type<pt>();
+}
+
+MODULE_RUN(CoulombsLaw) {
+    const auto& [qlhs, qrhs] = pt::unwrap_inputs(inputs);
+
+    // This algorithm only works as long as qlhs and qrhs are disjoint or
+    // identical
+    double e = 0.0;
+    for(const auto&& qi : qlhs) {
+        for(const auto&& qj : qrhs) {
+            if(qi == qj) break;
+            auto dx = qi.x() - qj.x();
+            auto dy = qi.y() - qj.y();
+            auto dz = qi.z() - qj.z();
+
+            e += qi.charge() * qj.charge() /
+                 std::sqrt(dx * dx + dy * dy + dz * dz);
+        }
+    }
+
+    auto rv = results();
+    return pt::wrap_results(rv, e);
+}
+
+} // namespace scf

src/scf/scf_mm.cpp

@@ -22,6 +22,7 @@ namespace scf {
 
 void load_modules(pluginplay::ModuleManager& mm) {
     fock_operator::load_modules(mm);
+    mm.add_module<CoulombsLaw>("Coulomb's Law");
 #ifdef BUILD_TAMM_SCF
     mm.add_module<TAMMEnergy>("SCF Energy via TAMM");
 #endif

src/scf/scf_modules.hpp

@@ -23,4 +23,6 @@ namespace scf {
 DECLARE_MODULE(TAMMEnergy);
 #endif
 
+DECLARE_MODULE(CoulombsLaw);
+
 } // namespace scf

tests/cxx/unit_tests/coulombs_law.cpp

@@ -0,0 +1,59 @@
+/*
+ * Copyright 2024 NWChemEx-Project
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+#include "../test_scf.hpp"
+#include <catch2/matchers/catch_matchers_floating_point.hpp>
+#include <scf/scf.hpp>
+#include <simde/simde.hpp>
+
+using pt = simde::charge_charge_interaction;
+using Catch::Matchers::WithinAbs;
+
+TEST_CASE("CoulombsLaw") {
+    pluginplay::ModuleManager mm;
+    scf::load_modules(mm);
+    auto& mod = mm.at("Coulomb's Law");
+
+    using charge_type = typename simde::type::charges::value_type;
+    charge_type q0(1.0, 2.0, 3.0, 4.0);
+    charge_type q1(-1.0, 3.0, 4.0, 5.0);
+    charge_type q2(1.5, 4.0, 5.0, 6.0);
+
+    simde::type::charges empty;
+    simde::type::charges qs{q0, q1, q2};
+
+    SECTION("empty points") {
+        auto e = mod.run_as<pt>(empty, empty);
+        REQUIRE(e == 0.0);
+    }
+
+    SECTION("charges w/ itself") {
+        auto e = mod.run_as<pt>(qs, qs);
+        REQUIRE_THAT(e, WithinAbs(-1.0103629710818451, 1E-6));
+    }
+
+    SECTION("charges w/ empty") {
+        auto e = mod.run_as<pt>(qs, empty);
+        REQUIRE_THAT(e, WithinAbs(0.0, 1E-6));
+    }
+
+    SECTION("charges w/ different charges") {
+        simde::type::charges qs0{q0};
+        simde::type::charges qs12{q1, q2};
+        auto e = mod.run_as<pt>(qs0, qs12);
+        REQUIRE_THAT(e, WithinAbs(-0.1443375672974065, 1E-6));
+    }
+}