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

  static unsigned int FindRingAtomsAndBonds2(OBMol &mol)
  {
    mol.SetRingAtomsAndBondsPerceived(); // mol.SetFlag(OB_RINGFLAGS_MOL);
    mol.SetClosureBondsPerceived();      // mol.SetFlag(OB_CLOSURE_MOL);

    unsigned int bsize = mol.NumBonds()+1;
    unsigned char *bvisit = (unsigned char*)malloc(bsize);
    memset(bvisit,0,bsize);

    unsigned int acount = mol.NumAtoms();
    unsigned int asize = (unsigned int)((acount+1)*sizeof(int));
    int *avisit = (int*)malloc(asize);
    memset(avisit,0,asize);

    unsigned int frj = 0;
    for(unsigned int i=1; i<=acount; i++ )
      if(avisit[i] == 0) {
        avisit[i] = 1;
        OBAtom *atom = mol.GetAtom(i);
        FindRings(atom,avisit,bvisit,frj,1);
      }
    free(avisit);
    free(bvisit);
    return frj;
  }

  static int DetermineFRJ(OBMol &mol)
  {
    vector<vector<int> >::iterator i;
    vector<vector<int> > cfl;
    //find all continuous graphs in the mol area
    mol.ContigFragList(cfl);

    if (cfl.empty())
      return(0);
    if (cfl.size() == 1)
      return(mol.NumBonds() - mol.NumAtoms() + 1);

    //count up the atoms and bonds belonging to each graph
    OBBond *bond;
    vector<OBBond*>::iterator j;
    int numatoms,numbonds,frj=0;
    OBBitVec frag;
    for (i = cfl.begin();i != cfl.end();++i)
      {
        frag.Clear();
        frag.FromVecInt(*i);
        numatoms = (*i).size();
        numbonds=0;
        for (bond = mol.BeginBond(j);bond;bond = mol.NextBond(j))
          if (frag.BitIsOn(bond->GetBeginAtomIdx()) &&
              frag.BitIsOn(bond->GetEndAtomIdx()))
            numbonds++;
        frj += numbonds - numatoms + 1;
      }

    return(frj);
  }

  void OutputUnions(ostream &ofs, OBMol &mol, string prefix)
  {
    /* ---- Build union of all atoms ---- */
    ofs << endl << "//All atoms of molecule " << prefix << endl;
    ofs << "#ifdef (TRANS)" << endl;
    ofs << "#declare " << prefix << "_atoms = merge {" << endl;
    ofs << "#else" << endl;
    ofs << "#declare " << prefix << "_atoms = union {" << endl;
    ofs << "#end //(End of TRANS)" << endl;

    /* ---- Write definition of all atoms ---- */
    for(unsigned int i = 1; i <= mol.NumAtoms(); ++i)
      {

        /* ---- Write definition of atom i ---- */
        ofs << "\t  object{" << prefix << "_atom" << i << "}" << endl;

      }
    ofs << "\t }" << endl << endl;

    /* ---- Check for number of bonds ---- */
    if(mol.NumBonds() > 0)
      {

        /* ---- Do a BAS or CST model ? ---- */
        ofs << "//Bonds only needed for ball and sticks or capped sticks models" << endl;
        ofs << "#if (BAS | CST)" << endl;
        ofs << "#declare " << prefix <<"_bonds = union {" << endl;

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

            /* ---- Write Definition of Bond i ---- */
            ofs << "\t  object{" << prefix << "_bond" << i << "}" << endl;

          }

        /* ---- End of povray-conditional for ball and sticks ---- */
        ofs << "\t }" << endl << "#end" << endl << endl;

      }

  }

