doxygen_documentation/git-master/LJMultisiteFunctor_8h_source.html

#pragma once


#include "MoleculeLJ.h"

#include "MultisiteMoleculeLJ.h"

#include "ParticlePropertiesLibrary.h"

#include "autopas/baseFunctors/PairwiseFunctor.h"

#include "autopas/particles/OwnershipState.h"

#include "autopas/utils/AlignedAllocator.h"

#include "autopas/utils/ArrayMath.h"

#include "autopas/utils/ExceptionHandler.h"

#include "autopas/utils/Math.h"

#include "autopas/utils/Quaternion.h"

#include "autopas/utils/SoA.h"

#include "autopas/utils/StaticBoolSelector.h"

#include "autopas/utils/WrapOpenMP.h"

#include "autopas/utils/inBox.h"


namespace mdLib {


template <class Particle_T, bool applyShift = false, bool useMixing = false,

          autopas::FunctorN3Modes useNewton3 = autopas::FunctorN3Modes::Both, bool calculateGlobals = false,

          bool countFLOPs = false, bool relevantForTuning = true>

class LJMultisiteFunctor

    : public autopas::PairwiseFunctor<Particle_T, LJMultisiteFunctor<Particle_T, applyShift, useMixing, useNewton3,

                                                                     calculateGlobals, relevantForTuning, countFLOPs>> {

  using SoAArraysType = typename Particle_T::SoAArraysType;


  using SoAFloatPrecision = typename Particle_T::ParticleSoAFloatPrecision;


  const double _cutoffSquared;


  double _epsilon24;


  double _sigmaSquared;


  double _shift6 = 0;


  const std::vector<std::array<double, 3>> _sitePositionsLJ{};


  ParticlePropertiesLibrary<SoAFloatPrecision, size_t> *_PPLibrary = nullptr;


  double _potentialEnergySum;


  std::array<double, 3> _virialSum;


  bool _postProcessed;


 public:

  LJMultisiteFunctor() = delete;


 private:

  explicit LJMultisiteFunctor(SoAFloatPrecision cutoff, void * /*dummy*/)

      : autopas::PairwiseFunctor<Particle_T, LJMultisiteFunctor<Particle_T, applyShift, useMixing, useNewton3,

                                                                calculateGlobals, relevantForTuning, countFLOPs>>(

            cutoff),

        _cutoffSquared{cutoff * cutoff},

        _potentialEnergySum{0.},

        _virialSum{0., 0., 0.},

        _aosThreadData(),

        _postProcessed{false} {

    if constexpr (calculateGlobals) {

      _aosThreadData.resize(autopas::autopas_get_max_threads());

    }

    if constexpr (countFLOPs) {

      AutoPasLog(DEBUG, "FLOP counting is enabled but is not supported for multi-site functors yet!");

    }

  }


 public:

  explicit LJMultisiteFunctor(double cutoff) : LJMultisiteFunctor(cutoff, nullptr) {

    static_assert(not useMixing,

                  "Mixing without a ParticlePropertiesLibrary is not possible! Use a different constructor or set "

                  "mixing to false.");

    AutoPasLog(WARN, "Using LJMultisiteFunctor with mixing disabled is untested!");

  }


  explicit LJMultisiteFunctor(double cutoff, ParticlePropertiesLibrary<double, size_t> &particlePropertiesLibrary)

      : LJMultisiteFunctor(cutoff, nullptr) {

    static_assert(useMixing,

                  "Not using Mixing but using a ParticlePropertiesLibrary is not allowed! Use a different constructor "

                  "or set mixing to true.");

    _PPLibrary = &particlePropertiesLibrary;

  }


  std::string getName() final { return "LJMultisiteFunctor"; }


  bool isRelevantForTuning() final { return relevantForTuning; }


  bool allowsNewton3() final {

    return useNewton3 == autopas::FunctorN3Modes::Newton3Only or useNewton3 == autopas::FunctorN3Modes::Both;

  }


  bool allowsNonNewton3() final {

    return useNewton3 == autopas::FunctorN3Modes::Newton3Off or useNewton3 == autopas::FunctorN3Modes::Both;

  }


  void AoSFunctor(Particle_T &particleA, Particle_T &particleB, bool newton3) final {

    using namespace autopas::utils::ArrayMath::literals;

    if (particleA.isDummy() or particleB.isDummy()) {

      return;

    }


    // Don't calculate force if particleB outside cutoff of particleA

    const auto displacementCoM = autopas::utils::ArrayMath::sub(particleA.getR(), particleB.getR());

    const auto distanceSquaredCoM = autopas::utils::ArrayMath::dot(displacementCoM, displacementCoM);


    if (distanceSquaredCoM > _cutoffSquared) {

      return;

    }


    // get number of sites

    const size_t numSitesA = useMixing ? _PPLibrary->getNumSites(particleA.getTypeId()) : _sitePositionsLJ.size();

    const size_t numSitesB = useMixing ? _PPLibrary->getNumSites(particleB.getTypeId()) : _sitePositionsLJ.size();


    // get siteIds

    const std::vector<size_t> siteIdsA =

        useMixing ? _PPLibrary->getSiteTypes(particleA.getTypeId()) : std::vector<unsigned long>();

    const std::vector<size_t> siteIdsB =

        useMixing ? _PPLibrary->getSiteTypes(particleB.getTypeId()) : std::vector<unsigned long>();


    // get unrotated relative site positions

    const std::vector<std::array<double, 3>> unrotatedSitePositionsA =

        useMixing ? _PPLibrary->getSitePositions(particleA.getTypeId()) : _sitePositionsLJ;

    const std::vector<std::array<double, 3>> unrotatedSitePositionsB =

        useMixing ? _PPLibrary->getSitePositions(particleB.getTypeId()) : _sitePositionsLJ;


    // calculate correctly rotated relative site positions

    const auto rotatedSitePositionsA =

        autopas::utils::quaternion::rotateVectorOfPositions(particleA.getQuaternion(), unrotatedSitePositionsA);

    const auto rotatedSitePositionsB =

        autopas::utils::quaternion::rotateVectorOfPositions(particleB.getQuaternion(), unrotatedSitePositionsB);


    for (int i = 0; i < numSitesA; i++) {

      for (int j = 0; j < numSitesB; j++) {

        const auto displacement = autopas::utils::ArrayMath::add(

            autopas::utils::ArrayMath::sub(displacementCoM, rotatedSitePositionsB[j]), rotatedSitePositionsA[i]);

        const auto distanceSquared = autopas::utils::ArrayMath::dot(displacement, displacement);


        const auto sigmaSquared =

            useMixing ? _PPLibrary->getMixingSigmaSquared(siteIdsA[i], siteIdsB[j]) : _sigmaSquared;

        const auto epsilon24 = useMixing ? _PPLibrary->getMixing24Epsilon(siteIdsA[i], siteIdsB[j]) : _epsilon24;

        const auto shift6 =

            applyShift ? (useMixing ? _PPLibrary->getMixingShift6(siteIdsA[i], siteIdsB[j]) : _shift6) : 0;


        // clang-format off

        // Calculate potential between sites and thus force

        // Force = 24 * epsilon * (2*(sigma/distance)^12 - (sigma/distance)^6) * (1/distance)^2 * [x_displacement, y_displacement, z_displacement]

        //         {                         scalarMultiple                                   } * {                     displacement             }

        // clang-format on

        const auto invDistSquared = 1. / distanceSquared;

        const auto lj2 = sigmaSquared * invDistSquared;

        const auto lj6 = lj2 * lj2 * lj2;

        const auto lj12 = lj6 * lj6;

        const auto lj12m6 = lj12 - lj6;  // = LJ potential / (4x epsilon)

        const auto scalarMultiple = epsilon24 * (lj12 + lj12m6) * invDistSquared;

        const auto force = autopas::utils::ArrayMath::mulScalar(displacement, scalarMultiple);


        // Add force on site to net force

        particleA.addF(force);

        if (newton3) {

          particleB.subF(force);

        }


        // Add torque applied by force

        particleA.addTorque(autopas::utils::ArrayMath::cross(rotatedSitePositionsA[i], force));

        if (newton3) {

          particleB.subTorque(autopas::utils::ArrayMath::cross(rotatedSitePositionsB[j], force));

        }


        if (calculateGlobals) {

          // We always add the full contribution for each owned particle and divide the sums by 2 in endTraversal().

          // Potential energy has an additional factor of 6, which is also handled in endTraversal().

          const auto potentialEnergy6 = epsilon24 * lj12m6 + shift6;

          const auto virial = displacement * force;


          const auto threadNum = autopas::autopas_get_thread_num();


          if (particleA.isOwned()) {

            _aosThreadData[threadNum].potentialEnergySum += potentialEnergy6;

            _aosThreadData[threadNum].virialSum += virial;

          }

          // for non-newton3 the second particle will be considered in a separate calculation

          if (newton3 and particleB.isOwned()) {

            _aosThreadData[threadNum].potentialEnergySum += potentialEnergy6;

            _aosThreadData[threadNum].virialSum += virial;

          }

        }

      }

    }

  }


  void SoAFunctorSingle(autopas::SoAView<SoAArraysType> soa, bool newton3) final {

    if (soa.size() == 0) return;


    const auto *const __restrict xptr = soa.template begin<Particle_T::AttributeNames::posX>();

    const auto *const __restrict yptr = soa.template begin<Particle_T::AttributeNames::posY>();

    const auto *const __restrict zptr = soa.template begin<Particle_T::AttributeNames::posZ>();


    const auto *const __restrict ownedStatePtr = soa.template begin<Particle_T::AttributeNames::ownershipState>();


    const auto *const __restrict q0ptr = soa.template begin<Particle_T::AttributeNames::quaternion0>();

    const auto *const __restrict q1ptr = soa.template begin<Particle_T::AttributeNames::quaternion1>();

    const auto *const __restrict q2ptr = soa.template begin<Particle_T::AttributeNames::quaternion2>();

    const auto *const __restrict q3ptr = soa.template begin<Particle_T::AttributeNames::quaternion3>();


    SoAFloatPrecision *const __restrict fxptr = soa.template begin<Particle_T::AttributeNames::forceX>();

    SoAFloatPrecision *const __restrict fyptr = soa.template begin<Particle_T::AttributeNames::forceY>();

    SoAFloatPrecision *const __restrict fzptr = soa.template begin<Particle_T::AttributeNames::forceZ>();


    SoAFloatPrecision *const __restrict txptr = soa.template begin<Particle_T::AttributeNames::torqueX>();

    SoAFloatPrecision *const __restrict typtr = soa.template begin<Particle_T::AttributeNames::torqueY>();

    SoAFloatPrecision *const __restrict tzptr = soa.template begin<Particle_T::AttributeNames::torqueZ>();


    [[maybe_unused]] auto *const __restrict typeptr = soa.template begin<Particle_T::AttributeNames::typeId>();


    SoAFloatPrecision potentialEnergySum = 0.;

    SoAFloatPrecision virialSumX = 0.;

    SoAFloatPrecision virialSumY = 0.;

    SoAFloatPrecision virialSumZ = 0.;


    // the local redeclaration of the following values helps the SoAFloatPrecision-generation of various compilers.

    const SoAFloatPrecision cutoffSquared = _cutoffSquared;


    std::vector<SoAFloatPrecision, autopas::AlignedAllocator<SoAFloatPrecision>> sigmaSquareds;

    std::vector<SoAFloatPrecision, autopas::AlignedAllocator<SoAFloatPrecision>> epsilon24s;

    std::vector<SoAFloatPrecision, autopas::AlignedAllocator<SoAFloatPrecision>> shift6s;


    std::vector<SoAFloatPrecision, autopas::AlignedAllocator<SoAFloatPrecision>> exactSitePositionX;

    std::vector<SoAFloatPrecision, autopas::AlignedAllocator<SoAFloatPrecision>> exactSitePositionY;

    std::vector<SoAFloatPrecision, autopas::AlignedAllocator<SoAFloatPrecision>> exactSitePositionZ;


    // we require arrays for forces for sites to maintain SIMD in site-site calculations

    std::vector<SoAFloatPrecision, autopas::AlignedAllocator<SoAFloatPrecision>> siteForceX;

    std::vector<SoAFloatPrecision, autopas::AlignedAllocator<SoAFloatPrecision>> siteForceY;

    std::vector<SoAFloatPrecision, autopas::AlignedAllocator<SoAFloatPrecision>> siteForceZ;


    std::vector<size_t, autopas::AlignedAllocator<size_t>> siteTypes;

    std::vector<char, autopas::AlignedAllocator<char>> isSiteOwned;


    const SoAFloatPrecision const_sigmaSquared = _sigmaSquared;

    const SoAFloatPrecision const_epsilon24 = _epsilon24;

    const SoAFloatPrecision const_shift6 = _shift6;


    const auto const_unrotatedSitePositions = _sitePositionsLJ;


    // count number of sites in SoA

    size_t siteCount = 0;

    if constexpr (useMixing) {

      for (size_t mol = 0; mol < soa.size(); ++mol) {

        siteCount += _PPLibrary->getNumSites(typeptr[mol]);

      }

    } else {

      siteCount = const_unrotatedSitePositions.size() * soa.size();

    }


    // pre-reserve site std::vectors

    exactSitePositionX.reserve(siteCount);

    exactSitePositionY.reserve(siteCount);

    exactSitePositionZ.reserve(siteCount);


    if constexpr (useMixing) {

      siteTypes.reserve(siteCount);

    }


    siteForceX.reserve((siteCount));

    siteForceY.reserve((siteCount));

    siteForceZ.reserve((siteCount));


    if constexpr (calculateGlobals) {

      // this is only needed for vectorization when calculating globals

      isSiteOwned.reserve(siteCount);

    }


    // Fill site-wise std::vectors for SIMD

    if constexpr (useMixing) {

      size_t siteIndex = 0;

      for (size_t mol = 0; mol < soa.size(); ++mol) {

        const auto rotatedSitePositions = autopas::utils::quaternion::rotateVectorOfPositions(

            {q0ptr[mol], q1ptr[mol], q2ptr[mol], q3ptr[mol]}, _PPLibrary->getSitePositions(typeptr[mol]));

        const auto siteTypesOfMol = _PPLibrary->getSiteTypes(typeptr[mol]);


        for (size_t site = 0; site < _PPLibrary->getNumSites(typeptr[mol]); ++site) {

          exactSitePositionX[siteIndex] = rotatedSitePositions[site][0] + xptr[mol];

          exactSitePositionY[siteIndex] = rotatedSitePositions[site][1] + yptr[mol];

          exactSitePositionZ[siteIndex] = rotatedSitePositions[site][2] + zptr[mol];

          siteTypes[siteIndex] = siteTypesOfMol[site];

          siteForceX[siteIndex] = 0.;

          siteForceY[siteIndex] = 0.;

          siteForceZ[siteIndex] = 0.;

          if (calculateGlobals) {

            isSiteOwned[siteIndex] = ownedStatePtr[mol] == autopas::OwnershipState::owned;

          }

          ++siteIndex;

        }

      }

    } else {

      size_t siteIndex = 0;

      for (size_t mol = 0; mol < soa.size(); mol++) {

        const auto rotatedSitePositions = autopas::utils::quaternion::rotateVectorOfPositions(

            {q0ptr[mol], q1ptr[mol], q2ptr[mol], q3ptr[mol]}, const_unrotatedSitePositions);

        for (size_t site = 0; site < const_unrotatedSitePositions.size(); ++site) {

          exactSitePositionX[siteIndex] = rotatedSitePositions[site][0] + xptr[mol];

          exactSitePositionY[siteIndex] = rotatedSitePositions[site][1] + yptr[mol];

          exactSitePositionZ[siteIndex] = rotatedSitePositions[site][2] + zptr[mol];

          siteForceX[siteIndex] = 0.;

          siteForceY[siteIndex] = 0.;

          siteForceZ[siteIndex] = 0.;

          if (calculateGlobals) {

            isSiteOwned[siteIndex] = ownedStatePtr[mol] == autopas::OwnershipState::owned;

          }

          ++siteIndex;

        }

      }

    }


    // main force calculation loop

    size_t siteIndexMolA = 0;  // index of first site in molA

    for (size_t molA = 0; molA < soa.size(); ++molA) {

      const size_t noSitesInMolA = useMixing ? _PPLibrary->getNumSites(typeptr[molA])

                                             : const_unrotatedSitePositions.size();  // Number of sites in molecule A


      const auto ownedStateA = ownedStatePtr[molA];

      if (ownedStateA == autopas::OwnershipState::dummy) {

        siteIndexMolA += noSitesInMolA;

        continue;

      }


      const size_t siteIndexMolB = siteIndexMolA + noSitesInMolA;  // index of first site in molB

      const size_t noSitesB = (siteCount - siteIndexMolB);         // Number of sites in molecules that A interacts with


      // create mask over every mol 'above' molA  (char to keep arrays aligned)

      std::vector<char, autopas::AlignedAllocator<char>> molMask;

      molMask.reserve(soa.size() - (molA + 1));


#pragma omp simd

      for (size_t molB = molA + 1; molB < soa.size(); ++molB) {

        const auto ownedStateB = ownedStatePtr[molB];


        const auto displacementCoMX = xptr[molA] - xptr[molB];

        const auto displacementCoMY = yptr[molA] - yptr[molB];

        const auto displacementCoMZ = zptr[molA] - zptr[molB];


        const auto distanceSquaredCoMX = displacementCoMX * displacementCoMX;

        const auto distanceSquaredCoMY = displacementCoMY * displacementCoMY;

        const auto distanceSquaredCoMZ = displacementCoMZ * displacementCoMZ;


        const auto distanceSquaredCoM = distanceSquaredCoMX + distanceSquaredCoMY + distanceSquaredCoMZ;


        // mask sites of molecules beyond cutoff or if molecule is a dummy

        molMask[molB - (molA + 1)] =

            distanceSquaredCoM <= cutoffSquared and ownedStateB != autopas::OwnershipState::dummy;

      }


      // generate mask for each site in the mols 'above' molA from molecular mask

      std::vector<char, autopas::AlignedAllocator<char>> siteMask;

      siteMask.reserve(noSitesB);


      for (size_t molB = molA + 1; molB < soa.size(); ++molB) {

        for (size_t siteB = 0; siteB < _PPLibrary->getNumSites(typeptr[molB]); ++siteB) {

          siteMask.emplace_back(molMask[molB - (molA + 1)]);

        }

      }


      // calculate LJ forces

      for (size_t siteA = siteIndexMolA; siteA < siteIndexMolB; ++siteA) {

        if (useMixing) {

          // preload sigmas, epsilons, and shifts

          sigmaSquareds.reserve(noSitesB);

          epsilon24s.reserve(noSitesB);

          if constexpr (applyShift) {

            shift6s.reserve(noSitesB);

          }


          for (size_t siteB = 0; siteB < siteCount - (siteIndexMolB); ++siteB) {

            const auto mixingData = _PPLibrary->getLJMixingData(siteTypes[siteA], siteTypes[siteIndexMolB + siteB]);

            sigmaSquareds[siteB] = mixingData.sigmaSquared;

            epsilon24s[siteB] = mixingData.epsilon24;

            if (applyShift) {

              shift6s[siteB] = mixingData.shift6;

            }

          }

        }

        // sums used for siteA

        SoAFloatPrecision forceSumX = 0.;

        SoAFloatPrecision forceSumY = 0.;

        SoAFloatPrecision forceSumZ = 0.;

        SoAFloatPrecision torqueSumX = 0.;

        SoAFloatPrecision torqueSumY = 0.;

        SoAFloatPrecision torqueSumZ = 0.;


#pragma omp simd reduction (+ : forceSumX, forceSumY, forceSumZ, torqueSumX, torqueSumY, torqueSumZ, potentialEnergySum, virialSumX, virialSumY, virialSumZ)

        for (size_t siteB = 0; siteB < noSitesB; ++siteB) {

          const size_t globalSiteBIndex = siteB + siteIndexMolB;


          const SoAFloatPrecision sigmaSquared = useMixing ? sigmaSquareds[siteB] : const_sigmaSquared;

          const SoAFloatPrecision epsilon24 = useMixing ? epsilon24s[siteB] : const_epsilon24;

          const SoAFloatPrecision shift6 = applyShift ? (useMixing ? shift6s[siteB] : const_shift6) : 0;


          const auto isSiteBOwned = !calculateGlobals || isSiteOwned[globalSiteBIndex];


          const auto displacementX = exactSitePositionX[siteA] - exactSitePositionX[globalSiteBIndex];

          const auto displacementY = exactSitePositionY[siteA] - exactSitePositionY[globalSiteBIndex];

          const auto displacementZ = exactSitePositionZ[siteA] - exactSitePositionZ[globalSiteBIndex];


          const auto distanceSquaredX = displacementX * displacementX;

          const auto distanceSquaredY = displacementY * displacementY;

          const auto distanceSquaredZ = displacementZ * displacementZ;


          const auto distanceSquared = distanceSquaredX + distanceSquaredY + distanceSquaredZ;


          const auto invDistSquared = 1. / distanceSquared;

          const auto lj2 = sigmaSquared * invDistSquared;

          const auto lj6 = lj2 * lj2 * lj2;

          const auto lj12 = lj6 * lj6;

          const auto lj12m6 = lj12 - lj6;

          const auto scalarMultiple = siteMask[siteB] ? epsilon24 * (lj12 + lj12m6) * invDistSquared : 0.;


          // calculate forces

          const auto forceX = scalarMultiple * displacementX;

          const auto forceY = scalarMultiple * displacementY;

          const auto forceZ = scalarMultiple * displacementZ;


          forceSumX += forceX;

          forceSumY += forceY;

          forceSumZ += forceZ;


          // newton's third law

          siteForceX[globalSiteBIndex] -= forceX;

          siteForceY[globalSiteBIndex] -= forceY;

          siteForceZ[globalSiteBIndex] -= forceZ;


          if constexpr (calculateGlobals) {

            const auto virialX = displacementX * forceX;

            const auto virialY = displacementY * forceY;

            const auto virialZ = displacementZ * forceZ;

            const auto potentialEnergy6 = siteMask[siteB] ? (epsilon24 * lj12m6 + shift6) : 0.;


            // We add 6 times the potential energy for each owned particle. The total sum is corrected in

            // endTraversal().

            const auto ownershipMask =

                (ownedStateA == autopas::OwnershipState::owned ? 1. : 0.) + (isSiteBOwned ? 1. : 0.);

            potentialEnergySum += potentialEnergy6 * ownershipMask;

            virialSumX += virialX * ownershipMask;

            virialSumY += virialY * ownershipMask;

            virialSumZ += virialZ * ownershipMask;

          }

        }

        // sum forces on single site in mol A

        siteForceX[siteA] += forceSumX;

        siteForceY[siteA] += forceSumY;

        siteForceZ[siteA] += forceSumZ;

      }

      siteIndexMolA += noSitesInMolA;

    }


    // reduce the forces on individual sites to forces & torques on whole molecules.

    if constexpr (useMixing) {

      size_t siteIndex = 0;

      for (size_t mol = 0; mol < soa.size(); ++mol) {

        const auto rotatedSitePositions = autopas::utils::quaternion::rotateVectorOfPositions(

            {q0ptr[mol], q1ptr[mol], q2ptr[mol], q3ptr[mol]}, _PPLibrary->getSitePositions(typeptr[mol]));

        for (size_t site = 0; site < _PPLibrary->getNumSites(typeptr[mol]); ++site) {

          fxptr[mol] += siteForceX[siteIndex];

          fyptr[mol] += siteForceY[siteIndex];

          fzptr[mol] += siteForceZ[siteIndex];

          txptr[mol] += rotatedSitePositions[site][1] * siteForceZ[siteIndex] -

                        rotatedSitePositions[site][2] * siteForceY[siteIndex];

          typtr[mol] += rotatedSitePositions[site][2] * siteForceX[siteIndex] -

                        rotatedSitePositions[site][0] * siteForceZ[siteIndex];

          tzptr[mol] += rotatedSitePositions[site][0] * siteForceY[siteIndex] -

                        rotatedSitePositions[site][1] * siteForceX[siteIndex];

          ++siteIndex;

        }

      }

    } else {

      size_t siteIndex = 0;

      for (size_t mol = 0; mol < soa.size(); mol++) {

        const auto rotatedSitePositions = autopas::utils::quaternion::rotateVectorOfPositions(

            {q0ptr[mol], q1ptr[mol], q2ptr[mol], q3ptr[mol]}, const_unrotatedSitePositions);

        for (size_t site = 0; site < const_unrotatedSitePositions.size(); ++site) {

          fxptr[mol] += siteForceX[siteIndex];

          fyptr[mol] += siteForceY[siteIndex];

          fzptr[mol] += siteForceZ[siteIndex];

          txptr[mol] += rotatedSitePositions[site][1] * siteForceZ[siteIndex] -

                        rotatedSitePositions[site][2] * siteForceY[siteIndex];

          typtr[mol] += rotatedSitePositions[site][2] * siteForceX[siteIndex] -

                        rotatedSitePositions[site][0] * siteForceZ[siteIndex];

          tzptr[mol] += rotatedSitePositions[site][0] * siteForceY[siteIndex] -

                        rotatedSitePositions[site][1] * siteForceX[siteIndex];

          ++siteIndex;

        }

      }

    }


    if constexpr (calculateGlobals) {

      const auto threadNum = autopas::autopas_get_thread_num();


      _aosThreadData[threadNum].potentialEnergySum += potentialEnergySum;

      _aosThreadData[threadNum].virialSum[0] += virialSumX;

      _aosThreadData[threadNum].virialSum[1] += virialSumY;

      _aosThreadData[threadNum].virialSum[2] += virialSumZ;

    }

  }

  void SoAFunctorPair(autopas::SoAView<SoAArraysType> soa1, autopas::SoAView<SoAArraysType> soa2,

                      const bool newton3) final {

    if (newton3) {

      SoAFunctorPairImpl<true>(soa1, soa2);

    } else {

      SoAFunctorPairImpl<false>(soa1, soa2);

    }

  }


  // clang-format off

  // clang-format on

  void SoAFunctorVerlet(autopas::SoAView<SoAArraysType> soa, const size_t indexFirst,

                        const std::vector<size_t, autopas::AlignedAllocator<size_t>> &neighborList,

                        bool newton3) final {

    if (soa.size() == 0 or neighborList.empty()) return;

    if (newton3) {

      SoAFunctorVerletImpl<true>(soa, indexFirst, neighborList);

    } else {

      SoAFunctorVerletImpl<false>(soa, indexFirst, neighborList);

    }

  }


  void setParticleProperties(SoAFloatPrecision epsilon24, SoAFloatPrecision sigmaSquared,

                             std::vector<std::array<SoAFloatPrecision, 3>> sitePositionsLJ) {

    _epsilon24 = epsilon24;

    _sigmaSquared = sigmaSquared;

    if (applyShift) {

      _shift6 = ParticlePropertiesLibrary<double, size_t>::calcShift6(_epsilon24, _sigmaSquared, _cutoffSquared);

    } else {

      _shift6 = 0;

    }

    _sitePositionsLJ = sitePositionsLJ;

  }


  constexpr static auto getNeededAttr() {

    return std::array<typename Particle_T::AttributeNames, 16>{

        Particle_T::AttributeNames::id,          Particle_T::AttributeNames::posX,

        Particle_T::AttributeNames::posY,        Particle_T::AttributeNames::posZ,

        Particle_T::AttributeNames::forceX,      Particle_T::AttributeNames::forceY,

        Particle_T::AttributeNames::forceZ,      Particle_T::AttributeNames::quaternion0,

        Particle_T::AttributeNames::quaternion1, Particle_T::AttributeNames::quaternion2,

        Particle_T::AttributeNames::quaternion3, Particle_T::AttributeNames::torqueX,

        Particle_T::AttributeNames::torqueY,     Particle_T::AttributeNames::torqueZ,

        Particle_T::AttributeNames::typeId,      Particle_T::AttributeNames::ownershipState};

  }


  constexpr static auto getNeededAttr(std::false_type) {

    return std::array<typename Particle_T::AttributeNames, 16>{

        Particle_T::AttributeNames::id,          Particle_T::AttributeNames::posX,

        Particle_T::AttributeNames::posY,        Particle_T::AttributeNames::posZ,

        Particle_T::AttributeNames::forceX,      Particle_T::AttributeNames::forceY,

        Particle_T::AttributeNames::forceZ,      Particle_T::AttributeNames::quaternion0,

        Particle_T::AttributeNames::quaternion1, Particle_T::AttributeNames::quaternion2,

        Particle_T::AttributeNames::quaternion3, Particle_T::AttributeNames::torqueX,

        Particle_T::AttributeNames::torqueY,     Particle_T::AttributeNames::torqueZ,

        Particle_T::AttributeNames::typeId,      Particle_T::AttributeNames::ownershipState};

  }


  constexpr static auto getComputedAttr() {

    return std::array<typename Particle_T::AttributeNames, 6>{

        Particle_T::AttributeNames::forceX,  Particle_T::AttributeNames::forceY,  Particle_T::AttributeNames::forceZ,

        Particle_T::AttributeNames::torqueX, Particle_T::AttributeNames::torqueY, Particle_T::AttributeNames::torqueZ};

  }


  constexpr static bool getMixing() { return useMixing; }


  unsigned long getNumFlopsPerKernelCall(size_t molAType, size_t molBType, bool newton3) {

    // Site-to-site displacement: 6 (3 in the SoA case, but this requires O(N) precomputing site positions)

    // Site-to-site distance squared: 4

    // Compute scale: 9

    // Apply scale to force: With newton3: 6, Without: 3

    // Apply scale to torque: With newton3 18, Without: 9 (0 in SoA case, with O(N) post computing)

    // Site-to-site total: With newton3: 33, Without: 26

    // (SoA total: With N3L: 22, Without N3L: 19)

    // Above multiplied by number sites of i * number sites of j

    const unsigned long siteToSiteFlops = newton3 ? 33ul : 26ul;

    return _PPLibrary->getNumSites(molAType) * _PPLibrary->getNumSites(molBType) * siteToSiteFlops;

  }


  void initTraversal() final {

    _potentialEnergySum = 0;

    _virialSum = {0., 0., 0.};

    _postProcessed = false;

    for (size_t i = 0; i < _aosThreadData.size(); i++) {

      _aosThreadData[i].setZero();

    }

  }


  void endTraversal(bool newton3) final {

    using namespace autopas::utils::ArrayMath::literals;


    if (_postProcessed) {

      throw autopas::utils::ExceptionHandler::AutoPasException(

          "Already postprocessed, endTraversal(bool newton3) was called twice without calling initTraversal().");

    }

    if (calculateGlobals) {

      for (size_t i = 0; i < _aosThreadData.size(); ++i) {

        _potentialEnergySum += _aosThreadData[i].potentialEnergySum;

        _virialSum += _aosThreadData[i].virialSum;

      }


      // For each interaction, we added the full contribution for both particles. Divide by 2 here, so that each

      // contribution is only counted once per pair.

      _potentialEnergySum *= 0.5;

      _virialSum *= 0.5;


      // We have always calculated 6*potentialEnergy, so we divide by 6 here!

      _potentialEnergySum /= 6.;

      _postProcessed = true;


      AutoPasLog(DEBUG, "Final potential energy {}", _potentialEnergySum);

      AutoPasLog(DEBUG, "Final virial           {}", _virialSum[0] + _virialSum[1] + _virialSum[2]);

    }

  }


  double getPotentialEnergy() {

    if (not calculateGlobals) {

      throw autopas::utils::ExceptionHandler::AutoPasException(

          "Trying to get potential energy even though calculateGlobals is false. If you want this functor to calculate "

          "global "

          "values, please specify calculateGlobals to be true.");

    }

    if (not _postProcessed) {

      throw autopas::utils::ExceptionHandler::AutoPasException(

          "Cannot get potential energy, because endTraversal was not called.");

    }

    return _potentialEnergySum;

  }


  double getVirial() {

    if (not calculateGlobals) {

      throw autopas::utils::ExceptionHandler::AutoPasException(

          "Trying to get virial even though calculateGlobals is false. If you want this functor to calculate global "

          "values, please specify calculateGlobals to be true.");

    }

    if (not _postProcessed) {

      throw autopas::utils::ExceptionHandler::AutoPasException(

          "Cannot get virial, because endTraversal was not called.");

    }

    return _virialSum[0] + _virialSum[1] + _virialSum[2];

  }


 private:

  template <bool newton3>

  void SoAFunctorPairImpl(autopas::SoAView<SoAArraysType> soaA, autopas::SoAView<SoAArraysType> soaB) {

    if (soaA.size() == 0 || soaB.size() == 0) return;


    const auto *const __restrict xAptr = soaA.template begin<Particle_T::AttributeNames::posX>();

    const auto *const __restrict yAptr = soaA.template begin<Particle_T::AttributeNames::posY>();

    const auto *const __restrict zAptr = soaA.template begin<Particle_T::AttributeNames::posZ>();

    const auto *const __restrict xBptr = soaB.template begin<Particle_T::AttributeNames::posX>();

    const auto *const __restrict yBptr = soaB.template begin<Particle_T::AttributeNames::posY>();

    const auto *const __restrict zBptr = soaB.template begin<Particle_T::AttributeNames::posZ>();


    const auto *const __restrict ownedStatePtrA = soaA.template begin<Particle_T::AttributeNames::ownershipState>();

    const auto *const __restrict ownedStatePtrB = soaB.template begin<Particle_T::AttributeNames::ownershipState>();


    const auto *const __restrict q0Aptr = soaA.template begin<Particle_T::AttributeNames::quaternion0>();

    const auto *const __restrict q1Aptr = soaA.template begin<Particle_T::AttributeNames::quaternion1>();

    const auto *const __restrict q2Aptr = soaA.template begin<Particle_T::AttributeNames::quaternion2>();

    const auto *const __restrict q3Aptr = soaA.template begin<Particle_T::AttributeNames::quaternion3>();

    const auto *const __restrict q0Bptr = soaB.template begin<Particle_T::AttributeNames::quaternion0>();

    const auto *const __restrict q1Bptr = soaB.template begin<Particle_T::AttributeNames::quaternion1>();

    const auto *const __restrict q2Bptr = soaB.template begin<Particle_T::AttributeNames::quaternion2>();

    const auto *const __restrict q3Bptr = soaB.template begin<Particle_T::AttributeNames::quaternion3>();


    SoAFloatPrecision *const __restrict fxAptr = soaA.template begin<Particle_T::AttributeNames::forceX>();

    SoAFloatPrecision *const __restrict fyAptr = soaA.template begin<Particle_T::AttributeNames::forceY>();

    SoAFloatPrecision *const __restrict fzAptr = soaA.template begin<Particle_T::AttributeNames::forceZ>();

    SoAFloatPrecision *const __restrict fxBptr = soaB.template begin<Particle_T::AttributeNames::forceX>();

    SoAFloatPrecision *const __restrict fyBptr = soaB.template begin<Particle_T::AttributeNames::forceY>();

    SoAFloatPrecision *const __restrict fzBptr = soaB.template begin<Particle_T::AttributeNames::forceZ>();


    SoAFloatPrecision *const __restrict txAptr = soaA.template begin<Particle_T::AttributeNames::torqueX>();

    SoAFloatPrecision *const __restrict tyAptr = soaA.template begin<Particle_T::AttributeNames::torqueY>();

    SoAFloatPrecision *const __restrict tzAptr = soaA.template begin<Particle_T::AttributeNames::torqueZ>();

    SoAFloatPrecision *const __restrict txBptr = soaB.template begin<Particle_T::AttributeNames::torqueX>();

    SoAFloatPrecision *const __restrict tyBptr = soaB.template begin<Particle_T::AttributeNames::torqueY>();

    SoAFloatPrecision *const __restrict tzBptr = soaB.template begin<Particle_T::AttributeNames::torqueZ>();


    [[maybe_unused]] auto *const __restrict typeptrA = soaA.template begin<Particle_T::AttributeNames::typeId>();

    [[maybe_unused]] auto *const __restrict typeptrB = soaB.template begin<Particle_T::AttributeNames::typeId>();


    SoAFloatPrecision potentialEnergySum = 0.;

    SoAFloatPrecision virialSumX = 0.;

    SoAFloatPrecision virialSumY = 0.;

    SoAFloatPrecision virialSumZ = 0.;


    // local redeclarations to help compilers

    const SoAFloatPrecision cutoffSquared = _cutoffSquared;


    std::vector<SoAFloatPrecision, autopas::AlignedAllocator<SoAFloatPrecision>> sigmaSquareds;

    std::vector<SoAFloatPrecision, autopas::AlignedAllocator<SoAFloatPrecision>> epsilon24s;

    std::vector<SoAFloatPrecision, autopas::AlignedAllocator<SoAFloatPrecision>> shift6s;


    std::vector<SoAFloatPrecision, autopas::AlignedAllocator<SoAFloatPrecision>> exactSitePositionBx;

    std::vector<SoAFloatPrecision, autopas::AlignedAllocator<SoAFloatPrecision>> exactSitePositionBy;

    std::vector<SoAFloatPrecision, autopas::AlignedAllocator<SoAFloatPrecision>> exactSitePositionBz;


    // we require arrays for forces for sites to maintain SIMD in site-site calculations

    std::vector<SoAFloatPrecision, autopas::AlignedAllocator<SoAFloatPrecision>> siteForceBx;

    std::vector<SoAFloatPrecision, autopas::AlignedAllocator<SoAFloatPrecision>> siteForceBy;

    std::vector<SoAFloatPrecision, autopas::AlignedAllocator<SoAFloatPrecision>> siteForceBz;


    std::vector<size_t, autopas::AlignedAllocator<size_t>> siteTypesB;

    std::vector<char, autopas::AlignedAllocator<char>> isSiteOwnedBArr;


    const SoAFloatPrecision const_sigmaSquared = _sigmaSquared;

    const SoAFloatPrecision const_epsilon24 = _epsilon24;

    const SoAFloatPrecision const_shift6 = _shift6;


    const auto const_unrotatedSitePositions = _sitePositionsLJ;


    // count number of sites in both SoAs

    size_t siteCountB = 0;

    if constexpr (useMixing) {

      for (size_t mol = 0; mol < soaB.size(); ++mol) {

        siteCountB += _PPLibrary->getNumSites(typeptrB[mol]);

      }

    } else {

      siteCountB = const_unrotatedSitePositions.size() * soaB.size();

    }


    // pre-reserve std::vectors

    exactSitePositionBx.reserve(siteCountB);

    exactSitePositionBy.reserve(siteCountB);

    exactSitePositionBz.reserve(siteCountB);


    if constexpr (useMixing) {

      siteTypesB.reserve(siteCountB);

    }


    siteForceBx.reserve(siteCountB);

    siteForceBy.reserve(siteCountB);

    siteForceBz.reserve(siteCountB);


    if constexpr (calculateGlobals) {

      // this is only needed for vectorization when calculating globals

      isSiteOwnedBArr.reserve(siteCountB);

    }


    if constexpr (useMixing) {

      siteTypesB.reserve(siteCountB);

      sigmaSquareds.reserve(siteCountB);

      epsilon24s.reserve(siteCountB);

      if constexpr (applyShift) {

        shift6s.reserve(siteCountB);

      }

    }


    // Fill site-wise std::vectors for SIMD

    if constexpr (useMixing) {

      size_t siteIndex = 0;

      for (size_t mol = 0; mol < soaB.size(); ++mol) {

        const auto rotatedSitePositions = autopas::utils::quaternion::rotateVectorOfPositions(

            {q0Bptr[mol], q1Bptr[mol], q2Bptr[mol], q3Bptr[mol]}, _PPLibrary->getSitePositions(typeptrB[mol]));

        const auto siteTypesOfMol = _PPLibrary->getSiteTypes(typeptrB[mol]);


        for (size_t site = 0; site < _PPLibrary->getNumSites(typeptrB[mol]); ++site) {

          exactSitePositionBx[siteIndex] = rotatedSitePositions[site][0] + xBptr[mol];

          exactSitePositionBy[siteIndex] = rotatedSitePositions[site][1] + yBptr[mol];

          exactSitePositionBz[siteIndex] = rotatedSitePositions[site][2] + zBptr[mol];

          siteTypesB[siteIndex] = siteTypesOfMol[site];

          siteForceBx[siteIndex] = 0.;

          siteForceBy[siteIndex] = 0.;

          siteForceBz[siteIndex] = 0.;

          if (calculateGlobals) {

            isSiteOwnedBArr[siteIndex] = ownedStatePtrB[mol] == autopas::OwnershipState::owned;

          }

          ++siteIndex;

        }

      }

    } else {

      size_t siteIndex = 0;

      for (size_t mol = 0; mol < soaB.size(); mol++) {

        const auto rotatedSitePositions = autopas::utils::quaternion::rotateVectorOfPositions(

            {q0Bptr[mol], q1Bptr[mol], q2Bptr[mol], q3Bptr[mol]}, const_unrotatedSitePositions);

        for (size_t site = 0; site < const_unrotatedSitePositions.size(); ++site) {

          exactSitePositionBx[siteIndex] = rotatedSitePositions[site][0] + xBptr[mol];

          exactSitePositionBy[siteIndex] = rotatedSitePositions[site][1] + yBptr[mol];

          exactSitePositionBz[siteIndex] = rotatedSitePositions[site][2] + zBptr[mol];

          siteForceBx[siteIndex] = 0.;

          siteForceBy[siteIndex] = 0.;

          siteForceBz[siteIndex] = 0.;

          if (calculateGlobals) {

            isSiteOwnedBArr[siteIndex] = ownedStatePtrB[mol] == autopas::OwnershipState::owned;

          }

          ++siteIndex;

        }

      }

    }


    // main force calculation loop

    for (size_t molA = 0; molA < soaA.size(); ++molA) {

      const auto ownedStateA = ownedStatePtrA[molA];

      if (ownedStateA == autopas::OwnershipState::dummy) {

        continue;

      }


      const auto noSitesInMolA =

          useMixing ? _PPLibrary->getNumSites(typeptrA[molA]) : const_unrotatedSitePositions.size();

      const auto unrotatedSitePositionsA =

          useMixing ? _PPLibrary->getSitePositions(typeptrA[molA]) : const_unrotatedSitePositions;


      const auto rotatedSitePositionsA = autopas::utils::quaternion::rotateVectorOfPositions(

          {q0Aptr[molA], q1Aptr[molA], q2Aptr[molA], q3Aptr[molA]}, unrotatedSitePositionsA);


      // create mask over every mol in cell B (char to keep arrays aligned)

      std::vector<char, autopas::AlignedAllocator<char>> molMask;

      molMask.reserve(soaB.size());


#pragma omp simd

      for (size_t molB = 0; molB < soaB.size(); ++molB) {

        const auto ownedStateB = ownedStatePtrB[molB];


        const auto displacementCoMX = xAptr[molA] - xBptr[molB];

        const auto displacementCoMY = yAptr[molA] - yBptr[molB];

        const auto displacementCoMZ = zAptr[molA] - zBptr[molB];


        const auto distanceSquaredCoMX = displacementCoMX * displacementCoMX;

        const auto distanceSquaredCoMY = displacementCoMY * displacementCoMY;

        const auto distanceSquaredCoMZ = displacementCoMZ * displacementCoMZ;


        const auto distanceSquaredCoM = distanceSquaredCoMX + distanceSquaredCoMY + distanceSquaredCoMZ;


        // mask sites of molecules beyond cutoff or if molecule is a dummy

        molMask[molB] = distanceSquaredCoM <= cutoffSquared and ownedStateB != autopas::OwnershipState::dummy;

      }


      // generate mask for each site in cell B from molecular mask

      std::vector<char, autopas::AlignedAllocator<char>> siteMask;

      siteMask.reserve(siteCountB);


      for (size_t molB = 0; molB < soaB.size(); ++molB) {

        for (size_t siteB = 0; siteB < _PPLibrary->getNumSites(typeptrB[molB]); ++siteB) {

          siteMask.emplace_back(molMask[molB]);

        }

      }


      // sums used for molA

      SoAFloatPrecision forceSumX = 0.;

      SoAFloatPrecision forceSumY = 0.;

      SoAFloatPrecision forceSumZ = 0.;

      SoAFloatPrecision torqueSumX = 0.;

      SoAFloatPrecision torqueSumY = 0.;

      SoAFloatPrecision torqueSumZ = 0.;


      for (size_t siteA = 0; siteA < noSitesInMolA; ++siteA) {

        if (useMixing) {

          // preload sigmas, epsilons, and shifts

          for (size_t siteB = 0; siteB < siteCountB; ++siteB) {

            const auto mixingData =

                _PPLibrary->getLJMixingData(_PPLibrary->getSiteTypes(typeptrA[molA])[siteA], siteTypesB[siteB]);

            sigmaSquareds[siteB] = mixingData.sigmaSquared;

            epsilon24s[siteB] = mixingData.epsilon24;

            if (applyShift) {

              shift6s[siteB] = mixingData.shift6;

            }

          }

        }


        const auto rotatedSitePositionAx = rotatedSitePositionsA[siteA][0];

        const auto rotatedSitePositionAy = rotatedSitePositionsA[siteA][1];

        const auto rotatedSitePositionAz = rotatedSitePositionsA[siteA][2];


        const auto exactSitePositionAx = rotatedSitePositionAx + xAptr[molA];

        const auto exactSitePositionAy = rotatedSitePositionAy + yAptr[molA];

        const auto exactSitePositionAz = rotatedSitePositionAz + zAptr[molA];


#pragma omp simd reduction (+ : forceSumX, forceSumY, forceSumZ, torqueSumX, torqueSumY, torqueSumZ, potentialEnergySum, virialSumX, virialSumY, virialSumZ)

        for (size_t siteB = 0; siteB < siteCountB; ++siteB) {

          const SoAFloatPrecision sigmaSquared = useMixing ? sigmaSquareds[siteB] : const_sigmaSquared;

          const SoAFloatPrecision epsilon24 = useMixing ? epsilon24s[siteB] : const_epsilon24;

          const SoAFloatPrecision shift6 = applyShift ? (useMixing ? shift6s[siteB] : const_shift6) : 0;


          const auto isSiteOwnedB = !calculateGlobals || isSiteOwnedBArr[siteB];


          const auto displacementX = exactSitePositionAx - exactSitePositionBx[siteB];

          const auto displacementY = exactSitePositionAy - exactSitePositionBy[siteB];

          const auto displacementZ = exactSitePositionAz - exactSitePositionBz[siteB];


          const auto distanceSquaredX = displacementX * displacementX;

          const auto distanceSquaredY = displacementY * displacementY;

          const auto distanceSquaredZ = displacementZ * displacementZ;


          const auto distanceSquared = distanceSquaredX + distanceSquaredY + distanceSquaredZ;


          const auto invDistSquared = 1. / distanceSquared;

          const auto lj2 = sigmaSquared * invDistSquared;

          const auto lj6 = lj2 * lj2 * lj2;

          const auto lj12 = lj6 * lj6;

          const auto lj12m6 = lj12 - lj6;

          const auto scalarMultiple = siteMask[siteB] ? epsilon24 * (lj12 + lj12m6) * invDistSquared : 0.;


          // calculate forces

          const auto forceX = scalarMultiple * displacementX;

          const auto forceY = scalarMultiple * displacementY;

          const auto forceZ = scalarMultiple * displacementZ;


          forceSumX += forceX;

          forceSumY += forceY;

          forceSumZ += forceZ;


          torqueSumX += rotatedSitePositionAy * forceZ - rotatedSitePositionAz * forceY;

          torqueSumY += rotatedSitePositionAz * forceX - rotatedSitePositionAx * forceZ;

          torqueSumZ += rotatedSitePositionAx * forceY - rotatedSitePositionAy * forceX;


          // N3L ( total molecular forces + torques to be determined later )

          if constexpr (newton3) {

            siteForceBx[siteB] -= forceX;

            siteForceBy[siteB] -= forceY;

            siteForceBz[siteB] -= forceZ;

          }


          // globals

          if constexpr (calculateGlobals) {

            const auto potentialEnergy6 = siteMask[siteB] ? (epsilon24 * lj12m6 + shift6) : 0.;

            const auto virialX = displacementX * forceX;

            const auto virialY = displacementY * forceY;

            const auto virialZ = displacementZ * forceZ;


            // We add 6 times the potential energy for each owned particle. The total sum is corrected in

            // endTraversal().

            const auto ownershipFactor =

                newton3 ? (ownedStateA == autopas::OwnershipState::owned ? 1. : 0.) + (isSiteOwnedB ? 1. : 0.)

                        : (ownedStateA == autopas::OwnershipState::owned ? 1. : 0.);

            potentialEnergySum += potentialEnergy6 * ownershipFactor;

            virialSumX += virialX * ownershipFactor;

            virialSumY += virialY * ownershipFactor;

            virialSumZ += virialZ * ownershipFactor;

          }

        }

      }

      fxAptr[molA] += forceSumX;

      fyAptr[molA] += forceSumY;

      fzAptr[molA] += forceSumZ;

      txAptr[molA] += torqueSumX;

      tyAptr[molA] += torqueSumY;

      tzAptr[molA] += torqueSumZ;

    }


    // reduce the forces on individual sites in SoA B to total forces & torques on whole molecules

    if constexpr (useMixing) {

      size_t siteIndex = 0;

      for (size_t mol = 0; mol < soaB.size(); ++mol) {

        const auto rotatedSitePositions = autopas::utils::quaternion::rotateVectorOfPositions(

            {q0Bptr[mol], q1Bptr[mol], q2Bptr[mol], q3Bptr[mol]}, _PPLibrary->getSitePositions(typeptrB[mol]));

        for (size_t site = 0; site < _PPLibrary->getNumSites(typeptrB[mol]); ++site) {

          fxBptr[mol] += siteForceBx[siteIndex];

          fyBptr[mol] += siteForceBy[siteIndex];

          fzBptr[mol] += siteForceBz[siteIndex];

          txBptr[mol] += rotatedSitePositions[site][1] * siteForceBz[siteIndex] -

                         rotatedSitePositions[site][2] * siteForceBy[siteIndex];

          tyBptr[mol] += rotatedSitePositions[site][2] * siteForceBx[siteIndex] -

                         rotatedSitePositions[site][0] * siteForceBz[siteIndex];

          tzBptr[mol] += rotatedSitePositions[site][0] * siteForceBy[siteIndex] -

                         rotatedSitePositions[site][1] * siteForceBx[siteIndex];

          ++siteIndex;

        }

      }

    } else {

      size_t siteIndex = 0;

      for (size_t mol = 0; mol < soaB.size(); ++mol) {

        const auto rotatedSitePositions = autopas::utils::quaternion::rotateVectorOfPositions(

            {q0Bptr[mol], q1Bptr[mol], q2Bptr[mol], q3Bptr[mol]}, const_unrotatedSitePositions);

        for (size_t site = 0; site < const_unrotatedSitePositions.size(); ++site) {

          fxBptr[mol] += siteForceBx[siteIndex];

          fyBptr[mol] += siteForceBy[siteIndex];

          fzBptr[mol] += siteForceBz[siteIndex];

          txBptr[mol] += rotatedSitePositions[site][1] * siteForceBz[siteIndex] -

                         rotatedSitePositions[site][2] * siteForceBy[siteIndex];

          tyBptr[mol] += rotatedSitePositions[site][2] * siteForceBx[siteIndex] -

                         rotatedSitePositions[site][0] * siteForceBz[siteIndex];

          tzBptr[mol] += rotatedSitePositions[site][0] * siteForceBy[siteIndex] -

                         rotatedSitePositions[site][1] * siteForceBx[siteIndex];

          ++siteIndex;

        }

      }

    }

    if constexpr (calculateGlobals) {

      const auto threadNum = autopas::autopas_get_thread_num();

      // SoAFunctorPairImpl obtains the potential energy * 12. For non-newton3, this sum is divided by 12 in

      // post-processing. For newton3, this sum is only divided by 6 in post-processing, so must be divided by 2 here.

      _aosThreadData[threadNum].potentialEnergySum += potentialEnergySum;

      _aosThreadData[threadNum].virialSum[0] += virialSumX;

      _aosThreadData[threadNum].virialSum[1] += virialSumY;

      _aosThreadData[threadNum].virialSum[2] += virialSumZ;

    }

  }


  template <bool newton3>

  void SoAFunctorVerletImpl(autopas::SoAView<SoAArraysType> soa, const size_t indexPrime,

                            const std::vector<size_t, autopas::AlignedAllocator<size_t>> &neighborList) {

    const auto *const __restrict ownedStatePtr = soa.template begin<Particle_T::AttributeNames::ownershipState>();


    // Skip if primary particle is dummy

    const auto ownedStatePrime = ownedStatePtr[indexPrime];

    if (ownedStatePrime == autopas::OwnershipState::dummy) {

      return;

    }


    const auto *const __restrict xptr = soa.template begin<Particle_T::AttributeNames::posX>();

    const auto *const __restrict yptr = soa.template begin<Particle_T::AttributeNames::posY>();

    const auto *const __restrict zptr = soa.template begin<Particle_T::AttributeNames::posZ>();


    const auto *const __restrict q0ptr = soa.template begin<Particle_T::AttributeNames::quaternion0>();

    const auto *const __restrict q1ptr = soa.template begin<Particle_T::AttributeNames::quaternion1>();

    const auto *const __restrict q2ptr = soa.template begin<Particle_T::AttributeNames::quaternion2>();

    const auto *const __restrict q3ptr = soa.template begin<Particle_T::AttributeNames::quaternion3>();


    SoAFloatPrecision *const __restrict fxptr = soa.template begin<Particle_T::AttributeNames::forceX>();

    SoAFloatPrecision *const __restrict fyptr = soa.template begin<Particle_T::AttributeNames::forceY>();

    SoAFloatPrecision *const __restrict fzptr = soa.template begin<Particle_T::AttributeNames::forceZ>();


    SoAFloatPrecision *const __restrict txptr = soa.template begin<Particle_T::AttributeNames::torqueX>();

    SoAFloatPrecision *const __restrict typtr = soa.template begin<Particle_T::AttributeNames::torqueY>();

    SoAFloatPrecision *const __restrict tzptr = soa.template begin<Particle_T::AttributeNames::torqueZ>();


    [[maybe_unused]] auto *const __restrict typeptr = soa.template begin<Particle_T::AttributeNames::typeId>();


    SoAFloatPrecision potentialEnergySum = 0.;

    SoAFloatPrecision virialSumX = 0.;

    SoAFloatPrecision virialSumY = 0.;

    SoAFloatPrecision virialSumZ = 0.;


    // the local redeclaration of the following values helps the SoAFloatPrecision-generation of various compilers.

    const SoAFloatPrecision cutoffSquared = _cutoffSquared;

    const auto const_unrotatedSitePositions = _sitePositionsLJ;


    std::vector<SoAFloatPrecision, autopas::AlignedAllocator<SoAFloatPrecision>> sigmaSquareds;

    std::vector<SoAFloatPrecision, autopas::AlignedAllocator<SoAFloatPrecision>> epsilon24s;

    std::vector<SoAFloatPrecision, autopas::AlignedAllocator<SoAFloatPrecision>> shift6s;


    const auto const_sigmaSquared = _sigmaSquared;

    const auto const_epsilon24 = _epsilon24;

    const auto const_shift6 = _shift6;


    const size_t neighborListSize = neighborList.size();

    const size_t *const __restrict neighborListPtr = neighborList.data();


    // Count sites

    const size_t siteCountMolPrime =

        useMixing ? _PPLibrary->getNumSites(typeptr[indexPrime]) : const_unrotatedSitePositions.size();


    size_t siteCountNeighbors = 0;  // site count of neighbours of primary molecule

    if constexpr (useMixing) {

      for (size_t neighborMol = 0; neighborMol < neighborListSize; ++neighborMol) {

        siteCountNeighbors += _PPLibrary->getNumSites(typeptr[neighborList[neighborMol]]);

      }

    } else {

      siteCountNeighbors = const_unrotatedSitePositions.size() * neighborListSize;

    }


    // initialize site-wise arrays for neighbors

    std::vector<SoAFloatPrecision, autopas::AlignedAllocator<SoAFloatPrecision>> exactNeighborSitePositionX;

    std::vector<SoAFloatPrecision, autopas::AlignedAllocator<SoAFloatPrecision>> exactNeighborSitePositionY;

    std::vector<SoAFloatPrecision, autopas::AlignedAllocator<SoAFloatPrecision>> exactNeighborSitePositionZ;


    std::vector<size_t, autopas::AlignedAllocator<size_t>> siteTypesNeighbors;

    std::vector<char, autopas::AlignedAllocator<char>> isNeighborSiteOwnedArr;


    // we require arrays for forces for sites to maintain SIMD in site-site calculations

    std::vector<SoAFloatPrecision, autopas::AlignedAllocator<SoAFloatPrecision>> siteForceX;

    std::vector<SoAFloatPrecision, autopas::AlignedAllocator<SoAFloatPrecision>> siteForceY;

    std::vector<SoAFloatPrecision, autopas::AlignedAllocator<SoAFloatPrecision>> siteForceZ;


    // pre-reserve arrays

    exactNeighborSitePositionX.reserve(siteCountNeighbors);

    exactNeighborSitePositionY.reserve(siteCountNeighbors);

    exactNeighborSitePositionZ.reserve(siteCountNeighbors);


    if constexpr (useMixing) {

      siteTypesNeighbors.reserve(siteCountNeighbors);

    }


    siteForceX.reserve(siteCountNeighbors);

    siteForceY.reserve(siteCountNeighbors);

    siteForceZ.reserve(siteCountNeighbors);


    if constexpr (calculateGlobals) {

      isNeighborSiteOwnedArr.reserve(siteCountNeighbors);

    }


    if constexpr (useMixing) {

      sigmaSquareds.reserve(siteCountNeighbors);

      epsilon24s.reserve(siteCountNeighbors);

      if constexpr (applyShift) {

        shift6s.reserve(siteCountNeighbors);

      }

    }


    const auto rotatedSitePositionsPrime =

        useMixing ? autopas::utils::quaternion::rotateVectorOfPositions(

                        {q0ptr[indexPrime], q1ptr[indexPrime], q2ptr[indexPrime], q3ptr[indexPrime]},

                        _PPLibrary->getSitePositions(typeptr[indexPrime]))

                  : autopas::utils::quaternion::rotateVectorOfPositions(

                        {q0ptr[indexPrime], q1ptr[indexPrime], q2ptr[indexPrime], q3ptr[indexPrime]},

                        const_unrotatedSitePositions);


    const auto siteTypesPrime = _PPLibrary->getSiteTypes(typeptr[indexPrime]);  // this is not used if non-mixing


    // generate site-wise arrays for neighbors of primary mol

    if constexpr (useMixing) {

      size_t siteIndex = 0;

      for (size_t neighborMol = 0; neighborMol < neighborListSize; ++neighborMol) {

        const auto neighborMolIndex = neighborList[neighborMol];

        const auto rotatedSitePositions = autopas::utils::quaternion::rotateVectorOfPositions(

            {q0ptr[neighborMolIndex], q1ptr[neighborMolIndex], q2ptr[neighborMolIndex], q3ptr[neighborMolIndex]},

            _PPLibrary->getSitePositions(typeptr[neighborMolIndex]));

        const auto siteTypesOfMol = _PPLibrary->getSiteTypes(typeptr[neighborMolIndex]);


        for (size_t site = 0; site < _PPLibrary->getNumSites(typeptr[neighborMolIndex]); ++site) {

          exactNeighborSitePositionX[siteIndex] = rotatedSitePositions[site][0] + xptr[neighborMolIndex];

          exactNeighborSitePositionY[siteIndex] = rotatedSitePositions[site][1] + yptr[neighborMolIndex];

          exactNeighborSitePositionZ[siteIndex] = rotatedSitePositions[site][2] + zptr[neighborMolIndex];

          siteTypesNeighbors[siteIndex] = siteTypesOfMol[site];

          siteForceX[siteIndex] = 0.;

          siteForceY[siteIndex] = 0.;

          siteForceZ[siteIndex] = 0.;

          if (calculateGlobals) {

            isNeighborSiteOwnedArr[siteIndex] = ownedStatePtr[neighborMolIndex] == autopas::OwnershipState::owned;

          }

          ++siteIndex;

        }

      }

    } else {

      size_t siteIndex = 0;

      for (size_t neighborMol = 0; neighborMol < neighborListSize; ++neighborMol) {

        const auto neighborMolIndex = neighborList[neighborMol];

        const auto rotatedSitePositions = autopas::utils::quaternion::rotateVectorOfPositions(

            {q0ptr[neighborMolIndex], q1ptr[neighborMolIndex], q2ptr[neighborMolIndex], q3ptr[neighborMolIndex]},

            const_unrotatedSitePositions);

        for (size_t site = 0; site < const_unrotatedSitePositions.size(); ++site) {

          exactNeighborSitePositionX[siteIndex] = rotatedSitePositions[site][0] + xptr[neighborMolIndex];

          exactNeighborSitePositionY[siteIndex] = rotatedSitePositions[site][1] + yptr[neighborMolIndex];

          exactNeighborSitePositionZ[siteIndex] = rotatedSitePositions[site][2] + zptr[neighborMolIndex];

          siteForceX[siteIndex] = 0.;

          siteForceY[siteIndex] = 0.;

          siteForceZ[siteIndex] = 0.;

          if (calculateGlobals) {

            isNeighborSiteOwnedArr[siteIndex] = ownedStatePtr[neighborMolIndex] == autopas::OwnershipState::owned;

          }

          ++siteIndex;

        }

      }

    }


    // -- main force calculation --


    // - calculate mol mask -

    std::vector<char, autopas::AlignedAllocator<char>> molMask;

    molMask.reserve(neighborListSize);


#pragma omp simd

    for (size_t neighborMol = 0; neighborMol < neighborListSize; ++neighborMol) {

      const auto neighborMolIndex = neighborList[neighborMol];  // index of neighbor mol in soa


      const auto ownedState = ownedStatePtr[neighborMolIndex];


      const auto displacementCoMX = xptr[indexPrime] - xptr[neighborMolIndex];

      const auto displacementCoMY = yptr[indexPrime] - yptr[neighborMolIndex];

      const auto displacementCoMZ = zptr[indexPrime] - zptr[neighborMolIndex];


      const auto distanceSquaredCoMX = displacementCoMX * displacementCoMX;

      const auto distanceSquaredCoMY = displacementCoMY * displacementCoMY;

      const auto distanceSquaredCoMZ = displacementCoMZ * displacementCoMZ;


      const auto distanceSquaredCoM = distanceSquaredCoMX + distanceSquaredCoMY + distanceSquaredCoMZ;


      // mask molecules beyond cutoff or if molecule is a dummy

      molMask[neighborMol] = distanceSquaredCoM <= cutoffSquared and ownedState != autopas::OwnershipState::dummy;

    }


    // generate mask for each site from molecular mask

    std::vector<char, autopas::AlignedAllocator<char>> siteMask;

    siteMask.reserve(siteCountNeighbors);


    for (size_t neighborMol = 0; neighborMol < neighborListSize; ++neighborMol) {

      const auto neighborMolIndex = neighborList[neighborMol];  // index of neighbor mol in soa

      for (size_t siteB = 0; siteB < _PPLibrary->getNumSites(typeptr[neighborMolIndex]); ++siteB) {

        siteMask.emplace_back(molMask[neighborMol]);

      }

    }


    // sums used for prime mol

    SoAFloatPrecision forceSumX = 0.;

    SoAFloatPrecision forceSumY = 0.;

    SoAFloatPrecision forceSumZ = 0.;

    SoAFloatPrecision torqueSumX = 0.;

    SoAFloatPrecision torqueSumY = 0.;

    SoAFloatPrecision torqueSumZ = 0.;


    // - actual LJ calculation -


    for (size_t primeSite = 0; primeSite < siteCountMolPrime; ++primeSite) {

      const auto rotatedPrimeSitePositionX = rotatedSitePositionsPrime[primeSite][0];

      const auto rotatedPrimeSitePositionY = rotatedSitePositionsPrime[primeSite][1];

      const auto rotatedPrimeSitePositionZ = rotatedSitePositionsPrime[primeSite][2];


      const auto exactPrimeSitePositionX = rotatedPrimeSitePositionX + xptr[indexPrime];

      const auto exactPrimeSitePositionY = rotatedPrimeSitePositionY + yptr[indexPrime];

      const auto exactPrimeSitePositionZ = rotatedPrimeSitePositionZ + zptr[indexPrime];


      // generate parameter data for chosen site

      if constexpr (useMixing) {

        const auto primeSiteType = siteTypesPrime[primeSite];


        size_t siteIndex = 0;

        for (size_t neighborMol = 0; neighborMol < neighborListSize; ++neighborMol) {

          const auto neighborMolIndex = neighborList[neighborMol];

          const auto siteTypesOfNeighborMol = _PPLibrary->getSiteTypes(typeptr[neighborMolIndex]);


          for (size_t site = 0; site < _PPLibrary->getNumSites(typeptr[neighborMolIndex]); ++site) {

            const auto mixingData = _PPLibrary->getLJMixingData(primeSiteType, siteTypesOfNeighborMol[site]);

            sigmaSquareds[siteIndex] = mixingData.sigmaSquared;

            epsilon24s[siteIndex] = mixingData.epsilon24;

            if constexpr (applyShift) {

              shift6s[siteIndex] = mixingData.shift6;

            }

            ++siteIndex;

          }

        }

      }


#pragma omp simd reduction(+ : forceSumX, forceSumY, forceSumZ, torqueSumX, torqueSumY, torqueSumZ)

      for (size_t neighborSite = 0; neighborSite < siteCountNeighbors; ++neighborSite) {

        const SoAFloatPrecision sigmaSquared = useMixing ? sigmaSquareds[neighborSite] : const_sigmaSquared;

        const SoAFloatPrecision epsilon24 = useMixing ? epsilon24s[neighborSite] : const_epsilon24;

        const SoAFloatPrecision shift6 = applyShift ? (useMixing ? shift6s[neighborSite] : const_shift6) : 0;


        const bool isNeighborSiteOwned = !calculateGlobals || isNeighborSiteOwnedArr[neighborSite];


        const auto displacementX = exactPrimeSitePositionX - exactNeighborSitePositionX[neighborSite];

        const auto displacementY = exactPrimeSitePositionY - exactNeighborSitePositionY[neighborSite];

        const auto displacementZ = exactPrimeSitePositionZ - exactNeighborSitePositionZ[neighborSite];


        const auto distanceSquaredX = displacementX * displacementX;

        const auto distanceSquaredY = displacementY * displacementY;

        const auto distanceSquaredZ = displacementZ * displacementZ;


        const auto distanceSquared = distanceSquaredX + distanceSquaredY + distanceSquaredZ;


        const auto invDistSquared = 1. / distanceSquared;

        const auto lj2 = sigmaSquared * invDistSquared;

        const auto lj6 = lj2 * lj2 * lj2;

        const auto lj12 = lj6 * lj6;

        const auto lj12m6 = lj12 - lj6;

        const auto scalarMultiple = siteMask[neighborSite] ? epsilon24 * (lj12 + lj12m6) * invDistSquared : 0.;


        // calculate forces

        const auto forceX = scalarMultiple * displacementX;

        const auto forceY = scalarMultiple * displacementY;

        const auto forceZ = scalarMultiple * displacementZ;


        forceSumX += forceX;

        forceSumY += forceY;

        forceSumZ += forceZ;


        torqueSumX += rotatedPrimeSitePositionY * forceZ - rotatedPrimeSitePositionZ * forceY;

        torqueSumY += rotatedPrimeSitePositionZ * forceX - rotatedPrimeSitePositionX * forceZ;

        torqueSumZ += rotatedPrimeSitePositionX * forceY - rotatedPrimeSitePositionY * forceX;


        // N3L

        if (newton3) {

          siteForceX[neighborSite] -= forceX;

          siteForceY[neighborSite] -= forceY;

          siteForceZ[neighborSite] -= forceZ;

        }


        // calculate globals

        if constexpr (calculateGlobals) {

          const auto potentialEnergy6 = siteMask[neighborSite] ? (epsilon24 * lj12m6 + shift6) : 0.;

          const auto virialX = displacementX * forceX;

          const auto virialY = displacementY * forceY;

          const auto virialZ = displacementZ * forceZ;


          // We add 6 times the potential energy for each owned particle. The total sum is corrected in endTraversal().

          const auto ownershipFactor =

              newton3 ? (ownedStatePrime == autopas::OwnershipState::owned ? 1. : 0.) + (isNeighborSiteOwned ? 1. : 0.)

                      : (ownedStatePrime == autopas::OwnershipState::owned ? 1. : 0.);

          potentialEnergySum += potentialEnergy6 * ownershipFactor;

          virialSumX += virialX * ownershipFactor;

          virialSumY += virialY * ownershipFactor;

          virialSumZ += virialZ * ownershipFactor;

        }

      }

    }

    // Add forces to prime mol

    fxptr[indexPrime] += forceSumX;

    fyptr[indexPrime] += forceSumY;

    fzptr[indexPrime] += forceSumZ;

    txptr[indexPrime] += torqueSumX;

    typtr[indexPrime] += torqueSumY;

    tzptr[indexPrime] += torqueSumZ;


    // Reduce forces on individual neighbor sites to molecular forces & torques if newton3=true

    if constexpr (newton3) {

      if constexpr (useMixing) {

        size_t siteIndex = 0;

        for (size_t neighborMol = 0; neighborMol < neighborListSize; ++neighborMol) {

          const auto neighborMolIndex = neighborList[neighborMol];

          const auto rotatedSitePositions = autopas::utils::quaternion::rotateVectorOfPositions(

              {q0ptr[neighborMolIndex], q1ptr[neighborMolIndex], q2ptr[neighborMolIndex], q3ptr[neighborMolIndex]},

              _PPLibrary->getSitePositions(typeptr[neighborMolIndex]));

          for (size_t site = 0; site < _PPLibrary->getNumSites(typeptr[neighborMolIndex]); ++site) {

            fxptr[neighborMolIndex] += siteForceX[siteIndex];

            fyptr[neighborMolIndex] += siteForceY[siteIndex];

            fzptr[neighborMolIndex] += siteForceZ[siteIndex];

            txptr[neighborMolIndex] += rotatedSitePositions[site][1] * siteForceZ[siteIndex] -

                                       rotatedSitePositions[site][2] * siteForceY[siteIndex];

            typtr[neighborMolIndex] += rotatedSitePositions[site][2] * siteForceX[siteIndex] -

                                       rotatedSitePositions[site][0] * siteForceZ[siteIndex];

            tzptr[neighborMolIndex] += rotatedSitePositions[site][0] * siteForceY[siteIndex] -

                                       rotatedSitePositions[site][1] * siteForceX[siteIndex];

            ++siteIndex;

          }

        }

      } else {

        size_t siteIndex = 0;

        for (size_t neighborMol = 0; neighborMol < neighborListSize; ++neighborMol) {

          const auto neighborMolIndex = neighborList[neighborMol];

          const auto rotatedSitePositions = autopas::utils::quaternion::rotateVectorOfPositions(

              {q0ptr[neighborMolIndex], q1ptr[neighborMolIndex], q2ptr[neighborMolIndex], q3ptr[neighborMolIndex]},

              const_unrotatedSitePositions);

          for (size_t site = 0; site < const_unrotatedSitePositions.size(); ++site) {

            fxptr[neighborMolIndex] += siteForceX[siteIndex];

            fyptr[neighborMolIndex] += siteForceY[siteIndex];

            fzptr[neighborMolIndex] += siteForceZ[siteIndex];

            txptr[neighborMolIndex] += rotatedSitePositions[site][1] * siteForceZ[siteIndex] -

                                       rotatedSitePositions[site][2] * siteForceY[siteIndex];

            typtr[neighborMolIndex] += rotatedSitePositions[site][2] * siteForceX[siteIndex] -

                                       rotatedSitePositions[site][0] * siteForceZ[siteIndex];

            txptr[neighborMolIndex] += rotatedSitePositions[site][0] * siteForceY[siteIndex] -

                                       rotatedSitePositions[site][1] * siteForceX[siteIndex];

            ++siteIndex;

          }

        }

      }

    }


    if constexpr (calculateGlobals) {

      const auto threadNum = autopas::autopas_get_thread_num();


      _aosThreadData[threadNum].potentialEnergySum += potentialEnergySum;

      _aosThreadData[threadNum].virialSum[0] += virialSumX;

      _aosThreadData[threadNum].virialSum[1] += virialSumY;

      _aosThreadData[threadNum].virialSum[2] += virialSumZ;

    }

  }


  class AoSThreadData {

   public:

    AoSThreadData() : virialSum{0., 0., 0.}, potentialEnergySum{0.}, __remainingTo64{} {}

    void setZero() {

      virialSum = {0., 0., 0.};

      potentialEnergySum = 0.;

    }


    // variables

    std::array<double, 3> virialSum;

    double potentialEnergySum;


   private:

    // dummy parameter to get the right size (64 bytes)

    double __remainingTo64[(64 - 4 * sizeof(double)) / sizeof(double)];

  };


  static_assert(sizeof(AoSThreadData) % 64 == 0, "AoSThreadData has wrong size (should be multiple of 64)");


  std::vector<AoSThreadData> _aosThreadData;

};

}  // namespace mdLib

