C++ OBAtom::SetVector方法代码示例

  bool CHEM3D1Format::ReadChem3d(istream &ifs,OBMol &mol,bool mmads,const char *type_key)
  {
    char buffer[BUFF_SIZE];
    int natoms,i;
    char tmp[16],tmp1[16];
    char atomic_type[16];
    double exponent = 0.0;
    double divisor = 1.0;
    double Alpha,Beta,Gamma,A,B,C;
    bool has_fractional = false, has_divisor = false;
    matrix3x3 m;

    vector<string> vs;
    ifs.getline(buffer,BUFF_SIZE);
    tokenize(vs,buffer);

    if (mmads)
      {
        if (vs.empty())
          return(false);
        natoms = atoi((char*)vs[0].c_str());
        if (vs.size() == 2)
          mol.SetTitle(vs[1]);
      }
    else
      {
        switch(vs.size())
          {
          case 7 :
            sscanf(buffer,"%d%lf%lf%lf%lf%lf%lf",
                   &natoms,&Alpha,&Beta,&Gamma,&A,&B,&C);
            m.FillOrth(Alpha,Beta,Gamma,A,B,C);
            has_fractional = true;
            break;
          case 8 :
            sscanf(buffer,"%d%lf%lf%lf%lf%lf%lf%lf",
                   &natoms,&Alpha,&Beta,&Gamma,&A,&B,&C,&exponent);
            m.FillOrth(Alpha,Beta,Gamma,A,B,C);
            has_fractional = true;
            has_divisor = true;
            break;
          default :
            sscanf(buffer,"%d",&natoms);
            break;
          }
      }

    if (!natoms)
      return(false);
    divisor = pow(10.0,exponent);
    mol.ReserveAtoms(natoms);

    ttab.SetToType("INT");
    ttab.SetFromType(type_key);

    OBAtom *atom;
    double x,y,z;
    vector3 v;

    unsigned int k;
    for (i = 1; i <= natoms; i++)
      {
        ifs.getline(buffer,BUFF_SIZE);
        sscanf(buffer,"%15s%*d%lf%lf%lf%15s",
               atomic_type,
               &x,
               &y,
               &z,
               tmp);
        v.Set(x,y,z);
        if (has_fractional)
          v *= m;
        if (has_divisor)
          v/= divisor;

        tokenize(vs,buffer);
        if (vs.empty())
          return(false);

        atom = mol.NewAtom();
        ttab.Translate(tmp1,tmp);
        atom->SetType(tmp1);
        atom->SetVector(v);
        atom->SetAtomicNum(etab.GetAtomicNum(atomic_type));

        for (k = 6;k < vs.size(); k++)
          mol.AddBond(atom->GetIdx(),atoi((char*)vs[k].c_str()),1);
      }

    // clean out remaining blank lines
    while(ifs.peek() != EOF && ifs.good() &&
          (ifs.peek() == '\n' || ifs.peek() == '\r'))
      ifs.getline(buffer,BUFF_SIZE);

    mol.PerceiveBondOrders();

    return(true);
  }

