31 siteMasses_(
List<scalar>(0)),
32 siteCharges_(
List<scalar>(0)),
33 siteIds_(
List<label>(0)),
34 pairPotentialSites_(
List<bool>(false)),
35 electrostaticSites_(
List<bool>(false)),
46 siteReferencePositions_(dict.
lookup(
"siteReferencePositions")),
47 siteMasses_(dict.
lookup(
"siteMasses")),
48 siteCharges_(dict.
lookup(
"siteCharges")),
50 pairPotentialSites_(),
51 electrostaticSites_(),
57 setInteracionSiteBools
63 mass_ =
sum(siteMasses_);
70 forAll(siteReferencePositions_, i)
72 centreOfMass += siteReferencePositions_[i]*siteMasses_[i];
75 centreOfMass /= mass_;
77 siteReferencePositions_ -= centreOfMass;
79 if(siteIds_.size() == 1)
87 else if(linearMoleculeTest())
93 vector dir = siteReferencePositions_[1] - siteReferencePositions_[0];
99 siteReferencePositions_ = (Q & siteReferencePositions_);
106 forAll(siteReferencePositions_, i)
108 centreOfMass += siteReferencePositions_[i]*siteMasses_[i];
111 centreOfMass /= mass_;
113 siteReferencePositions_ -= centreOfMass;
117 forAll(siteReferencePositions_, i)
119 const vector&
p(siteReferencePositions_[i]);
140 forAll(siteReferencePositions_, i)
142 const vector&
p(siteReferencePositions_[i]);
144 momOfInertia += siteMasses_[i]*
tensor
146 p.
y()*p.
y() + p.
z()*p.
z(), -p.
x()*p.
y(), -p.
x()*p.
z(),
147 -p.
y()*p.
x(), p.
x()*p.
x() + p.
z()*p.
z(), -p.
y()*p.
z(),
148 -p.
z()*p.
x(), -p.
z()*p.
y(), p.
x()*p.
x() + p.
y()*p.
y()
154 FatalErrorIn(
"molecule::constantProperties::constantProperties")
155 <<
"An eigenvalue of the inertia tensor is zero, the molecule "
156 << dict.
name() <<
" is not a valid 6DOF shape."
175 siteReferencePositions_ = (Q & siteReferencePositions_);
184 forAll(siteReferencePositions_, i)
186 centreOfMass += siteReferencePositions_[i]*siteMasses_[i];
189 centreOfMass /= mass_;
191 siteReferencePositions_ -= centreOfMass;
198 forAll(siteReferencePositions_, i)
200 const vector&
p(siteReferencePositions_[i]);
202 momOfInertia += siteMasses_[i]*
tensor
204 p.
y()*p.
y() + p.
z()*p.
z(), -p.
x()*p.
y(), -p.
x()*p.
z(),
205 -p.
y()*p.
x(), p.
x()*p.
x() + p.
z()*p.
z(), -p.
y()*p.
z(),
206 -p.
z()*p.
x(), -p.
z()*p.
y(), p.
x()*p.
x() + p.
y()*p.
y()
243 specialPosition_(specialPosition),
244 potentialEnergy_(0.0),
249 sitePositions_(constProps.
nSites())
258 inline void Foam::molecule::constantProperties::checkSiteListSizes()
const
262 siteIds_.size() != siteReferencePositions_.size()
263 || siteIds_.size() != siteCharges_.size()
266 FatalErrorIn(
"molecule::constantProperties::checkSiteListSizes")
267 <<
"Sizes of site id, charge and "
268 <<
"referencePositions are not the same. "
274 inline void Foam::molecule::constantProperties::setInteracionSiteBools
276 const List<word>& siteIds,
277 const List<word>& pairPotSiteIds
280 pairPotentialSites_.setSize(siteIds_.size());
282 electrostaticSites_.setSize(siteIds_.size());
286 const word&
id(siteIds[i]);
288 pairPotentialSites_[i] = (
findIndex(pairPotSiteIds,
id) > -1);
290 electrostaticSites_[i] = (
mag(siteCharges_[i]) > VSMALL);
295 inline bool Foam::molecule::constantProperties::linearMoleculeTest()
const
297 if (siteIds_.size() == 2)
302 vector refDir = siteReferencePositions_[1] - siteReferencePositions_[0];
304 refDir /=
mag(refDir);
309 i < siteReferencePositions_.size();
313 vector dir = siteReferencePositions_[i] - siteReferencePositions_[i-1];
317 if (
mag(refDir & dir) < 1 - SMALL)
332 return siteReferencePositions_;
367 return pairPotentialSites_;
381 << sId <<
" site not found."
385 return pairPotentialSites_[s];
392 return electrostaticSites_;
406 "molecule::constantProperties::electrostaticSite(label)"
407 ) << sId <<
" site not found."
411 return electrostaticSites_[s];
418 return momentOfInertia_;
424 return ((momentOfInertia_.xx() < 0) && (momentOfInertia_.yy() > 0));
430 return (momentOfInertia_.zz() < 0);
436 if (linearMolecule())
440 else if (pointMolecule())
459 return siteIds_.size();
553 return sitePositions_;
559 return sitePositions_;
565 return specialPosition_;
571 return specialPosition_;
577 return potentialEnergy_;
583 return potentialEnergy_;