ParticlePropertiesLibrary.h

Go to the documentation of this file.

 
#pragma once
 
#include <algorithm>
#include <cmath>
#include <map>
#include <set>
#include <vector>
 
#include "autopas/utils/AlignedAllocator.h"
#include "autopas/utils/ExceptionHandler.h"
 
template <typename floatType = double, typename intType = unsigned long>
class ParticlePropertiesLibrary {
 public:
  explicit ParticlePropertiesLibrary(const double cutoff) : _cutoff(cutoff) {}
 
  ParticlePropertiesLibrary(const ParticlePropertiesLibrary &particlePropertiesLibrary) = default;
 
  ParticlePropertiesLibrary &operator=(const ParticlePropertiesLibrary &particlePropertiesLibrary) = default;
 
  void addSiteType(const intType siteId, const floatType mass);
 
  void addLJParametersToSite(const intType siteId, const floatType epsilon, const floatType sigma);
 
  void addATParametersToSite(const intType siteId, const floatType nu);
 
  void addMolType(const intType molId, const std::vector<intType> siteIds,
                  const std::vector<std::array<floatType, 3>> relPos, const std::array<floatType, 3> momentOfInertia);
 
  void calculateMixingCoefficients();
 
  ~ParticlePropertiesLibrary() = default;
 
  [[nodiscard]] int getNumberRegisteredSiteTypes() const { return _numRegisteredSiteTypes; }
 
  [[nodiscard]] int getNumberRegisteredMolTypes() const {
#if MD_FLEXIBLE_MODE == SINGLESITE
    AutoPasLog(
        WARN,
        "ParticlePropertiesLibrary::getNumberRegisteredMolTypes(): trying to get the number of registered multi-site"
        "molecule types when md-flexible has been compiled without support for multi-site molecules. Please compile "
        "with the CMake argument '-DMD_FLEXIBLE_MODE=MULTISITE'.");
#endif
    return _numRegisteredMolTypes;
  }
 
  floatType getEpsilon(intType i) const;
 
  floatType getSigma(intType i) const;
 
  floatType getNu(intType i) const;
 
  floatType getSiteMass(intType i) const;
 
  floatType getMolMass(intType i) const;
 
  std::array<floatType, 3> getMomentOfInertia(intType i) const;
 
  std::vector<std::array<floatType, 3>> getSitePositions(intType i) const;
 
  std::vector<intType> getSiteTypes(intType i) const;
 
  intType getNumSites(intType i) const;
 
  floatType getMoleculesLargestSigma(intType i) const;
 
  floatType getMixing24Epsilon(intType i, intType j) const {
    return _computedLJMixingData[i * _numRegisteredSiteTypes + j].epsilon24;
  }
 
  auto getLJMixingData(intType i, intType j) const { return _computedLJMixingData[i * _numRegisteredSiteTypes + j]; }
 
  const double *getLJMixingDataPtr(intType i, intType j) {
    return reinterpret_cast<const double *>(&_computedLJMixingData[i * _numRegisteredSiteTypes + j]);
  }
 
  floatType getMixingSigmaSquared(intType i, intType j) const {
    return _computedLJMixingData[i * _numRegisteredSiteTypes + j].sigmaSquared;
  }
 
  floatType getMixingShift6(intType i, intType j) const {
    return _computedLJMixingData[i * _numRegisteredSiteTypes + j].shift6;
  }
 
  static double calcShift6(double epsilon24, double sigmaSquared, double cutoffSquared);
 
  floatType getMixingNu(intType i, intType j, intType k) const {
    return _computedATMixingData[i * _numRegisteredSiteTypes * _numRegisteredSiteTypes + j * _numRegisteredSiteTypes +
                                 k]
        .nu;
  }
 
  auto getATMixingData(intType i, intType j, intType k) const {
    return _computedATMixingData[i * _numRegisteredSiteTypes * _numRegisteredSiteTypes + j * _numRegisteredSiteTypes +
                                 k];
  }
 
 private:
  intType _numRegisteredSiteTypes{0};
  intType _numRegisteredMolTypes{0};
  const double _cutoff;
 
  std::vector<floatType> _epsilons;
  std::vector<floatType> _sigmas;
  std::vector<floatType> _siteMasses;
  std::vector<floatType> _nus;  // Factor for AxilrodTeller potential
 
