Detect all halogen bonds of the type Y-O...X-C
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)
1324510
def show_cationpi(self):
        """"""
        grp = self.getPseudoBondGroup("Cation-Pi-%i" % self.tid, associateWith=[self.model])
        grp.lineWidth = 3
        grp.lineType = self.chimera.Dash
        for i, cat in enumerate(self.plcomplex.pication):

            m = self.model
            r = m.newResidue("pseudoatoms", " ", 1, " ")
            chargecenter = m.newAtom("CHARGE", self.chimera.Element("CHARGE"))
            x, y, z = cat.charge_center
            chargecenter.setCoord(self.chimera.Coord(x, y, z))
            r.addAtom(chargecenter)

            centroid = m.newAtom("CENTROID", self.chimera.Element("CENTROID"))
            x, y, z = cat.ring_center
            centroid.setCoord(self.chimera.Coord(x, y, z))
            r.addAtom(centroid)

            b = grp.newPseudoBond(centroid, chargecenter)
            b.color = self.colorbyname('orange')

            if cat.protcharged:
                self.bs_res_ids += cat.charge_atoms
            else:
                self.bs_res_ids += cat.ring_atoms
1075233
def find_metal_binding(self, lig_atoms, water_oxygens):
        
        a_set = []
        data = namedtuple('metal_binding', 'atom orig_atom atom_orig_idx type fgroup restype resnr reschain location')
        for oxygen in water_oxygens:
            a_set.append(data(atom=oxygen.oxy, atom_orig_idx=oxygen.oxy_orig_idx, type='O', fgroup='water',
                              restype=whichrestype(oxygen.oxy), resnr=whichresnumber(oxygen.oxy),
                              reschain=whichchain(oxygen.oxy), location='water', orig_atom=self.Mapper.id_to_atom(oxygen.oxy_orig_idx)))
        # #@todo Refactor code
        for a in lig_atoms:
            a_orig_idx = self.Mapper.mapid(a.idx, mtype='ligand', bsid=self.bsid)
            n_atoms = pybel.ob.OBAtomAtomIter(a.OBAtom)  # Neighboring atoms
            # All atomic numbers of neighboring atoms
            n_atoms_atomicnum = [n.GetAtomicNum() for n in pybel.ob.OBAtomAtomIter(a.OBAtom)]
            if a.atomicnum == 8:  # Oxygen
                if n_atoms_atomicnum.count('1') == 1 and len(n_atoms_atomicnum) == 2:  # Oxygen in alcohol (R-[O]-H)
                    a_set.append(data(atom=a, atom_orig_idx=a_orig_idx, type='O', fgroup='alcohol',
                                      restype=self.hetid, resnr=self.position, reschain=self.chain,
                                      location='ligand', orig_atom=self.Mapper.id_to_atom(a_orig_idx)))
                if True in [n.IsAromatic() for n in n_atoms] and not a.OBAtom.IsAromatic():  # Phenolate oxygen
                    a_set.append(data(atom=a, atom_orig_idx=a_orig_idx, type='O', fgroup='phenolate',
                                      restype=self.hetid, resnr=self.position, reschain=self.chain,
                                      location='ligand', orig_atom=self.Mapper.id_to_atom(a_orig_idx)))
            if a.atomicnum == 6:  # It's a carbon atom
                if n_atoms_atomicnum.count(8) == 2 and n_atoms_atomicnum.count(6) == 1:  # It's a carboxylate group
                    for neighbor in [n for n in n_atoms if n.GetAtomicNum() == 8]:
                        neighbor_orig_idx = self.Mapper.mapid(neighbor.GetIdx(), mtype='ligand', bsid=self.bsid)
                        a_set.append(data(atom=pybel.Atom(neighbor), atom_orig_idx=neighbor_orig_idx, type='O',
                                          fgroup='carboxylate',
                                          restype=self.hetid,
                                          resnr=self.position, reschain=self.chain,
                                          location='ligand', orig_atom=self.Mapper.id_to_atom(a_orig_idx)))
            if a.atomicnum == 15:  # It's a phosphor atom
                if n_atoms_atomicnum.count(8) >= 3:  # It's a phosphoryl
                    for neighbor in [n for n in n_atoms if n.GetAtomicNum() == 8]:
                        neighbor_orig_idx = self.Mapper.mapid(neighbor.GetIdx(), mtype='ligand', bsid=self.bsid)
                        a_set.append(data(atom=pybel.Atom(neighbor), atom_orig_idx=neighbor_orig_idx, type='O',
                                          fgroup='phosphoryl',
                                          restype=self.hetid,
                                          resnr=self.position, reschain=self.chain,
                                          location='ligand', orig_atom=self.Mapper.id_to_atom(a_orig_idx)))
                if n_atoms_atomicnum.count(8) == 2:  # It's another phosphor-containing group #@todo (correct name?)
                    for neighbor in [n for n in n_atoms if n.GetAtomicNum() == 8]:
                        neighbor_orig_idx = self.Mapper.mapid(neighbor.GetIdx(), mtype='ligand', bsid=self.bsid)
                        a_set.append(data(atom=pybel.Atom(neighbor), atom_orig_idx=neighbor_orig_idx, type='O',
                                          fgroup='phosphor.other', restype=self.hetid,
                                          resnr=self.position,
                                          reschain=self.chain, location='ligand', orig_atom=self.Mapper.id_to_atom(a_orig_idx)))
            if a.atomicnum == 7:  # It's a nitrogen atom
                if n_atoms_atomicnum.count(6) == 2:  # It's imidazole/pyrrole or similar
                    a_set.append(data(atom=a, atom_orig_idx=a_orig_idx, type='N', fgroup='imidazole/pyrrole',
                                      restype=self.hetid, resnr=self.position, reschain=self.chain,
                                      location='ligand', orig_atom=self.Mapper.id_to_atom(a_orig_idx)))
            if a.atomicnum == 16:  # It's a sulfur atom
                if True in [n.IsAromatic() for n in n_atoms] and not a.OBAtom.IsAromatic():  # Thiolate
                    a_set.append(data(atom=a, atom_orig_idx=a_orig_idx, type='S', fgroup='thiolate',
                                      restype=self.hetid, resnr=self.position, reschain=self.chain,
                                      location='ligand', orig_atom=self.Mapper.id_to_atom(a_orig_idx)))
                if set(n_atoms_atomicnum) == {26}:  # Sulfur in Iron sulfur cluster
                    a_set.append(data(atom=a, atom_orig_idx=a_orig_idx, type='S', fgroup='iron-sulfur.cluster',
                                      restype=self.hetid, resnr=self.position, reschain=self.chain,
                                      location='ligand', orig_atom=self.Mapper.id_to_atom(a_orig_idx)))

        return a_set