AlignedAllocator.h

ArrayMath.h

ExceptionHandler.h

AutoPasLog
#define AutoPasLog(lvl, fmt,...)
Macro for logging providing common meta information without filename.
Definition: Logger.h:24

Math.h

MoleculeLJ.h

MultisiteMoleculeLJ.h

OwnershipState.h

PairwiseFunctor.h

ParticlePropertiesLibrary.h

Quaternion.h

SoA.h

StaticBoolSelector.h

WrapOpenMP.h

ParticlePropertiesLibrary
This class stores the (physical) properties of molecule types, and, in the case of multi-site molecul...
Definition: ParticlePropertiesLibrary.h:28

ParticlePropertiesLibrary::getMixingShift6
floatType getMixingShift6(intType i, intType j) const
Returns precomputed mixed shift * 6 for one pair of site types.
Definition: ParticlePropertiesLibrary.h:262

ParticlePropertiesLibrary::getNumSites
intType getNumSites(intType i) const
Get number of sites of a multi-site molecule.
Definition: ParticlePropertiesLibrary.h:554

ParticlePropertiesLibrary::getSitePositions
std::vector< std::array< floatType, 3 > > getSitePositions(intType i) const
Get relative site positions to a multi-site molecule's center-of-mass.
Definition: ParticlePropertiesLibrary.h:517

