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

 void OBChargeModel::FillChargeVectors(OBMol &mol)
 {
   OBAtom *atom;
   vector<OBAtom*>::iterator itr;
   m_partialCharges.clear();
   m_partialCharges.reserve(mol.NumAtoms());
   m_formalCharges.clear();
   m_formalCharges.reserve(mol.NumAtoms());
   
   for (atom = mol.BeginAtom(itr);atom;atom = mol.NextAtom(itr))
     {
       m_partialCharges.push_back(atom->GetPartialCharge());
       m_formalCharges.push_back((double)(atom->GetFormalCharge()));
     }
 }

//.........这里部分代码省略.........
                  {
                    char tmp[10];
                    strncpy(tmp, type_name, 9); // make sure to null-terminate tmp
                    snprintf(padded_name, sizeof(padded_name), " %-3s", tmp);
                    strncpy(type_name,padded_name,4);
                    type_name[4] = '\0';
                  }
                else
                  {
                    /*
                      type_name[4] = type_name[3];
                      type_name[3] = type_name[2];
                      type_name[2] = type_name[1];
                      type_name[1] = type_name[0];
                      type_name[0] = type_name[4];
                    */
                    type_name[4] = '\0';
                  }
              }
            res_num = res->GetNum();
          }
        else
          {
            strcpy(the_res,"UNK");
            the_res[3] = '\0';
            snprintf(padded_name,sizeof(padded_name), "%s",type_name);
            strncpy(type_name,padded_name,4);
            type_name[4] = '\0';
            res_num = 1;
          }

        element_name = etab.GetSymbol(atom->GetAtomicNum());

        int charge = atom->GetFormalCharge();
        char scharge[3] = { ' ', ' ', '\0' };
        if(0 != charge)
          {
            snprintf(scharge, 3, "%+d", charge);
            char tmp = scharge[1];
            scharge[1] = scharge[0];
            scharge[0] = tmp;
          }
        snprintf(buffer, BUFF_SIZE, "%s%5d %-4s %-3s %c%4d    %8.3f%8.3f%8.3f  1.00  0.00          %2s%2s\n",
                 het?"HETATM":"ATOM  ",
                 i,
                 type_name,
                 the_res,
                 the_chain,
                 res_num,
                 atom->GetX(),
                 atom->GetY(),
                 atom->GetZ(),
                 element_name,
                 scharge);
        ofs << buffer;
      }

    OBAtom *nbr;
    vector<OBBond*>::iterator k;
    for (i = 1; i <= mol.NumAtoms(); i ++)
      {
        atom = mol.GetAtom(i);
        if (atom->GetValence() == 0)
          continue; // no need to write a CONECT record -- no bonds

        snprintf(buffer, BUFF_SIZE, "CONECT%5d", i);

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    xmlTextWriterStartElement(writer(), C_MOLECULE);

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

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

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

    xmlTextWriterEndElement(writer());//molecule

    //TODO: Writing property block

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

//.........这里部分代码省略.........
        else
          d->painter->SetPenColor(d->bondColor);

        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:
          rightAligned = true;
        default:
          break;
      }

      if((d->options & internalColor) && atom->HasData("color"))
        d->painter->SetPenColor(OBColor(atom->GetData("color")->GetValue()));
      else if(d->options & bwAtoms)
        d->painter->SetPenColor(d->bondColor);
      else
        d->painter->SetPenColor(OBColor(etab.GetRGB(atom->GetAtomicNum())));

      //charge and radical
      int charge = atom->GetFormalCharge();
      int spin = atom->GetSpinMultiplicity();
      if(charge || spin) {
        OBFontMetrics metrics = d->painter->GetFontMetrics("N");
        double yoffset = d->HasLabel(atom) ? 0.4 * metrics.height : 0.0;
        switch (GetLabelAlignment(atom)) {
          case TopCenter:
          case TopRight:
          case TopLeft:
          case CenterLeft:
          case CenterRight:
            yoffset = - 1.2 * metrics.height;
        }
        stringstream ss;
        if(charge) {
          if(abs(charge)!=1)
            ss << abs(charge);
          ss << (charge>0 ? "+" : "-") ;
        }
        if(spin) {
          ss << (spin==2 ? "." : "..");
          yoffset += 0.5 * metrics.height;
        }
        if(spin || charge<0)
          d->painter->SetFontSize(2 * metrics.fontSize);
        d->painter->DrawText(x-0.4*metrics.width, y-yoffset, ss.str());
        d->painter->SetFontSize(metrics.fontSize);//restore
      }
 
      if (atom->IsCarbon()) { 
        if(!(d->options & drawAllC))
        {
          if (atom->GetValence() > 1)

//.........这里部分代码省略.........
        _paramFileLoaded = true;
      } else
      {
        return false;
      }
    }

    // Calculate atomic properties based around their ionic charge
    for (i = 0; i < N; i++)
    {
      atom = mol.GetAtom(i + 1);
      a = atom->GetAtomicNum();
      c = _chargeCenter[a];

      // Fail if ionization data is missing for any atom in the molecule
      if (_ionizations[a][c + 1] == -1 || _ionizations[a][c] == -1 || a > TABLE_OF_ELEMENTS_SIZE)
      {
        obErrorLog.ThrowError(__FUNCTION__, "Insufficient ionization data for atoms in the given molecule. Update `data/eqeqIonizations.txt` with missing information and re-run this function.", obError);
        return false;
      }

      J(i) = _ionizations[a][c + 1] - _ionizations[a][c];
      chi(i) = 0.5 * (_ionizations[a][c + 1] + _ionizations[a][c]) - (a == 1? 0 : c * J(i));
    }

    // If a unit cell is defined, use the periodic Ewald calculation
    if (mol.HasData(OBGenericDataType::UnitCell))
    {
      // Get unit cell and calculate its Fourier transform + volume
      obuc = (OBUnitCell *) mol.GetData(OBGenericDataType::UnitCell);
      unitcell = obuc->GetCellMatrix();
      fourier = (2 * PI * unitcell.inverse()).transpose();
      cellVolume = obuc->GetCellVolume();

      // Get the number of radial unit cells to use in x, y, and z
      numNeighbors[0] = int(ceil(minCellLength / (2.0 * (obuc->GetA())))) - 1;
      numNeighbors[1] = int(ceil(minCellLength / (2.0 * (obuc->GetB())))) - 1;
      numNeighbors[2] = int(ceil(minCellLength / (2.0 * (obuc->GetC())))) - 1;

      for (i = 0; i < N; i++)
      {
        atom = mol.GetAtom(i + 1);
        for (j = 0; j < N; j++)
        {
          dx = atom->GetVector() - (mol.GetAtom(j + 1))->GetVector();
          J_ij(i, j) = GetPeriodicEwaldJij(J(i), J(j), dx, (i == j), unitcell, fourier, cellVolume, numNeighbors);
        }
      }
    // If no unit cell, use the simplified nonperiodic calculation
    } else
    {
      for (i = 0; i < N; i++)
      {
        atom = mol.GetAtom(i + 1);
        for (j = 0; j < N; j++)
        {
          J_ij(i, j) = GetNonperiodicJij(J(i), J(j), atom->GetDistance(j + 1), (i == j));
        }
        return false;
      }
    }

    // Formulate problem as A x = b, where x is the calculated partial charges
    // First equation is a simple overall balance: sum(Q) = 0
    A.row(0) = VectorXf::Ones(N);
    b(0) = 0;

    // Remaining equations are based off of the fact that, at equilibrium, the
    // energy of the system changes equally for a change in any charge:
    //     dE/dQ_1 = dE/dQ_2 = ... = dE/dQ_N
    A.block(1, 0, N - 1, N) = J_ij.block(0, 0, N - 1, N) - J_ij.block(1, 0, N - 1, N);
    b.tail(N - 1) = chi.tail(N - 1) - chi.head(N - 1);

    // The solution is a list of charges in the system
    x = A.colPivHouseholderQr().solve(b);

    // Now we are done calculating, pass all this back to OpenBabel molecule
    mol.SetPartialChargesPerceived();
    OBPairData *dp = new OBPairData;
    dp->SetAttribute("PartialCharges");
    dp->SetValue("EQEq");
    dp->SetOrigin(perceived);
    mol.SetData(dp);

    m_partialCharges.clear();
    m_partialCharges.reserve(N);
    m_formalCharges.clear();
    m_formalCharges.reserve(N);

    for (i = 0; i < N; i ++)
    {
      atom = mol.GetAtom(i + 1);
      atom->SetPartialCharge(x(i));
      m_partialCharges.push_back(x(i));
      m_formalCharges.push_back(atom->GetFormalCharge());
    }

    obErrorLog.ThrowError(__FUNCTION__, "EQEq charges successfully assigned.", obInfo);
    return true;
  }

