src/lagrangian/molecularDynamics/potential/potential/potential.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 "potential.H"
  28
  29 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  30
  31 void Foam::potential::setSiteIdList(const IOdictionary& moleculePropertiesDict)
  32 {
  33     DynamicList<word> siteIdList;
  34
  35     DynamicList<word> pairPotentialSiteIdList;
  36
  37     forAll(idList_, i)
  38     {
  39         const word& id(idList_[i]);
  40
  41         if (!moleculePropertiesDict.found(id))
  42         {
  43             FatalErrorIn("potential.C") << nl
  44                 << id << " molecule subDict not found"
  45                 << nl << abort(FatalError);
  46         }
  47
  48         const dictionary& molDict(moleculePropertiesDict.subDict(id));
  49
  50         List<word> siteIdNames = molDict.lookup("siteIds");
  51
  52         forAll(siteIdNames, sI)
  53         {
  54             const word& siteId = siteIdNames[sI];
  55
  56             if(findIndex(siteIdList, siteId) == -1)
  57             {
  58                 siteIdList.append(siteId);
  59             }
  60         }
  61
  62         List<word> pairPotSiteIds = molDict.lookup("pairPotentialSiteIds");
  63
  64         forAll(pairPotSiteIds, sI)
  65         {
  66             const word& siteId = pairPotSiteIds[sI];
  67
  68             if(findIndex(siteIdNames, siteId) == -1)
  69             {
  70                 FatalErrorIn("potential.C") << nl
  71                     << siteId << " in pairPotentialSiteIds is not in siteIds: "
  72                     << siteIdNames << nl << abort(FatalError);
  73             }
  74
  75             if(findIndex(pairPotentialSiteIdList, siteId) == -1)
  76             {
  77                 pairPotentialSiteIdList.append(siteId);
  78             }
  79         }
  80     }
  81
  82     nPairPotIds_ = pairPotentialSiteIdList.size();
  83
  84     forAll(siteIdList, aSIN)
  85     {
  86         const word& siteId = siteIdList[aSIN];
  87
  88         if(findIndex(pairPotentialSiteIdList, siteId) == -1)
  89         {
  90             pairPotentialSiteIdList.append(siteId);
  91         }
  92     }
  93
  94     siteIdList_.transfer(pairPotentialSiteIdList.shrink());
  95 }
  96
  97
  98 void Foam::potential::potential::readPotentialDict()
  99 {
 100     Info<< nl <<  "Reading potential dictionary:" << endl;
 101
 102     IOdictionary idListDict
 103     (
 104         IOobject
 105         (
 106             "idList",
 107             mesh_.time().constant(),
 108             mesh_,
 109             IOobject::MUST_READ,
 110             IOobject::NO_WRITE
 111         )
 112     );
 113
 114     idList_ = List<word>(idListDict.lookup("idList"));
 115
 116     IOdictionary moleculePropertiesDict
 117     (
 118         IOobject
 119         (
 120             "moleculeProperties",
 121             mesh_.time().constant(),
 122             mesh_,
 123             IOobject::MUST_READ,
 124             IOobject::NO_WRITE,
 125             false
 126         )
 127     );
 128
 129     setSiteIdList(moleculePropertiesDict);
 130
 131     List<word> pairPotentialSiteIdList
 132     (
 133         SubList<word>(siteIdList_, nPairPotIds_)
 134     );
 135
 136     Info<< nl << "Unique site ids found: " << siteIdList_
 137         << nl << "Site Ids requiring a pair potential: "
 138         << pairPotentialSiteIdList
 139         << endl;
 140
 141     List<word> tetherSiteIdList(0);
 142
 143     if (idListDict.found("tetherSiteIdList"))
 144     {
 145         tetherSiteIdList = List<word>(idListDict.lookup("tetherSiteIdList"));
 146     }
 147
 148     IOdictionary potentialDict
 149     (
 150         IOobject
 151         (
 152             "potentialDict",
 153             mesh_.time().system(),
 154             mesh_,
 155             IOobject::MUST_READ,
 156             IOobject::NO_WRITE
 157         )
 158     );
 159
 160     potentialEnergyLimit_ = readScalar
 161     (
 162         potentialDict.lookup("potentialEnergyLimit")
 163     );
 164
 165     if (potentialDict.found("removalOrder"))
 166     {
 167         List<word> remOrd = potentialDict.lookup("removalOrder");
 168
 169         removalOrder_.setSize(remOrd.size());
 170
 171         forAll(removalOrder_, rO)
 172         {
 173             removalOrder_[rO] = findIndex(idList_, remOrd[rO]);
 174
 175             if (removalOrder_[rO] == -1)
 176             {
 177                 FatalErrorIn("potentials.C") << nl
 178                     << "removalOrder entry: " << remOrd[rO]
 179                     << " not found in idList."
 180                     << nl << abort(FatalError);
 181             }
 182         }
 183     }
 184
 185     // *************************************************************************
 186     // Pair potentials
 187
 188     if (!potentialDict.found("pair"))
 189     {
 190         FatalErrorIn("potentials.C") << nl
 191             << "pair potential specification subDict not found"
 192             << abort(FatalError);
 193     }
 194
 195     const dictionary& pairDict = potentialDict.subDict("pair");
 196
 197     pairPotentials_.buildPotentials
 198     (
 199         pairPotentialSiteIdList,
 200         pairDict,
 201         mesh_
 202     );
 203
 204     // *************************************************************************
 205     // Tether potentials
 206
 207     if (tetherSiteIdList.size())
 208     {
 209         if (!potentialDict.found("tether"))
 210         {
 211             FatalErrorIn("potential.C") << nl
 212                 << "tether potential specification subDict not found"
 213                 << abort(FatalError);
 214         }
 215
 216         const dictionary& tetherDict = potentialDict.subDict("tether");
 217
 218         tetherPotentials_.buildPotentials
 219         (
 220             siteIdList_,
 221             tetherDict,
 222             tetherSiteIdList
 223         );
 224     }
 225
 226     // *************************************************************************
 227     // External Forces
 228
 229     gravity_ = vector::zero;
 230
 231     if (potentialDict.found("external"))
 232     {
 233
 234         Info << nl << "Reading external forces:" << endl;
 235
 236         const dictionary& externalDict = potentialDict.subDict("external");
 237
 238         // *********************************************************************
 239         // gravity
 240
 241         if (externalDict.found("gravity"))
 242         {
 243             gravity_ = externalDict.lookup("gravity");
 244         }
 245     }
 246
 247     Info << nl << tab << "gravity = " << gravity_ << endl;
 248 }
 249
 250
 251 void Foam::potential::potential::readMdInitialiseDict
 252 (
 253     const IOdictionary& mdInitialiseDict,
 254     IOdictionary& idListDict
 255 )
 256 {
 257     IOdictionary moleculePropertiesDict
 258     (
 259         IOobject
 260         (
 261             "moleculeProperties",
 262             mesh_.time().constant(),
 263             mesh_,
 264             IOobject::MUST_READ,
 265             IOobject::NO_WRITE,
 266             false
 267         )
 268     );
 269
 270     DynamicList<word> idList;
 271
 272     DynamicList<word> tetherSiteIdList;
 273
 274     forAll(mdInitialiseDict.toc(), zone)
 275     {
 276         const dictionary& zoneDict = mdInitialiseDict.subDict
 277         (
 278             mdInitialiseDict.toc()[zone]
 279         );
 280
 281         List<word> latticeIds
 282         (
 283             zoneDict.lookup("latticeIds")
 284         );
 285
 286         forAll(latticeIds, i)
 287         {
 288             const word& id = latticeIds[i];
 289
 290             if(!moleculePropertiesDict.found(id))
 291             {
 292                 FatalErrorIn("potential.C") << nl
 293                     << "Molecule type "
 294                     << id
 295                     << " not found in moleculeProperties dictionary."
 296                     << nl
 297                     << abort(FatalError);
 298             }
 299
 300             if (findIndex(idList,id) == -1)
 301             {
 302                 idList.append(id);
 303             }
 304         }
 305
 306         List<word> tetherSiteIds
 307         (
 308             zoneDict.lookup("tetherSiteIds")
 309         );
 310
 311         forAll(tetherSiteIds, t)
 312         {
 313             const word& tetherSiteId = tetherSiteIds[t];
 314
 315             bool idFound = false;
 316
 317             forAll(latticeIds, i)
 318             {
 319                 if (idFound)
 320                 {
 321                     break;
 322                 }
 323
 324                 const word& id = latticeIds[i];
 325
 326                 List<word> siteIds
 327                 (
 328                     moleculePropertiesDict.subDict(id).lookup("siteIds")
 329                 );
 330
 331                 if (findIndex(siteIds, tetherSiteId) != -1)
 332                 {
 333                     idFound = true;
 334                 }
 335             }
 336
 337             if (idFound)
 338             {
 339                 tetherSiteIdList.append(tetherSiteId);
 340             }
 341             else
 342             {
 343                 FatalErrorIn("potential.C") << nl
 344                     << "Tether id  "
 345                     << tetherSiteId
 346                     << " not found as a site of any molecule in zone."
 347                     << nl
 348                     << abort(FatalError);
 349             }
 350         }
 351     }
 352
 353     idList_.transfer(idList);
 354
 355     tetherSiteIdList.shrink();
 356
 357     idListDict.add("idList", idList_);
 358
 359     idListDict.add("tetherSiteIdList", tetherSiteIdList);
 360
 361     setSiteIdList(moleculePropertiesDict);
 362 }
 363
 364 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 365
 366 Foam::potential::potential(const polyMesh& mesh)
 367 :
 368     mesh_(mesh)
 369 {
 370     readPotentialDict();
 371 }
 372
 373
 374 Foam::potential::potential
 375 (
 376     const polyMesh& mesh,
 377     const IOdictionary& mdInitialiseDict,
 378     IOdictionary& idListDict
 379 )
 380 :
 381     mesh_(mesh)
 382 {
 383     readMdInitialiseDict(mdInitialiseDict, idListDict);
 384 }
 385
 386 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 387
 388 Foam::potential::~potential()
 389 {}
 390
 391
 392 // ************************************************************************* //