ParticlePropertiesLibrary::getMixingSigmaSquared
floatType getMixingSigmaSquared(intType i, intType j) const
Returns precomputed mixed squared sigma for one pair of site types.
Definition: ParticlePropertiesLibrary.h:252

ParticlePropertiesLibrary::getMixing24Epsilon
floatType getMixing24Epsilon(intType i, intType j) const
Returns the precomputed mixed epsilon * 24.
Definition: ParticlePropertiesLibrary.h:224

ParticlePropertiesLibrary::calcShift6
static double calcShift6(double epsilon24, double sigmaSquared, double cutoffSquared)
Calculate the shift multiplied 6 of the lennard jones potential from given cutoff,...
Definition: ParticlePropertiesLibrary.h:576

ParticlePropertiesLibrary::getLJMixingData
auto getLJMixingData(intType i, intType j) const
Get complete mixing data for one pair of LJ site types.
Definition: ParticlePropertiesLibrary.h:234

ParticlePropertiesLibrary::getSiteTypes
std::vector< intType > getSiteTypes(intType i) const
Get site types of a multi-site molecule.
Definition: ParticlePropertiesLibrary.h:528

autopas::AlignedAllocator
AlignedAllocator class.
Definition: AlignedAllocator.h:29

autopas::PairwiseFunctor
PairwiseFunctor class.
Definition: PairwiseFunctor.h:31