975665
def molecule(lines):
    
    count_line = lines[3]
    num_atoms = int(count_line[0:3])
    num_bonds = int(count_line[3:6])
    # chiral_flag = int(count_line[12:15])  # Not used
    # num_prop = int(count_line[30:33])  # "No longer supported"
    compound = Compound()
    compound.graph._node = atoms(lines[4: num_atoms+4])
    compound.graph._adj = bonds(lines[num_atoms+4: num_atoms+num_bonds+4],
                                compound.graph._node.keys())
    props = properties(lines[num_atoms+num_bonds+4:])
    add_properties(props, compound)
    return compound
1191855
private int getDoubleBondedOxygenCount(IAtomContainer ac, IAtom atom) {
        List<IAtom> neighbours = ac.getConnectedAtomsList(atom);
        IBond bond;
        int odbcounter = 0;
        boolean chargeFlag = false;
        if (atom.getFormalCharge() >= 1) {
            chargeFlag = true;
        }
        for (IAtom neighbour : neighbours) {
            if (neighbour.getSymbol().equals("O")) {
                bond = ac.getBond(neighbour, atom);
                if (chargeFlag && neighbour.getFormalCharge() == -1 && bond.getOrder() == IBond.Order.SINGLE) {
                    odbcounter += 1;
                }
                if (!neighbour.getFlag(CDKConstants.ISAROMATIC)) {
                    if (bond.getOrder() == IBond.Order.DOUBLE) {
                        odbcounter += 1;
                    }
                }
            }
        }
        return odbcounter;
    }
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');
  }
