33 siteReferencePositions_(),
37 pairPotentialSites_(),
38 electrostaticSites_(),
39 momentOfInertia_(
Zero),
49 siteReferencePositions_(
dict.
lookup(
"siteReferencePositions")),
53 pairPotentialSites_(),
54 electrostaticSites_(),
60 setInteracionSiteBools
66 mass_ =
sum(siteMasses_);
73 forAll(siteReferencePositions_, i)
75 centreOfMass += siteReferencePositions_[i]*siteMasses_[i];
78 centreOfMass /= mass_;
80 siteReferencePositions_ -= centreOfMass;
82 if (siteIds_.
size() == 1)
86 siteReferencePositions_[0] =
Zero;
90 else if (linearMoleculeTest())
99 siteReferencePositions_[1] - siteReferencePositions_[0]
104 siteReferencePositions_ = (
Q & siteReferencePositions_);
111 forAll(siteReferencePositions_, i)
113 centreOfMass += siteReferencePositions_[i]*siteMasses_[i];
116 centreOfMass /= mass_;
118 siteReferencePositions_ -= centreOfMass;
122 forAll(siteReferencePositions_, i)
124 const vector&
p(siteReferencePositions_[i]);
145 forAll(siteReferencePositions_, i)
147 const vector&
p(siteReferencePositions_[i]);
151 p.y()*
p.y() +
p.z()*
p.z(), -
p.x()*
p.y(), -
p.x()*
p.z(),
152 p.x()*
p.x() +
p.z()*
p.z(), -
p.y()*
p.z(),
153 p.x()*
p.x() +
p.y()*
p.y()
160 <<
"An eigenvalue of the inertia tensor is zero, the molecule " 161 <<
dict.name() <<
" is not a valid 6DOF shape." 180 siteReferencePositions_ = (
Q & siteReferencePositions_);
189 forAll(siteReferencePositions_, i)
191 centreOfMass += siteReferencePositions_[i]*siteMasses_[i];
194 centreOfMass /= mass_;
196 siteReferencePositions_ -= centreOfMass;
203 forAll(siteReferencePositions_, i)
205 const vector&
p(siteReferencePositions_[i]);
209 p.y()*
p.y() +
p.z()*
p.z(), -
p.x()*
p.y(), -
p.x()*
p.z(),
210 p.x()*
p.x() +
p.z()*
p.z(), -
p.y()*
p.z(),
211 p.x()*
p.x() +
p.y()*
p.y()
227 const polyMesh&
mesh,
230 const label tetFacei,
238 const constantProperties& constProps,
251 potentialEnergy_(0.0),
255 siteForces_(constProps.nSites(),
Zero),
256 sitePositions_(constProps.nSites())
273 const constantProperties& constProps,
286 potentialEnergy_(0.0),
290 siteForces_(constProps.nSites(),
Zero),
291 sitePositions_(constProps.nSites())
299 inline void Foam::molecule::constantProperties::checkSiteListSizes()
const 303 siteIds_.size() != siteReferencePositions_.size()
304 || siteIds_.size() != siteCharges_.size()
308 <<
"Sizes of site id, charge and " 309 <<
"referencePositions are not the same. " 315 inline void Foam::molecule::constantProperties::setInteracionSiteBools
317 const List<word>& siteIds,
318 const List<word>& pairPotSiteIds
321 pairPotentialSites_.setSize(siteIds_.size());
323 electrostaticSites_.setSize(siteIds_.size());
327 const word&
id = siteIds[i];
329 pairPotentialSites_[i] = pairPotSiteIds.found(
id);
330 electrostaticSites_[i] = (
mag(siteCharges_[i]) > VSMALL);
335 inline bool Foam::molecule::constantProperties::linearMoleculeTest()
const 337 if (siteIds_.size() == 2)
345 siteReferencePositions_[1] - siteReferencePositions_[0]
351 i < siteReferencePositions_.size();
358 siteReferencePositions_[i] - siteReferencePositions_[i-1]
361 if (
mag(refDir & dir) < 1 - SMALL)
376 return siteReferencePositions_;
411 return pairPotentialSites_;
420 label
s = siteIds_.find(sId);
425 << sId <<
" site not found." 429 return pairPotentialSites_[
s];
436 return electrostaticSites_;
445 label
s = siteIds_.find(sId);
450 << sId <<
" site not found." 454 return electrostaticSites_[
s];
461 return momentOfInertia_;
467 return ((momentOfInertia_.xx() < 0) && (momentOfInertia_.yy() > 0));
473 return (momentOfInertia_.zz() < 0);
479 if (linearMolecule())
483 else if (pointMolecule())
502 return siteIds_.size();
583 return sitePositions_;
589 return sitePositions_;
595 return specialPosition_;
601 return specialPosition_;
607 return potentialEnergy_;
613 return potentialEnergy_;
637 return special_ == SPECIAL_TETHERED;
scalar potentialEnergy() const
const List< bool > & electrostaticSites() const
void size(const label n)
Older name for setAddressableSize.
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &f1)
const List< vector > & sitePositions() const
DiagTensor< scalar > diagTensor
DiagTensor of scalars, i.e. DiagTensor<scalar>.
const tensor & rf() const
const diagTensor & momentOfInertia() const
void setSitePositions(const constantProperties &constProps)
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
const barycentric & coordinates() const noexcept
Return current particle coordinates.
error FatalError
Error stream (stdout output on all processes), with additional 'FOAM FATAL ERROR' header text and sta...
A list of keyword definitions, which are a keyword followed by a number of values (eg...
const List< vector > & siteForces() const
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
A 1D array of objects of type <T>, where the size of the vector is known and used for subscript bound...
constexpr char nl
The newline '\n' character (0x0a)
const List< label > & siteIds() const
Ostream & endl(Ostream &os)
Add newline and flush stream.
const vector & pi() const
Barycentric< scalar > barycentric
A scalar version of the templated Barycentric.
quaternion normalised(const quaternion &q)
Return the normalised (unit) quaternion of the given quaternion.
tensor rotationTensor(const vector &n1, const vector &n2)
Rotational transformation tensor from vector n1 to n2.
Lookup type of boundary radiation properties.
dimensionedTensor eigenVectors(const dimensionedSymmTensor &dt)
#define forAll(list, i)
Loop across all elements in list.
const Field< vector > & siteReferencePositions() const
const dimensionedScalar e
Elementary charge.
molecule(const polyMesh &mesh, const barycentric &coordinates, const label celli, const label tetFacei, const label tetPti, const tensor &Q, const vector &v, const vector &a, const vector &pi, const vector &tau, const vector &specialPosition, const constantProperties &constProps, const label special, const label id)
Construct from components.
SymmTensor< scalar > symmTensor
SymmTensor of scalars, i.e. SymmTensor<scalar>.
A class for handling words, derived from Foam::string.
particle(const polyMesh &mesh, const barycentric &coordinates, const label celli, const label tetFacei, const label tetPti)
Construct from components.
errorManip< error > abort(error &err)
const List< bool > & pairPotentialSites() const
A Vector of values with scalar precision, where scalar is float/double depending on the compilation f...
bool pointMolecule() const
bool electrostaticSite(label sId) const
const polyMesh & mesh() const noexcept
Return the mesh database.
label degreesOfFreedom() const
messageStream Info
Information stream (stdout output on master, null elsewhere)
dimensionedVector eigenValues(const dimensionedSymmTensor &dt)
Mesh consisting of general polyhedral cells.
bool linearMolecule() const
bool pairPotentialSite(label sId) const
gmvFile<< "tracers "<< particles.size()<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().x()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().y()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().z()<< " ";}gmvFile<< nl;forAll(lagrangianScalarNames, i){ word name=lagrangianScalarNames[i];IOField< scalar > s(IOobject(name, runTime.timeName(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
Tensor of scalars, i.e. Tensor<scalar>.
const List< scalar > & siteMasses() const
const vector & tau() const
vector position() const
Return current particle position.
const List< scalar > & siteCharges() const
const vector & specialPosition() const
static constexpr const zero Zero
Global zero (0)