autopas::PairwiseFunctor::PairwiseFunctor
PairwiseFunctor(double cutoff)
Constructor.
Definition: PairwiseFunctor.h:42

autopas::SoAView
View on a fixed part of a SoA between a start index and an end index.
Definition: SoAView.h:23

autopas::SoAView::size
size_t size() const
Returns the number of particles in the view.
Definition: SoAView.h:83

autopas::utils::ExceptionHandler::AutoPasException
Default exception class for autopas exceptions.
Definition: ExceptionHandler.h:116

mdLib::LJMultisiteFunctor
A functor to handle Lennard-Jones interactions between two Multisite Molecules.
Definition: LJMultisiteFunctor.h:49

mdLib::LJMultisiteFunctor::SoAFunctorSingle
void SoAFunctorSingle(autopas::SoAView< SoAArraysType > soa, bool newton3) final
PairwiseFunctor for structure of arrays (SoA)
Definition: LJMultisiteFunctor.h:281

mdLib::LJMultisiteFunctor::AoSFunctor
void AoSFunctor(Particle_T &particleA, Particle_T &particleB, bool newton3) final
Functor for arrays of structures (AoS).
Definition: LJMultisiteFunctor.h:181

mdLib::LJMultisiteFunctor::LJMultisiteFunctor
LJMultisiteFunctor()=delete
Delete Default constructor.

