Electrostatics now implemented using pair potential object, stored in pairPotentialList. Using damped, shifted force coulomb interactions as per J.Chem.Phys. 124, 234104.

src/lagrangian/molecularDynamics/molecule/mdTools/temperatureAndPressure.H

@@ -84,12 +84,9 @@ if (runTime.outputTime())
         Info << "----------------------------------------" << nl
             << "Averaged properties" << nl
             << "Average |velocity| = "
-                << mag(accumulatedTotalLinearMomentum)/accumulatedTotalMass
-                << " m/s" << nl
-            << "Average temperature = "
-                << averageTemperature << " K" << nl
-            << "Average pressure = "
-                << averagePressure << " N/m^2" << nl
+            << mag(accumulatedTotalLinearMomentum)/accumulatedTotalMass << nl
+            << "Average temperature = " << averageTemperature << nl
+            << "Average pressure = " << averagePressure << nl
             << "----------------------------------------" << endl;
     }
     else

src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudI.H

@@ -34,7 +34,7 @@ inline void Foam::moleculeCloud::evaluatePair
 {
     const pairPotentialList& pairPot(pot_.pairPotentials());
 
-    const electrostaticPotential& electrostatic(pot_.electrostatic());
+    const pairPotential& electrostatic(pairPot.electrostatic());
 
     label idI = molI->id();
 
@@ -143,7 +143,7 @@ inline void Foam::moleculeCloud::evaluatePair
 {
     const pairPotentialList& pairPot(pot_.pairPotentials());
 
-    const electrostaticPotential& electrostatic(pot_.electrostatic());
+    const pairPotential& electrostatic(pairPot.electrostatic());
 
     label idReal = molReal->id();
 
@@ -243,7 +243,7 @@ inline bool Foam::moleculeCloud::evaluatePotentialLimit
 {
     const pairPotentialList& pairPot(pot_.pairPotentials());
 
-    const electrostaticPotential& electrostatic(pot_.electrostatic());
+    const pairPotential& electrostatic(pairPot.electrostatic());
 
     label idI = molI->id();
 
@@ -381,7 +381,7 @@ inline bool Foam::moleculeCloud::evaluatePotentialLimit
 {
     const pairPotentialList& pairPot(pot_.pairPotentials());
 
-    const electrostaticPotential& electrostatic(pot_.electrostatic());
+    const pairPotential& electrostatic(pairPot.electrostatic());
 
     label idReal = molReal->id();
 

src/lagrangian/molecularDynamics/potential/Make/files

@@ -15,6 +15,7 @@ $(pairPotential)/derived/maitlandSmith/maitlandSmith.C
 $(pairPotential)/derived/azizChen/azizChen.C
 $(pairPotential)/derived/exponentialRepulsion/exponentialRepulsion.C
 $(pairPotential)/derived/coulomb/coulomb.C
+$(pairPotential)/derived/dampedCoulomb/dampedCoulomb.C
 $(pairPotential)/derived/noInteraction/noInteraction.C
 
 energyScalingFunction = energyScalingFunction
@@ -43,4 +44,4 @@ electrostaticPotential = electrostaticPotential
 
 $(electrostaticPotential)/electrostaticPotential.C
 
-LIB = $(FOAM_LIBBIN)/libpotential
+LIB = $(FOAM_LIBBIN)/libpotential

src/lagrangian/molecularDynamics/potential/energyScalingFunction/derived/doubleSigmoid/doubleSigmoid.H

@@ -98,7 +98,7 @@ public:
 
         void scaleEnergy(scalar& e, const scalar r) const;
 
-        //- Read transportProperties dictionary
+        //- Read dictionary
         bool read(const dictionary& energyScalingFunctionProperties);
 };
 

src/lagrangian/molecularDynamics/potential/energyScalingFunction/derived/noScaling/noScaling.H

@@ -80,7 +80,7 @@ public:
 
         void scaleEnergy(scalar& e, const scalar r) const;
 
-        //- Read transportProperties dictionary
+        //- Read dictionary
         bool read(const dictionary& energyScalingFunctionProperties);
 };
 

src/lagrangian/molecularDynamics/potential/energyScalingFunction/derived/shifted/shifted.H

@@ -83,7 +83,7 @@ public:
 
         void scaleEnergy(scalar& e, const scalar r) const;
 
-        //- Read transportProperties dictionary
+        //- Read dictionary
         bool read(const dictionary& energyScalingFunctionProperties);
 };
 

src/lagrangian/molecularDynamics/potential/energyScalingFunction/derived/shiftedForce/shiftedForce.H

@@ -89,7 +89,7 @@ public:
 
         void scaleEnergy(scalar& e, const scalar r) const;
 
-        //- Read transportProperties dictionary
+        //- Read dictionary
         bool read(const dictionary& energyScalingFunctionProperties);
 };
 

