src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloud.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 "moleculeCloud.H"
  27 #include "fvMesh.H"
  28 #include "mathematicalConstants.H"
  29
  30 using namespace Foam::constant::mathematical;
  31
  32 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  33
  34 namespace Foam
  35 {
  36     defineTemplateTypeNameAndDebug(Cloud<molecule>, 0);
  37 }
  38
  39 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  40
  41 void Foam::moleculeCloud::buildConstProps()
  42 {
  43     Info<< nl << "Reading moleculeProperties dictionary." << endl;
  44
  45     const List<word>& idList(pot_.idList());
  46
  47     constPropList_.setSize(idList.size());
  48
  49     const List<word>& siteIdList(pot_.siteIdList());
  50
  51     IOdictionary moleculePropertiesDict
  52     (
  53         IOobject
  54         (
  55             "moleculeProperties",
  56             mesh_.time().constant(),
  57             mesh_,
  58             IOobject::MUST_READ_IF_MODIFIED,
  59             IOobject::NO_WRITE,
  60             false
  61         )
  62     );
  63
  64     forAll(idList, i)
  65     {
  66         const word& id = idList[i];
  67         const dictionary& molDict = moleculePropertiesDict.subDict(id);
  68
  69         List<word> siteIdNames = molDict.lookup("siteIds");
  70
  71         List<label> siteIds(siteIdNames.size());
  72
  73         forAll(siteIdNames, sI)
  74         {
  75             const word& siteId = siteIdNames[sI];
  76
  77             siteIds[sI] = findIndex(siteIdList, siteId);
  78
  79             if (siteIds[sI] == -1)
  80             {
  81                 FatalErrorIn("moleculeCloud::buildConstProps()")
  82                     << siteId << " site not found."
  83                     << nl << abort(FatalError);
  84             }
  85         }
  86
  87         molecule::constantProperties& constProp = constPropList_[i];
  88
  89         constProp = molecule::constantProperties(molDict);
  90
  91         constProp.siteIds() = siteIds;
  92     }
  93 }
  94
  95
  96 void Foam::moleculeCloud::setSiteSizesAndPositions()
  97 {
  98     forAllIter(moleculeCloud, *this, mol)
  99     {
 100         const molecule::constantProperties& cP = constProps(mol().id());
 101
 102         mol().setSiteSizes(cP.nSites());
 103
 104         mol().setSitePositions(cP);
 105     }
 106 }
 107
 108
 109 void Foam::moleculeCloud::buildCellOccupancy()
 110 {
 111     forAll(cellOccupancy_, cO)
 112     {
 113         cellOccupancy_[cO].clear();
 114     }
 115
 116     forAllIter(moleculeCloud, *this, mol)
 117     {
 118         cellOccupancy_[mol().cell()].append(&mol());
 119     }
 120
 121     forAll(cellOccupancy_, cO)
 122     {
 123         cellOccupancy_[cO].shrink();
 124     }
 125 }
 126
 127
 128 void Foam::moleculeCloud::calculatePairForce()
 129 {
 130     PstreamBuffers pBufs(Pstream::nonBlocking);
 131
 132     // Start sending referred data
 133     il_.sendReferredData(cellOccupancy(), pBufs);
 134
 135     molecule* molI = NULL;
 136     molecule* molJ = NULL;
 137
 138     {
 139         // Real-Real interactions
 140
 141         const labelListList& dil = il_.dil();
 142
 143         forAll(dil, d)
 144         {
 145             forAll(cellOccupancy_[d],cellIMols)
 146             {
 147                 molI = cellOccupancy_[d][cellIMols];
 148
 149                 forAll(dil[d], interactingCells)
 150                 {
 151                     List<molecule*> cellJ =
 152                         cellOccupancy_[dil[d][interactingCells]];
 153
 154                     forAll(cellJ, cellJMols)
 155                     {
 156                         molJ = cellJ[cellJMols];
 157
 158                         evaluatePair(*molI, *molJ);
 159                     }
 160                 }
 161
 162                 forAll(cellOccupancy_[d], cellIOtherMols)
 163                 {
 164                     molJ = cellOccupancy_[d][cellIOtherMols];
 165
 166                     if (molJ > molI)
 167                     {
 168                         evaluatePair(*molI, *molJ);
 169                     }
 170                 }
 171             }
 172         }
 173     }
 174
 175     // Receive referred data
 176     il_.receiveReferredData(pBufs);
 177
 178     {
 179         // Real-Referred interactions
 180
 181         const labelListList& ril = il_.ril();
 182
 183         List<IDLList<molecule> >& referredMols = il_.referredParticles();
 184
 185         forAll(ril, r)
 186         {
 187             const List<label>& realCells = ril[r];
 188
 189             IDLList<molecule>& refMols = referredMols[r];
 190
 191             forAllIter
 192             (
 193                 IDLList<molecule>,
 194                 refMols,
 195                 refMol
 196             )
 197             {
 198                 forAll(realCells, rC)
 199                 {
 200                     List<molecule*> cellI = cellOccupancy_[realCells[rC]];
 201
 202                     forAll(cellI, cellIMols)
 203                     {
 204                         molI = cellI[cellIMols];
 205
 206                         evaluatePair(*molI, refMol());
 207                     }
 208                 }
 209             }
 210         }
 211     }
 212 }
 213
 214
 215 void Foam::moleculeCloud::calculateTetherForce()
 216 {
 217     const tetherPotentialList& tetherPot(pot_.tetherPotentials());
 218
 219     forAllIter(moleculeCloud, *this, mol)
 220     {
 221         if (mol().tethered())
 222         {
 223             vector rIT = mol().position() - mol().specialPosition();
 224
 225             label idI = mol().id();
 226
 227             scalar massI = constProps(idI).mass();
 228
 229             vector fIT = tetherPot.force(idI, rIT);
 230
 231             mol().a() += fIT/massI;
 232
 233             mol().potentialEnergy() += tetherPot.energy(idI, rIT);
 234
 235             // What to do here?
 236             // mol().rf() += rIT*fIT;
 237         }
 238     }
 239 }
 240
 241
 242 void Foam::moleculeCloud::calculateExternalForce()
 243 {
 244     forAllIter(moleculeCloud, *this, mol)
 245     {
 246         mol().a() += pot_.gravity();
 247     }
 248 }
 249
 250
 251 void Foam::moleculeCloud::removeHighEnergyOverlaps()
 252 {
 253     Info<< nl << "Removing high energy overlaps, limit = "
 254         << pot_.potentialEnergyLimit()
 255         << nl << "Removal order:";
 256
 257     forAll(pot_.removalOrder(), rO)
 258     {
 259         Info<< ' ' << pot_.idList()[pot_.removalOrder()[rO]];
 260     }
 261
 262     Info<< nl ;
 263
 264     label initialSize = this->size();
 265
 266     buildCellOccupancy();
 267
 268     // Real-Real interaction
 269
 270     molecule* molI = NULL;
 271     molecule* molJ = NULL;
 272
 273     {
 274         DynamicList<molecule*> molsToDelete;
 275
 276         const labelListList& dil(il_.dil());
 277
 278         forAll(dil, d)
 279         {
 280             forAll(cellOccupancy_[d],cellIMols)
 281             {
 282                 molI = cellOccupancy_[d][cellIMols];
 283
 284                 forAll(dil[d], interactingCells)
 285                 {
 286                     List<molecule*> cellJ =
 287                         cellOccupancy_[dil[d][interactingCells]];
 288
 289                     forAll(cellJ, cellJMols)
 290                     {
 291                         molJ = cellJ[cellJMols];
 292
 293                         if (evaluatePotentialLimit(*molI, *molJ))
 294                         {
 295                             label idI = molI->id();
 296
 297                             label idJ = molJ->id();
 298
 299                             if
 300                             (
 301                                 idI == idJ
 302                              || findIndex(pot_.removalOrder(), idJ)
 303                               < findIndex(pot_.removalOrder(), idI)
 304                             )
 305                             {
 306                                 if (findIndex(molsToDelete, molJ) == -1)
 307                                 {
 308                                     molsToDelete.append(molJ);
 309                                 }
 310                             }
 311                             else if (findIndex(molsToDelete, molI) == -1)
 312                             {
 313                                 molsToDelete.append(molI);
 314                             }
 315                         }
 316                     }
 317                 }
 318             }
 319
 320             forAll(cellOccupancy_[d], cellIOtherMols)
 321             {
 322                 molJ = cellOccupancy_[d][cellIOtherMols];
 323
 324                 if (molJ > molI)
 325                 {
 326                     if (evaluatePotentialLimit(*molI, *molJ))
 327                     {
 328                         label idI = molI->id();
 329
 330                         label idJ = molJ->id();
 331
 332                         if
 333                         (
 334                             idI == idJ
 335                          || findIndex(pot_.removalOrder(), idJ)
 336                           < findIndex(pot_.removalOrder(), idI)
 337                         )
 338                         {
 339                             if (findIndex(molsToDelete, molJ) == -1)
 340                             {
 341                                 molsToDelete.append(molJ);
 342                             }
 343                         }
 344                         else if (findIndex(molsToDelete, molI) == -1)
 345                         {
 346                             molsToDelete.append(molI);
 347                         }
 348                     }
 349                 }
 350             }
 351         }
 352
 353         forAll(molsToDelete, mTD)
 354         {
 355             deleteParticle(*(molsToDelete[mTD]));
 356         }
 357     }
 358
 359     buildCellOccupancy();
 360
 361     PstreamBuffers pBufs(Pstream::nonBlocking);
 362
 363     // Start sending referred data
 364     il_.sendReferredData(cellOccupancy(), pBufs);
 365
 366         // Receive referred data
 367     il_.receiveReferredData(pBufs);
 368
 369     // Real-Referred interaction
 370
 371     {
 372         DynamicList<molecule*> molsToDelete;
 373
 374         const labelListList& ril(il_.ril());
 375
 376         List<IDLList<molecule> >& referredMols = il_.referredParticles();
 377
 378         forAll(ril, r)
 379         {
 380             IDLList<molecule>& refMols = referredMols[r];
 381
 382             forAllIter
 383             (
 384                 IDLList<molecule>,
 385                 refMols,
 386                 refMol
 387             )
 388             {
 389                 molJ = &refMol();
 390
 391                 const List<label>& realCells = ril[r];
 392
 393                 forAll(realCells, rC)
 394                 {
 395                     label cellI = realCells[rC];
 396
 397                     List<molecule*> cellIMols = cellOccupancy_[cellI];
 398
 399                     forAll(cellIMols, cIM)
 400                     {
 401                         molI = cellIMols[cIM];
 402
 403                         if (evaluatePotentialLimit(*molI, *molJ))
 404                         {
 405                             label idI = molI->id();
 406
 407                             label idJ = molJ->id();
 408
 409                             if
 410                             (
 411                                 findIndex(pot_.removalOrder(), idI)
 412                               < findIndex(pot_.removalOrder(), idJ)
 413                             )
 414                             {
 415                                 if (findIndex(molsToDelete, molI) == -1)
 416                                 {
 417                                     molsToDelete.append(molI);
 418                                 }
 419                             }
 420                             else if
 421                             (
 422                                 findIndex(pot_.removalOrder(), idI)
 423                              == findIndex(pot_.removalOrder(), idJ)
 424                             )
 425                             {
 426                                 // Remove one of the molecules
 427                                 // arbitrarily, assuring that a
 428                                 // consistent decision is made for
 429                                 // both real-referred pairs.
 430
 431                                 if (molI->origId() > molJ->origId())
 432                                 {
 433                                     if (findIndex(molsToDelete, molI) == -1)
 434                                     {
 435                                         molsToDelete.append(molI);
 436                                     }
 437                                 }
 438                             }
 439                         }
 440                     }
 441                 }
 442             }
 443         }
 444
 445         forAll(molsToDelete, mTD)
 446         {
 447             deleteParticle(*(molsToDelete[mTD]));
 448         }
 449     }
 450
 451     buildCellOccupancy();
 452
 453     // Start sending referred data
 454     il_.sendReferredData(cellOccupancy(), pBufs);
 455
 456     // Receive referred data
 457     il_.receiveReferredData(pBufs);
 458
 459     label molsRemoved = initialSize - this->size();
 460
 461     if (Pstream::parRun())
 462     {
 463         reduce(molsRemoved, sumOp<label>());
 464     }
 465
 466     Info<< tab << molsRemoved << " molecules removed" << endl;
 467 }
 468
 469
 470 void Foam::moleculeCloud::initialiseMolecules
 471 (
 472     const IOdictionary& mdInitialiseDict
 473 )
 474 {
 475     Info<< nl
 476         << "Initialising molecules in each zone specified in mdInitialiseDict."
 477         << endl;
 478
 479     const cellZoneMesh& cellZones = mesh_.cellZones();
 480
 481     if (!cellZones.size())
 482     {
 483         FatalErrorIn("void Foam::moleculeCloud::initialiseMolecules")
 484             << "No cellZones found in the mesh."
 485             << abort(FatalError);
 486     }
 487
 488     forAll(cellZones, z)
 489     {
 490         const cellZone& zone(cellZones[z]);
 491
 492         if (zone.size())
 493         {
 494             if (!mdInitialiseDict.found(zone.name()))
 495             {
 496                 Info<< "No specification subDictionary for zone "
 497                     << zone.name() << " found, skipping." << endl;
 498             }
 499             else
 500             {
 501                 const dictionary& zoneDict =
 502                     mdInitialiseDict.subDict(zone.name());
 503
 504                 const scalar temperature
 505                 (
 506                     readScalar(zoneDict.lookup("temperature"))
 507                 );
 508
 509                 const vector bulkVelocity(zoneDict.lookup("bulkVelocity"));
 510
 511                 List<word> latticeIds
 512                 (
 513                     zoneDict.lookup("latticeIds")
 514                 );
 515
 516                 List<vector> latticePositions
 517                 (
 518                     zoneDict.lookup("latticePositions")
 519                 );
 520
 521                 if (latticeIds.size() != latticePositions.size())
 522                 {
 523                     FatalErrorIn("Foam::moleculeCloud::initialiseMolecules")
 524                         << "latticeIds and latticePositions must be the same "
 525                         << " size." << nl
 526                         << abort(FatalError);
 527                 }
 528
 529                 diagTensor latticeCellShape
 530                 (
 531                     zoneDict.lookup("latticeCellShape")
 532                 );
 533
 534                 scalar latticeCellScale = 0.0;
 535
 536                 if (zoneDict.found("numberDensity"))
 537                 {
 538                     scalar numberDensity = readScalar
 539                     (
 540                         zoneDict.lookup("numberDensity")
 541                     );
 542
 543                     if (numberDensity < VSMALL)
 544                     {
 545                         WarningIn("moleculeCloud::initialiseMolecules")
 546                             << "numberDensity too small, not filling zone "
 547                             << zone.name() << endl;
 548
 549                         continue;
 550                     }
 551
 552                     latticeCellScale = pow
 553                     (
 554                         latticeIds.size()/(det(latticeCellShape)*numberDensity),
 555                         (1.0/3.0)
 556                     );
 557                 }
 558                 else if (zoneDict.found("massDensity"))
 559                 {
 560                     scalar unitCellMass = 0.0;
 561
 562                     forAll(latticeIds, i)
 563                     {
 564                         label id = findIndex(pot_.idList(), latticeIds[i]);
 565
 566                         const molecule::constantProperties& cP(constProps(id));
 567
 568                         unitCellMass += cP.mass();
 569                     }
 570
 571                     scalar massDensity = readScalar
 572                     (
 573                         zoneDict.lookup("massDensity")
 574                     );
 575
 576                     if (massDensity < VSMALL)
 577                     {
 578                         WarningIn("moleculeCloud::initialiseMolecules")
 579                             << "massDensity too small, not filling zone "
 580                             << zone.name() << endl;
 581
 582                         continue;
 583                     }
 584
 585
 586                     latticeCellScale = pow
 587                     (
 588                         unitCellMass/(det(latticeCellShape)*massDensity),
 589                         (1.0/3.0)
 590                     );
 591                 }
 592                 else
 593                 {
 594                     FatalErrorIn("Foam::moleculeCloud::initialiseMolecules")
 595                         << "massDensity or numberDensity not specified " << nl
 596                         << abort(FatalError);
 597                 }
 598
 599                 latticeCellShape *= latticeCellScale;
 600
 601                 vector anchor(zoneDict.lookup("anchor"));
 602
 603                 bool tethered = false;
 604
 605                 if (zoneDict.found("tetherSiteIds"))
 606                 {
 607                     tethered = bool
 608                     (
 609                         List<word>(zoneDict.lookup("tetherSiteIds")).size()
 610                     );
 611                 }
 612
 613                 const vector orientationAngles
 614                 (
 615                     zoneDict.lookup("orientationAngles")
 616                 );
 617
 618                 scalar phi(orientationAngles.x()*pi/180.0);
 619
 620                 scalar theta(orientationAngles.y()*pi/180.0);
 621
 622                 scalar psi(orientationAngles.z()*pi/180.0);
 623
 624                 const tensor R
 625                 (
 626                     cos(psi)*cos(phi) - cos(theta)*sin(phi)*sin(psi),
 627                     cos(psi)*sin(phi) + cos(theta)*cos(phi)*sin(psi),
 628                     sin(psi)*sin(theta),
 629                     - sin(psi)*cos(phi) - cos(theta)*sin(phi)*cos(psi),
 630                     - sin(psi)*sin(phi) + cos(theta)*cos(phi)*cos(psi),
 631                     cos(psi)*sin(theta),
 632                     sin(theta)*sin(phi),
 633                     - sin(theta)*cos(phi),
 634                     cos(theta)
 635                 );
 636
 637                 // Find the optimal anchor position.  Finding the approximate
 638                 // mid-point of the zone of cells and snapping to the nearest
 639                 // lattice location.
 640
 641                 vector zoneMin = VGREAT*vector::one;
 642
 643                 vector zoneMax = -VGREAT*vector::one;
 644
 645                 forAll(zone, cell)
 646                 {
 647                     const point cellCentre = mesh_.cellCentres()[zone[cell]];
 648
 649                     if (cellCentre.x() > zoneMax.x())
 650                     {
 651                         zoneMax.x() = cellCentre.x();
 652                     }
 653                     if (cellCentre.x() < zoneMin.x())
 654                     {
 655                         zoneMin.x() = cellCentre.x();
 656                     }
 657                     if (cellCentre.y() > zoneMax.y())
 658                     {
 659                         zoneMax.y() = cellCentre.y();
 660                     }
 661                     if (cellCentre.y() < zoneMin.y())
 662                     {
 663                         zoneMin.y() = cellCentre.y();
 664                     }
 665                     if (cellCentre.z() > zoneMax.z())
 666                     {
 667                         zoneMax.z() = cellCentre.z();
 668                     }
 669                     if (cellCentre.z() < zoneMin.z())
 670                     {
 671                         zoneMin.z() = cellCentre.z();
 672                     }
 673                 }
 674
 675                 point zoneMid = 0.5*(zoneMin + zoneMax);
 676
 677                 point latticeMid = inv(latticeCellShape) & (R.T()
 678                 & (zoneMid - anchor));
 679
 680                 point latticeAnchor
 681                 (
 682                     label(latticeMid.x() + 0.5*sign(latticeMid.x())),
 683                     label(latticeMid.y() + 0.5*sign(latticeMid.y())),
 684                     label(latticeMid.z() + 0.5*sign(latticeMid.z()))
 685                 );
 686
 687                 anchor += (R & (latticeCellShape & latticeAnchor));
 688
 689                 // Continue trying to place molecules as long as at
 690                 // least one molecule is placed in each iteration.
 691                 // The "|| totalZoneMols == 0" condition means that the
 692                 // algorithm will continue if the origin is outside the
 693                 // zone.
 694
 695                 label n = 0;
 696
 697                 label totalZoneMols = 0;
 698
 699                 label molsPlacedThisIteration = 0;
 700
 701                 while
 702                 (
 703                     molsPlacedThisIteration != 0
 704                  || totalZoneMols == 0
 705                 )
 706                 {
 707                     label sizeBeforeIteration = this->size();
 708
 709                     bool partOfLayerInBounds = false;
 710
 711                     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 712                     // start of placement of molecules
 713                     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 714
 715                     if (n == 0)
 716                     {
 717                         // Special treatment is required for the first position,
 718                         // i.e. iteration zero.
 719
 720                         labelVector unitCellLatticePosition(0,0,0);
 721
 722                         forAll(latticePositions, p)
 723                         {
 724                             label id = findIndex(pot_.idList(), latticeIds[p]);
 725
 726                             const vector& latticePosition =
 727                                 vector
 728                                 (
 729                                     unitCellLatticePosition.x(),
 730                                     unitCellLatticePosition.y(),
 731                                     unitCellLatticePosition.z()
 732                                 )
 733                               + latticePositions[p];
 734
 735                             point globalPosition =
 736                                 anchor
 737                               + (R & (latticeCellShape & latticePosition));
 738
 739                             partOfLayerInBounds = mesh_.bounds().contains
 740                             (
 741                                 globalPosition
 742                             );
 743
 744                             label cell = -1;
 745                             label tetFace = -1;
 746                             label tetPt = -1;
 747
 748                             mesh_.findCellFacePt
 749                             (
 750                                 globalPosition,
 751                                 cell,
 752                                 tetFace,
 753                                 tetPt
 754                             );
 755
 756                             if (findIndex(zone, cell) != -1)
 757                             {
 758                                 createMolecule
 759                                 (
 760                                     globalPosition,
 761                                     cell,
 762                                     tetFace,
 763                                     tetPt,
 764                                     id,
 765                                     tethered,
 766                                     temperature,
 767                                     bulkVelocity
 768                                 );
 769                             }
 770                         }
 771                     }
 772                     else
 773                     {
 774                         // Place top and bottom caps.
 775
 776                         labelVector unitCellLatticePosition(0,0,0);
 777
 778                         for
 779                         (
 780                             unitCellLatticePosition.z() = -n;
 781                             unitCellLatticePosition.z() <= n;
 782                             unitCellLatticePosition.z() += 2*n
 783                         )
 784                         {
 785                             for
 786                             (
 787                                 unitCellLatticePosition.y() = -n;
 788                                 unitCellLatticePosition.y() <= n;
 789                                 unitCellLatticePosition.y()++
 790                             )
 791                             {
 792                                 for
 793                                 (
 794                                     unitCellLatticePosition.x() = -n;
 795                                     unitCellLatticePosition.x() <= n;
 796                                     unitCellLatticePosition.x()++
 797                                 )
 798                                 {
 799                                     forAll(latticePositions, p)
 800                                     {
 801                                         label id = findIndex
 802                                         (
 803                                             pot_.idList(),
 804                                             latticeIds[p]
 805                                         );
 806
 807                                         const vector& latticePosition =
 808                                             vector
 809                                             (
 810                                                 unitCellLatticePosition.x(),
 811                                                 unitCellLatticePosition.y(),
 812                                                 unitCellLatticePosition.z()
 813                                             )
 814                                           + latticePositions[p];
 815
 816                                         point globalPosition =
 817                                             anchor
 818                                           + (
 819                                                 R
 820                                               & (
 821                                                     latticeCellShape
 822                                                   & latticePosition
 823                                                 )
 824                                             );
 825
 826                                         partOfLayerInBounds =
 827                                             mesh_.bounds().contains
 828                                             (
 829                                                 globalPosition
 830                                             );
 831
 832                                         label cell = -1;
 833                                         label tetFace = -1;
 834                                         label tetPt = -1;
 835
 836                                         mesh_.findCellFacePt
 837                                         (
 838                                             globalPosition,
 839                                             cell,
 840                                             tetFace,
 841                                             tetPt
 842                                         );
 843
 844                                         if (findIndex(zone, cell) != -1)
 845                                         {
 846                                             createMolecule
 847                                             (
 848                                                 globalPosition,
 849                                                 cell,
 850                                                 tetFace,
 851                                                 tetPt,
 852                                                 id,
 853                                                 tethered,
 854                                                 temperature,
 855                                                 bulkVelocity
 856                                             );
 857                                         }
 858                                     }
 859                                 }
 860                             }
 861                         }
 862
 863                         for
 864                         (
 865                             unitCellLatticePosition.z() = -(n-1);
 866                             unitCellLatticePosition.z() <= (n-1);
 867                             unitCellLatticePosition.z()++
 868                         )
 869                         {
 870                             for (label iR = 0; iR <= 2*n -1; iR++)
 871                             {
 872                                 unitCellLatticePosition.x() = n;
 873
 874                                 unitCellLatticePosition.y() = -n + (iR + 1);
 875
 876                                 for (label iK = 0; iK < 4; iK++)
 877                                 {
 878                                     forAll(latticePositions, p)
 879                                     {
 880                                         label id = findIndex
 881                                         (
 882                                             pot_.idList(),
 883                                             latticeIds[p]
 884                                         );
 885
 886                                         const vector& latticePosition =
 887                                             vector
 888                                             (
 889                                                 unitCellLatticePosition.x(),
 890                                                 unitCellLatticePosition.y(),
 891                                                 unitCellLatticePosition.z()
 892                                             )
 893                                           + latticePositions[p];
 894
 895                                         point globalPosition =
 896                                             anchor
 897                                           + (
 898                                                 R
 899                                               & (
 900                                                    latticeCellShape
 901                                                  & latticePosition
 902                                                 )
 903                                             );
 904
 905                                         partOfLayerInBounds =
 906                                             mesh_.bounds().contains
 907                                             (
 908                                                 globalPosition
 909                                             );
 910
 911                                         label cell = -1;
 912                                         label tetFace = -1;
 913                                         label tetPt = -1;
 914
 915                                         mesh_.findCellFacePt
 916                                         (
 917                                             globalPosition,
 918                                             cell,
 919                                             tetFace,
 920                                             tetPt
 921                                         );
 922
 923                                         if (findIndex(zone, cell) != -1)
 924                                         {
 925                                             createMolecule
 926                                             (
 927                                                 globalPosition,
 928                                                 cell,
 929                                                 tetFace,
 930                                                 tetPt,
 931                                                 id,
 932                                                 tethered,
 933                                                 temperature,
 934                                                 bulkVelocity
 935                                             );
 936                                         }
 937                                     }
 938
 939                                     unitCellLatticePosition =
 940                                         labelVector
 941                                         (
 942                                           - unitCellLatticePosition.y(),
 943                                             unitCellLatticePosition.x(),
 944                                             unitCellLatticePosition.z()
 945                                         );
 946                                 }
 947                             }
 948                         }
 949                     }
 950
 951                     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 952                     // end of placement of molecules
 953                     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 954
 955                     if
 956                     (
 957                         totalZoneMols == 0
 958                      && !partOfLayerInBounds
 959                     )
 960                     {
 961                         WarningIn("Foam::moleculeCloud::initialiseMolecules()")
 962                             << "A whole layer of unit cells was placed "
 963                             << "outside the bounds of the mesh, but no "
 964                             << "molecules have been placed in zone '"
 965                             << zone.name()
 966                             << "'.  This is likely to be because the zone "
 967                             << "has few cells ("
 968                             << zone.size()
 969                             << " in this case) and no lattice position "
 970                             << "fell inside them.  "
 971                             << "Aborting filling this zone."
 972                             << endl;
 973
 974                         break;
 975                     }
 976
 977                     molsPlacedThisIteration =
 978                         this->size() - sizeBeforeIteration;
 979
 980                     totalZoneMols += molsPlacedThisIteration;
 981
 982                     n++;
 983                 }
 984             }
 985         }
 986     }
 987 }
 988
 989
 990 void Foam::moleculeCloud::createMolecule
 991 (
 992     const point& position,
 993     label cell,
 994     label tetFace,
 995     label tetPt,
 996     label id,
 997     bool tethered,
 998     scalar temperature,
 999     const vector& bulkVelocity
1000 )
1001 {
1002     if (cell == -1)
1003     {
1004         mesh_.findCellFacePt(position, cell, tetFace, tetPt);
1005     }
1006
1007     if (cell == -1)
1008     {
1009         FatalErrorIn("Foam::moleculeCloud::createMolecule")
1010             << "Position specified does not correspond to a mesh cell." << nl
1011             << abort(FatalError);
1012     }
1013
1014     point specialPosition(vector::zero);
1015
1016     label special = 0;
1017
1018     if (tethered)
1019     {
1020         specialPosition = position;
1021
1022         special = molecule::SPECIAL_TETHERED;
1023     }
1024
1025     const molecule::constantProperties& cP(constProps(id));
1026
1027     vector v = equipartitionLinearVelocity(temperature, cP.mass());
1028
1029     v += bulkVelocity;
1030
1031     vector pi = vector::zero;
1032
1033     tensor Q = I;
1034
1035     if (!cP.pointMolecule())
1036     {
1037         pi = equipartitionAngularMomentum(temperature, cP);
1038
1039         scalar phi(rndGen_.scalar01()*twoPi);
1040
1041         scalar theta(rndGen_.scalar01()*twoPi);
1042
1043         scalar psi(rndGen_.scalar01()*twoPi);
1044
1045         Q = tensor
1046         (
1047             cos(psi)*cos(phi) - cos(theta)*sin(phi)*sin(psi),
1048             cos(psi)*sin(phi) + cos(theta)*cos(phi)*sin(psi),
1049             sin(psi)*sin(theta),
1050             - sin(psi)*cos(phi) - cos(theta)*sin(phi)*cos(psi),
1051             - sin(psi)*sin(phi) + cos(theta)*cos(phi)*cos(psi),
1052             cos(psi)*sin(theta),
1053             sin(theta)*sin(phi),
1054             - sin(theta)*cos(phi),
1055             cos(theta)
1056         );
1057     }
1058
1059     addParticle
1060     (
1061         new molecule
1062         (
1063             mesh_,
1064             position,
1065             cell,
1066             tetFace,
1067             tetPt,
1068             Q,
1069             v,
1070             vector::zero,
1071             pi,
1072             vector::zero,
1073             specialPosition,
1074             constProps(id),
1075             special,
1076             id
1077         )
1078     );
1079 }
1080
1081
1082 Foam::label Foam::moleculeCloud::nSites() const
1083 {
1084     label n = 0;
1085
1086     forAllConstIter(moleculeCloud, *this, mol)
1087     {
1088         n += constProps(mol().id()).nSites();
1089     }
1090
1091     return n;
1092 }
1093
1094
1095 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
1096
1097 Foam::moleculeCloud::moleculeCloud
1098 (
1099     const polyMesh& mesh,
1100     const potential& pot,
1101     bool readFields
1102 )
1103 :
1104     Cloud<molecule>(mesh, "moleculeCloud", false),
1105     mesh_(mesh),
1106     pot_(pot),
1107     cellOccupancy_(mesh_.nCells()),
1108     il_(mesh_, pot_.pairPotentials().rCutMax(), false),
1109     constPropList_(),
1110     rndGen_(clock::getTime())
1111 {
1112     if (readFields)
1113     {
1114         molecule::readFields(*this);
1115     }
1116
1117     buildConstProps();
1118
1119     setSiteSizesAndPositions();
1120
1121     removeHighEnergyOverlaps();
1122
1123     calculateForce();
1124 }
1125
1126
1127 Foam::moleculeCloud::moleculeCloud
1128 (
1129     const polyMesh& mesh,
1130     const potential& pot,
1131     const IOdictionary& mdInitialiseDict,
1132     bool readFields
1133 )
1134 :
1135     Cloud<molecule>(mesh, "moleculeCloud", false),
1136     mesh_(mesh),
1137     pot_(pot),
1138     il_(mesh_, 0.0, false),
1139     constPropList_(),
1140     rndGen_(clock::getTime())
1141 {
1142     if (readFields)
1143     {
1144         molecule::readFields(*this);
1145     }
1146
1147     clear();
1148
1149     buildConstProps();
1150
1151     initialiseMolecules(mdInitialiseDict);
1152 }
1153
1154
1155 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
1156
1157 void Foam::moleculeCloud::evolve()
1158 {
1159     molecule::trackingData td0(*this, 0);
1160     Cloud<molecule>::move(td0, mesh_.time().deltaTValue());
1161
1162     molecule::trackingData td1(*this, 1);
1163     Cloud<molecule>::move(td1, mesh_.time().deltaTValue());
1164
1165     molecule::trackingData td2(*this, 2);
1166     Cloud<molecule>::move(td2, mesh_.time().deltaTValue());
1167
1168     calculateForce();
1169
1170     molecule::trackingData td3(*this, 3);
1171     Cloud<molecule>::move(td3, mesh_.time().deltaTValue());
1172 }
1173
1174
1175 void Foam::moleculeCloud::calculateForce()
1176 {
1177     buildCellOccupancy();
1178
1179     // Set accumulated quantities to zero
1180     forAllIter(moleculeCloud, *this, mol)
1181     {
1182         mol().siteForces() = vector::zero;
1183
1184         mol().potentialEnergy() = 0.0;
1185
1186         mol().rf() = tensor::zero;
1187     }
1188
1189     calculatePairForce();
1190
1191     calculateTetherForce();
1192
1193     calculateExternalForce();
1194 }
1195
1196
1197 void Foam::moleculeCloud::applyConstraintsAndThermostats
1198 (
1199     const scalar targetTemperature,
1200     const scalar measuredTemperature
1201 )
1202 {
1203     scalar temperatureCorrectionFactor =
1204         sqrt(targetTemperature/max(VSMALL, measuredTemperature));
1205
1206     Info<< "----------------------------------------" << nl
1207         << "Temperature equilibration" << nl
1208         << "Target temperature = "
1209         << targetTemperature << nl
1210         << "Measured temperature = "
1211         << measuredTemperature << nl
1212         << "Temperature correction factor = "
1213         << temperatureCorrectionFactor << nl
1214         << "----------------------------------------"
1215         << endl;
1216
1217     forAllIter(moleculeCloud, *this, mol)
1218     {
1219         mol().v() *= temperatureCorrectionFactor;
1220
1221         mol().pi() *= temperatureCorrectionFactor;
1222     }
1223 }
1224
1225
1226 void Foam::moleculeCloud::writeXYZ(const fileName& fName) const
1227 {
1228     OFstream os(fName);
1229
1230     os  << nSites() << nl << "moleculeCloud site positions in angstroms" << nl;
1231
1232     forAllConstIter(moleculeCloud, *this, mol)
1233     {
1234         const molecule::constantProperties& cP = constProps(mol().id());
1235
1236         forAll(mol().sitePositions(), i)
1237         {
1238             const point& sP = mol().sitePositions()[i];
1239
1240             os  << pot_.siteIdList()[cP.siteIds()[i]]
1241                 << ' ' << sP.x()*1e10
1242                 << ' ' << sP.y()*1e10
1243                 << ' ' << sP.z()*1e10
1244                 << nl;
1245         }
1246     }
1247 }
1248
1249
1250 // ************************************************************************* //