615144
def check_reaction(reactants, products):
    
    reactant_list = [reactant.split('@')[0].strip(
        'star').strip('gas') for reactant in reactants]
    product_list = [product.split('@')[0].strip(
        'star').strip('gas') for product in products]

    reactant_atoms = [extract_atoms(reactant) for reactant in reactant_list]
    product_atoms = [extract_atoms(product) for product in product_list]

    reactants = add_atoms(reactant_atoms)
    products = add_atoms(product_atoms)

    r_stars = 0
    p_stars = 0

    for i, a in enumerate(reactant_atoms):
        if a == '' or 'star' in reactant_list[i]:
            r_stars += 1
        elif isinstance(a, float):
            r_stars += a
    for a in product_atoms:
        if a == '':
            p_stars += 1
        elif isinstance(a, float):
            p_stars += a
    assert ''.join(sorted(reactants)) == ''.join(sorted(products))
101237
def roc_calculator(screened_molecules, status_field, active_label, decoy_label):
    

    P = 0  # Total no. of actives
    N = 0  # Total no. of decoys
    tpf = [];
    tpf.append(0)  # true positive fraction list
    fpf = [];
    fpf.append(0)  # false positive fraction list
    fpindex = []  # indeces where decoys are found are labeled '1'

    # Tally the # of positives & negatives at each threshold & in total
    for index in range(len(screened_molecules)):
        if screened_molecules[index].GetProp(status_field) == active_label and index == 0:
            tpf[index] = float(1)
            P = P + 1
            fpindex.append(0)
        elif screened_molecules[index].GetProp(status_field) == active_label and index > 0:
            tpf.append(float(tpf[index - 1] + 1))
            fpf.append(float(fpf[index - 1]))
            P = P + 1
            fpindex.append(0)
        elif screened_molecules[index].GetProp(status_field) == decoy_label and index == 0:
            fpf[index] = float(1)
            N = N + 1
            fpindex.append(1)
        elif screened_molecules[index].GetProp(status_field) == decoy_label and index > 0:
            fpf.append(float(fpf[index - 1] + 1))
            tpf.append(float(tpf[index - 1]))
            N = N + 1
            fpindex.append(1)

    # calculate TPF & FPF
    for index in range(len(tpf)):
        tpf[index] = tpf[index] / P
        fpf[index] = fpf[index] / N

    return tpf, fpf, P, N
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
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;
    }
646148
def _is_from_chemical_system(chemical_system, struct):
        
        chemsys = list(set([sp.symbol for sp in struct.composition]))
        if len(chemsys) != len(chemical_system):
            return False
        for el in chemsys:
            if not el in chemical_system:
                return False
        return True
573931
private int getHalogenCount(IAtomContainer ac, IAtom atom) {
        List<IAtom> neighbours = ac.getConnectedAtomsList(atom);
        int acounter = 0;
        for (IAtom neighbour : neighbours) {
            if (neighbour.getSymbol().equals("F") || neighbour.getSymbol().equals("I")
                    || neighbour.getSymbol().equals("Cl") || neighbour.getSymbol().equals("Br")) {
                acounter += 1;
            }
        }
        return acounter;
    }
31140
private void identifyBonds() {
        IAtomContainer spt = getSpanningTree();
        IRing ring;
        int nBasicRings = 0;
        for (int i = 0; i < totalEdgeCount; i++) {
            if (!bondsInTree[i]) {
                ring = getRing(spt, molecule.getBond(i));
                for (int b = 0; b < ring.getBondCount(); b++) {
                    int m = molecule.indexOf(ring.getBond(b));
                    cb[nBasicRings][m] = 1;
                }
                nBasicRings++;
            }
        }
        bondsAcyclicCount = 0;
        bondsCyclicCount = 0;
        for (int i = 0; i < totalEdgeCount; i++) {
            int s = 0;
            for (int j = 0; j < nBasicRings; j++) {
                s += cb[j][i];
            }
            switch (s) {
                case (0): {
                    bondsAcyclicCount++;
                    break;
                }
                case (1): {
                    bondsCyclicCount++;
                    break;
                }
                default: {
                    bondsCyclicCount++;
                }
            }
        }
        identifiedBonds = true;
    }
