Potential.h

Go to the documentation of this file.

/*
** This file is part of eOn.
**
** SPDX-License-Identifier: BSD-3-Clause
**
** Copyright (c) 2010--present, eOn Development Team
** All rights reserved.
**
** Repo:
** https://github.com/TheochemUI/eOn
*/
#pragma once
#include "EonLogger.h"
 
#include "Eigen.h"
#include "PotRegistry.h"
#include <atomic>
#include <memory>
#include <span>
#include <stdexcept>
 
namespace eonc {
 
class Parameters;
class Runtime;

class PotentialConstructionScope {
public:
  explicit PotentialConstructionScope(IPotRegistry &registry);
  ~PotentialConstructionScope();
  PotentialConstructionScope(const PotentialConstructionScope &) = delete;
  PotentialConstructionScope &
  operator=(const PotentialConstructionScope &) = delete;
 
private:
  IPotRegistry *prev_;
};
 
class Potential {
protected:
  PotType ptype;
 
private:
  IPotRegistry &registry_;
  uint64_t m_registry_id;
  IPotRegistry::TimePoint m_created_at;
  bool force_serial_{false};
 
public:
  std::atomic<size_t> forceCallCounter;

  explicit Potential(PotType a_ptype);

  Potential(PotType a_ptype, IPotRegistry &registry)
      : ptype{a_ptype}, registry_{registry},
        m_registry_id{registry.on_created(a_ptype)},
        m_created_at{IPotRegistry::Clock::now()}, forceCallCounter{0} {}
 
  // Out-of-line in eoncbase (PotentialParams.cpp) so shared plugins
  // that link only eoncbase get the symbols. This header only
  // forward-declares Parameters (header-surface contract).
  Potential(PotType a_ptype, const Parameters &p);
  explicit Potential(const Parameters &a_params);
 
  virtual ~Potential() {
    registry_.on_destroyed(m_registry_id, ptype, forceCallCounter,
                           m_created_at);
  }
 
  // Does not take into account the fixed / free atoms
  // Variance here is null when not needed and that's OK
  void virtual force(long nAtoms, const double *positions, const int *atomicNrs,
                     double *forces, double *energy, double *variance,
                     const double *box) = 0;

  void force(std::span<const double> positions, std::span<const int> atomicNrs,
             std::span<double> forces, double *energy, double *variance,
             std::span<const double> box) {
    if (positions.size() % 3 != 0 || positions.size() / 3 != atomicNrs.size() ||
        forces.size() != positions.size() || box.size() != 9) {
      throw std::invalid_argument("Potential::force span size mismatch");
    }
    force(static_cast<long>(atomicNrs.size()), positions.data(),
          atomicNrs.data(), forces.data(), energy, variance, box.data());
  }

  virtual void setFixedMask(long nAtoms, const double *isFixed) {
    (void)nAtoms;
    (void)isFixed;
  }
 
  std::tuple<double, AtomMatrix>
  get_ef(const AtomMatrix &pos, const VectorXi &atmnrs, const Matrix3d &box);
 
  [[nodiscard]] PotType getType() const { return this->ptype; }

  [[nodiscard]] virtual double finiteCutoff() const noexcept { return 0.0; }

  [[nodiscard]] virtual bool isSurrogate() const noexcept { return false; }

  [[nodiscard]] virtual bool requiresIsolatedMoleculeLayout() const noexcept {
    return false;
  }

  [[nodiscard]] virtual bool isThreadSafe() const noexcept { return true; }

  enum class PotLayout : unsigned {
    InProcess = 1u << 0,
    NeedsWorkingDirectory = 1u << 1,
    Subprocess = 1u << 2,
  };
  [[nodiscard]] virtual unsigned layoutFlags() const noexcept {
    return static_cast<unsigned>(PotLayout::InProcess);
  }

  [[nodiscard]] virtual bool isSharedInstanceThreadSafe() const noexcept {
    if (force_serial_)
      return false;
    return isThreadSafe();
  }

  [[nodiscard]] virtual bool needsPerImageInstance() const noexcept {
    return false;
  }

  [[nodiscard]] virtual std::shared_ptr<Potential> clonePotential() const {
    return nullptr;
  }

  [[nodiscard]] virtual bool supportsBatchEvaluation() const noexcept {
    return false;
  }

  [[nodiscard]] virtual bool computesStress() const noexcept { return false; }

  [[nodiscard]] virtual Matrix3d cauchyStress() const {
    throw std::logic_error("Potential::cauchyStress is not implemented");
  }

  virtual void forceBatchOwned(long nSystems, long nAtoms,
                               const double *const *positions,
                               const int *const *atomicNrs,
                               double *const *forces, double *energies,
                               double *variances, const double *const *boxes,
                               const long *owners) {
    (void)owners;
    forceBatch(nSystems, nAtoms, positions, atomicNrs, forces, energies,
               variances, boxes);
  }
 
  virtual 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) {
    for (long i = 0; i < nSystems; i++) {
      double var = 0;
      force(nAtoms, positions[i], atomicNrs[i], forces[i], &energies[i], &var,
            boxes[i]);
      if (variances)
        variances[i] = var;
      forceCallCounter++;
      registry_.on_force_call(ptype);
    }
  }
};
 
namespace helpers {
std::shared_ptr<Potential> makePotential(const Parameters &params);
std::shared_ptr<Potential> makePotential(PotType ptype,
                                         const Parameters &params);
std::shared_ptr<Potential>
makePotential(PotType ptype, const Parameters &params, Runtime &runtime);
} // namespace helpers
 
} // namespace eonc