34inline void Foam::moleculeCloud::evaluatePair
44 label idI = molI.
id();
46 label idJ = molJ.
id();
56 List<bool> pairPotentialSitesI = constPropI.pairPotentialSites();
58 List<bool> electrostaticSitesI = constPropI.electrostaticSites();
60 List<bool> pairPotentialSitesJ = constPropJ.pairPotentialSites();
62 List<bool> electrostaticSitesJ = constPropJ.electrostaticSites();
66 label idsI(siteIdsI[sI]);
70 label idsJ(siteIdsJ[sJ]);
72 if (pairPotentialSitesI[sI] && pairPotentialSitesJ[sJ])
77 scalar rsIsJMagSq =
magSqr(rsIsJ);
79 if (pairPot.
rCutSqr(idsI, idsJ, rsIsJMagSq))
81 scalar rsIsJMag =
mag(rsIsJ);
85 *pairPot.
force(idsI, idsJ, rsIsJMag);
91 scalar potentialEnergy
93 pairPot.
energy(idsI, idsJ, rsIsJMag)
102 tensor virialContribution =
103 (rsIsJ*fsIsJ)*(rsIsJ & rIJ)/rsIsJMagSq;
105 molI.
rf() += virialContribution;
107 molJ.
rf() += virialContribution;
111 if (electrostaticSitesI[sI] && electrostaticSitesJ[sJ])
116 scalar rsIsJMagSq =
magSqr(rsIsJ);
118 if (rsIsJMagSq <= electrostatic.
rCutSqr())
120 scalar rsIsJMag =
mag(rsIsJ);
122 scalar chargeI = constPropI.siteCharges()[sI];
124 scalar chargeJ = constPropJ.siteCharges()[sJ];
128 *chargeI*chargeJ*electrostatic.
force(rsIsJMag);
134 scalar potentialEnergy =
136 *electrostatic.
energy(rsIsJMag);
144 tensor virialContribution =
145 (rsIsJ*fsIsJ)*(rsIsJ & rIJ)/rsIsJMagSq;
147 molI.
rf() += virialContribution;
149 molJ.
rf() += virialContribution;
157inline bool Foam::moleculeCloud::evaluatePotentialLimit
167 label idI = molI.
id();
169 label idJ = molJ.
id();
179 List<bool> pairPotentialSitesI = constPropI.pairPotentialSites();
181 List<bool> electrostaticSitesI = constPropI.electrostaticSites();
183 List<bool> pairPotentialSitesJ = constPropJ.pairPotentialSites();
185 List<bool> electrostaticSitesJ = constPropJ.electrostaticSites();
189 label idsI(siteIdsI[sI]);
193 label idsJ(siteIdsJ[sJ]);
195 if (pairPotentialSitesI[sI] && pairPotentialSitesJ[sJ])
200 scalar rsIsJMagSq =
magSqr(rsIsJ);
202 if (pairPot.
rCutSqr(idsI, idsJ, rsIsJMagSq))
204 scalar rsIsJMag =
mag(rsIsJ);
210 if (rsIsJMag < SMALL)
213 <<
"Molecule site pair closer than "
215 <<
": mag separation = " << rsIsJMag
216 <<
". These may have been placed on top of each"
217 <<
" other by a rounding error in mdInitialise in"
218 <<
" parallel or a block filled with molecules"
219 <<
" twice. Removing one of the molecules."
228 if (rsIsJMag < pairPot.
rMin(idsI, idsJ))
235 mag(pairPot.
energy(idsI, idsJ, rsIsJMag))
236 > pot_.potentialEnergyLimit()
244 if (electrostaticSitesI[sI] && electrostaticSitesJ[sJ])
249 scalar rsIsJMagSq =
magSqr(rsIsJ);
253 scalar rsIsJMag =
mag(rsIsJ);
259 if (rsIsJMag < SMALL)
262 <<
"Molecule site pair closer than "
264 <<
": mag separation = " << rsIsJMag
265 <<
". These may have been placed on top of each"
266 <<
" other by a rounding error in mdInitialise in"
267 <<
" parallel or a block filled with molecules"
268 <<
" twice. Removing one of the molecules."
274 if (rsIsJMag < electrostatic.
rMin())
279 scalar chargeI = constPropI.siteCharges()[sI];
281 scalar chargeJ = constPropJ.siteCharges()[sJ];
285 mag(chargeI*chargeJ*electrostatic.
energy(rsIsJMag))
286 > pot_.potentialEnergyLimit()
300inline Foam::vector Foam::moleculeCloud::equipartitionLinearVelocity
308 *rndGen_.GaussNormal<
vector>();
312inline Foam::vector Foam::moleculeCloud::equipartitionAngularMomentum
358 return cellOccupancy_;
372 return constPropList_;
379 return constPropList_[id];
Builds direct interaction list, specifying which local (real) cells are potentially in range of each ...
A 1D array of objects of type <T>, where the size of the vector is known and used for subscript bound...
const potential & pot() const
const List< molecule::constantProperties > constProps() const
const List< DynamicList< molecule * > > & cellOccupancy() const
const polyMesh & mesh() const
const InteractionLists< molecule > & il() const
Class to hold molecule constant properties.
const diagTensor & momentOfInertia() const
bool linearMolecule() const
const List< vector > & sitePositions() const
const List< vector > & siteForces() const
const tensor & rf() const
scalar potentialEnergy() const
scalar energy(const label a, const label b, const scalar rIJMag) const
scalar force(const label a, const label b, const scalar rIJMag) const
scalar rCutMaxSqr() const
const pairPotential & electrostatic() const
bool rCutSqr(const label a, const label b, const scalar rIJMagSqr) const
scalar rMin(const label a, const label b) const
scalar energy(const scalar r) const
scalar force(const scalar r) const
vector position() const
Return current particle position.
Mesh consisting of general polyhedral cells.
const pairPotentialList & pairPotentials() 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 WarningInFunction
Report a warning using Foam::Warning.
const dimensionedScalar k
Boltzmann constant.
Different types of constants.
Ostream & endl(Ostream &os)
Add newline and flush stream.
dimensionedScalar sqrt(const dimensionedScalar &ds)
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
dimensioned< typename typeOfMag< Type >::type > magSqr(const dimensioned< Type > &dt)
#define forAll(list, i)
Loop across all elements in list.