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

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