  // Note: this is a vector of site type Ids for the sites of a certain molecular Id
  std::vector<std::vector<intType>> _siteIds;
  // This is a vector (indexed by mol ID) of vectors of site positions (which are 3D arrays)
  std::vector<std::vector<std::array<floatType, 3>>> _relativeSitePositions;
  std::vector<floatType> _molMasses;
  std::vector<std::array<floatType, 3>> _momentOfInertias;
  std::vector<size_t> _numSites;
  std::vector<floatType> _moleculesLargestSigma;
 
  // Allocate memory for the respective parameters
  bool _storeLJData{false};
  bool _storeATData{false};
 
  struct PackedLJMixingData {
    floatType epsilon24;
    floatType sigmaSquared;
    floatType shift6;
  };
 
  struct PackedATMixingData {
    floatType nu;
  };
 
  std::vector<PackedLJMixingData, autopas::AlignedAllocator<PackedLJMixingData>> _computedLJMixingData;
  std::vector<PackedATMixingData, autopas::AlignedAllocator<PackedATMixingData>> _computedATMixingData;
};
 
template <typename floatType, typename intType>
void ParticlePropertiesLibrary<floatType, intType>::addSiteType(intType siteID, floatType mass) {
  if (_numRegisteredSiteTypes != siteID) {
    autopas::utils::ExceptionHandler::exception(
        "ParticlePropertiesLibrary::addSiteType(): trying to register a site type with id {}. Please "
        "register types "
        "consecutively, starting at id 0. Currently there are {} registered types.",
        siteID, _numRegisteredSiteTypes);
  }
  ++_numRegisteredSiteTypes;
  _siteMasses.emplace_back(mass);
 
  // Allocate memory for all parameters of used models
  if (_storeLJData) {
    _sigmas.emplace_back(0.0);
    _epsilons.emplace_back(0.0);
  }
  if (_storeATData) {
    _nus.emplace_back(0.0);
  }
}
 
template <typename floatType, typename intType>
void ParticlePropertiesLibrary<floatType, intType>::addLJParametersToSite(intType siteID, floatType epsilon,
                                                                          floatType sigma) {
  if (siteID >= _numRegisteredSiteTypes) {
    autopas::utils::ExceptionHandler::exception(
        "ParticlePropertiesLibrary::addLJParametersToSite(): Trying to set lennard-jones parameters for a site type "
        "with id {},"
        " which has not been registered yet. Currently there are {} registered types.",
        siteID, _numRegisteredSiteTypes);
  }
  _storeLJData = true;
  if (_epsilons.size() != _numRegisteredSiteTypes) {
    _epsilons.resize(_numRegisteredSiteTypes);
    _sigmas.resize(_numRegisteredSiteTypes);
  }
  _epsilons[siteID] = epsilon;
  _sigmas[siteID] = sigma;
}
 
template <typename floatType, typename intType>
void ParticlePropertiesLibrary<floatType, intType>::addATParametersToSite(intType siteID, floatType nu) {
  if (siteID >= _numRegisteredSiteTypes) {
    autopas::utils::ExceptionHandler::exception(
        "ParticlePropertiesLibrary::addATParametersToSite(): Trying to set the axilrod-teller parameter for a site "
        "type with id {},"
        " which has not been registered yet. Currently there are {} registered types.",
        siteID, _numRegisteredSiteTypes);
  }
  _storeATData = true;
  if (_nus.size() != _numRegisteredSiteTypes) {
    _nus.resize(_numRegisteredSiteTypes);
  }
  _nus[siteID] = nu;
}
 
