28 inline void Foam::moleculeCloud::evaluatePair
38 label idI = molI->id();
40 label idJ = molJ->id();
42 const molecule::constantProperties& constPropI(
constProps(idI));
44 const molecule::constantProperties& constPropJ(
constProps(idJ));
46 List<label> siteIdsI = constPropI.siteIds();
48 List<label> siteIdsJ = constPropJ.siteIds();
50 List<bool> pairPotentialSitesI = constPropI.pairPotentialSites();
52 List<bool> electrostaticSitesI = constPropI.electrostaticSites();
54 List<bool> pairPotentialSitesJ = constPropJ.pairPotentialSites();
56 List<bool> electrostaticSitesJ = constPropJ.electrostaticSites();
60 label idsI(siteIdsI[sI]);
64 label idsJ(siteIdsJ[sJ]);
66 if (pairPotentialSitesI[sI] && pairPotentialSitesJ[sJ])
69 molI->sitePositions()[sI] - molJ->sitePositions()[sJ];
71 scalar rsIsJMagSq =
magSqr(rsIsJ);
73 if(pairPot.rCutSqr(idsI, idsJ, rsIsJMagSq))
75 scalar rsIsJMag =
mag(rsIsJ);
79 *pairPot.force(idsI, idsJ, rsIsJMag);
81 molI->siteForces()[sI] += fsIsJ;
83 molJ->siteForces()[sJ] += -fsIsJ;
85 scalar potentialEnergy
87 pairPot.energy(idsI, idsJ, rsIsJMag)
90 molI->potentialEnergy() += 0.5*potentialEnergy;
92 molJ->potentialEnergy() += 0.5*potentialEnergy;
94 vector rIJ = molI->position() - molJ->position();
96 tensor virialContribution =
97 (rsIsJ*fsIsJ)*(rsIsJ & rIJ)/rsIsJMagSq;
99 molI->rf() += virialContribution;
101 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;
159 inline void Foam::moleculeCloud::evaluatePair
162 referredMolecule* molRef
165 const pairPotentialList& pairPot = pot_.pairPotentials();
167 const pairPotential& electrostatic = pairPot.electrostatic();
169 label idReal = molReal->id();
171 label idRef = molRef->id();
173 const molecule::constantProperties& constPropReal(constProps(idReal));
175 const molecule::constantProperties& constPropRef(constProps(idRef));
177 List<label> siteIdsReal = constPropReal.siteIds();
179 List<label> siteIdsRef = constPropRef.siteIds();
181 List<bool> pairPotentialSitesReal = constPropReal.pairPotentialSites();
183 List<bool> electrostaticSitesReal = constPropReal.electrostaticSites();
185 List<bool> pairPotentialSitesRef = constPropRef.pairPotentialSites();
187 List<bool> electrostaticSitesRef = constPropRef.electrostaticSites();
189 forAll(siteIdsReal, sReal)
191 label idsReal(siteIdsReal[sReal]);
195 label idsRef(siteIdsRef[sRef]);
197 if (pairPotentialSitesReal[sReal] && pairPotentialSitesRef[sRef])
200 molReal->sitePositions()[sReal]
201 - molRef->sitePositions()[sRef];
203 scalar rsRealsRefMagSq =
magSqr(rsRealsRef);
205 if (pairPot.rCutSqr(idsReal, idsRef, rsRealsRefMagSq))
207 scalar rsRealsRefMag =
mag(rsRealsRef);
210 (rsRealsRef/rsRealsRefMag)
211 *pairPot.force(idsReal, idsRef, rsRealsRefMag);
213 molReal->siteForces()[sReal] += fsRealsRef;
215 scalar potentialEnergy
217 pairPot.energy(idsReal, idsRef, rsRealsRefMag)
220 molReal->potentialEnergy() += 0.5*potentialEnergy;
222 vector rRealRef = molReal->position() - molRef->position();
225 (rsRealsRef*fsRealsRef)
226 *(rsRealsRef & rRealRef)
234 if (electrostaticSitesReal[sReal] && electrostaticSitesRef[sRef])
237 molReal->sitePositions()[sReal]
238 - molRef->sitePositions()[sRef];
240 scalar rsRealsRefMagSq =
magSqr(rsRealsRef);
242 if (rsRealsRefMagSq <= electrostatic.rCutSqr())
244 scalar rsRealsRefMag =
mag(rsRealsRef);
246 scalar chargeReal = constPropReal.siteCharges()[sReal];
248 scalar chargeRef = constPropRef.siteCharges()[sRef];
251 (rsRealsRef/rsRealsRefMag)
252 *chargeReal*chargeRef
253 *electrostatic.force(rsRealsRefMag);
255 molReal->siteForces()[sReal] += fsRealsRef;
257 scalar potentialEnergy =
259 *electrostatic.energy(rsRealsRefMag);
261 molReal->potentialEnergy() += 0.5*potentialEnergy;
263 vector rRealRef = molReal->position() - molRef->position();
266 (rsRealsRef*fsRealsRef)
267 *(rsRealsRef & rRealRef)
278 inline bool Foam::moleculeCloud::evaluatePotentialLimit
284 const pairPotentialList& pairPot = pot_.pairPotentials();
286 const pairPotential& electrostatic = pairPot.electrostatic();
288 label idI = molI->id();
290 label idJ = molJ->id();
292 const molecule::constantProperties& constPropI(constProps(idI));
294 const molecule::constantProperties& constPropJ(constProps(idJ));
296 List<label> siteIdsI = constPropI.siteIds();
298 List<label> siteIdsJ = constPropJ.siteIds();
300 List<bool> pairPotentialSitesI = constPropI.pairPotentialSites();
302 List<bool> electrostaticSitesI = constPropI.electrostaticSites();
304 List<bool> pairPotentialSitesJ = constPropJ.pairPotentialSites();
306 List<bool> electrostaticSitesJ = constPropJ.electrostaticSites();
310 label idsI(siteIdsI[sI]);
314 label idsJ(siteIdsJ[sJ]);
316 if (pairPotentialSitesI[sI] && pairPotentialSitesJ[sJ])
319 molI->sitePositions()[sI] - molJ->sitePositions()[sJ];
321 scalar rsIsJMagSq =
magSqr(rsIsJ);
323 if (pairPot.rCutSqr(idsI, idsJ, rsIsJMagSq))
325 scalar rsIsJMag =
mag(rsIsJ);
331 if (rsIsJMag < SMALL)
333 WarningIn(
"moleculeCloud::removeHighEnergyOverlaps()")
334 <<
"Molecule site pair closer than "
336 <<
": mag separation = " << rsIsJMag
337 <<
". These may have been placed on top of each"
338 <<
" other by a rounding error in mdInitialise in"
339 <<
" parallel or a block filled with molecules"
340 <<
" twice. Removing one of the molecules."
349 if (rsIsJMag < pairPot.rMin(idsI, idsJ))
356 mag(pairPot.energy(idsI, idsJ, rsIsJMag))
357 > pot_.potentialEnergyLimit()
365 if (electrostaticSitesI[sI] && electrostaticSitesJ[sJ])
368 molI->sitePositions()[sI] - molJ->sitePositions()[sJ];
370 scalar rsIsJMagSq =
magSqr(rsIsJ);
372 if (pairPot.rCutMaxSqr(rsIsJMagSq))
374 scalar rsIsJMag =
mag(rsIsJ);
380 if (rsIsJMag < SMALL)
382 WarningIn(
"moleculeCloud::removeHighEnergyOverlaps()")
383 <<
"Molecule site pair closer than "
385 <<
": mag separation = " << rsIsJMag
386 <<
". These may have been placed on top of each"
387 <<
" other by a rounding error in mdInitialise in"
388 <<
" parallel or a block filled with molecules"
389 <<
" twice. Removing one of the molecules."
395 if (rsIsJMag < electrostatic.rMin())
400 scalar chargeI = constPropI.siteCharges()[sI];
402 scalar chargeJ = constPropJ.siteCharges()[sJ];
406 mag(chargeI*chargeJ*electrostatic.energy(rsIsJMag))
407 > pot_.potentialEnergyLimit()
421 inline bool Foam::moleculeCloud::evaluatePotentialLimit
424 referredMolecule* molRef
427 const pairPotentialList& pairPot = pot_.pairPotentials();
429 const pairPotential& electrostatic = pairPot.electrostatic();
431 label idReal = molReal->id();
433 label idRef = molRef->id();
435 const molecule::constantProperties& constPropReal(constProps(idReal));
437 const molecule::constantProperties& constPropRef(constProps(idRef));
439 List<label> siteIdsReal = constPropReal.siteIds();
441 List<label> siteIdsRef = constPropRef.siteIds();
443 List<bool> pairPotentialSitesReal = constPropReal.pairPotentialSites();
445 List<bool> electrostaticSitesReal = constPropReal.electrostaticSites();
447 List<bool> pairPotentialSitesRef = constPropRef.pairPotentialSites();
449 List<bool> electrostaticSitesRef = constPropRef.electrostaticSites();
451 forAll(siteIdsReal, sReal)
453 label idsReal(siteIdsReal[sReal]);
457 label idsRef(siteIdsRef[sRef]);
459 if (pairPotentialSitesReal[sReal] && pairPotentialSitesRef[sRef])
462 molReal->sitePositions()[sReal]
463 - molRef->sitePositions()[sRef];
465 scalar rsRealsRefMagSq =
magSqr(rsRealsRef);
467 if (pairPot.rCutSqr(idsReal, idsRef, rsRealsRefMagSq))
469 scalar rsRealsRefMag =
mag(rsRealsRef);
475 if (rsRealsRefMag < SMALL)
477 WarningIn(
"moleculeCloud::removeHighEnergyOverlaps()")
478 <<
"Molecule site pair closer than "
480 <<
": mag separation = " << rsRealsRefMag
481 <<
". These may have been placed on top of each"
482 <<
" other by a rounding error in mdInitialise in"
483 <<
" parallel or a block filled with molecules"
484 <<
" twice. Removing one of the molecules."
494 if (rsRealsRefMag < pairPot.rMin(idsReal, idsRef))
501 mag(pairPot.energy(idsReal, idsRef, rsRealsRefMag))
502 > pot_.potentialEnergyLimit()
510 if (electrostaticSitesReal[sReal] && electrostaticSitesRef[sRef])
513 molReal->sitePositions()[sReal]
514 - molRef->sitePositions()[sRef];
516 scalar rsRealsRefMagSq =
magSqr(rsRealsRef);
518 if (pairPot.rCutMaxSqr(rsRealsRefMagSq))
520 scalar rsRealsRefMag =
mag(rsRealsRef);
526 if (rsRealsRefMag < SMALL)
528 WarningIn(
"moleculeCloud::removeHighEnergyOverlaps()")
529 <<
"Molecule site pair closer than "
531 <<
": mag separation = " << rsRealsRefMag
532 <<
". These may have been placed on top of each"
533 <<
" other by a rounding error in mdInitialise in"
534 <<
" parallel or a block filled with molecules"
535 <<
" twice. Removing one of the molecules."
541 if (rsRealsRefMag < electrostatic.rMin())
546 scalar chargeReal = constPropReal.siteCharges()[sReal];
548 scalar chargeRef = constPropRef.siteCharges()[sRef];
556 *electrostatic.energy(rsRealsRefMag)
558 > pot_.potentialEnergyLimit()
572 inline Foam::vector Foam::moleculeCloud::equipartitionLinearVelocity
580 rndGen_.GaussNormal(),
581 rndGen_.GaussNormal(),
582 rndGen_.GaussNormal()
587 inline Foam::vector Foam::moleculeCloud::equipartitionAngularMomentum
590 const molecule::constantProperties& cP
593 scalar sqrtKbT =
sqrt(kb*temperature);
595 if (cP.linearMolecule())
600 sqrt(cP.momentOfInertia().yy())*rndGen_.GaussNormal(),
601 sqrt(cP.momentOfInertia().zz())*rndGen_.GaussNormal()
608 sqrt(cP.momentOfInertia().xx())*rndGen_.GaussNormal(),
609 sqrt(cP.momentOfInertia().yy())*rndGen_.GaussNormal(),
610 sqrt(cP.momentOfInertia().zz())*rndGen_.GaussNormal()
633 return cellOccupancy_;
647 return constPropList_;
654 return constPropList_[id];