Searches S-S bonds based on distances
            between atoms in the structure (first model only).
            Average distance is 2.05A. Threshold is 3A default.
            Returns iterator with tuples of residues.
517777
def search_ss_bonds(self, threshold=3.0):
        
        
        # Taken from http://docs.python.org/library/itertools.html
        # Python 2.4 does not include itertools.combinations

        def combinations(iterable, r):
            # combinations('ABCD', 2) --> AB AC AD BC BD CD
            # combinations(range(4), 3) --> 012 013 023 123
            pool = tuple(iterable)
            n = len(pool)
            if r > n:
                return
            indices = range(r)
            yield tuple(pool[i] for i in indices)
            while True:
                for i in reversed(range(r)):
                    if indices[i] != i + n - r:
                        break
                else:
                    return
                indices[i] += 1
                for j in range(i+1, r):
                    indices[j] = indices[j-1] + 1
                yield tuple(pool[i] for i in indices)

        
        
        model = self.child_list[0]
        cysteines = [r for r in model.get_residues() if r.get_resname() == 'CYS']

        pairs = combinations(cysteines, 2) # Iterator with pairs

        for cys_pair in pairs:
            if cys_pair[0]['SG'] - cys_pair[1]['SG'] < threshold:
                yield cys_pair
887616
def get_molecule(self, index=0):
        
        return Molecule(self.numbers, self.geometries[index], self.titles[index], symbols=self.symbols)
452598
def bonds(lines, atoms):
    
    # Convert sdf style stereobond (see chem.model.bond.Bond)
    conv_stereo_table = {0: 0, 1: 1, 3: 3, 4: 3, 6: 2}
    results = {a: {} for a in atoms}
    for line in lines:
        bond = Bond()
        first = int(line[0:3])
        second = int(line[3:6])
        if first > second:
            bond.is_lower_first = 0
        order = int(line[6:9])
        if order < 4:
            bond.order = order
        bond.type = conv_stereo_table[int(line[9:12])]
        results[first][second] = {"bond": bond}
        results[second][first] = {"bond": bond}
    return results
37133
private int getAtomTypeXCount(IAtomContainer ac, IAtom atom) {
        List<IAtom> neighbours = ac.getConnectedAtomsList(atom);
        int nocounter = 0;
        IBond bond;
        for (IAtom neighbour : neighbours) {
            if ((neighbour.getSymbol().equals("N") || neighbour.getSymbol().equals("O"))
                    && !(Boolean) neighbour.getProperty("IS_IN_AROMATIC_RING")) {
                //if (ac.getMaximumBondOrder(neighbours[i]) == 1.0) {
                bond = ac.getBond(neighbour, atom);
                if (bond.getOrder() != IBond.Order.DOUBLE) {
                    nocounter += 1;
                }
            }
        }
        return nocounter;
    }
189291
def from_structure(structure):
    

    try:
        from Bio import PDB
        if isinstance(structure, str):
            p = PDB.PDBParser()
            structure = p.get_structure('S', structure)
        if isinstance(structure, PDB.Structure.Structure):
            for _ in structure.get_chains():
                atoms = [np.array(atom.get_coord())
                         for atom in structure.get_atoms()]
    except ImportError:
        print("Biopython not found.")
        raise

    atoms = np.array(structure)
    try:
        import scipy
        D = scipy.spatial.distance.pdist(atoms, 'euclidean')
        D = scipy.spatial.distance.squareform(D)
    except ImportError:
        print("Scipy not found.")
        raise
    m = np.max(1 / D[D != 0])
    M = np.zeros(D.shape)
    M[D != 0] = 1 / D[D != 0]
    M[D == 0] = m
    return M
325924
def getCentroid(attribute_variants, comparator):
    

    n = len(attribute_variants)

    distance_matrix = numpy.zeros([n, n])

    # populate distance matrix by looping through elements of matrix triangle
    for i in range(0, n):
        for j in range(0, i):
            distance = comparator(attribute_variants[i], attribute_variants[j])
            distance_matrix[i, j] = distance_matrix[j, i] = distance

    average_distance = distance_matrix.mean(0)

    # there can be ties for minimum, average distance string
    min_dist_indices = numpy.where(
        average_distance == average_distance.min())[0]

    if len(min_dist_indices) > 1:
        centroid = breakCentroidTie(attribute_variants, min_dist_indices)
    else:
        centroid_index = min_dist_indices[0]
        centroid = attribute_variants[centroid_index]

    return centroid
