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