A donor atom already forming a hydrogen bond is not allowed to form a water bridge. Each water molecule
        can only be donor for two water bridges, selecting the constellation with the omega angle closest to 110 deg.
430188
def refine_water_bridges(self, wbridges, hbonds_ldon, hbonds_pdon):
        """"""
        donor_atoms_hbonds = [hb.d.idx for hb in hbonds_ldon + hbonds_pdon]
        wb_dict = {}
        wb_dict2 = {}
        omega = 110.0

        # Just one hydrogen bond per donor atom
        for wbridge in [wb for wb in wbridges if wb.d.idx not in donor_atoms_hbonds]:
            if (wbridge.water.idx, wbridge.a.idx) not in wb_dict:
                wb_dict[(wbridge.water.idx, wbridge.a.idx)] = wbridge
            else:
                if abs(omega - wb_dict[(wbridge.water.idx, wbridge.a.idx)].w_angle) < abs(omega - wbridge.w_angle):
                    wb_dict[(wbridge.water.idx, wbridge.a.idx)] = wbridge
        for wb_tuple in wb_dict:
            water, acceptor = wb_tuple
            if water not in wb_dict2:
                wb_dict2[water] = [(abs(omega - wb_dict[wb_tuple].w_angle), wb_dict[wb_tuple]), ]
            elif len(wb_dict2[water]) == 1:
                wb_dict2[water].append((abs(omega - wb_dict[wb_tuple].w_angle), wb_dict[wb_tuple]))
                wb_dict2[water] = sorted(wb_dict2[water])
            else:
                if wb_dict2[water][1][0] < abs(omega - wb_dict[wb_tuple].w_angle):
                    wb_dict2[water] = [wb_dict2[water][0], (wb_dict[wb_tuple].w_angle, wb_dict[wb_tuple])]

        filtered_wb = []
        for fwbridges in wb_dict2.values():
            [filtered_wb.append(fwb[1]) for fwb in fwbridges]
        return filtered_wb
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');
  }
278376
boolean matchNext() {
        if (numAtoms == 0)
            return false;
        if (sptr > 1)
            backtrack();
        main:
        while (sptr != 0) {
            final int bidx = currBondIdx();

            if (bidx == numBonds) {

                // done
                if (numMapped == numAtoms)
                    return true;

                // handle disconnected atoms
                for (IAtom qatom : query.atoms()) {
                    if (amap[qatom.getIndex()] == UNMAPPED) {
                        Iterator<IAtom> iter = atoms();
                        while (iter.hasNext()) {
                            IAtom atom = iter.next();
                            if (feasible(bidx, (IQueryAtom) qatom, atom))
                                continue main;
                        }
                        break;
                    }
                }
                backtrack();
                continue;
            }

            IQueryBond qbond = qbonds[bidx];
            IQueryAtom qbeg  = (IQueryAtom) qbond.getBegin();
            IQueryAtom qend  = (IQueryAtom) qbond.getEnd();

            int begIdx = amap[qbeg.getIndex()];
            int endIdx = amap[qend.getIndex()];

            // both atoms matched, there must be a bond between them
            if (begIdx != UNMAPPED && endIdx != UNMAPPED) {
                IBond bond = mol.getAtom(begIdx).getBond(mol.getAtom(endIdx));
                if (feasible(qbond, bond))
                    continue;
            }
            // 'beg' is mapped, find a feasible 'end' from it's neighbor list
            else if (begIdx != UNMAPPED) {
                IAtom           beg   = mol.getAtom(begIdx);
                Iterator<IBond> biter = bonds(beg);
                while (biter.hasNext()) {
                    IBond bond = biter.next();
                    IAtom end  = bond.getOther(beg);
                    if (qbond.matches(bond) && feasible(bidx + 1, qend, end))
                        continue main;
                }
            }
            // 'end' is mapped, find a feasible 'beg' from it's neighbor list
            else if (endIdx != UNMAPPED) {
                IAtom           end   = mol.getAtom(endIdx);
                Iterator<IBond> biter = bonds(end);
                while (biter.hasNext()) {
                    IBond bond = biter.next();
                    IAtom beg  = bond.getOther(end);
                    if (qbond.matches(bond) && feasible(bidx + 1, qbeg, beg))
                        continue main;
                }
            }
            // 'beg' nor 'end' matched, find a feasible mapping from
            // any atom in the molecule
            else {
                Iterator<IAtom> aiter = atoms();
                while (aiter.hasNext()) {
                    if (feasible(bidx, qbeg, aiter.next()))
                        continue main;
                }
            }
            backtrack();
        }
        return false;
    }
525667
def select_atoms(indices):
    '''

    '''
    rep = current_representation()
    rep.select({'atoms': Selection(indices, current_system().n_atoms)})
    return rep.selection_state
1072690
def H_donor_count(mol):
    
    mol.require("Valence")
    return sum(1 for _, a in mol.atoms_iter() if a.H_donor)
