Extract the ligand by copying atoms and bonds and assign all information necessary for later steps.
515832
def extract_ligand(self, kmer):
        """"""
        data = namedtuple('ligand', 'mol hetid chain position water members longname type atomorder can_to_pdb')
        members = [(res.GetName(), res.GetChain(), int32_to_negative(res.GetNum())) for res in kmer]
        members = sort_members_by_importance(members)
        rname, rchain, rnum = members[0]
        write_message("Finalizing extraction for ligand %s:%s:%s with %i elements\n" %
                      (rname, rchain, rnum, len(kmer)), mtype='debug')
        names = [x[0] for x in members]
        longname = '-'.join([x[0] for x in members])

        if config.PEPTIDES != []:
            ligtype = 'PEPTIDE'
        elif config.INTRA is not None:
            ligtype = 'INTRA'
        else:
            # Classify a ligand by its HETID(s)
            ligtype = classify_by_name(names)
        write_message("Ligand classified as {}\n".format(ligtype), mtype='debug')

        hetatoms = set()
        for obresidue in kmer:
            hetatoms_res = set([(obatom.GetIdx(), obatom) for obatom in pybel.ob.OBResidueAtomIter(obresidue)
                                if obatom.GetAtomicNum() != 1])

            if not config.ALTLOC:
                # Remove alternative conformations (standard -> True)
                hetatoms_res = set([atm for atm in hetatoms_res
                                    if not self.mapper.mapid(atm[0], mtype='protein',
                                                             to='internal') in self.altconformations])
            hetatoms.update(hetatoms_res)
        write_message("Hetero atoms determined (n={})\n".format(len(hetatoms)), mtype='debug')

        hetatoms = dict(hetatoms)  # make it a dict with idx as key and OBAtom as value
        lig = pybel.ob.OBMol()  # new ligand mol
        neighbours = dict()
        for obatom in hetatoms.values():  # iterate over atom objects
            idx = obatom.GetIdx()
            lig.AddAtom(obatom)
            # ids of all neighbours of obatom
            neighbours[idx] = set([neighbour_atom.GetIdx() for neighbour_atom
                                   in pybel.ob.OBAtomAtomIter(obatom)]) & set(hetatoms.keys())
        write_message("Atom neighbours mapped\n", mtype='debug')

        ##############################################################
        # map the old atom idx of OBMol to the new idx of the ligand #
        ##############################################################

        newidx = dict(zip(hetatoms.keys(), [obatom.GetIdx() for obatom in pybel.ob.OBMolAtomIter(lig)]))
        mapold = dict(zip(newidx.values(), newidx))
        # copy the bonds
        for obatom in hetatoms:
            for neighbour_atom in neighbours[obatom]:
                bond = hetatoms[obatom].GetBond(hetatoms[neighbour_atom])
                lig.AddBond(newidx[obatom], newidx[neighbour_atom], bond.GetBondOrder())
        lig = pybel.Molecule(lig)

        # For kmers, the representative ids are chosen (first residue of kmer)
        lig.data.update({'Name': rname, 'Chain': rchain, 'ResNr': rnum})

        # Check if a negative residue number is represented as a 32 bit integer
        if rnum > 10 ** 5:
            rnum = int32_to_negative(rnum)

        lig.title = ':'.join((rname, rchain, str(rnum)))
        self.mapper.ligandmaps[lig.title] = mapold

        write_message("Renumerated molecule generated\n", mtype='debug')

        if not config.NOPDBCANMAP:
            atomorder = canonicalize(lig)
        else:
            atomorder = None

        can_to_pdb = {}
        if atomorder is not None:
            can_to_pdb = {atomorder[key-1]: mapold[key] for key in mapold}

        ligand = data(mol=lig, hetid=rname, chain=rchain, position=rnum, water=self.water,
                      members=members, longname=longname, type=ligtype, atomorder=atomorder,
                      can_to_pdb=can_to_pdb)
        return ligand
19759
def copy(self):
        """"""
        rv = object.__new__(self.__class__)
        rv.__dict__.update(self.__dict__)
        rv.symbols = self.symbols.copy()
        return rv
