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) 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 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  28
  29 inline void Foam::moleculeCloud::evaluatePair
  30 (
  31     molecule* molI,
  32     molecule* molJ
  33 )
  34 {
  35     const pairPotentialList& pairPot = pot_.pairPotentials();
  36
  37     const pairPotential& electrostatic = pairPot.electrostatic();
  38
  39     label idI = molI->id();
  40
  41     label idJ = molJ->id();
  42
  43     const molecule::constantProperties& constPropI(constProps(idI));
  44
  45     const molecule::constantProperties& constPropJ(constProps(idJ));
  46
  47     List<label> siteIdsI = constPropI.siteIds();
  48
  49     List<label> siteIdsJ = constPropJ.siteIds();
  50
  51     List<bool> pairPotentialSitesI = constPropI.pairPotentialSites();
  52
  53     List<bool> electrostaticSitesI = constPropI.electrostaticSites();
  54
  55     List<bool> pairPotentialSitesJ = constPropJ.pairPotentialSites();
  56
  57     List<bool> electrostaticSitesJ = constPropJ.electrostaticSites();
  58
  59     forAll(siteIdsI, sI)
  60     {
  61         label idsI(siteIdsI[sI]);
  62
  63         forAll(siteIdsJ, sJ)
  64         {
  65             label idsJ(siteIdsJ[sJ]);
  66
  67             if (pairPotentialSitesI[sI] && pairPotentialSitesJ[sJ])
  68             {
  69                 vector rsIsJ =
  70                     molI->sitePositions()[sI] - molJ->sitePositions()[sJ];
  71
  72                 scalar rsIsJMagSq = magSqr(rsIsJ);
  73
  74                 if(pairPot.rCutSqr(idsI, idsJ, rsIsJMagSq))
  75                 {
  76                     scalar rsIsJMag = mag(rsIsJ);
  77
  78                     vector fsIsJ =
  79                         (rsIsJ/rsIsJMag)
  80                        *pairPot.force(idsI, idsJ, rsIsJMag);
  81
  82                     molI->siteForces()[sI] += fsIsJ;
  83
  84                     molJ->siteForces()[sJ] += -fsIsJ;
  85
  86                     scalar potentialEnergy
  87                     (
  88                         pairPot.energy(idsI, idsJ, rsIsJMag)
  89                     );
  90
  91                     molI->potentialEnergy() += 0.5*potentialEnergy;
  92
  93                     molJ->potentialEnergy() += 0.5*potentialEnergy;
  94
  95                     vector rIJ = molI->position() - molJ->position();
  96
  97                     tensor virialContribution =
  98                         (rsIsJ*fsIsJ)*(rsIsJ & rIJ)/rsIsJMagSq;
  99
 100                     molI->rf() += virialContribution;
 101
 102                     molJ->rf() += virialContribution;
 103
 104                     // molI->rf() += rsIsJ * fsIsJ;
 105
 106                     // molJ->rf() += rsIsJ * fsIsJ;
 107                 }
 108             }
 109
 110             if (electrostaticSitesI[sI] && electrostaticSitesJ[sJ])
 111             {
 112                 vector rsIsJ =
 113                 molI->sitePositions()[sI] - molJ->sitePositions()[sJ];
 114
 115                 scalar rsIsJMagSq = magSqr(rsIsJ);
 116
 117                 if(rsIsJMagSq <= electrostatic.rCutSqr())
 118                 {
 119                     scalar rsIsJMag = mag(rsIsJ);
 120
 121                     scalar chargeI = constPropI.siteCharges()[sI];
 122
 123                     scalar chargeJ = constPropJ.siteCharges()[sJ];
 124
 125                     vector fsIsJ =
 126                         (rsIsJ/rsIsJMag)
 127                        *chargeI*chargeJ*electrostatic.force(rsIsJMag);
 128
 129                     molI->siteForces()[sI] += fsIsJ;
 130
 131                     molJ->siteForces()[sJ] += -fsIsJ;
 132
 133                     scalar potentialEnergy =
 134                         chargeI*chargeJ
 135                        *electrostatic.energy(rsIsJMag);
 136
 137                     molI->potentialEnergy() += 0.5*potentialEnergy;
 138
 139                     molJ->potentialEnergy() += 0.5*potentialEnergy;
 140
 141                     vector rIJ = molI->position() - molJ->position();
 142
 143                     tensor virialContribution =
 144                         (rsIsJ*fsIsJ)*(rsIsJ & rIJ)/rsIsJMagSq;
 145
 146                     molI->rf() += virialContribution;
 147
 148                     molJ->rf() += virialContribution;
 149
 150                     // molI->rf() += rsIsJ * fsIsJ;
 151
 152                     // molJ->rf() += rsIsJ * fsIsJ;
 153                 }
 154             }
 155         }
 156     }
 157 }
 158
 159
 160 inline void Foam::moleculeCloud::evaluatePair
 161 (
 162     molecule* molReal,
 163     referredMolecule* molRef
 164 )
 165 {
 166     const pairPotentialList& pairPot = pot_.pairPotentials();
 167
 168     const pairPotential& electrostatic = pairPot.electrostatic();
 169
 170     label idReal = molReal->id();
 171
 172     label idRef = molRef->id();
 173
 174     const molecule::constantProperties& constPropReal(constProps(idReal));
 175
 176     const molecule::constantProperties& constPropRef(constProps(idRef));
 177
 178     List<label> siteIdsReal = constPropReal.siteIds();
 179
 180     List<label> siteIdsRef = constPropRef.siteIds();
 181
 182     List<bool> pairPotentialSitesReal = constPropReal.pairPotentialSites();
 183
 184     List<bool> electrostaticSitesReal = constPropReal.electrostaticSites();
 185
 186     List<bool> pairPotentialSitesRef = constPropRef.pairPotentialSites();
 187
 188     List<bool> electrostaticSitesRef = constPropRef.electrostaticSites();
 189
 190     forAll(siteIdsReal, sReal)
 191     {
 192         label idsReal(siteIdsReal[sReal]);
 193
 194         forAll(siteIdsRef, sRef)
 195         {
 196             label idsRef(siteIdsRef[sRef]);
 197
 198             if (pairPotentialSitesReal[sReal] && pairPotentialSitesRef[sRef])
 199             {
 200                 vector rsRealsRef =
 201                     molReal->sitePositions()[sReal]
 202                   - molRef->sitePositions()[sRef];
 203
 204                 scalar rsRealsRefMagSq = magSqr(rsRealsRef);
 205
 206                 if (pairPot.rCutSqr(idsReal, idsRef, rsRealsRefMagSq))
 207                 {
 208                     scalar rsRealsRefMag = mag(rsRealsRef);
 209
 210                     vector fsRealsRef =
 211                         (rsRealsRef/rsRealsRefMag)
 212                        *pairPot.force(idsReal, idsRef, rsRealsRefMag);
 213
 214                     molReal->siteForces()[sReal] += fsRealsRef;
 215
 216                     scalar potentialEnergy
 217                     (
 218                         pairPot.energy(idsReal, idsRef, rsRealsRefMag)
 219                     );
 220
 221                     molReal->potentialEnergy() += 0.5*potentialEnergy;
 222
 223                     vector rRealRef = molReal->position() - molRef->position();
 224
 225                     molReal->rf() +=
 226                         (rsRealsRef*fsRealsRef)
 227                        *(rsRealsRef & rRealRef)
 228                        /rsRealsRefMagSq;
 229
 230                     // molReal->rf() += rsRealsRef * fsRealsRef;
 231
 232                 }
 233             }
 234
 235             if (electrostaticSitesReal[sReal] && electrostaticSitesRef[sRef])
 236             {
 237                 vector rsRealsRef =
 238                     molReal->sitePositions()[sReal]
 239                   - molRef->sitePositions()[sRef];
 240
 241                 scalar rsRealsRefMagSq = magSqr(rsRealsRef);
 242
 243                 if (rsRealsRefMagSq <= electrostatic.rCutSqr())
 244                 {
 245                     scalar rsRealsRefMag = mag(rsRealsRef);
 246
 247                     scalar chargeReal = constPropReal.siteCharges()[sReal];
 248
 249                     scalar chargeRef = constPropRef.siteCharges()[sRef];
 250
 251                     vector fsRealsRef =
 252                         (rsRealsRef/rsRealsRefMag)
 253                        *chargeReal*chargeRef
 254                        *electrostatic.force(rsRealsRefMag);
 255
 256                     molReal->siteForces()[sReal] += fsRealsRef;
 257
 258                     scalar potentialEnergy =
 259                         chargeReal*chargeRef
 260                        *electrostatic.energy(rsRealsRefMag);
 261
 262                     molReal->potentialEnergy() += 0.5*potentialEnergy;
 263
 264                     vector rRealRef = molReal->position() - molRef->position();
 265
 266                     molReal->rf() +=
 267                         (rsRealsRef*fsRealsRef)
 268                        *(rsRealsRef & rRealRef)
 269                        /rsRealsRefMagSq;
 270
 271                     // molReal->rf() += rsRealsRef * fsRealsRef;
 272                 }
 273             }
 274         }
 275     }
 276 }
 277
 278
 279 inline bool Foam::moleculeCloud::evaluatePotentialLimit
 280 (
 281     molecule* molI,
 282     molecule* molJ
 283 ) const
 284 {
 285     const pairPotentialList& pairPot = pot_.pairPotentials();
 286
 287     const pairPotential& electrostatic = pairPot.electrostatic();
 288
 289     label idI = molI->id();
 290
 291     label idJ = molJ->id();
 292
 293     const molecule::constantProperties& constPropI(constProps(idI));
 294
 295     const molecule::constantProperties& constPropJ(constProps(idJ));
 296
 297     List<label> siteIdsI = constPropI.siteIds();
 298
 299     List<label> siteIdsJ = constPropJ.siteIds();
 300
 301     List<bool> pairPotentialSitesI = constPropI.pairPotentialSites();
 302
 303     List<bool> electrostaticSitesI = constPropI.electrostaticSites();
 304
 305     List<bool> pairPotentialSitesJ = constPropJ.pairPotentialSites();
 306
 307     List<bool> electrostaticSitesJ = constPropJ.electrostaticSites();
 308
 309     forAll(siteIdsI, sI)
 310     {
 311         label idsI(siteIdsI[sI]);
 312
 313         forAll(siteIdsJ, sJ)
 314         {
 315             label idsJ(siteIdsJ[sJ]);
 316
 317             if (pairPotentialSitesI[sI] && pairPotentialSitesJ[sJ])
 318             {
 319                 vector rsIsJ =
 320                     molI->sitePositions()[sI] - molJ->sitePositions()[sJ];
 321
 322                 scalar rsIsJMagSq = magSqr(rsIsJ);
 323
 324                 if (pairPot.rCutSqr(idsI, idsJ, rsIsJMagSq))
 325                 {
 326                     scalar rsIsJMag = mag(rsIsJ);
 327
 328                     // Guard against pairPot.energy being evaluated
 329                     // if rIJMag < SMALL. A floating point exception will
 330                     // happen otherwise.
 331
 332                     if (rsIsJMag < SMALL)
 333                     {
 334                         WarningIn("moleculeCloud::removeHighEnergyOverlaps()")
 335                             << "Molecule site pair closer than "
 336                             << SMALL
 337                             << ": mag separation = " << rsIsJMag
 338                             << ". These may have been placed on top of each"
 339                             << " other by a rounding error in mdInitialise in"
 340                             << " parallel or a block filled with molecules"
 341                             << " twice. Removing one of the molecules."
 342                             << endl;
 343
 344                         return true;
 345                     }
 346
 347                     // Guard against pairPot.energy being evaluated if rIJMag <
 348                     // rMin.  A tabulation lookup error will occur otherwise.
 349
 350                     if (rsIsJMag < pairPot.rMin(idsI, idsJ))
 351                     {
 352                         return true;
 353                     }
 354
 355                     if
 356                     (
 357                         mag(pairPot.energy(idsI, idsJ, rsIsJMag))
 358                       > pot_.potentialEnergyLimit()
 359                     )
 360                     {
 361                         return true;
 362                     };
 363                 }
 364             }
 365
 366             if (electrostaticSitesI[sI] && electrostaticSitesJ[sJ])
 367             {
 368                 vector rsIsJ =
 369                     molI->sitePositions()[sI] - molJ->sitePositions()[sJ];
 370
 371                 scalar rsIsJMagSq = magSqr(rsIsJ);
 372
 373                 if (pairPot.rCutMaxSqr(rsIsJMagSq))
 374                 {
 375                     scalar rsIsJMag = mag(rsIsJ);
 376
 377                     // Guard against pairPot.energy being evaluated
 378                     // if rIJMag < SMALL. A floating point exception will
 379                     // happen otherwise.
 380
 381                     if (rsIsJMag < SMALL)
 382                     {
 383                         WarningIn("moleculeCloud::removeHighEnergyOverlaps()")
 384                             << "Molecule site pair closer than "
 385                             << SMALL
 386                             << ": mag separation = " << rsIsJMag
 387                             << ". These may have been placed on top of each"
 388                             << " other by a rounding error in mdInitialise in"
 389                             << " parallel or a block filled with molecules"
 390                             << " twice. Removing one of the molecules."
 391                             << endl;
 392
 393                         return true;
 394                     }
 395
 396                     if (rsIsJMag < electrostatic.rMin())
 397                     {
 398                         return true;
 399                     }
 400
 401                     scalar chargeI = constPropI.siteCharges()[sI];
 402
 403                     scalar chargeJ = constPropJ.siteCharges()[sJ];
 404
 405                     if
 406                     (
 407                         mag(chargeI*chargeJ*electrostatic.energy(rsIsJMag))
 408                       > pot_.potentialEnergyLimit()
 409                     )
 410                     {
 411                         return true;
 412                     };
 413                 }
 414             }
 415         }
 416     }
 417
 418     return false;
 419 }
 420
 421
 422 inline bool Foam::moleculeCloud::evaluatePotentialLimit
 423 (
 424     molecule* molReal,
 425     referredMolecule* molRef
 426 ) const
 427 {
 428     const pairPotentialList& pairPot = pot_.pairPotentials();
 429
 430     const pairPotential& electrostatic = pairPot.electrostatic();
 431
 432     label idReal = molReal->id();
 433
 434     label idRef = molRef->id();
 435
 436     const molecule::constantProperties& constPropReal(constProps(idReal));
 437
 438     const molecule::constantProperties& constPropRef(constProps(idRef));
 439
 440     List<label> siteIdsReal = constPropReal.siteIds();
 441
 442     List<label> siteIdsRef = constPropRef.siteIds();
 443
 444     List<bool> pairPotentialSitesReal = constPropReal.pairPotentialSites();
 445
 446     List<bool> electrostaticSitesReal = constPropReal.electrostaticSites();
 447
 448     List<bool> pairPotentialSitesRef = constPropRef.pairPotentialSites();
 449
 450     List<bool> electrostaticSitesRef = constPropRef.electrostaticSites();
 451
 452     forAll(siteIdsReal, sReal)
 453     {
 454         label idsReal(siteIdsReal[sReal]);
 455
 456         forAll(siteIdsRef, sRef)
 457         {
 458             label idsRef(siteIdsRef[sRef]);
 459
 460             if (pairPotentialSitesReal[sReal] && pairPotentialSitesRef[sRef])
 461             {
 462                 vector rsRealsRef =
 463                     molReal->sitePositions()[sReal]
 464                   - molRef->sitePositions()[sRef];
 465
 466                 scalar rsRealsRefMagSq = magSqr(rsRealsRef);
 467
 468                 if (pairPot.rCutSqr(idsReal, idsRef, rsRealsRefMagSq))
 469                 {
 470                     scalar rsRealsRefMag = mag(rsRealsRef);
 471
 472                     // Guard against pairPot.energy being evaluated
 473                     // if rRealRefMag < SMALL. A floating point exception will
 474                     // happen otherwise.
 475
 476                     if (rsRealsRefMag < SMALL)
 477                     {
 478                         WarningIn("moleculeCloud::removeHighEnergyOverlaps()")
 479                             << "Molecule site pair closer than "
 480                             << SMALL
 481                             << ": mag separation = " << rsRealsRefMag
 482                             << ". These may have been placed on top of each"
 483                             << " other by a rounding error in mdInitialise in"
 484                             << " parallel or a block filled with molecules"
 485                             << " twice. Removing one of the molecules."
 486                             << endl;
 487
 488                         return true;
 489                     }
 490
 491                     // Guard against pairPot.energy being evaluated if
 492                     // rRealRefMag < rMin.  A tabulation lookup error will occur
 493                     // otherwise.
 494
 495                     if (rsRealsRefMag < pairPot.rMin(idsReal, idsRef))
 496                     {
 497                         return true;
 498                     }
 499
 500                     if
 501                     (
 502                         mag(pairPot.energy(idsReal, idsRef, rsRealsRefMag))
 503                       > pot_.potentialEnergyLimit()
 504                     )
 505                     {
 506                         return true;
 507                     };
 508                 }
 509             }
 510
 511             if (electrostaticSitesReal[sReal] && electrostaticSitesRef[sRef])
 512             {
 513                 vector rsRealsRef =
 514                     molReal->sitePositions()[sReal]
 515                   - molRef->sitePositions()[sRef];
 516
 517                 scalar rsRealsRefMagSq = magSqr(rsRealsRef);
 518
 519                 if (pairPot.rCutMaxSqr(rsRealsRefMagSq))
 520                 {
 521                     scalar rsRealsRefMag = mag(rsRealsRef);
 522
 523                     // Guard against pairPot.energy being evaluated
 524                     // if rRealRefMag < SMALL. A floating point exception will
 525                     // happen otherwise.
 526
 527                     if (rsRealsRefMag < SMALL)
 528                     {
 529                         WarningIn("moleculeCloud::removeHighEnergyOverlaps()")
 530                             << "Molecule site pair closer than "
 531                             << SMALL
 532                             << ": mag separation = " << rsRealsRefMag
 533                             << ". These may have been placed on top of each"
 534                             << " other by a rounding error in mdInitialise in"
 535                             << " parallel or a block filled with molecules"
 536                             << " twice. Removing one of the molecules."
 537                             << endl;
 538
 539                         return true;
 540                     }
 541
 542                     if (rsRealsRefMag < electrostatic.rMin())
 543                     {
 544                         return true;
 545                     }
 546
 547                     scalar chargeReal = constPropReal.siteCharges()[sReal];
 548
 549                     scalar chargeRef = constPropRef.siteCharges()[sRef];
 550
 551                     if
 552                     (
 553                         mag
 554                         (
 555                             chargeReal
 556                            *chargeRef
 557                            *electrostatic.energy(rsRealsRefMag)
 558                         )
 559                       > pot_.potentialEnergyLimit()
 560                     )
 561                     {
 562                         return true;
 563                     };
 564                 }
 565             }
 566         }
 567     }
 568
 569     return false;
 570 }
 571
 572
 573 inline Foam::vector Foam::moleculeCloud::equipartitionLinearVelocity
 574 (
 575     scalar temperature,
 576     scalar mass
 577 )
 578 {
 579     return sqrt(kb*temperature/mass)*vector
 580     (
 581         rndGen_.GaussNormal(),
 582         rndGen_.GaussNormal(),
 583         rndGen_.GaussNormal()
 584     );
 585 }
 586
 587
 588 inline Foam::vector Foam::moleculeCloud::equipartitionAngularMomentum
 589 (
 590     scalar temperature,
 591     const molecule::constantProperties& cP
 592 )
 593 {
 594     scalar sqrtKbT = sqrt(kb*temperature);
 595
 596     if (cP.linearMolecule())
 597     {
 598         return sqrtKbT*vector
 599         (
 600             0.0,
 601             sqrt(cP.momentOfInertia().yy())*rndGen_.GaussNormal(),
 602             sqrt(cP.momentOfInertia().zz())*rndGen_.GaussNormal()
 603         );
 604     }
 605     else
 606     {
 607         return sqrtKbT*vector
 608         (
 609             sqrt(cP.momentOfInertia().xx())*rndGen_.GaussNormal(),
 610             sqrt(cP.momentOfInertia().yy())*rndGen_.GaussNormal(),
 611             sqrt(cP.momentOfInertia().zz())*rndGen_.GaussNormal()
 612         );
 613     }
 614 }
 615
 616
 617 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 618
 619 inline const Foam::polyMesh& Foam::moleculeCloud::mesh() const
 620 {
 621     return mesh_;
 622 }
 623
 624
 625 inline const Foam::potential& Foam::moleculeCloud::pot() const
 626 {
 627     return pot_;
 628 }
 629
 630
 631 inline const Foam::List<Foam::DynamicList<Foam::molecule*> >&
 632     Foam::moleculeCloud::cellOccupancy() const
 633 {
 634     return cellOccupancy_;
 635 }
 636
 637
 638 inline const Foam::interactionLists&
 639     Foam::moleculeCloud::il() const
 640 {
 641     return il_;
 642 }
 643
 644
 645 inline const Foam::List<Foam::molecule::constantProperties>
 646     Foam::moleculeCloud::constProps() const
 647 {
 648     return constPropList_;
 649 }
 650
 651
 652 inline const Foam::molecule::constantProperties&
 653     Foam::moleculeCloud::constProps(label id) const
 654 {
 655     return constPropList_[id];
 656 }
 657
 658
 659 inline Foam::Random& Foam::moleculeCloud::rndGen()
 660 {
 661     return rndGen_;
 662 }
 663
 664
 665 // ************************************************************************* //