80235
def find_unpaired_ligand(self):
        
        unpaired_hba, unpaired_hbd, unpaired_hal = [], [], []
        # Unpaired hydrogen bond acceptors/donors in ligand (not used for hydrogen bonds/water, salt bridges/mcomplex)
        involved_atoms = [hbond.a.idx for hbond in self.hbonds_pdon] + [hbond.d.idx for hbond in self.hbonds_ldon]
        [[involved_atoms.append(atom.idx) for atom in sb.negative.atoms] for sb in self.saltbridge_lneg]
        [[involved_atoms.append(atom.idx) for atom in sb.positive.atoms] for sb in self.saltbridge_pneg]
        [involved_atoms.append(wb.a.idx) for wb in self.water_bridges if wb.protisdon]
        [involved_atoms.append(wb.d.idx) for wb in self.water_bridges if not wb.protisdon]
        [involved_atoms.append(mcomplex.target.atom.idx) for mcomplex in self.metal_complexes
         if mcomplex.location == 'ligand']
        for atom in [hba.a for hba in self.ligand.get_hba()]:
            if atom.idx not in involved_atoms:
                unpaired_hba.append(atom)
        for atom in [hbd.d for hbd in self.ligand.get_hbd()]:
            if atom.idx not in involved_atoms:
                unpaired_hbd.append(atom)

        # unpaired halogen bond donors in ligand (not used for the previous + halogen bonds)
        [involved_atoms.append(atom.don.x.idx) for atom in self.halogen_bonds]
        for atom in [haldon.x for haldon in self.ligand.halogenbond_don]:
            if atom.idx not in involved_atoms:
                unpaired_hal.append(atom)

        return unpaired_hba, unpaired_hbd, unpaired_hal
1220538
def parse_results_mol2(mol2_outpath):
    
    docked_ligands = pd.DataFrame()

    lines = [line.strip() for line in open(mol2_outpath, 'r')]
    props = {}

    for i, line in enumerate(lines):
        if line.startswith('########## Name:'):
            ligand = line.strip().strip('##########').replace(' ', '').replace('\t', '').split(':')[1]
            line = lines[i + 1]
            props = {}
            props['Ligand'] = ligand
        if line.startswith('##########'):
            splitter = line.strip().strip('##########').replace(' ', '').replace('\t', '').split(':')
            props[splitter[0]] = float(splitter[1])
        if line.startswith('@<TRIPOS>MOLECULE'):
            if props:
                docked_ligands = docked_ligands.append(props, ignore_index=True)

    return docked_ligands
40469
def calculate_descriptors(self,mol):
		""""""
		#make dictionary
		self.ligand_atoms = {index:{"name":x.name} for index,x in enumerate(self.topology_data.universe.ligand_noH.atoms)}

		#Calculate logP and MR
		contribs = self.calculate_logP(mol)

		#Calculate Gasteiger charges
		self.calculate_Gasteiger_charges(mol)

		#Calculate formal charges
		fcharges = self.calculate_formal_charge(mol)

		for atom in self.ligand_atoms.keys():
			self.ligand_atoms[atom]["logP"]=contribs[atom][0]
			self.ligand_atoms[atom]["MR"]=contribs[atom][1]
			self.ligand_atoms[atom]["Gasteiger_ch"]=mol.GetAtomWithIdx(atom).GetProp("_GasteigerCharge")
			self.ligand_atoms[atom]["Formal charges"]=fcharges[atom]

		#Determine rotatable bonds
		self.rot_bonds=self.get_rotatable_bonds(mol)
1265365
def from_indigo(idgmol, assign_descriptor=True):
    """"""
    mol = Compound()
    for atom in idgmol.iterateAtoms():
        key = atom.index()
        a = Atom(atom.symbol())
        a.coords = list(atom.xyz())
        mol.add_atom(key, a)
    for bond in idgmol.iterateBonds():
        u = bond.source()
        v = bond.destination()
        b = Bond()
        b.order = int(bond.bondOrder())
        mol.add_bond(u.index(), v.index(), b)
    if assign_descriptor:
        molutil.assign_descriptors(mol)
    return mol