Add batch calculations for energies and forces

avogadro/calc/energycalculator.cpp

@@ -66,6 +66,28 @@ void EnergyCalculator::gradient(const Eigen::VectorXd& x, Eigen::VectorXd& grad)
   constraintGradients(x, grad);
 }
 
+std::vector<Real> EnergyCalculator::valueBatch(
+  const std::vector<Eigen::VectorXd>& coords)
+{
+  // Default: loop one coordinate set at a time. Derived classes with a true
+  // batch transport (e.g. ScriptEnergy over binary-v1) override this.
+  std::vector<Real> energies;
+  energies.reserve(coords.size());
+  for (const auto& x : coords)
+    energies.push_back(value(x));
+  return energies;
+}
+
+void EnergyCalculator::gradientBatch(const std::vector<Eigen::VectorXd>& coords,
+                                     std::vector<Eigen::VectorXd>& grads)
+{
+  // Default: loop one coordinate set at a time. gradient() already cleans and
+  // applies constraints, so the fallback matches the single-shot path.
+  grads.resize(coords.size());
+  for (std::size_t i = 0; i < coords.size(); ++i)
+    gradient(coords[i], grads[i]);
+}
+
 void EnergyCalculator::hessian(const Eigen::VectorXd& x, Eigen::MatrixXd& hess)
 {
   finiteHessian(x, hess);

avogadro/calc/energycalculator.h

@@ -8,6 +8,7 @@
 
 #include "avogadrocalcexport.h"
 
+#include <avogadro/core/avogadrocore.h>
 #include <avogadro/core/constraint.h>
 #include <avogadro/core/molecule.h>
 #include <avogadro/core/variantmap.h>
@@ -20,7 +21,8 @@ class Molecule;
 
 namespace Calc {
 
-constexpr Real KCAL_TO_KJ = 4.184;
+// Alias the shared constant in avogadrocore.h so the value lives in one place.
+constexpr Real KCAL_TO_KJ = Avogadro::KCAL_TO_KJ;
 
 class AVOGADROCALC_EXPORT EnergyCalculator
 {
@@ -117,6 +119,41 @@ class AVOGADROCALC_EXPORT EnergyCalculator
   virtual Real evaluate(const Eigen::VectorXd& x,
                         Eigen::VectorXd* grad = nullptr);
 
+  /**
+   * @brief Indicate if this method can evaluate many coordinate sets in a
+   * single call (e.g. an external script that batches frames).
+   *
+   * When false (the default), valueBatch()/gradientBatch() simply loop over
+   * the supplied coordinate sets, calling value()/gradient() once each.
+   */
+  virtual bool supportsBatch() const { return false; }
+
+  /**
+   * Evaluate the energy for several coordinate sets (e.g. conformers or
+   * trajectory frames).
+   *
+   * The default implementation loops over @p coords calling value(). Derived
+   * classes with a true batch transport should override this.
+   *
+   * @param coords A vector of coordinate vectors, each of length 3 * atomCount.
+   * @return One energy per coordinate set, in the same order as @p coords.
+   */
+  virtual std::vector<Real> valueBatch(
+    const std::vector<Eigen::VectorXd>& coords);
+
+  /**
+   * Evaluate the gradients for several coordinate sets.
+   *
+   * The default implementation loops over @p coords calling gradient().
+   * Derived classes with a true batch transport should override this.
+   *
+   * @param coords A vector of coordinate vectors, each of length 3 * atomCount.
+   * @param grads Filled with one gradient vector per coordinate set, in the
+   * same order as @p coords (resized as needed).
+   */
+  virtual void gradientBatch(const std::vector<Eigen::VectorXd>& coords,
+                             std::vector<Eigen::VectorXd>& grads);
+
   /**
    * Calculate the Hessian matrix for this method, defaulting to numerical
    * finite-difference methods.

avogadro/core/avogadrocore.h

@@ -91,6 +91,11 @@ constexpr Real ANGSTROM_TO_BOHR = static_cast<Real>(ANGSTROM_TO_BOHR_D);
 constexpr double HARTREE_TO_EV_D = 27.211386245981;
 constexpr float HARTREE_TO_EV_F = static_cast<float>(HARTREE_TO_EV_D);
 constexpr Real HARTREE_TO_EV = static_cast<Real>(HARTREE_TO_EV_D);
+
+// thermochemical calorie
+constexpr double KCAL_TO_KJ_D = 4.184;
+constexpr float KCAL_TO_KJ_F = static_cast<float>(KCAL_TO_KJ_D);
+constexpr Real KCAL_TO_KJ = static_cast<Real>(KCAL_TO_KJ_D);
 /** @} */
 
 } // namespace Avogadro

avogadro/core/molecule.cpp

@@ -1442,15 +1442,22 @@ bool Molecule::setCoordinate3d(int coord)
 {
   if (coord >= 0 && coord < static_cast<int>(m_coordinates3d.size())) {
     m_positions3d = m_coordinates3d[coord];
+    m_coordinate3dIndex = coord;
     return true;
   }
   return false;
 }
 
+int Molecule::coordinate3d() const
+{
+  return m_coordinate3dIndex;
+}
+
 void Molecule::clearCoordinate3d()
 {
   m_coordinates3d.clear();
   m_conformerProperties.clear();
+  m_coordinate3dIndex = 0;
 }
 
 Array<Vector3> Molecule::coordinate3d(size_t index) const

avogadro/core/molecule.h

@@ -682,6 +682,8 @@ class AVOGADROCORE_EXPORT Molecule
 
   size_t coordinate3dCount() const;
   bool setCoordinate3d(int coord);
+  /** @return the index of the currently active coordinate set. */
+  int coordinate3d() const;
   Array<Vector3> coordinate3d(size_t index) const;
   bool setCoordinate3d(const Array<Vector3>& coords, size_t index);
 
@@ -1004,6 +1006,7 @@ class AVOGADROCORE_EXPORT Molecule
   Array<std::string> m_bondLabels;
   Array<std::string> m_residueLabels;
   Array<Array<Vector3>> m_coordinates3d; //!< Store conformers/trajectories.
+  int m_coordinate3dIndex = 0;           //!< Active coordinate set index.
   Array<Array<Vector3>> m_velocities;    //!< Store velocities.
   Array<double> m_timesteps;
   Array<AtomHybridization> m_hybridizations;

avogadro/qtgui/calcworker.cpp

@@ -8,6 +8,9 @@
 #include <avogadro/calc/energycalculator.h>
 #include <avogadro/calc/energyoptimizer.h>
 
+#include <algorithm>
+#include <iterator>
+
 namespace Avogadro::QtGui {
 
 CalcWorker::CalcWorker(QObject* parent) : QObject(parent)
@@ -20,6 +23,9 @@ CalcWorker::CalcWorker(QObject* parent) : QObject(parent)
     "std::vector<Avogadro::Core::Constraint>");
   qRegisterMetaType<Calc::EnergyCalculator*>(
     "Avogadro::Calc::EnergyCalculator*");
+  qRegisterMetaType<std::vector<Eigen::VectorXd>>(
+    "std::vector<Eigen::VectorXd>");
+  qRegisterMetaType<std::vector<double>>("std::vector<double>");
 }
 
 CalcWorker::~CalcWorker() = default;
@@ -103,6 +109,69 @@ void CalcWorker::runGradient(Eigen::VectorXd positions)
   emit evaluateFinished(gradient, energy);
 }
 
