Commit d3a863e7 authored by Axel Kohlmeyer's avatar Axel Kohlmeyer
Browse files

when identifying molecules/clusters fall back to unfiltered coordinates for ghost atoms

parent 64e80007
Loading
Loading
Loading
Loading
+69 −55
Original line number Diff line number Diff line
@@ -520,6 +520,8 @@ void FixReaxCSpecies::FindMolecule ()
  int *ilist;
  double bo_tmp,bo_cut;
  double **spec_atom = f_SPECBOND->array_atom;
  const double * const * const x = atom->x;
  const int nlocal = atom->nlocal;

  inum = reaxc->list->inum;
  ilist = reaxc->list->ilist;
@@ -553,11 +555,14 @@ void FixReaxCSpecies::FindMolecule ()
        for (jj = 0; jj < MAXSPECBOND; jj++) {
          j = reaxc->tmpid[i][jj];

      	  if (j < i) continue;
          if ((j == 0) || (j < i)) continue;
          if (!(mask[j] & groupbit)) continue;

      	  if (clusterID[i] == clusterID[j] && PBCconnected[i] == PBCconnected[j]
	    && x0[i].x == x0[j].x && x0[i].y == x0[j].y && x0[i].z == x0[j].z) continue;
          if (clusterID[i] == clusterID[j]
              && PBCconnected[i] == PBCconnected[j]
              && x0[i].x == x0[j].x
              && x0[i].y == x0[j].y
              && x0[i].z == x0[j].z) continue;

          jtype = atom->type[j];
          bo_cut = BOCut[itype][jtype];
@@ -567,10 +572,19 @@ void FixReaxCSpecies::FindMolecule ()
            clusterID[i] = clusterID[j] = MIN(clusterID[i], clusterID[j]);
            PBCconnected[i] = PBCconnected[j] = MAX(PBCconnected[i], PBCconnected[j]);
            x0[i] = x0[j] = chAnchor(x0[i], x0[j]);
            // spec_atom[][] contains filtered coordinates only for local atoms,
            // so we have to use unfiltered ones for ghost atoms.
            if (j < nlocal) {
              if ((fabs(spec_atom[i][1] - spec_atom[j][1]) > reaxc->control->bond_cut)
                  || (fabs(spec_atom[i][2] - spec_atom[j][2]) > reaxc->control->bond_cut)
                  || (fabs(spec_atom[i][3] - spec_atom[j][3]) > reaxc->control->bond_cut))
                PBCconnected[i] = PBCconnected[j] = 1;
            } else {
              if ((fabs(spec_atom[i][1] - x[j][1]) > reaxc->control->bond_cut)
                  || (fabs(spec_atom[i][2] - x[j][2]) > reaxc->control->bond_cut)
                  || (fabs(spec_atom[i][3] - x[j][3]) > reaxc->control->bond_cut))
                PBCconnected[i] = PBCconnected[j] = 1;
            }
            done = 0;
          }
        }