41 void Foam::moleculeCloud::buildConstProps()
43 Info<< nl <<
"Reading moleculeProperties dictionary." <<
endl;
45 const List<word>& idList(pot_.idList());
47 constPropList_.setSize(idList.size());
49 const List<word>& siteIdList(pot_.siteIdList());
51 IOdictionary moleculePropertiesDict
56 mesh_.time().constant(),
66 const word&
id = idList[i];
67 const dictionary& molDict = moleculePropertiesDict.subDict(
id);
69 List<word> siteIdNames = molDict.lookup(
"siteIds");
71 List<label> siteIds(siteIdNames.size());
75 const word& siteId = siteIdNames[sI];
77 siteIds[sI] =
findIndex(siteIdList, siteId);
79 if (siteIds[sI] == -1)
82 << siteId <<
" site not found." 87 molecule::constantProperties&
constProp = constPropList_[i];
89 constProp = molecule::constantProperties(molDict);
91 constProp.siteIds() = siteIds;
96 void Foam::moleculeCloud::setSiteSizesAndPositions()
100 const molecule::constantProperties& cP = constProps(mol().
id());
102 mol().setSiteSizes(cP.nSites());
104 mol().setSitePositions(cP);
109 void Foam::moleculeCloud::buildCellOccupancy()
111 forAll(cellOccupancy_, cO)
113 cellOccupancy_[cO].clear();
118 cellOccupancy_[mol().cell()].append(&mol());
121 forAll(cellOccupancy_, cO)
123 cellOccupancy_[cO].shrink();
128 void Foam::moleculeCloud::calculatePairForce()
136 molecule* molI =
nullptr;
137 molecule* molJ =
nullptr;
146 forAll(cellOccupancy_[d],cellIMols)
148 molI = cellOccupancy_[d][cellIMols];
150 forAll(dil[d], interactingCells)
152 List<molecule*> cellJ =
153 cellOccupancy_[dil[d][interactingCells]];
157 molJ = cellJ[cellJMols];
159 evaluatePair(*molI, *molJ);
163 forAll(cellOccupancy_[d], cellIOtherMols)
165 molJ = cellOccupancy_[d][cellIOtherMols];
169 evaluatePair(*molI, *molJ);
177 il_.receiveReferredData(pBufs, startOfRequests);
184 List<IDLList<molecule>>& referredMols = il_.referredParticles();
188 const List<label>& realCells = ril[r];
190 IDLList<molecule>& refMols = referredMols[r];
201 List<molecule*> celli = cellOccupancy_[realCells[rC]];
205 molI = celli[cellIMols];
207 evaluatePair(*molI, refMol());
216 void Foam::moleculeCloud::calculateTetherForce()
218 const tetherPotentialList& tetherPot(pot_.tetherPotentials());
222 if (mol().tethered())
224 vector rIT = mol().position() - mol().specialPosition();
226 label idI = mol().id();
228 scalar massI = constProps(idI).mass();
230 vector fIT = tetherPot.force(idI, rIT);
232 mol().a() += fIT/massI;
234 mol().potentialEnergy() += tetherPot.energy(idI, rIT);
243 void Foam::moleculeCloud::calculateExternalForce()
247 mol().a() += pot_.gravity();
252 void Foam::moleculeCloud::removeHighEnergyOverlaps()
254 Info<< nl <<
"Removing high energy overlaps, limit = " 255 << pot_.potentialEnergyLimit()
256 << nl <<
"Removal order:";
258 forAll(pot_.removalOrder(), rO)
260 Info<<
' ' << pot_.idList()[pot_.removalOrder()[rO]];
265 label initialSize = this->size();
267 buildCellOccupancy();
271 molecule* molI =
nullptr;
272 molecule* molJ =
nullptr;
275 DynamicList<molecule*> molsToDelete;
281 forAll(cellOccupancy_[d],cellIMols)
283 molI = cellOccupancy_[d][cellIMols];
285 forAll(dil[d], interactingCells)
287 List<molecule*> cellJ =
288 cellOccupancy_[dil[d][interactingCells]];
292 molJ = cellJ[cellJMols];
294 if (evaluatePotentialLimit(*molI, *molJ))
296 label idI = molI->id();
298 label idJ = molJ->id();
309 molsToDelete.append(molJ);
312 else if (
findIndex(molsToDelete, molI) == -1)
314 molsToDelete.append(molI);
321 forAll(cellOccupancy_[d], cellIOtherMols)
323 molJ = cellOccupancy_[d][cellIOtherMols];
327 if (evaluatePotentialLimit(*molI, *molJ))
329 label idI = molI->id();
331 label idJ = molJ->id();
342 molsToDelete.append(molJ);
345 else if (
findIndex(molsToDelete, molI) == -1)
347 molsToDelete.append(molI);
356 deleteParticle(*(molsToDelete[mTD]));
360 buildCellOccupancy();
370 il_.receiveReferredData(pBufs, startOfRequests);
375 DynamicList<molecule*> molsToDelete;
379 List<IDLList<molecule>>& referredMols = il_.referredParticles();
383 IDLList<molecule>& refMols = referredMols[r];
394 const List<label>& realCells = ril[r];
398 label celli = realCells[rC];
400 List<molecule*> cellIMols = cellOccupancy_[celli];
404 molI = cellIMols[cIM];
406 if (evaluatePotentialLimit(*molI, *molJ))
408 label idI = molI->id();
410 label idJ = molJ->id();
420 molsToDelete.append(molI);
434 if (molI->origId() > molJ->origId())
438 molsToDelete.append(molI);
450 deleteParticle(*(molsToDelete[mTD]));
454 buildCellOccupancy();
462 il_.receiveReferredData(pBufs, startOfRequests);
464 label molsRemoved = initialSize - this->size();
468 reduce(molsRemoved, sumOp<label>());
471 Info<<
tab << molsRemoved <<
" molecules removed" <<
endl;
475 void Foam::moleculeCloud::initialiseMolecules
477 const IOdictionary& mdInitialiseDict
481 <<
"Initialising molecules in each zone specified in mdInitialiseDict." 486 if (!cellZones.size())
489 <<
"No cellZones found in the mesh." 495 const cellZone& zone(cellZones[z]);
499 if (!mdInitialiseDict.found(zone.name()))
501 Info<<
"No specification subDictionary for zone " 502 << zone.name() <<
" found, skipping." <<
endl;
506 const dictionary& zoneDict =
507 mdInitialiseDict.subDict(zone.name());
509 const scalar temperature
511 zoneDict.template lookup<scalar>(
"temperature")
514 const vector bulkVelocity(zoneDict.lookup(
"bulkVelocity"));
516 List<word> latticeIds
518 zoneDict.lookup(
"latticeIds")
521 List<vector> latticePositions
523 zoneDict.lookup(
"latticePositions")
526 if (latticeIds.size() != latticePositions.size())
529 <<
"latticeIds and latticePositions must be the same " 536 zoneDict.lookup(
"latticeCellShape")
539 scalar latticeCellScale = 0.0;
541 if (zoneDict.found(
"numberDensity"))
543 scalar numberDensity =
544 zoneDict.lookup<scalar>(
"numberDensity");
546 if (numberDensity < vSmall)
549 <<
"numberDensity too small, not filling zone " 550 << zone.name() <<
endl;
555 latticeCellScale =
pow 557 latticeIds.size()/(
det(latticeCellShape)*numberDensity),
561 else if (zoneDict.found(
"massDensity"))
563 scalar unitCellMass = 0.0;
569 const molecule::constantProperties& cP(constProps(
id));
571 unitCellMass += cP.mass();
574 scalar massDensity = zoneDict.lookup<scalar>(
"massDensity");
576 if (massDensity < vSmall)
579 <<
"massDensity too small, not filling zone " 580 << zone.name() <<
endl;
586 latticeCellScale =
pow 588 unitCellMass/(
det(latticeCellShape)*massDensity),
595 <<
"massDensity or numberDensity not specified " << nl
599 latticeCellShape *= latticeCellScale;
601 vector anchor(zoneDict.lookup(
"anchor"));
603 bool tethered =
false;
605 if (zoneDict.found(
"tetherSiteIds"))
609 List<word>(zoneDict.lookup(
"tetherSiteIds")).size()
613 const vector orientationAngles
615 zoneDict.lookup(
"orientationAngles")
618 scalar
phi(orientationAngles.x()*pi/180.0);
620 scalar theta(orientationAngles.y()*pi/180.0);
622 scalar
psi(orientationAngles.z()*pi/180.0);
647 const point cellCentre = mesh_.cellCentres()[zone[cell]];
649 if (cellCentre.x() > zoneMax.x())
651 zoneMax.
x() = cellCentre.x();
653 if (cellCentre.x() < zoneMin.x())
655 zoneMin.x() = cellCentre.x();
657 if (cellCentre.y() > zoneMax.y())
659 zoneMax.y() = cellCentre.y();
661 if (cellCentre.y() < zoneMin.y())
663 zoneMin.y() = cellCentre.y();
665 if (cellCentre.z() > zoneMax.z())
667 zoneMax.z() = cellCentre.z();
669 if (cellCentre.z() < zoneMin.z())
671 zoneMin.z() = cellCentre.z();
675 point zoneMid = 0.5*(zoneMin + zoneMax);
677 point latticeMid =
inv(latticeCellShape) & (
R.T()
678 & (zoneMid - anchor));
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()))
687 anchor += (
R & (latticeCellShape & latticeAnchor));
697 label totalZoneMols = 0;
699 label molsPlacedThisIteration = 0;
703 molsPlacedThisIteration != 0
704 || totalZoneMols == 0
707 label sizeBeforeIteration = this->size();
709 bool partOfLayerInBounds =
false;
726 const vector& latticePosition =
729 unitCellLatticePosition.x(),
730 unitCellLatticePosition.y(),
731 unitCellLatticePosition.z()
733 + latticePositions[
p];
735 point globalPosition =
737 + (
R & (latticeCellShape & latticePosition));
739 partOfLayerInBounds = mesh_.bounds().contains
745 mesh_.cellTree().findInside(globalPosition);
769 unitCellLatticePosition.z() = -
n;
770 unitCellLatticePosition.z() <=
n;
771 unitCellLatticePosition.z() += 2*
n 776 unitCellLatticePosition.y() = -
n;
777 unitCellLatticePosition.y() <=
n;
778 unitCellLatticePosition.y()++
783 unitCellLatticePosition.x() = -
n;
784 unitCellLatticePosition.x() <=
n;
785 unitCellLatticePosition.x()++
796 const vector& latticePosition =
799 unitCellLatticePosition.x(),
800 unitCellLatticePosition.y(),
801 unitCellLatticePosition.z()
803 + latticePositions[
p];
805 point globalPosition =
815 partOfLayerInBounds =
816 mesh_.bounds().contains
822 mesh_.cellTree().findInside
846 unitCellLatticePosition.z() = -(n-1);
847 unitCellLatticePosition.z() <= (n-1);
848 unitCellLatticePosition.z()++
851 for (
label iR = 0; iR <= 2*n -1; iR++)
853 unitCellLatticePosition.x() =
n;
855 unitCellLatticePosition.y() = -n + (iR + 1);
857 for (
label iK = 0; iK < 4; iK++)
867 const vector& latticePosition =
870 unitCellLatticePosition.x(),
871 unitCellLatticePosition.y(),
872 unitCellLatticePosition.z()
874 + latticePositions[
p];
876 point globalPosition =
886 partOfLayerInBounds =
887 mesh_.bounds().contains
893 mesh_.cellTree().findInside
912 unitCellLatticePosition =
915 - unitCellLatticePosition.y(),
916 unitCellLatticePosition.x(),
917 unitCellLatticePosition.z()
931 && !partOfLayerInBounds
935 <<
"A whole layer of unit cells was placed " 936 <<
"outside the bounds of the mesh, but no " 937 <<
"molecules have been placed in zone '" 939 <<
"'. This is likely to be because the zone " 942 <<
" in this case) and no lattice position " 943 <<
"fell inside them. " 944 <<
"Aborting filling this zone." 950 molsPlacedThisIteration =
951 this->size() - sizeBeforeIteration;
953 totalZoneMols += molsPlacedThisIteration;
963 void Foam::moleculeCloud::createMolecule
965 const point& position,
970 const vector& bulkVelocity
979 specialPosition = position;
984 const molecule::constantProperties& cP(constProps(
id));
986 vector v = equipartitionLinearVelocity(temperature, cP.mass());
994 if (!cP.pointMolecule())
996 pi = equipartitionAngularMomentum(temperature, cP);
998 scalar
phi(rndGen_.scalar01()*
twoPi);
1000 scalar theta(rndGen_.scalar01()*
twoPi);
1002 scalar
psi(rndGen_.scalar01()*
twoPi);
1045 n += constProps(mol().
id()).nSites();
1064 cellOccupancy_(mesh_.nCells()),
1065 il_(mesh_, pot_.pairPotentials().rCutMax(),
false),
1076 setSiteSizesAndPositions();
1078 removeHighEnergyOverlaps();
1095 il_(mesh_, 0.0,
false),
1108 initialiseMolecules(mdInitialiseDict);
1136 buildCellOccupancy();
1141 mol().siteForces() =
Zero;
1143 mol().potentialEnergy() = 0.0;
1148 calculatePairForce();
1150 calculateTetherForce();
1152 calculateExternalForce();
1158 const scalar targetTemperature,
1159 const scalar measuredTemperature
1162 scalar temperatureCorrectionFactor =
1163 sqrt(targetTemperature/
max(vSmall, measuredTemperature));
1165 Info<<
"----------------------------------------" << nl
1166 <<
"Temperature equilibration" << nl
1167 <<
"Target temperature = " 1168 << targetTemperature << nl
1169 <<
"Measured temperature = " 1170 << measuredTemperature << nl
1171 <<
"Temperature correction factor = " 1172 << temperatureCorrectionFactor << nl
1173 <<
"----------------------------------------" 1178 mol().v() *= temperatureCorrectionFactor;
1180 mol().pi() *= temperatureCorrectionFactor;
1189 os << nSites() << nl <<
"moleculeCloud site positions in angstroms" <<
nl;
1195 forAll(mol().sitePositions(), i)
1197 const point& sP = mol().sitePositions()[i];
1199 os << pot_.siteIdList()[cP.
siteIds()[i]]
1200 <<
' ' << sP.
x()*1e10
1201 <<
' ' << sP.
y()*1e10
1202 <<
' ' << sP.
z()*1e10
void writeXYZ(const fileName &fName) const
Write molecule sites in XYZ format.
dimensionedScalar sign(const dimensionedScalar &ds)
List< labelList > labelListList
A List of labelList.
#define forAll(list, i)
Loop across all elements in list.
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
DiagTensor< scalar > diagTensor
A scalar version of the templated DiagTensor.
A class for handling file names.
Class to hold molecule constant properties.
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
void evolve()
Evolve the molecules (move, calculate forces, control state etc)
#define forAllIter(Container, container, iter)
Iterate across all elements in the container object of type.
dimensionedSphericalTensor inv(const dimensionedSphericalTensor &dt)
const List< label > & siteIds() const
dimensionedScalar sqrt(const dimensionedScalar &ds)
Ostream & endl(Ostream &os)
Add newline and flush stream.
static time_t getTime()
Get the current clock time in seconds.
static label nRequests()
Get number of outstanding requests.
Vector< scalar > vector
A scalar version of the templated Vector.
dimensionedScalar det(const dimensionedSphericalTensor &dt)
static void readFields(Cloud< molecule > &mC)
IOdictionary is derived from dictionary and IOobject to give the dictionary automatic IO functionalit...
dimensionedScalar cos(const dimensionedScalar &ds)
static const Identity< scalar > I
void applyConstraintsAndThermostats(const scalar targetTemperature, const scalar measuredTemperature)
trackPart part() const
Return the part of the tracking operation taking place.
forAllConstIter(PtrDictionary< phaseModel >, mixture.phases(), phase)
errorManip< error > abort(error &err)
defineTemplateTypeNameAndDebug(IOPtrList< ensightPart >, 0)
const scalar twoPi(2 *pi)
dimensionedScalar sin(const dimensionedScalar &ds)
moleculeCloud(const polyMesh &mesh, const potential &pot, bool readFields=true)
Construct given mesh and potential references.
const word constProp(initialConditions.lookup("constantProperty"))
label findIndex(const ListType &, typename ListType::const_reference, const label start=0)
Find first occurrence of given element and return index,.
void reduce(const List< UPstream::commsStruct > &comms, T &Value, const BinaryOp &bop, const int tag, const label comm)
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
void move(TrackCloudType &cloud, typename ParticleType::trackingData &td, const scalar trackTime)
Move the particles.
static bool & parRun()
Is this a parallel run?
const List< DynamicList< molecule * > > & cellOccupancy
vector point
Point is a vector.
#define R(A, B, C, D, E, F, K, M)
#define WarningInFunction
Report a warning using Foam::Warning.
Class used to pass tracking data to the trackToFace function.
static const Vector< scalar > one
const volScalarField & psi
Vector< label > labelVector
Vector of labels.
Mesh consisting of general polyhedral cells.
Tensor< scalar > tensor
Tensor of scalars.
ZoneMesh< cellZone, polyMesh > cellZoneMesh
A ZoneMesh with the type cellZone.