C++ OBBond::GetBeginAtom方法代码示例

  void ReportFormat::WriteTorsions(ostream &ofs,OBMol &mol)
  {
    vector<OBBond*>::iterator bi1,bi2,bi3;
    OBBond* bond;
    OBAtom *a,*b,*c,*d;
    char buffer[BUFF_SIZE];

    //loop through all bonds generating torsions
    for(bond = mol.BeginBond(bi1); bond; bond = mol.NextBond(bi1))
      {
        b = bond->GetBeginAtom();
        c = bond->GetEndAtom();

        for(a = b->BeginNbrAtom(bi2);a;a = b->NextNbrAtom(bi2))
          {
            if(a == c)
              continue;

            for(d = c->BeginNbrAtom(bi3);d;d = c->NextNbrAtom(bi3))
              {
                if(d == b)
                  continue;

                snprintf(buffer, BUFF_SIZE, "%4d %4d %4d %4d %10.3f",
                        a->GetIdx(), b->GetIdx(),c->GetIdx(),d->GetIdx(),
                        CalcTorsionAngle(a->GetVector(), b->GetVector(),
                                         c->GetVector(), d->GetVector()));
                ofs << buffer << "\n";
              }
          }
      }
  }

  //! Sets atom->IsChiral() to true for chiral atoms
  void OBMol::FindChiralCenters()
  {
    if (HasChiralityPerceived())
      return;
    SetChiralityPerceived();

    obErrorLog.ThrowError(__FUNCTION__,
                          "Ran OpenBabel::FindChiralCenters", obAuditMsg);

    //do quick test to see if there are any possible chiral centers
    bool mayHaveChiralCenter=false;
    OBAtom *atom,*nbr;
    vector<OBAtom*>::iterator i;
    for (atom = BeginAtom(i);atom;atom = NextAtom(i))
      if (atom->GetHyb() == 3 && atom->GetHvyValence() >= 3)
        {
          mayHaveChiralCenter=true;
          break;
        }

    if (!mayHaveChiralCenter)
      return;

    OBBond *bond;
    vector<OBBond*>::iterator j;
    for (bond = BeginBond(j);bond;bond = NextBond(j))
      if (bond->IsWedge() || bond->IsHash())
        (bond->GetBeginAtom())->SetChiral();

    vector<unsigned int> vgid;
    GetGIDVector(vgid);
    vector<unsigned int> tlist;
    vector<unsigned int>::iterator k;

    bool ischiral;
    for (atom = BeginAtom(i);atom;atom = NextAtom(i))
      if (atom->GetHyb() == 3 && atom->GetHvyValence() >= 3 && !atom->IsChiral())
        {
          tlist.clear();
          ischiral = true;

          for (nbr = atom->BeginNbrAtom(j);nbr;nbr = atom->NextNbrAtom(j))
            {
              for (k = tlist.begin();k != tlist.end();++k)
                if (vgid[nbr->GetIdx()-1] == *k)
                  ischiral = false;

              if (ischiral)
                tlist.push_back(vgid[nbr->GetIdx()-1]);

              else
                break;
            }

          if (ischiral)
            atom->SetChiral();
        }
  }

 void OBRotorList::RemoveSymVals(OBMol &mol)
  {
    OBGraphSym gs(&mol);
    vector<unsigned int> sym_classes;
    gs.GetSymmetry(sym_classes);

    OBRotor *rotor;
    vector<OBRotor*>::iterator i;
    std::set<unsigned int> syms;
    for (rotor = BeginRotor(i);rotor;rotor = NextRotor(i)) {
      OBBond* bond = rotor->GetBond();
      OBAtom* end = bond->GetEndAtom();
      OBAtom* begin = bond->GetBeginAtom();
      int N_fold_symmetry = 1;
      for (int here=0; here <= 1; ++here) { // Try each side of the bond in turn
        
        OBAtom *this_side, *other_side;
        if (here == 0) {
          this_side = begin; other_side = end;
        }
        else {
          this_side = end; other_side = begin;
        }

        for (int hyb=2; hyb<=3; ++hyb) { // sp2 and sp3 carbons, with explicit Hs
          if (this_side->GetAtomicNum() == 6 && this_side->GetHyb() == hyb && this_side->GetValence() == (hyb + 1) ) {
            syms.clear();
            FOR_NBORS_OF_ATOM(nbr, this_side) {
              if ( &(*nbr) == other_side ) continue;
              syms.insert(sym_classes[nbr->GetIdx() - 1]);
            }
            if (syms.size() == 1) // All of the rotated atoms have the same symmetry class
              N_fold_symmetry *= hyb;
          }
        }
      }

      if (N_fold_symmetry  > 1) {
        size_t old_size = rotor->Size();
        rotor->RemoveSymTorsionValues(N_fold_symmetry);
        if (!_quiet) {
          cout << "...." << N_fold_symmetry << "-fold symmetry at rotor between " << 
                 begin->GetIdx() << " and " << end->GetIdx();
          cout << " - reduced from " << old_size << " to " << rotor->Size() << endl;
        }
      }
    }

  static void FindRings(OBMol &mol,vector<int> &path,OBBitVec &avisit,
                        OBBitVec &bvisit, int natom,int depth )
  {
    OBAtom *atom;
    OBBond *bond;
    vector<OBBond*>::iterator k;

    // don't return if all atoms are visited
    // (For example, some atoms are in multiple rings!) -GRH
      
    if (avisit[natom])
      {
        int j = depth-1;
        bond=mol.GetBond(path[j--]);
        bond->SetInRing();
        while( j >= 0 )
          {
            bond=mol.GetBond(path[j--]);
            bond->SetInRing();
            (bond->GetBeginAtom())->SetInRing();
            (bond->GetEndAtom())->SetInRing();
            if(bond->GetBeginAtomIdx()==static_cast<unsigned int>(natom) || bond->
               GetEndAtomIdx()==static_cast<unsigned int>(natom))
              break;
          }
      }
    else
      {
        avisit.SetBitOn(natom);
        atom = mol.GetAtom(natom);
        for(bond = atom->BeginBond(k);bond;bond=atom->NextBond(k))
          if( !bvisit[bond->GetIdx()])
            {
              path[depth] = bond->GetIdx();
              bvisit.SetBitOn(bond->GetIdx());
              FindRings(mol,path,avisit,bvisit,bond->GetNbrAtomIdx(atom),
                        depth+1);
            }
      }
  }

//.........这里部分代码省略.........
              {
                vector<OBResidue*>::iterator ri;
                for (res = mol.BeginResidue(ri) ; res ; res = mol.NextResidue(ri))
                  if (res->GetName() == resname &&
                      static_cast<int>(res->GetNum())
                      == resnum)
                    break;

                if (res == NULL)
                  {
                    res = mol.NewResidue();
                    res->SetName(resname);
                    res->SetNum(resnum);
                  }
              }
            OBAtom *atomPtr = mol.GetAtom(mol.NumAtoms());
            res->AddAtom(atomPtr);
            res->SetAtomID(atomPtr, atmid);
          } // end adding residue info
      }

    for (;;)
      {
        if (!ifs.getline(buffer,BUFF_SIZE))
          return(false);
        str = buffer;
        if (!strncmp(buffer,"@<TRIPOS>BOND",13))
          break;
      }

    int start,end,order;
    for (i = 0; i < nbonds; i++)
      {
        if (!ifs.getline(buffer,BUFF_SIZE))
          return(false);

        sscanf(buffer,"%*d %d %d %1024s",&start,&end,temp_type);
        str = temp_type;
        order = 1;
        if (str == "ar" || str == "AR" || str == "Ar")
          order = 5;
        else if (str == "AM" || str == "am" || str == "Am")
          order = 1;
        else
          order = atoi(str.c_str());

        mol.AddBond(start,end,order);
      }

 
    // update neighbour bonds information for each atom.
    vector<OBAtom*>::iterator apos;
    vector<OBBond*>::iterator bpos;
    OBAtom* patom;
    OBBond* pbond;
    
    for (patom = mol.BeginAtom(apos); patom; patom = mol.NextAtom(apos))
      {
        patom->ClearBond();
        for (pbond = mol.BeginBond(bpos); pbond; pbond = mol.NextBond(bpos))
          {
            if (patom == pbond->GetBeginAtom() || patom == pbond->GetEndAtom())
              {
                patom->AddBond(pbond);
              }
          }
      }

    // Suggestion by Liu Zhiguo 2008-01-26
    // Mol2 files define atom types -- there is no need to re-perceive
    mol.SetAtomTypesPerceived();
    mol.EndModify();

    //must add generic data after end modify - otherwise it will be blown away
    if (comment)
      {
        OBCommentData *cd = new OBCommentData;
        cd->SetData(comment);
        cd->SetOrigin(fileformatInput);
        mol.SetData(cd);
        delete [] comment;
        comment = NULL;
      }
    if (hasPartialCharges)
      mol.SetPartialChargesPerceived();

    // continue untill EOF or untill next molecule record
    streampos pos;
    for(;;)
      {
        pos = ifs.tellg();
        if (!ifs.getline(buffer,BUFF_SIZE))
          break;
        if (EQn(buffer,"@<TRIPOS>MOLECULE",17))
          break;
      }

    ifs.seekg(pos); // go back to the end of the molecule
    return(true);
  }

  bool OBChemTsfm::Apply(OBMol &mol)
  {
    if (!_bgn.Match(mol))
      return(false);
    mol.BeginModify();
    vector<vector<int> > mlist = _bgn.GetUMapList();

    obErrorLog.ThrowError(__FUNCTION__,
                          "Ran OpenBabel::OBChemTransform", obAuditMsg);

    if (!_vchrg.empty()) //modify charges
      {
        vector<vector<int> >::iterator i;
        vector<pair<int,int> >::iterator j;

        for (i = mlist.begin();i != mlist.end();++i)
          for (j = _vchrg.begin();j != _vchrg.end();++j)
            if (j->first < (signed)i->size()) { //goof proofing
              OBAtom *atom = mol.GetAtom((*i)[j->first]);
              int old_charge = atom->GetFormalCharge();
              atom->SetFormalCharge(j->second);
              int new_hcount = atom->GetImplicitHCount() + (j->second - old_charge);
              if (new_hcount < 0)
                new_hcount = 0;
              atom->SetImplicitHCount(new_hcount);
            }
      }

    if (!_vbond.empty()) //modify bond orders
      {
        OBBond *bond;
        vector<vector<int> >::iterator i;
        vector<pair<pair<int,int>,int> >::iterator j;
        for (i = mlist.begin();i != mlist.end();++i)
          for (j = _vbond.begin();j != _vbond.end();++j)
            {
              bond = mol.GetBond((*i)[j->first.first],(*i)[j->first.second]);
              if (!bond)
                {
                  obErrorLog.ThrowError(__FUNCTION__, "unable to find bond", obDebug);
                  continue;
                }
              unsigned int old_bond_order = bond->GetBondOrder();
              bond->SetBondOrder(j->second);
              for (int k = 0; k < 2; ++k) {
                OBAtom* atom = k == 0 ? bond->GetBeginAtom() : bond->GetEndAtom();
                int new_hcount = atom->GetImplicitHCount() - (j->second - old_bond_order);
                if (new_hcount < 0)
                  new_hcount = 0;
                atom->SetImplicitHCount(new_hcount);
              }
            }
      }

    if (!_vadel.empty() || !_vele.empty()) //delete atoms and change elements
      {
        vector<int>::iterator j;
        vector<vector<int> >::iterator i;

        if (!_vele.empty())
          {
            vector<pair<int,int> >::iterator k;
            for (i = mlist.begin();i != mlist.end();++i)
              for (k = _vele.begin();k != _vele.end();++k)
                mol.GetAtom((*i)[k->first])->SetAtomicNum(k->second);
          }

        //make sure same atom isn't deleted twice
        vector<bool> vda;
        vector<OBAtom*> vdel;
        vda.resize(mol.NumAtoms()+1,false);
        for (i = mlist.begin();i != mlist.end();++i)
          for (j = _vadel.begin();j != _vadel.end();++j)
            if (!vda[(*i)[*j]])
              {
                vda[(*i)[*j]] = true;
                vdel.push_back(mol.GetAtom((*i)[*j]));
              }

        vector<OBAtom*>::iterator k;
        for (k = vdel.begin();k != vdel.end();++k)
          mol.DeleteAtom((OBAtom*)*k);
      }

    mol.EndModify();
    return(true);
  }

