C++ Molecule::getEdge方法代码示例

bool MoleculeAutomorphismSearch::_isCisTransBondMappedRigid (Molecule &mol, int i, const int *mapping)
{
   int parity = mol.cis_trans.getParity(i);
   int parity2 = MoleculeCisTrans::applyMapping(parity, mol.cis_trans.getSubstituents(i), mapping, false);

   const Edge &edge = mol.getEdge(i);
   int i2 = mol.findEdgeIndex(mapping[edge.beg], mapping[edge.end]);
   if (mol.cis_trans.getParity(i2) != parity2)
      return false;

   return true;
}

void MoleculeAutomorphismSearch::_findAllPossibleCisTransOneStep (Molecule &mol)
{
   _approximation_orbits_saved.copy(_approximation_orbits);

   // Uniquify such bonds and mark them temorary as cis-trans
   int mark = mol.vertexEnd();
   for (int i = 0; i < possible_cis_trans_to_check.size(); i++)
   {
      int bond = possible_cis_trans_to_check[i];

      int subst[4];
      if (!MoleculeCisTrans::isGeomStereoBond(mol, bond, subst, false))
      {
         possible_cis_trans_to_check.remove(i);
         i--;
         continue;
      }

      if (mol.cis_trans.getParity(bond))
         throw Error("Possible cis-trans check allowed only for non cis-trans bonds");

      mol.cis_trans.add(bond, subst, MoleculeCisTrans::CIS); // Parity doesn't matter

      int validity = _validCisTransBond(bond, _approximation_orbits);
      _cistrans_bond_state[bond] = validity;

      // Uniqufy this bond
      const Edge &e = mol.getEdge(bond);
      _approximation_orbits[e.beg] = mark++;
   }

   _findInvalidStereoCisTrans(mol);

   for (int i = 0; i < possible_cis_trans_to_check.size(); i++)
   {
      int bond = possible_cis_trans_to_check[i];

      int state = _cistrans_bond_state[bond];

      // Restore state
      _cistrans_bond_state[bond] = _NO_STEREO;
      mol.cis_trans.setParity(bond, 0);

      if (state == _INVALID)
      {
         possible_cis_trans_to_check.remove(i);
         i--;
      }
   }

   // Restore orbits
   _approximation_orbits.copy(_approximation_orbits_saved);
}

bool MoleculeAutomorphismSearch::_checkCisTransInvalid (Molecule &mol, int bond_idx)
{
   _target_bond = bond_idx;
   _target_bond_parity_inv = false;
   _fixed_atom = mol.getEdge(bond_idx).beg;

   AutomorphismSearch::process(mol);

   _target_bond = -1;
   _fixed_atom = -1;

   // If _target_bond_parity_inv is true then stereobond is invalid
   return _target_bond_parity_inv;
}

bool edge_intersection (const Molecule &mol, int edge1_idx, int edge2_idx, Vec2f &p)
{
   const Edge &edge1 = mol.getEdge(edge1_idx);
   const Edge &edge2 = mol.getEdge(edge2_idx);

   if (edge1.beg == edge2.beg || edge1.beg == edge2.end || edge1.end == edge2.beg || edge1.end == edge2.end)
      return false;

   Vec2f v1_1(mol.getAtomPos(edge1.beg).x, mol.getAtomPos(edge1.beg).y);
   Vec2f v1_2(mol.getAtomPos(edge1.end).x, mol.getAtomPos(edge1.end).y);
   Vec2f v2_1(mol.getAtomPos(edge2.beg).x, mol.getAtomPos(edge2.beg).y);
   Vec2f v2_2(mol.getAtomPos(edge2.end).x, mol.getAtomPos(edge2.end).y);

   return Vec2f::intersection(v1_1, v1_2, v2_1, v2_2, p);
}

void MoleculePiSystemsMatcher::_calcConnectivity (Molecule &mol, Array<int> &conn)
{
   conn.clear_resize(mol.vertexEnd());
   conn.zerofill();
   for (int e = mol.edgeBegin(); e != mol.edgeEnd(); e = mol.edgeNext(e))
   {
      int bond_order = mol.getBondOrder(e);

      const Edge &edge = mol.getEdge(e);
      conn[edge.beg] += bond_order;
      conn[edge.end] += bond_order;
   }
   for (int v = mol.vertexBegin(); v != mol.vertexEnd(); v = mol.vertexNext(v))
      if (!mol.isPseudoAtom(v) && !mol.isRSite(v))
         conn[v] += mol.getImplicitH(v);
}