//.........这里部分代码省略.........
          } else
            if (resname == "ADR" || resname == "COA" || resname == "FAD" ||
                resname == "GPG" || resname == "NAD" || resname == "NAL" ||
                resname == "NDP" || resname == "ABA") {
              if (type.size() > 1)
                type = type.substr(0,1);
              //type.erase(1,type.size()-1);
            } else // other residues
              if (isdigit(type[0])){
                type = type.substr(1,1);
              }
              else
                if (type.size() > 1 && isdigit(type[1]))
                  type = type.substr(0,1);
                else
                  if (type.size() > 1 && isalpha(type[1])) {
                    if (type[0] == 'O' && type[1] == 'H')
                      type = type.substr(0,1); // no "Oh" element (e.g. 1MBN)
                    else if(isupper(type[1])) {
                      type[1] = tolower(type[1]);
                    }
                  }
        }

      } // HETATM records
    } // no element column to use

    OBAtom atom;
    /* X, Y, Z */
    string xstr = sbuf.substr(24,8);
    string ystr = sbuf.substr(32,8);
    string zstr = sbuf.substr(40,8);
    vector3 v(atof(xstr.c_str()),atof(ystr.c_str()),atof(zstr.c_str()));
    atom.SetVector(v);
    atom.ForceImplH();

    // useful for debugging unknown atom types (e.g., PR#1577238)
    //    cout << mol.NumAtoms() + 1  << " : '" << element << "'" << " " << etab.GetAtomicNum(element.c_str()) << endl;
    if (elementFound)
      atom.SetAtomicNum(etab.GetAtomicNum(element.c_str()));
    else // use our old-style guess from athe atom type
      atom.SetAtomicNum(etab.GetAtomicNum(type.c_str()));

    if ( (! scharge.empty()) && "  " != scharge )
      {
        if ( isdigit(scharge[0]) && ('+' == scharge[1] || '-' == scharge[1]) )
          {
            const char reorderCharge[3] = { scharge[1], scharge[0], '\0' };
            const int charge = atoi(reorderCharge);
            atom.SetFormalCharge(charge);
          }
        else
          {
            stringstream errorMsg;
            errorMsg << "WARNING: Problems reading a PDB file\n"
                     << "  Problems reading a HETATM or ATOM record.\n"
                     << "  According to the PDB specification,\n"
                     << "  columns 79-80 should contain charge of the atom\n"
                     << "  but OpenBabel found '" << scharge << "' (atom " << mol.NumAtoms()+1 << ").";
            obErrorLog.ThrowError(__FUNCTION__, errorMsg.str(), obWarning);
          }
      }
    else {
      atom.SetFormalCharge(0);
    }

  //! 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)
