42 void 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);
100 void 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);
113 void Foam::moleculeCloud::buildCellOccupancy()
115 for (
auto& list : cellOccupancy_)
120 for (molecule& mol : *
this)
122 cellOccupancy_[mol.cell()].append(&mol);
125 for (
auto& list : cellOccupancy_)
132 void Foam::moleculeCloud::calculatePairForce()
134 PstreamBuffers pBufs;
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);
215 void 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);
242 void Foam::moleculeCloud::calculateExternalForce()
244 for (molecule& mol : *
this)
246 mol.a() += pot_.gravity();
251 void 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();
359 PstreamBuffers pBufs;
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;
466 void 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;
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;
952 void 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);
1026 Foam::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;
1176 for (
const molecule& mol : *
this)
1178 const molecule::constantProperties& cP = constProps(mol.id());
1180 forAll(mol.sitePositions(), i)
1182 const point& sP = mol.sitePositions()[i];
1184 os << pot_.siteIdList()[cP.siteIds()[i]]
1185 <<
' ' << sP.x()*1e10
1186 <<
' ' << sP.y()*1e10
1187 <<
' ' << sP.z()*1e10
scalar potentialEnergy() const
void writeXYZ(const fileName &fName) const
Write molecule sites in XYZ format.
dimensionedScalar sign(const dimensionedScalar &ds)
DiagTensor< scalar > diagTensor
DiagTensor of scalars, i.e. DiagTensor<scalar>.
const tensor & rf() const
error FatalError
Error stream (stdout output on all processes), with additional 'FOAM FATAL ERROR' header text and sta...
const List< vector > & siteForces() const
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
void evolve()
Evolve the molecules (move, calculate forces, control state etc)
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
Unit conversion functions.
static label nRequests() noexcept
Number of outstanding requests (on the internal list of requests)
dimensionedSphericalTensor inv(const dimensionedSphericalTensor &dt)
constexpr char nl
The newline '\n' character (0x0a)
Read access to the system clock with formatting.
dimensionedScalar sqrt(const dimensionedScalar &ds)
Ostream & endl(Ostream &os)
Add newline and flush stream.
static bool & parRun() noexcept
Test if this a parallel run.
constexpr char tab
The tab '\t' character(0x09)
const vector & pi() const
dimensionedScalar det(const dimensionedSphericalTensor &dt)
Ignore writing from objectRegistry::writeObject()
List< labelList > labelListList
List of labelList.
static void readFields(Cloud< molecule > &mC)
IOdictionary is derived from dictionary and IOobject to give the dictionary automatic IO functionalit...
#define forAll(list, i)
Loop across all elements in list.
dimensionedScalar cos(const dimensionedScalar &ds)
static const Identity< scalar > I
void readFields(const typename GeoFieldType::Mesh &mesh, const IOobjectList &objects, const NameMatchPredicate &selectedFields, DynamicList< regIOobject *> &storedObjects)
Read the selected GeometricFields of the templated type and store on the objectRegistry.
constexpr scalar twoPi(2 *M_PI)
const Time & time() const noexcept
Return time registry.
void clear()
Clear the contents of the list.
void applyConstraintsAndThermostats(const scalar targetTemperature, const scalar measuredTemperature)
constexpr scalar pi(M_PI)
moleculeCloud(const moleculeCloud &)=delete
No copy construct.
errorManip< error > abort(error &err)
word constProp(initialConditions.get< word >("constantProperty"))
Base cloud calls templated on particle type.
const List< word > & idList() const
dimensionedScalar sin(const dimensionedScalar &ds)
const word & constant() const noexcept
Return constant name.
OBJstream os(runTime.globalPath()/outputName)
defineTemplateTypeNameAndDebug(faScalarMatrix, 0)
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
void move(TrackCloudType &cloud, typename ParticleType::trackingData &td, const scalar trackTime)
Move the particles.
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.
messageStream Info
Information stream (stdout output on master, null elsewhere)
const volScalarField & psi
Vector< label > labelVector
Vector of labels.
Mesh consisting of general polyhedral cells.
constexpr scalar degToRad(const scalar deg) noexcept
Conversion from degrees to radians.
Do not request registration (bool: false)
ZoneMesh< cellZone, polyMesh > cellZoneMesh
A ZoneMesh with the type cellZone.
void reduce(T &value, const BinaryOp &bop, const int tag=UPstream::msgType(), const label comm=UPstream::worldComm)
Reduce inplace (cf. MPI Allreduce) using linear/tree communication schedule.
const List< word > & siteIdList() const
static constexpr const zero Zero
Global zero (0)