344895
def get_structure(atoms, cls=None):
        
        symbols = atoms.get_chemical_symbols()
        positions = atoms.get_positions()
        lattice = atoms.get_cell()

        cls = Structure if cls is None else cls
        return cls(lattice, symbols, positions,
                   coords_are_cartesian=True)
406492
private IBond getNextBondWithUnplacedRingAtom() {
        for (IBond bond : molecule.bonds()) {
            IAtom beg = bond.getBegin();
            IAtom end = bond.getEnd();
            if (beg.getPoint2d() != null && end.getPoint2d() != null) {
                if (end.getFlag(CDKConstants.ISPLACED) && !beg.getFlag(CDKConstants.ISPLACED) && beg.isInRing()) {
                    return bond;
                }
                if (beg.getFlag(CDKConstants.ISPLACED) && !end.getFlag(CDKConstants.ISPLACED) && end.isInRing()) {
                    return bond;
                }
            }
        }
        return null;
    }
819269
def calc2dcoords(mol):
    
    topology.recognize(mol)
    g = set(i for i, _ in mol.atoms_iter())
    # 1: get nodes in scaffolds
    scaffolds = []
    belongs = {}
    for i, rkeys in enumerate(sorted(mol.scaffolds, key=len)):
        scf = []
        for rkey in rkeys:
            ring = mol.rings[rkey]
            for r in ring:
                belongs[r] = i
            scf.append(ring)
            g -= set(ring)
        scaffolds.append(scf)
    # 2: traverse nodes and scaffolds
    # the node and scaffold graph should be a tree (no cycles)
    f = True
    # print(scaffolds)
    coords = {}
    while g:
        if f and scaffolds:  # largest scaffold is first
            stack = [scaffolds[-1][0][0]]
            f = False
        else:
            stack = [g.pop()]
        pred = {}
        branch = {}
        while stack:
            # print("stack: {}".format(stack))
            tail = stack.pop()
            # print("tail: {}".format(tail))
            if tail in belongs:  # scaffolds
                scf = scaffold_coords(scaffolds[belongs[tail]])
                # print(scf.keys())
                # rotate and translate
                if not coords:
                    coords = scf
                else:
                    u = coords[pred[tail]]
                    v = scf[tail]
                    op = [u[0] + math.cos(u[2]), u[1] + math.sin(u[2])]
                    translate(scf, gm.vector(v[:2], op))
                    rotate(scf, op, gm.rad(u[2] + math.pi - v[2]))
                    coords.update(scf)
                # stack nbrs of scaffold
                for k in scf.keys():
                    pred[k] = None
                    for nbr in mol.neighbors(k):
                        if nbr not in scf.keys():
                            stack.append(nbr)
                            pred[nbr] = k
            else:  # append linker
                if tail not in pred:  # isolated
                    coords[tail] = [0, 0, 0, 1]
                    continue
                p = pred[tail]
                x, y, ang, d = coords[p]
                # TODO: ring configuration
                coords[tail] = [x + math.cos(ang), y + math.sin(ang),
                                ang + d * math.pi / 3, d * -1]
                if p not in branch:
                    coords[p][2] = gm.rad(coords[p][2] + math.pi * 2 / 3 * d)
                    branch[p] = 1
                elif branch[p] == 1:
                    coords[p][2] = gm.rad(coords[p][2] + math.pi * d)
                    branch[p] += 1
                elif branch[p] == 2:
                    coords[p][2] = gm.rad(coords[p][2] + math.pi * 2 / 3 * d)
                    branch[p] += 1
                for nbr in mol.neighbors(tail):
                    if nbr not in pred:
                        stack.append(nbr)
                        pred[nbr] = tail
        g -= set(pred)
    resolve_overlap(coords)
    for i, a in mol.atoms_iter():
        mol.atom(i).coords = coords[i][:2]