int MoleculeInChIUtils::getParityInChI (Molecule &mol, int bond)
{
   if (mol.cis_trans.getParity(bond) == 0)
      throw Error("Specified bond ins't stereogenic");

   const Edge &edge = mol.getEdge(bond);
   
   const int *subst = mol.cis_trans.getSubstituents(bond);
   // Find substituents with maximal indices
   int max_first = __max(subst[0], subst[1]);
   int max_second = __max(subst[2], subst[3]);

   int value = MoleculeCisTrans::sameside(
      mol.getAtomXyz(edge.beg), mol.getAtomXyz(edge.end),
      mol.getAtomXyz(max_first), mol.getAtomXyz(max_second));
   if (value > 0)
      return -1;
   return 1;
}

bool MoleculeAutomorphismSearch::_findInvalidStereoCisTrans (Molecule &mol)
{
   _treat_undef_as = _VALID;

   QS_DEF(Array<int>, invalid_stereo);
   invalid_stereo.clear();

   bool invalid_found = false;

   for (int i = mol.edgeBegin(); i != mol.edgeEnd(); i = mol.edgeNext(i))
   {
      if (_cistrans_bond_state[i] != _UNDEF || mol.cis_trans.getParity(i) == 0)
         continue;

      if (ignored_vertices != 0)
      {
         const Edge &edge = mol.getEdge(i);
         if (ignored_vertices[edge.beg] || ignored_vertices[edge.end])
            continue;
      }

      _cistrans_bond_state[i] = _INVALID;

      if (_checkCisTransInvalid(mol, i))
         // Stereobond is invalid
         invalid_stereo.push(i);

      _cistrans_bond_state[i] = _UNDEF;
   }

   // Mark invalid stereocenters
   for (int i = 0; i < invalid_stereo.size(); i++)
   {
      int bond_index = invalid_stereo[i];
      _cistrans_bond_state[bond_index] = _INVALID;
      invalid_found = true;
   }

   return invalid_found;
}

void MoleculeAutomorphismSearch::_findCisTransStereoBondParirties (Molecule &mol)
{
   const MoleculeStereocenters &stereocenters = mol.stereocenters;

   // Mark edges that connects two stereocenters and that common parity can be detected
   for (int i = mol.edgeBegin(); i != mol.edgeEnd(); i = mol.edgeNext(i))
   {
      const Edge &edge = mol.getEdge(i);
      if (!stereocenters.exists(edge.beg) || !stereocenters.exists(edge.end))
         continue;
      if (stereocenters.getGroup(edge.beg) != stereocenters.getGroup(edge.end))
         continue;
      if (stereocenters.getType(edge.beg) != stereocenters.getType(edge.end))
         continue;

      // Both stereocenters exists and groups are the same

      int orb_beg = _approximation_orbits[edge.beg];
      int orb_end = _approximation_orbits[edge.end];

      _approximation_orbits[edge.beg] = -2; // Some value different from -1 and other orbits
      _approximation_orbits[edge.end] = -2;

      int parity_beg, parity_end;
      if (_validStereocenterByAtom(edge.beg, _approximation_orbits, &parity_beg) == _VALID &&
         _validStereocenterByAtom(edge.end, _approximation_orbits, &parity_end) == _VALID)
      {
         // 1 - means that both stereocenters have the same parity
         _cistrans_stereo_bond_parity[i] = -parity_beg * parity_end;
      }

      // Restore orbits
      _approximation_orbits[edge.beg] = orb_beg;
      _approximation_orbits[edge.end] = orb_end;
   }
}

