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;

Junwei

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.

Regards,

Graham