654835
public static int getVertexCountAtDistance(IAtomContainer atomContainer, int distance) {
        int natom = atomContainer.getAtomCount();

        int[][] admat = AdjacencyMatrix.getMatrix(atomContainer);
        int[][] distanceMatrix = computeFloydAPSP(admat);

        int matches = 0;

        for (int i = 0; i < natom; i++) {
            for (int j = 0; j < natom; j++) {
                if (distanceMatrix[i][j] == distance) matches++;
            }
        }
        return matches / 2;
    }
1149754
public IAtom getNextUnplacedHeavyAtomWithAliphaticPlacedNeighbour(IAtomContainer molecule) {
        Iterator<IBond> bonds = molecule.bonds().iterator();
        while (bonds.hasNext()) {
            IBond bond = bonds.next();
            if (bond.getBegin().getFlag(CDKConstants.ISPLACED) && !(bond.getEnd().getFlag(CDKConstants.ISPLACED))) {
                if (isAliphaticHeavyAtom(bond.getEnd())) {
                    return bond.getEnd();
                }
            }
            if (bond.getEnd().getFlag(CDKConstants.ISPLACED) && !(bond.getBegin().getFlag(CDKConstants.ISPLACED))) {
                if (isAliphaticHeavyAtom(bond.getBegin())) {
                    return bond.getBegin();
                }
            }
        }
        return null;
    }
1237288
public IRingSet findSSSR(IAtomContainer mol) {
        IBond brokenBond = null;
        IChemObjectBuilder builder = mol.getBuilder();
        IRingSet sssr = builder.newInstance(IRingSet.class);
        IAtomContainer molecule = builder.newInstance(IAtomContainer.class);
        molecule.add(mol);
        IAtom smallest;
        int smallestDegree, nodesToBreakCounter, degree;
        IAtom[] rememberNodes;
        IRing ring;

        //Two Vectors - as defined in the article. One to hold the
        //full set of atoms in the structure and on to store the numbers
        //of the nodes that have been trimmed away.
        //Furhter there is a Vector nodesN2 to store the number of N2 nodes
        List<IAtom> fullSet = new ArrayList<IAtom>();
        List<IAtom> trimSet = new ArrayList<IAtom>();
        List<IAtom> nodesN2 = new ArrayList<IAtom>();

        initPath(molecule);
        logger.debug("molecule.getAtomCount(): " + molecule.getAtomCount());
        // load fullSet with the numbers of our atoms
        for (int f = 0; f < molecule.getAtomCount(); f++) {
            fullSet.add(molecule.getAtom(f));
        }
        logger.debug("fullSet.size(): " + fullSet.size());

        do {
            //Add nodes of degree zero to trimset.
            //Also add nodes of degree 2 to nodesN2.
            //In the same run, check, which node has the lowest degree
            //greater than zero.
            smallestDegree = 7;
            smallest = null;
            nodesN2.clear();
            for (int f = 0; f < molecule.getAtomCount(); f++) {
                IAtom atom = molecule.getAtom(f);
                degree = molecule.getConnectedBondsCount(atom);
                if (degree == 0) {
                    if (!trimSet.contains(atom)) {
                        logger.debug("Atom of degree 0");
                        trimSet.add(atom);
                    }
                }
                if (degree == 2) {
                    nodesN2.add(atom);
                }
                if (degree < smallestDegree && degree > 0) {
                    smallest = atom;
                    smallestDegree = degree;
                }
            }
            if (smallest == null) break;

            // If there are nodes of degree 1, trim them away
            if (smallestDegree == 1) {
                trimCounter++;
                trim(smallest, molecule);
                trimSet.add(smallest);
            }

            // if there are nodes of degree 2, find out of which rings
            // they are part of.
            else if (smallestDegree == 2) {
                rememberNodes = new IAtom[nodesN2.size()];
                nodesToBreakCounter = 0;
                for (int f = 0; f < nodesN2.size(); f++) {
                    ring = getRing((IAtom) nodesN2.get(f), molecule);
                    if (ring != null) {
                        // check, if this ring already is in SSSR
                        if (!RingSetManipulator.ringAlreadyInSet(ring, sssr)) {
                            sssr.addAtomContainer(ring);
                            rememberNodes[nodesToBreakCounter] = (IAtom) nodesN2.get(f);
                            nodesToBreakCounter++;
                        }
                    }
                }
                if (nodesToBreakCounter == 0) {
                    nodesToBreakCounter = 1;
                    rememberNodes[0] = (IAtom) nodesN2.get(0);
                }
                for (int f = 0; f < nodesToBreakCounter; f++) {
                    breakBond(rememberNodes[f], molecule);
                }
                if (brokenBond != null) {
                    molecule.addBond(brokenBond);
                    brokenBond = null;
                }
            }
            // if there are nodes of degree 3
            else if (smallestDegree == 3) {
                ring = getRing(smallest, molecule);
                if (ring != null) {

                    // check, if this ring already is in SSSR
                    if (!RingSetManipulator.ringAlreadyInSet(ring, sssr)) {
                        sssr.addAtomContainer(ring);
                    }
                    brokenBond = checkEdges(ring, molecule);
                    molecule.removeElectronContainer(brokenBond);
                }
            }
        } while (trimSet.size() < fullSet.size());
        logger.debug("fullSet.size(): " + fullSet.size());
        logger.debug("trimSet.size(): " + trimSet.size());
        logger.debug("trimCounter: " + trimCounter);
        //		molecule.setProperty(CDKConstants.SMALLEST_RINGS, sssr);
        return sssr;
    }