//.........这里部分代码省略.........
                if (strlen(type_name) < 4)
                  {
                    char tmp[16];
                    strncpy(tmp, type_name, 16);
                    snprintf(padded_name, sizeof(padded_name), " %-3s", tmp);
                    strncpy(type_name,padded_name,4);
                    type_name[4] = '\0';
                  }
                else
                  {
                    /*
                      type_name[4] = type_name[3];
                      type_name[3] = type_name[2];
                      type_name[2] = type_name[1];
                      type_name[1] = type_name[0];
                      type_name[0] = type_name[4];
                    */
                    type_name[4] = '\0';
                  }
              }
            res_num = res->GetNum();
          }
        else
          {
            strcpy(the_res,"UNK");
            snprintf(padded_name,sizeof(padded_name), "%s",type_name);
            strncpy(type_name,padded_name,4);
            type_name[4] = '\0';
            res_num = 1;
          }

        element_name = etab.GetSymbol(atom->GetAtomicNum());
        
        int charge = atom->GetFormalCharge();
        char scharge[3] = { ' ', ' ', '\0' };
        if(0 != charge)
          {
            snprintf(scharge, 3, "%+d", charge);
            char tmp = scharge[1];
            scharge[1] = scharge[0];
            scharge[0] = tmp;
          }
        snprintf(buffer, BUFF_SIZE, "%s%5d %-4s %-3s %c%4d    %8.3f%8.3f%8.3f  1.00  0.00          %2s%2s\n",
                 het?"HETATM":"ATOM  ",
                 i,
                 type_name,
                 the_res,
                 the_chain,
                 res_num,
                 atom->GetX(),
                 atom->GetY(),
                 atom->GetZ(),
                 element_name,
                 scharge);
        ofs << buffer;
      }

    OBAtom *nbr;
    vector<OBBond*>::iterator k;
    for (i = 1; i <= mol.NumAtoms(); i ++)
      {
        atom = mol.GetAtom(i);
        if (atom->GetValence() == 0)
          continue; // no need to write a CONECT record -- no bonds

        snprintf(buffer, BUFF_SIZE, "CONECT%5d", i);