src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudRemoveHighEnergyOverlapsRealCellsCalculationStep.H

   1 idI = molI->id();
   2
   3 idJ = molJ->id();
   4
   5 rIJ = molI->position() - molJ->position();
   6
   7 rIJMagSq = magSqr(rIJ);
   8
   9 if (pairPotentials_.rCutSqr(idI, idJ, rIJMagSq))
  10 {
  11     rIJMag = mag(rIJ);
  12
  13     bool remove = false;
  14
  15     // Guard against pairPotentials_.energy being evaluated
  16     // if rIJMag < SMALL. A floating point exception will
  17     // happen otherwise.
  18
  19     if (rIJMag < SMALL)
  20     {
  21         WarningIn("moleculeCloud::removeHighEnergyOverlaps()")
  22             << "Real molecule pair "
  23             << " idI = " << idI
  24             << " at position " << molI->position()
  25             << " idJ = " << idJ
  26             << " at position " << molJ->position()
  27             << " are closer than " << SMALL
  28             << ": mag separation = " << rIJMag
  29             << ". These may have been placed on top of each"
  30             << " other by a rounding error in molConfig in parallel"
  31             << " or a block filled with molecules twice."
  32             << " Removing one of the molecules."
  33             << endl;
  34
  35         remove = true;
  36     }
  37
  38     // Guard against pairPotentials_.energy being evaluated
  39     // if rIJMag < rMin. A tabulation lookup error will occur otherwise.
  40
  41     if (rIJMag < pairPotentials_.rMin(idI, idJ))
  42     {
  43         remove = true;
  44     }
  45
  46     if (!remove)
  47     {
  48         if
  49         (
  50             pairPotentials_.energy(idI, idJ, rIJMag) > potentialEnergyLimit_
  51         )
  52         {
  53
  54             remove = true;
  55         }
  56     }
  57
  58     if (remove)
  59     {
  60         if
  61         (
  62             idI == idJ
  63          || findIndex(removalOrder_, idJ) < findIndex(removalOrder_, idI)
  64         )
  65         {
  66             if (findIndex(molsToDelete, molJ) == -1)
  67             {
  68                 molsToDelete.append(molJ);
  69             }
  70         }
  71         else
  72         {
  73             if (findIndex(molsToDelete, molI) == -1)
  74             {
  75                 molsToDelete.append(molI);
  76             }
  77         }
  78     }
  79 }