1269110
public static Set<IBond> findMappedBonds(IReaction reaction) {
        Set<IBond> mapped = new HashSet<>();

        // first we collect the occurrance of mapped bonds from reacants then products
        Set<IntTuple> mappedReactantBonds = new HashSet<>();
        Set<IntTuple> mappedProductBonds  = new HashSet<>();
        for (IAtomContainer reactant : reaction.getReactants().atomContainers()) {
            for (IBond bond : reactant.bonds()) {
                Integer begidx = bond.getBegin().getProperty(CDKConstants.ATOM_ATOM_MAPPING);
                Integer endidx = bond.getEnd().getProperty(CDKConstants.ATOM_ATOM_MAPPING);
                if (begidx != null && endidx != null)
                    mappedReactantBonds.add(new IntTuple(begidx, endidx));
            }
        }
        // fail fast
        if (mappedReactantBonds.isEmpty())
            return Collections.emptySet();

        for (IAtomContainer product : reaction.getProducts().atomContainers()) {
            for (IBond bond : product.bonds()) {
                Integer begidx = bond.getBegin().getProperty(CDKConstants.ATOM_ATOM_MAPPING);
                Integer endidx = bond.getEnd().getProperty(CDKConstants.ATOM_ATOM_MAPPING);
                if (begidx != null && endidx != null)
                    mappedProductBonds.add(new IntTuple(begidx, endidx));
            }
        }
        // fail fast
        if (mappedProductBonds.isEmpty())
            return Collections.emptySet();

        // repeat above but now store any that are different or unmapped as being mapped
        for (IAtomContainer reactant : reaction.getReactants().atomContainers()) {
            for (IBond bond : reactant.bonds()) {
                Integer begidx = bond.getBegin().getProperty(CDKConstants.ATOM_ATOM_MAPPING);
                Integer endidx = bond.getEnd().getProperty(CDKConstants.ATOM_ATOM_MAPPING);
                if (begidx != null && endidx != null && mappedProductBonds.contains(new IntTuple(begidx, endidx)))
                    mapped.add(bond);
            }
        }
        for (IAtomContainer product : reaction.getProducts().atomContainers()) {
            for (IBond bond : product.bonds()) {
                Integer begidx = bond.getBegin().getProperty(CDKConstants.ATOM_ATOM_MAPPING);
                Integer endidx = bond.getEnd().getProperty(CDKConstants.ATOM_ATOM_MAPPING);
                if (begidx != null && endidx != null && mappedReactantBonds.contains(new IntTuple(begidx, endidx)))
                    mapped.add(bond);
            }
        }
        return mapped;
    }