src/lagrangian/molecularDynamics/potential/energyScalingFunction/derived/sigmoid/sigmoid.H

@@ -98,7 +98,7 @@ public:
 
         void scaleEnergy(scalar& e, const scalar r) const;
 
-        //- Read transportProperties dictionary
+        //- Read dictionary
         bool read(const dictionary& energyScalingFunctionProperties);
 };
 

src/lagrangian/molecularDynamics/potential/pairPotential/basic/pairPotential.C

@@ -93,7 +93,7 @@ void Foam::pairPotential::setLookupTables()
 }
 
 
-Foam::scalar Foam::pairPotential::forceLookup(const scalar r) const
+Foam::scalar Foam::pairPotential::force(const scalar r) const
 {
     scalar k_rIJ = (r - rMin_)/dr_;
 
@@ -102,7 +102,7 @@ Foam::scalar Foam::pairPotential::forceLookup(const scalar r) const
     if (k < 0)
     {
         FatalErrorIn("pairPotential.C") << nl
-            << "r less than rMin" << nl
+            << "r less than rMin in pair potential " << name_ << nl
             << abort(FatalError);
     }
 
@@ -130,7 +130,7 @@ Foam::pairPotential::forceTable() const
 }
 
 
-Foam::scalar Foam::pairPotential::energyLookup(const scalar r) const
+Foam::scalar Foam::pairPotential::energy(const scalar r) const
 {
     scalar k_rIJ = (r - rMin_)/dr_;
 
@@ -139,7 +139,7 @@ Foam::scalar Foam::pairPotential::energyLookup(const scalar r) const
     if (k < 0)
     {
         FatalErrorIn("pairPotential.C") << nl
-            << "r less than rMin" << nl
+            << "r less than rMin in pair potential " << name_ << nl
             << abort(FatalError);
     }
 

src/lagrangian/molecularDynamics/potential/pairPotential/basic/pairPotential.H

@@ -94,7 +94,7 @@ protected:
 public:
 
     //- Runtime type information
-    TypeName("pairPotential");
+        TypeName("pairPotential");
 
 
     // Declare run-time constructor selection table
@@ -150,9 +150,9 @@ public:
 
         inline scalar rCutSqr() const;
 
-        scalar energyLookup (const scalar r) const;
+        scalar energy (const scalar r) const;
 
-        scalar forceLookup (const scalar r) const;
+        scalar force (const scalar r) const;
 
         List<Pair<scalar> > energyTable() const;
 

src/lagrangian/molecularDynamics/potential/pairPotential/derived/azizChen/azizChen.H

@@ -111,7 +111,7 @@ public:
 
         scalar unscaledEnergy(const scalar r) const;
 
-        //- Read transportProperties dictionary
+        //- Read dictionary
         bool read(const dictionary& pairPotentialProperties);
 };
 

src/lagrangian/molecularDynamics/potential/pairPotential/derived/coulomb/coulomb.C

