FreeFOAM The Cross-Platform CFD Toolkit
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 | Copyright (C) 2008-2010 OpenCFD Ltd.
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 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
27 
28 inline void Foam::moleculeCloud::evaluatePair
29 (
30  molecule* molI,
31  molecule* molJ
32 )
33 {
34  const pairPotentialList& pairPot = pot_.pairPotentials();
35 
36  const pairPotential& electrostatic = pairPot.electrostatic();
37 
38  label idI = molI->id();
39 
40  label idJ = molJ->id();
41 
42  const molecule::constantProperties& constPropI(constProps(idI));
43 
44  const molecule::constantProperties& constPropJ(constProps(idJ));
45 
46  List<label> siteIdsI = constPropI.siteIds();
47 
48  List<label> siteIdsJ = constPropJ.siteIds();
49 
50  List<bool> pairPotentialSitesI = constPropI.pairPotentialSites();
51 
52  List<bool> electrostaticSitesI = constPropI.electrostaticSites();
53 
54  List<bool> pairPotentialSitesJ = constPropJ.pairPotentialSites();
55 
56  List<bool> electrostaticSitesJ = constPropJ.electrostaticSites();
57 
58  forAll(siteIdsI, sI)
59  {
60  label idsI(siteIdsI[sI]);
61 
62  forAll(siteIdsJ, sJ)
63  {
64  label idsJ(siteIdsJ[sJ]);
65 
66  if (pairPotentialSitesI[sI] && pairPotentialSitesJ[sJ])
67  {
68  vector rsIsJ =
69  molI->sitePositions()[sI] - molJ->sitePositions()[sJ];
70 
71  scalar rsIsJMagSq = magSqr(rsIsJ);
72 
73  if(pairPot.rCutSqr(idsI, idsJ, rsIsJMagSq))
74  {
75  scalar rsIsJMag = mag(rsIsJ);
76 
77  vector fsIsJ =
78  (rsIsJ/rsIsJMag)
79  *pairPot.force(idsI, idsJ, rsIsJMag);
80 
81  molI->siteForces()[sI] += fsIsJ;
82 
83  molJ->siteForces()[sJ] += -fsIsJ;
84 
85  scalar potentialEnergy
86  (
87  pairPot.energy(idsI, idsJ, rsIsJMag)
88  );
89 
90  molI->potentialEnergy() += 0.5*potentialEnergy;
91 
92  molJ->potentialEnergy() += 0.5*potentialEnergy;
93 
94  vector rIJ = molI->position() - molJ->position();
95 
96  tensor virialContribution =
97  (rsIsJ*fsIsJ)*(rsIsJ & rIJ)/rsIsJMagSq;
98 
99  molI->rf() += virialContribution;
100 
101  molJ->rf() += virialContribution;
102 
103  // molI->rf() += rsIsJ * fsIsJ;
104 
105  // molJ->rf() += rsIsJ * fsIsJ;
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() - molJ->position();
141 
142  tensor virialContribution =
143  (rsIsJ*fsIsJ)*(rsIsJ & rIJ)/rsIsJMagSq;
144 
145  molI->rf() += virialContribution;
146 
147  molJ->rf() += virialContribution;
148 
149  // molI->rf() += rsIsJ * fsIsJ;
150 
151  // molJ->rf() += rsIsJ * fsIsJ;
152  }
153  }
154  }
155  }
156 }
157 
158 
159 inline void Foam::moleculeCloud::evaluatePair
160 (
161  molecule* molReal,
162  referredMolecule* molRef
163 )
164 {
165  const pairPotentialList& pairPot = pot_.pairPotentials();
166 
167  const pairPotential& electrostatic = pairPot.electrostatic();
168 
169  label idReal = molReal->id();
170 
171  label idRef = molRef->id();
172 
173  const molecule::constantProperties& constPropReal(constProps(idReal));
174 
175  const molecule::constantProperties& constPropRef(constProps(idRef));
176 
177  List<label> siteIdsReal = constPropReal.siteIds();
178 
179  List<label> siteIdsRef = constPropRef.siteIds();
180 
181  List<bool> pairPotentialSitesReal = constPropReal.pairPotentialSites();
182 
183  List<bool> electrostaticSitesReal = constPropReal.electrostaticSites();
184 
185  List<bool> pairPotentialSitesRef = constPropRef.pairPotentialSites();
186 
187  List<bool> electrostaticSitesRef = constPropRef.electrostaticSites();
188 
189  forAll(siteIdsReal, sReal)
190  {
191  label idsReal(siteIdsReal[sReal]);
192 
193  forAll(siteIdsRef, sRef)
194  {
195  label idsRef(siteIdsRef[sRef]);
196 
197  if (pairPotentialSitesReal[sReal] && pairPotentialSitesRef[sRef])
198  {
199  vector rsRealsRef =
200  molReal->sitePositions()[sReal]
201  - molRef->sitePositions()[sRef];
202 
203  scalar rsRealsRefMagSq = magSqr(rsRealsRef);
204 
205  if (pairPot.rCutSqr(idsReal, idsRef, rsRealsRefMagSq))
206  {
207  scalar rsRealsRefMag = mag(rsRealsRef);
208 
209  vector fsRealsRef =
210  (rsRealsRef/rsRealsRefMag)
211  *pairPot.force(idsReal, idsRef, rsRealsRefMag);
212 
213  molReal->siteForces()[sReal] += fsRealsRef;
214 
215  scalar potentialEnergy
216  (
217  pairPot.energy(idsReal, idsRef, rsRealsRefMag)
218  );
219 
220  molReal->potentialEnergy() += 0.5*potentialEnergy;
221 
222  vector rRealRef = molReal->position() - molRef->position();
223 
224  molReal->rf() +=
225  (rsRealsRef*fsRealsRef)
226  *(rsRealsRef & rRealRef)
227  /rsRealsRefMagSq;
228 
229  // molReal->rf() += rsRealsRef * fsRealsRef;
230 
231  }
232  }
233 
234  if (electrostaticSitesReal[sReal] && electrostaticSitesRef[sRef])
235  {
236  vector rsRealsRef =
237  molReal->sitePositions()[sReal]
238  - molRef->sitePositions()[sRef];
239 
240  scalar rsRealsRefMagSq = magSqr(rsRealsRef);
241 
242  if (rsRealsRefMagSq <= electrostatic.rCutSqr())
243  {
244  scalar rsRealsRefMag = mag(rsRealsRef);
245 
246  scalar chargeReal = constPropReal.siteCharges()[sReal];
247 
248  scalar chargeRef = constPropRef.siteCharges()[sRef];
249 
250  vector fsRealsRef =
251  (rsRealsRef/rsRealsRefMag)
252  *chargeReal*chargeRef
253  *electrostatic.force(rsRealsRefMag);
254 
255  molReal->siteForces()[sReal] += fsRealsRef;
256 
257  scalar potentialEnergy =
258  chargeReal*chargeRef
259  *electrostatic.energy(rsRealsRefMag);
260 
261  molReal->potentialEnergy() += 0.5*potentialEnergy;
262 
263  vector rRealRef = molReal->position() - molRef->position();
264 
265  molReal->rf() +=
266  (rsRealsRef*fsRealsRef)
267  *(rsRealsRef & rRealRef)
268  /rsRealsRefMagSq;
269 
270  // molReal->rf() += rsRealsRef * fsRealsRef;
271  }
272  }
273  }
274  }
275 }
276 
277 
278 inline bool Foam::moleculeCloud::evaluatePotentialLimit
279 (
280  molecule* molI,
281  molecule* molJ
282 ) const
283 {
284  const pairPotentialList& pairPot = pot_.pairPotentials();
285 
286  const pairPotential& electrostatic = pairPot.electrostatic();
287 
288  label idI = molI->id();
289 
290  label idJ = molJ->id();
291 
292  const molecule::constantProperties& constPropI(constProps(idI));
293 
294  const molecule::constantProperties& constPropJ(constProps(idJ));
295 
296  List<label> siteIdsI = constPropI.siteIds();
297 
298  List<label> siteIdsJ = constPropJ.siteIds();
299 
300  List<bool> pairPotentialSitesI = constPropI.pairPotentialSites();
301 
302  List<bool> electrostaticSitesI = constPropI.electrostaticSites();
303 
304  List<bool> pairPotentialSitesJ = constPropJ.pairPotentialSites();
305 
306  List<bool> electrostaticSitesJ = constPropJ.electrostaticSites();
307 
308  forAll(siteIdsI, sI)
309  {
310  label idsI(siteIdsI[sI]);
311 
312  forAll(siteIdsJ, sJ)
313  {
314  label idsJ(siteIdsJ[sJ]);
315 
316  if (pairPotentialSitesI[sI] && pairPotentialSitesJ[sJ])
317  {
318  vector rsIsJ =
319  molI->sitePositions()[sI] - molJ->sitePositions()[sJ];
320 
321  scalar rsIsJMagSq = magSqr(rsIsJ);
322 
323  if (pairPot.rCutSqr(idsI, idsJ, rsIsJMagSq))
324  {
325  scalar rsIsJMag = mag(rsIsJ);
326 
327  // Guard against pairPot.energy being evaluated
328  // if rIJMag < SMALL. A floating point exception will
329  // happen otherwise.
330 
331  if (rsIsJMag < SMALL)
332  {
333  WarningIn("moleculeCloud::removeHighEnergyOverlaps()")
334  << "Molecule site pair closer than "
335  << SMALL
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."
341  << endl;
342 
343  return true;
344  }
345 
346  // Guard against pairPot.energy being evaluated if rIJMag <
347  // rMin. A tabulation lookup error will occur otherwise.
348 
349  if (rsIsJMag < pairPot.rMin(idsI, idsJ))
350  {
351  return true;
352  }
353 
354  if
355  (
356  mag(pairPot.energy(idsI, idsJ, rsIsJMag))
357  > pot_.potentialEnergyLimit()
358  )
359  {
360  return true;
361  };
362  }
363  }
364 
365  if (electrostaticSitesI[sI] && electrostaticSitesJ[sJ])
366  {
367  vector rsIsJ =
368  molI->sitePositions()[sI] - molJ->sitePositions()[sJ];
369 
370  scalar rsIsJMagSq = magSqr(rsIsJ);
371 
372  if (pairPot.rCutMaxSqr(rsIsJMagSq))
373  {
374  scalar rsIsJMag = mag(rsIsJ);
375 
376  // Guard against pairPot.energy being evaluated
377  // if rIJMag < SMALL. A floating point exception will
378  // happen otherwise.
379 
380  if (rsIsJMag < SMALL)
381  {
382  WarningIn("moleculeCloud::removeHighEnergyOverlaps()")
383  << "Molecule site pair closer than "
384  << SMALL
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."
390  << endl;
391 
392  return true;
393  }
394 
395  if (rsIsJMag < electrostatic.rMin())
396  {
397  return true;
398  }
399 
400  scalar chargeI = constPropI.siteCharges()[sI];
401 
402  scalar chargeJ = constPropJ.siteCharges()[sJ];
403 
404  if
405  (
406  mag(chargeI*chargeJ*electrostatic.energy(rsIsJMag))
407  > pot_.potentialEnergyLimit()
408  )
409  {
410  return true;
411  };
412  }
413  }
414  }
415  }
416 
417  return false;
418 }
419 
420 
421 inline bool Foam::moleculeCloud::evaluatePotentialLimit
422 (
423  molecule* molReal,
424  referredMolecule* molRef
425 ) const
426 {
427  const pairPotentialList& pairPot = pot_.pairPotentials();
428 
429  const pairPotential& electrostatic = pairPot.electrostatic();
430 
431  label idReal = molReal->id();
432 
433  label idRef = molRef->id();
434 
435  const molecule::constantProperties& constPropReal(constProps(idReal));
436 
437  const molecule::constantProperties& constPropRef(constProps(idRef));
438 
439  List<label> siteIdsReal = constPropReal.siteIds();
440 
441  List<label> siteIdsRef = constPropRef.siteIds();
442 
443  List<bool> pairPotentialSitesReal = constPropReal.pairPotentialSites();
444 
445  List<bool> electrostaticSitesReal = constPropReal.electrostaticSites();
446 
447  List<bool> pairPotentialSitesRef = constPropRef.pairPotentialSites();
448 
449  List<bool> electrostaticSitesRef = constPropRef.electrostaticSites();
450 
451  forAll(siteIdsReal, sReal)
452  {
453  label idsReal(siteIdsReal[sReal]);
454 
455  forAll(siteIdsRef, sRef)
456  {
457  label idsRef(siteIdsRef[sRef]);
458 
459  if (pairPotentialSitesReal[sReal] && pairPotentialSitesRef[sRef])
460  {
461  vector rsRealsRef =
462  molReal->sitePositions()[sReal]
463  - molRef->sitePositions()[sRef];
464 
465  scalar rsRealsRefMagSq = magSqr(rsRealsRef);
466 
467  if (pairPot.rCutSqr(idsReal, idsRef, rsRealsRefMagSq))
468  {
469  scalar rsRealsRefMag = mag(rsRealsRef);
470 
471  // Guard against pairPot.energy being evaluated
472  // if rRealRefMag < SMALL. A floating point exception will
473  // happen otherwise.
474 
475  if (rsRealsRefMag < SMALL)
476  {
477  WarningIn("moleculeCloud::removeHighEnergyOverlaps()")
478  << "Molecule site pair closer than "
479  << SMALL
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."
485  << endl;
486 
487  return true;
488  }
489 
490  // Guard against pairPot.energy being evaluated if
491  // rRealRefMag < rMin. A tabulation lookup error will occur
492  // otherwise.
493 
494  if (rsRealsRefMag < pairPot.rMin(idsReal, idsRef))
495  {
496  return true;
497  }
498 
499  if
500  (
501  mag(pairPot.energy(idsReal, idsRef, rsRealsRefMag))
502  > pot_.potentialEnergyLimit()
503  )
504  {
505  return true;
506  };
507  }
508  }
509 
510  if (electrostaticSitesReal[sReal] && electrostaticSitesRef[sRef])
511  {
512  vector rsRealsRef =
513  molReal->sitePositions()[sReal]
514  - molRef->sitePositions()[sRef];
515 
516  scalar rsRealsRefMagSq = magSqr(rsRealsRef);
517 
518  if (pairPot.rCutMaxSqr(rsRealsRefMagSq))
519  {
520  scalar rsRealsRefMag = mag(rsRealsRef);
521 
522  // Guard against pairPot.energy being evaluated
523  // if rRealRefMag < SMALL. A floating point exception will
524  // happen otherwise.
525 
526  if (rsRealsRefMag < SMALL)
527  {
528  WarningIn("moleculeCloud::removeHighEnergyOverlaps()")
529  << "Molecule site pair closer than "
530  << SMALL
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."
536  << endl;
537 
538  return true;
539  }
540 
541  if (rsRealsRefMag < electrostatic.rMin())
542  {
543  return true;
544  }
545 
546  scalar chargeReal = constPropReal.siteCharges()[sReal];
547 
548  scalar chargeRef = constPropRef.siteCharges()[sRef];
549 
550  if
551  (
552  mag
553  (
554  chargeReal
555  *chargeRef
556  *electrostatic.energy(rsRealsRefMag)
557  )
558  > pot_.potentialEnergyLimit()
559  )
560  {
561  return true;
562  };
563  }
564  }
565  }
566  }
567 
568  return false;
569 }
570 
571 
572 inline Foam::vector Foam::moleculeCloud::equipartitionLinearVelocity
573 (
574  scalar temperature,
575  scalar mass
576 )
577 {
578  return sqrt(kb*temperature/mass)*vector
579  (
580  rndGen_.GaussNormal(),
581  rndGen_.GaussNormal(),
582  rndGen_.GaussNormal()
583  );
584 }
585 
586 
587 inline Foam::vector Foam::moleculeCloud::equipartitionAngularMomentum
588 (
589  scalar temperature,
590  const molecule::constantProperties& cP
591 )
592 {
593  scalar sqrtKbT = sqrt(kb*temperature);
594 
595  if (cP.linearMolecule())
596  {
597  return sqrtKbT*vector
598  (
599  0.0,
600  sqrt(cP.momentOfInertia().yy())*rndGen_.GaussNormal(),
601  sqrt(cP.momentOfInertia().zz())*rndGen_.GaussNormal()
602  );
603  }
604  else
605  {
606  return sqrtKbT*vector
607  (
608  sqrt(cP.momentOfInertia().xx())*rndGen_.GaussNormal(),
609  sqrt(cP.momentOfInertia().yy())*rndGen_.GaussNormal(),
610  sqrt(cP.momentOfInertia().zz())*rndGen_.GaussNormal()
611  );
612  }
613 }
614 
615 
616 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
617 
619 {
620  return mesh_;
621 }
622 
623 
625 {
626  return pot_;
627 }
628 
629 
632 {
633  return cellOccupancy_;
634 }
635 
636 
637 inline const Foam::interactionLists&
639 {
640  return il_;
641 }
642 
643 
646 {
647  return constPropList_;
648 }
649 
650 
653 {
654  return constPropList_[id];
655 }
656 
657 
659 {
660  return rndGen_;
661 }
662 
663 
664 // ************************ vim: set sw=4 sts=4 et: ************************ //