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

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | foam-extend: Open Source CFD
   4    \\    /   O peration     | Version:     3.2
   5     \\  /    A nd           | Web:         http://www.foam-extend.org
   6      \\/     M anipulation  | For copyright notice see file Copyright
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of foam-extend.
  10
  11     foam-extend is free software: you can redistribute it and/or modify it
  12     under the terms of the GNU General Public License as published by the
  13     Free Software Foundation, either version 3 of the License, or (at your
  14     option) any later version.
  15
  16     foam-extend is distributed in the hope that it will be useful, but
  17     WITHOUT ANY WARRANTY; without even the implied warranty of
  18     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19     General Public License for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 \*---------------------------------------------------------------------------*/
  25
  26 #include "pairPotential.H"
  27 #include "energyScalingFunction.H"
  28
  29 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  30
  31 namespace Foam
  32 {
  33     defineTypeNameAndDebug(pairPotential, 0);
  34     defineRunTimeSelectionTable(pairPotential, dictionary);
  35 }
  36
  37
  38 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  39
  40
  41 void Foam::pairPotential::scaleEnergy(scalar& e, const scalar r) const
  42 {
  43     if (!esfPtr_)
  44     {
  45         esfPtr_ = energyScalingFunction::New
  46         (
  47             name_, pairPotentialProperties_, *this
  48         ).ptr();
  49     }
  50
  51     esfPtr_->scaleEnergy(e, r);
  52 }
  53
  54
  55 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  56
  57 Foam::pairPotential::pairPotential
  58 (
  59     const word& name,
  60     const dictionary& pairPotentialProperties
  61 )
  62 :
  63     name_(name),
  64     pairPotentialProperties_(pairPotentialProperties),
  65     rCut_(readScalar(pairPotentialProperties_.lookup("rCut"))),
  66     rCutSqr_(rCut_*rCut_),
  67     rMin_(readScalar(pairPotentialProperties_.lookup("rMin"))),
  68     dr_(readScalar(pairPotentialProperties_.lookup("dr"))),
  69     forceLookup_(0),
  70     energyLookup_(0),
  71     esfPtr_(NULL),
  72     writeTables_(Switch(pairPotentialProperties_.lookup("writeTables")))
  73 {}
  74
  75
  76 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
  77
  78 void Foam::pairPotential::setLookupTables()
  79 {
  80     label N = label((rCut_ - rMin_)/dr_) + 1;
  81
  82     forceLookup_.setSize(N);
  83
  84     energyLookup_.setSize(N);
  85
  86     forAll(forceLookup_, k)
  87     {
  88         energyLookup_[k] = scaledEnergy(k*dr_ + rMin_);
  89
  90         forceLookup_[k] = -energyDerivative((k*dr_ + rMin_), true);
  91     }
  92 }
  93
  94
  95 Foam::scalar Foam::pairPotential::force(const scalar r) const
  96 {
  97     scalar k_rIJ = (r - rMin_)/dr_;
  98
  99     label k = label(k_rIJ);
 100
 101     if (k < 0)
 102     {
 103         FatalErrorIn("pairPotential.C") << nl
 104             << "r less than rMin in pair potential " << name_ << nl
 105             << abort(FatalError);
 106     }
 107
 108     scalar f =
 109         (k_rIJ - k)*forceLookup_[k+1]
 110       + (k + 1 - k_rIJ)*forceLookup_[k];
 111
 112     return f;
 113 }
 114
 115
 116 Foam::List< Foam::Pair< Foam::scalar > >
 117 Foam::pairPotential::forceTable() const
 118 {
 119     List<Pair<scalar> > forceTab(forceLookup_.size());
 120
 121     forAll(forceLookup_,k)
 122     {
 123         forceTab[k].first() = rMin_ + k*dr_;
 124
 125         forceTab[k].second() = forceLookup_[k];
 126     }
 127
 128     return forceTab;
 129 }
 130
 131
 132 Foam::scalar Foam::pairPotential::energy(const scalar r) const
 133 {
 134     scalar k_rIJ = (r - rMin_)/dr_;
 135
 136     label k = label(k_rIJ);
 137
 138     if (k < 0)
 139     {
 140         FatalErrorIn("pairPotential.C") << nl
 141             << "r less than rMin in pair potential " << name_ << nl
 142             << abort(FatalError);
 143     }
 144
 145     scalar e =
 146         (k_rIJ - k)*energyLookup_[k+1]
 147       + (k + 1 - k_rIJ)*energyLookup_[k];
 148
 149     return e;
 150 }
 151
 152
 153 Foam::List< Foam::Pair< Foam::scalar > >
 154     Foam::pairPotential::energyTable() const
 155 {
 156     List<Pair<scalar> > energyTab(energyLookup_.size());
 157
 158     forAll(energyLookup_,k)
 159     {
 160         energyTab[k].first() = rMin_ + k*dr_;
 161
 162         energyTab[k].second() = energyLookup_[k];
 163     }
 164
 165     return energyTab;
 166 }
 167
 168
 169 Foam::scalar Foam::pairPotential::scaledEnergy(const scalar r) const
 170 {
 171     scalar e = unscaledEnergy(r);
 172
 173     scaleEnergy(e, r);
 174
 175     return e;
 176 }
 177
 178
 179 Foam::scalar Foam::pairPotential::energyDerivative
 180 (
 181     const scalar r,
 182     const bool scaledEnergyDerivative
 183 ) const
 184 {
 185     // Local quadratic fit to energy: E = a0 + a1*r + a2*r^2
 186     // Differentiate to give f = -dE/dr = -a1 - 2*a2*r
 187
 188     scalar ra = r - dr_;
 189     scalar rf = r;
 190     scalar rb = r + dr_;
 191
 192     scalar Ea, Ef, Eb;
 193
 194     if (scaledEnergyDerivative)
 195     {
 196         Ea = scaledEnergy(ra);
 197         Ef = scaledEnergy(rf);
 198         Eb = scaledEnergy(rb);
 199     }
 200     else
 201     {
 202         Ea = unscaledEnergy(ra);
 203         Ef = unscaledEnergy(rf);
 204         Eb = unscaledEnergy(rb);
 205     }
 206
 207     scalar denominator = (ra - rf)*(ra - rb)*(rf - rb);
 208
 209     scalar a1 =
 210     (
 211         rb*rb*(Ea - Ef) + ra*ra*(Ef - Eb) + rf*rf*(Eb - Ea)
 212     ) / denominator;
 213
 214     scalar a2 =
 215     (
 216         rb*(Ef - Ea) + rf*(Ea - Eb) + ra*(Eb - Ef)
 217     ) / denominator;
 218
 219     return a1 + 2.0*a2*r;
 220 }
 221
 222
 223 bool Foam::pairPotential::read(const dictionary& pairPotentialProperties)
 224 {
 225     pairPotentialProperties_ = pairPotentialProperties;
 226
 227     return true;
 228 }
 229
 230
 231 // ************************************************************************* //