@@ -25,6 +25,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "coulomb.H"
+#include "mathematicalConstants.H"
 #include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -45,6 +46,9 @@ addToRunTimeSelectionTable
     dictionary
 );
 
+scalar coulomb::oneOverFourPiEps0 =
+1.0/(4.0 * mathematicalConstant::pi * 8.854187817e-12);
+
 // * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
 
 
@@ -65,7 +69,7 @@ coulomb::coulomb
 
 scalar coulomb::unscaledEnergy(const scalar r) const
 {
-    return 1.0/(4.0 * mathematicalConstant::pi * 8.854187817e-12 * r);
+    return oneOverFourPiEps0/r;
 }
 
 

src/lagrangian/molecularDynamics/potential/pairPotential/derived/coulomb/coulomb.H

@@ -23,7 +23,7 @@ License
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
 Class
-    Foam::pairPotentials::coulomb
+    Foam::pairPotentials::electrostatic
 
 Description
 
@@ -53,12 +53,14 @@ class coulomb
 :
     public pairPotential
 {
-
 public:
 
     //- Runtime type information
-    TypeName("coulomb");
+        TypeName("coulomb");
 
+    // Static data members
+
+        static scalar oneOverFourPiEps0;
 
     // Constructors
 
@@ -80,7 +82,7 @@ public:
 
         scalar unscaledEnergy(const scalar r) const;
 
-        //- Read transportProperties dictionary
+        //- Read dictionary
         bool read(const dictionary& pairPotentialProperties);
 };
 

src/lagrangian/molecularDynamics/potential/pairPotential/derived/dampedCoulomb/dampedCoulomb.C

@@ -0,0 +1,99 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "dampedCoulomb.H"
+#include "mathematicalConstants.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace pairPotentials
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(dampedCoulomb, 0);
+
+addToRunTimeSelectionTable
+(
+    pairPotential,
+    dampedCoulomb,
+    dictionary
+);
+
+scalar dampedCoulomb::oneOverFourPiEps0 =
+1.0/(4.0 * mathematicalConstant::pi * 8.854187817e-12);
+
+// * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+dampedCoulomb::dampedCoulomb
+(
+    const word& name,
+    const dictionary& pairPotentialProperties
+)
+:
+    pairPotential(name, pairPotentialProperties),
+    dampedCoulombCoeffs_
+    (
+        pairPotentialProperties.subDict(typeName + "Coeffs")
+    ),
+    alpha_(readScalar(dampedCoulombCoeffs_.lookup("alpha")))
+{
+    setLookupTables();
+}
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+scalar dampedCoulomb::unscaledEnergy(const scalar r) const
+{
+    return oneOverFourPiEps0*erfc(alpha_*r)/r;
+}
+
+
+bool dampedCoulomb::read(const dictionary& pairPotentialProperties)
+{
+    pairPotential::read(pairPotentialProperties);
+
+    dampedCoulombCoeffs_ =
+    pairPotentialProperties.subDict(typeName + "Coeffs");
+
+    dampedCoulombCoeffs_.lookup("alpha") >> alpha_;
+
+    return true;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace pairPotentials
+} // End namespace Foam
+
+// ************************************************************************* //

src/lagrangian/molecularDynamics/potential/pairPotential/derived/dampedCoulomb/dampedCoulomb.H

