determine connectivity structure. for simplicity only connects atoms of the same residue, peptide bonds and nucleotides
9835
function() {
    console.time('Mol.deriveConnectivity');
    var thisStructure = this;
    var prevResidue = null;
    this.eachResidue(function(res) {
      var sqrDist;
      var atoms = res.atoms();
      var numAtoms = atoms.length;
      for (var i = 0; i < numAtoms; i+=1) {
        var atomI = atoms[i];
        var posI = atomI.pos();
        var covalentI = covalentRadius(atomI.element());
        for (var j = 0; j < i; j+=1) {
          var atomJ = atoms[j];
          var covalentJ = covalentRadius(atomJ.element());
          sqrDist = vec3.sqrDist(posI, atomJ.pos());
          var lower = covalentI+covalentJ-0.30;
          var upper = covalentI+covalentJ+0.30;
          if (sqrDist < upper*upper && sqrDist > lower*lower) {
            thisStructure.connect(atomI, atomJ);
          }
        }
      }
      res._deduceType();
      if (prevResidue !== null) {
        if (res.isAminoacid() && prevResidue.isAminoacid()) {
          connectPeptides(thisStructure, prevResidue, res);
        }
        if (res.isNucleotide() && prevResidue.isNucleotide()) {
          connectNucleotides(thisStructure, prevResidue, res);
        }
      }
      prevResidue = res;
    });
    console.timeEnd('Mol.deriveConnectivity');
  }
626144
def environment(self, atom):
        
        return tuple((bond, self._node[n]) for n, bond in self._adj[atom].items())
447386
private int getPiSystemsCount(IAtomContainer ac, IAtom atom) {
        List neighbours = ac.getConnectedAtomsList(atom);
        int picounter = 0;
        List bonds = null;
        for (int i = 0; i < neighbours.size(); i++) {
            IAtom neighbour = (IAtom) neighbours.get(i);
            bonds = ac.getConnectedBondsList(neighbour);
            for (int j = 0; j < bonds.size(); j++) {
                IBond bond = (IBond) bonds.get(j);
                if (bond.getOrder() != IBond.Order.SINGLE && !bond.getOther(neighbour).equals(atom)
                        && !neighbour.getSymbol().equals("P") && !neighbour.getSymbol().equals("S")) {
                    picounter += 1;
                }/*
                  * else if (bonds[j].getOther(neighbours[i])!=atom &&
                  * !neighbours[i].getSymbol().equals("P") &&
                  * !neighbours[i].getSymbol().equals("S") &&
                  * bonds[j].getOther
                  * (neighbours[i]).getFlag(CDKConstants.ISAROMATIC)){ picounter
                  * += 1; }
                  */
            }
        }
        return picounter;
    }
364155
def cylinder_and_strand(self):
        '''

        '''
        top = self.topology
        # We build a  mini-state machine to find the
        # start end of helices and such
        in_helix = False
        helices_starts = []
        helices_ends = []
        coils = []

        coil = []
        for i, typ in enumerate(top['secondary_structure']):
            if typ == 'H':
                if in_helix == False:
                    # We become helices
                    helices_starts.append(top['residue_indices'][i][0])
                    in_helix = True

                    # We end the previous coil
                    coil.append(top['residue_indices'][i][0])
            else:
                if in_helix == True:
                    # We stop being helices
                    helices_ends.append(top['residue_indices'][i][0])

                    # We start a new coil
                    coil = []
                    coils.append(coil)
                    in_helix = False

                # We add control points
                coil.append(top['residue_indices'][i][0])
                [coil.append(j) for j in top['residue_indices'][i] if top['atom_names'][j] == 'CA']

        # We add the coils
        coil_representations = []
        for control_points in coils:
            rid = self.add_representation('smoothtube', {'coordinates': self.coordinates[control_points],
                                                   'radius': 0.05,
                                                   'resolution': 4,
                                                   'color': 0xffffff})
            coil_representations.append(rid)


        start_idx, end_idx = helices_starts, helices_ends
        cylinders = self.add_representation('cylinders', {'startCoords': self.coordinates[list(start_idx)],
                                              'endCoords': self.coordinates[list(end_idx)],
                                              'colors': [0xffff00] * len(self.coordinates),
                                              'radii': [0.15] * len(self.coordinates)})
        def update(self=self, cylinders=cylinders, coils=coils,
                   coil_representations=coil_representations,
                   start_idx=start_idx, end_idx=end_idx, control_points=control_points):
            for i, control_points in enumerate(coils):
                rid = self.update_representation(coil_representations[i],
                                                 {'coordinates': self.coordinates[control_points]})

            self.update_representation(cylinders, {'startCoords': self.coordinates[list(start_idx)],
                                                  'endCoords': self.coordinates[list(end_idx)]})

        self.update_callbacks.append(update)
        self.autozoom(self.coordinates)