mdLib::LJMultisiteFunctor::isRelevantForTuning
bool isRelevantForTuning() final
Specifies whether the functor should be considered for the auto-tuning process.
Definition: LJMultisiteFunctor.h:164

mdLib::LJMultisiteFunctor::LJMultisiteFunctor
LJMultisiteFunctor(double cutoff, ParticlePropertiesLibrary< double, size_t > &particlePropertiesLibrary)
Constructor for Functor with particle mixing enabled.
Definition: LJMultisiteFunctor.h:154

mdLib::LJMultisiteFunctor::LJMultisiteFunctor
LJMultisiteFunctor(double cutoff)
Constructor for Functor with particle mixing disabled.
Definition: LJMultisiteFunctor.h:140

mdLib::LJMultisiteFunctor::initTraversal
void initTraversal() final
Reset the global values.
Definition: LJMultisiteFunctor.h:712

mdLib::LJMultisiteFunctor::getMixing
static constexpr bool getMixing()
Definition: LJMultisiteFunctor.h:683

mdLib::LJMultisiteFunctor::getComputedAttr
static constexpr auto getComputedAttr()
Get attributes computed by this functor.
Definition: LJMultisiteFunctor.h:674

mdLib::LJMultisiteFunctor::setParticleProperties
void setParticleProperties(SoAFloatPrecision epsilon24, SoAFloatPrecision sigmaSquared, std::vector< std::array< SoAFloatPrecision, 3 > > sitePositionsLJ)
Sets the molecule properties constants for this functor.
Definition: LJMultisiteFunctor.h:629

