32 inline void Foam::moleculeCloud::evaluatePair
38 const pairPotentialList& pairPot = pot_.pairPotentials();
40 const pairPotential& electrostatic = pairPot.electrostatic();
42 label idI = molI.id();
44 label idJ = molJ.id();
46 const molecule::constantProperties& constPropI(constProps(idI));
48 const molecule::constantProperties& constPropJ(constProps(idJ));
50 List<label> siteIdsI = constPropI.siteIds();
52 List<label> siteIdsJ = constPropJ.siteIds();
54 List<bool> pairPotentialSitesI = constPropI.pairPotentialSites();
56 List<bool> electrostaticSitesI = constPropI.electrostaticSites();
58 List<bool> pairPotentialSitesJ = constPropJ.pairPotentialSites();
60 List<bool> electrostaticSitesJ = constPropJ.electrostaticSites();
64 label idsI(siteIdsI[sI]);
68 label idsJ(siteIdsJ[sJ]);
70 if (pairPotentialSitesI[sI] && pairPotentialSitesJ[sJ])
73 molI.sitePositions()[sI] - molJ.sitePositions()[sJ];
75 scalar rsIsJMagSq =
magSqr(rsIsJ);
77 if (pairPot.rCutSqr(idsI, idsJ, rsIsJMagSq))
79 scalar rsIsJMag =
mag(rsIsJ);
83 *pairPot.force(idsI, idsJ, rsIsJMag);
85 molI.siteForces()[sI] += fsIsJ;
87 molJ.siteForces()[sJ] += -fsIsJ;
89 scalar potentialEnergy
91 pairPot.energy(idsI, idsJ, rsIsJMag)
94 molI.potentialEnergy() += 0.5*potentialEnergy;
96 molJ.potentialEnergy() += 0.5*potentialEnergy;
98 vector rIJ = molI.position() - molJ.position();
100 tensor virialContribution =
101 (rsIsJ*fsIsJ)*(rsIsJ & rIJ)/rsIsJMagSq;
103 molI.rf() += virialContribution;
105 molJ.rf() += virialContribution;
109 if (electrostaticSitesI[sI] && electrostaticSitesJ[sJ])
112 molI.sitePositions()[sI] - molJ.sitePositions()[sJ];
114 scalar rsIsJMagSq =
magSqr(rsIsJ);
116 if (rsIsJMagSq <= electrostatic.rCutSqr())
118 scalar rsIsJMag =
mag(rsIsJ);
120 scalar chargeI = constPropI.siteCharges()[sI];
122 scalar chargeJ = constPropJ.siteCharges()[sJ];
126 *chargeI*chargeJ*electrostatic.force(rsIsJMag);
128 molI.siteForces()[sI] += fsIsJ;
130 molJ.siteForces()[sJ] += -fsIsJ;
132 scalar potentialEnergy =
134 *electrostatic.energy(rsIsJMag);
136 molI.potentialEnergy() += 0.5*potentialEnergy;
138 molJ.potentialEnergy() += 0.5*potentialEnergy;
140 vector rIJ = molI.position() - molJ.position();
142 tensor virialContribution =
143 (rsIsJ*fsIsJ)*(rsIsJ & rIJ)/rsIsJMagSq;
145 molI.rf() += virialContribution;
147 molJ.rf() += virialContribution;
155 inline bool Foam::moleculeCloud::evaluatePotentialLimit
161 const pairPotentialList& pairPot = pot_.pairPotentials();
163 const pairPotential& electrostatic = pairPot.electrostatic();
165 label idI = molI.id();
167 label idJ = molJ.id();
169 const molecule::constantProperties& constPropI(constProps(idI));
171 const molecule::constantProperties& constPropJ(constProps(idJ));
173 List<label> siteIdsI = constPropI.siteIds();
175 List<label> siteIdsJ = constPropJ.siteIds();
177 List<bool> pairPotentialSitesI = constPropI.pairPotentialSites();
179 List<bool> electrostaticSitesI = constPropI.electrostaticSites();
181 List<bool> pairPotentialSitesJ = constPropJ.pairPotentialSites();
183 List<bool> electrostaticSitesJ = constPropJ.electrostaticSites();
187 label idsI(siteIdsI[sI]);
191 label idsJ(siteIdsJ[sJ]);
193 if (pairPotentialSitesI[sI] && pairPotentialSitesJ[sJ])
196 molI.sitePositions()[sI] - molJ.sitePositions()[sJ];
198 scalar rsIsJMagSq =
magSqr(rsIsJ);
200 if (pairPot.rCutSqr(idsI, idsJ, rsIsJMagSq))
202 scalar rsIsJMag =
mag(rsIsJ);
208 if (rsIsJMag < small)
211 <<
"Molecule site pair closer than " 213 <<
": mag separation = " << rsIsJMag
214 <<
". These may have been placed on top of each" 215 <<
" other by a rounding error in mdInitialise in" 216 <<
" parallel or a block filled with molecules" 217 <<
" twice. Removing one of the molecules." 226 if (rsIsJMag < pairPot.rMin(idsI, idsJ))
233 mag(pairPot.energy(idsI, idsJ, rsIsJMag))
234 > pot_.potentialEnergyLimit()
242 if (electrostaticSitesI[sI] && electrostaticSitesJ[sJ])
245 molI.sitePositions()[sI] - molJ.sitePositions()[sJ];
247 scalar rsIsJMagSq =
magSqr(rsIsJ);
249 if (pairPot.rCutMaxSqr(rsIsJMagSq))
251 scalar rsIsJMag =
mag(rsIsJ);
257 if (rsIsJMag < small)
260 <<
"Molecule site pair closer than " 262 <<
": mag separation = " << rsIsJMag
263 <<
". These may have been placed on top of each" 264 <<
" other by a rounding error in mdInitialise in" 265 <<
" parallel or a block filled with molecules" 266 <<
" twice. Removing one of the molecules." 272 if (rsIsJMag < electrostatic.rMin())
277 scalar chargeI = constPropI.siteCharges()[sI];
279 scalar chargeJ = constPropJ.siteCharges()[sJ];
283 mag(chargeI*chargeJ*electrostatic.energy(rsIsJMag))
284 > pot_.potentialEnergyLimit()
298 inline Foam::vector Foam::moleculeCloud::equipartitionLinearVelocity
306 *rndGen_.sampleNormal<
vector>();
310 inline Foam::vector Foam::moleculeCloud::equipartitionAngularMomentum
313 const molecule::constantProperties& cP
318 if (cP.linearMolecule())
323 sqrt(cP.momentOfInertia().yy())*rndGen_.scalarNormal(),
324 sqrt(cP.momentOfInertia().zz())*rndGen_.scalarNormal()
331 sqrt(cP.momentOfInertia().xx())*rndGen_.scalarNormal(),
332 sqrt(cP.momentOfInertia().yy())*rndGen_.scalarNormal(),
333 sqrt(cP.momentOfInertia().zz())*rndGen_.scalarNormal()
356 return cellOccupancy_;
370 return constPropList_;
377 return constPropList_[id];
#define forAll(list, i)
Loop across all elements in list.
const InteractionLists< molecule > & il() const
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Class to hold molecule constant properties.
A 1D array of objects of type <T>, where the size of the vector is known and used for subscript bound...
dimensionedScalar sqrt(const dimensionedScalar &ds)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Vector< scalar > vector
A scalar version of the templated Vector.
const List< DynamicList< molecule * > > & cellOccupancy() const
dimensioned< scalar > magSqr(const dimensioned< Type > &)
const potential & pot() const
const dimensionedScalar k
Boltzmann constant.
#define WarningInFunction
Report a warning using Foam::Warning.
dimensioned< scalar > mag(const dimensioned< Type > &)
const List< molecule::constantProperties > constProps() const
Mesh consisting of general polyhedral cells.
const polyMesh & mesh() const
Tensor< scalar > tensor
Tensor of scalars.