bool WriteAlchemy(ostream &ofs,OBMol &mol)
{
  unsigned int i;
  char buffer[BUFF_SIZE];
  char bond_string[10];
  
  snprintf(buffer, BUFF_SIZE, "%5d ATOMS, %5d BONDS,     0 CHARGES",
	  mol.NumAtoms(),
	  mol.NumBonds());
  ofs << buffer << endl;
  ttab.SetFromType("INT"); ttab.SetToType("ALC");

  OBAtom *atom;
  string str,str1;
  for(i = 1;i <= mol.NumAtoms(); i++)
  {
    atom = mol.GetAtom(i);
    str = atom->GetType();
    ttab.Translate(str1,str);
    snprintf(buffer, BUFF_SIZE, "%5d %-6s%8.4f %8.4f %8.4f     0.0000",
	    i,
	    (char*)str1.c_str(),
	    atom->GetX(),
	    atom->GetY(),
	    atom->GetZ());
    ofs << buffer << endl;
  }

  OBBond *bond;
  vector<OBEdgeBase*>::iterator j;

  for (bond = mol.BeginBond(j);bond;bond = mol.NextBond(j))
  {
    switch(bond->GetBO())
    {
    case 1 :  strcpy(bond_string,"SINGLE");  break;
    case 2 :  strcpy(bond_string,"DOUBLE");  break;
    case 3 :  strcpy(bond_string,"TRIPLE");  break;
    case 5 :  strcpy(bond_string,"AROMATIC"); break;
    default : strcpy(bond_string,"SINGLE");
    }
    snprintf(buffer, BUFF_SIZE, "%5d  %4d  %4d  %s",
	    bond->GetIdx()+1,
	    bond->GetBeginAtomIdx(),
	    bond->GetEndAtomIdx(),
	    bond_string);
    ofs << buffer << endl;
  }
  return(true);
}

  void OutputHeader(ostream &ofs, OBMol &mol, string prefix)
  {
    time_t akttime;                              /* Systemtime                        */
    char timestr[TIME_STR_SIZE + 1] = "";        /* Timestring                        */
    size_t time_res;                             /* Result of strftime                */

    /* ---- Get the system-time ---- */
    akttime = time((time_t *) NULL);
    time_res = strftime(timestr,
                        TIME_STR_SIZE,
                        "%a %b %d %H:%M:%S %Z %Y",
                        localtime((time_t *) &akttime)
                        );

    /* ---- Write some header information ---- */
    ofs << "//Povray V3.1 code generated by Open Babel" << endl;
    ofs << "//Author: Steffen Reith <[email protected]>" << endl;

    /* ---- Include timestamp in header ---- */
    ofs << "//Date: " << timestr << endl << endl;

    /* ---- Include header statement for babel ---- */
    ofs << "//Include header for povray" << endl;
    ofs << "#include \"babel31.inc\"" << endl << endl;

    /* ---- You should do a spacefill model for molecules without bonds ---- */
    if (mol.NumBonds() == 0)
      {

        /* ---- Check if a spacefill-model is selected ---- */
        ofs << "#if (BAS | CST)\"" << endl;
        ofs << "#warning \"Molecule without bonds!\"" << endl;
        ofs << "#warning \"You should do a spacefill-model\"" << endl;
        ofs << "#end" << endl << endl;

      }

    /* ---- Set version ---- */
    ofs << "//Use PovRay3.1" << endl;
    ofs << "#version 3.1;" << endl << endl;

    /* ---- Print of name of molecule (#\b depends on size of babel.inc!) ---- */
    ofs << "//Print name of molecule while rendering" << endl;
    ofs << "#render \"\\b\\b " << mol.GetTitle() << "\\n\\n\"" << endl << endl;

  }

extern "C" unsigned int
ob_num_bonds (char *smiles)
{
  OBMol mol;
  OBConversion conv;
  string tmpStr (smiles);
  istringstream molstream (tmpStr);
  unsigned int numbonds = 0;

  conv.SetInAndOutFormats ("SMI", "SMI");

  conv.Read (&mol, &molstream);
  //mol.AddHydrogens (false, false);

  numbonds = mol.NumBonds ();

  return (numbonds);
}

