Method MakeBonds

IO/Chemistry/vtkMoleculeReaderBase.cxx:276–384

}

unsigned int vtkMoleculeReaderBase::MakeBonds(
  vtkPoints* points, vtkIdTypeArray* atomTypes, vtkCellArray* newBonds)
{
  double X[3], Y[3];

  vtkSmartPointer<vtkPolyData> dataset = vtkSmartPointer<vtkPolyData>::New();
  dataset->SetPoints(points);

  vtkSmartPointer<vtkIdList> neighborAtoms = vtkSmartPointer<vtkIdList>::New();

  // Add atoms to the molecule first because an atom must
  // must be declared before bonds involving it.
  if (this->Molecule)
  {
    for (vtkIdType i = 0; i < this->NumberOfAtoms; ++i)
    {
      points->GetPoint(i, X);
      this->Molecule->AppendAtom(atomTypes->GetValue(i), X[0], X[1], X[2]);
    }
  }

  double max, dist, radius;
  vtkIdType bond[2];
  unsigned int numberOfBonds = 0;

  vtkSmartPointer<vtkPointLocator> locator = vtkSmartPointer<vtkPointLocator>::New();
  locator->SetDataSet(dataset);

  for (vtkIdType atomId = this->NumberOfAtoms - 1; atomId > 0; --atomId)
  {
    bond[0] = atomId;
    points->GetPoint(atomId, X);

    vtkIdType atom1Type = atomTypes->GetValue(atomId);

    /* Find all the atoms in the neighborhood at the max acceptable
     * bond distance
     */
    radius =
      (this->PeriodicTable->GetCovalentRadius(atom1Type) + 2.0 + 0.56) * std::max(BScale, HBScale);
    locator->FindPointsWithinRadius(radius, X, neighborAtoms);

    for (vtkIdType k = neighborAtoms->GetNumberOfIds() - 1; k >= 0; --k)
    {
      vtkIdType neighborAtomId = neighborAtoms->GetId(k);
      vtkIdType atom2Type = atomTypes->GetValue(neighborAtomId);

      // Skip points with which a bond may have already been created
      if (neighborAtomId >= atomId)
      {
        continue;
      }

      /*
       * The outer loop index 'atomId' is AFTER the inner loop 'neighborAtomId': 'atomId'
       * leads 'neighborAtomId' in the list: since hydrogens traditionally follow
       * the heavy atom they're bonded to, this makes it easy to quit
       * bonding to hydrogens after one bond is made by breaking out of