template <typename floatType, typename intType>
void ParticlePropertiesLibrary<floatType, intType>::addMolType(const intType molId, const std::vector<intType> siteIds,
                                                               const std::vector<std::array<floatType, 3>> relPos,
                                                               const std::array<floatType, 3> momentOfInertia) {
  // Error handling
#if MD_FLEXIBLE_MODE == SINGLESITE
  AutoPasLog(WARN,
             "ParticlePropertiesLibrary::addMolType(): trying to register a multi-site molecule type when md-flexible "
             "has been compiled without support for multi-site molecules. Please compile with the CMake argument "
             "'-DMD_FLEXIBLE_MODE=MULTISITE'.");
#endif
  if (_numRegisteredMolTypes != molId) {
    autopas::utils::ExceptionHandler::exception(
        "ParticlePropertiesLibrary::addMolType(): trying to register a molecule type with id {}. Please register types "
        "consecutively, starting at id 0. Currently there are {} registered types.",
        molId, _numRegisteredSiteTypes);
  }
  if (std::any_of(siteIds.cbegin(), siteIds.cend(), [this](intType i) { return i >= this->_numRegisteredSiteTypes; })) {
    autopas::utils::ExceptionHandler::exception(
        "ParticlePropertiesLibrary::addMolType(): trying to register a molecule type with an unregistered site type "
        "Id.");
  }
  if (siteIds.size() != relPos.size()) {
    autopas::utils::ExceptionHandler::exception(
        "ParticlePropertiesLibrary::addMolType(): trying to register a molecule type with vectors of site IDs and site"
        "positions that do not match in size.");
  }
 
  // Add molecule type if there are no errors.
  ++_numRegisteredMolTypes;
  _siteIds.emplace_back(siteIds);
  _relativeSitePositions.emplace_back(relPos);
  _numSites.emplace_back(siteIds.size());
  floatType molMass{0.};
  for (intType site = 0; site < siteIds.size(); ++site) {
    molMass += _siteMasses[siteIds[site]];
  }
  _molMasses.emplace_back(molMass);
  _momentOfInertias.emplace_back(momentOfInertia);
 
  floatType molLargestSigma{0.};
  for (size_t site = 0; site < siteIds.size(); site++) {
    molLargestSigma = std::max(molLargestSigma, _sigmas[siteIds[site]]);
  }
  _moleculesLargestSigma.emplace_back(molLargestSigma);
}
 
template <typename floatType, typename intType>
void ParticlePropertiesLibrary<floatType, intType>::calculateMixingCoefficients() {
  if (_numRegisteredSiteTypes == 0) {
    autopas::utils::ExceptionHandler::AutoPasException(
        "ParticlePropertiesLibrary::calculateMixingCoefficients was called without any site types being registered!");
  }
 
  // There are Lennard-Jones Sites
  if (_storeLJData) {
    const auto cutoffSquared = _cutoff * _cutoff;
    _computedLJMixingData.resize(_numRegisteredSiteTypes * _numRegisteredSiteTypes);
 
    for (size_t firstIndex = 0ul; firstIndex < _numRegisteredSiteTypes; ++firstIndex) {
      for (size_t secondIndex = 0ul; secondIndex < _numRegisteredSiteTypes; ++secondIndex) {
        auto globalIndex = _numRegisteredSiteTypes * firstIndex + secondIndex;
 
        // epsilon
        const floatType epsilon24 = 24 * sqrt(_epsilons[firstIndex] * _epsilons[secondIndex]);
        _computedLJMixingData[globalIndex].epsilon24 = epsilon24;
 
        // sigma
        const floatType sigma = (_sigmas[firstIndex] + _sigmas[secondIndex]) / 2.0;
        const floatType sigmaSquared = sigma * sigma;
        _computedLJMixingData[globalIndex].sigmaSquared = sigmaSquared;
 
        // shift6
        const floatType shift6 = calcShift6(epsilon24, sigmaSquared, cutoffSquared);
        _computedLJMixingData[globalIndex].shift6 = shift6;
      }
    }
  }
 
  if (_storeATData) {
    _computedATMixingData.resize(_numRegisteredSiteTypes * _numRegisteredSiteTypes * _numRegisteredSiteTypes);
    for (size_t firstIndex = 0ul; firstIndex < _numRegisteredSiteTypes; ++firstIndex) {
      for (size_t secondIndex = 0ul; secondIndex < _numRegisteredSiteTypes; ++secondIndex) {
        for (size_t thirdIndex = 0ul; thirdIndex < _numRegisteredSiteTypes; ++thirdIndex) {
          const auto globalIndex3B = _numRegisteredSiteTypes * _numRegisteredSiteTypes * firstIndex +
                                     _numRegisteredSiteTypes * secondIndex + thirdIndex;
          // geometric mixing as used in e.g. https://doi.org/10.1063/1.3567308
          const floatType mixedNu = cbrt(_nus[firstIndex] * _nus[secondIndex] * _nus[thirdIndex]);
          _computedATMixingData[globalIndex3B].nu = mixedNu;
        }
      }
    }
  }
}
 
template <typename floatType, typename intType>
floatType ParticlePropertiesLibrary<floatType, intType>::getSiteMass(intType i) const {
  return _siteMasses[i];
}
 