mdLib::LJMultisiteFunctor::allowsNewton3
bool allowsNewton3() final
Specifies whether the functor is capable of Newton3-like functors.
Definition: LJMultisiteFunctor.h:166

mdLib::LJMultisiteFunctor::allowsNonNewton3
bool allowsNonNewton3() final
Specifies whether the functor is capable of non-Newton3-like functors.
Definition: LJMultisiteFunctor.h:170

mdLib::LJMultisiteFunctor::getPotentialEnergy
double getPotentialEnergy()
Get the potential energy.
Definition: LJMultisiteFunctor.h:757

mdLib::LJMultisiteFunctor::getName
std::string getName() final
Returns name of functor.
Definition: LJMultisiteFunctor.h:162

mdLib::LJMultisiteFunctor::getNumFlopsPerKernelCall
unsigned long getNumFlopsPerKernelCall(size_t molAType, size_t molBType, bool newton3)
Get the number of flops used per kernel call - i.e.
Definition: LJMultisiteFunctor.h:695

mdLib::LJMultisiteFunctor::SoAFunctorPair
void SoAFunctorPair(autopas::SoAView< SoAArraysType > soa1, autopas::SoAView< SoAArraysType > soa2, const bool newton3) final
PairwiseFunctor for structure of arrays (SoA)
Definition: LJMultisiteFunctor.h:596