1282083
def detect_aromatic_rings_in_ligand(self):
        """"""
        self.ligrings = {}
        try:
            ring_info = self.topology_data.mol.GetRingInfo()
            self.ligand_ring_num = ring_info.NumRings()
        except Exception as e:
            m = Chem.MolFromPDBFile("lig.pdb")
            ring_info = m.GetRingInfo()
            self.ligand_ring_num = ring_info.NumRings()
        i=0
        for ring in range(self.ligand_ring_num):
            if 4 < len(ring_info.AtomRings()[ring]) <= 6 and  False not in [self.topology_data.mol.GetAtomWithIdx(x).GetIsAromatic() for x in ring_info.AtomRings()[ring]]: #narrow ring definition
                atom_ids_in_ring = []
                for atom in ring_info.AtomRings()[ring]:
                    atom_ids_in_ring.append(self.topology_data.universe.ligand.atoms[atom].name)
                self.ligrings[i]=atom_ids_in_ring
                i+=1
427823
def clone_subgraphs(self, g):
        if not isinstance(g, CGRContainer):
            raise InvalidData('only CGRContainer acceptable')

        r_group = []
        x_group = {}
        r_group_clones = []
        newcomponents = []

        ''' 
        '''
        components, lost_bonds, term_atoms = self.__split_graph(g)
        lost_map = {x: y for x, y in lost_bonds}
        ''' extract subgraphs and sort by group type (R or X)
        '''
        x_terminals = set(lost_map.values())
        r_terminals = set(lost_map)

        for i in components:
            x_terminal_atom = x_terminals.intersection(i)
            if x_terminal_atom:
                x_group[x_terminal_atom.pop()] = i
                continue

            r_terminal_atom = r_terminals.intersection(i)
            if r_terminal_atom:
                r_group.append([r_terminal_atom, i])
                continue

            newcomponents.append(i)
        ''' search similar R groups and patch.
        '''
        tmp = g
        for i in newcomponents:
            for k, j in r_group:
                gm = GraphMatcher(j, i, node_match=self.__node_match_products,
                                  edge_match=self.__edge_match_products)
                ''' search for similar R-groups started from bond breaks.
                '''
                mapping = next((x for x in gm.subgraph_isomorphisms_iter() if k.issubset(x) and
                                all(x[y] in term_atoms for y in k)), None)
                if mapping:
                    r_group_clones.append([k, mapping])
                    tmp = compose(tmp, self.__remap_group(j, tmp, mapping)[0])
                    break

        ''' add lose X groups to R groups
        '''
        for i, j in r_group_clones:
            for k in i:
                remappedgroup, mapping = self.__remap_group(x_group[lost_map[k]], tmp, {})
                tmp = CGRcore.union(tmp, remappedgroup)
                tmp.add_edge(j[k], mapping[lost_map[k]], s_bond=1, sp_bond=(1, None))

        if r_group_clones:
            tmp.meta.update(g.meta)
            return tmp

        return tmp.copy()
1220013
def copy(self):
        """"""
        result = Vector3(self.size, self.deriv)
        result.x.v = self.x.v
        result.y.v = self.y.v
        result.z.v = self.z.v
        if self.deriv > 0:
            result.x.d[:] = self.x.d
            result.y.d[:] = self.y.d
            result.z.d[:] = self.z.d
        if self.deriv > 1:
            result.x.dd[:] = self.x.dd
            result.y.dd[:] = self.y.dd
            result.z.dd[:] = self.z.dd
        return result
