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