42void Foam::moleculeCloud::buildConstProps()
44 Info<<
nl <<
"Reading moleculeProperties dictionary." <<
endl;
46 const List<word>& idList(pot_.
idList());
48 constPropList_.setSize(idList.size());
50 const List<word>& siteIdList(pot_.
siteIdList());
52 IOdictionary moleculePropertiesDict
67 const word&
id = idList[i];
68 const dictionary& molDict = moleculePropertiesDict.subDict(
id);
70 List<word> siteIdNames
72 molDict.lookup(
"siteIds")
75 List<label> siteIds(siteIdNames.size());
79 const word& siteId = siteIdNames[sI];
81 siteIds[sI] = siteIdList.find(siteId);
83 if (siteIds[sI] == -1)
86 << siteId <<
" site not found."
91 molecule::constantProperties&
constProp = constPropList_[i];
93 constProp = molecule::constantProperties(molDict);
100void Foam::moleculeCloud::setSiteSizesAndPositions()
102 for (molecule& mol : *
this)
104 const molecule::constantProperties& cP = constProps(mol.id());
106 mol.setSiteSizes(cP.nSites());
108 mol.setSitePositions(cP);
113void Foam::moleculeCloud::buildCellOccupancy()
115 for (
auto& list : cellOccupancy_)
120 for (molecule& mol : *
this)
122 cellOccupancy_[mol.cell()].append(&mol);
125 for (
auto& list : cellOccupancy_)
132void Foam::moleculeCloud::calculatePairForce()
140 molecule* molI =
nullptr;
141 molecule* molJ =
nullptr;
150 forAll(cellOccupancy_[d],cellIMols)
152 molI = cellOccupancy_[d][cellIMols];
154 forAll(dil[d], interactingCells)
156 List<molecule*> cellJ =
157 cellOccupancy_[dil[d][interactingCells]];
161 molJ = cellJ[cellJMols];
163 evaluatePair(*molI, *molJ);
167 forAll(cellOccupancy_[d], cellIOtherMols)
169 molJ = cellOccupancy_[d][cellIOtherMols];
173 evaluatePair(*molI, *molJ);
181 il_.receiveReferredData(pBufs, startOfRequests);
188 List<IDLList<molecule>>& referredMols = il_.referredParticles();
192 const List<label>& realCells = ril[r];
194 IDLList<molecule>& refMols = referredMols[r];
196 for (molecule& refMol : refMols)
200 List<molecule*> celli = cellOccupancy_[realCells[rC]];
204 molI = celli[cellIMols];
206 evaluatePair(*molI, refMol);
215void Foam::moleculeCloud::calculateTetherForce()
217 const tetherPotentialList& tetherPot(pot_.tetherPotentials());
219 for (molecule& mol : *
this)
223 vector rIT = mol.position() - mol.specialPosition();
225 label idI = mol.id();
227 scalar massI = constProps(idI).mass();
229 vector fIT = tetherPot.force(idI, rIT);
231 mol.a() += fIT/massI;
233 mol.potentialEnergy() += tetherPot.energy(idI, rIT);
242void Foam::moleculeCloud::calculateExternalForce()
244 for (molecule& mol : *
this)
246 mol.a() += pot_.gravity();
251void Foam::moleculeCloud::removeHighEnergyOverlaps()
253 Info<<
nl <<
"Removing high energy overlaps, limit = "
254 << pot_.potentialEnergyLimit()
255 <<
nl <<
"Removal order:";
257 forAll(pot_.removalOrder(), rO)
259 Info<<
' ' << pot_.idList()[pot_.removalOrder()[rO]];
264 label initialSize = this->size();
266 buildCellOccupancy();
270 molecule* molI =
nullptr;
271 molecule* molJ =
nullptr;
274 DynamicList<molecule*> molsToDelete;
280 forAll(cellOccupancy_[d],cellIMols)
282 molI = cellOccupancy_[d][cellIMols];
284 forAll(dil[d], interactingCells)
286 List<molecule*> cellJ =
287 cellOccupancy_[dil[d][interactingCells]];
291 molJ = cellJ[cellJMols];
293 if (evaluatePotentialLimit(*molI, *molJ))
295 const label idI = molI->id();
296 const label idJ = molJ->id();
301 || pot_.removalOrder().find(idJ)
302 < pot_.removalOrder().find(idI)
305 if (!molsToDelete.found(molJ))
307 molsToDelete.append(molJ);
310 else if (!molsToDelete.found(molI))
312 molsToDelete.append(molI);
319 forAll(cellOccupancy_[d], cellIOtherMols)
321 molJ = cellOccupancy_[d][cellIOtherMols];
325 if (evaluatePotentialLimit(*molI, *molJ))
327 const label idI = molI->id();
328 const label idJ = molJ->id();
333 || pot_.removalOrder().find(idJ)
334 < pot_.removalOrder().find(idI)
337 if (!molsToDelete.found(molJ))
339 molsToDelete.append(molJ);
342 else if (!molsToDelete.found(molI))
344 molsToDelete.append(molI);
353 deleteParticle(*(molsToDelete[mTD]));
357 buildCellOccupancy();
367 il_.receiveReferredData(pBufs, startOfRequests);
372 DynamicList<molecule*> molsToDelete;
376 List<IDLList<molecule>>& referredMols = il_.referredParticles();
380 IDLList<molecule>& refMols = referredMols[r];
382 for (molecule& refMol : refMols)
386 const List<label>& realCells = ril[r];
390 label celli = realCells[rC];
392 List<molecule*> cellIMols = cellOccupancy_[celli];
396 molI = cellIMols[cIM];
398 if (evaluatePotentialLimit(*molI, *molJ))
400 const label idI = molI->id();
401 const label idJ = molJ->id();
405 pot_.removalOrder().find(idI)
406 < pot_.removalOrder().find(idJ)
409 if (!molsToDelete.found(molI))
411 molsToDelete.append(molI);
416 pot_.removalOrder().find(idI)
417 == pot_.removalOrder().find(idJ)
425 if (molI->origId() > molJ->origId())
427 if (!molsToDelete.found(molI))
429 molsToDelete.append(molI);
441 deleteParticle(*(molsToDelete[mTD]));
445 buildCellOccupancy();
453 il_.receiveReferredData(pBufs, startOfRequests);
455 label molsRemoved = initialSize - this->size();
459 reduce(molsRemoved, sumOp<label>());
462 Info<<
tab << molsRemoved <<
" molecules removed" <<
endl;
466void Foam::moleculeCloud::initialiseMolecules
468 const IOdictionary& mdInitialiseDict
472 <<
"Initialising molecules in each zone specified in mdInitialiseDict."
477 if (!cellZones.size())
480 <<
"No cellZones found in the mesh."
484 for (
const cellZone& zone : cellZones)
488 if (!mdInitialiseDict.found(zone.name()))
490 Info<<
"No specification subDictionary for zone "
491 << zone.name() <<
" found, skipping." <<
endl;
495 const dictionary& zoneDict =
496 mdInitialiseDict.subDict(zone.name());
498 const scalar temperature
500 zoneDict.get<scalar>(
"temperature")
505 zoneDict.get<
vector>(
"bulkVelocity")
508 List<word> latticeIds
510 zoneDict.lookup(
"latticeIds")
513 List<vector> latticePositions
515 zoneDict.lookup(
"latticePositions")
518 if (latticeIds.size() != latticePositions.size())
521 <<
"latticeIds and latticePositions must be the same "
528 zoneDict.lookup(
"latticeCellShape")
531 scalar latticeCellScale = 0.0;
533 if (zoneDict.found(
"numberDensity"))
535 const scalar numberDensity
537 zoneDict.get<scalar>(
"numberDensity")
540 if (numberDensity < VSMALL)
543 <<
"numberDensity too small, not filling zone "
544 << zone.name() <<
endl;
549 latticeCellScale =
pow
551 latticeIds.size()/(
det(latticeCellShape)*numberDensity),
555 else if (zoneDict.found(
"massDensity"))
557 scalar unitCellMass = 0.0;
561 label
id = pot_.idList().find(latticeIds[i]);
563 const molecule::constantProperties& cP(constProps(
id));
565 unitCellMass += cP.mass();
568 const scalar massDensity
570 zoneDict.get<scalar>(
"massDensity")
573 if (massDensity < VSMALL)
576 <<
"massDensity too small, not filling zone "
577 << zone.name() <<
endl;
583 latticeCellScale =
pow
585 unitCellMass/(
det(latticeCellShape)*massDensity),
592 <<
"massDensity or numberDensity not specified " <<
nl
596 latticeCellShape *= latticeCellScale;
600 bool tethered =
false;
602 if (zoneDict.found(
"tetherSiteIds"))
606 List<word>(zoneDict.lookup(
"tetherSiteIds")).size()
610 const vector orientationAngles
612 zoneDict.get<
vector>(
"orientationAngles")
616 const scalar theta(
degToRad(orientationAngles.y()));
636 vector zoneMin = VGREAT*vector::one;
638 vector zoneMax = -VGREAT*vector::one;
642 const point cellCentre = mesh_.cellCentres()[zone[cell]];
644 if (cellCentre.x() > zoneMax.x())
646 zoneMax.x() = cellCentre.x();
648 if (cellCentre.x() < zoneMin.x())
650 zoneMin.x() = cellCentre.x();
652 if (cellCentre.y() > zoneMax.y())
654 zoneMax.y() = cellCentre.y();
656 if (cellCentre.y() < zoneMin.y())
658 zoneMin.y() = cellCentre.y();
660 if (cellCentre.z() > zoneMax.z())
662 zoneMax.z() = cellCentre.z();
664 if (cellCentre.z() < zoneMin.z())
666 zoneMin.z() = cellCentre.z();
670 point zoneMid = 0.5*(zoneMin + zoneMax);
672 point latticeMid =
inv(latticeCellShape) & (
R.T()
673 & (zoneMid - anchor));
677 label(latticeMid.x() + 0.5*
sign(latticeMid.x())),
678 label(latticeMid.y() + 0.5*
sign(latticeMid.y())),
679 label(latticeMid.z() + 0.5*
sign(latticeMid.z()))
682 anchor += (
R & (latticeCellShape & latticeAnchor));
692 label totalZoneMols = 0;
694 label molsPlacedThisIteration = 0;
698 molsPlacedThisIteration != 0
699 || totalZoneMols == 0
702 label sizeBeforeIteration = this->size();
704 bool partOfLayerInBounds =
false;
719 label
id = pot_.idList().find(latticeIds[
p]);
721 const vector& latticePosition =
724 unitCellLatticePosition.x(),
725 unitCellLatticePosition.y(),
726 unitCellLatticePosition.z()
728 + latticePositions[
p];
730 point globalPosition =
732 + (
R & (latticeCellShape & latticePosition));
734 partOfLayerInBounds = mesh_.bounds().contains
740 mesh_.cellTree().findInside(globalPosition);
742 if (zone.found(cell))
764 unitCellLatticePosition.z() = -
n;
765 unitCellLatticePosition.z() <=
n;
766 unitCellLatticePosition.z() += 2*
n
771 unitCellLatticePosition.y() = -
n;
772 unitCellLatticePosition.y() <=
n;
773 unitCellLatticePosition.y()++
778 unitCellLatticePosition.x() = -
n;
779 unitCellLatticePosition.x() <=
n;
780 unitCellLatticePosition.x()++
786 pot_.idList().find(latticeIds[
p]);
788 const vector& latticePosition =
791 unitCellLatticePosition.x(),
792 unitCellLatticePosition.y(),
793 unitCellLatticePosition.z()
795 + latticePositions[
p];
797 point globalPosition =
807 partOfLayerInBounds =
808 mesh_.bounds().contains
814 mesh_.cellTree().findInside
819 if (zone.found(cell))
838 unitCellLatticePosition.z() = -(
n-1);
839 unitCellLatticePosition.z() <= (
n-1);
840 unitCellLatticePosition.z()++
843 for (label iR = 0; iR <= 2*
n -1; iR++)
845 unitCellLatticePosition.x() =
n;
847 unitCellLatticePosition.y() = -
n + (iR + 1);
849 for (label iK = 0; iK < 4; iK++)
854 pot_.idList().find(latticeIds[
p]);
856 const vector& latticePosition =
859 unitCellLatticePosition.x(),
860 unitCellLatticePosition.y(),
861 unitCellLatticePosition.z()
863 + latticePositions[
p];
865 point globalPosition =
875 partOfLayerInBounds =
876 mesh_.bounds().contains
882 mesh_.cellTree().findInside
887 if (zone.found(cell))
901 unitCellLatticePosition =
904 - unitCellLatticePosition.y(),
905 unitCellLatticePosition.x(),
906 unitCellLatticePosition.z()
920 && !partOfLayerInBounds
924 <<
"A whole layer of unit cells was placed "
925 <<
"outside the bounds of the mesh, but no "
926 <<
"molecules have been placed in zone '"
928 <<
"'. This is likely to be because the zone "
931 <<
" in this case) and no lattice position "
932 <<
"fell inside them. "
933 <<
"Aborting filling this zone."
939 molsPlacedThisIteration =
940 this->size() - sizeBeforeIteration;
942 totalZoneMols += molsPlacedThisIteration;
952void Foam::moleculeCloud::createMolecule
954 const point& position,
959 const vector& bulkVelocity
968 specialPosition = position;
973 const molecule::constantProperties& cP(constProps(
id));
975 vector v = equipartitionLinearVelocity(temperature, cP.mass());
983 if (!cP.pointMolecule())
985 pi = equipartitionAngularMomentum(temperature, cP);
1026Foam::label Foam::moleculeCloud::nSites()
const
1030 for (
const molecule& mol : *
this)
1032 n += constProps(mol.id()).nSites();
1051 cellOccupancy_(mesh_.nCells()),
1052 il_(mesh_, pot_.pairPotentials().rCutMax(), false),
1054 rndGen_(
clock::getTime())
1063 setSiteSizesAndPositions();
1065 removeHighEnergyOverlaps();
1082 il_(mesh_, 0.0, false),
1084 rndGen_(
clock::getTime())
1095 initialiseMolecules(mdInitialiseDict);
1121 buildCellOccupancy();
1133 calculatePairForce();
1135 calculateTetherForce();
1137 calculateExternalForce();
1143 const scalar targetTemperature,
1144 const scalar measuredTemperature
1147 scalar temperatureCorrectionFactor =
1148 sqrt(targetTemperature/
max(VSMALL, measuredTemperature));
1150 Info<<
"----------------------------------------" <<
nl
1151 <<
"Temperature equilibration" <<
nl
1152 <<
"Target temperature = "
1153 << targetTemperature <<
nl
1154 <<
"Measured temperature = "
1155 << measuredTemperature <<
nl
1156 <<
"Temperature correction factor = "
1157 << temperatureCorrectionFactor <<
nl
1158 <<
"----------------------------------------"
1163 mol.
v() *= temperatureCorrectionFactor;
1165 mol.
pi() *= temperatureCorrectionFactor;
1174 os << nSites() <<
nl <<
"moleculeCloud site positions in angstroms" <<
nl;
1184 os << pot_.siteIdList()[cP.
siteIds()[i]]
1185 <<
' ' << sP.
x()*1e10
1186 <<
' ' << sP.
y()*1e10
1187 <<
' ' << sP.
z()*1e10
#define R(A, B, C, D, E, F, K, M)
const List< DynamicList< molecule * > > & cellOccupancy
Base cloud calls templated on particle type.
void clear()
Clear the contents of the list.
IOdictionary is derived from dictionary and IOobject to give the dictionary automatic IO functionalit...
Output to file stream, using an OSstream.
const word & constant() const
Return constant name.
@ nonBlocking
"nonBlocking"
static label nRequests()
Get number of outstanding requests.
static bool & parRun() noexcept
Test if this a parallel run.
const Cmpt & z() const
Access to the vector z component.
const Cmpt & y() const
Access to the vector y component.
const Cmpt & x() const
Access to the vector x component.
Read access to the system clock with formatting.
A class for handling file names.
Reads fields from the time directories and adds them to the mesh database for further post-processing...
void writeXYZ(const fileName &fName) const
Write molecule sites in XYZ format.
void evolve()
Evolve the molecules (move, calculate forces, control state etc)
void applyConstraintsAndThermostats(const scalar targetTemperature, const scalar measuredTemperature)
Class to hold molecule constant properties.
const List< label > & siteIds() const
Class used to pass tracking data to the trackToFace function.
const List< vector > & sitePositions() const
const List< vector > & siteForces() const
const tensor & rf() const
scalar potentialEnergy() const
const vector & pi() const
const Time & time() const noexcept
Return time registry.
Mesh consisting of general polyhedral cells.
const List< word > & siteIdList() const
const List< word > & idList() const
Tensor of scalars, i.e. Tensor<scalar>.
A Vector of values with scalar precision, where scalar is float/double depending on the compilation f...
#define defineTemplateTypeNameAndDebug(Type, DebugSwitch)
Define the typeName and debug information for templates, useful.
const volScalarField & psi
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
OBJstream os(runTime.globalPath()/outputName)
#define WarningInFunction
Report a warning using Foam::Warning.
constexpr scalar pi(M_PI)
constexpr scalar twoPi(2 *M_PI)
dimensionedScalar det(const dimensionedSphericalTensor &dt)
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
dimensionedScalar sign(const dimensionedScalar &ds)
ZoneMesh< cellZone, polyMesh > cellZoneMesh
A ZoneMesh with the type cellZone.
Vector< label > labelVector
Vector of labels.
dimensionedScalar sin(const dimensionedScalar &ds)
messageStream Info
Information stream (stdout output on master, null elsewhere)
static const Identity< scalar > I
void readFields(const typename GeoFieldType::Mesh &mesh, const IOobjectList &objects, const wordHashSet &selectedFields, LIFOStack< regIOobject * > &storedObjects)
Read the selected GeometricFields of the templated type.
vector point
Point is a vector.
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
DiagTensor< scalar > diagTensor
DiagTensor of scalars, i.e. DiagTensor<scalar>.
Ostream & endl(Ostream &os)
Add newline and flush stream.
constexpr scalar degToRad() noexcept
Multiplication factor for degrees to radians conversion.
dimensionedScalar sqrt(const dimensionedScalar &ds)
List< labelList > labelListList
A List of labelList.
void reduce(const List< UPstream::commsStruct > &comms, T &value, const BinaryOp &bop, const int tag, const label comm)
errorManip< error > abort(error &err)
static constexpr const zero Zero
Global zero (0)
dimensionedSphericalTensor inv(const dimensionedSphericalTensor &dt)
dimensionedScalar cos(const dimensionedScalar &ds)
constexpr char nl
The newline '\n' character (0x0a)
constexpr char tab
The tab '\t' character(0x09)
word constProp(initialConditions.get< word >("constantProperty"))
#define forAll(list, i)
Loop across all elements in list.
Unit conversion functions.