104383
private void addDoubleBondStereochemistry(Graph g, IAtomContainer ac) {

        for (final Edge e : g.edges()) {

            if (e.bond() != Bond.DOUBLE) continue;

            int u = e.either();
            int v = e.other(u);

            // find a directional bond for either end
            Edge first  = null;
            Edge second = null;

            // if either atom is not incident to a directional label there
            // is no configuration
            if ((first = findDirectionalEdge(g, u)) != null) {
                if ((second = findDirectionalEdge(g, v)) != null) {

                    // if the directions (relative to the double bond) are the
                    // same then they are on the same side - otherwise they
                    // are opposite
                    Conformation conformation = first.bond(u) == second.bond(v) ?
                                                Conformation.TOGETHER : Conformation.OPPOSITE;

                    // get the stereo bond and build up the ligands for the
                    // stereo-element - linear search could be improved with
                    // map or API change to double bond element
                    IBond db = ac.getBond(ac.getAtom(u), ac.getAtom(v));

                    IBond[] ligands = new IBond[]{
                        ac.getBond(ac.getAtom(u), ac.getAtom(first.other(u))),
                        ac.getBond(ac.getAtom(v), ac.getAtom(second.other(v)))};

                    ac.addStereoElement(new DoubleBondStereochemistry(db, ligands, conformation));
                } else if (g.degree(v) == 2) {
                    List<Edge> edges = new ArrayList<>();
                    edges.add(e);
                    Edge f = findCumulatedEdge(g, v, e);
                    while (f != null) {
                        edges.add(f);
                        v = f.other(v);
                        f = findCumulatedEdge(g, v, f);
                    }
                    // only odd number of cumulated bonds here, otherwise is
                    // extended tetrahedral
                    if ((edges.size() & 0x1) == 0)
                        continue;
                    second = findDirectionalEdge(g, v);
                    if (second != null) {
                        int   cfg        = first.bond(u) == second.bond(v)
                                           ? IStereoElement.TOGETHER
                                           : IStereoElement.OPPOSITE;
                        Edge  middleEdge = edges.get(edges.size()/2);
                        IBond middleBond = ac.getBond(ac.getAtom(middleEdge.either()),
                                                      ac.getAtom(middleEdge.other(middleEdge.either())));
                        IBond[] ligands = new IBond[]{
                            ac.getBond(ac.getAtom(u), ac.getAtom(first.other(u))),
                            ac.getBond(ac.getAtom(v), ac.getAtom(second.other(v)))};

                        ac.addStereoElement(new ExtendedCisTrans(middleBond, ligands, cfg));
                    }
                }
            }
            // extension F[C@]=[C@@]F
            else {
                Configuration uConf = g.configurationOf(u);
                Configuration vConf = g.configurationOf(v);
                if (uConf.type() == Configuration.Type.DoubleBond &&
                    vConf.type() == Configuration.Type.DoubleBond) {

                    int[] nbrs = new int[6];
                    int[] uNbrs = g.neighbors(u);
                    int[] vNbrs = g.neighbors(v);

                    if (uNbrs.length < 2 || uNbrs.length > 3)
                        continue;
                    if (vNbrs.length < 2 || vNbrs.length > 3)
                        continue;

                    int idx   = 0;
                    System.arraycopy(uNbrs, 0, nbrs, idx, uNbrs.length);
                    idx += uNbrs.length;
                    if (uNbrs.length == 2) nbrs[idx++] = u;
                    System.arraycopy(vNbrs, 0, nbrs, idx, vNbrs.length);
                    idx += vNbrs.length;
                    if (vNbrs.length == 2) nbrs[idx] = v;
                    Arrays.sort(nbrs, 0, 3);
                    Arrays.sort(nbrs, 3, 6);

                    int vPos = Arrays.binarySearch(nbrs, 0, 3, v);
                    int uPos = Arrays.binarySearch(nbrs, 3, 6, u);

                    int uhi = 0, ulo = 0;
                    int vhi = 0, vlo = 0;

                    uhi = nbrs[(vPos + 1) % 3];
                    ulo = nbrs[(vPos + 2) % 3];
                    vhi = nbrs[3 + ((uPos + 1) % 3)];
                    vlo = nbrs[3 + ((uPos + 2) % 3)];

                    if (uConf.shorthand() == Configuration.CLOCKWISE) {
                        int tmp = uhi;
                        uhi = ulo;
                        ulo = tmp;
                    }
                    if (vConf.shorthand() == Configuration.ANTI_CLOCKWISE) {
                        int tmp = vhi;
                        vhi = vlo;
                        vlo = tmp;
                    }

                    DoubleBondStereochemistry.Conformation conf = null;
                    IBond[] bonds = new IBond[2];

                    if (uhi != u) {
                        bonds[0] = ac.getBond(ac.getAtom(u), ac.getAtom(uhi));
                        if (vhi != v) {
                            // System.err.println(uhi + "\\=/" + vhi);
                            conf = Conformation.TOGETHER;
                            bonds[1] = ac.getBond(ac.getAtom(v), ac.getAtom(vhi));
                        } else if (vlo != v) {
                            // System.err.println(uhi + "\\=\\" + vlo);
                            conf = Conformation.OPPOSITE;
                            bonds[1] = ac.getBond(ac.getAtom(v), ac.getAtom(vlo));
                        }
                    } else if (ulo != u) {
                        bonds[0] = ac.getBond(ac.getAtom(u), ac.getAtom(ulo));
                        if (vhi != v) {
                            // System.err.println(ulo + "/=/" + vhi);
                            conf = Conformation.OPPOSITE;
                            bonds[1] = ac.getBond(ac.getAtom(v), ac.getAtom(vhi));
                        } else if (vlo != v) {
                            // System.err.println(ulo + "/=\\" + vlo);
                            conf = Conformation.TOGETHER;
                            bonds[1] = ac.getBond(ac.getAtom(v), ac.getAtom(vlo));
                        }
                    }

                    ac.addStereoElement(new DoubleBondStereochemistry(ac.getBond(ac.getAtom(u), ac.getAtom(v)),
                                                                      bonds,
                                                                      conf));
                }
            }
        }
    }
