Add reciprocal ewald force

src/energy.cpp

@@ -375,6 +375,42 @@ double Ewald::energy(Change &change) {
     return u;
 }
 
+/**
+ * @param forces Destination force vector
+ *
+ * Calculate forces from reciprocal space. Note that
+ * the destination force vector will *not* be zeroed
+ * before addition.
+ */
+void Ewald::force(std::vector<Point> &forces) {
+    assert(forces.size() == spc.p.size());
+    const double volume = spc.geo.getVolume();
+
+    // Surface contribution
+    Point total_dipole_moment = {0.0, 0.0, 0.0};
+    for (auto &particle : spc.p) {
+        auto mu = particle.hasExtension() ? particle.getExt().mu * particle.getExt().mulen : Point(0, 0, 0);
+        total_dipole_moment += particle.pos * particle.charge + mu;
+    }
+
+    auto force = forces.begin(); // iterator to force vector on first particle
+
+    for (auto &particle : spc.p) { // loop over particles
+        (*force) = total_dipole_moment * particle.charge / (2.0 * data.surface_dielectric_constant + 1.0);
+        double mu_scalar = particle.hasExtension() ? particle.getExt().mulen : 0.0;
+        std::complex<double> qmu(mu_scalar, particle.charge);
+        for (size_t i = 0; i < data.k_vectors.cols(); i++) { // loop over k vectors
+            std::complex<double> Q = data.Q_ion[i] + data.Q_dipole[i];
+            double k_dot_r = data.k_vectors.col(i).dot(particle.pos);
+            std::complex<double> expKri(std::cos(k_dot_r), std::sin(k_dot_r));
+            std::complex<double> repart = expKri * qmu * std::conj(Q);
+            (*force) += std::real(repart) * data.k_vectors.col(i) * data.Aks[i];
+        }
+        (*force) *= -4.0 * pc::pi / volume * data.bjerrum_length; // to units of kT/Angstrom^2
+        force++;                                                  // advance to next force vector
+    }
+}
+
 /**
  * @todo Implement a sync() function in EwaldData to selectively copy information
  */

src/energy.h

@@ -205,6 +205,7 @@ class Ewald : public Energybase {
               Change &) override; //!< Called after a move is rejected/accepted
                                   //! as well as before simulation
     void to_json(json &) const override;
+    void force(std::vector<Point> &) override; // update forces on all particles
 };
 
 class Isobaric : public Energybase {