34 inline void Foam::moleculeCloud::evaluatePair
40 const pairPotentialList& pairPot = pot_.pairPotentials();
42 const pairPotential& electrostatic = pairPot.electrostatic();
44 label idI = molI.id();
46 label idJ = molJ.id();
48 const molecule::constantProperties& constPropI(constProps(idI));
50 const molecule::constantProperties& constPropJ(constProps(idJ));
52 List<label> siteIdsI = constPropI.siteIds();
54 List<label> siteIdsJ = constPropJ.siteIds();
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])
75 molI.sitePositions()[sI] - molJ.sitePositions()[sJ];
77 scalar rsIsJMagSq =
magSqr(rsIsJ);
79 if (pairPot.rCutSqr(idsI, idsJ, rsIsJMagSq))
81 scalar rsIsJMag =
mag(rsIsJ);
85 *pairPot.force(idsI, idsJ, rsIsJMag);
87 molI.siteForces()[sI] += fsIsJ;
89 molJ.siteForces()[sJ] += -fsIsJ;
91 scalar potentialEnergy
93 pairPot.energy(idsI, idsJ, rsIsJMag)
96 molI.potentialEnergy() += 0.5*potentialEnergy;
98 molJ.potentialEnergy() += 0.5*potentialEnergy;
100 vector rIJ = molI.position() - molJ.position();
102 tensor virialContribution =
103 (rsIsJ*fsIsJ)*(rsIsJ & rIJ)/rsIsJMagSq;
105 molI.rf() += virialContribution;
107 molJ.rf() += virialContribution;
111 if (electrostaticSitesI[sI] && electrostaticSitesJ[sJ])
114 molI.sitePositions()[sI] - molJ.sitePositions()[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);
130 molI.siteForces()[sI] += fsIsJ;
132 molJ.siteForces()[sJ] += -fsIsJ;
134 scalar potentialEnergy =
136 *electrostatic.energy(rsIsJMag);
138 molI.potentialEnergy() += 0.5*potentialEnergy;
140 molJ.potentialEnergy() += 0.5*potentialEnergy;
142 vector rIJ = molI.position() - molJ.position();
144 tensor virialContribution =
145 (rsIsJ*fsIsJ)*(rsIsJ & rIJ)/rsIsJMagSq;
147 molI.rf() += virialContribution;
149 molJ.rf() += virialContribution;
157 inline bool Foam::moleculeCloud::evaluatePotentialLimit
163 const pairPotentialList& pairPot = pot_.pairPotentials();
165 const pairPotential& electrostatic = pairPot.electrostatic();
167 label idI = molI.id();
169 label idJ = molJ.id();
171 const molecule::constantProperties& constPropI(constProps(idI));
173 const molecule::constantProperties& constPropJ(constProps(idJ));
175 List<label> siteIdsI = constPropI.siteIds();
177 List<label> siteIdsJ = constPropJ.siteIds();
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])
198 molI.sitePositions()[sI] - molJ.sitePositions()[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])
247 molI.sitePositions()[sI] - molJ.sitePositions()[sJ];
249 scalar rsIsJMagSq =
magSqr(rsIsJ);
251 if (pairPot.rCutMaxSqr(rsIsJMagSq))
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()
300 inline Foam::vector Foam::moleculeCloud::equipartitionLinearVelocity
308 *rndGen_.GaussNormal<
vector>();
312 inline Foam::vector Foam::moleculeCloud::equipartitionAngularMomentum
315 const molecule::constantProperties& cP
320 if (cP.linearMolecule())
325 sqrt(cP.momentOfInertia().yy())*rndGen_.GaussNormal<scalar>(),
326 sqrt(cP.momentOfInertia().zz())*rndGen_.GaussNormal<scalar>()
333 sqrt(cP.momentOfInertia().xx())*rndGen_.GaussNormal<scalar>(),
334 sqrt(cP.momentOfInertia().yy())*rndGen_.GaussNormal<scalar>(),
335 sqrt(cP.momentOfInertia().zz())*rndGen_.GaussNormal<scalar>()
358 return cellOccupancy_;
372 return constPropList_;
379 return constPropList_[id];
Different types of constants.
const InteractionLists< molecule > & il() const
Class to hold molecule constant properties.
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
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.
#define forAll(list, i)
Loop across all elements in list.
const List< DynamicList< molecule * > > & cellOccupancy() const
A Vector of values with scalar precision, where scalar is float/double depending on the compilation f...
const potential & pot() const
const dimensionedScalar k
Boltzmann constant.
#define WarningInFunction
Report a warning using Foam::Warning.
const List< molecule::constantProperties > constProps() const
Mesh consisting of general polyhedral cells.
const polyMesh & mesh() const
dimensioned< typename typeOfMag< Type >::type > magSqr(const dimensioned< Type > &dt)