615432
def copy_from(self, other):
        
        # Need to copy all attributes, fields, relations
        self.__attributes__ = {k: v.copy() for k, v in other.__attributes__.items()}
        self.__fields__ = {k: v.copy() for k, v in other.__fields__.items()}
        self.__relations__ = {k: v.copy() for k, v in other.__relations__.items()}
        self.maps = {k: m.copy() for k, m in other.maps.items()}
        self.dimensions = other.dimensions.copy()
99892
def mol_from_text(text, assign_descriptors=True):
    
    cg = mols_from_text(text, False, assign_descriptors)
    return next(cg)
162548
def from_filepath(filepath, strict = True, parse_ligands = False):
        ''''''
        return PDB(read_file(filepath), strict = strict, parse_ligands = parse_ligands)
889608
def copy(self, newdata=None):
        ''''''
        if newdata is None:
            newdata = self.data.copy()
        return self.__class__(
            self.molecule, self.origin.copy(), self.axes.copy(),
            self.nrep.copy(), newdata, self.subtitle, self.nuclear_charges
        )
1014186
def to_rdkit_molecule(data):
    
    mol = RWMol()
    conf = Conformer()
    mapping = {}
    is_3d = False
    for n, a in data.atoms():
        ra = Atom(a.number)
        ra.SetAtomMapNum(n)
        if a.charge:
            ra.SetFormalCharge(a.charge)
        if a.isotope != a.common_isotope:
            ra.SetIsotope(a.isotope)
        if a.radical:
            ra.SetNumRadicalElectrons(a.radical)
        mapping[n] = m = mol.AddAtom(ra)
        conf.SetAtomPosition(m, (a.x, a.y, a.z))
        if a.z:
            is_3d = True
    if not is_3d:
        conf.Set3D(False)

    for n, m, b in data.bonds():
        mol.AddBond(mapping[n], mapping[m], _bond_map[b.order])

    mol.AddConformer(conf)
    SanitizeMol(mol)
    return mol
640456
def assign_valence(mol):
    """"""
    for u, v, bond in mol.bonds_iter():
        if bond.order == 2:
            mol.atom(u).pi = 1
            mol.atom(v).pi = 1
            if mol.atom(u).symbol == "O" and not mol.atom(u).charge:
                mol.atom(v).carbonyl_C = 1
            if mol.atom(v).symbol == "O" and not mol.atom(v).charge:
                mol.atom(u).carbonyl_C = 1
        elif bond.order == 3:
            mol.atom(u).pi = mol.atom(v).pi = 2
    max_nbr = {"C": 4, "Si": 4, "N": 3, "P": 3, "As": 3,
               "O": 2, "S": 2, "Se": 2, "F": 1, "Cl": 1, "Br": 1, "I": 1}
    for i, nbrs in mol.neighbors_iter():
        atom = mol.atom(i)
        if len(nbrs) == 2 and all(bond.order == 2 for bond in nbrs.values()):
            atom.pi = 2  # sp (allene, ketene)
        if atom.symbol in max_nbr:
            h_cnt = max_nbr[atom.symbol] - len(nbrs) - atom.pi + atom.charge
            if h_cnt > 0:
                mol.atom(i).add_hydrogen(h_cnt)
    mol.descriptors.add("Valence")
