RgpotAdapter.h

Go to the documentation of this file.

/*
** This file is part of eOn.
**
** SPDX-License-Identifier: BSD-3-Clause
**
** Adapter over rgpot potential kernels. Migrated pair potentials live in
** librgpot (github.com/OmniPotentRPC/rgpot); eOn's config surface,
** PotRegistry, force counters, and caching stay on this side.
*/
#pragma once
 
#include <memory>
#include <stdexcept>
#include <type_traits>
#include <vector>
 
#include "eon/Parameters.h"
#include "eon/Potential.h"
#include "eon/potentials/PluginLoader.h"
 
#include "rgpot/ForceStructs.hpp"
#include "rgpot/pot_caps.hpp"

template <class RPot> class RgpotAdapter final : public eonc::Potential {
public:
  template <class Cfg>
    requires(!std::is_base_of_v<eonc::IPluginLoader, std::remove_cvref_t<Cfg>>)
  RgpotAdapter(eonc::PotType ptype, const eonc::Parameters &params,
               const Cfg &cfg)
      : RgpotAdapter(ptype, params, cfg, eonc::PluginLoader::instance()) {}

  template <class Cfg>
    requires(!std::is_base_of_v<eonc::IPluginLoader, std::remove_cvref_t<Cfg>>)
  RgpotAdapter(eonc::PotType ptype, const eonc::Parameters &params,
               const Cfg &cfg, eonc::IPluginLoader &loader)
      : eonc::Potential(ptype, params), pot_(cfg), loader_(loader) {
    finishConstruct(params);
  }

  RgpotAdapter(eonc::PotType ptype, const eonc::Parameters &params)
      : eonc::Potential(ptype, params), pot_(),
        loader_(eonc::PluginLoader::instance()) {
    finishConstruct(params);
  }

  RgpotAdapter(eonc::PotType ptype, const eonc::Parameters &params,
               eonc::IPluginLoader &loader)
      : eonc::Potential(ptype, params), pot_(), loader_(loader) {
    finishConstruct(params);
  }

  RgpotAdapter(eonc::PotType ptype, const eonc::Parameters &params,
               RPot &&kernel)
      : RgpotAdapter(ptype, params, std::move(kernel),
                     eonc::PluginLoader::instance()) {}

  RgpotAdapter(eonc::PotType ptype, const eonc::Parameters &params,
               RPot &&kernel, eonc::IPluginLoader &loader)
      : eonc::Potential(ptype, params), pot_(std::move(kernel)),
        loader_(loader) {
    finishConstruct(params);
  }
 
  void force(long N, const double *R, const int *atomicNrs, double *F,
             double *U, double *variance, const double *box) override {
    rgpot::ForceInput in{static_cast<size_t>(N), R, atomicNrs, box};
    rgpot::ForceOut out{F, 0.0, 0.0};
    pot_.forceImpl(in, &out);
    *U = out.energy;
    if (variance != nullptr) {
      *variance = out.variance;
    }
    storeStress(out);
  }
 
  [[nodiscard]] bool computesStress() const noexcept override {
    return haveStress_;
  }

  [[nodiscard]] Matrix3d cauchyStress() const override {
    if (!haveStress_) {
      throw std::logic_error("rgpot kernel did not return a stress tensor");
    }
    return stress_;
  }

  [[nodiscard]] bool supportsBatchEvaluation() const noexcept override {
    return pot_.caps().batched;
  }

  void forceBatch(long nSystems, long nAtoms, const double *const *positions,
                  const int *const *atomicNrs, double *const *forces,
                  double *energies, double *variances,
                  const double *const *boxes) override {
    const auto n = static_cast<size_t>(nSystems);
    std::vector<rgpot::ForceInput> in;
    std::vector<rgpot::ForceOut> out;
    in.reserve(n);
    out.reserve(n);
    for (size_t i = 0; i < n; ++i) {
      in.push_back(rgpot::ForceInput{static_cast<size_t>(nAtoms), positions[i],
                                     atomicNrs[i], boxes[i]});
      out.push_back(rgpot::ForceOut{forces[i], 0.0, 0.0});
    }
 
    pot_.forceBatchImpl(
        rgpot::ForceBatch{.nSystems = n, .in = in.data(), .out = out.data()});
 
    haveStress_ = false;
    for (size_t i = 0; i < n; ++i) {
      energies[i] = out[i].energy;
      if (variances != nullptr) {
        variances[i] = out[i].variance;
      }
      if (out[i].has_stress) {
        storeStress(out[i]);
      }
      forceCallCounter++;
      eonc::PotRegistry::get().on_force_call(ptype);
    }
  }
 
  [[nodiscard]] bool isThreadSafe() const noexcept override {
    return pot_.caps().reentrancy == rgpot::Reentrancy::SharedInstance;
  }

  [[nodiscard]] bool requiresIsolatedMoleculeLayout() const noexcept override {
    return !pot_.caps().periodic;
  }
 
  [[nodiscard]] bool needsPerImageInstance() const noexcept override {
    const auto caps = pot_.caps();
    return caps.perImageInstances ||
           caps.reentrancy == rgpot::Reentrancy::PerInstance;
  }
 
  [[nodiscard]] const RPot &kernel() const noexcept { return pot_; }

  [[nodiscard]] double finiteCutoff() const noexcept override {
    if constexpr (requires(const RPot &k) { k.config().cutoff; }) {
      return pot_.config().cutoff;
    } else {
      return 0.0;
    }
  }
 
private:
  void finishConstruct(const eonc::Parameters &params) {
    loader_.add_config_paths(params.potential_options().potentialsPath);
    validateCaps();
  }
 
  void validateCaps() const {
    const auto caps = pot_.caps();
    using R = rgpot::Reentrancy;
    if (caps.reentrancy != R::SharedInstance &&
        caps.reentrancy != R::PerInstance &&
        caps.reentrancy != R::ProcessSerial) {
      throw std::runtime_error("rgpot kernel advertised an unknown reentrancy");
    }
  }
 
  void storeStress(const rgpot::ForceOut &out) {
    haveStress_ = out.has_stress != 0;
    if (!haveStress_) {
      return;
    }
    for (int row = 0; row < 3; ++row) {
      for (int col = 0; col < 3; ++col) {
        stress_(row, col) = out.stress[row * 3 + col];
      }
    }
  }
 
  RPot pot_;
  eonc::IPluginLoader &loader_;
  Matrix3d stress_{Matrix3d::Zero()};
  bool haveStress_{false};
};

template <class RPot>
std::shared_ptr<eonc::Potential>
makeRgpotDefault(eonc::PotType ptype, const eonc::Parameters &params) {
  return std::make_shared<RgpotAdapter<RPot>>(ptype, params);
}

template <class RPot, class Cfg>
std::shared_ptr<eonc::Potential>
makeRgpot(eonc::PotType ptype, const eonc::Parameters &params, const Cfg &cfg) {
  return std::make_shared<RgpotAdapter<RPot>>(ptype, params, cfg);
}