ENH: Using site lists in force calculation.

applications/solvers/discreteMethods/molecularDynamics/mdFoam/Make/options

@@ -1,7 +1,6 @@
 EXE_INC = \
     -I$(LIB_SRC)/lagrangian/molecularDynamics/molecule/lnInclude \
     -I$(LIB_SRC)/lagrangian/molecularDynamics/potential/lnInclude \
-    -I$(LIB_SRC)/lagrangian/molecularDynamics/molecularMeasurements/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/lagrangian/basic/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude
@@ -11,6 +10,4 @@ EXE_LIBS = \
     -lfiniteVolume \
     -llagrangian \
     -lmolecule \
-    -lpotential \
+    -lpotential
-    -lmolecularMeasurements
-

src/lagrangian/molecularDynamics/molecule/polyatomic/polyatomic/polyatomic.C

@@ -156,13 +156,15 @@ bool Foam::polyatomic::move(polyatomic::trackingData& td, const scalar trackTime
 
         tau_ = vector::zero;
 
-        forAll(siteForces_, s)
+        forAll(siteForces_, sI)
         {
-            const vector& f = siteForces_[s];
+            const vector& f = siteForces_[sI];
 
             a_ += f/m;
 
-            tau_ += (constProps.siteReferencePositions()[s] ^ (Q_.T() & f));
+            tau_ +=
+                constProps.sites()[sI].siteReferencePosition()
+              ^ (Q_.T() & f);
         }
 
         v_ += 0.5*trackTime*a_;
@@ -231,7 +233,12 @@ void Foam::polyatomic::transformProperties(const vector& separation)
 
 void Foam::polyatomic::setSitePositions(const constantProperties& constProps)
 {
-    sitePositions_ = position_ + (Q_ & constProps.siteReferencePositions());
+    forAll(constProps.sites(), sI)
+    {
+        sitePositions_[sI] =
+            position_
+          + (Q_ & constProps.sites()[sI].siteReferencePosition());
+    }
 }
 
 

src/lagrangian/molecularDynamics/molecule/polyatomic/polyatomic/polyatomic.H

@@ -41,6 +41,7 @@ SourceFiles
 #include "IOstream.H"
 #include "autoPtr.H"
 #include "diagTensor.H"
+#include "constPropSite.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -78,18 +79,16 @@ public:
     {
         // Private data
 
-            Field<vector> siteReferencePositions_;
+            //- Sites of mass, charge or interaction
+            List<constPropSite> sites_;
 
-            List<scalar> siteMasses_;
+            //- Which sites require electrostatic interactions
+            List<label> electrostaticSites_;
 
-            List<scalar> siteCharges_;
+            //- Which sites require pair interactions
-
+            List<label> pairPotSites_;
-            List<label> siteIds_;
-
-            List<bool> pairPotentialSites_;
-
-            List<bool> electrostaticSites_;
 
+            //- Moment of intertia (in principal axis configiration)
             diagTensor momentOfInertia_;
 
             scalar mass_;
@@ -97,8 +96,6 @@ public:
 
         // Private Member Functions
 
-            void checkSiteListSizes() const;
-
             void setInteracionSiteBools
             (
                 const List<word>& siteIds,
@@ -113,27 +110,19 @@ public:
         inline constantProperties();
 
         //- Construct from dictionary
-        inline constantProperties(const dictionary& dict);
+        inline constantProperties
+        (
+            const dictionary& dict,
+            const labelList& siteIds
+        );
 
         // Member functions
 
-        inline const Field<vector>& siteReferencePositions() const;
+        inline const List<constPropSite>& sites() const;
 
-        inline const List<scalar>& siteMasses() const;
+        inline const List<label>& pairPotSites() const;
 
-        inline const List<scalar>& siteCharges() const;
+        inline const List<label>& electrostaticSites() const;
-
-        inline const List<label>& siteIds() const;
-
-        inline List<label>& siteIds();
-
-        inline const List<bool>& pairPotentialSites() const;
-
-        inline bool pairPotentialSite(label sId) const;
-
-        inline const List<bool>& electrostaticSites() const;
-
-        inline bool electrostaticSite(label sId) const;
 
         inline const diagTensor& momentOfInertia() const;
 
@@ -181,14 +170,22 @@ private:
 
     // Private data
 
+        //- Orientation, stored as the rotation tensor to transform
+        //  from the polyatomic to the global reference frame, i.e.:
+        //  globalVector = Q_ & bodyLocalVector
+        //  bodyLocalVector = Q_.T() & globalVector
         tensor Q_;
 
+        // Linear velocity of polyatomic
         vector v_;
 
+        // Total linear acceleration of polyatomic
         vector a_;
 
+        //- Angular momentum of polyatomic, in body local reference frame
         vector pi_;
 
+        //- Total torque on polyatomic, in body local reference frame
         vector tau_;
 
         vector specialPosition_;
@@ -198,6 +195,9 @@ private:
         // - r_ij f_ij, stress dyad
         tensor rf_;
 
+//         // - r_ij outer product f_ij: virial contribution
+//         tensor rDotf_;
+
         label special_;
 
         label id_;

src/lagrangian/molecularDynamics/molecule/polyatomic/polyatomic/polyatomicI.H

@@ -27,12 +27,9 @@ License
 
 inline Foam::polyatomic::constantProperties::constantProperties()
 :
-    siteReferencePositions_(Field<vector>(0)),
+    sites_(),
-    siteMasses_(List<scalar>(0)),
+    electrostaticSites_(),
-    siteCharges_(List<scalar>(0)),
+    pairPotSites_(),
-    siteIds_(List<label>(0)),
-    pairPotentialSites_(List<bool>(false)),
-    electrostaticSites_(List<bool>(false)),
     momentOfInertia_(diagTensor(0, 0, 0)),
     mass_(0)
 {}
@@ -40,47 +37,117 @@ inline Foam::polyatomic::constantProperties::constantProperties()
 
 inline Foam::polyatomic::constantProperties::constantProperties
 (
-    const dictionary& dict
+    const dictionary& dict,
+    const labelList& siteIds
 )
 :
-    siteReferencePositions_(dict.lookup("siteReferencePositions")),
+    sites_(),
-    siteMasses_(dict.lookup("siteMasses")),
-    siteCharges_(dict.lookup("siteCharges")),
-    siteIds_(List<word>(dict.lookup("siteIds")).size(), -1),
-    pairPotentialSites_(),
     electrostaticSites_(),
+    pairPotSites_(),
     momentOfInertia_(),
-    mass_()
+    mass_(0)
 {
-    checkSiteListSizes();
+    List<word> siteIdNames = dict.lookup("siteIds");
 
-    setInteracionSiteBools
+    List<vector> siteReferencePositions(dict.lookup("siteReferencePositions"));
+
+    List<scalar> siteMasses(dict.lookup("siteMasses"));
+
+    List<scalar> siteCharges(dict.lookup("siteCharges"));
+
+    List<word> pairPotentialIds(dict.lookup("pairPotentialSiteIds"));
+
+    if
     (
-        List<word>(dict.lookup("siteIds")),
+        ( siteIdNames.size() != sites_.size() )
-        List<word>(dict.lookup("pairPotentialSiteIds"))
+     || ( siteReferencePositions.size() != sites_.size() )
-    );
+     || ( siteCharges.size() != sites_.size() )
+    )
+    {
+            FatalErrorIn
+            (
+                "Foam::polyatomic::constantProperties::constantProperties"
+                "("
+                    "const dictionary& dict"
+                ")"
+            )
+            << "Sizes of site id, charge and "
+            << "referencePositions are not the same: " << sites_.size()
+            << nl << abort(FatalError);
+    }
 
-    mass_ = sum(siteMasses_);
+    sites_.setSize(siteIdNames.size());
 
-    vector centreOfMass(vector::zero);
+    electrostaticSites_.setSize(sites_.size(), -1);
+    pairPotSites_.setSize(sites_.size(), -1);
+
+    label pairI = 0;
+    label electrostaticI = 0;
+
+    forAll(sites_, sI)
+    {
+        sites_[sI] = constPropSite
+        (
+            siteReferencePositions[sI],
+            siteMasses[sI],
+            siteCharges[sI],
+            siteIds[sI],
+            siteIdNames[sI],
+            (findIndex(pairPotentialIds, siteIdNames[sI]) != -1),   // pair
+            (mag(siteCharges[sI]) > VSMALL)                         // charge
+        );
+
+        if (sites_[sI].pairPotentialSite())
+        {
+            pairPotSites_[pairI++] = sI;
+        }
+
+        if (sites_[sI].electrostaticSite())
+        {
+            electrostaticSites_[electrostaticI++] = sI;
+        }
+
+        if (sites_[sI].pairPotentialSite() && !sites_[sI].electrostaticSite())
+        {
+            WarningIn
+            (
+                "Foam::polyatomic::constantProperties::constantProperties"
+                "("
+                    "const dictionary& dict"
+                ")"
+            )
+                << siteIdNames[sI] << " is a non-interacting site." << endl;
+        }
+    }
+
+    pairPotSites_.setSize(pairI);
+    electrostaticSites_.setSize(electrostaticI);
 
     // Calculate the centre of mass of the body and subtract it from each
     // position
 
-    forAll(siteReferencePositions_, i)
+    vector centreOfMass(vector::zero);
+
+    forAll(sites_, sI)
     {
-        centreOfMass += siteReferencePositions_[i]*siteMasses_[i];
+        mass_ += sites_[sI].siteMass();
+
+        centreOfMass +=
+            sites_[sI].siteReferencePosition()*sites_[sI].siteMass();
     }
 
     centreOfMass /= mass_;
 
-    siteReferencePositions_ -= centreOfMass;
+    forAll(sites_, sI)
+    {
+        sites_[sI].siteReferencePosition() -= centreOfMass;
+    }
 
-    if (siteIds_.size() == 1)
+    if (sites_.size() == 1)
     {
         // Single site polyatomic - no rotational motion.
 
-        siteReferencePositions_[0] = vector::zero;
+        sites_[0].siteReferencePosition() = vector::zero;
 
         momentOfInertia_ = diagTensor(-1, -1, -1);
     }
@@ -90,36 +157,46 @@ inline Foam::polyatomic::constantProperties::constantProperties
 
         Info<< nl << "Linear polyatomic." << endl;
 
-        vector dir = siteReferencePositions_[1] - siteReferencePositions_[0];
+        vector dir =
+            sites_[1].siteReferencePosition()
+          - sites_[0].siteReferencePosition();
 
         dir /= mag(dir);
 
         tensor Q = rotationTensor(dir, vector(1,0,0));
 
-        siteReferencePositions_ = (Q & siteReferencePositions_);
+        forAll(sites_, sI)
+        {
+            sites_[sI].siteReferencePosition() =
+                (Q & sites_[sI].siteReferencePosition());
+        }
 
         // The rotation was around the centre of mass but remove any
         // components that have crept in due to floating point errors
 
         centreOfMass = vector::zero;
 
-        forAll(siteReferencePositions_, i)
+        forAll(sites_, sI)
         {
-            centreOfMass += siteReferencePositions_[i]*siteMasses_[i];
+            centreOfMass +=
+                sites_[sI].siteReferencePosition()*sites_[sI].siteMass();
         }
 
         centreOfMass /= mass_;
 
-        siteReferencePositions_ -= centreOfMass;
+        forAll(sites_, sI)
+        {
+            sites_[sI].siteReferencePosition() -= centreOfMass;
+        }
 
         diagTensor momOfInertia = diagTensor::zero;
 
-        forAll(siteReferencePositions_, i)
+        forAll(sites_, sI)
         {
-            const vector& p(siteReferencePositions_[i]);
+            const vector& p(sites_[sI].siteReferencePosition());
 
             momOfInertia +=
-                siteMasses_[i]*diagTensor(0, p.x()*p.x(), p.x()*p.x());
+                sites_[sI].siteMass()*diagTensor(0, p.x()*p.x(), p.x()*p.x());
         }
 
         momentOfInertia_ = diagTensor
@@ -137,11 +214,11 @@ inline Foam::polyatomic::constantProperties::constantProperties
 
         tensor momOfInertia(tensor::zero);
 
-        forAll(siteReferencePositions_, i)
+        forAll(sites_, sI)
         {
-            const vector& p(siteReferencePositions_[i]);
+            const vector& p(sites_[sI].siteReferencePosition());
 
-            momOfInertia += siteMasses_[i]*tensor
+            momOfInertia += sites_[sI].siteMass()*tensor
             (
                 p.y()*p.y() + p.z()*p.z(), -p.x()*p.y(), -p.x()*p.z(),
                 -p.y()*p.x(), p.x()*p.x() + p.z()*p.z(), -p.y()*p.z(),
@@ -152,8 +229,8 @@ inline Foam::polyatomic::constantProperties::constantProperties
         if (eigenValues(momOfInertia).x() < VSMALL)
         {
             FatalErrorIn("polyatomic::constantProperties::constantProperties")
-                << "An eigenvalue of the inertia tensor is zero, the polyatomic "
+                << "An eigenvalue of the inertia tensor is zero, the "
-                << dict.name() << " is not a valid 6DOF shape."
+                << "polyatomic " << dict.name() << " is not a valid 6DOF shape."
                 << nl << abort(FatalError);
         }
 
@@ -172,7 +249,11 @@ inline Foam::polyatomic::constantProperties::constantProperties
 
         // Transform the site positions
 
-        siteReferencePositions_ = (Q & siteReferencePositions_);
+        forAll(sites_, sI)
+        {
+            sites_[sI].siteReferencePosition() =
+                (Q & sites_[sI].siteReferencePosition());
+        }
 
         // Recalculating the moment of inertia with the new site positions
 
@@ -181,25 +262,29 @@ inline Foam::polyatomic::constantProperties::constantProperties
 
         centreOfMass = vector::zero;
 
-        forAll(siteReferencePositions_, i)
+        forAll(sites_, sI)
         {
-            centreOfMass += siteReferencePositions_[i]*siteMasses_[i];
+            centreOfMass +=
+                sites_[sI].siteReferencePosition()*sites_[sI].siteMass();
         }
 
         centreOfMass /= mass_;
 
-        siteReferencePositions_ -= centreOfMass;
+        forAll(sites_, sI)
+        {
+            sites_[sI].siteReferencePosition() -= centreOfMass;
+        }
 
         // Calculate the moment of inertia in the principle component
         // reference frame
 
         momOfInertia = tensor::zero;
 
-        forAll(siteReferencePositions_, i)
+        forAll(sites_, sI)
         {
-            const vector& p(siteReferencePositions_[i]);
+            const vector& p(sites_[sI].siteReferencePosition());
 
-            momOfInertia += siteMasses_[i]*tensor
+            momOfInertia += sites_[sI].siteMass()*tensor
             (
                 p.y()*p.y() + p.z()*p.z(), -p.x()*p.y(), -p.x()*p.z(),
                 -p.y()*p.x(), p.x()*p.x() + p.z()*p.z(), -p.y()*p.z(),
@@ -256,62 +341,29 @@ inline Foam::polyatomic::polyatomic
 
 // * * * * * * * constantProperties Private Member Functions * * * * * * * * //
 
-inline void Foam::polyatomic::constantProperties::checkSiteListSizes() const
-{
-    if
-    (
-        siteIds_.size() != siteReferencePositions_.size()
-     || siteIds_.size() != siteCharges_.size()
-    )
-    {
-        FatalErrorIn("polyatomic::constantProperties::checkSiteListSizes")
-            << "Sizes of site id, charge and "
-            << "referencePositions are not the same. "
-            << nl << abort(FatalError);
-    }
-}
-
-
-inline void Foam::polyatomic::constantProperties::setInteracionSiteBools
-(
-    const List<word>& siteIds,
-    const List<word>& pairPotSiteIds
-)
-{
-    pairPotentialSites_.setSize(siteIds_.size());
-
-    electrostaticSites_.setSize(siteIds_.size());
-
-    forAll(siteIds_, i)
-    {
-        const word& id(siteIds[i]);
-
-        pairPotentialSites_[i] = (findIndex(pairPotSiteIds, id) > -1);
-
-        electrostaticSites_[i] = (mag(siteCharges_[i]) > VSMALL);
-    }
-}
-
-
 inline bool Foam::polyatomic::constantProperties::linearMoleculeTest() const
 {
-    if (siteIds_.size() == 2)
+    if (sites_.size() == 2)
     {
         return true;
     }
 
-    vector refDir = siteReferencePositions_[1] - siteReferencePositions_[0];
+    vector refDir =
+        sites_[1].siteReferencePosition()
+      - sites_[0].siteReferencePosition();
 
     refDir /= mag(refDir);
 
     for
     (
         label i = 2;
-        i < siteReferencePositions_.size();
+        i < sites_.size();
         i++
     )
     {
-        vector dir = siteReferencePositions_[i] - siteReferencePositions_[i-1];
+        vector dir =
+            sites_[i].siteReferencePosition()
+          - sites_[i - 1].siteReferencePosition();
 
         dir /= mag(dir);
 
@@ -327,93 +379,27 @@ inline bool Foam::polyatomic::constantProperties::linearMoleculeTest() const
 
 // * * * * * * * constantProperties Member Functions * * * * * * * * * * * * //
 
-inline const Foam::Field<Foam::vector>&
+inline const Foam::List<Foam::constPropSite>&
-Foam::polyatomic::constantProperties::siteReferencePositions() const
+Foam::polyatomic::constantProperties::sites() const
 {
-    return siteReferencePositions_;
+    return sites_;
-}
-
-
-inline const Foam::List<Foam::scalar>&
-Foam::polyatomic::constantProperties::siteMasses() const
-{
-    return siteMasses_;
-}
-
-
-inline const Foam::List<Foam::scalar>&
-Foam::polyatomic::constantProperties::siteCharges() const
-{
-    return siteCharges_;
 }
 
 
 inline const Foam::List<Foam::label>&
-Foam::polyatomic::constantProperties::siteIds() const
+Foam::polyatomic::constantProperties::pairPotSites() const
 {
-    return siteIds_;
+    return pairPotSites_;
 }
 
 
-inline Foam::List<Foam::label>&
+inline const Foam::List<Foam::label>&
-Foam::polyatomic::constantProperties::siteIds()
-{
-    return siteIds_;
-}
-
-
-inline const Foam::List<bool>&
-Foam::polyatomic::constantProperties::pairPotentialSites() const
-{
-    return pairPotentialSites_;
-}
-
-
-inline bool Foam::polyatomic::constantProperties::pairPotentialSite
-(
-    label sId
-) const
-{
-    label s = findIndex(siteIds_, sId);
-
-    if (s == -1)
-    {
-        FatalErrorIn("polyatomicI.H") << nl
-            << sId << " site not found."
-            << nl << abort(FatalError);
-    }
-
-    return pairPotentialSites_[s];
-}
-
-
-inline const Foam::List<bool>&
 Foam::polyatomic::constantProperties::electrostaticSites() const
 {
     return electrostaticSites_;
 }
 
 
-inline bool Foam::polyatomic::constantProperties::electrostaticSite
-(
-    label sId
-) const
-{
-    label s = findIndex(siteIds_, sId);
-
-    if (s == -1)
-    {
-        FatalErrorIn
-        (
-            "polyatomic::constantProperties::electrostaticSite(label)"
-        )   << sId << " site not found."
-            << nl << abort(FatalError);
-    }
-
-    return electrostaticSites_[s];
-}
-
-
 inline const Foam::diagTensor&
 Foam::polyatomic::constantProperties::momentOfInertia() const
 {
@@ -433,7 +419,8 @@ inline bool Foam::polyatomic::constantProperties::pointMolecule() const
 }
 
 
-inline Foam::label Foam::polyatomic::constantProperties::degreesOfFreedom() const
+inline Foam::label
+Foam::polyatomic::constantProperties::degreesOfFreedom() const
 {
     if (linearMolecule())
     {
@@ -458,8 +445,7 @@ inline Foam::scalar Foam::polyatomic::constantProperties::mass() const
 
 inline Foam::label Foam::polyatomic::constantProperties::nSites() const
 {
-    return siteIds_.size();
+    return sites_.size();
-
 }
 
 

src/lagrangian/molecularDynamics/molecule/polyatomic/polyatomicCloud/polyatomicCloud.C

@@ -64,7 +64,8 @@ void Foam::polyatomicCloud::buildConstProps()
     forAll(idList, i)
     {
         const word& id = idList[i];
-        const dictionary& molDict = polyatomicPropertiesDict.subDict(id);
+
+        const dictionary& molDict(polyatomicPropertiesDict.subDict(id));
 
         List<word> siteIdNames = molDict.lookup("siteIds");
 
@@ -78,17 +79,23 @@ void Foam::polyatomicCloud::buildConstProps()
 
             if (siteIds[sI] == -1)
             {
-                FatalErrorIn("polyatomicCloud::buildConstProps()")
+                FatalErrorIn
+                (
+                    "Foam::polyatomic::constantProperties::constantProperties"
+                    "("
+                    "const dictionary& dict"
+                    ")"
+                )
                     << siteId << " site not found."
                     << nl << abort(FatalError);
             }
         }
 
-        polyatomic::constantProperties& constProp = constPropList_[i];
+        constPropList_[i] = polyatomic::constantProperties
-
+        (
-        constProp = polyatomic::constantProperties(molDict);
+            molDict,
-
+            siteIds
-        constProp.siteIds() = siteIds;
+        );
     }
 }
 
@@ -1237,7 +1244,7 @@ void Foam::polyatomicCloud::writeXYZ(const fileName& fName) const
         {
             const point& sP = mol().sitePositions()[i];
 
-            os  << pot_.siteIdList()[cP.siteIds()[i]]
+            os  << pot_.siteIdList()[cP.sites()[i].siteId()]
                 << ' ' << sP.x()*1e10
                 << ' ' << sP.y()*1e10
                 << ' ' << sP.z()*1e10

src/lagrangian/molecularDynamics/molecule/polyatomic/polyatomicCloud/polyatomicCloudI.H

@@ -47,105 +47,102 @@ inline void Foam::polyatomicCloud::evaluatePair
 
     const polyatomic::constantProperties& constPropJ(constProps(idJ));
 
-    List<label> siteIdsI = constPropI.siteIds();
+    forAll(constPropI.pairPotSites(), pI)
-
-    List<label> siteIdsJ = constPropJ.siteIds();
-
-    List<bool> pairPotentialSitesI = constPropI.pairPotentialSites();
-
-    List<bool> electrostaticSitesI = constPropI.electrostaticSites();
-
-    List<bool> pairPotentialSitesJ = constPropJ.pairPotentialSites();
-
-    List<bool> electrostaticSitesJ = constPropJ.electrostaticSites();
-
-    forAll(siteIdsI, sI)
     {
-        label idsI(siteIdsI[sI]);
+        label sI = constPropI.pairPotSites()[pI];
 
-        forAll(siteIdsJ, sJ)
+        label idsI = constPropI.sites()[sI].siteId();
+
+        forAll(constPropJ.pairPotSites(), pJ)
         {
-            label idsJ(siteIdsJ[sJ]);
+            label sJ = constPropJ.pairPotSites()[pJ];
 
-            if (pairPotentialSitesI[sI] && pairPotentialSitesJ[sJ])
+            label idsJ = constPropJ.sites()[sJ].siteId();
-            {
-                vector rsIsJ =
-                    molI.sitePositions()[sI] - molJ.sitePositions()[sJ];
 
-                scalar rsIsJMagSq = magSqr(rsIsJ);
+            vector rsIsJ =
-
-                if (pairPot.rCutSqr(idsI, idsJ, rsIsJMagSq))
-                {
-                    scalar rsIsJMag = mag(rsIsJ);
-
-                    vector fsIsJ =
-                        (rsIsJ/rsIsJMag)
-                       *pairPot.force(idsI, idsJ, rsIsJMag);
-
-                    molI.siteForces()[sI] += fsIsJ;
-
-                    molJ.siteForces()[sJ] += -fsIsJ;
-
-                    scalar potentialEnergy
-                    (
-                        pairPot.energy(idsI, idsJ, rsIsJMag)
-                    );
-
-                    molI.potentialEnergy() += 0.5*potentialEnergy;
-
-                    molJ.potentialEnergy() += 0.5*potentialEnergy;
-
-                    vector rIJ = molI.position() - molJ.position();
-
-                    tensor virialContribution =
-                        (rsIsJ*fsIsJ)*(rsIsJ & rIJ)/rsIsJMagSq;
-
-                    molI.rf() += virialContribution;
-
-                    molJ.rf() += virialContribution;
-                }
-            }
-
-            if (electrostaticSitesI[sI] && electrostaticSitesJ[sJ])
-            {
-                vector rsIsJ =
                 molI.sitePositions()[sI] - molJ.sitePositions()[sJ];
 
-                scalar rsIsJMagSq = magSqr(rsIsJ);
+            scalar rsIsJMagSq = magSqr(rsIsJ);
 
-                if (rsIsJMagSq <= electrostatic.rCutSqr())
+            if (pairPot.rCutSqr(idsI, idsJ, rsIsJMagSq))
-                {
+            {
-                    scalar rsIsJMag = mag(rsIsJ);
+                scalar rsIsJMag = mag(rsIsJ);
 
-                    scalar chargeI = constPropI.siteCharges()[sI];
+                vector fsIsJ =
+                (rsIsJ/rsIsJMag)
+               *pairPot.force(idsI, idsJ, rsIsJMag);
 
-                    scalar chargeJ = constPropJ.siteCharges()[sJ];
+                molI.siteForces()[sI] += fsIsJ;
 
-                    vector fsIsJ =
+                molJ.siteForces()[sJ] += -fsIsJ;
-                        (rsIsJ/rsIsJMag)
-                       *chargeI*chargeJ*electrostatic.force(rsIsJMag);
 
-                    molI.siteForces()[sI] += fsIsJ;
+                scalar potentialEnergy
+                (
+                    pairPot.energy(idsI, idsJ, rsIsJMag)
+                );
 
-                    molJ.siteForces()[sJ] += -fsIsJ;
+                molI.potentialEnergy() += 0.5*potentialEnergy;
 
-                    scalar potentialEnergy =
+                molJ.potentialEnergy() += 0.5*potentialEnergy;
-                        chargeI*chargeJ
-                       *electrostatic.energy(rsIsJMag);
 
-                    molI.potentialEnergy() += 0.5*potentialEnergy;
+                vector rIJ = molI.position() - molJ.position();
 
-                    molJ.potentialEnergy() += 0.5*potentialEnergy;
+                tensor virialContribution =
+                (rsIsJ*fsIsJ)*(rsIsJ & rIJ)/rsIsJMagSq;
 
-                    vector rIJ = molI.position() - molJ.position();
+                molI.rf() += virialContribution;
 
-                    tensor virialContribution =
+                molJ.rf() += virialContribution;
-                        (rsIsJ*fsIsJ)*(rsIsJ & rIJ)/rsIsJMagSq;
+            }
+        }
+    }
 
-                    molI.rf() += virialContribution;
 
-                    molJ.rf() += virialContribution;
+    forAll(constPropI.electrostaticSites(), pI)
-                }
+    {
+        label sI = constPropI.electrostaticSites()[pI];
+
+        forAll(constPropJ.electrostaticSites(), pJ)
+        {
+            label sJ = constPropJ.electrostaticSites()[pJ];
+
+            vector rsIsJ =
+                molI.sitePositions()[sI] - molJ.sitePositions()[sJ];
+
+            scalar rsIsJMagSq = magSqr(rsIsJ);
+
+            if (rsIsJMagSq <= electrostatic.rCutSqr())
+            {
+                scalar rsIsJMag = mag(rsIsJ);
+
+                scalar chargeI = constPropI.sites()[sI].siteCharge();
+
+                scalar chargeJ = constPropJ.sites()[sJ].siteCharge();
+
+                vector fsIsJ =
+                    (rsIsJ/rsIsJMag)
+                   *chargeI*chargeJ*electrostatic.force(rsIsJMag);
+
+                molI.siteForces()[sI] += fsIsJ;
+
+                molJ.siteForces()[sJ] += -fsIsJ;
+
+                scalar potentialEnergy =
+                    chargeI*chargeJ
+                   *electrostatic.energy(rsIsJMag);
+
+                molI.potentialEnergy() += 0.5*potentialEnergy;
+
+                molJ.potentialEnergy() += 0.5*potentialEnergy;
+
+                vector rIJ = molI.position() - molJ.position();
+
+                tensor virialContribution =
+                    (rsIsJ*fsIsJ)*(rsIsJ & rIJ)/rsIsJMagSq;
+
+                molI.rf() += virialContribution;
+
+                molJ.rf() += virialContribution;
             }
         }
     }
@@ -170,123 +167,119 @@ inline bool Foam::polyatomicCloud::evaluatePotentialLimit
 
     const polyatomic::constantProperties& constPropJ(constProps(idJ));
 
-    List<label> siteIdsI = constPropI.siteIds();
+    forAll(constPropI.pairPotSites(), pI)
-
-    List<label> siteIdsJ = constPropJ.siteIds();
-
-    List<bool> pairPotentialSitesI = constPropI.pairPotentialSites();
-
-    List<bool> electrostaticSitesI = constPropI.electrostaticSites();
-
-    List<bool> pairPotentialSitesJ = constPropJ.pairPotentialSites();
-
-    List<bool> electrostaticSitesJ = constPropJ.electrostaticSites();
-
-    forAll(siteIdsI, sI)
     {
-        label idsI(siteIdsI[sI]);
+        label sI = constPropI.pairPotSites()[pI];
 
-        forAll(siteIdsJ, sJ)
+        label idsI = constPropI.sites()[sI].siteId();
+
+        forAll(constPropJ.pairPotSites(), pJ)
         {
-            label idsJ(siteIdsJ[sJ]);
+            label sJ = constPropJ.pairPotSites()[pJ];
 
-            if (pairPotentialSitesI[sI] && pairPotentialSitesJ[sJ])
+            label idsJ = constPropJ.sites()[sJ].siteId();
+
+            vector rsIsJ =
+                molI.sitePositions()[sI] - molJ.sitePositions()[sJ];
+
+            scalar rsIsJMagSq = magSqr(rsIsJ);
+
+            if (pairPot.rCutSqr(idsI, idsJ, rsIsJMagSq))
             {
-                vector rsIsJ =
+                scalar rsIsJMag = mag(rsIsJ);
-                    molI.sitePositions()[sI] - molJ.sitePositions()[sJ];
 
-                scalar rsIsJMagSq = magSqr(rsIsJ);
+                // Guard against pairPot.energy being evaluated
+                // if rIJMag < SMALL. A floating point exception will
+                // happen otherwise.
 
-                if (pairPot.rCutSqr(idsI, idsJ, rsIsJMagSq))
+                if (rsIsJMag < SMALL)
                 {
-                    scalar rsIsJMag = mag(rsIsJ);
+                    WarningIn("polyatomicCloud::removeHighEnergyOverlaps()")
+                        << "Molecule site pair closer than "
+                        << SMALL
+                        << ": mag separation = " << rsIsJMag
+                        << ". These may have been placed on top of each"
+                        << " other by a rounding error in mdInitialise in"
+                        << " parallel or a block filled with polyatomics"
+                        << " twice. Removing one of the polyatomics."
+                        << endl;
 
-                    // Guard against pairPot.energy being evaluated
+                    return true;
-                    // if rIJMag < SMALL. A floating point exception will
-                    // happen otherwise.
-
-                    if (rsIsJMag < SMALL)
-                    {
-                        WarningIn("polyatomicCloud::removeHighEnergyOverlaps()")
-                            << "Molecule site pair closer than "
-                            << SMALL
-                            << ": mag separation = " << rsIsJMag
-                            << ". These may have been placed on top of each"
-                            << " other by a rounding error in mdInitialise in"
-                            << " parallel or a block filled with polyatomics"
-                            << " twice. Removing one of the polyatomics."
-                            << endl;
-
-                        return true;
-                    }
-
-                    // Guard against pairPot.energy being evaluated if rIJMag <
-                    // rMin.  A tabulation lookup error will occur otherwise.
-
-                    if (rsIsJMag < pairPot.rMin(idsI, idsJ))
-                    {
-                        return true;
-                    }
-
-                    if
-                    (
-                        mag(pairPot.energy(idsI, idsJ, rsIsJMag))
-                      > pot_.potentialEnergyLimit()
-                    )
-                    {
-                        return true;
-                    };
                 }
+
+                // Guard against pairPot.energy being evaluated if rIJMag <
+                // rMin.  A tabulation lookup error will occur otherwise.
+
+                if (rsIsJMag < pairPot.rMin(idsI, idsJ))
+                {
+                    return true;
+                }
+
+                if
+                (
+                    mag(pairPot.energy(idsI, idsJ, rsIsJMag))
+                  > pot_.potentialEnergyLimit()
+                )
+                {
+                    return true;
+                };
             }
+        }
+    }
 
-            if (electrostaticSitesI[sI] && electrostaticSitesJ[sJ])
+    forAll(constPropI.electrostaticSites(), pI)
+    {
+        label sI = constPropI.electrostaticSites()[pI];
+
+        forAll(constPropJ.electrostaticSites(), pJ)
+        {
+            label sJ = constPropJ.electrostaticSites()[pJ];
+
+            vector rsIsJ =
+                molI.sitePositions()[sI] - molJ.sitePositions()[sJ];
+
+            scalar rsIsJMagSq = magSqr(rsIsJ);
+
+            if (pairPot.rCutMaxSqr(rsIsJMagSq))
             {
-                vector rsIsJ =
+                scalar rsIsJMag = mag(rsIsJ);
-                    molI.sitePositions()[sI] - molJ.sitePositions()[sJ];
 
-                scalar rsIsJMagSq = magSqr(rsIsJ);
+                // Guard against pairPot.energy being evaluated
+                // if rIJMag < SMALL. A floating point exception will
+                // happen otherwise.
 
-                if (pairPot.rCutMaxSqr(rsIsJMagSq))
+                if (rsIsJMag < SMALL)
                 {
-                    scalar rsIsJMag = mag(rsIsJ);
+                    WarningIn("polyatomicCloud::removeHighEnergyOverlaps()")
+                        << "Molecule site pair closer than "
+                        << SMALL
+                        << ": mag separation = " << rsIsJMag
+                        << ". These may have been placed on top of each"
+                        << " other by a rounding error in mdInitialise in"
+                        << " parallel or a block filled with polyatomics"
+                        << " twice. Removing one of the polyatomics."
+                        << endl;
 
-                    // Guard against pairPot.energy being evaluated
+                    return true;
-                    // if rIJMag < SMALL. A floating point exception will
-                    // happen otherwise.
-
-                    if (rsIsJMag < SMALL)
-                    {
-                        WarningIn("polyatomicCloud::removeHighEnergyOverlaps()")
-                            << "Molecule site pair closer than "
-                            << SMALL
-                            << ": mag separation = " << rsIsJMag
-                            << ". These may have been placed on top of each"
-                            << " other by a rounding error in mdInitialise in"
-                            << " parallel or a block filled with polyatomics"
-                            << " twice. Removing one of the polyatomics."
-                            << endl;
-
-                        return true;
-                    }
-
-                    if (rsIsJMag < electrostatic.rMin())
-                    {
-                        return true;
-                    }
-
-                    scalar chargeI = constPropI.siteCharges()[sI];
-
-                    scalar chargeJ = constPropJ.siteCharges()[sJ];
-
-                    if
-                    (
-                        mag(chargeI*chargeJ*electrostatic.energy(rsIsJMag))
-                      > pot_.potentialEnergyLimit()
-                    )
-                    {
-                        return true;
-                    };
                 }
+
+                if (rsIsJMag < electrostatic.rMin())
+                {
+                    return true;
+                }
+
+                scalar chargeI = constPropI.sites()[sI].siteCharge();
+
+                scalar chargeJ = constPropJ.sites()[sJ].siteCharge();
+
+                if
+                (
+                    mag(chargeI*chargeJ*electrostatic.energy(rsIsJMag))
+                  > pot_.potentialEnergyLimit()
+                )
+                {
+                    return true;
+                };
             }
         }
     }