@@ -0,0 +1,105 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::pairPotentials::electrostatic
+
+Description
+
+
+SourceFiles
+    dampedCoulomb.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef dampedCoulomb_H
+#define dampedCoulomb_H
+
+#include "pairPotential.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace pairPotentials
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class dampedCoulomb Declaration
+\*---------------------------------------------------------------------------*/
+
+class dampedCoulomb
+:
+    public pairPotential
+{
+    // Private data
+
+        dictionary dampedCoulombCoeffs_;
+
+        scalar alpha_;
+
+public:
+
+    //- Runtime type information
+        TypeName("dampedCoulomb");
+
+    // Static data members
+
+        static scalar oneOverFourPiEps0;
+
+    // Constructors
+
+        //- Construct from components
+        dampedCoulomb
+        (
+            const word& name,
+            const dictionary& pairPotentialProperties
+        );
+
+
+    // Destructor
+
+        ~dampedCoulomb()
+        {}
+
+
+    // Member Functions
+
+        scalar unscaledEnergy(const scalar r) const;
+
+        //- Read dictionary
+        bool read(const dictionary& pairPotentialProperties);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace pairPotentials
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/lagrangian/molecularDynamics/potential/pairPotential/derived/exponentialRepulsion/exponentialRepulsion.H

@@ -86,7 +86,7 @@ public:
 
         scalar unscaledEnergy(const scalar r) const;
 
-        //- Read transportProperties dictionary
+        //- Read dictionary
         bool read(const dictionary& pairPotentialProperties);
 };
 

src/lagrangian/molecularDynamics/potential/pairPotential/derived/lennardJones/lennardJones.H

@@ -86,7 +86,7 @@ public:
 
         scalar unscaledEnergy(const scalar r) const;
 
-        //- Read transportProperties dictionary
+        //- Read dictionary
         bool read(const dictionary& pairPotentialProperties);
 };
 

src/lagrangian/molecularDynamics/potential/pairPotential/derived/maitlandSmith/maitlandSmith.H

@@ -114,7 +114,7 @@ public:
 
         scalar unscaledEnergy(const scalar r) const;
 
-        //- Read transportProperties dictionary
+        //- Read dictionary
         bool read(const dictionary& pairPotentialProperties);
 };
 

src/lagrangian/molecularDynamics/potential/pairPotential/derived/noInteraction/noInteraction.H

@@ -80,7 +80,7 @@ public:
 
         scalar unscaledEnergy(const scalar r) const;
 
-        //- Read transportProperties dictionary
+        //- Read dictionary
         bool read(const dictionary& pairPotentialProperties);
 };
 

src/lagrangian/molecularDynamics/potential/pairPotential/pairPotentialList/pairPotentialList.C

@@ -136,6 +136,37 @@ void Foam::pairPotentialList::readPairPotentialDict
         }
     }
 
+    if (!pairPotentialDict.found("electrostatic"))
+    {
+        FatalErrorIn("pairPotentialList::buildPotentials") << nl
+            << "Pair pairPotential specification subDict electrostatic"
+            << nl << abort(FatalError);
+    }
+
+    electrostaticPotential_ = pairPotential::New
+    (
+        "electrostatic",
+        pairPotentialDict.subDict("electrostatic")
+    );
+
+    if (electrostaticPotential_->rCut() > rCutMax_)
+    {
+        rCutMax_ = electrostaticPotential_->rCut();
+    }
+
+    if (electrostaticPotential_->writeTables())
+    {
+        OFstream ppTabFile(mesh.time().path()/"electrostatic");
+
+        if(!electrostaticPotential_->writeEnergyAndForceTables(ppTabFile))
+        {
+            FatalErrorIn("pairPotentialList::readPairPotentialDict")
+                << "Failed writing to "
+                << ppTabFile.name() << nl
+                << abort(FatalError);
+        }
+    }
+
     rCutMaxSqr_ = rCutMax_*rCutMax_;
 }
 
@@ -270,7 +301,7 @@ Foam::scalar Foam::pairPotentialList::force
     const scalar rIJMag
 ) const
 {
-    scalar f = (*this)[pairPotentialIndex (a, b)].forceLookup(rIJMag);
+    scalar f = (*this)[pairPotentialIndex (a, b)].force(rIJMag);
 
     return f;
 }
@@ -283,7 +314,7 @@ Foam::scalar Foam::pairPotentialList::energy
     const scalar rIJMag
 ) const
 {
-    scalar e = (*this)[pairPotentialIndex (a, b)].energyLookup(rIJMag);
+    scalar e = (*this)[pairPotentialIndex (a, b)].energy(rIJMag);
 
     return e;
 }