void CmfSaver::saveMolecule (Molecule &mol)
{
   /* Walk molecule */
   DfsWalk walk(mol);
   QS_DEF(Array<int>, mapping);

   if (_ext_encoder != 0)
      _ext_encoder->start();

   walk.walk();

   /* Get walking sequence */
   const Array<DfsWalk::SeqElem> &v_seq = walk.getSequence();

   /* Calculate mapping to the encoded molecule */
   walk.calcMapping(mapping);
   
   QS_DEF(Array<int>, branch_counters);
   QS_DEF(Array<int>, cycle_numbers);

   branch_counters.clear_resize(mol.vertexEnd());
   branch_counters.zerofill();
   cycle_numbers.clear();

   _atom_sequence.clear();

   QS_DEF(Array<int>, bond_mapping);
   bond_mapping.clear_resize(mol.edgeEnd());
   bond_mapping.fffill();
   int bond_index = 0;

   /* Encode first atom */
   if (v_seq.size() > 0)
   {
      _encodeAtom(mol, v_seq[0].idx, mapping.ptr());
      _atom_sequence.push(v_seq[0].idx);
      
      int j, openings = walk.numOpenings(v_seq[0].idx);

      for (j = 0; j < openings; j++)
      {
         cycle_numbers.push(v_seq[0].idx);
         _encodeCycleNumer(j);
      }
   }

   /* Main cycle */
   int i, j, k;

   for (i = 1; i < v_seq.size(); i++)
   {
      int v_idx = v_seq[i].idx;
      int e_idx = v_seq[i].parent_edge;
      int v_prev_idx = v_seq[i].parent_vertex;
      bool write_atom = true;

      if (v_prev_idx >= 0)
      {
         if (walk.numBranches(v_prev_idx) > 1)
            if (branch_counters[v_prev_idx] > 0)
               _encode(CMF_CLOSE_BRACKET);

         int branches = walk.numBranches(v_prev_idx);

         if (branches > 1)
            if (branch_counters[v_prev_idx] < branches - 1)
               _encode(CMF_OPEN_BRACKET);

         branch_counters[v_prev_idx]++;

         if (branch_counters[v_prev_idx] > branches)
            throw Error("unexpected branch");

         _encodeBond(mol, e_idx, mapping.ptr());
         bond_mapping[e_idx] = bond_index++;

         if (save_bond_dirs)
         {
            int dir = mol.getBondDirection(e_idx);

            if (dir != 0)
            {
               if (dir == BOND_UP)
                  dir = CMF_BOND_UP;
               else if (dir == BOND_DOWN)
                  dir = CMF_BOND_DOWN;
               else
                  dir = CMF_BOND_EITHER;

               const Edge &edge = mol.getEdge(e_idx);

               if (edge.beg == v_prev_idx && edge.end == v_idx)
                  ;
               else if (edge.beg == v_idx && edge.end == v_prev_idx)
                  _encode(CMF_BOND_SWAP_ENDS);
               else
                  throw Error("internal");

               _encode(dir);
            }
//.........这里部分代码省略.........

void IndigoInchi::generateInchiInput (Molecule &mol, inchi_Input &input, 
   Array<inchi_Atom> &atoms, Array<inchi_Stereo0D> &stereo)
{
   QS_DEF(Array<int>, mapping);
   mapping.clear_resize(mol.vertexEnd());
   mapping.fffill();
   int index = 0;
   for (int v = mol.vertexBegin(); v != mol.vertexEnd(); v = mol.vertexNext(v))
      mapping[v] = index++;
   atoms.clear_resize(index);
   atoms.zerofill();

   stereo.clear();
   for (int v = mol.vertexBegin(); v != mol.vertexEnd(); v = mol.vertexNext(v))
   {
      inchi_Atom &atom = atoms[mapping[v]];
      
      int atom_number = mol.getAtomNumber(v);
      if (atom_number == ELEM_PSEUDO)
         throw IndigoError("Molecule with pseudoatom (%s) cannot be converted into InChI", mol.getPseudoAtom(v));
      if (atom_number == ELEM_RSITE)
         throw IndigoError("Molecule with RGroups cannot be converted into InChI");
      strncpy(atom.elname, Element::toString(atom_number), ATOM_EL_LEN);

      Vec3f &c = mol.getAtomXyz(v);
      atom.x = c.x;
      atom.y = c.y;
      atom.z = c.z;
                              
      // connectivity
      const Vertex &vtx = mol.getVertex(v);
      int nei_idx = 0;
      for (int nei = vtx.neiBegin(); nei != vtx.neiEnd(); nei = vtx.neiNext(nei))
      {
         int v_nei = vtx.neiVertex(nei);
         atom.neighbor[nei_idx] = mapping[v_nei];
         int edge_idx = vtx.neiEdge(nei);
         atom.bond_type[nei_idx] = getInchiBondType(mol.getBondOrder(edge_idx));

         int bond_stereo = INCHI_BOND_STEREO_NONE;
         if (mol.cis_trans.isIgnored(edge_idx))
            bond_stereo = INCHI_BOND_STEREO_DOUBLE_EITHER;
         else
         {
            int dir = mol.getBondDirection2(v, v_nei);
            if (mol.getBondDirection2(v, v_nei) == BOND_EITHER)
               bond_stereo = INCHI_BOND_STEREO_SINGLE_1EITHER;
            else if (mol.getBondDirection2(v_nei, v) == BOND_EITHER)
               bond_stereo = INCHI_BOND_STEREO_SINGLE_2EITHER;
         }
         atom.bond_stereo[nei_idx] = bond_stereo;
         nei_idx++;
      }
      atom.num_bonds = vtx.degree();

      // Other properties
      atom.isotopic_mass = mol.getAtomIsotope(v);
      atom.radical = mol.getAtomRadical(v);
      atom.charge = mol.getAtomCharge(v);

      // Hydrogens
      int hcount = -1;
      if (Molecule::shouldWriteHCount(mol, v) || mol.isExplicitValenceSet(v) || mol.isImplicitHSet(v))
      {
         if (mol.getAtomAromaticity(v) == ATOM_AROMATIC &&
            atom_number == ELEM_C && atom.charge == 0 && atom.radical == 0)
         {
            // Do not set number of implicit hydrogens here as InChI throws an exception on
            // the molecule B1=CB=c2cc3B=CC=c3cc12
            ;
         }
         else
            // set -1 to tell InChI add implicit hydrogens automatically
            hcount = mol.getImplicitH_NoThrow(v, -1); 
      }
      atom.num_iso_H[0] = hcount;
   }
  
   // Process cis-trans bonds
   for (int e = mol.edgeBegin(); e != mol.edgeEnd(); e = mol.edgeNext(e))
   {
      if (mol.cis_trans.getParity(e) == 0)
         continue;

      int subst[4];
      mol.cis_trans.getSubstituents_All(e, subst);

      const Edge &edge = mol.getEdge(e);

      inchi_Stereo0D &st = stereo.push();

      // Write it as
      // #0 - #1 = #2 - #3
      st.neighbor[0] = mapping[subst[0]];
      st.neighbor[1] = mapping[edge.beg];
      st.neighbor[2] = mapping[edge.end];
      st.neighbor[3] = mapping[subst[2]];

      if (mol.cis_trans.getParity(e) == MoleculeCisTrans::CIS)
         st.parity = INCHI_PARITY_ODD;
//.........这里部分代码省略.........

void IndigoInchi::parseInchiOutput (const inchi_OutputStruct &inchi_output, Molecule &mol)
{
   mol.clear();

   Array<int> atom_indices;
   atom_indices.clear();
   
   // Add atoms
   for (AT_NUM i = 0; i < inchi_output.num_atoms; i ++)
   {
      const inchi_Atom &inchi_atom = inchi_output.atom[i];

      int idx = mol.addAtom(Element::fromString(inchi_atom.elname));
      atom_indices.push(idx);
   }

   // Add bonds
   for (AT_NUM i = 0; i < inchi_output.num_atoms; i ++)
   {
      const inchi_Atom &inchi_atom = inchi_output.atom[i];
      for (AT_NUM bi = 0; bi < inchi_atom.num_bonds; bi++)
      {
         AT_NUM nei = inchi_atom.neighbor[bi];
         if (i > nei)
            // Add bond only once
            continue;
         int bond_order = inchi_atom.bond_type[bi];
         if (bond_order == INCHI_BOND_TYPE_NONE)
            throw Molecule::Error("Indigo-InChI: NONE-typed bonds are not supported");
         if (bond_order >= INCHI_BOND_TYPE_ALTERN)
            throw Molecule::Error("Indigo-InChI: ALTERN-typed bonds are not supported");
         int bond = mol.addBond(atom_indices[i], atom_indices[nei], bond_order);
      }
   }

   // Add Hydrogen isotope atoms at the end to preserver 
   // the same atom ordering
   for (AT_NUM i = 0; i < inchi_output.num_atoms; i ++)
   {
      const inchi_Atom &inchi_atom = inchi_output.atom[i];

      int root_atom = atom_indices[i];
      for (int iso = 1; iso <= NUM_H_ISOTOPES; iso++)
      {
         int count = inchi_atom.num_iso_H[iso];
         while (count-- > 0)
         {
            int h = mol.addAtom(ELEM_H);
            mol.setAtomIsotope(h, iso);
            mol.addBond(root_atom, h, BOND_SINGLE);
         }
      }
   }

   // Set atom charges, radicals and etc.
   for (int i = 0; i < inchi_output.num_atoms; i++)
   {
      const inchi_Atom &inchi_atom = inchi_output.atom[i];

      int idx = atom_indices[i];
      mol.setAtomCharge(idx, inchi_atom.charge);
      if (inchi_atom.isotopic_mass)
         mol.setAtomIsotope(idx, inchi_atom.isotopic_mass);
      if (inchi_atom.radical)
         mol.setAtomRadical(idx, inchi_atom.radical);
      mol.setImplicitH(idx, inchi_atom.num_iso_H[0]);
   }

   neutralizeV5Nitrogen(mol);

   // Process stereoconfiguration
   for (int i = 0; i < inchi_output.num_stereo0D; i++)
   {
      inchi_Stereo0D &stereo0D = inchi_output.stereo0D[i];
      if (stereo0D.type == INCHI_StereoType_DoubleBond)
      {
         if (stereo0D.parity != INCHI_PARITY_ODD && stereo0D.parity != INCHI_PARITY_EVEN)
            continue;

         int bond = mol.findEdgeIndex(stereo0D.neighbor[1], stereo0D.neighbor[2]);

         bool valid = mol.cis_trans.registerBondAndSubstituents(bond);
         if (!valid)
            throw IndigoError("Indigo-InChI: Unsupported cis-trans configuration for "
               "bond %d (atoms %d-%d-%d-%d)", bond, stereo0D.neighbor[0], stereo0D.neighbor[1], 
               stereo0D.neighbor[2], stereo0D.neighbor[3]);

         int vb, ve;
         const Edge &edge = mol.getEdge(bond);
         if (edge.beg == stereo0D.neighbor[1])
         {
            vb = stereo0D.neighbor[0];
            ve = stereo0D.neighbor[3];
         }
         else if (edge.beg == stereo0D.neighbor[2])
         {
            vb = stereo0D.neighbor[3];
            ve = stereo0D.neighbor[0];
         }
         else
//.........这里部分代码省略.........

TautomerSuperStructure::TautomerSuperStructure (Molecule &mol) :
CP_INIT,
TL_CP_GET(_atomsEmitBond),
TL_CP_GET(_atomsAcceptBond),
TL_CP_GET(_isBondAttachedArray),
TL_CP_GET(_mapping),
TL_CP_GET(_inv_mapping),
TL_CP_GET(_edge_mapping),
TL_CP_GET(_total_h)
{
   int i;

   _inside_ctor = true;

   clone(mol, &_inv_mapping, &_mapping);

   _edge_mapping.clear_resize(edgeEnd());
   _edge_mapping.fffill();

   for (i = mol.edgeBegin(); i != mol.edgeEnd(); i = mol.edgeNext(i))
   {
      const Edge &edge = mol.getEdge(i);

      _edge_mapping[findEdgeIndex(_inv_mapping[edge.beg],
              _inv_mapping[edge.end])] = i;
   }

   // Collect atom properties
   _collectAtomProperties();

   // Detect distances from _atomsEmitBond elements to _atomsAcceptBond elements
   QS_DEF(Array<int>, distancesMatrix);
   distancesMatrix.resize(_atomsEmitBond.size() * _atomsAcceptBond.size());
   for (int i = 0; i < _atomsEmitBond.size(); i++)
   {
      int *result = distancesMatrix.ptr() + _atomsAcceptBond.size() * i;
      _findMinDistance(_atomsEmitBond[i], 6, _atomsAcceptBond, result);
   }


   QS_DEF(Array<int>, attachedBonds);
   attachedBonds.clear();
   for (int i = 0; i < _atomsEmitBond.size(); i++)
      for (int j = 0; j < _atomsAcceptBond.size(); j++)
      {
         int v1 = _atomsEmitBond[i];
         int v2 = _atomsAcceptBond[j];
         if (findEdgeIndex(v1, v2) != -1)
            continue;
         // Check new loop size: 5 or 6
         int size = distancesMatrix[_atomsAcceptBond.size() * i + j];
         if (size != 4 && size != 5)
            continue;

         attachedBonds.push(addEdge(v1, v2));
      }

   _isBondAttachedArray.resize(edgeEnd());
   _isBondAttachedArray.zerofill();
   for (int i = 0; i < attachedBonds.size(); i++)
      _isBondAttachedArray[attachedBonds[i]] = true;

   _inside_ctor = false;
}

TautomerSuperStructure::~TautomerSuperStructure ()
{
}


void TautomerSuperStructure::clear ()
{
   if (_inside_ctor)
      Molecule::clear();
   else
      throw Exception("clear(): not supported");
}

int TautomerSuperStructure::getBondOrder (int idx)
{
   if (!_inside_ctor && _isBondAttachedArray[idx])
      return -1;
   return Molecule::getBondOrder(idx);
}

int TautomerSuperStructure::getBondTopology (int idx)
{
   if (!_inside_ctor &&_isBondAttachedArray[idx])
      return -1;
   return Molecule::getBondTopology(idx);
}

bool TautomerSuperStructure::possibleBondOrder (int idx, int order)
{
   if (!_inside_ctor &&_isBondAttachedArray[idx])
      return order == 0 || order == BOND_SINGLE;

   return Molecule::possibleBondOrder(idx, order);
}

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

void convertMolfile (char *path, char *filename, FileOutput &cpp_file)
{
   FileScanner molfile("%s\\%s", path, filename);
   MolfileLoader mf_loader(molfile);
   Molecule mol;
   QS_DEF(Array<int>, edges);

   printf("%s\n", filename);

   mf_loader.loadMolecule(mol, true);

   BiconnectedDecomposer bd(mol); 

   if (bd.decompose() != 1)
   {
      printf("Error: %s is not biconnected\n", filename);
      return;
   }

   int i, j;

   edges.clear_reserve(mol.edgeCount());

   for (i = mol.edgeBegin() ; i < mol.edgeEnd(); i = mol.edgeNext(i))
      edges.push(i);

   edges.qsort(edge_cmp, &mol);

   const Edge &edge = mol.getEdge(edges[edges.size() / 2]);

   Vec3f v1 = mol.getAtomPos(edge.beg);
   Vec3f v2 = mol.getAtomPos(edge.end);

   v1.z = 0.f;
   v2.z = 0.f;

   float scale = Vec3f::dist(v1, v2);

   if (scale < 0.0001f)
   {
      printf("Error: %s has zero bond\n", filename);
      return;
   }

   scale = 1.f / scale;

   int first_idx = mol.vertexBegin();
   Vec3f pos = mol.getAtomPos(first_idx);

   for (i = mol.vertexNext(first_idx); i < mol.vertexEnd(); i = mol.vertexNext(i))
   {
      if (mol.getAtomPos(i).y < pos.y)
      {
         pos = mol.getAtomPos(i);
         first_idx = i;
      }
   }


   for (i = mol.vertexBegin() ; i < mol.vertexEnd(); i = mol.vertexNext(i))
   {
      mol.getAtom2(i).pos.sub(pos);
      mol.getAtom2(i).pos.scale(scale);
   }

   char buf[1024];

   sprintf_s(buf, "BEGIN_PATTERN(\"%s\")", filename);
   cpp_file.writeStringCR(buf);   

   for (i = mol.vertexBegin(); i < mol.vertexEnd(); i = mol.vertexNext(i))
   {
      sprintf_s(buf, "   ADD_ATOM(%d, %ff, %ff)", i, mol.getAtomPos(i).x, mol.getAtomPos(i).y);
      cpp_file.writeStringCR(buf);   
   }

   for (i = mol.edgeBegin(); i < mol.edgeEnd(); i = mol.edgeNext(i))
   {
      const Edge &edge = mol.getEdge(i);
      int type = mol.getBond(i).type;
      int qtype = mol.getQueryBond(i).type;

      sprintf_s(buf, "   ADD_BOND(%d, %d, %d)", edge.beg, edge.end, qtype != 0 ? qtype : type);
      cpp_file.writeStringCR(buf);   
   }

   Vec2f v, inter;
   Vec2f pos_i;
   int idx = mol.vertexCount();

   i = first_idx;

   float max_angle, cur_angle;
   float i_angle = 0;
   int next_nei = 0;
   int point_idx = 0;

   pos_i.set(mol.getAtomPos(i).x, mol.getAtomPos(i).y);

   while (true)
//.........这里部分代码省略.........

bool LayeredMolecules::addLayerFromMolecule(const Molecule &molecule, Array<int> &aam)
{
   Array<int> aam_inverse;
   aam_inverse.expandFill(aam.size(), -1);
   for(int i = 0; i < aam.size(); ++i)
   {
      if(aam[i] != -1)
         aam_inverse[aam[i]] = i;
   }

   unsigned newTautomerIndex = layers;
   _resizeLayers(newTautomerIndex + 1);
   for(auto e1_idx : edges())
   {
      _bond_masks[BOND_ZERO][e1_idx].reset(newTautomerIndex);
      _bond_masks[BOND_SINGLE][e1_idx].reset(newTautomerIndex);
      _bond_masks[BOND_DOUBLE][e1_idx].reset(newTautomerIndex);
      _bond_masks[BOND_TRIPLE][e1_idx].reset(newTautomerIndex);
      _bond_masks[BOND_AROMATIC][e1_idx].reset(newTautomerIndex);
   }

   unsigned node = _trie.getRoot();
   bool unique = false;

   for(auto e2_idx : const_cast<Molecule&>(molecule).edges())
   {
      auto e2 = molecule.getEdge(e2_idx);
      int u2 = e2.beg;
      int v2 = e2.end;
      int u1 = aam_inverse[u2];
      int v1 = aam_inverse[v2];
      if(u1 == -1 || v1 == -1)
         continue;
      int e1_idx = findEdgeIndex(u1, v1);
      if(e1_idx == -1)
      {
         e1_idx = addEdge(u1, v1);
         _proto.addEdge(u1, v1);
         _proto.setBondOrder(e1_idx, BOND_ZERO, false);
         _bond_masks[BOND_ZERO].resize(e1_idx + 1);
         _bond_masks[BOND_SINGLE].resize(e1_idx + 1);
         _bond_masks[BOND_DOUBLE].resize(e1_idx + 1);
         _bond_masks[BOND_TRIPLE].resize(e1_idx + 1);
         _bond_masks[BOND_AROMATIC].resize(e1_idx + 1);
      }
      int order = const_cast<Molecule&>(molecule).getBondOrder(e2_idx);
      _bond_masks[order][e1_idx].set(newTautomerIndex);
   }

   for (auto e1_idx : edges())
   {
      int order = BOND_ZERO;
      if(_bond_masks[BOND_SINGLE][e1_idx].get(newTautomerIndex))
         order = BOND_SINGLE;
      else if(_bond_masks[BOND_DOUBLE][e1_idx].get(newTautomerIndex))
         order = BOND_DOUBLE;
      else if(_bond_masks[BOND_TRIPLE][e1_idx].get(newTautomerIndex))
         order = BOND_TRIPLE;

      bool newlyAdded;
      node = _trie.add(node, order, newlyAdded);
      unique = (newlyAdded ? true : unique);
   }

   if(unique)
   {
      ++layers;
      return true;
   }

   // This means that we avoided adding non-unique layer, and we need to reduce the size of bitsets.
   _resizeLayers(layers);
   return false;
}