404842
private static int checkAtom(IAtomContainer ac, IAtom currentAtom) {
        int check = -1;
        List<IAtom> atoms = ac.getConnectedAtomsList(currentAtom);
        List<IBond> bonds = ac.getConnectedBondsList(currentAtom);
        if (currentAtom.getFlag(CDKConstants.ISAROMATIC)) {
            check = 0;
        } else if (currentAtom.getFormalCharge() == 1 /*
                                                       * &&
                                                       * currentAtom.getSymbol
                                                       * ().equals("C")
                                                       */) {
            check = 0;
        } else if (currentAtom.getFormalCharge() == -1) {
            //// NEGATIVE CHARGES WITH A NEIGHBOOR PI BOND //////////////
            int counterOfPi = 0;
            for (IAtom atom : atoms) {
                if (ac.getMaximumBondOrder(atom) != IBond.Order.SINGLE) {
                    counterOfPi++;
                }
            }
            if (counterOfPi > 0) check = 0;
        } else {
            int se = ac.getConnectedSingleElectronsCount(currentAtom);
            if (se == 1) {
                check = 0; //// DETECTION of radicals
            } else if (ac.getConnectedLonePairsCount(currentAtom) > 0
            /* && (currentAtom.getSymbol().equals("N") */) {
                check = 0; //// DETECTION of  lone pair
            } else {
                int highOrderBondCount = 0;
                for (int j = 0; j < atoms.size(); j++) {
                    IBond bond = bonds.get(j);
                    if (bond == null || bond.getOrder() != IBond.Order.SINGLE) {
                        highOrderBondCount++;
                    } else {
                    }
                }
                if (highOrderBondCount == 1) {
                    check = 0;
                } else if (highOrderBondCount > 1) {
                    check = 1;
                }
            }
        }
        return check;
    }
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)
1064025
private void assembleLinkages(List<List<Atom[]>> matchedAtomsOfLinkages,
			ProteinModification mod, List<ModifiedCompound> ret) {
		ModificationCondition condition = mod.getCondition();
		List<ModificationLinkage> modLinks = condition.getLinkages();

		int nLink = matchedAtomsOfLinkages.size();
		int[] indices = new int[nLink];
		Set<ModifiedCompound> identifiedCompounds = new HashSet<ModifiedCompound>();
		while (indices[0]<matchedAtomsOfLinkages.get(0).size()) {
			List<Atom[]> atomLinkages = new ArrayList<Atom[]>(nLink);
			for (int iLink=0; iLink<nLink; iLink++) {
				Atom[] atoms = matchedAtomsOfLinkages.get(iLink).get(indices[iLink]);
				atomLinkages.add(atoms);
			}
			if (matchLinkages(modLinks, atomLinkages)) {
				// matched

				int n = atomLinkages.size();
				List<StructureAtomLinkage> linkages = new ArrayList<StructureAtomLinkage>(n);
				for (int i=0; i<n; i++) {
					Atom[] linkage = atomLinkages.get(i);
					StructureAtomLinkage link = StructureUtil.getStructureAtomLinkage(
							linkage[0], residues.contains(linkage[0].getGroup()),
							linkage[1], residues.contains(linkage[1].getGroup()));
					linkages.add(link);
				}

				ModifiedCompound mc = new ModifiedCompoundImpl(mod, linkages);
				if (!identifiedCompounds.contains(mc)) {
					ret.add(mc);
					identifiedCompounds.add(mc);
				}
			}

			// indices++ (e.g. [0,0,1]=>[0,0,2]=>[1,2,0])
			int i = nLink-1;
			while (i>=0) {
				if (i==0 || indices[i]<matchedAtomsOfLinkages.get(i).size()-1) {
					indices[i]++;
					break;
				} else {
					indices[i] = 0;
					i--;
				}
			}
		}
	}