1100766
def formula(mol):
    
    mol.require("Valence")
    mol.require("Topology")
    total_cntr = Counter()
    for m in sorted(mols_iter(mol), key=len, reverse=True):
        cntr = Counter()
        for i in m:
            cntr += mol.atom(i).composition()
        text = []
        Cs = cntr.pop("C", 0)
        if Cs:
            text.append("C")
            if Cs > 1:
                text.append(str(Cs))
        Hs = cntr.pop("H", 0)
        if Hs:
            text.append("H")
            if Hs > 1:
                text.append(str(Hs))
        heteros = sorted(cntr.items(), key=lambda x: atom_number(x[0]))
        for k, v in heteros:
            text.append(k)
            if v > 1:
                text.append(str(v))
        total_cntr["".join(text)] += 1
    total = sorted(total_cntr.items(), key=lambda x: len(x[0]), reverse=True)
    total_text = []
    for k, v in total:
        if v > 1:
            total_text.append(str(v) + k)
        else:
            total_text.append(k)
    return ".".join(total_text)
1222381
def ready_to_draw(mol):
    
    copied = molutil.clone(mol)
    # display_terminal_carbon(mol)
    equalize_terminal_double_bond(copied)
    # spine_to_terminal_wedge(copied)
    scale_and_center(copied)
    format_ring_double_bond(copied)
    return copied
1175133
def from_rdkit_molecule(data):
    
    m = MoleculeContainer()
    atoms, mapping = [], []
    for a in data.GetAtoms():
        atom = {'element': a.GetSymbol(), 'charge': a.GetFormalCharge()}
        atoms.append(atom)
        mapping.append(a.GetAtomMapNum())

        isotope = a.GetIsotope()
        if isotope:
            atom['isotope'] = isotope
        radical = a.GetNumRadicalElectrons()
        if radical:
            atom['multiplicity'] = radical + 1

    conformers = data.GetConformers()
    if conformers:
        for atom, (x, y, z) in zip(atoms, conformers[0].GetPositions()):
            atom['x'] = x
            atom['y'] = y
            atom['z'] = z

    for atom, mapping in zip(atoms, mapping):
        a = m.add_atom(atom)
        if mapping:
            m.atom(a)._parsed_mapping = mapping

    for bond in data.GetBonds():
        m.add_bond(bond.GetBeginAtomIdx() + 1, bond.GetEndAtomIdx() + 1, _rdkit_bond_map[bond.GetBondType()])

    return m
709214
def _get_bonds(self, mol):
        
        num_atoms = len(mol)
        # index starting from 0
        if self.ignore_ionic_bond:
            covalent_atoms = [i for i in range(num_atoms) if mol.species[i].symbol not in self.ionic_element_list]
        else:
            covalent_atoms = list(range(num_atoms))
        all_pairs = list(itertools.combinations(covalent_atoms, 2))
        pair_dists = [mol.get_distance(*p) for p in all_pairs]
        elements = mol.composition.as_dict().keys()
        unavailable_elements = list(set(elements) -
                                    set(self.covalent_radius.keys()))
        if len(unavailable_elements) > 0:
            raise ValueError("The covalent radius for element {} is not "
                             "available".format(unavailable_elements))
        bond_13 = self.get_13_bonds(self.priority_bonds)
        max_length = [(self.covalent_radius[mol.sites[p[0]].specie.symbol] +
                       self.covalent_radius[mol.sites[p[1]].specie.symbol]) *
                      (1 + (self.priority_cap if p in self.priority_bonds
                            else (self.bond_length_cap if p not in bond_13
                                  else self.bond_13_cap))) *
                      (0.1 if (self.ignore_halogen_self_bond and p not in self.priority_bonds and
                               mol.sites[p[0]].specie.symbol in self.halogen_list and
                               mol.sites[p[1]].specie.symbol in self.halogen_list)
                       else 1.0)
                      for p in all_pairs]

        bonds = [bond
                 for bond, dist, cap in zip(all_pairs, pair_dists, max_length)
                 if dist <= cap]
        return bonds
1229480
def from_lines(pdb_file_lines, strict = True, parse_ligands = False):
        ''''''
        return PDB("\n".join(pdb_file_lines), strict = strict, parse_ligands = parse_ligands)