//.........这里部分代码省略.........

  void SuperCellExtension::fillCell()
  {
    /* Change coords back to inverse space, apply the space group transforms
     *  then change coords back to real space
     */
    if (!m_molecule)
      return;

    OBUnitCell *uc = m_molecule->OBUnitCell();
    if (!uc) {
      qDebug() << "No unit cell found - fillCell() returning...";
      return;
    }

    const SpaceGroup *sg = uc->GetSpaceGroup(); // the actual space group and transformations for this unit cell
    if (sg) {
      qDebug() << "Space group:" << sg->GetId();// << sg->GetHMName();
      // We operate on a copy of the Avogadro molecule
      // For each atom, we loop through:
      // * convert the coords back to inverse space
      // * apply the transformations
      // * create new (duplicate) atoms
      OBMol mol = m_molecule->OBMol();
      vector3 uniqueV, newV;
      list<vector3> transformedVectors; // list of symmetry-defined copies of the atom
      list<vector3>::iterator transformIterator, duplicateIterator;
      vector3 updatedCoordinate;
      bool foundDuplicate;

      OBAtom *addAtom;
      QList<OBAtom*> atoms; // keep the current list of unique atoms -- don't double-create
      list<vector3> coordinates; // all coordinates to prevent duplicates
      FOR_ATOMS_OF_MOL(atom, mol)
        atoms.push_back(&(*atom));

      foreach(OBAtom *atom, atoms) {
        uniqueV = atom->GetVector();
        // Assert: won't crash because we already ensure uc != NULL
        uniqueV = uc->CartesianToFractional(uniqueV);
        uniqueV = transformedFractionalCoordinate(uniqueV);
        coordinates.push_back(uniqueV);

        transformedVectors = sg->Transform(uniqueV);
        for (transformIterator = transformedVectors.begin();
             transformIterator != transformedVectors.end(); ++transformIterator) {
          // coordinates are in reciprocal space -- check if it's in the unit cell
          // if not, transform it in place
          updatedCoordinate = transformedFractionalCoordinate(*transformIterator);
          foundDuplicate = false;

          // Check if the transformed coordinate is a duplicate of an atom
          for (duplicateIterator = coordinates.begin();
               duplicateIterator != coordinates.end(); ++duplicateIterator) {
            if (duplicateIterator->distSq(updatedCoordinate) < 1.0e-4) {
              foundDuplicate = true;
              break;
            }
          }
          if (foundDuplicate)
            continue;

          coordinates.push_back(updatedCoordinate); // make sure to check the new atom for dupes
          addAtom = mol.NewAtom();
          addAtom->Duplicate(atom);
          addAtom->SetVector(uc->FractionalToCartesian(updatedCoordinate));
        } // end loop of transformed atoms

        // Put the original atom into the proper space in the unit cell too
        atom->SetVector(uc->FractionalToCartesian(uniqueV));
      } // end loop of atoms

      m_molecule->setOBMol(&mol);
      qDebug() << "Spacegroups done...";

      // Need a fresh pointer to the new unit cell - setOBMol is invalidating
      // the old one. This should be cleaned up to use a more permanent data
      // structure.
      uc = m_molecule->OBUnitCell();
      uc->SetSpaceGroup(1);
    }

  //! 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)
//.........这里部分代码省略.........

  bool BGFFormat::ReadMolecule(OBBase* pOb, OBConversion* pConv)
  {

    OBMol* pmol = pOb->CastAndClear<OBMol>();
    if(pmol==NULL)
      return false;

    //Define some references so we can use the old parameter names
    istream &ifs = *pConv->GetInStream();
    OBMol &mol = *pmol;
    mol.SetTitle( pConv->GetTitle()); //default title is the filename
    mol.BeginModify();

    char buffer[BUFF_SIZE];
    char tmp[16],tmptyp[16];

    while (ifs.getline(buffer,BUFF_SIZE))
      if (EQn(buffer,"FORMAT",6))
        break;

    ttab.SetFromType("DRE");
    ttab.SetToType("INT");
    OBAtom *atom;
    double x,y,z,chrg;
    for (;;)
      {
        if (!ifs.getline(buffer,BUFF_SIZE))
          break;
        if (EQn(buffer,"FORMAT",6))
          break;

        sscanf(buffer,"%*s %*s %*s %*s %*s %*s %lf %lf %lf %15s %*s %*s %lf",
               &x,&y,&z,
               tmptyp,
               &chrg);
        atom = mol.NewAtom();

        ttab.Translate(tmp,tmptyp);
        atom->SetType(tmp);

        CleanAtomType(tmptyp);
        atom->SetAtomicNum(etab.GetAtomicNum(tmptyp));

        atom->SetVector(x,y,z);
      }
    unsigned int i;
    vector<int> vtmp;
    vector<vector<int> > vcon;
    vector<vector<int> > vord;

    for (i = 0; i < mol.NumAtoms();i++)
      {
        vcon.push_back(vtmp);
        vord.push_back(vtmp);
      }

    unsigned int bgn;
    vector<string> vs;
    for (;;)
      {
        if (!ifs.getline(buffer,BUFF_SIZE) || EQn(buffer,"END",3))
          break;

        tokenize(vs,buffer);
        if (vs.empty() || vs.size() < 3 || vs.size() > 10)
          continue;

        if (EQn(buffer,"CONECT",6))
          {
            bgn = atoi((char*)vs[1].c_str()) - 1;
            if (bgn < 1 || bgn > mol.NumAtoms())
              continue;
            for (i = 2;i < vs.size();i++)
              {
                vcon[bgn].push_back(atoi((char*)vs[i].c_str()));
                vord[bgn].push_back(1);
              }
          }
        else
          if (EQn(buffer,"ORDER",5))
            {
              bgn = atoi((char*)vs[1].c_str()) - 1;
              if (bgn < 1 || bgn > mol.NumAtoms())
                continue;
              if (vs.size() > vord[bgn].size()+2)
                continue;
              for (i = 2;i < vs.size();i++)
                vord[bgn][i-2] = atoi((char*)vs[i].c_str());
            }
      }

    unsigned int j;
    for (i = 1;i <= mol.NumAtoms();i++)
      if (!vcon[i - 1].empty())
        for (j = 0;j < vcon[i - 1].size();j++)
          {
            mol.AddBond(i,vcon[i - 1][j],vord[i - 1][j]);
          }

    //load up the next line after the END marker
//.........这里部分代码省略.........

bool CCCFormat::ReadMolecule(OBBase* pOb, OBConversion* pConv)
{

    OBMol* pmol = pOb->CastAndClear<OBMol>();
    if(pmol==NULL)
        return false;

    //Define some references so we can use the old parameter names
    istream &ifs = *pConv->GetInStream();
    OBMol &mol = *pmol;
    mol.SetTitle( pConv->GetTitle()); //default title is the filename

    char buffer[BUFF_SIZE];
    ifs.getline(buffer,BUFF_SIZE);

    if (strlen(buffer) > 5)
        mol.SetTitle(&buffer[5]);
    mol.SetEnergy(0.0);

    int natoms;
    ifs.getline(buffer,BUFF_SIZE);
    sscanf(buffer,"%*s%d",&natoms);
    mol.ReserveAtoms(natoms);
    mol.BeginModify();

    int end,order;
    double x,y,z;
    OBAtom atom;
    vector3 v;
    vector<string> vs;
    char element[3];
    element[2] = '\0';

    for (int i = 1;i <= natoms;i++)
    {
        if (!ifs.getline(buffer,BUFF_SIZE))
            return(false);
        atom.Clear();
        element[0] = buffer[0];
        element[1] = (buffer[1] != ' ') ? buffer[1]:'\0';
        atom.SetAtomicNum(etab.GetAtomicNum(element));
        sscanf(&buffer[15],"%lf%lf%lf",&x,&y,&z);
        v.Set(x,y,z);
        atom.SetVector(v);

        if (!mol.AddAtom(atom))
            return(false);
        tokenize(vs,&buffer[60]);
        vector<string>::iterator j;

        for (j = vs.begin();j != vs.end();j++)
            if (!j->empty())
            {
                //get the bond order
                switch((char)(*j)[j->size()-1])
                {
                case 'S':
                    order = 1;
                    break;
                case 'D':
                    order = 2;
                    break;
                case 'T':
                    order = 3;
                    break;
                default:
                    order = 1;
                }
                (*j)[j->size()-1] = ' ';
                end = atoi(j->c_str());
                if (i>end)
                    mol.AddBond(i,end,order);
            }
    }

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

//.........这里部分代码省略.........
        FOR_ATOMS_OF_MOL(a, *pmol)
          toDelete.push_back(&*a);
        for (size_t i = 0; i < toDelete.size(); i++)
          pmol->DeleteAtom(toDelete.at(i));


        // Discover units
        matrix3x3 conv (1);
        tokenize(vs, buffer);

        if (strstr(vs[1].c_str(), "alat")) {
          conv *= (alat * BOHR_TO_ANGSTROM);
        }
        else if (strstr(vs[1].c_str(), "crystal")) {
          // Set to the zero matrix and test below.
          conv = matrix3x3 (0.0);
        }
        // Add others if needed

        // Load new atoms from molecule
        ifs.getline(buffer,BUFF_SIZE); // First entry
        tokenize(vs, buffer);
        int size = vs.size();
        while (size == 4) {
          atomicNum = OBElements::GetAtomicNum(vs[0].c_str());
          x = atof((char*)vs[1].c_str());
          y = atof((char*)vs[2].c_str());
          z = atof((char*)vs[3].c_str());
          // Add atom
          OBAtom *atom = pmol->NewAtom();
          atom->SetAtomicNum(atomicNum);
          vector3 coords (x,y,z);
          if (conv.determinant() == 0.0) { // Fractional coords
            atom->SetVector(cell->FractionalToCartesian(coords));
          }
          else {
            atom->SetVector(conv * coords);
          }

          // Reset vars
          ifs.getline(buffer,BUFF_SIZE); // First entry
          tokenize(vs, buffer);
          size = vs.size();
        }
      }

      // Free energy
      if (strstr(buffer, "Final energy =")) {
        tokenize(vs, buffer);
        pmol->SetEnergy(atof(vs[3].c_str()) * RYDBERG_TO_KCAL_PER_MOL);
      }

      // H - PV = U energy
      if (strstr(buffer, "!    total energy              =")) {
        tokenize(vs, buffer);
        pmol->SetEnergy(atof(vs[4].c_str()) * RYDBERG_TO_KCAL_PER_MOL);
      }

      // Enthalphy
      if (strstr(buffer, "Final enthalpy =")) {
        tokenize(vs, buffer);

        hasEnthalpy = true;
        enthalpy = atof(vs.at(3).c_str()) * RYDBERG_TO_KCAL_PER_MOL;
        pv = enthalpy - pmol->GetEnergy();
      }

//.........这里部分代码省略.........
                if (comment != NULL)
                  delete [] comment;
                comment = new char [len];
                memcpy(comment,buffer,len);
              }
          }
      }

    if (!foundAtomLine)
      {
        mol.EndModify();
        mol.Clear();
        obErrorLog.ThrowError(__FUNCTION__, "Unable to read Mol2 format file. No atoms found.", obWarning);
        return(false);
      }

    mol.ReserveAtoms(natoms);

    int i;
    vector3 v;
    OBAtom atom;
    double x,y,z,pcharge;
    char temp_type[BUFF_SIZE], resname[BUFF_SIZE], atmid[BUFF_SIZE];
    int elemno, resnum = -1;

    ttab.SetFromType("SYB");
    for (i = 0;i < natoms;i++)
      {
        if (!ifs.getline(buffer,BUFF_SIZE))
          return(false);
        sscanf(buffer," %*s %1024s %lf %lf %lf %1024s %d %1024s %lf",
               atmid, &x,&y,&z, temp_type, &resnum, resname, &pcharge);

        atom.SetVector(x, y, z);

        // Handle "CL" and "BR" and other mis-typed atoms
        str = temp_type;
        if (strncmp(temp_type, "CL", 2) == 0) {
          str = "Cl";
        } else  if (strncmp(temp_type,"BR",2) == 0) {
          str = "Br";
        } else if (strncmp(temp_type,"S.o2", 4) == 02) {
          str = "S.O2";
        } else if (strncmp(temp_type,"S.o", 3) == 0) {
          str = "S.O";
        } else if (strncmp(temp_type,"SI", 2) == 0) {
          str = "Si";
        // The following cases are entries which are not in openbabel/data/types.txt
        // and should probably be added there
        } else if (strncmp(temp_type,"S.1", 3) == 0) {
          str = "S.2"; // no idea what the best type might be here
        } else if (strncmp(temp_type,"P.", 2) == 0) {
          str = "P.3";
        } else if (strncasecmp(temp_type,"Ti.", 3) == 0) { // e.g. Ti.th
          str = "Ti";
        } else if (strncasecmp(temp_type,"Ru.", 3) == 0) { // e.g. Ru.oh
          str = "Ru";
        }

        ttab.SetToType("ATN");
        ttab.Translate(str1,str);
        elemno = atoi(str1.c_str());
        ttab.SetToType("IDX");

        // We might have missed some SI or FE type things above, so here's
        // another check

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

        if (strstr(buffer, "SpaceGroup") != NULL) {
          tokenize(vs, buffer);
          if (vs.size() != 5)
            continue; // invalid space group
          setSpaceGroup = true;
          sg = SpaceGroup::GetSpaceGroup(vs[4]); // remove the initial " character
        }

        // atom information
        if (atomRecord) {
          if (strstr(buffer, "ACL") != NULL) {
            tokenize(vs, buffer);
            // size should be 5 -- need a test here
            if (vs.size() != 5) return false; // timvdm 18/06/2008
            vs[3].erase(0,1); // "6 => remove the first " character
            unsigned int atomicNum = atoi(vs[3].c_str());
            if (atomicNum == 0)
              atomicNum = 1; // hydrogen ?

            // valid element, so create the atom
            atom = mol.NewAtom();
            atom->SetAtomicNum(atomicNum);
            continue;
          }
          else if (strstr(buffer, "XYZ") != NULL) {
            tokenize(vs, buffer);
            // size should be 6 -- need a test here
            if (vs.size() != 6) return false; // timvdm 18/06/2008
            vs[3].erase(0,1); // remove ( character
            vs[5].erase(vs[5].length()-2, 2); // remove trailing )) characters
            atom->SetVector(atof(vs[3].c_str()),
                            atof(vs[4].c_str()),
                            atof(vs[5].c_str()));
            continue;
          }
        } // end of atom records

        // bond information
        if (bondRecord) {
          if (strstr(buffer, "Atom1") != NULL) {
            tokenize(vs, buffer);
            if (vs.size() < 4) return false; // timvdm 18/06/2008
            vs[3].erase(vs[3].length()-1,1);
            startBondAtom = atoi(vs[3].c_str());
            continue;
          }
          else if (strstr(buffer, "Atom2") != NULL) {
            tokenize(vs, buffer);
            if (vs.size() < 4) return false; // timvdm 18/06/2008
            vs[3].erase(vs[3].length()-1,1);
            endBondAtom = atoi(vs[3].c_str());
            continue;
          }
          else if (strstr(buffer, "Type") != NULL) {
            tokenize(vs, buffer);
            if (vs.size() < 4) return false; // timvdm 18/06/2008
            vs[3].erase(vs[3].length()-1,1);
            bondOrder = atoi(vs[3].c_str());
            if (bondOrder == 4) // triple bond?
              bondOrder = 3;
            else if (bondOrder == 8) // aromatic?
              bondOrder = 5;
            else if (bondOrder != 2) // 1 OK, 2 OK, others unknown

  bool MacroModFormat::ReadMolecule(OBBase* pOb, OBConversion* pConv)
  {

    OBMol* pmol = pOb->CastAndClear<OBMol>();
    if(pmol==NULL)
      return false;

    //Define some references so we can use the old parameter names
    istream &ifs = *pConv->GetInStream();
    OBMol &mol = *pmol;
    const char* defaultTitle = pConv->GetTitle();

    // Get Title
    char buffer[BUFF_SIZE];
    int natoms;
    vector<vector<pair<int,int> > > connections;

    if (ifs.getline(buffer,BUFF_SIZE))
      {
        vector<string> vs;
        tokenize(vs,buffer," \n");

        if ( !vs.empty() && vs.size() > 0)
          sscanf(buffer,"%i%*s",&natoms);

        if (natoms == 0)
          return false;

        if ( !vs.empty() && vs.size() > 1)
          mol.SetTitle(vs[1]);
        else
          {
            string s = defaultTitle;
            mol.SetTitle(defaultTitle);
          }
      }
    else
      return(false);

    mol.BeginModify();
    mol.ReserveAtoms(natoms);
    connections.resize(natoms+1);

    /***********************************************************************/

    // Get Type Bonds, BondOrder, X, Y, Z

    double x,y,z;
    vector3 v;
    char temp_type[10];
    int i,j;
    double charge;
    OBAtom atom;

    ttab.SetFromType("MMD");
    for (i = 1; i <= natoms; i++)
      {
        if (!ifs.getline(buffer,BUFF_SIZE))
          break;

        int end[6], order[6];

        sscanf(buffer,"%9s%d%d%d%d%d%d%d%d%d%d%d%d%lf%lf%lf",
               temp_type,&end[0],&order[0],&end[1],&order[1],&end[2],&order[2],
               &end[3], &order[3], &end[4], &order[4], &end[5], &order[5],
               &x, &y, &z);

        pair<int,int> tmp;
        for ( j = 0 ; j <=5 ; j++ )
          {
            if ( end[j] > 0  && end[j] > i)
              {
                tmp.first = end[j];
                tmp.second = order[j];
                connections[i].push_back(tmp);
              }
          }

        v.SetX(x);
        v.SetY(y);
        v.SetZ(z);
        atom.SetVector(v);

        string str = temp_type,str1;
        ttab.SetToType("ATN");
        ttab.Translate(str1,str);
        atom.SetAtomicNum(atoi(str1.c_str()));
        ttab.SetToType("INT");
        ttab.Translate(str1,str);
        atom.SetType(str1);

        // stuff for optional fields

        buffer[109]='\0';
        sscanf(&buffer[101],"%lf", &charge);
        atom.SetPartialCharge(charge);
        mol.AddAtom(atom);
      }

    for (i = 1; i <= natoms; i++)
//.........这里部分代码省略.........

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

    //Define some references so we can use the old parameter names
    istream &ifs = *pConv->GetInStream();
    OBMol &mol = *pmol;
    double UnitConv=AAU;
    if(pConv->IsOption("a", OBConversion::INOPTIONS))
      UnitConv=1;


    char buffer[BUFF_SIZE];
    do
    {
        ifs.getline(buffer,BUFF_SIZE);
	if (ifs.peek() == EOF || !ifs.good())
	  return false;
    }
    while(strncmp(buffer,"$coord",6));

    mol.BeginModify();
    OBAtom atom;
    while(!(!ifs))
    {
        ifs.getline(buffer,BUFF_SIZE);
        if(*buffer=='$')
            break;
        if(*buffer=='#')
            continue;
        float x,y,z;
        char atomtype[8];
        if(sscanf(buffer,"%f %f %f %7s",&x,&y,&z,atomtype)!=4)
            return false;

        atom.SetVector(x*UnitConv, y*UnitConv, z*UnitConv);
        atom.SetAtomicNum(OBElements::GetAtomicNum(atomtype));
        atom.SetType(atomtype);

        if(!mol.AddAtom(atom))
            return false;
        atom.Clear();
    }
    while(!(!ifs) && strncmp(buffer,"$end",4))
        ifs.getline(buffer,BUFF_SIZE);

    if (!pConv->IsOption("b",OBConversion::INOPTIONS))
      mol.ConnectTheDots();
    if (!pConv->IsOption("s",OBConversion::INOPTIONS) && !pConv->IsOption("b",OBConversion::INOPTIONS))
      mol.PerceiveBondOrders();

    // clean out remaining blank lines
    std::streampos ipos;
    do
    {
      ipos = ifs.tellg();
      ifs.getline(buffer,BUFF_SIZE);
    }
    while(strlen(buffer) == 0 && !ifs.eof() );
    ifs.seekg(ipos);

    mol.EndModify();
    return true;
}

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

    d->mol = mol;

    OBAtom *atom;
    OBBondIterator j;
    OBAtomIterator i;

    // 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);

    double width = max_x - min_x + 2*margin;
    double 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
    for (OBBond *bond = mol->BeginBond(j); bond; bond = mol->NextBond(j)) {
      OBAtom *begin = bond->GetBeginAtom();
      OBAtom *end = bond->GetEndAtom();

      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;
        d->DrawRingBond(begin, end, center, ringBond->GetBO());
        drawnBonds.SetBitOn(ringBond->GetId());
      }

    }

    // draw atom labels
    for (atom = mol->BeginAtom(i); atom; atom = mol->NextAtom(i)) {
      double x = atom->GetX();
      double y = atom->GetY();

      int alignment = GetLabelAlignment(atom);
      bool rightAligned = false;
      switch (alignment) {
        case TopRight:
        case CenterRight:
        case BottomRight:
//.........这里部分代码省略.........

int mol(int argc, char* argv[])
{
  int defaultchoice = 1;
  
  int choice = defaultchoice;

  if (argc > 1) {
    if(sscanf(argv[1], "%d", &choice) != 1) {
      printf("Couldn't parse that input as a number\n");
      return -1;
    }
  }


  // Define location of file formats for testing
  #ifdef FORMATDIR
    char env[BUFF_SIZE];
    snprintf(env, BUFF_SIZE, "BABEL_LIBDIR=%s", FORMATDIR);
    putenv(env);
  #endif

  cout << "# Unit tests for OBMol \n";

  cout << "ok 1\n"; // for loading tests

  OBMol emptyMol, testMol1;
  cout << "ok 2\n"; // ctor works

  testMol1.ReserveAtoms(-1);
  testMol1.ReserveAtoms(0);
  testMol1.ReserveAtoms(2);
  cout << "ok 3\n";

  // atom component tests

  if (testMol1.NumAtoms() == 0) {
    cout << "ok 4\n";
  } else {
    cout << "not ok 4\n";
  }

  testMol1.NewAtom();
  if (testMol1.NumAtoms() == 1) {
    cout << "ok 5\n";
  } else {
    cout << "not ok 5\n";
  }

  testMol1.NewAtom();
  testMol1.AddBond(1, 2, 1);
  if (testMol1.NumBonds() == 1) {
    cout << "ok 6\n";
  } else {
    cout << "not ok 6\n";
  }

  testMol1.Clear();
  if (testMol1.NumAtoms() == 0) {
    cout << "ok 7\n";
  } else {
    cout << "not ok 7\n";
  }

  ifstream ifs1(kd3file.c_str());
  if (!ifs1)
    {
      cout << "Bail out! Cannot read input file!" << endl;
      return(-1);
    }
  OBConversion conv(&ifs1, &cout);
  OBFormat* pFormat;
  
  pFormat = conv.FindFormat("XYZ");
  if ( pFormat == NULL )
    {
      cout << "Bail out! Cannot read file format!" << endl;
      return(-1);
    }
  if (! conv.SetInAndOutFormats(pFormat, pFormat))
    {
      cout << "Bail out! File format isn't loaded" << endl;
      return (-1);
    }

  OBMol testMol2D, testMol3D;
  if (conv.Read(&testMol3D))
    cout << "ok 8\n";
  else
    cout << "not ok 8\n";
  testMol3D.Center();
  
  // test bond insertion (PR#1665649)
  OBMol doubleBondMol;
  OBAtom *a1, *a2;
  OBBond *b;
  doubleBondMol.BeginModify();
  a1 = doubleBondMol.NewAtom();
  a1->SetVector(0.0, 0.0, 0.0);
  a1->SetAtomicNum(6);
  a2 = doubleBondMol.NewAtom();
//.........这里部分代码省略.........

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

    //Define some references so we can use the old parameter names
    istream &ifs = *pConv->GetInStream();
    OBMol &mol = *pmol;
    const char* title = pConv->GetTitle();

    int natoms;
    char buffer[BUFF_SIZE];
    vector<string> vs;

    ifs.getline(buffer, BUFF_SIZE);
    tokenize(vs,buffer);
    if (vs.size() < 2)
        return false;
    natoms = atoi(vs[0].c_str());

    // title is 2nd token (usually add tokens for the atom types)
    mol.SetTitle(vs[1]);

    mol.ReserveAtoms(natoms);
    mol.BeginModify();

    string str;
    double x,y,z;
    OBAtom *atom;
    int atomicNum;

    for (int i = 1; i <= natoms; ++i)
    {
        if (!ifs.getline(buffer,BUFF_SIZE))
            return(false);
        tokenize(vs,buffer);
        // e.g. "2  C      2.476285    0.121331   -0.001070     2     1     3    14"
        if (vs.size() < 5)
            return(false);

        atom = mol.NewAtom();
        x = atof((char*)vs[2].c_str());
        y = atof((char*)vs[3].c_str());
        z = atof((char*)vs[4].c_str());
        atom->SetVector(x,y,z); //set coordinates

        //set atomic number
        atom->SetAtomicNum(etab.GetAtomicNum(vs[1].c_str()));

        // add bonding
        if (vs.size() > 6)
            for (int j = 6; j < vs.size(); ++j)
                mol.AddBond(mol.NumAtoms(), atoi((char *)vs[j].c_str()), 1); // we don't know the bond order

    }
    if (!pConv->IsOption("s",OBConversion::INOPTIONS))
        mol.PerceiveBondOrders();

    // clean out remaining blank lines
    while(ifs.peek() != EOF && ifs.good() &&
            (ifs.peek() == '\n' || ifs.peek() == '\r'))
        ifs.getline(buffer,BUFF_SIZE);

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