src/lagrangian/molecularDynamics/potential/pairPotential/pairPotentialList/pairPotentialList.H

@@ -62,6 +62,8 @@ class pairPotentialList
 
         scalar rCutMaxSqr_;
 
+        autoPtr<pairPotential> electrostaticPotential_;
+
     // Private Member Functions
 
         inline label pairPotentialIndex
@@ -153,6 +155,8 @@ public:
                 const label b,
                 const scalar rIJMag
             ) const;
+
+            inline const pairPotential& electrostatic() const;
 };
 
 

src/lagrangian/molecularDynamics/potential/pairPotential/pairPotentialList/pairPotentialListI.H

@@ -70,4 +70,10 @@ inline Foam::scalar Foam::pairPotentialList::rCutMaxSqr() const
 }
 
 
+inline const Foam::pairPotential& Foam::pairPotentialList::electrostatic() const
+{
+    return electrostaticPotential_;
+}
+
+
 // ************************************************************************* //

src/lagrangian/molecularDynamics/potential/potential/potential.C

@@ -365,8 +365,7 @@ void Foam::potential::potential::readMdInitialiseDict
 
 Foam::potential::potential(const polyMesh& mesh)
 :
-    mesh_(mesh),
-    electrostaticPotential_()
+    mesh_(mesh)
 {
     readPotentialDict();
 }
@@ -379,8 +378,7 @@ Foam::potential::potential
     IOdictionary& idListDict
 )
 :
-    mesh_(mesh),
-    electrostaticPotential_()
+    mesh_(mesh)
 {
     readMdInitialiseDict(mdInitialiseDict, idListDict);
 }

src/lagrangian/molecularDynamics/potential/potential/potential.H

@@ -70,8 +70,6 @@ class potential
 
         pairPotentialList pairPotentials_;
 
-        electrostaticPotential electrostaticPotential_;
-
         tetherPotentialList tetherPotentials_;
 
         vector gravity_;
@@ -133,8 +131,6 @@ public:
 
             inline const pairPotentialList& pairPotentials() const;
 
-            inline const electrostaticPotential& electrostatic() const;
-
             inline const tetherPotentialList& tetherPotentials() const;
 
             inline const vector& gravity() const;

src/lagrangian/molecularDynamics/potential/potential/potentialI.H

@@ -70,13 +70,6 @@ inline const Foam::pairPotentialList& Foam::potential::pairPotentials() const
 }
 
 
-inline const Foam::electrostaticPotential&
-Foam::potential::electrostatic() const
-{
-    return electrostaticPotential_;
-}
-
-
 inline const Foam::tetherPotentialList&
 Foam::potential::tetherPotentials() const
 {

src/lagrangian/molecularDynamics/potential/tetherPotential/derived/harmonicSpring/harmonicSpring.H

@@ -88,7 +88,7 @@ public:
 
         vector force(const vector r) const;
 
-        //- Read transportProperties dictionary
+        //- Read dictionary
         bool read(const dictionary& tetherPotentialProperties);
 };
 

src/lagrangian/molecularDynamics/potential/tetherPotential/derived/pitchForkRing/pitchForkRing.H

@@ -90,7 +90,7 @@ public:
 
         vector force(const vector r) const;
 
-        //- Read transportProperties dictionary
+        //- Read dictionary
         bool read(const dictionary& tetherPotentialProperties);
 };
 

src/lagrangian/molecularDynamics/potential/tetherPotential/derived/restrainedHarmonicSpring/restrainedHarmonicSpring.H

@@ -90,7 +90,7 @@ public:
 
         vector force(const vector r) const;
 
-        //- Read transportProperties dictionary
+        //- Read dictionary
         bool read(const dictionary& tetherPotentialProperties);
 };