From 30c9efca7d63b657ccf5e49b1076f25ddc1a6df8 Mon Sep 17 00:00:00 2001 From: graham Date: Thu, 15 Jan 2009 20:00:51 +0000 Subject: [PATCH 01/21] Electrostatics now implemented using pair potential object, stored in pairPotentialList. Using damped, shifted force coulomb interactions as per J.Chem.Phys. 124, 234104. --- .../molecule/mdTools/temperatureAndPressure.H | 9 +- .../molecule/moleculeCloud/moleculeCloudI.H | 8 +- .../molecularDynamics/potential/Make/files | 3 +- .../derived/doubleSigmoid/doubleSigmoid.H | 2 +- .../derived/noScaling/noScaling.H | 2 +- .../derived/shifted/shifted.H | 2 +- .../derived/shiftedForce/shiftedForce.H | 2 +- .../derived/sigmoid/sigmoid.H | 2 +- .../pairPotential/basic/pairPotential.C | 8 +- .../pairPotential/basic/pairPotential.H | 6 +- .../pairPotential/derived/azizChen/azizChen.H | 2 +- .../pairPotential/derived/coulomb/coulomb.C | 6 +- .../pairPotential/derived/coulomb/coulomb.H | 10 +- .../derived/dampedCoulomb/dampedCoulomb.C | 99 +++++++++++++++++ .../derived/dampedCoulomb/dampedCoulomb.H | 105 ++++++++++++++++++ .../exponentialRepulsion.H | 2 +- .../derived/lennardJones/lennardJones.H | 2 +- .../derived/maitlandSmith/maitlandSmith.H | 2 +- .../derived/noInteraction/noInteraction.H | 2 +- .../pairPotentialList/pairPotentialList.C | 35 +++++- .../pairPotentialList/pairPotentialList.H | 4 + .../pairPotentialList/pairPotentialListI.H | 6 + .../potential/potential/potential.C | 6 +- .../potential/potential/potential.H | 4 - .../potential/potential/potentialI.H | 7 -- .../derived/harmonicSpring/harmonicSpring.H | 2 +- .../derived/pitchForkRing/pitchForkRing.H | 2 +- .../restrainedHarmonicSpring.H | 2 +- 28 files changed, 289 insertions(+), 53 deletions(-) create mode 100644 src/lagrangian/molecularDynamics/potential/pairPotential/derived/dampedCoulomb/dampedCoulomb.C create mode 100644 src/lagrangian/molecularDynamics/potential/pairPotential/derived/dampedCoulomb/dampedCoulomb.H diff --git a/src/lagrangian/molecularDynamics/molecule/mdTools/temperatureAndPressure.H b/src/lagrangian/molecularDynamics/molecule/mdTools/temperatureAndPressure.H index 4ba2bc1bc4..250fa1e2ee 100644 --- a/src/lagrangian/molecularDynamics/molecule/mdTools/temperatureAndPressure.H +++ b/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 diff --git a/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudI.H b/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudI.H index 55813e6658..9dd9ee06ae 100644 --- a/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudI.H +++ b/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(); diff --git a/src/lagrangian/molecularDynamics/potential/Make/files b/src/lagrangian/molecularDynamics/potential/Make/files index 7b16b6e2b6..4afe4de3d9 100644 --- a/src/lagrangian/molecularDynamics/potential/Make/files +++ b/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 \ No newline at end of file diff --git a/src/lagrangian/molecularDynamics/potential/energyScalingFunction/derived/doubleSigmoid/doubleSigmoid.H b/src/lagrangian/molecularDynamics/potential/energyScalingFunction/derived/doubleSigmoid/doubleSigmoid.H index 07c4dfbd49..6144bc388d 100644 --- a/src/lagrangian/molecularDynamics/potential/energyScalingFunction/derived/doubleSigmoid/doubleSigmoid.H +++ b/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); }; diff --git a/src/lagrangian/molecularDynamics/potential/energyScalingFunction/derived/noScaling/noScaling.H b/src/lagrangian/molecularDynamics/potential/energyScalingFunction/derived/noScaling/noScaling.H index 192c862fb4..845634d7be 100644 --- a/src/lagrangian/molecularDynamics/potential/energyScalingFunction/derived/noScaling/noScaling.H +++ b/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); }; diff --git a/src/lagrangian/molecularDynamics/potential/energyScalingFunction/derived/shifted/shifted.H b/src/lagrangian/molecularDynamics/potential/energyScalingFunction/derived/shifted/shifted.H index 6f95adbdb9..1aeca9e87f 100644 --- a/src/lagrangian/molecularDynamics/potential/energyScalingFunction/derived/shifted/shifted.H +++ b/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); }; diff --git a/src/lagrangian/molecularDynamics/potential/energyScalingFunction/derived/shiftedForce/shiftedForce.H b/src/lagrangian/molecularDynamics/potential/energyScalingFunction/derived/shiftedForce/shiftedForce.H index d67164c64e..5f90885481 100644 --- a/src/lagrangian/molecularDynamics/potential/energyScalingFunction/derived/shiftedForce/shiftedForce.H +++ b/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); }; diff --git a/src/lagrangian/molecularDynamics/potential/energyScalingFunction/derived/sigmoid/sigmoid.H b/src/lagrangian/molecularDynamics/potential/energyScalingFunction/derived/sigmoid/sigmoid.H index 42d0ac85b8..81a6e31261 100644 --- a/src/lagrangian/molecularDynamics/potential/energyScalingFunction/derived/sigmoid/sigmoid.H +++ b/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); }; diff --git a/src/lagrangian/molecularDynamics/potential/pairPotential/basic/pairPotential.C b/src/lagrangian/molecularDynamics/potential/pairPotential/basic/pairPotential.C index d4a807f352..0d42b0b2f9 100644 --- a/src/lagrangian/molecularDynamics/potential/pairPotential/basic/pairPotential.C +++ b/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); } diff --git a/src/lagrangian/molecularDynamics/potential/pairPotential/basic/pairPotential.H b/src/lagrangian/molecularDynamics/potential/pairPotential/basic/pairPotential.H index 6d6701a106..c804030307 100644 --- a/src/lagrangian/molecularDynamics/potential/pairPotential/basic/pairPotential.H +++ b/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 > energyTable() const; diff --git a/src/lagrangian/molecularDynamics/potential/pairPotential/derived/azizChen/azizChen.H b/src/lagrangian/molecularDynamics/potential/pairPotential/derived/azizChen/azizChen.H index 43a1ee9dca..e7904674dd 100644 --- a/src/lagrangian/molecularDynamics/potential/pairPotential/derived/azizChen/azizChen.H +++ b/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); }; diff --git a/src/lagrangian/molecularDynamics/potential/pairPotential/derived/coulomb/coulomb.C b/src/lagrangian/molecularDynamics/potential/pairPotential/derived/coulomb/coulomb.C index 0476e378da..d7a87238f0 100644 --- a/src/lagrangian/molecularDynamics/potential/pairPotential/derived/coulomb/coulomb.C +++ b/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; } diff --git a/src/lagrangian/molecularDynamics/potential/pairPotential/derived/coulomb/coulomb.H b/src/lagrangian/molecularDynamics/potential/pairPotential/derived/coulomb/coulomb.H index 44e5e2cc5c..115884990b 100644 --- a/src/lagrangian/molecularDynamics/potential/pairPotential/derived/coulomb/coulomb.H +++ b/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); }; diff --git a/src/lagrangian/molecularDynamics/potential/pairPotential/derived/dampedCoulomb/dampedCoulomb.C b/src/lagrangian/molecularDynamics/potential/pairPotential/derived/dampedCoulomb/dampedCoulomb.C new file mode 100644 index 0000000000..ca1038350c --- /dev/null +++ b/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 + +// ************************************************************************* // diff --git a/src/lagrangian/molecularDynamics/potential/pairPotential/derived/dampedCoulomb/dampedCoulomb.H b/src/lagrangian/molecularDynamics/potential/pairPotential/derived/dampedCoulomb/dampedCoulomb.H new file mode 100644 index 0000000000..625d9a9081 --- /dev/null +++ b/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 + +// ************************************************************************* // diff --git a/src/lagrangian/molecularDynamics/potential/pairPotential/derived/exponentialRepulsion/exponentialRepulsion.H b/src/lagrangian/molecularDynamics/potential/pairPotential/derived/exponentialRepulsion/exponentialRepulsion.H index 7026710427..895bb05706 100644 --- a/src/lagrangian/molecularDynamics/potential/pairPotential/derived/exponentialRepulsion/exponentialRepulsion.H +++ b/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); }; diff --git a/src/lagrangian/molecularDynamics/potential/pairPotential/derived/lennardJones/lennardJones.H b/src/lagrangian/molecularDynamics/potential/pairPotential/derived/lennardJones/lennardJones.H index f12027d115..829827120b 100644 --- a/src/lagrangian/molecularDynamics/potential/pairPotential/derived/lennardJones/lennardJones.H +++ b/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); }; diff --git a/src/lagrangian/molecularDynamics/potential/pairPotential/derived/maitlandSmith/maitlandSmith.H b/src/lagrangian/molecularDynamics/potential/pairPotential/derived/maitlandSmith/maitlandSmith.H index af4612644b..d2097b42fe 100644 --- a/src/lagrangian/molecularDynamics/potential/pairPotential/derived/maitlandSmith/maitlandSmith.H +++ b/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); }; diff --git a/src/lagrangian/molecularDynamics/potential/pairPotential/derived/noInteraction/noInteraction.H b/src/lagrangian/molecularDynamics/potential/pairPotential/derived/noInteraction/noInteraction.H index d09564ef88..d4d7a9f208 100644 --- a/src/lagrangian/molecularDynamics/potential/pairPotential/derived/noInteraction/noInteraction.H +++ b/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); }; diff --git a/src/lagrangian/molecularDynamics/potential/pairPotential/pairPotentialList/pairPotentialList.C b/src/lagrangian/molecularDynamics/potential/pairPotential/pairPotentialList/pairPotentialList.C index ef856bc505..c9ff27fe1d 100644 --- a/src/lagrangian/molecularDynamics/potential/pairPotential/pairPotentialList/pairPotentialList.C +++ b/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; } diff --git a/src/lagrangian/molecularDynamics/potential/pairPotential/pairPotentialList/pairPotentialList.H b/src/lagrangian/molecularDynamics/potential/pairPotential/pairPotentialList/pairPotentialList.H index a7e8190a01..a70179fc78 100644 --- a/src/lagrangian/molecularDynamics/potential/pairPotential/pairPotentialList/pairPotentialList.H +++ b/src/lagrangian/molecularDynamics/potential/pairPotential/pairPotentialList/pairPotentialList.H @@ -62,6 +62,8 @@ class pairPotentialList scalar rCutMaxSqr_; + autoPtr electrostaticPotential_; + // Private Member Functions inline label pairPotentialIndex @@ -153,6 +155,8 @@ public: const label b, const scalar rIJMag ) const; + + inline const pairPotential& electrostatic() const; }; diff --git a/src/lagrangian/molecularDynamics/potential/pairPotential/pairPotentialList/pairPotentialListI.H b/src/lagrangian/molecularDynamics/potential/pairPotential/pairPotentialList/pairPotentialListI.H index 413ffea05e..844f641ca1 100644 --- a/src/lagrangian/molecularDynamics/potential/pairPotential/pairPotentialList/pairPotentialListI.H +++ b/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_; +} + + // ************************************************************************* // diff --git a/src/lagrangian/molecularDynamics/potential/potential/potential.C b/src/lagrangian/molecularDynamics/potential/potential/potential.C index 97414fd7c1..84eb338ddd 100644 --- a/src/lagrangian/molecularDynamics/potential/potential/potential.C +++ b/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); } diff --git a/src/lagrangian/molecularDynamics/potential/potential/potential.H b/src/lagrangian/molecularDynamics/potential/potential/potential.H index 6f6d47d9ee..890d179238 100644 --- a/src/lagrangian/molecularDynamics/potential/potential/potential.H +++ b/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; diff --git a/src/lagrangian/molecularDynamics/potential/potential/potentialI.H b/src/lagrangian/molecularDynamics/potential/potential/potentialI.H index eddd2b047c..f0cd7de3da 100644 --- a/src/lagrangian/molecularDynamics/potential/potential/potentialI.H +++ b/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 { diff --git a/src/lagrangian/molecularDynamics/potential/tetherPotential/derived/harmonicSpring/harmonicSpring.H b/src/lagrangian/molecularDynamics/potential/tetherPotential/derived/harmonicSpring/harmonicSpring.H index 6508adaec9..5ecb91023c 100644 --- a/src/lagrangian/molecularDynamics/potential/tetherPotential/derived/harmonicSpring/harmonicSpring.H +++ b/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); }; diff --git a/src/lagrangian/molecularDynamics/potential/tetherPotential/derived/pitchForkRing/pitchForkRing.H b/src/lagrangian/molecularDynamics/potential/tetherPotential/derived/pitchForkRing/pitchForkRing.H index 82b93b087b..f1c734a14f 100644 --- a/src/lagrangian/molecularDynamics/potential/tetherPotential/derived/pitchForkRing/pitchForkRing.H +++ b/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); }; diff --git a/src/lagrangian/molecularDynamics/potential/tetherPotential/derived/restrainedHarmonicSpring/restrainedHarmonicSpring.H b/src/lagrangian/molecularDynamics/potential/tetherPotential/derived/restrainedHarmonicSpring/restrainedHarmonicSpring.H index 62b9261859..5680437bd3 100644 --- a/src/lagrangian/molecularDynamics/potential/tetherPotential/derived/restrainedHarmonicSpring/restrainedHarmonicSpring.H +++ b/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); }; From 3ec52af20f4b251a8fa35623e7f0a6b33088ff80 Mon Sep 17 00:00:00 2001 From: graham Date: Thu, 15 Jan 2009 20:10:53 +0000 Subject: [PATCH 02/21] Adding rMin check to electrostatic evaluation in removeHighEnergyOverlaps --- .../molecule/moleculeCloud/moleculeCloudI.H | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudI.H b/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudI.H index 9dd9ee06ae..e5320a4f3b 100644 --- a/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudI.H +++ b/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudI.H @@ -352,6 +352,11 @@ inline bool Foam::moleculeCloud::evaluatePotentialLimit return true; } + if (rsIsJMag < electrostatic.rMin()) + { + return true; + } + scalar chargeI = constPropI.siteCharges()[sI]; scalar chargeJ = constPropJ.siteCharges()[sJ]; @@ -491,6 +496,11 @@ inline bool Foam::moleculeCloud::evaluatePotentialLimit return true; } + if (rsRealsRefMag < electrostatic.rMin()) + { + return true; + } + scalar chargeReal = constPropReal.siteCharges()[sReal]; scalar chargeRef = constPropRef.siteCharges()[sRef]; From 64391aabbc1ed6596febef2a95306c730e45bf4f Mon Sep 17 00:00:00 2001 From: graham Date: Thu, 15 Jan 2009 20:16:12 +0000 Subject: [PATCH 03/21] Warning message update: molConfig -> mdInitialise --- .../molecule/moleculeCloud/moleculeCloudI.H | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudI.H b/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudI.H index e5320a4f3b..9af2e04f55 100644 --- a/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudI.H +++ b/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudI.H @@ -344,8 +344,8 @@ inline bool Foam::moleculeCloud::evaluatePotentialLimit << SMALL << ": mag separation = " << rsIsJMag << ". These may have been placed on top of each" - << " other by a rounding error in molConfig in" - << " parallel or a block filled with molecules " + << " other by a rounding error in mdInitialise in" + << " parallel or a block filled with molecules" << " twice. Removing one of the molecules." << endl; @@ -438,8 +438,8 @@ inline bool Foam::moleculeCloud::evaluatePotentialLimit << SMALL << ": mag separation = " << rsRealsRefMag << ". These may have been placed on top of each" - << " other by a rounding error in molConfig in" - << " parallel or a block filled with molecules " + << " other by a rounding error in mdInitialise in" + << " parallel or a block filled with molecules" << " twice. Removing one of the molecules." << endl; @@ -488,8 +488,8 @@ inline bool Foam::moleculeCloud::evaluatePotentialLimit << SMALL << ": mag separation = " << rsRealsRefMag << ". These may have been placed on top of each" - << " other by a rounding error in molConfig in" - << " parallel or a block filled with molecules " + << " other by a rounding error in mdInitialise in" + << " parallel or a block filled with molecules" << " twice. Removing one of the molecules." << endl; From 399e6a3baba08df0a4c13faefcb788579a1fecef Mon Sep 17 00:00:00 2001 From: graham Date: Fri, 16 Jan 2009 17:43:21 +0000 Subject: [PATCH 04/21] Modifying site virial contribution. --- .../molecule/moleculeCloud/moleculeCloudI.H | 46 +++++++++++++++---- 1 file changed, 38 insertions(+), 8 deletions(-) diff --git a/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudI.H b/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudI.H index 9af2e04f55..c4b8e5d89c 100644 --- a/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudI.H +++ b/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudI.H @@ -91,9 +91,18 @@ inline void Foam::moleculeCloud::evaluatePair molJ->potentialEnergy() += 0.5*potentialEnergy; - molI->rf() += rsIsJ * fsIsJ; + vector rIJ = molI->position() - molJ->position(); - molJ->rf() += rsIsJ * fsIsJ; + tensor virialContribution = + (rsIsJ * fsIsJ)*(rsIsJ & rIJ)/rsIsJMagSq; + + molI->rf() += virialContribution; + + molJ->rf() += virialContribution; + + // molI->rf() += rsIsJ * fsIsJ; + + // molJ->rf() += rsIsJ * fsIsJ; } } @@ -104,7 +113,7 @@ inline void Foam::moleculeCloud::evaluatePair scalar rsIsJMagSq = magSqr(rsIsJ); - if(pairPot.rCutMaxSqr(rsIsJMagSq)) + if(rsIsJMagSq <= electrostatic.rCutSqr()) { scalar rsIsJMag = mag(rsIsJ); @@ -126,9 +135,18 @@ inline void Foam::moleculeCloud::evaluatePair molJ->potentialEnergy() += 0.5*potentialEnergy; - molI->rf() += rsIsJ * fsIsJ; + vector rIJ = molI->position() - molJ->position(); - molJ->rf() += rsIsJ * fsIsJ; + tensor virialContribution = + (rsIsJ * fsIsJ)*(rsIsJ & rIJ)/rsIsJMagSq; + + molI->rf() += virialContribution; + + molJ->rf() += virialContribution; + + // molI->rf() += rsIsJ * fsIsJ; + + // molJ->rf() += rsIsJ * fsIsJ; } } } @@ -196,7 +214,13 @@ inline void Foam::moleculeCloud::evaluatePair molReal->potentialEnergy() += 0.5*potentialEnergy; - molReal->rf() += rsRealsRef * fsRealsRef; + vector rRealRef = molReal->position() - molRef->position(); + + molReal->rf() += (rsRealsRef * fsRealsRef) + *(rsRealsRef & rRealRef) + /rsRealsRefMagSq; + + // molReal->rf() += rsRealsRef * fsRealsRef; } } @@ -208,7 +232,7 @@ inline void Foam::moleculeCloud::evaluatePair scalar rsRealsRefMagSq = magSqr(rsRealsRef); - if(pairPot.rCutMaxSqr(rsRealsRefMagSq)) + if(rsRealsRefMagSq <= electrostatic.rCutSqr()) { scalar rsRealsRefMag = mag(rsRealsRef); @@ -227,7 +251,13 @@ inline void Foam::moleculeCloud::evaluatePair molReal->potentialEnergy() += 0.5*potentialEnergy; - molReal->rf() += rsRealsRef * fsRealsRef; + vector rRealRef = molReal->position() - molRef->position(); + + molReal->rf() += (rsRealsRef * fsRealsRef) + *(rsRealsRef & rRealRef) + /rsRealsRefMagSq; + + // molReal->rf() += rsRealsRef * fsRealsRef; } } } From 4f87d6053088835d3022f453f9016bf22b12efb7 Mon Sep 17 00:00:00 2001 From: graham Date: Tue, 20 Jan 2009 16:34:10 +0000 Subject: [PATCH 05/21] rename dictionary variable --- .../mdEquilibrationFoam/readmdEquilibrationDict.H | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/applications/solvers/molecularDynamics/mdEquilibrationFoam/readmdEquilibrationDict.H b/applications/solvers/molecularDynamics/mdEquilibrationFoam/readmdEquilibrationDict.H index 8289d014b0..b94745d9df 100644 --- a/applications/solvers/molecularDynamics/mdEquilibrationFoam/readmdEquilibrationDict.H +++ b/applications/solvers/molecularDynamics/mdEquilibrationFoam/readmdEquilibrationDict.H @@ -14,5 +14,5 @@ IOdictionary mdEquilibrationDict scalar targetTemperature = readScalar ( - mdEquilibrationDict.lookup("equilibrationTargetTemperature") + mdEquilibrationDict.lookup("targetTemperature") ); From 1b9bef5dd3aaaf9dbf8b33ab661c0d950f90125c Mon Sep 17 00:00:00 2001 From: graham Date: Thu, 22 Jan 2009 17:20:48 +0000 Subject: [PATCH 06/21] Adding vectors for global pi and tau as well as a triplet of orientation vectors to visualise rotational motion to written fields for post-processing. --- .../molecule/molecule/moleculeIO.C | 46 +++++++++++++++++++ 1 file changed, 46 insertions(+) diff --git a/src/lagrangian/molecularDynamics/molecule/molecule/moleculeIO.C b/src/lagrangian/molecularDynamics/molecule/molecule/moleculeIO.C index c1e99d8a0b..e01e5ec328 100644 --- a/src/lagrangian/molecularDynamics/molecule/molecule/moleculeIO.C +++ b/src/lagrangian/molecularDynamics/molecule/molecule/moleculeIO.C @@ -169,6 +169,37 @@ void Foam::molecule::writeFields(const moleculeCloud& mC) IOField