moleculeCloudI.H
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2011-2024 OpenFOAM Foundation
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8 License
9  This file is part of OpenFOAM.
10 
11  OpenFOAM is free software: you can redistribute it and/or modify it
12  under the terms of the GNU General Public License as published by
13  the Free Software Foundation, either version 3 of the License, or
14  (at your option) any later version.
15 
16  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
19  for more details.
20 
21  You should have received a copy of the GNU General Public License
22  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
23 
24 \*---------------------------------------------------------------------------*/
25 
26 #include "constants.H"
27 
28 using namespace Foam::constant;
29 
30 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
31 
32 inline void Foam::moleculeCloud::evaluatePair
33 (
34  molecule& molI,
35  molecule& molJ
36 )
37 {
38  const pairPotentialList& pairPot = pot_.pairPotentials();
39 
40  const pairPotential& electrostatic = pairPot.electrostatic();
41 
42  label idI = molI.id();
43 
44  label idJ = molJ.id();
45 
46  const molecule::constantProperties& constPropI(constProps(idI));
47 
48  const molecule::constantProperties& constPropJ(constProps(idJ));
49 
50  List<label> siteIdsI = constPropI.siteIds();
51 
52  List<label> siteIdsJ = constPropJ.siteIds();
53 
54  List<bool> pairPotentialSitesI = constPropI.pairPotentialSites();
55 
56  List<bool> electrostaticSitesI = constPropI.electrostaticSites();
57 
58  List<bool> pairPotentialSitesJ = constPropJ.pairPotentialSites();
59 
60  List<bool> electrostaticSitesJ = constPropJ.electrostaticSites();
61 
62  forAll(siteIdsI, sI)
63  {
64  label idsI(siteIdsI[sI]);
65 
66  forAll(siteIdsJ, sJ)
67  {
68  label idsJ(siteIdsJ[sJ]);
69 
70  if (pairPotentialSitesI[sI] && pairPotentialSitesJ[sJ])
71  {
72  vector rsIsJ =
73  molI.sitePositions()[sI] - molJ.sitePositions()[sJ];
74 
75  scalar rsIsJMagSq = magSqr(rsIsJ);
76 
77  if (pairPot.rCutSqr(idsI, idsJ, rsIsJMagSq))
78  {
79  scalar rsIsJMag = mag(rsIsJ);
80 
81  vector fsIsJ =
82  (rsIsJ/rsIsJMag)
83  *pairPot.force(idsI, idsJ, rsIsJMag);
84 
85  molI.siteForces()[sI] += fsIsJ;
86 
87  molJ.siteForces()[sJ] += -fsIsJ;
88 
89  scalar potentialEnergy
90  (
91  pairPot.energy(idsI, idsJ, rsIsJMag)
92  );
93 
94  molI.potentialEnergy() += 0.5*potentialEnergy;
95 
96  molJ.potentialEnergy() += 0.5*potentialEnergy;
97 
98  vector rIJ = molI.position(mesh()) - molJ.position(mesh());
99 
100  tensor virialContribution =
101  (rsIsJ*fsIsJ)*(rsIsJ & rIJ)/rsIsJMagSq;
102 
103  molI.rf() += virialContribution;
104 
105  molJ.rf() += virialContribution;
106  }
107  }
108 
109  if (electrostaticSitesI[sI] && electrostaticSitesJ[sJ])
110  {
111  vector rsIsJ =
112  molI.sitePositions()[sI] - molJ.sitePositions()[sJ];
113 
114  scalar rsIsJMagSq = magSqr(rsIsJ);
115 
116  if (rsIsJMagSq <= electrostatic.rCutSqr())
117  {
118  scalar rsIsJMag = mag(rsIsJ);
119 
120  scalar chargeI = constPropI.siteCharges()[sI];
121 
122  scalar chargeJ = constPropJ.siteCharges()[sJ];
123 
124  vector fsIsJ =
125  (rsIsJ/rsIsJMag)
126  *chargeI*chargeJ*electrostatic.force(rsIsJMag);
127 
128  molI.siteForces()[sI] += fsIsJ;
129 
130  molJ.siteForces()[sJ] += -fsIsJ;
131 
132  scalar potentialEnergy =
133  chargeI*chargeJ
134  *electrostatic.energy(rsIsJMag);
135 
136  molI.potentialEnergy() += 0.5*potentialEnergy;
137 
138  molJ.potentialEnergy() += 0.5*potentialEnergy;
139 
140  vector rIJ = molI.position(mesh()) - molJ.position(mesh());
141 
142  tensor virialContribution =
143  (rsIsJ*fsIsJ)*(rsIsJ & rIJ)/rsIsJMagSq;
144 
145  molI.rf() += virialContribution;
146 
147  molJ.rf() += virialContribution;
148  }
149  }
150  }
151  }
152 }
153 
154 
155 inline bool Foam::moleculeCloud::evaluatePotentialLimit
156 (
157  molecule& molI,
158  molecule& molJ
159 ) const
160 {
161  const pairPotentialList& pairPot = pot_.pairPotentials();
162 
163  const pairPotential& electrostatic = pairPot.electrostatic();
164 
165  label idI = molI.id();
166 
167  label idJ = molJ.id();
168 
169  const molecule::constantProperties& constPropI(constProps(idI));
170 
171  const molecule::constantProperties& constPropJ(constProps(idJ));
172 
173  List<label> siteIdsI = constPropI.siteIds();
174 
175  List<label> siteIdsJ = constPropJ.siteIds();
176 
177  List<bool> pairPotentialSitesI = constPropI.pairPotentialSites();
178 
179  List<bool> electrostaticSitesI = constPropI.electrostaticSites();
180 
181  List<bool> pairPotentialSitesJ = constPropJ.pairPotentialSites();
182 
183  List<bool> electrostaticSitesJ = constPropJ.electrostaticSites();
184 
185  forAll(siteIdsI, sI)
186  {
187  label idsI(siteIdsI[sI]);
188 
189  forAll(siteIdsJ, sJ)
190  {
191  label idsJ(siteIdsJ[sJ]);
192 
193  if (pairPotentialSitesI[sI] && pairPotentialSitesJ[sJ])
194  {
195  vector rsIsJ =
196  molI.sitePositions()[sI] - molJ.sitePositions()[sJ];
197 
198  scalar rsIsJMagSq = magSqr(rsIsJ);
199 
200  if (pairPot.rCutSqr(idsI, idsJ, rsIsJMagSq))
201  {
202  scalar rsIsJMag = mag(rsIsJ);
203 
204  // Guard against pairPot.energy being evaluated
205  // if rIJMag < small. A floating point exception will
206  // happen otherwise.
207 
208  if (rsIsJMag < small)
209  {
211  << "Molecule site pair closer than "
212  << small
213  << ": mag separation = " << rsIsJMag
214  << ". These may have been placed on top of each"
215  << " other by a rounding error in mdInitialise in"
216  << " parallel or a block filled with molecules"
217  << " twice. Removing one of the molecules."
218  << endl;
219 
220  return true;
221  }
222 
223  // Guard against pairPot.energy being evaluated if rIJMag <
224  // rMin. A tabulation lookup error will occur otherwise.
225 
226  if (rsIsJMag < pairPot.rMin(idsI, idsJ))
227  {
228  return true;
229  }
230 
231  if
232  (
233  mag(pairPot.energy(idsI, idsJ, rsIsJMag))
234  > pot_.potentialEnergyLimit()
235  )
236  {
237  return true;
238  };
239  }
240  }
241 
242  if (electrostaticSitesI[sI] && electrostaticSitesJ[sJ])
243  {
244  vector rsIsJ =
245  molI.sitePositions()[sI] - molJ.sitePositions()[sJ];
246 
247  scalar rsIsJMagSq = magSqr(rsIsJ);
248 
249  if (pairPot.rCutMaxSqr(rsIsJMagSq))
250  {
251  scalar rsIsJMag = mag(rsIsJ);
252 
253  // Guard against pairPot.energy being evaluated
254  // if rIJMag < small. A floating point exception will
255  // happen otherwise.
256 
257  if (rsIsJMag < small)
258  {
260  << "Molecule site pair closer than "
261  << small
262  << ": mag separation = " << rsIsJMag
263  << ". These may have been placed on top of each"
264  << " other by a rounding error in mdInitialise in"
265  << " parallel or a block filled with molecules"
266  << " twice. Removing one of the molecules."
267  << endl;
268 
269  return true;
270  }
271 
272  if (rsIsJMag < electrostatic.rMin())
273  {
274  return true;
275  }
276 
277  scalar chargeI = constPropI.siteCharges()[sI];
278 
279  scalar chargeJ = constPropJ.siteCharges()[sJ];
280 
281  if
282  (
283  mag(chargeI*chargeJ*electrostatic.energy(rsIsJMag))
284  > pot_.potentialEnergyLimit()
285  )
286  {
287  return true;
288  };
289  }
290  }
291  }
292  }
293 
294  return false;
295 }
296 
297 
298 inline Foam::vector Foam::moleculeCloud::equipartitionLinearVelocity
299 (
300  scalar temperature,
301  scalar mass
302 )
303 {
304  return
305  sqrt(physicoChemical::k.value()*temperature/mass)
306  *stdNormal_.sample<vector>();
307 }
308 
309 
310 inline Foam::vector Foam::moleculeCloud::equipartitionAngularMomentum
311 (
312  scalar temperature,
313  const molecule::constantProperties& cP
314 )
315 {
316  scalar sqrtKbT = sqrt(physicoChemical::k.value()*temperature);
317 
318  if (cP.linearMolecule())
319  {
320  return sqrtKbT*vector
321  (
322  0.0,
323  sqrt(cP.momentOfInertia().yy())*stdNormal_.sample(),
324  sqrt(cP.momentOfInertia().zz())*stdNormal_.sample()
325  );
326  }
327  else
328  {
329  return sqrtKbT*vector
330  (
331  sqrt(cP.momentOfInertia().xx())*stdNormal_.sample(),
332  sqrt(cP.momentOfInertia().yy())*stdNormal_.sample(),
333  sqrt(cP.momentOfInertia().zz())*stdNormal_.sample()
334  );
335  }
336 }
337 
338 
339 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
340 
342 {
343  return mesh_;
344 }
345 
346 
348 {
349  return pot_;
350 }
351 
352 
355 {
356  return cellOccupancy_;
357 }
358 
359 
362 {
363  return il_;
364 }
365 
366 
369 {
370  return constPropList_;
371 }
372 
373 
376 {
377  return constPropList_[id];
378 }
379 
380 
382 {
383  return rndGen_;
384 }
385 
386 
388 {
389  return stdNormal_;
390 }
391 
392 
393 // ************************************************************************* //
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
A 1D array of objects of type <T>, where the size of the vector is known and used for subscript bound...
Definition: List.H:91
Standard normal distribution. Not selectable.
const potential & pot() const
const List< molecule::constantProperties > constProps() const
const List< DynamicList< molecule * > > & cellOccupancy() const
const polyMesh & mesh() const
distributions::standardNormal & stdNormal()
const InteractionLists< molecule > & il() const
randomGenerator & rndGen()
Class to hold molecule constant properties.
Definition: molecule.H:91
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:80
Random number generator.
#define WarningInFunction
Report a warning using Foam::Warning.
const dimensionedScalar k
Boltzmann constant.
Collection of constants.
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Definition: label.H:59
Tensor< scalar > tensor
Tensor of scalars.
Definition: tensor.H:51
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:257
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
dimensionedScalar sqrt(const dimensionedScalar &ds)
dimensioned< scalar > mag(const dimensioned< Type > &)
dimensioned< scalar > magSqr(const dimensioned< Type > &)