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