src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudI.H

   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 "constants.H"
  27
  28 using namespace Foam::constant;
  29
  30 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  31
  32 inline void Foam::moleculeCloud::evaluatePair
  33 (
  34     molecule& molI,
  35     molecule& molJ
  36 )
  37 {
  38     const pairPotentialList& pairPot = pot_.pairPotentials();
  39
  40     const pairPotential& electrostatic = pairPot.electrostatic();
  41
  42     label idI = molI.id();
  43
  44     label idJ = molJ.id();
  45
  46     const molecule::constantProperties& constPropI(constProps(idI));
  47
  48     const molecule::constantProperties& constPropJ(constProps(idJ));
  49
  50     List<label> siteIdsI = constPropI.siteIds();
  51
  52     List<label> siteIdsJ = constPropJ.siteIds();
  53
  54     List<bool> pairPotentialSitesI = constPropI.pairPotentialSites();
  55
  56     List<bool> electrostaticSitesI = constPropI.electrostaticSites();
  57
  58     List<bool> pairPotentialSitesJ = constPropJ.pairPotentialSites();
  59
  60     List<bool> electrostaticSitesJ = constPropJ.electrostaticSites();
  61
  62     forAll(siteIdsI, sI)
  63     {
  64         label idsI(siteIdsI[sI]);
  65
  66         forAll(siteIdsJ, sJ)
  67         {
  68             label idsJ(siteIdsJ[sJ]);
  69
  70             if (pairPotentialSitesI[sI] && pairPotentialSitesJ[sJ])
  71             {
  72                 vector rsIsJ =
  73                     molI.sitePositions()[sI] - molJ.sitePositions()[sJ];
  74
  75                 scalar rsIsJMagSq = magSqr(rsIsJ);
  76
  77                 if (pairPot.rCutSqr(idsI, idsJ, rsIsJMagSq))
  78                 {
  79                     scalar rsIsJMag = mag(rsIsJ);
  80
  81                     vector fsIsJ =
  82                         (rsIsJ/rsIsJMag)
  83                        *pairPot.force(idsI, idsJ, rsIsJMag);
  84
  85                     molI.siteForces()[sI] += fsIsJ;
  86
  87                     molJ.siteForces()[sJ] += -fsIsJ;
  88
  89                     scalar potentialEnergy
  90                     (
  91                         pairPot.energy(idsI, idsJ, rsIsJMag)
  92                     );
  93
  94                     molI.potentialEnergy() += 0.5*potentialEnergy;
  95
  96                     molJ.potentialEnergy() += 0.5*potentialEnergy;
  97
  98                     vector rIJ = molI.position() - molJ.position();
  99
 100                     tensor virialContribution =
 101                         (rsIsJ*fsIsJ)*(rsIsJ & rIJ)/rsIsJMagSq;
 102
 103                     molI.rf() += virialContribution;
 104
 105                     molJ.rf() += virialContribution;
 106                 }
 107             }
 108
 109             if (electrostaticSitesI[sI] && electrostaticSitesJ[sJ])
 110             {
 111                 vector rsIsJ =
 112                 molI.sitePositions()[sI] - molJ.sitePositions()[sJ];
 113
 114                 scalar rsIsJMagSq = magSqr(rsIsJ);
 115
 116                 if (rsIsJMagSq <= electrostatic.rCutSqr())
 117                 {
 118                     scalar rsIsJMag = mag(rsIsJ);
 119
 120                     scalar chargeI = constPropI.siteCharges()[sI];
 121
 122                     scalar chargeJ = constPropJ.siteCharges()[sJ];
 123
 124                     vector fsIsJ =
 125                         (rsIsJ/rsIsJMag)
 126                        *chargeI*chargeJ*electrostatic.force(rsIsJMag);
 127
 128                     molI.siteForces()[sI] += fsIsJ;
 129
 130                     molJ.siteForces()[sJ] += -fsIsJ;
 131
 132                     scalar potentialEnergy =
 133                         chargeI*chargeJ
 134                        *electrostatic.energy(rsIsJMag);
 135
 136                     molI.potentialEnergy() += 0.5*potentialEnergy;
 137
 138                     molJ.potentialEnergy() += 0.5*potentialEnergy;
 139
 140                     vector rIJ = molI.position() - molJ.position();
 141
 142                     tensor virialContribution =
 143                         (rsIsJ*fsIsJ)*(rsIsJ & rIJ)/rsIsJMagSq;
 144
 145                     molI.rf() += virialContribution;
 146
 147                     molJ.rf() += virialContribution;
 148                 }
 149             }
 150         }
 151     }
 152 }
 153
 154
 155 inline bool Foam::moleculeCloud::evaluatePotentialLimit
 156 (
 157     molecule& molI,
 158     molecule& molJ
 159 ) const
 160 {
 161     const pairPotentialList& pairPot = pot_.pairPotentials();
 162
 163     const pairPotential& electrostatic = pairPot.electrostatic();
 164
 165     label idI = molI.id();
 166
 167     label idJ = molJ.id();
 168
 169     const molecule::constantProperties& constPropI(constProps(idI));
 170
 171     const molecule::constantProperties& constPropJ(constProps(idJ));
 172
 173     List<label> siteIdsI = constPropI.siteIds();
 174
 175     List<label> siteIdsJ = constPropJ.siteIds();
 176
 177     List<bool> pairPotentialSitesI = constPropI.pairPotentialSites();
 178
 179     List<bool> electrostaticSitesI = constPropI.electrostaticSites();
 180
 181     List<bool> pairPotentialSitesJ = constPropJ.pairPotentialSites();
 182
 183     List<bool> electrostaticSitesJ = constPropJ.electrostaticSites();
 184
 185     forAll(siteIdsI, sI)
 186     {
 187         label idsI(siteIdsI[sI]);
 188
 189         forAll(siteIdsJ, sJ)
 190         {
 191             label idsJ(siteIdsJ[sJ]);
 192
 193             if (pairPotentialSitesI[sI] && pairPotentialSitesJ[sJ])
 194             {
 195                 vector rsIsJ =
 196                     molI.sitePositions()[sI] - molJ.sitePositions()[sJ];
 197
 198                 scalar rsIsJMagSq = magSqr(rsIsJ);
 199
 200                 if (pairPot.rCutSqr(idsI, idsJ, rsIsJMagSq))
 201                 {
 202                     scalar rsIsJMag = mag(rsIsJ);
 203
 204                     // Guard against pairPot.energy being evaluated
 205                     // if rIJMag < SMALL. A floating point exception will
 206                     // happen otherwise.
 207
 208                     if (rsIsJMag < SMALL)
 209                     {
 210                         WarningIn("moleculeCloud::removeHighEnergyOverlaps()")
 211                             << "Molecule site pair closer than "
 212                             << SMALL
 213                             << ": mag separation = " << rsIsJMag
 214                             << ". These may have been placed on top of each"
 215                             << " other by a rounding error in mdInitialise in"
 216                             << " parallel or a block filled with molecules"
 217                             << " twice. Removing one of the molecules."
 218                             << endl;
 219
 220                         return true;
 221                     }
 222
 223                     // Guard against pairPot.energy being evaluated if rIJMag <
 224                     // rMin.  A tabulation lookup error will occur otherwise.
 225
 226                     if (rsIsJMag < pairPot.rMin(idsI, idsJ))
 227                     {
 228                         return true;
 229                     }
 230
 231                     if
 232                     (
 233                         mag(pairPot.energy(idsI, idsJ, rsIsJMag))
 234                       > pot_.potentialEnergyLimit()
 235                     )
 236                     {
 237                         return true;
 238                     };
 239                 }
 240             }
 241
 242             if (electrostaticSitesI[sI] && electrostaticSitesJ[sJ])
 243             {
 244                 vector rsIsJ =
 245                     molI.sitePositions()[sI] - molJ.sitePositions()[sJ];
 246
 247                 scalar rsIsJMagSq = magSqr(rsIsJ);
 248
 249                 if (pairPot.rCutMaxSqr(rsIsJMagSq))
 250                 {
 251                     scalar rsIsJMag = mag(rsIsJ);
 252
 253                     // Guard against pairPot.energy being evaluated
 254                     // if rIJMag < SMALL. A floating point exception will
 255                     // happen otherwise.
 256
 257                     if (rsIsJMag < SMALL)
 258                     {
 259                         WarningIn("moleculeCloud::removeHighEnergyOverlaps()")
 260                             << "Molecule site pair closer than "
 261                             << SMALL
 262                             << ": mag separation = " << rsIsJMag
 263                             << ". These may have been placed on top of each"
 264                             << " other by a rounding error in mdInitialise in"
 265                             << " parallel or a block filled with molecules"
 266                             << " twice. Removing one of the molecules."
 267                             << endl;
 268
 269                         return true;
 270                     }
 271
 272                     if (rsIsJMag < electrostatic.rMin())
 273                     {
 274                         return true;
 275                     }
 276
 277                     scalar chargeI = constPropI.siteCharges()[sI];
 278
 279                     scalar chargeJ = constPropJ.siteCharges()[sJ];
 280
 281                     if
 282                     (
 283                         mag(chargeI*chargeJ*electrostatic.energy(rsIsJMag))
 284                       > pot_.potentialEnergyLimit()
 285                     )
 286                     {
 287                         return true;
 288                     };
 289                 }
 290             }
 291         }
 292     }
 293
 294     return false;
 295 }
 296
 297
 298 inline Foam::vector Foam::moleculeCloud::equipartitionLinearVelocity
 299 (
 300     scalar temperature,
 301     scalar mass
 302 )
 303 {
 304     return sqrt(physicoChemical::k.value()*temperature/mass)*vector
 305     (
 306         rndGen_.GaussNormal(),
 307         rndGen_.GaussNormal(),
 308         rndGen_.GaussNormal()
 309     );
 310 }
 311
 312
 313 inline Foam::vector Foam::moleculeCloud::equipartitionAngularMomentum
 314 (
 315     scalar temperature,
 316     const molecule::constantProperties& cP
 317 )
 318 {
 319     scalar sqrtKbT = sqrt(physicoChemical::k.value()*temperature);
 320
 321     if (cP.linearMolecule())
 322     {
 323         return sqrtKbT*vector
 324         (
 325             0.0,
 326             sqrt(cP.momentOfInertia().yy())*rndGen_.GaussNormal(),
 327             sqrt(cP.momentOfInertia().zz())*rndGen_.GaussNormal()
 328         );
 329     }
 330     else
 331     {
 332         return sqrtKbT*vector
 333         (
 334             sqrt(cP.momentOfInertia().xx())*rndGen_.GaussNormal(),
 335             sqrt(cP.momentOfInertia().yy())*rndGen_.GaussNormal(),
 336             sqrt(cP.momentOfInertia().zz())*rndGen_.GaussNormal()
 337         );
 338     }
 339 }
 340
 341
 342 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 343
 344 inline const Foam::polyMesh& Foam::moleculeCloud::mesh() const
 345 {
 346     return mesh_;
 347 }
 348
 349
 350 inline const Foam::potential& Foam::moleculeCloud::pot() const
 351 {
 352     return pot_;
 353 }
 354
 355
 356 inline const Foam::List<Foam::DynamicList<Foam::molecule*> >&
 357     Foam::moleculeCloud::cellOccupancy() const
 358 {
 359     return cellOccupancy_;
 360 }
 361
 362
 363 inline const Foam::InteractionLists<Foam::molecule>&
 364     Foam::moleculeCloud::il() const
 365 {
 366     return il_;
 367 }
 368
 369
 370 inline const Foam::List<Foam::molecule::constantProperties>
 371     Foam::moleculeCloud::constProps() const
 372 {
 373     return constPropList_;
 374 }
 375
 376
 377 inline const Foam::molecule::constantProperties&
 378     Foam::moleculeCloud::constProps(label id) const
 379 {
 380     return constPropList_[id];
 381 }
 382
 383
 384 inline Foam::Random& Foam::moleculeCloud::rndGen()
 385 {
 386     return rndGen_;
 387 }
 388
 389
 390 // ************************************************************************* //