198374
def side_chain(self):
        
        side_chain_atoms = []
        if self.mol_code != 'GLY':
            covalent_bond_graph = generate_covalent_bond_graph(
                find_covalent_bonds(self))
            try:
                subgraphs = generate_bond_subgraphs_from_break(
                    covalent_bond_graph, self['CA'], self['CB'])
                if len(subgraphs) == 1:
                    subgraphs = generate_bond_subgraphs_from_break(
                        subgraphs[0], self['CD'], self['N'])
                if len(subgraphs) == 2:
                    for g in subgraphs:
                        if self['CB'] in g:
                            side_chain_atoms = g.nodes()
                            break
            except:
                warning_message = "Malformed PDB for Residue {0}: {1}.".format(
                    self.id, self)
                if 'CB' in self.atoms.keys():
                    side_chain_atoms.append(self['CB'])
                    warning_message += " Side-chain is just the CB atom."
                else:
                    warning_message += " Empty side-chain."
                warnings.warn(warning_message, MalformedPDBWarning)
        return side_chain_atoms
1231111
def select_connected_bonds():
    ''''''
    s = current_system()
    start, end = s.bonds.transpose()
    selected = np.zeros(s.n_bonds, 'bool')
    for i in selected_atoms():
        selected |= (i == start) | (i == end)

    csel = current_selection()
    bsel = csel['bonds'].add(
        Selection(selected.nonzero()[0], s.n_bonds))
    
    ret = csel.copy()
    ret['bonds'] = bsel

    return select_selection(ret)
293210
def select_molecules(name):
    ''''''
    mol_formula = current_system().get_derived_molecule_array('formula')
    mask = mol_formula == name
    ind = current_system().mol_to_atom_indices(mask.nonzero()[0])

    selection = {'atoms': Selection(ind, current_system().n_atoms)}

    # Need to find the bonds between the atoms
    b = current_system().bonds
    if len(b) == 0:
        selection['bonds'] = Selection([], 0)
    else:
        molbonds = np.zeros(len(b), 'bool')
        for i in ind:
            matching = (i == b).sum(axis=1).astype('bool')
            molbonds[matching] = True
        
        selection['bonds'] = Selection(molbonds.nonzero()[0], len(b))

    return select_selection(selection)
933648
int[][] findNeighborIndices() {

		// looking at a typical protein case, number of neighbours are from ~10 to ~50, with an average of ~30
		int initialCapacity = 60;

		int[][] nbsIndices = new int[atomCoords.length][];

		for (int k=0; k<atomCoords.length; k++) {
			double radius = radii[k] + probe + probe;

			List<Integer> thisNbIndices = new ArrayList<>(initialCapacity);

			for (int i = 0; i < atomCoords.length; i++) {
				if (i == k) continue;

				double dist = atomCoords[i].distance(atomCoords[k]);

				if (dist < radius + radii[i]) {
					thisNbIndices.add(i);
				}
			}

			int[] indicesArray = new int[thisNbIndices.size()];
			for (int i=0;i<thisNbIndices.size();i++) indicesArray[i] = thisNbIndices.get(i);
			nbsIndices[k] = indicesArray;
		}
		return nbsIndices;
	}
1371333
private static BitSet available(int[][] graph, IAtom[] atoms, EdgeToBondMap bonds) {

        final BitSet available = new BitSet();

        // for all atoms, select those that require a double-bond
        ATOMS: for (int i = 0; i < atoms.length; i++) {

            final IAtom atom = atoms[i];

            // preconditions
            if (atom.getAtomicNumber() == null)
                throw new IllegalArgumentException("atom " + (i + 1) + " had unset atomic number");
            if (atom.getFormalCharge() == null)
                throw new IllegalArgumentException("atom " + (i + 1) + " had unset formal charge");
            if (atom.getImplicitHydrogenCount() == null)
                throw new IllegalArgumentException("atom " + (i + 1) + " had unset implicit hydrogen count");

            if (!atom.getFlag(ISAROMATIC)) continue;

            // count preexisting pi-bonds, a higher bond order causes a skip
            int nPiBonds = 0;
            for (final int w : graph[i]) {
                IBond.Order order = bonds.get(i, w).getOrder();
                if (order == DOUBLE) {
                    nPiBonds++;
                } else if (order.numeric() > 2) {
                    continue ATOMS;
                }
            }

            // check if a pi bond can be assigned
            final int element = atom.getAtomicNumber();
            final int charge = atom.getFormalCharge();
            final int valence = graph[i].length + atom.getImplicitHydrogenCount() + nPiBonds;

            if (available(element, charge, valence)) {
                available.set(i);
            }
        }

        return available;
    }
