C++ OBMol::NextAtom方法代码示例

void ring_test()
{
  ostringstream os;

#ifdef TESTDATADIR
  string testdatadir = TESTDATADIR;
  string results_file = testdatadir + "ringresults.txt";
  string smilestypes_file = testdatadir + "attype.00.smi";
#else
  string results_file = "files/ringresults.txt";
  string smilestypes_file = "files/attype.00.smi";
#endif

  cout << "# Testing ring perception..." << endl;

  std::ifstream mifs;
  os << "Bail out! Cannot read test data " << smilestypes_file.c_str();
  BOOST_REQUIRE_MESSAGE( SafeOpen(mifs, smilestypes_file.c_str()), os.str().c_str() );

  std::ifstream rifs;
  os.str("");
  os << "Bail out! Cannot read test data " << results_file.c_str();
  BOOST_REQUIRE_MESSAGE( SafeOpen(rifs, results_file.c_str()), os.str().c_str() );

  unsigned int size;
  OBBond *bond;
  OBAtom *atom;
  int count;
  char buffer[BUFF_SIZE];
  vector<string> vs;
  vector<OBRing*> vr;
  vector<bool> vb;
  vector<int> vi;
  OBMol mol;
  vector<string>::iterator i;
  vector<OBBond*>::iterator j;
  vector<OBAtom*>::iterator k;
  vector<OBRing*>::iterator m;
  OBConversion conv(&mifs, &cout);
  unsigned int currentTest = 0;

  BOOST_REQUIRE_MESSAGE( conv.SetInAndOutFormats("SMI","SMI"), "Bail out! SMILES format is not loaded" );

  for (;mifs;)
    {
      mol.Clear();
      conv.Read(&mol);
      if (mol.Empty())
        continue;
      
      BOOST_REQUIRE_MESSAGE( rifs.getline(buffer,BUFF_SIZE), "Bail out! error reading reference data" );

      vb.clear();
      vb.resize(mol.NumBonds(),false);
      //check ring bonds
      tokenize(vs,buffer);
      for (i = vs.begin();i != vs.end();i++)
        vb[atoi(i->c_str())] = true;

      for (bond = mol.BeginBond(j);bond;bond = mol.NextBond(j))
        {
          os.str("");
	  os << "ring bond data different than reference # Molecule: " << mol.GetTitle();
          BOOST_CHECK_MESSAGE( vb[bond->GetIdx()] == bond->IsInRing(), os.str().c_str() );
        }

      vr = mol.GetSSSR();
      BOOST_REQUIRE_MESSAGE( rifs.getline(buffer,BUFF_SIZE), "Bail out! error reading reference data" );
      sscanf(buffer,"%d",&size);

      os.str("");
      os << "SSSR size different than reference # Molecule: " << mol.GetTitle();
      BOOST_CHECK_MESSAGE( vr.size() == size, os.str().c_str() ); //check SSSR size

      BOOST_REQUIRE_MESSAGE( rifs.getline(buffer,BUFF_SIZE), "Bail out! error reading reference data" );

      tokenize(vs,buffer);
      i = vs.begin();
      for (atom = mol.BeginAtom(k);atom;atom = mol.NextAtom(k))
        {
          os.str("");
          os << "error in SSSR count # Molecule: " << mol.GetTitle();
          BOOST_CHECK_MESSAGE( i != vs.end(), os.str().c_str() ); //check SSSR size

          count = 0;
          for (m = vr.begin();m != vr.end();m++)
            if ((*m)->_pathset[atom->GetIdx()])
              count++;

	  os.str("");
          os << "ring membership test failed # Molecule: " << mol.GetTitle();
          BOOST_CHECK_MESSAGE( atoi(i->c_str()) == count, os.str().c_str() ); 

          i++;
        }
    }

}

  bool OBResidueData::AssignBonds(OBMol &mol,OBBitVec &bv)
  {
    if (!_init)
      Init();

    OBAtom *a1,*a2;
    OBResidue *r1,*r2;
    vector<OBAtom*>::iterator i,j;
    vector3 v;

    int bo;
    string skipres = ""; // Residue Number to skip
    string rname = "";
    //assign residue bonds
    for (a1 = mol.BeginAtom(i);a1;a1 = mol.NextAtom(i))
      {
        r1 = a1->GetResidue();
        if (r1 == NULL) // atoms may not have residues
          continue;

        if (skipres.length() && r1->GetNumString() == skipres)
          continue;

        if (r1->GetName() != rname)
          {
            skipres = SetResName(r1->GetName()) ? "" : r1->GetNumString();
            rname = r1->GetName();
          }
        //assign bonds for each atom
        for (j=i,a2 = mol.NextAtom(j);a2;a2 = mol.NextAtom(j))
          {
            r2 = a2->GetResidue();
            if (r2 == NULL) // atoms may not have residues
              continue;

            if (r1->GetNumString() != r2->GetNumString())
              break;
            if (r1->GetName() != r2->GetName())
              break;
            if (r1->GetChain() != r2->GetChain())
              break; // Fixes PR#2889763 - Fabian

            if ((bo = LookupBO(r1->GetAtomID(a1),r2->GetAtomID(a2))))
              {
                // Suggested by Liu Zhiguo 2007-08-13
                // for predefined residues, don't perceive connection
                // by distance
                //                v = a1->GetVector() - a2->GetVector();
                //                if (v.length_2() < 3.5) //check by distance
                  mol.AddBond(a1->GetIdx(),a2->GetIdx(),bo);
              }
          }
      }

    int hyb;
    string type;

    //types and hybridization
    rname = ""; // name of current residue
    skipres = ""; // don't skip any residues right now
    for (a1 = mol.BeginAtom(i);a1;a1 = mol.NextAtom(i))
      {
        if (a1->GetAtomicNum() == OBElements::Oxygen && !a1->GetValence())
          {
            a1->SetType("O3");
            continue;
          }
        if (a1->GetAtomicNum() == OBElements::Hydrogen)
          {
            a1->SetType("H");
            continue;
          }

        //***valence rule for O-
        if (a1->GetAtomicNum() == OBElements::Oxygen && a1->GetValence() == 1)
          {
            OBBond *bond;
            bond = (OBBond*)*(a1->BeginBonds());
            if (bond->GetBO() == 2)
              {
                a1->SetType("O2");
                a1->SetHyb(2);
              }
            else if (bond->GetBO() == 1)
              {
                // Leave the protonation/deprotonation to phmodel.txt
                a1->SetType("O3");
                a1->SetHyb(3);
                // PR#3203039 -- Fix from Magnus Lundborg
                //                a1->SetFormalCharge(0);
              }
            continue;
          }

        r1 = a1->GetResidue();
        if (r1 == NULL) continue; // atoms may not have residues
        if (skipres.length() && r1->GetNumString() == skipres)
          continue;

        if (r1->GetName() != rname)
//.........这里部分代码省略.........

  bool parseConectRecord(char *buffer,OBMol &mol)
  {
    stringstream errorMsg;
    string clearError;

    // Setup strings and string buffers
    vector<string> vs;
    buffer[70] = '\0';
    if (strlen(buffer) < 70)
      {
        errorMsg << "WARNING: Problems reading a PDB file\n"
                 << "  Problems reading a CONECT record.\n"
                 << "  According to the PDB specification,\n"
                 << "  the record should have 70 columns, but OpenBabel found "
                 << strlen(buffer) << " columns." << endl;
        obErrorLog.ThrowError(__FUNCTION__, errorMsg.str() , obInfo);
        errorMsg.str(clearError);
      }

    // Serial number of the first atom, read from column 7-11 of the
    // connect record, to which the other atoms connect to.
    long int startAtomSerialNumber;
    // A pointer to the first atom.
    OBAtom *firstAtom = NULL;
    // Serial numbers of the atoms which bind to firstAtom, read from
    // columns 12-16, 17-21, 22-27 and 27-31 of the connect record. Note
    // that we reserve space for 5 integers, but read only four of
    // them. This is to simplify the determination of the bond order;
    // see below.
    long int boundedAtomsSerialNumbers[5]  = {0,0,0,0,0};
    // Bools which tell us which of the serial numbers in
    // boundedAtomsSerialNumbers are read from the file, and which are
    // invalid
    bool boundedAtomsSerialNumbersValid[5] = {false, false, false, false, false};

    // Pragmatic approach -- too many non-standard PDB files out there
    // (including some old ones from us)
    // So if we have a small number of atoms, then try to break by spaces
    // Otherwise (i.e., NumAtoms() > 9,999 we need to go by position)
    // We'll switch back and forth a few times to save duplicating common code

    if (mol.NumAtoms() <= 9999)
      {
        // make sure we don't look at salt bridges or whatever, so cut the buffer short
        buffer[32] = '\0';
        tokenize(vs,buffer);
        if( vs.empty() || vs.size() < 2)
          return false;
        vs.erase(vs.begin()); // remove "CONECT"

        startAtomSerialNumber = atoi(vs[0].c_str());
      }
    else
      {
        if (readIntegerFromRecord(buffer, 7, &startAtomSerialNumber) == false)
          {
            errorMsg << "WARNING: Problems reading a PDB file\n"
                     << "  Problems reading a CONECT record.\n"
                     << "  According to the PDB specification,\n"
                     << "  columns 7-11 should contain the serial number of an atom.\n"
                     << "  THIS CONECT RECORD WILL BE IGNORED." << endl;
            obErrorLog.ThrowError(__FUNCTION__, errorMsg.str() , obWarning);
            return(false);
          }
      }

    vector<OBAtom*>::iterator i;
    for (OBAtom *a1 = mol.BeginAtom(i);a1;a1 = mol.NextAtom(i)) {
      // atoms may not have residue information, but if they do,
      // check serial numbers
      if (a1->GetResidue() != NULL &&
          static_cast<long int>(a1->GetResidue()->
                                GetSerialNum(a1)) == startAtomSerialNumber)
        {
          firstAtom = a1;
          break;
        }
    }

    if (firstAtom == NULL)
      {
        errorMsg << "WARNING: Problems reading a PDB file:\n"
                 << "  Problems reading a CONECT record.\n"
                 << "  According to the PDB specification,\n"
                 << "  columns 7-11 should contain the serial number of an atom.\n"
                 << "  No atom was found with this serial number.\n"
                 << "  THIS CONECT RECORD WILL BE IGNORED." << endl;
        obErrorLog.ThrowError(__FUNCTION__, errorMsg.str() , obWarning);
        return(false);
      }

    if (mol.NumAtoms() < 9999)
      {
        if (vs.size() > 1) boundedAtomsSerialNumbers[0] = atoi(vs[1].c_str());
        if (vs.size() > 2) boundedAtomsSerialNumbers[1] = atoi(vs[2].c_str());
        if (vs.size() > 3) boundedAtomsSerialNumbers[2] = atoi(vs[3].c_str());
        if (vs.size() > 4) boundedAtomsSerialNumbers[3] = atoi(vs[4].c_str());

        unsigned int limit = 4;
        if (vs.size() <= 4)
//.........这里部分代码省略.........

  void OBAromaticTyper::AssignAromaticFlags(OBMol &mol)
  {
    if (!_init)
      Init();

    if (mol.HasAromaticPerceived())
      return;
    mol.SetAromaticPerceived();
    obErrorLog.ThrowError(__FUNCTION__,
                          "Ran OpenBabel::AssignAromaticFlags", obAuditMsg);

    _vpa.clear();
    _vpa.resize(mol.NumAtoms()+1);
    _velec.clear();
    _velec.resize(mol.NumAtoms()+1);
    _root.clear();
    _root.resize(mol.NumAtoms()+1);

    OBBond *bond;
    OBAtom *atom;
    vector<OBAtom*>::iterator i;
    vector<OBBond*>::iterator j;

    //unset all aromatic flags
    for (atom = mol.BeginAtom(i);atom;atom = mol.NextAtom(i))
      atom->UnsetAromatic();
    for (bond = mol.BeginBond(j);bond;bond = mol.NextBond(j))
      bond->UnsetAromatic();

    int idx;
    vector<vector<int> >::iterator m;
    vector<OBSmartsPattern*>::iterator k;

    //mark atoms as potentially aromatic
    for (idx=0,k = _vsp.begin();k != _vsp.end();++k,++idx)
      if ((*k)->Match(mol))
        {
          _mlist = (*k)->GetMapList();
          for (m = _mlist.begin();m != _mlist.end();++m)
            {
              _vpa[(*m)[0]] = true;
              _velec[(*m)[0]] = _verange[idx];
            }
        }

    //sanity check - exclude all 4 substituted atoms and sp centers
    for (atom = mol.BeginAtom(i);atom;atom = mol.NextAtom(i))
      {
        if (atom->GetImplicitValence() > 3)
          {
            _vpa[atom->GetIdx()] = false;
            continue;
          }

        switch(atom->GetAtomicNum())
          {
            //phosphorus and sulfur may be initially typed as sp3
          case 6:
          case 7:
          case 8:
            if (atom->GetHyb() != 2)
              _vpa[atom->GetIdx()] = false;
            break;
          }
      }

    //propagate potentially aromatic atoms
    for (atom = mol.BeginAtom(i);atom;atom = mol.NextAtom(i))
      if (_vpa[atom->GetIdx()])
        PropagatePotentialAromatic(atom);

    //select root atoms
    SelectRootAtoms(mol);

    ExcludeSmallRing(mol); //remove 3 membered rings from consideration

    //loop over root atoms and look for aromatic rings
    _visit.clear();
    _visit.resize(mol.NumAtoms()+1);
    for (atom = mol.BeginAtom(i);atom;atom = mol.NextAtom(i))
      if (_root[atom->GetIdx()])
        CheckAromaticity(atom,14);

    //for (atom = mol.BeginAtom(i);atom;atom = mol.NextAtom(i))
    //	  if (atom->IsAromatic())
    //		  cerr << "aro = " <<atom->GetIdx()  << endl;

    //for (bond = mol.BeginBond(j);bond;bond = mol.NextBond(j))
    //if (bond->IsAromatic())
    //cerr << bond->GetIdx() << ' ' << bond->IsAromatic() << endl;
  }

  //! \todo Make this method bulletproof. Currently it causes segfaults sometimes
  void CacaoFormat::SetHilderbrandt(OBMol &mol,vector<OBInternalCoord*> &vit)
  {
    // Roundtrip testing shows that some atoms are NULL
    //  which causes segfaults when dereferencing later
    //   (e.g. in the last "segment" of this routine)
    double sum,r;

    OBAtom dummy1,dummy2;
    dummy1.SetVector(0.0,0.0,1.0);
    dummy2.SetVector(1.0,0.0,0.0);

    OBAtom *atom,*a1,*a2,*ref;
    vector<OBAtom*>::iterator ai;

    vit.push_back((OBInternalCoord*)NULL);
    for (atom = mol.BeginAtom(ai);atom;atom = mol.NextAtom(ai))
      vit.push_back(new OBInternalCoord (atom));

    vit[1]->_a = &dummy1;
    vit[1]->_b = &dummy2;
    if (vit.size() > 2) {
      vit[2]->_b = &dummy1;
      vit[2]->_c = &dummy2;
      if  (vit.size() > 3) {
        vit[3]->_c = &dummy1;
      }
    }

    unsigned int i,j;
    for (i = 2;i <= mol.NumAtoms();i++)
      {
        ref = (OBAtom*)NULL;
        a1 = mol.GetAtom(i);
        sum = 100.0;
        for (j = 1;j < i;j++)
          {
            a2 = mol.GetAtom(j);
            r = (a1->GetVector()-a2->GetVector()).length_2();
            if ((r < sum) && (vit[j]->_a != a2) && (vit[j]->_b != a2))
              {
                sum = r;
                ref = a2;
              }
          }
        vit[i]->_a = ref;
      }

    for (i = 3;i <= mol.NumAtoms();i++)
      vit[i]->_b = vit[vit[i]->_a->GetIdx()]->_a;

    for (i = 4;i <= mol.NumAtoms();i++)
      {
        if (vit[i]->_b && vit[i]->_b->GetIdx())
          vit[i]->_c = vit[vit[i]->_b->GetIdx()]->_b;
        else
          vit[i]->_c = &dummy1;
      }

    OBAtom *a,*b,*c;
    vector3 v1,v2;
    for (i = 2;i <= mol.NumAtoms();i++)
      {
        atom = mol.GetAtom(i);
        a = vit[i]->_a;
        b = vit[i]->_b;
        c = vit[i]->_c;
        v1 = atom->GetVector() - a->GetVector();
        v2 = b->GetVector() - a->GetVector();
        vit[i]->_ang = vectorAngle(v1,v2);
        vit[i]->_tor = CalcTorsionAngle(atom->GetVector(),
                                        a->GetVector(),
                                        b->GetVector(),
                                        c->GetVector());
        vit[i]->_dst = (vit[i]->_a->GetVector() - atom->GetVector()).length();
      }

  }

bool Database::readTreeSmi (string smi, Tid tid, Tid orig_tid, int line_nr) {

    OBMol mol;

    istringstream iss (smi, istringstream::in);
    OBConversion conv(&iss,&cout);

    // read the molecule
    conv.SetInAndOutFormats("SMI","SDF");

    if (!conv.ReadString(&mol,smi)) {
        cerr << "Error during conversion" << endl;
        return(0);
    }   

    if (!mol.DeleteHydrogens()) {
        cerr << "Unable to delete hydrogens" << endl;
        return(0);
    }

    // create and store new tree object
    DatabaseTreePtr tree = new DatabaseTree ( tid , orig_tid , line_nr );


    // SEGFAULT
    trees.push_back(tree);
    trees_map[orig_tid] = tree;

    int nodessize = 0, edgessize = 0;
    static vector<DatabaseTreeNode> nodes;
    static vector<vector<DatabaseTreeEdge> > edges;
    nodes.resize ( 0 );

//    cerr << "Atoms are (Type(ID)):" << endl;
    OBAtomIterator atom;
    mol.BeginAtom(atom);

	///////////
	// NODES //
	///////////

//    cerr << endl;
    do {

        //cerr << " " << (*atom)->GetType() << " (idx " << (*atom)->GetIdx() << ")" << endl;

        InputNodeLabel inputnodelabel=0;

        // set atom type as label
        // code for 'c' is set to -1 (aromatic carbon).
        if (fm::aromatic) {
            (*atom)->IsAromatic() ? inputnodelabel = (*atom)->GetAtomicNum()+150 : inputnodelabel = (*atom)->GetAtomicNum();
        }
        else inputnodelabel = (*atom)->GetAtomicNum();
        nodessize++;

        // Insert into map, using subsequent numbering for internal labels:
    	// node nr. 1, node nr. 2, ...
	    // Store direction inputlabel -> internal label
    	//cerr << "NodeLabelMap: insert " << inputnodelabel << " --> " << nodelabels.size() << endl;
        map_insert_pair ( nodelabelmap ) p = nodelabelmap.insert ( make_pair ( inputnodelabel, (nodelabels.size()) ) );

        // Store direction internal label -> node
        // if node has NOT been present, label it and set frequency to 1
        if ( p.second ) {
          vector_push_back ( DatabaseNodeLabel, nodelabels, nodelabel );
          nodelabel.inputlabel = inputnodelabel;
          nodelabel.occurrences.parent = NULL;
          nodelabel.occurrences.number = 1;
          nodelabel.lasttid = tid;
          //cerr << "Created node label " << nodelabel.inputlabel << " (freq " << nodelabel.frequency <<  ")" << endl;
        }
        // if node has been present, increase frequency
        else {
          DatabaseNodeLabel &nodelabel = nodelabels[p.first->second];
          if ( nodelabel.lasttid != tid ) nodelabel.frequency++;
          nodelabel.lasttid = tid;
          //cerr << "Updated node label " << nodelabel.inputlabel << " (freq " << nodelabel.frequency << ")" << endl;
        }

        // Tree node
        vector_push_back ( DatabaseTreeNode, nodes, node );
        node.nodelabel = p.first->second;                           // refer to nodelabel
    //    node.atom = (*atom);                                        // attach OB atom
        node.incycle = false;
        //cerr << "Created tree node for OB index " << node.atom->GetIdx() << " (nodelabel " << (int) node.nodelabel << ", nodes[] size " << nodessize << ")" << endl;

    } while (mol.NextAtom(atom));


    // copy nodes to tree and prepare edges storage size
    // edges[nodeid] gives the edges going out of node with id 'nodeid'
    tree->nodes.reserve ( nodessize );
    if ( edges.size () < (unsigned int) nodessize )
        edges.resize ( nodessize );					// edges stored per node
    for ( int i = 0; i < nodessize; i++ ) {
        edges[i].resize ( 0 );						// no edges yet
        tree->nodes.push_back ( nodes[i] );
    }
    //cerr << endl;
//.........这里部分代码省略.........

int main(int argc,char **argv)
{
  char *program_name= argv[0];
  char *FileIn = NULL;

  if (argc != 2)
    {
      cout << "Usage: " << program_name << " <filename>" << endl;
      exit(-1);
    }
  else
    {
      FileIn  = argv[1];
      //   const char* p = strrchr(FileIn,'.');
    }

  // Find Input filetype
  OBConversion conv;
  OBFormat *format = conv.FormatFromExt(FileIn);
    
  if (!format || !conv.SetInAndOutFormats(format, format))
    {
      cerr << program_name << ": cannot read input format!" << endl;
      exit (-1);
    }

  ifstream ifs;

  // Read the file
  ifs.open(FileIn);
  if (!ifs)
    {
      cerr << program_name << ": cannot read input file!" << endl;
      exit (-1);
    }

  OBMol mol;
  OBAtom *atom;

  for (int c=1;;++c) // big for loop (replace with do while?)
    {
      mol.Clear();
      conv.Read(&mol, &ifs);
      if (mol.Empty())
        break;
      cout << "Molecule "<< c << ": " << mol.GetTitle() << endl;
      //mol.FindChiralCenters(); // labels all chiral atoms
      vector<OBAtom*>::iterator i; // iterate over all atoms
      for (atom = mol.BeginAtom(i);atom;atom = mol.NextAtom(i))
        {
          if(!atom->IsChiral())continue; // aborts if atom isn't chiral
          cout << "Atom " << atom->GetIdx() << " Is Chiral ";
          cout << atom->GetType()<<endl;
        
          OBChiralData* cd = (OBChiralData*)atom->GetData(OBGenericDataType::ChiralData);
        
          if (cd){
            vector<unsigned int> x=cd->GetAtom4Refs(input);
            size_t n=0;
            cout <<"Atom4refs:";
            for (n=0;n<x.size();++n)
              cout <<" "<<x[n];
            cout <<endl;
          }
          else{cd=new OBChiralData;atom->SetData(cd);}
          vector<unsigned int> _output;
          unsigned int n;
          for(n=1;n<5;++n) _output.push_back(n);
          cd->SetAtom4Refs(_output,output);
          /* // MOLV3000 uses 1234 unless an H then 123H
             if (atom->GetHvyValence()==3)
             {
             OBAtom *nbr;
             int Hid=1000;// max Atom ID +1 should be used here
             vector<unsigned int> nbr_atms;
             vector<OBBond*>::iterator i;
             for (nbr = atom->BeginNbrAtom(i);nbr;nbr = atom->NextNbrAtom(i))
             {
             if (nbr->IsHydrogen()){Hid=nbr->GetIdx();continue;}
             nbr_atms.push_back(nbr->GetIdx());
             }
             sort(nbr_atms.begin(),nbr_atms.end());
             nbr_atms.push_back(Hid);
             OBChiralData* cd=(OBChiralData*)atom->GetData(OBGenericDataType::ChiralData);
             cd->SetAtom4Refs(nbr_atms,output);   
             } 
             else if (atom->GetHvyValence()==4)
             {
             OBChiralData* cd=(OBChiralData*)atom->GetData(OBGenericDataType::ChiralData);
             vector<unsigned int> nbr_atms;
             int n;
             for(n=1;n<5;++n)nbr_atms.push_back(n);
             cd->SetAtom4Refs(nbr_atms,output); 
             } */
    /* FIXME          
          if (!mol.HasNonZeroCoords())
            {
              cout << "Calcing 0D chirality "<< CorrectChirality(mol,atom)<<endl;
            }
          else {
//.........这里部分代码省略.........

//.........这里部分代码省略.........
              {
                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)
      {
        cd = new OBChiralData;
        cd->SetOrigin(perceived);
        atm->SetData(cd);
      }
    cd->SetAtom4Refs(nbr_atms,calcvolume);
    
    //re-zero psuedo-coords
    if (is2D && ReZeroZ)
      {
        vector3 v;
        OBAtom *atom;
        vector<OBAtom*>::iterator k;
        for (atom = mol.BeginAtom(k);atom;atom = mol.NextAtom(k))
          {
            v = atom->GetVector();
            v.SetZ(0.0);
            atom->SetVector(v);
          }
      }
    
    return(signed_volume(nbr_crds[0],nbr_crds[1],nbr_crds[2],nbr_crds[3]));
  }

  //! Use the supplied OBMol and its Cartesian coordinates to generate
  //! a set of internal (z-matrix) coordinates as supplied in the
  //! vector<OBInternalCoord*> argument.
  //! Implements <a href="http://qsar.sourceforge.net/dicts/blue-obelisk/index.xhtml#cartesianCoordinatesIntoZmatrixCoordinates">blue-obelisk:cartesianCoordinatesIntoZmatrixCoordinates</a>.
  void CartesianToInternal(std::vector<OBInternalCoord*> &vic,OBMol &mol)
  {
    double r,sum;
    OBAtom *atom,*nbr,*ref;
    vector<OBAtom*>::iterator i,j,m;

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

    //set reference atoms
    for (atom = mol.BeginAtom(i);atom;atom = mol.NextAtom(i))
      {
        if      (atom->GetIdx() == 1)
          continue;
        else if (atom->GetIdx() == 2)
          {
            vic[atom->GetIdx()]->_a = mol.GetAtom(1);
            continue;
          }
        else if (atom->GetIdx() == 3)
          {
            if( (atom->GetVector()-mol.GetAtom(2)->GetVector()).length_2()
                <(atom->GetVector()-mol.GetAtom(1)->GetVector()).length_2())
              {
                vic[atom->GetIdx()]->_a = mol.GetAtom(2);
                vic[atom->GetIdx()]->_b = mol.GetAtom(1);
              }
            else
              {
                vic[atom->GetIdx()]->_a = mol.GetAtom(1);
                vic[atom->GetIdx()]->_b = mol.GetAtom(2);
              }
            continue;
          }
        sum=1.0E10;
        ref = mol.GetAtom(1);
        for(nbr = mol.BeginAtom(j);nbr && (i != j);nbr = mol.NextAtom(j))
          {
            r = (atom->GetVector()-nbr->GetVector()).length_2();
            if((r < sum) && (vic[nbr->GetIdx()]->_a != nbr) &&
               (vic[nbr->GetIdx()]->_b != nbr))
              {
                sum = r;
                ref = nbr;
              }
          }

        vic[atom->GetIdx()]->_a = ref;
        if (ref->GetIdx() >= 3)
          {
            vic[atom->GetIdx()]->_b = vic[ref->GetIdx()]->_a;
            vic[atom->GetIdx()]->_c = vic[ref->GetIdx()]->_b;
          }
        else
          {
            if(ref->GetIdx()== 1)
              {
                vic[atom->GetIdx()]->_b = mol.GetAtom(2);
                vic[atom->GetIdx()]->_c = mol.GetAtom(3);
              }
            else
              {//ref->GetIdx()== 2
                vic[atom->GetIdx()]->_b = mol.GetAtom(1);
                vic[atom->GetIdx()]->_c = mol.GetAtom(3);
              }
          }
      }

    //fill in geometries
    unsigned int k;
    vector3 v1,v2;
    OBAtom *a,*b,*c;
    for (k = 2;k <= mol.NumAtoms();++k)
      {
        atom = mol.GetAtom(k);
        a = vic[k]->_a;
        b = vic[k]->_b;
        c = vic[k]->_c;
        v1 = atom->GetVector() - a->GetVector();
        vic[k]->_dst = v1.length();
        if (k == 2)
          continue;

        v2 = b->GetVector()    - a->GetVector();
        vic[k]->_ang = vectorAngle(v1,v2);
        if (k == 3)
          continue;

        vic[k]->_tor = CalcTorsionAngle(atom->GetVector(),
                                        a->GetVector(),
                                        b->GetVector(),
                                        c->GetVector());
      }

    //check for linear geometries and try to correct if possible
    bool done;
//.........这里部分代码省略.........

  //! Transform the supplied vector<OBInternalCoord*> into cartesian and update
  //! the OBMol accordingly. The size of supplied internal coordinate vector
  //! has to be the same as the number of atoms in molecule (+ NULL in the
  //! beginning).
  //! Implements <a href="http://qsar.sourceforge.net/dicts/blue-obelisk/index.xhtml#zmatrixCoordinatesIntoCartesianCoordinates">blue-obelisk:zmatrixCoordinatesIntoCartesianCoordinates</a>
  void InternalToCartesian(std::vector<OBInternalCoord*> &vic,OBMol &mol)
  {
    vector3 n,nn,v1,v2,v3,avec,bvec,cvec;
    double dst = 0.0, ang = 0.0, tor = 0.0;
    OBAtom *atom;
    vector<OBAtom*>::iterator i;
    unsigned int index;

    if (vic.empty())
      return;

    if (vic[0] != NULL) {
      std::vector<OBInternalCoord*>::iterator it = vic.begin();
      vic.insert(it, static_cast<OBInternalCoord*>(NULL));
    }

    if (vic.size() != mol.NumAtoms() + 1) {
      string error = "Number of internal coordinates is not the same as";
      error += " the number of atoms in molecule";
      obErrorLog.ThrowError(__FUNCTION__, error, obError);
      return;
    }

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

    for (atom = mol.BeginAtom(i);atom;atom = mol.NextAtom(i))
      {
        index = atom->GetIdx();

        // make sure we always have valid pointers
        if (index >= vic.size() || !vic[index])
          return;

        if (vic[index]->_a) // make sure we have a valid ptr
          {
            avec = vic[index]->_a->GetVector();
            dst = vic[index]->_dst;
          }
        else
          {
            // atom 1
            atom->SetVector(0.0, 0.0, 0.0);
            continue;
          }

        if (vic[index]->_b)
          {
            bvec = vic[index]->_b->GetVector();
            ang = vic[index]->_ang * DEG_TO_RAD;
          }
        else
          {
            // atom 2
            atom->SetVector(dst, 0.0, 0.0);
            continue;
          }

        if (vic[index]->_c)
          {
            cvec = vic[index]->_c->GetVector();
            tor = vic[index]->_tor * DEG_TO_RAD;
          }
        else
          {
            // atom 3
            cvec = VY;
            tor = 90. * DEG_TO_RAD;
          }

        v1 = avec - bvec;
        v2 = avec - cvec;
        n = cross(v1,v2);
        nn = cross(v1,n);
        n.normalize();
        nn.normalize();

        n  *= -sin(tor);
        nn *= cos(tor);
        v3 = n + nn;
        v3.normalize();
        v3 *= dst * sin(ang);
        v1.normalize();
        v1 *= dst * cos(ang);
        v2 = avec + v3 - v1;

        atom->SetVector(v2);
      }

    // Delete dummy atoms
    vector<OBAtom*> for_deletion;
    FOR_ATOMS_OF_MOL(a, mol)
      if (a->GetAtomicNum() == 0)
        for_deletion.push_back(&(*a));
    for(vector<OBAtom*>::iterator a_it=for_deletion.begin(); a_it!=for_deletion.end(); ++a_it)
//.........这里部分代码省略.........