template <typename floatType, typename intType>
floatType ParticlePropertiesLibrary<floatType, intType>::getMolMass(intType i) const {
#if MD_FLEXIBLE_MODE == MULTISITE
  return _molMasses[i];
#else
  return _siteMasses[i];
#endif
}
 
template <typename floatType, typename intType>
std::array<floatType, 3> ParticlePropertiesLibrary<floatType, intType>::getMomentOfInertia(intType i) const {
#if MD_FLEXIBLE_MODE == SINGLESITE
  AutoPasLog(WARN,
             "ParticlePropertiesLibrary::getMomentOfInertia(): trying to get the Moment of Inertia of a multi-site "
             "molecule type when md-flexible has been compiled without support for multi-site molecules. Please "
             "compile with the CMake argument '-DMD_FLEXIBLE_MODE=MULTISITE'.");
#endif
  return _momentOfInertias[i];
}
 
template <typename floatType, typename intType>
std::vector<std::array<floatType, 3>> ParticlePropertiesLibrary<floatType, intType>::getSitePositions(intType i) const {
#if MD_FLEXIBLE_MODE == SINGLESITE
  AutoPasLog(WARN,
             "ParticlePropertiesLibrary::getSitePositions(): trying to get the site positions of a multi-site molecule "
             "type when md-flexible has been compiled without support for multi-site molecules. Please compile with "
             "the CMake argument '-DMD_FLEXIBLE_MODE=MULTISITE'.");
#endif
  return _relativeSitePositions[i];
}
 
template <typename floatType, typename intType>
std::vector<intType> ParticlePropertiesLibrary<floatType, intType>::getSiteTypes(intType i) const {
#if MD_FLEXIBLE_MODE == SINGLESITE
  AutoPasLog(WARN,
             "ParticlePropertiesLibrary::getSiteTypes(): trying to get the site types of a multi-site molecule type "
             "when md-flexible has been compiled without support for multi-site molecules. Please compile with the "
             "CMake argument '-DMD_FLEXIBLE_MODE=MULTISITE'.");
#endif
  return _siteIds[i];
}
 
template <typename floatType, typename intType>
floatType ParticlePropertiesLibrary<floatType, intType>::getEpsilon(intType i) const {
  return _epsilons[i];
}
 
template <typename floatType, typename intType>
floatType ParticlePropertiesLibrary<floatType, intType>::getSigma(intType i) const {
  return _sigmas[i];
}
 
template <typename floatType, typename intType>
floatType ParticlePropertiesLibrary<floatType, intType>::getNu(intType i) const {
  return _nus[i];
}
 
template <typename floatType, typename intType>
intType ParticlePropertiesLibrary<floatType, intType>::getNumSites(intType i) const {
#if MD_FLEXIBLE_MODE == SINGLESITE
  AutoPasLog(WARN,
             "ParticlePropertiesLibrary::getNumSites(): trying to get the number of sites of a multi-site molecule "
             "type when md-flexible has been compiled without support for multi-site molecules. Please compile with "
             "the CMake argument '-DMD_FLEXIBLE_MODE=MULTISITE'.");
#endif
  return _numSites[i];
}
 
template <typename floatType, typename intType>
floatType ParticlePropertiesLibrary<floatType, intType>::getMoleculesLargestSigma(intType i) const {
#if MD_FLEXIBLE_MODE == SINGLESITE
  AutoPasLog(WARN,
             "ParticlePropertiesLibrary::getNumSites(): trying to get the number of sites of a multi-site molecule "
             "type when md-flexible has been compiled without support for multi-site molecules. Please compile with "
             "the CMake argument '-DMD_FLEXIBLE_MODE=MULTISITE'.");
#endif
  return _moleculesLargestSigma[i];
}
 
template <typename floatType, typename intType>
double ParticlePropertiesLibrary<floatType, intType>::calcShift6(double epsilon24, double sigmaSquared,
                                                                 double cutoffSquared) {
  const auto sigmaDivCutoffPow2 = sigmaSquared / cutoffSquared;
  const auto sigmaDivCutoffPow6 = sigmaDivCutoffPow2 * sigmaDivCutoffPow2 * sigmaDivCutoffPow2;
  const auto shift6 = epsilon24 * (sigmaDivCutoffPow6 - sigmaDivCutoffPow6 * sigmaDivCutoffPow6);
  return shift6;
}