void testIdsAddBond()
{
  OBMol mol;
  // add 5 bonds
  for (int i = 0; i < 5; ++i)
    mol.NewBond();

  OBBond bond;
  OBAtom *a = mol.NewAtom();
  OBAtom *b = mol.NewAtom();
  bond.SetBegin(a);
  bond.SetEnd(b);
  // add a sixth bond
  mol.AddBond(bond);

  OB_REQUIRE( mol.NumBonds() == 6 );
  OB_REQUIRE( mol.GetBondById(5) );
  OB_REQUIRE( mol.GetBondById(5)->GetId() == 5 );
}

  void FastSearchFormat::AddPattern(vector<OBMol>& patternMols, OBMol patternMol, int idx)
  {
    //Recursive function to generate all combinations of aromatic/single bonds for each tilde bond
    //Copying an OBMol, which happens when adding it to a vector, kekulizes it,
    // changing aromatic (bo=5) bonds. So set order after adding. Should work here,
    // but is dangerous if the vector needs to be reallocated.

    if(idx>=patternMol.NumBonds())
      return;
    if(patternMol.GetBond(idx)->GetBO()==4)
    {
      patternMol.GetBond(idx)->SetBO(1);
      patternMols.push_back(patternMol);
      AddPattern(patternMols, patternMol,idx+1);

      patternMols.push_back(patternMol);
      patternMols.back().GetBond(idx)->SetBO(5);
    }
    AddPattern(patternMols, patternMol,idx+1);
  }

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++;
        }
    }

}

  void OutputCSTBonds(ostream &ofs, OBMol &mol, string prefix)
  {
    /* ---- 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 ---- */
        ofs << "\t    pigment{color Color_"
            << bond -> GetBeginAtom() -> GetType()
            << "}" << 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;

        /* ---- End of description of Start-Bond ---- */
        ofs << "\t   }" << endl;

        /* ---- Begin of End-Half of Bond i ---- */
//.........这里部分代码省略.........

int main(int argc,char *argv[])
{
  // turn off slow sync with C-style output (we don't use it anyway).
  std::ios::sync_with_stdio(false);

  OBConversion conv;
  OBFormat *inFormat, *canFormat;
  OBMol mol;
  ifstream ifs;
  vector<OBMol> fragments;
  unsigned int fragmentCount = 0; // track how many in library -- give a running count
  map<string, int> index; // index of cansmi
  string currentCAN;
  unsigned int size;
  OBAtom *atom;
  OBBond *bond;
  bool nonRingAtoms, nonRingBonds;
  char buffer[BUFF_SIZE];

  canFormat = conv.FindFormat("can");
  conv.SetOutFormat(canFormat);

  if (argc < 2)
    {
      cout << "Usage: obfragment <file>" << endl;
      return(-1);
    }

  for (int i = 1; i < argc; i++) {
    cerr << " Reading file " << argv[i] << endl;

    inFormat = conv.FormatFromExt(argv[i]);
    if(inFormat==NULL || !conv.SetInFormat(inFormat))
      {
        cerr << " Cannot read file format for " << argv[i] << endl;
        continue; // try next file
      }
    
    ifs.open(argv[i]);
    
    if (!ifs)
      {
        cerr << "Cannot read input file: " << argv[i] << endl;
        continue;
      }
    
    
    while(ifs.peek() != EOF && ifs.good())
      {
        conv.Read(&mol, &ifs);
        if (!mol.Has3D()) continue; // invalid coordinates!
        mol.DeleteHydrogens(); // remove these before we do anything else
        
        do {
          nonRingAtoms = false;
          size = mol.NumAtoms();
          for (unsigned int i = 1; i <= size; ++i)
            {
              atom = mol.GetAtom(i);
              if (!atom->IsInRing()) {
                mol.DeleteAtom(atom);
                nonRingAtoms = true;
                break; // don't know how many atoms there are
              } 
              // Previously, we changed atoms to carbon here.
              // Now we perform this alchemy in terms of string-rewriting
              // once the canonical SMILES is generated
            }
        } while (nonRingAtoms);
        
        if (mol.NumAtoms() < 3)
          continue;
        
        if (mol.NumBonds() == 0)
          continue;
        
        do {
          nonRingBonds = false;
          size = mol.NumBonds();
          for (unsigned int i = 0; i < size; ++i)
            {
              bond = mol.GetBond(i);
              if (!bond->IsInRing()) {
                mol.DeleteBond(bond);
                nonRingBonds = true;
                break; // don't know how many bonds there are
              }
            }        
        } while (nonRingBonds);

        fragments = mol.Separate();
        for (unsigned int i = 0; i < fragments.size(); ++i)
          {
            if (fragments[i].NumAtoms() < 3) // too small to care
              continue;
              
            currentCAN = conv.WriteString(&fragments[i], true);
            currentCAN = RewriteSMILES(currentCAN); // change elements to "a/A" for compression
            if (index.find(currentCAN) != index.end()) { // already got this
              index[currentCAN] += 1; // add to the count for bookkeeping
//.........这里部分代码省略.........

  bool ACRFormat::ReadMolecule(OBBase* pOb, OBConversion* pConv)
  {
    OBMol* pmol = pOb->CastAndClear<OBMol>();
    if(pmol==NULL)
      return false;

    istream& ifs = *pConv->GetInStream();

    pmol->BeginModify();

    /** Parse the input stream and use the OpenBabel API to populate the OBMol **/
    int id;
    char buf[BUFF_SIZE];
    int atoms, bonds, tmp;
    float scale, dtmp;
    bool atom_input = false, bond_input = false;
    string type;
    //int from, to;
    double X,Y,Z;
    vector<string> vs;

    // read in one at a time
    /* WARNING: Atom id starts from zero in Carine; not so in openbabel.
     * Solution: increment atom id's */

    while (true) {
      ifs.getline(buf, BUFF_SIZE);
      if (ifs.eof()) {
        break;
      }

      if (sscanf(buf, "General Scale=%f\n", &dtmp)) {
        scale = dtmp;
        continue;
      } else if (sscanf(buf, "Number of Atoms in Crystal=%d\n", &tmp)) {
        atoms = tmp;
        atom_input = true;

        // read table column names
        ifs.getline(buf, BUFF_SIZE);
        continue;
      } else if (sscanf(buf, "Number of Links in Crystal=%d\n", &tmp)) {
        atom_input = false;
        bond_input = true;
        bonds = tmp;

        // read table column names
        ifs.getline(buf, BUFF_SIZE);
        continue;
      } else if ( '#' == buf[0] || '\r' == buf[0] || '\n' == buf[0] ) {
        // between sections, in both windows and unix.
        continue;
      }
      tokenize(vs, buf, " \t\r\n");

      if (atom_input) {
	if (vs.size() < 9) return false; // timvdm 18/06/2008
        id = atoi((char*)vs[0].c_str()) + 1; // see warning above
        type = vs[1];
        X = atof((char*)vs[6].c_str())/scale;
        Y = atof((char*)vs[7].c_str())/scale;
        Z = atof((char*)vs[8].c_str())/scale;

        OBAtom* a = pmol->NewAtom();
        if (*(type.c_str()) != '*')
          a->SetAtomicNum(etab.GetAtomicNum(type.c_str()));
        a->SetVector(X,Y,Z);

      } else if (bond_input) {
	if (vs.size() < 2) return false; // timvdm 18/06/2008
        // add to pmol
        if (!pmol->AddBond(atoi((char*)vs[0].c_str()) + 1, atoi((char*)vs[1].c_str()) + 1,
                           1 /* bond order not specified in Carine, use PerceiveBondOrder later */))
          {
            obErrorLog.ThrowError(__FUNCTION__, "addition of bond between " + vs[0] + " and " + vs[1] + " failed", obError);
            return false;
          }
      }
    }

    /* got sanity? */
    if ( pmol->NumBonds() != bonds ) {
      // then we read a different number of bonds than those promised.
      obErrorLog.ThrowError(__FUNCTION__, "Number of bonds read does not match the number promised", obError);
      return false;
    } else if ( pmol->NumAtoms() != atoms ) {
      obErrorLog.ThrowError(__FUNCTION__, "Number of atoms read does not match the number promised", obError);
      return false;
    }

    pmol->PerceiveBondOrders();

    pmol->EndModify();

    return true;
  }

//.........这里部分代码省略.........
      ifstream ifs;
      OBConversion conv;
      OBFormat *format_in = conv.FormatFromExt(vs[1].c_str());
   
      if (!format_in || !conv.SetInFormat(format_in)) {
        cout << "could not detect format." << endl;
        continue;
      }
       
      ifs.open(vs[1].c_str());
      if (!ifs) {
        cout << "could not open '" << vs[1] << "'." <<endl;
        continue;
      }
      
      mol.Clear();
      if (!conv.Read(&mol, &ifs)) {
        cout << "could not read a molecule from '" << vs[1] << "'." <<endl;
        continue;
      }
      
      if (mol.Empty()) {
        cout << "this molecule is empty." <<endl;
        continue;
      }

      if (!pFF->Setup(mol)) {
        cout << "error while initializing the force field for this molecule." <<endl;
        continue;
      }

      cout << "molecule succesfully loaded." << endl;
      cout << "  " << mol.NumAtoms() << " atoms" << endl;
      cout << "  " << mol.NumBonds() << " bonds" << endl;

      ifs.close();
 
      continue;
    }
    
    // save <filename>
    if (EQn(commandline, "save", 4)) {
      if (vs.size() < 2) {
        cout << "no <filename> specified." << endl;
        continue;
      }
      
      ofstream ofs;
      OBConversion conv;
      OBFormat *format_out = conv.FormatFromExt(vs[1].c_str());
   
      if (!format_out || !conv.SetOutFormat(format_out)) {
        cout << "could not detect format." << endl;
        continue;
      }
       
      ofs.open(vs[1].c_str());
      if (!ofs) {
        cout << "could not open '" << vs[1] << "'." <<endl;
        continue;
      }
      
      if (!conv.Write(&mol, &ofs)) {
        cout << "could not read a molecule from '" << vs[1] << "'." <<endl;
        continue;
      }

//.........这里部分代码省略.........
    OpenBabel::vector3 centroid = my_center_coords(mol.GetCoordinates(), mol.NumAtoms());

    ofs << "#include \"colors.inc\"\n" << endl;

    ofs <<  "// create a regular point light source\n"
            "light_source {\n"
            "  <" << centroid.x()  + 2.0 << "," << centroid.y() + 3.0 << "," << centroid.z() - 8.0 << ">\n"
            "  color rgb <1,1,1>    // light's color\n"
            "}\n" << endl;

    if (sky) {
      ofs << "// Add some nice sky with clouds\n"
"sky_sphere {\n"
"    pigment {\n"
"      gradient y\n"
"      color_map {\n"
"        [0.0 1.0 color SkyBlue  color NavyBlue]\n"
"      }\n"
"      scale 2\n"
"      translate -1\n"
"    }\n"
"    pigment {\n"
"      bozo\n"
"      turbulence 0.65\n"
"      octaves 6\n"
"      omega 0.7\n"
"      lambda 2\n"
"      color_map {\n"
"          [0.0 0.1 color rgb <0.85, 0.85, 0.85>\n"
"                   color rgb <0.75, 0.75, 0.75>]\n"
"          [0.1 0.5 color rgb <0.75, 0.75, 0.75>\n"
"                   color rgbt <1, 1, 1, 1>]\n"
"          [0.5 1.0 color rgbt <1, 1, 1, 1>\n"
"                   color rgbt <1, 1, 1, 1>]\n"
"      }\n"
"      scale <0.2, 0.5, 0.2>\n"
"    }\n"
"    rotate -135*x\n"
"  }\n" << endl;
    }
    else { // Simple background
      ofs << "// set a color of the background (sky)" << endl;
      ofs << "background { color rgb <0.95 0.95 0.95> }\n" << endl;
    }

    ofs <<  "// perspective (default) camera\n"
            "camera {\n"
            "  location  <" << centroid.x() << "," << centroid.y() << "," << centroid.z() - 10.0 << ">\n"
            "  look_at   <" << centroid.x() << "," << centroid.y() << "," << centroid.z() << ">\n"
            "  right     x*image_width/image_height\n"
            "}\n" << endl;

    /* ---- Checkerboard and mirror sphere ---- */
    if (sphere) {
      ofs <<
"// a mirror sphere\n"
"sphere\n"
"{ <" << centroid.x()  + 8.0 << "," << centroid.y() - 4 << "," << centroid.z() + 8.0 << ">,4\n"
"  pigment { rgb <0,0,0> } // A perfect mirror with no color\n"
"  finish { reflection 1 } // It reflects all\n"
"}\n" << endl;
    }
    if (checkerboard) {
      ofs <<
"// simple Black on White checkerboard... it's a classic\n"
"plane {\n"
" -y, " << -(centroid.y()-8) << "\n"
" pigment {\n"
"  checker color Black color White\n"
"  scale 2\n"
" }\n"
"}\n" << endl;
    }


    /* ---- Include header statement for babel ---- */
    ofs << "//Include header for povray" << endl;
    ofs << "#include \"babel_povray3.inc\"" << endl << endl;

    /* ---- You should do a spacefill model for molecules without bonds ---- */
    if (mol.NumBonds() == 0)
      {

        /* ---- Check if a spacefill-model is selected ---- */
        ofs << "#if (BAS | CST)\"" << endl;
        ofs << "#warning \"Molecule without bonds!\"" << endl;
        ofs << "#warning \"You should do a spacefill-model\"" << endl;
        ofs << "#end" << endl << endl;

      }

    /* ---- Set version ---- */
    ofs << "//Use PovRay3.6" << endl;
    ofs << "#version 3.6;" << endl << endl;

    /* ---- Print of name of molecule (#\b depends on size of babel.inc!) ---- */
    ofs << "//Print name of molecule while rendering" << endl;
    ofs << "#render \"\\b\\b " << mol.GetTitle() << "\\n\\n\"" << endl << 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;
  OBFormat *canSMIFormat = conv.FindFormat("can");
  OBFormat *inchiFormat = conv.FindFormat("inchi");


  ////////////////////////////////////////////////////////////////////////////
  // List of properties
  // Name
  // Molecular weight (Standard molar mass given by IUPAC atomic masses)
  // Number of rings : the size of the smallest set of smallest rings (SSSR)
  
  //.....ADD YOURS HERE.....
  
  for (c = 1;; ++c)
    {
      mol.Clear();
      conv.Read(&mol, &ifs);
      if (mol.Empty())
        break;
      
      if (!mol.HasHydrogensAdded())
        mol.AddHydrogens();
      // Print the properties
      if (strlen(mol.GetTitle()) != 0)
        cout << "name             " << mol.GetTitle() << endl;
      else 
        cout << "name             " << FileIn << " " << c << endl;

      cout << "formula          " << mol.GetFormula() << endl;
      cout << "mol_weight       " << mol.GetMolWt() << endl;
      cout << "exact_mass       " << mol.GetExactMass() << endl;

      string smilesString = "-";
      if (canSMIFormat) {
        conv.SetOutFormat(canSMIFormat);
        smilesString = conv.WriteString(&mol);
        if ( smilesString.length() == 0 )
        {
          smilesString = "-";
        }
      }
      cout << "canonical_SMILES " << smilesString << endl;

      string inchiString = "-";
      if (inchiFormat) {
        conv.SetOutFormat(inchiFormat);
        inchiString = conv.WriteString(&mol);
        if ( inchiString.length() == 0 )
        {
          inchiString = "-";
        }
      }
      cout << "InChI            " << inchiString << endl;

      cout << "num_atoms        " << mol.NumAtoms() << endl;
      cout << "num_bonds        " << mol.NumBonds() << endl;
      cout << "num_residues     " << mol.NumResidues() << endl;
      cout << "num_rotors       " << mol.NumRotors() << endl;
      if (mol.NumResidues() > 0)
        cout << "sequence         " << sequence(mol) << endl;
      else
        cout << "sequence         " << "-" << endl;

      cout << "num_rings        " << nrings(mol) << endl;

      OBDescriptor* pDesc;
      pDesc= OBDescriptor::FindType("logP");
      if(pDesc)
        cout << "logP             " << pDesc->Predict(&mol) << endl;

      pDesc = OBDescriptor::FindType("TPSA");
      if(pDesc)
        cout << "PSA              " << pDesc->Predict(&mol) << endl;

      pDesc = OBDescriptor::FindType("MR");
      if(pDesc)
        cout << "MR               " << pDesc->Predict(&mol) << endl;

      cout << "$$$$" << endl; // SDF like end of compound descriptor list
      
      //Other OBDescriptors could be output here, even ones that were rarely
      // used. Since these are plugin classes, they may not be loaded, but
      // then with code like the above they are just ignored.
    } // end for loop
  
  return(0);
}