bool ChemDrawXMLFormat::WriteMolecule(OBBase* pOb, OBConversion* pConv)
{
    static const xmlChar C_MOLECULE[]         = "fragment";
    static const xmlChar C_CDXML[]            = "CDXML";
    static const xmlChar C_BONDLENGTH[]       = "BondLength";
    static const xmlChar C_PAGE[]             = "page";
    static const xmlChar C_ATOM[]             = "n";
    static const xmlChar C_BOND[]             = "b";
    static const xmlChar C_ID[]               = "id";

    static const xmlChar C_CHARGE[]           = "Charge";
    static const xmlChar C_COORDS[]           = "p";
    static const xmlChar C_ELEMENT[]          = "Element";
    static const xmlChar C_ORDER[]            = "Order";
    static const xmlChar C_BEGIN[]            = "B";
    static const xmlChar C_END[]              = "E";
    static const xmlChar C_DISPLAY[]          = "Display";

    _pxmlConv = XMLConversion::GetDerived(pConv,false);
    if(!_pxmlConv)
        return false;

    OBMol* pmol = dynamic_cast<OBMol*>(pOb);
    if(pmol==NULL)
        return false;
    OBMol &mol = *pmol;

    OBBond *pbond;
    vector<OBBond*>::iterator j;
    if(_pxmlConv->GetOutputIndex() == 1)
    {
        xmlTextWriterStartDocument(writer(), NULL, NULL, NULL);
        xmlTextWriterWriteDTD(writer(), BAD_CAST "CDXML", NULL, BAD_CAST "http://www.camsoft.com/xml/cdxml.dtd", NULL);
        xmlTextWriterStartElement(writer(), C_CDXML);
        xmlTextWriterWriteFormatAttribute(writer(), C_BONDLENGTH , "30");
        xmlTextWriterStartElement(writer(), C_PAGE); // put everything on one page
        // now guess the average bond size for the first molecule and scale to 30.
        _scale = 0.;
        if (mol.NumBonds())
        {
            for (pbond = mol.BeginBond(j); pbond; pbond = mol.NextBond(j))
                _scale += pbond->GetLength();
            _scale /= mol.NumBonds();
        }
        else
            _scale = 1.; // FIXME: what happens if the molecule has no bond?
        _scale = 30. / _scale;
        _offset = 0;
    }

    xmlTextWriterStartElement(writer(), C_MOLECULE);

    OBAtom *patom;
    vector<OBAtom*>::iterator i;
    int n;
    for (patom = mol.BeginAtom(i); patom; patom = mol.NextAtom(i))
    {
        xmlTextWriterStartElement(writer(), C_ATOM);

        xmlTextWriterWriteFormatAttribute(writer(), C_ID , "%d", patom->GetIdx() + _offset);
        xmlTextWriterWriteFormatAttribute(writer(), C_COORDS , "%f %f", patom->GetX() * _scale, patom->GetY() * _scale);
        n = patom->GetAtomicNum();
        if (n != 6)
        {
            xmlTextWriterWriteFormatAttribute(writer(), C_ELEMENT , "%d", n);
        }
        n = patom->GetFormalCharge();
        if (n != 0)
        {
            xmlTextWriterWriteFormatAttribute(writer(), C_CHARGE , "%d", n);
        }
        xmlTextWriterEndElement(writer());
    }

    for (pbond = mol.BeginBond(j); pbond; pbond = mol.NextBond(j))
    {
        xmlTextWriterStartElement(writer(), C_BOND);
        patom = pbond->GetBeginAtom();
        xmlTextWriterWriteFormatAttribute(writer(), C_BEGIN , "%d", patom->GetIdx() + _offset);
        patom = pbond->GetEndAtom();
        xmlTextWriterWriteFormatAttribute(writer(), C_END , "%d", patom->GetIdx() + _offset);
        n = pbond->GetBO();
        if (n != 1)
        {
            xmlTextWriterWriteFormatAttribute(writer(), C_ORDER , "%d", n);
        }
        if (pbond->IsHash())
            xmlTextWriterWriteFormatAttribute(writer(), C_DISPLAY , "WedgeBegin");
        else if (pbond->IsWedge())
            xmlTextWriterWriteFormatAttribute(writer(), C_DISPLAY , "WedgedHashEnd");
        xmlTextWriterEndElement(writer());
    }
    _offset += mol.NumAtoms ();

    xmlTextWriterEndElement(writer());//molecule

    //TODO: Writing property block

    if(_pxmlConv->IsLast())
    {
//.........这里部分代码省略.........

  bool OBDepict::DrawMolecule(OBMol *mol)
  {
    if (!d->painter)
      return false;

    d->mol = mol;

    double width=0.0, height=0.0;

    OBAtom *atom;
    OBBondIterator j;
    OBAtomIterator i;

    if(mol->NumAtoms()>0) {
      // scale bond lengths
      double bondLengthSum = 0.0;
      for (OBBond *bond = mol->BeginBond(j); bond; bond = mol->NextBond(j))
        bondLengthSum += bond->GetLength();
      const double averageBondLength = bondLengthSum / mol->NumBonds();
      const double f = mol->NumBonds() ? d->bondLength / averageBondLength : 1.0;
      for (atom = mol->BeginAtom(i); atom; atom = mol->NextAtom(i))
        atom->SetVector(atom->GetX() * f, atom->GetY() * f, 0.0);

      // find min/max values
      double min_x, max_x;
      double min_y, max_y;
      atom = mol->BeginAtom(i);
      min_x = max_x = atom->GetX();
      min_y = max_y = atom->GetY();
      for (atom = mol->NextAtom(i); atom; atom = mol->NextAtom(i)) {
        min_x = std::min(min_x, atom->GetX());
        max_x = std::max(max_x, atom->GetX());
        min_y = std::min(min_y, atom->GetY());
        max_y = std::max(max_y, atom->GetY());
      }

      const double margin = 40.0;
      // translate all atoms so the bottom-left atom is at margin,margin
      for (atom = mol->BeginAtom(i); atom; atom = mol->NextAtom(i))
        atom->SetVector(atom->GetX() - min_x + margin, atom->GetY() - min_y + margin, 0.0);

      width  = max_x - min_x + 2*margin;
      height = max_y - min_y + 2*margin;
      
      //d->painter->SetPenWidth(d->penWidth);
      //d->painter->SetPenColor(d->pen));
      //d->painter->SetFillColor(OBColor("black"));
    }

    d->painter->NewCanvas(width, height);
    
    // draw bonds
    if(d->options & genWedgeHash)
      d->SetWedgeAndHash(mol);
    for (OBBond *bond = mol->BeginBond(j); bond; bond = mol->NextBond(j)) {
      OBAtom *begin = bond->GetBeginAtom();
      OBAtom *end = bond->GetEndAtom();

      if((d->options & internalColor) && bond->HasData("color"))
        d->painter->SetPenColor(OBColor(bond->GetData("color")->GetValue()));
      else
        d->painter->SetPenColor(d->bondColor);

      if (bond->IsWedge()) {
        d->DrawWedge(begin, end);
      } else if (bond->IsHash()) {
        d->DrawHash(begin, end);
      } else if (!bond->IsInRing()) {
        d->DrawSimpleBond(begin, end, bond->GetBO());
      }
    }
    
    // draw ring bonds
    std::vector<OBRing*> rings(mol->GetSSSR());
    OBBitVec drawnBonds;
    for (std::vector<OBRing*>::iterator k = rings.begin(); k != rings.end(); ++k) {
      OBRing *ring = *k;
      std::vector<int> indexes = ring->_path;
      vector3 center(VZero);
      for (std::vector<int>::iterator l = indexes.begin(); l != indexes.end(); ++l) {
        center += mol->GetAtom(*l)->GetVector();        
      }
      center /= indexes.size();

      for (unsigned int l = 0; l < indexes.size(); ++l) {
        OBAtom *begin = mol->GetAtom(indexes[l]);
        OBAtom *end;
        if (l+1 < indexes.size())
          end = mol->GetAtom(indexes[l+1]);
        else
          end = mol->GetAtom(indexes[0]);

        OBBond *ringBond = mol->GetBond(begin, end);
        if (drawnBonds.BitIsSet(ringBond->GetId()))
          continue;

        if((d->options & internalColor) && ringBond->HasData("color"))
          d->painter->SetPenColor(OBColor(ringBond->GetData("color")->GetValue()));
        else
          d->painter->SetPenColor(d->bondColor);
//.........这里部分代码省略.........

  void OutputCSTBonds(ostream &ofs, OBMol &mol, string prefix)
  {
    string bond_type;

    /* ---- Write povray-description of all bonds---- */
    for(unsigned int i = 0; i < mol.NumBonds(); ++i)
      {

        double x1,y1,z1,x2,y2,z2; /* Start and stop coordinates of a bond       */
        double dist;              /* Distance between (x1|y1|z1) and (x2|y2|z2) */
        double phi,theta;         /* Angles between (x1|y1|z1) and (x2|y2|z2)   */
        double dy;                /* Distance between (x1|0|z1) and (x2|0|z2)   */

        /* ---- Get a pointer to ith atom ---- */
        OBBond *bond = mol.GetBond(i);

        /* ---- Assign start of bond i ---- */
        x1 = (bond -> GetBeginAtom()) -> GetX();
        y1 = (bond -> GetBeginAtom()) -> GetY();
        z1 = (bond -> GetBeginAtom()) -> GetZ();

        /* ---- Assign end of bond i ---- */
        x2 = (bond -> GetEndAtom()) -> GetX();
        y2 = (bond -> GetEndAtom()) -> GetY();
        z2 = (bond -> GetEndAtom()) -> GetZ();

        /* ---- Calculate length of bond and (x1|0|z1) - (x2|0|z2) ---- */
        dist = sqrt(SQUARE(x2-x1) + SQUARE(y2-y1) + SQUARE(z2-z1));
        dy = sqrt(SQUARE(x2-x1) + SQUARE(z2-z1));

        /* ---- Calculate Phi and Theta ---- */
        phi = (double) 0.0;
        theta = (double) 0.0;
        if (fabs(dist) >= EPSILON)
          phi = acos((y2-y1)/dist);
        if (fabs(dy) >= EPSILON)
          theta = acos((x2-x1)/dy);

        /* ---- Begin of description of bond i (for a capped sticks model) ---- */
        ofs << "#declare " << prefix << "_bond" << i << " = object {" << endl;
        ofs << "\t  union {" << endl;

        /* ---- Begin of Start-Half of Bond (i) ---- */
        ofs << "\t   object {" << endl << "\t    bond_" << bond -> GetBondOrder()  << "\n";

        /* ---- Add a pigment - statement for start-atom of bond ---- */
        bond_type = bond->GetBeginAtom() -> GetType();
        bond_type.erase(remove_if(bond_type.begin(), bond_type.end(), bind1st(equal_to<char>(), '.')), bond_type.end());
        ofs << "\t    pigment{color Color_"
            << bond_type
            << "}" << endl;

        /* ---- Scale bond if needed ---- */
        if (fabs((double) 2.0 * dist) >= EPSILON)
          {

            /* ---- Print povray scale-statement (x-Axis) ---- */
            ofs << "\t    scale <" << (double) 0.5 * dist << ",1.0000,1.0000>" << endl;

          }

        /* ---- Rotate (Phi) bond if needed ---- */
        if (fabs(RAD2DEG(-phi) + (double) 90.0) >= EPSILON)
          {

            /* ---- Rotate along z-axis ---- */
            ofs << "\t    rotate <0.0000,0.0000,"
                << RAD2DEG(-phi) + (double) 90.0
                << ">" << endl;

          }

        /* ---- Check angle between (x1|0|z1) and (x2|0|z2) ---- */
        if (theta >= EPSILON)
          {

            /* ---- Check direction ---- */
            if ((z2 - z1) >= (double) 0.0)
              {

                /* ---- Rotate along y-Axis (negative) ---- */
                ofs << "\t    rotate <0.0000,"
                    << RAD2DEG((double) -1.0 *theta) << ",0.0000>"
                    << endl;

              }
            else
              {

                /* ---- Rotate along y-Axis (positive) ---- */
                ofs << "\t    rotate <0.0000," << RAD2DEG(theta) << ",0.0000>" << endl;

              }

          }

        /* ---- Translate bond to start ---- */

        ofs << "\t    translate " << prefix << "_pos_" << bond -> GetBeginAtomIdx() << endl;

//.........这里部分代码省略.........

  //! Calculate the signed volume for an atom.  If the atom has a valence of 3
  //! the coordinates of an attached hydrogen are calculated
  //! Puts attached Hydrogen last at the moment, like mol V3000 format.
  //! If ReZero=false (the default is true) always make pseudo z coords and leave them in mol
  double CalcSignedVolume(OBMol &mol,OBAtom *atm, bool ReZeroZ)
  {
    vector3 tmp_crd;
    vector<unsigned int> nbr_atms;
    vector<vector3> nbr_crds;
    bool use_central_atom = false,is2D=false;
    //   double hbrad = etab.CorrectedBondRad(1,0);
           
    if (!ReZeroZ || !mol.Has3D()) //give pseudo Z coords if mol is 2D
      {
        vector3 v,vz(0.0,0.0,1.0);
        is2D = true;
        OBAtom *nbr;
        OBBond *bond;
        vector<OBBond*>::iterator i;
        for (bond = atm->BeginBond(i);bond;bond = atm->NextBond(i))
          {
            nbr = bond->GetEndAtom();
            if (nbr != atm)
              {
                v = nbr->GetVector();
                if (bond->IsWedge())
                  v += vz;
                else
                  if (bond->IsHash())
                    v -= vz;

                nbr->SetVector(v);
              }
            else
              {
                nbr = bond->GetBeginAtom();
                v = nbr->GetVector();
                if (bond->IsWedge())
                  v -= vz;
                else
                  if (bond->IsHash())
                    v += vz;

                nbr->SetVector(v);
              }
          }
      }
    
    if (atm->GetHvyValence() < 3)
      {
        stringstream errorMsg;
        errorMsg << "Cannot calculate a signed volume for an atom with a heavy atom valence of " << atm->GetHvyValence() << endl;
        obErrorLog.ThrowError(__FUNCTION__, errorMsg.str(), obInfo);
        return(0.0);
      }

    // Create a vector with the coordinates of the neighbor atoms
    // Also make a vector with Atom IDs
    OBAtom *nbr;
    vector<OBBond*>::iterator bint;
    for (nbr = atm->BeginNbrAtom(bint);nbr;nbr = atm->NextNbrAtom(bint))
      {
        nbr_atms.push_back(nbr->GetIdx());
      }
    // sort the neighbor atoms to insure a consistent ordering
    sort(nbr_atms.begin(),nbr_atms.end());
    for (unsigned int i = 0; i < nbr_atms.size(); ++i)
      {
        OBAtom *tmp_atm = mol.GetAtom(nbr_atms[i]);
        nbr_crds.push_back(tmp_atm->GetVector());
      }
    /*
    // If we have three heavy atoms we need to calculate the position of the fourth
    if (atm->GetHvyValence() == 3)
    {
    double bondlen = hbrad+etab.CorrectedBondRad(atm->GetAtomicNum(),atm->GetHyb());
    atm->GetNewBondVector(tmp_crd,bondlen);
    nbr_crds.push_back(tmp_crd);
    }
    */
    for(unsigned int j=0;j < nbr_crds.size();++j) // Checks for a neighbour having 0 co-ords (added hydrogen etc)
      {
        // are the coordinates zero to 6 or more significant figures
        if (nbr_crds[j].IsApprox(VZero, 1.0e-6) && use_central_atom==false)
          use_central_atom=true;
        else if (nbr_crds[j].IsApprox(VZero, 1.0e-6))
          {
            obErrorLog.ThrowError(__FUNCTION__, "More than 2 neighbours have 0 co-ords when attempting 3D chiral calculation", obInfo);
          }
      }

    // If we have three heavy atoms we can use the chiral center atom itself for the fourth
    // will always give same sign (for tetrahedron), magnitude will be smaller.
    if(nbr_atms.size()==3 || use_central_atom==true)
      {
        nbr_crds.push_back(atm->GetVector());
        nbr_atms.push_back(mol.NumAtoms()+1); // meed to add largest number on end to work
      }
    OBChiralData* cd=(OBChiralData*)atm->GetData(OBGenericDataType::ChiralData); //Set the output atom4refs to the ones used
    if(cd==NULL)
//.........这里部分代码省略.........