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) 1991-2008 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 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  31
  32 namespace Foam
  33 {
  34
  35 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  36
  37 defineTypeNameAndDebug(pairPotential, 0);
  38 defineRunTimeSelectionTable(pairPotential, dictionary);
  39
  40 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  41
  42
  43 void Foam::pairPotential::scaleEnergy(scalar& e, const scalar r) const
  44 {
  45     if (!esfPtr_)
  46     {
  47         esfPtr_ = energyScalingFunction::New
  48         (
  49             name_, pairPotentialProperties_, *this
  50         ).ptr();
  51     }
  52
  53     esfPtr_->scaleEnergy(e, r);
  54 }
  55
  56
  57 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  58
  59 Foam::pairPotential::pairPotential
  60 (
  61     const word& name,
  62     const dictionary& pairPotentialProperties
  63 )
  64 :
  65     name_(name),
  66     pairPotentialProperties_(pairPotentialProperties),
  67     rCut_(readScalar(pairPotentialProperties_.lookup("rCut"))),
  68     rCutSqr_(rCut_*rCut_),
  69     rMin_(readScalar(pairPotentialProperties_.lookup("rMin"))),
  70     dr_(readScalar(pairPotentialProperties_.lookup("dr"))),
  71     forceLookup_(0),
  72     energyLookup_(0),
  73     esfPtr_(NULL),
  74     writeTables_(Switch(pairPotentialProperties_.lookup("writeTables")))
  75 {}
  76
  77
  78 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
  79
  80 void Foam::pairPotential::setLookupTables()
  81 {
  82     label N = label((rCut_ - rMin_)/dr_) + 1;
  83
  84     forceLookup_.setSize(N);
  85
  86     energyLookup_.setSize(N);
  87
  88     forAll(forceLookup_, k)
  89     {
  90         energyLookup_[k] = scaledEnergy(k*dr_ + rMin_);
  91
  92         forceLookup_[k] = -energyDerivative((k*dr_ + rMin_), true);
  93     }
  94 }
  95
  96
  97 Foam::scalar Foam::pairPotential::forceLookup(const scalar r) const
  98 {
  99     scalar k_rIJ = (r - rMin_)/dr_;
 100
 101     label k(k_rIJ);
 102
 103     if (k < 0)
 104     {
 105         FatalErrorIn("pairPotential.C") << nl
 106             << "r less than rMin" << nl
 107             << abort(FatalError);
 108     }
 109
 110     scalar f =
 111         (k_rIJ - k)*forceLookup_[k+1]
 112       + (k + 1 - k_rIJ)*forceLookup_[k];
 113
 114     return f;
 115 }
 116
 117
 118 Foam::List< Foam::Pair< Foam::scalar > >
 119 Foam::pairPotential::forceTable() const
 120 {
 121     List<Pair<scalar> > forceTab(forceLookup_.size());
 122
 123     forAll(forceLookup_,k)
 124     {
 125         forceTab[k].first() = rMin_ + k*dr_;
 126
 127         forceTab[k].second() = forceLookup_[k];
 128     }
 129
 130     return forceTab;
 131 }
 132
 133
 134 Foam::scalar Foam::pairPotential::energyLookup(const scalar r) const
 135 {
 136     scalar k_rIJ = (r - rMin_)/dr_;
 137
 138     label k(k_rIJ);
 139
 140     if (k < 0)
 141     {
 142         FatalErrorIn("pairPotential.C") << nl
 143             << "r less than rMin" << nl
 144             << abort(FatalError);
 145     }
 146
 147     scalar e =
 148         (k_rIJ - k)*energyLookup_[k+1]
 149       + (k + 1 - k_rIJ)*energyLookup_[k];
 150
 151     return e;
 152 }
 153
 154
 155 Foam::List< Foam::Pair< Foam::scalar > >
 156     Foam::pairPotential::energyTable() const
 157 {
 158     List<Pair<scalar> > energyTab(energyLookup_.size());
 159
 160     forAll(energyLookup_,k)
 161     {
 162         energyTab[k].first() = rMin_ + k*dr_;
 163
 164         energyTab[k].second() = energyLookup_[k];
 165     }
 166
 167     return energyTab;
 168 }
 169
 170
 171 Foam::scalar Foam::pairPotential::scaledEnergy(const scalar r) const
 172 {
 173     scalar e = unscaledEnergy(r);
 174
 175     scaleEnergy(e, r);
 176
 177     return e;
 178 }
 179
 180
 181 Foam::scalar Foam::pairPotential::energyDerivative
 182 (
 183     const scalar r,
 184     const bool scaledEnergyDerivative
 185 ) const
 186 {
 187     // Local quadratic fit to energy: E = a0 + a1*r + a2*r^2
 188     // Differentiate to give f = -dE/dr = -a1 - 2*a2*r
 189
 190     scalar ra = r - dr_;
 191     scalar rf = r;
 192     scalar rb = r + dr_;
 193
 194     scalar Ea, Ef, Eb;
 195
 196     if (scaledEnergyDerivative)
 197     {
 198         Ea = scaledEnergy(ra);
 199         Ef = scaledEnergy(rf);
 200         Eb = scaledEnergy(rb);
 201     }
 202     else
 203     {
 204         Ea = unscaledEnergy(ra);
 205         Ef = unscaledEnergy(rf);
 206         Eb = unscaledEnergy(rb);
 207     }
 208
 209     scalar denominator = (ra - rf)*(ra - rb)*(rf - rb);
 210
 211     scalar a1 =
 212     (
 213         rb*rb*(Ea - Ef) + ra*ra*(Ef - Eb) + rf*rf*(Eb - Ea)
 214     ) / denominator;
 215
 216     scalar a2 =
 217     (
 218         rb*(Ef - Ea) + rf*(Ea - Eb) + ra*(Eb - Ef)
 219     ) / denominator;
 220
 221     return a1 + 2.0*a2*r;
 222 }
 223
 224
 225 bool Foam::pairPotential::read(const dictionary& pairPotentialProperties)
 226 {
 227     pairPotentialProperties_ = pairPotentialProperties;
 228
 229     return true;
 230 }
 231
 232
 233 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 234
 235 } // End namespace Foam
 236
 237 // ************************************************************************* //