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