1166957
private void addBonds(Atom atom, List<Atom> atomsInGroup, List<Atom> allAtoms) {
		if(atom.getBonds()==null){
			return;
		}
		for(Bond bond : atom.getBonds()) {
			// Now set the bonding information.
			Atom other = bond.getOther(atom);
			// If both atoms are in the group
			if (atomsInGroup.indexOf(other)!=-1){
				Integer firstBondIndex = atomsInGroup.indexOf(atom);
				Integer secondBondIndex = atomsInGroup.indexOf(other);
				// Don't add the same bond twice
				if(firstBondIndex>secondBondIndex){
					int bondOrder = bond.getBondOrder();
					mmtfDecoderInterface.setGroupBond(firstBondIndex, secondBondIndex, bondOrder);
				}
			}
			// Otherwise it's an inter group bond - so add it here
			else {
				Integer firstBondIndex = allAtoms.indexOf(atom);
				Integer secondBondIndex = allAtoms.indexOf(other);
				if(firstBondIndex>secondBondIndex){
					// Don't add the same bond twice
					int bondOrder = bond.getBondOrder();
					mmtfDecoderInterface.setInterGroupBond(firstBondIndex, secondBondIndex, bondOrder);
				}
			}
		}
	}
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)
777318
def halogen(acceptor, donor):
    """"""
    data = namedtuple('halogenbond', 'acc acc_orig_idx don don_orig_idx distance don_angle acc_angle restype '
                                     'resnr reschain restype_l resnr_l reschain_l donortype acctype sidechain')
    pairings = []
    for acc, don in itertools.product(acceptor, donor):
        dist = euclidean3d(acc.o.coords, don.x.coords)
        if not config.MIN_DIST < dist < config.HALOGEN_DIST_MAX:
            continue
        vec1, vec2 = vector(acc.o.coords, acc.y.coords), vector(acc.o.coords, don.x.coords)
        vec3, vec4 = vector(don.x.coords, acc.o.coords), vector(don.x.coords, don.c.coords)
        acc_angle, don_angle = vecangle(vec1, vec2), vecangle(vec3, vec4)
        is_sidechain_hal = acc.o.OBAtom.GetResidue().GetAtomProperty(acc.o.OBAtom, 8)  # Check if sidechain atom
        if not config.HALOGEN_ACC_ANGLE - config.HALOGEN_ANGLE_DEV < acc_angle \
                < config.HALOGEN_ACC_ANGLE + config.HALOGEN_ANGLE_DEV:
            continue
        if not config.HALOGEN_DON_ANGLE - config.HALOGEN_ANGLE_DEV < don_angle \
                < config.HALOGEN_DON_ANGLE + config.HALOGEN_ANGLE_DEV:
            continue
        restype, reschain, resnr = whichrestype(acc.o), whichchain(acc.o), whichresnumber(acc.o)
        restype_l, reschain_l, resnr_l = whichrestype(don.orig_x), whichchain(don.orig_x), whichresnumber(don.orig_x)
        contact = data(acc=acc, acc_orig_idx=acc.o_orig_idx, don=don, don_orig_idx=don.x_orig_idx,
                       distance=dist, don_angle=don_angle, acc_angle=acc_angle,
                       restype=restype, resnr=resnr,
                       reschain=reschain, restype_l=restype_l,
                       reschain_l=reschain_l, resnr_l=resnr_l, donortype=don.x.OBAtom.GetType(), acctype=acc.o.type,
                       sidechain=is_sidechain_hal)
        pairings.append(contact)
    return filter_contacts(pairings)
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;
    }
228239
def neighbor_atoms(self, atom_symbol=None):
        
        if not atom_symbol:
            return self.neighbors
        else:
            return [atom for atom in self.neighbors if atom['atom_symbol'] == atom_symbol]
