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

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
   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
  13     the Free Software Foundation, either version 3 of the License, or
  14     (at your 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, see <http://www.gnu.org/licenses/>.
  23
  24 \*----------------------------------------------------------------------------*/
  25
  26 #include "pairPotentialList.H"
  27 #include "OFstream.H"
  28 #include "Time.H"
  29
  30 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  31
  32 void Foam::pairPotentialList::readPairPotentialDict
  33 (
  34     const List<word>& idList,
  35     const dictionary& pairPotentialDict,
  36     const polyMesh& mesh
  37 )
  38 {
  39     Info<< nl << "Building pair potentials." << endl;
  40
  41     rCutMax_ = 0.0;
  42
  43     for (label a = 0; a < nIds_; ++a)
  44     {
  45         word idA = idList[a];
  46
  47         for (label b = a; b < nIds_; ++b)
  48         {
  49             word idB = idList[b];
  50
  51             word pairPotentialName;
  52
  53             if (a == b)
  54             {
  55                 if (pairPotentialDict.found(idA + "-" + idB))
  56                 {
  57                     pairPotentialName = idA + "-" + idB;
  58                 }
  59                 else
  60                 {
  61                     FatalErrorIn("pairPotentialList::buildPotentials") << nl
  62                         << "Pair pairPotential specification subDict "
  63                         << idA << "-" << idB << " not found"
  64                         << nl << abort(FatalError);
  65                 }
  66             }
  67             else
  68             {
  69                 if (pairPotentialDict.found(idA + "-" + idB))
  70                 {
  71                     pairPotentialName = idA + "-" + idB;
  72                 }
  73
  74                 else if (pairPotentialDict.found(idB + "-" + idA))
  75                 {
  76                     pairPotentialName = idB + "-" + idA;
  77                 }
  78
  79                 else
  80                 {
  81                     FatalErrorIn("pairPotentialList::buildPotentials") << nl
  82                         << "Pair pairPotential specification subDict "
  83                         << idA << "-" << idB << " or "
  84                         << idB << "-" << idA << " not found"
  85                         << nl << abort(FatalError);
  86                 }
  87
  88                 if
  89                 (
  90                     pairPotentialDict.found(idA+"-"+idB)
  91                  && pairPotentialDict.found(idB+"-"+idA)
  92                 )
  93                 {
  94                     FatalErrorIn("pairPotentialList::buildPotentials") << nl
  95                         << "Pair pairPotential specification subDict "
  96                         << idA << "-" << idB << " and "
  97                         << idB << "-" << idA << " found multiple definition"
  98                         << nl << abort(FatalError);
  99                 }
 100             }
 101
 102             (*this).set
 103             (
 104                 pairPotentialIndex(a, b),
 105                 pairPotential::New
 106                 (
 107                     pairPotentialName,
 108                     pairPotentialDict.subDict(pairPotentialName)
 109                 )
 110             );
 111
 112             if ((*this)[pairPotentialIndex(a, b)].rCut() > rCutMax_)
 113             {
 114                 rCutMax_ = (*this)[pairPotentialIndex(a, b)].rCut();
 115             }
 116
 117             if ((*this)[pairPotentialIndex(a, b)].writeTables())
 118             {
 119                 OFstream ppTabFile(mesh.time().path()/pairPotentialName);
 120
 121                 if
 122                 (
 123                     !(*this)[pairPotentialIndex(a, b)].writeEnergyAndForceTables
 124                     (
 125                         ppTabFile
 126                     )
 127                 )
 128                 {
 129                     FatalErrorIn("pairPotentialList::readPairPotentialDict")
 130                         << "Failed writing to "
 131                         << ppTabFile.name() << nl
 132                         << abort(FatalError);
 133                 }
 134             }
 135         }
 136     }
 137
 138     if (!pairPotentialDict.found("electrostatic"))
 139     {
 140         FatalErrorIn("pairPotentialList::buildPotentials") << nl
 141             << "Pair pairPotential specification subDict electrostatic"
 142             << nl << abort(FatalError);
 143     }
 144
 145     electrostaticPotential_ = pairPotential::New
 146     (
 147         "electrostatic",
 148         pairPotentialDict.subDict("electrostatic")
 149     );
 150
 151     if (electrostaticPotential_->rCut() > rCutMax_)
 152     {
 153         rCutMax_ = electrostaticPotential_->rCut();
 154     }
 155
 156     if (electrostaticPotential_->writeTables())
 157     {
 158         OFstream ppTabFile(mesh.time().path()/"electrostatic");
 159
 160         if (!electrostaticPotential_->writeEnergyAndForceTables(ppTabFile))
 161         {
 162             FatalErrorIn("pairPotentialList::readPairPotentialDict")
 163                 << "Failed writing to "
 164                 << ppTabFile.name() << nl
 165                 << abort(FatalError);
 166         }
 167     }
 168
 169     rCutMaxSqr_ = rCutMax_*rCutMax_;
 170 }
 171
 172
 173 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 174
 175 Foam::pairPotentialList::pairPotentialList()
 176 :
 177     PtrList<pairPotential>()
 178 {}
 179
 180
 181 Foam::pairPotentialList::pairPotentialList
 182 (
 183     const List<word>& idList,
 184     const dictionary& pairPotentialDict,
 185     const polyMesh& mesh
 186 )
 187 :
 188     PtrList<pairPotential>()
 189 {
 190     buildPotentials(idList, pairPotentialDict, mesh);
 191 }
 192
 193
 194 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 195
 196 Foam::pairPotentialList::~pairPotentialList()
 197 {}
 198
 199
 200 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 201
 202 void Foam::pairPotentialList::buildPotentials
 203 (
 204     const List<word>& idList,
 205     const dictionary& pairPotentialDict,
 206     const polyMesh& mesh
 207 )
 208 {
 209     setSize(((idList.size()*(idList.size() + 1))/2));
 210
 211     nIds_ = idList.size();
 212
 213     readPairPotentialDict(idList, pairPotentialDict, mesh);
 214 }
 215
 216
 217 const Foam::pairPotential& Foam::pairPotentialList::pairPotentialFunction
 218 (
 219     const label a,
 220     const label b
 221 ) const
 222 {
 223     return (*this)[pairPotentialIndex(a, b)];
 224 }
 225
 226
 227 bool Foam::pairPotentialList::rCutMaxSqr(const scalar rIJMagSqr) const
 228 {
 229     if (rIJMagSqr < rCutMaxSqr_)
 230     {
 231         return true;
 232     }
 233     else
 234     {
 235         return false;
 236     }
 237 }
 238
 239
 240 bool Foam::pairPotentialList::rCutSqr
 241 (
 242     const label a,
 243     const label b,
 244     const scalar rIJMagSqr
 245 ) const
 246 {
 247     if (rIJMagSqr < rCutSqr(a, b))
 248     {
 249         return true;
 250     }
 251     else
 252     {
 253         return false;
 254     }
 255 }
 256
 257
 258 Foam::scalar Foam::pairPotentialList::rMin
 259 (
 260     const label a,
 261     const label b
 262 ) const
 263 {
 264     return (*this)[pairPotentialIndex(a, b)].rMin();
 265 }
 266
 267
 268 Foam::scalar Foam::pairPotentialList::dr
 269 (
 270     const label a,
 271     const label b
 272 ) const
 273 {
 274     return (*this)[pairPotentialIndex(a, b)].dr();
 275 }
 276
 277
 278 Foam::scalar Foam::pairPotentialList::rCutSqr
 279 (
 280     const label a,
 281     const label b
 282 ) const
 283 {
 284     return (*this)[pairPotentialIndex(a, b)].rCutSqr();
 285 }
 286
 287
 288 Foam::scalar Foam::pairPotentialList::rCut
 289 (
 290     const label a,
 291     const label b
 292 ) const
 293 {
 294     return (*this)[pairPotentialIndex(a, b)].rCut();
 295 }
 296
 297
 298 Foam::scalar Foam::pairPotentialList::force
 299 (
 300     const label a,
 301     const label b,
 302     const scalar rIJMag
 303 ) const
 304 {
 305     scalar f = (*this)[pairPotentialIndex(a, b)].force(rIJMag);
 306
 307     return f;
 308 }
 309
 310
 311 Foam::scalar Foam::pairPotentialList::energy
 312 (
 313     const label a,
 314     const label b,
 315     const scalar rIJMag
 316 ) const
 317 {
 318     scalar e = (*this)[pairPotentialIndex(a, b)].energy(rIJMag);
 319
 320     return e;
 321 }
 322
 323
 324 // ************************************************************************* //