mdLib::LJMultisiteFunctor::getNeededAttr
static constexpr auto getNeededAttr(std::false_type)
Get attributes needed for computation without N3 optimization.
Definition: LJMultisiteFunctor.h:659

mdLib::LJMultisiteFunctor::SoAFunctorVerlet
void SoAFunctorVerlet(autopas::SoAView< SoAArraysType > soa, const size_t indexFirst, const std::vector< size_t, autopas::AlignedAllocator< size_t > > &neighborList, bool newton3) final
PairwiseFunctor for structure of arrays (SoA) for neighbor lists.
Definition: LJMultisiteFunctor.h:610

mdLib::LJMultisiteFunctor::getVirial
double getVirial()
Get the virial.
Definition: LJMultisiteFunctor.h:775

mdLib::LJMultisiteFunctor::endTraversal
void endTraversal(bool newton3) final
Postprocesses global values, e.g.
Definition: LJMultisiteFunctor.h:725

mdLib::LJMultisiteFunctor::getNeededAttr
static constexpr auto getNeededAttr()
Get attributes needed for computation.
Definition: LJMultisiteFunctor.h:644

inBox.h

autopas::utils::ArrayMath::dot
constexpr T dot(const std::array< T, SIZE > &a, const std::array< T, SIZE > &b)
Generates the dot product of two arrays.
Definition: ArrayMath.h:233