1073179
def getligs(self):
        

        if config.PEPTIDES == [] and config.INTRA is None:
            # Extract small molecule ligands (default)
            ligands = []

            # Filter for ligands using lists
            ligand_residues, self.lignames_all, self.water = self.filter_for_ligands()

            all_res_dict = {(a.GetName(), a.GetChain(), a.GetNum()): a for a in ligand_residues}
            self.lignames_kept = list(set([a.GetName() for a in ligand_residues]))

            if not config.BREAKCOMPOSITE:
                #  Update register of covalent links with those between DNA/RNA subunits
                self.covalent += nucleotide_linkage(all_res_dict)
                #  Find fragment linked by covalent bonds
                res_kmers = self.identify_kmers(all_res_dict)
            else:
                res_kmers = [[a, ] for a in ligand_residues]
            write_message("{} ligand kmer(s) detected for closer inspection.\n".format(len(res_kmers)), mtype='debug')
            for kmer in res_kmers:  # iterate over all ligands and extract molecules + information
                if len(kmer) > config.MAX_COMPOSITE_LENGTH:
                    write_message("Ligand kmer(s) filtered out with a length of {} fragments ({} allowed).\n".format(
                        len(kmer), config.MAX_COMPOSITE_LENGTH), mtype='debug')
                else:
                    ligands.append(self.extract_ligand(kmer))

        else:
            # Extract peptides from given chains
            self.water = [o for o in pybel.ob.OBResidueIter(self.proteincomplex.OBMol) if o.GetResidueProperty(9)]
            if config.PEPTIDES != []:
                peptide_ligands = [self.getpeptides(chain) for chain in config.PEPTIDES]
            elif config.INTRA is not None:
                peptide_ligands = [self.getpeptides(config.INTRA), ]

            ligands = [p for p in peptide_ligands if p is not None]
            self.covalent, self.lignames_kept, self.lignames_all = [], [], set()

        return [lig for lig in ligands if len(lig.mol.atoms) != 0]
376120
def substructure(mol, query, largest_only=True, ignore_hydrogen=True):
    
    def subset_filter(cnt1, cnt2):
        diff = cnt2
        diff.subtract(cnt1)
        if any(v < 0 for v in diff.values()):
            return True

    if not (len(mol) and len(query)):
        return False  # two blank molecules are not isomorphic
    m = molutil.clone(mol)
    q = molutil.clone(query)
    if largest_only:
        m = molutil.largest_graph(m)
        q = molutil.largest_graph(q)
    if ignore_hydrogen:
        m = molutil.make_Hs_implicit(m)
        q = molutil.make_Hs_implicit(q)
    if filter_(m, q, f=subset_filter):
        gm = GraphMatcher(q.graph, m.graph, node_match=atom_match)
        return gm.subgraph_is_isomorphic()
    return False
522223
def get_connected_components_as_subgraphs(graph):
    
    components = get_connected_components(graph)

    list_of_graphs = []

    for c in components:
        edge_ids = set()
        nodes = [graph.get_node(node) for node in c]
        for n in nodes:
            # --Loop through the edges in each node, to determine if it should be included
            for e in n['edges']:
                # --Only add the edge to the subgraph if both ends are in the subgraph
                edge = graph.get_edge(e)
                a, b = edge['vertices']
                if a in c and b in c:
                    edge_ids.add(e)
        # --Build the subgraph and add it to the list
        list_of_edges = list(edge_ids)
        subgraph = make_subgraph(graph, c, list_of_edges)
        list_of_graphs.append(subgraph)

    return list_of_graphs
810488
private static IAtom findPiBondedOxygen(IAtomContainer container, IAtom atom) {
        for (IBond bond : container.getConnectedBondsList(atom)) {
            if (bond.getOrder() == IBond.Order.DOUBLE) {
                IAtom neighbor = bond.getOther(atom);
                int charge = neighbor.getFormalCharge() == null ? 0 : neighbor.getFormalCharge();
                if (neighbor.getAtomicNumber() == 8 && charge == 0) return neighbor;
            }
        }
        return null;
    }
985556
static public SimpleGraph getMoleculeGraph(IAtomContainer molecule) {
        SimpleGraph graph = new SimpleGraph();
        for (int i = 0; i < molecule.getAtomCount(); i++) {
            IAtom atom = molecule.getAtom(i);
            graph.addVertex(atom);
        }

        for (int i = 0; i < molecule.getBondCount(); i++) {
            IBond bond = molecule.getBond(i);

            /*
             * int order = (int) bond.getOrder(); for (int j=0; j<order; j++) {
             * graph.addEdge(bond.getAtoms()[0], bond.getAtoms()[1]); }
             */
            graph.addEdge(bond.getBegin(), bond.getEnd());
        }
        return graph;
    }