CFD Online Discussion Forums

CFD Online Discussion Forums (
-   OpenFOAM Bugs (
-   -   bugs in moleculeDynamics model? (

su_junwei May 8, 2009 19:33

bugs in moleculeDynamics model?
Dear foamers
there seems two bugs in the following two files

(1)file : moleculeCloudBuildDirectInteractionList.H: line 13
location: src/lagrangian/molecularDynamics/molecule/moleculeCloud

if face f is a boundary face, it hasn't a neighbour face and the operation of mesh_.faceNeighbour()[f] will exceed the array bound.
A possible solution is moving this line of code to line 40 in that file.

(2) file: pairPotential.C : line 148 and line 111
location: /src/lagrangian/molecularDynamics/potential/pairPotential/basic/

In function of forceLookup, if r=rCut, k=N-1;
and energyLookup[k+1] or will exceed the array bound.

It seems that the size of lookupTable definition is wrong.

In line pairPotential.C:82
label N = label((rCut_ - rMin_)/dr_) + 1;
should be changed to
label N = label((rCut_ - rMin_)/dr_) + 2;


gmacpherson May 18, 2009 09:00

Hi Junwei,

Thanks for spotting these. Both fixes have been made in 1.5.x.

1) Correct, I've moved the the neighbour cell label declaration within

if (mesh.isInternalFace(f))

as you suggest.

2) This has been fixed in pairPotentialList.C by making the cut-off radius test

if (rIJMagSqr < rCutSqr(a, b))

instead of

if (rIJMagSqr <= rCutSqr(a, b))

so it is not possible for r=rCut when the potential is evaluated. This is now in line with what's in the MD books by Allen & Tildesley and Rapaport who use U(r) = 0 when r >= rCut.



All times are GMT -4. The time now is 16:51.