994924
def refine_hydrophobic(self, all_h, pistacks):
        """"""
        sel = {}
        #  1. Rings interacting via stacking can't have additional hydrophobic contacts between each other.
        for pistack, h in itertools.product(pistacks, all_h):
            h1, h2 = h.bsatom.idx, h.ligatom.idx
            brs, lrs = [p1.idx for p1 in pistack.proteinring.atoms], [p2.idx for p2 in pistack.ligandring.atoms]
            if h1 in brs and h2 in lrs:
                sel[(h1, h2)] = "EXCLUDE"
        hydroph = [h for h in all_h if not (h.bsatom.idx, h.ligatom.idx) in sel]
        sel2 = {}
        #  2. If a ligand atom interacts with several binding site atoms in the same residue,
        #  keep only the one with the closest distance
        for h in hydroph:
            if not (h.ligatom.idx, h.resnr) in sel2:
                sel2[(h.ligatom.idx, h.resnr)] = h
            else:
                if sel2[(h.ligatom.idx, h.resnr)].distance > h.distance:
                    sel2[(h.ligatom.idx, h.resnr)] = h
        hydroph = [h for h in sel2.values()]
        hydroph_final = []
        bsclust = {}
        #  3. If a protein atom interacts with several neighboring ligand atoms, just keep the one with the closest dist
        for h in hydroph:
            if h.bsatom.idx not in bsclust:
                bsclust[h.bsatom.idx] = [h, ]
            else:
                bsclust[h.bsatom.idx].append(h)

        idx_to_h = {}
        for bs in [a for a in bsclust if len(bsclust[a]) == 1]:
            hydroph_final.append(bsclust[bs][0])

        # A list of tuples with the idx of an atom and one of its neighbours is created
        for bs in [a for a in bsclust if not len(bsclust[a]) == 1]:
            tuples = []
            all_idx = [i.ligatom.idx for i in bsclust[bs]]
            for b in bsclust[bs]:
                idx = b.ligatom.idx
                neigh = [na for na in pybel.ob.OBAtomAtomIter(b.ligatom.OBAtom)]
                for n in neigh:
                    n_idx = n.GetIdx()
                    if n_idx in all_idx:
                        if n_idx < idx:
                            tuples.append((n_idx, idx))
                        else:
                            tuples.append((idx, n_idx))
                        idx_to_h[idx] = b

            tuples = list(set(tuples))
            tuples = sorted(tuples, key=itemgetter(1))
            clusters = cluster_doubles(tuples)  # Cluster connected atoms (i.e. find hydrophobic patches)

            for cluster in clusters:
                min_dist = float('inf')
                min_h = None
                for atm_idx in cluster:
                    h = idx_to_h[atm_idx]
                    if h.distance < min_dist:
                        min_dist = h.distance
                        min_h = h
                hydroph_final.append(min_h)
        before, reduced = len(all_h), len(hydroph_final)
        if not before == 0 and not before == reduced:
            write_message('Reduced number of hydrophobic contacts from %i to %i.\n' % (before, reduced), indent=True)
        return hydroph_final
1147667
private int getDoubleBondedSulfurCount(IAtomContainer ac, IAtom atom) {
        List<IAtom> neighbours = ac.getConnectedAtomsList(atom);
        IBond bond;
        int sdbcounter = 0;
        for (IAtom neighbour : neighbours) {
            if (neighbour.getSymbol().equals("S")) {
                if (atom.getFormalCharge() == 1 && neighbour.getFormalCharge() == -1) {
                    sdbcounter += 1;
                }
                bond = ac.getBond(neighbour, atom);
                if (!neighbour.getFlag(CDKConstants.ISAROMATIC)) {
                    if (bond.getOrder() == IBond.Order.DOUBLE) {
                        sdbcounter += 1;
                    }
                }
            }
        }
        return sdbcounter;
    }
802872
@Override
    public List<IBond> getConnectedBondsList(IAtom atom) {
        List<IBond> bondsList = new ArrayList<IBond>();
        for (int i = 0; i < bondCount; i++) {
            if (bonds[i].contains(atom)) bondsList.add(bonds[i]);
        }
        return bondsList;
    }
258855
def refine_hbonds_ldon(self, all_hbonds, salt_lneg, salt_pneg):
        """"""
        i_set = {}
        for hbond in all_hbonds:
            i_set[hbond] = False
            for salt in salt_pneg:
                protidx, ligidx = [at.idx for at in salt.negative.atoms], [at.idx for at in salt.positive.atoms]
                if hbond.d.idx in ligidx and hbond.a.idx in protidx:
                    i_set[hbond] = True
            for salt in salt_lneg:
                protidx, ligidx = [at.idx for at in salt.positive.atoms], [at.idx for at in salt.negative.atoms]
                if hbond.d.idx in ligidx and hbond.a.idx in protidx:
                    i_set[hbond] = True

        # Allow only one hydrogen bond per donor, select interaction with larger donor angle
        second_set = {}
        hbls = [k for k in i_set.keys() if not i_set[k]]
        for hbl in hbls:
            if hbl.d.idx not in second_set:
                second_set[hbl.d.idx] = (hbl.angle, hbl)
            else:
                if second_set[hbl.d.idx][0] < hbl.angle:
                    second_set[hbl.d.idx] = (hbl.angle, hbl)
        return [hb[1] for hb in second_set.values()]
800816
@Override
    public boolean isConnectedTo(IBond bond) {
        for (IAtom atom : atoms) {
            if (bond.contains(atom)) return true;
        }
        return false;
    }
995156
private void breakBond(IAtom atom, IAtomContainer molecule) {
        Iterator<IBond> bonds = molecule.bonds().iterator();
        while (bonds.hasNext()) {
            IBond bond = (IBond) bonds.next();
            if (bond.contains(atom)) {
                molecule.removeElectronContainer(bond);
                break;
            }
        }
    }
341700
public IAtom getOtherBondAtom(IAtom atom, IBond bond) {
        if (!bond.contains(atom)) return null;
        if (bond.getBegin().equals(atom))
            return bond.getEnd();
        else
            return bond.getBegin();
    }