+void CalcWorker::runEvaluateBatch(std::vector<Eigen::VectorXd> coordsList,
+                                  bool computeGradient, int chunkSize)
+{
+  std::vector<double> energies;
+  std::vector<Eigen::VectorXd> gradients;
+
+  if (!m_calc || coordsList.empty() || m_cancelled) {
+    emit evaluateBatchFinished(energies, gradients);
+    return;
+  }
+
+  const int total = static_cast<int>(coordsList.size());
+
+  // Memory-bounded auto chunk size: cap the per-chunk payload (which also
+  // bounds the equally-sized gradient response) to a sane budget.
+  if (chunkSize <= 0) {
+    constexpr std::size_t byteBudget = 32 * 1024 * 1024; // ~32 MB
+    constexpr std::size_t maxFramesCap = 256;
+    const std::size_t frameDoubles =
+      static_cast<std::size_t>(coordsList.front().size());
+    const std::size_t frameBytes =
+      std::max<std::size_t>(1, frameDoubles * sizeof(double));
+    std::size_t frames = byteBudget / frameBytes;
+    frames = std::clamp<std::size_t>(frames, 1, maxFramesCap);
+    chunkSize = static_cast<int>(frames);
+  }
+
+  energies.reserve(coordsList.size());
+  if (computeGradient)
+    gradients.reserve(coordsList.size());
+
+  for (std::size_t start = 0; start < coordsList.size();
+       start += static_cast<std::size_t>(chunkSize)) {
+    if (m_cancelled)
+      break;
+
+    const std::size_t end =
+      std::min(coordsList.size(), start + static_cast<std::size_t>(chunkSize));
+    std::vector<Eigen::VectorXd> chunk(coordsList.begin() + start,
+                                       coordsList.begin() + end);
+
+    if (computeGradient) {
+      std::vector<Eigen::VectorXd> chunkGrads;
+      m_calc->gradientBatch(chunk, chunkGrads);
+      // Energy alongside the gradient so callers always get both.
+      std::vector<double> chunkEnergies = m_calc->valueBatch(chunk);
+      energies.insert(energies.end(), chunkEnergies.begin(),
+                      chunkEnergies.end());
+      gradients.insert(gradients.end(),
+                       std::make_move_iterator(chunkGrads.begin()),
+                       std::make_move_iterator(chunkGrads.end()));
+    } else {
+      std::vector<double> chunkEnergies = m_calc->valueBatch(chunk);
+      energies.insert(energies.end(), chunkEnergies.begin(),
+                      chunkEnergies.end());
+    }
+
+    emit batchProgress(static_cast<int>(end), total);
+  }
+
+  emit evaluateBatchFinished(energies, gradients);
+}
+
 void CalcWorker::cancel()
 {
   m_cancelled = true;

avogadro/qtgui/calcworker.h

@@ -28,6 +28,8 @@ Q_DECLARE_METATYPE(Avogadro::Core::Molecule)
 Q_DECLARE_METATYPE(Avogadro::Calc::OptimizationOptions)
 Q_DECLARE_METATYPE(std::vector<Avogadro::Core::Constraint>)
 Q_DECLARE_METATYPE(Avogadro::Calc::EnergyCalculator*)
+Q_DECLARE_METATYPE(std::vector<Eigen::VectorXd>)
+Q_DECLARE_METATYPE(std::vector<double>)
 
 namespace Avogadro {
 
@@ -87,6 +89,23 @@ class AVOGADROQTGUI_EXPORT CalcWorker : public QObject
    */
   void evaluateFinished(Eigen::VectorXd gradient, double energy);
 
+  /**
+   * Emitted after each chunk of a batch evaluation completes.
+   * @param done Number of coordinate sets evaluated so far.
+   * @param total Total number of coordinate sets in the batch.
+   */
+  void batchProgress(int done, int total);
+
+  /**
+   * Emitted after a batch energy/gradient evaluation completes (or is
+   * cancelled, in which case the vectors hold the partial results computed
+   * before cancellation).
+   * @param energies One energy per coordinate set.
+   * @param gradients One gradient per coordinate set (empty if not requested).
+   */
+  void evaluateBatchFinished(std::vector<double> energies,
+                             std::vector<Eigen::VectorXd> gradients);
+
   /**
    * Emitted when the calculator has been initialized on the worker thread.
    */
@@ -122,6 +141,20 @@ public slots:
    */
   void runGradient(Eigen::VectorXd positions);
 
+  /**
+   * Evaluate energy (and optionally gradient) for many coordinate sets.
+   *
+   * The batch is split into chunks of @p chunkSize frames so a large
+   * trajectory never builds one giant wire buffer / model batch and so
+   * progress can be reported and cancellation honored between chunks.
+   *
+   * @param coordsList One coordinate vector per set (each 3 * atomCount long).
+   * @param computeGradient Request gradients in addition to energies.
+   * @param chunkSize Frames per chunk; <= 0 selects a memory-bounded default.
+   */
+  void runEvaluateBatch(std::vector<Eigen::VectorXd> coordsList,
+                        bool computeGradient, int chunkSize);
+
   /**
    * Request cancellation of the current computation.
    * Thread-safe: can be called from any thread.