184392
def find_hal(self, atoms):
        """"""
        data = namedtuple('hal_acceptor', 'o o_orig_idx y y_orig_idx')
        a_set = []
        # All oxygens, nitrogen, sulfurs with neighboring carbon, phosphor, nitrogen or sulfur
        for a in [at for at in atoms if at.atomicnum in [8, 7, 16]]:
            n_atoms = [na for na in pybel.ob.OBAtomAtomIter(a.OBAtom) if na.GetAtomicNum() in [6, 7, 15, 16]]
            if len(n_atoms) == 1:  # Proximal atom
                o_orig_idx = self.Mapper.mapid(a.idx, mtype=self.mtype, bsid=self.bsid)
                y_orig_idx = self.Mapper.mapid(n_atoms[0].GetIdx(), mtype=self.mtype, bsid=self.bsid)
                a_set.append(data(o=a, o_orig_idx=o_orig_idx, y=pybel.Atom(n_atoms[0]), y_orig_idx=y_orig_idx))
        return a_set
1311138
private static boolean queryAdjacency(IBond bondA1, IBond bondB1, IBond bondA2, IBond bondB2) {

        IAtom atom1 = null;
        IAtom atom2 = null;

        if (bondA1.contains(bondB1.getBegin())) {
            atom1 = bondB1.getBegin();
        } else if (bondA1.contains(bondB1.getEnd())) {
            atom1 = bondB1.getEnd();
        }

        if (bondA2.contains(bondB2.getBegin())) {
            atom2 = bondB2.getBegin();
        } else if (bondA2.contains(bondB2.getEnd())) {
            atom2 = bondB2.getEnd();
        }

        if (atom1 != null && atom2 != null) {
            // well, this looks fishy: the atom2 is not always atom IQueryAtom !
            return ((IQueryAtom) atom2).matches(atom1);
        } else {
            return atom1 == null && atom2 == null;
        }

    }
819792
def centers_list(self):
        
        center = set()
        adj = defaultdict(set)
        for n, atom in self.atoms():
            if atom._reactant != atom._product:
                center.add(n)

        for n, m, bond in self.bonds():
            if bond._reactant != bond._product:
                adj[n].add(m)
                adj[m].add(n)
                center.add(n)
                center.add(m)

        out = []
        while center:
            n = center.pop()
            if n in adj:
                c = set(self.__plain_bfs(adj, n))
                out.append(list(c))
                center.difference_update(c)
            else:
                out.append([n])

        return out
9881
def get_surface_subseq_3D(protein,
                              depth_prop='RES_DEPTH-msms', depth_cutoff=2.5, depth_condition='<',
                              acc_prop='RSA_ALL-freesasa_het', acc_cutoff=25, acc_condition='>'):
        """"""
        empty = {'surface_3D': {'subseq_len'    : 0, 'subseq': None,
                                   'subseq_resnums': []},
                    'notdeep_3D': {'subseq_len'    : 0, 'subseq': None,
                                   'subseq_resnums': []},
                    'acc_3D'    : {'subseq_len'    : 0, 'subseq': None,
                                   'subseq_resnums': []}}
        if not protein.representative_structure:
            log.error('{}: no representative structure, cannot search for subseq'.format(protein.id))
            return empty

        notdeep_subseq, notdeep_subseq_resnums = protein.get_seqprop_subsequence_from_structchain_property(
            property_key=depth_prop,
            property_value=depth_cutoff,
            condition=depth_condition,
            use_representatives=True,
            return_resnums=True) or (None, [])

        acc_subseq, acc_subseq_resnums = protein.get_seqprop_subsequence_from_structchain_property(
            property_key=acc_prop,
            property_value=acc_cutoff,
            condition=acc_condition,
            use_representatives=True,
            return_resnums=True) or (None, [])

        surface_subseq_resnums = list(set(notdeep_subseq_resnums).intersection(acc_subseq_resnums))
        surface_subseq = protein.representative_sequence.get_subsequence(surface_subseq_resnums)

        all_info = {'surface_3D': {'subseq_len'    : len(surface_subseq_resnums), 'subseq': surface_subseq,
                                   'subseq_resnums': surface_subseq_resnums},
                    'notdeep_3D': {'subseq_len'    : len(notdeep_subseq_resnums), 'subseq': notdeep_subseq,
                                   'subseq_resnums': notdeep_subseq_resnums},
                    'acc_3D'    : {'subseq_len'    : len(acc_subseq_resnums), 'subseq': acc_subseq,
                                   'subseq_resnums': acc_subseq_resnums}}
        return all_info