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

  bool OBDepict::AddAtomLabels(AtomLabelType type)
  {
    d->painter->SetPenColor(OBColor("red"));
    d->painter->SetFillColor(OBColor("red"));
    d->painter->SetFontSize((int)(GetFontSize() * 0.8));// smaller text
    OBAtomIterator i;
    for (OBAtom *atom = d->mol->BeginAtom(i); atom; atom = d->mol->NextAtom(i)) {
      vector3 pos(atom->GetVector());
      std::stringstream ss;
      switch (type) {
        case AtomId:
          ss << atom->GetId();
          d->painter->DrawText(pos.x(), pos.y(), ss.str());
          break;
        case AtomSymmetryClass:
          ss << GetAtomSymClass(atom);
          d->painter->DrawText(pos.x(), pos.y(), ss.str());
          break;
        case AtomIndex:
          ss << atom->GetIdx();
          d->painter->DrawText(pos.x(), pos.y(), ss.str());
          break;

        default:
          break;
      }
    }

    return true;    
  }

  int GetLabelAlignment(OBAtom *atom) 
  {
    // compute the sum of the bond vectors, this gives 
    vector3 direction(VZero);
    OBBondIterator i;
    for (OBAtom *nbr = atom->BeginNbrAtom(i); nbr; nbr = atom->NextNbrAtom(i))
      direction += atom->GetVector() - nbr->GetVector();
    
    const double bias = -0.1; //towards left-alignment, which is more natural
    int alignment = 0;
    if (direction.y() < 0.0) {
      if (direction.x() < bias)
        alignment = BottomRight;
      else 
        alignment = BottomLeft;
    } else if (direction.x() > 0.0) {
      if (direction.x() < bias)
        alignment = TopRight;
      else
        alignment = TopLeft;
    } else {
      if (direction.x() < bias)
         alignment = CenterRight;
      else 
        alignment = CenterLeft;
    }

    return alignment;
  }

bool BoxFormat::WriteMolecule(OBBase* pOb, OBConversion* pConv)
{
    OBMol* pmol = dynamic_cast<OBMol*>(pOb);
    if(pmol==NULL)
        return false;

    //Define some references so we can use the old parameter names
    ostream &ofs = *pConv->GetOutStream();
    OBMol &mol = *pmol;

    //margin hardwired in new framework. Also was in old fileformat
    double margin=1.0;

    char buffer[BUFF_SIZE];
    vector3 vcenter,vmin,vmax,vmid,vdim;

    OBAtom *atom;
    vector<OBAtom*>::iterator i;
    vmax.Set(-10E10,-10E10,-10E10);
    vmin.Set( 10E10, 10E10, 10E10);

    for (atom = mol.BeginAtom(i); atom; atom = mol.NextAtom(i))
    {
        vcenter += atom->GetVector();
        if (atom->x() < vmin.x())
            vmin.SetX(atom->x());
        if (atom->y() < vmin.y())
            vmin.SetY(atom->y());
        if (atom->z() < vmin.z())
            vmin.SetZ(atom->z());

        if (atom->x() > vmax.x())
            vmax.SetX(atom->x());
        if (atom->y() > vmax.y())
            vmax.SetY(atom->y());
        if (atom->z() > vmax.z())
            vmax.SetZ(atom->z());
    }
    vcenter /= (double)mol.NumAtoms();

    vector3 vmarg(margin,margin,margin);
    vmin -= vmarg;
    vmax += vmarg;
    vdim = vmax - vmin;
    vmid = vmin+vmax;
    vmid /= 2.0;

    ofs << "HEADER    CORNERS OF BOX" << endl;
    snprintf(buffer, BUFF_SIZE, "REMARK    CENTER (X Y Z)      %10.3f %10.3f %10.3f",
             vmid.x(),vmid.y(),vmid.z());
    ofs << buffer << endl;
    snprintf(buffer, BUFF_SIZE, "REMARK    DIMENSIONS (X Y Z)  %10.3f %10.3f %10.3f",
             vdim.x(),vdim.y(),vdim.z());
    ofs << buffer << endl;
    vdim /= 2.0;

    vector3 vtmp;
    int j;
    for (j = 1; j <= 8; j++)
    {
        switch(j)
        {
        case 1:
            vtmp = vmid-vdim;
            break;
        case 2:
            vtmp.SetX(vmid.x()+vdim.x());
            break;
        case 3:
            vtmp.SetZ(vmid.z()+vdim.z());
            break;
        case 4:
            vtmp.SetX(vmid.x()-vdim.x());
            break;
        case 5:
            vtmp = vmid-vdim;
            vtmp.SetY(vmid.y()+vdim.y());
            break;
        case 6:
            vtmp = vmid+vdim;
            vtmp.SetZ(vmid.z()-vdim.z());
            break;
        case 7:
            vtmp = vmid+vdim;
            break;
        case 8:
            vtmp.SetX(vmid.x()-vdim.x());
            break;
        }
        snprintf(buffer, BUFF_SIZE, "ATOM      %d  DUA BOX     1    %8.3f%8.3f%8.3f",
                 j,vtmp.x(),vtmp.y(),vtmp.z());
        ofs << buffer << endl;
    }

    ofs << "CONECT    1    2    4    5" << endl;
    ofs << "CONECT    2    1    3    6" << endl;
    ofs << "CONECT    3    2    4    7" << endl;
    ofs << "CONECT    4    1    3    8" << endl;
    ofs << "CONECT    5    1    6    8" << endl;
    ofs << "CONECT    6    2    5    7" << endl;
//.........这里部分代码省略.........

  //! \return whether partial charges were successfully assigned to this molecule
  bool EQEqCharges::ComputeCharges(OBMol &mol)
  {
    int i, j, a, c, N = mol.NumAtoms();
    double cellVolume;
    VectorXf chi(N), J(N), b(N), x(N);
    MatrixXf J_ij(N, N), A(N, N);
    OBUnitCell *obuc;
    matrix3x3 unitcell, fourier;
    vector3 dx;
    int numNeighbors[3];
    OBAtom *atom;

    // If parameters have not yet been loaded, do that
    if (!_paramFileLoaded)
    {
      if (ParseParamFile())
      {
        _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);
//.........这里部分代码省略.........

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