autopas::utils::ArrayMath::mulScalar
constexpr std::array< T, SIZE > mulScalar(const std::array< T, SIZE > &a, T s)
Multiplies a scalar s to each element of array a and returns the result.
Definition: ArrayMath.h:181

autopas::utils::ArrayMath::sub
constexpr std::array< T, SIZE > sub(const std::array< T, SIZE > &a, const std::array< T, SIZE > &b)
Subtracts array b from array a and returns the result.
Definition: ArrayMath.h:45

autopas::utils::ArrayMath::add
constexpr std::array< T, SIZE > add(const std::array< T, SIZE > &a, const std::array< T, SIZE > &b)
Adds two arrays, returns the result.
Definition: ArrayMath.h:28

autopas::utils::ArrayMath::cross
constexpr std::array< T, 3 > cross(const std::array< T, 3 > &a, const std::array< T, 3 > &b)
Generates the cross product of two arrays of 3 floats.
Definition: ArrayMath.h:249

autopas::utils::quaternion::rotateVectorOfPositions
std::vector< std::array< double, 3 > > rotateVectorOfPositions(const std::array< double, 4 > &q, const std::vector< std::array< double, 3 > > &positionVector)
Rotates a std::vector of 3D positions.
Definition: Quaternion.cpp:13

autopas
This is the main namespace of AutoPas.
Definition: AutoPasDecl.h:34

autopas::autopas_get_max_threads
int autopas_get_max_threads()
Dummy for omp_get_max_threads() when no OpenMP is available.
Definition: WrapOpenMP.h:144

autopas::OwnershipState
OwnershipState
Enum that specifies the state of ownership.
Definition: OwnershipState.h:20

autopas::OwnershipState::dummy
@ dummy
Dummy or deleted state, a particle with this state is not an actual particle!

autopas::OwnershipState::owned
@ owned
Owned state, a particle with this state is an actual particle and owned by the current AutoPas object...

autopas::FunctorN3Modes
FunctorN3Modes
Newton 3 modes for the Functor.
Definition: Functor.h:23

autopas::autopas_get_thread_num
int autopas_get_thread_num()
Dummy for omp_set_lock() when no OpenMP is available.
Definition: WrapOpenMP.h:132