Look for halogen bond acceptors (Y-{O|P|N|S}, with Y=C,P,S)
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
278004
private Set<IChemObject> getAromatics(IAtomContainer mol) {
        Set<IChemObject> oldArom = new HashSet<>();
        for (IAtom atom : mol.atoms())
            if (atom.getFlag(CDKConstants.ISAROMATIC))
                oldArom.add(atom);
        for (IBond bond : mol.bonds())
            if (bond.getFlag(CDKConstants.ISAROMATIC))
                oldArom.add(bond);
        return oldArom;
    }
250662
private List<List<Atom[]>> getMatchedAtomsOfLinkages(
			ModificationCondition condition, Map<Component, Set<Group>> mapCompGroups) {
		List<ModificationLinkage> linkages = condition.getLinkages();
		int nLink = linkages.size();

		List<List<Atom[]>> matchedAtomsOfLinkages =
				new ArrayList<List<Atom[]>>(nLink);

		for (int iLink=0; iLink<nLink; iLink++) {
			ModificationLinkage linkage = linkages.get(iLink);
			Component comp1 = linkage.getComponent1();
			Component comp2 = linkage.getComponent2();

//			boolean isAA1 = comp1.;
//			boolean isAA2 = comp2.getType()==true;

			Set<Group> groups1 = mapCompGroups.get(comp1);
			Set<Group> groups2 = mapCompGroups.get(comp2);

			List<Atom[]> list = new ArrayList<Atom[]>();

			List<String> potentialNamesOfAtomOnGroup1 = linkage.getPDBNameOfPotentialAtomsOnComponent1();
			for (String name : potentialNamesOfAtomOnGroup1) {
				if (name.equals("*")) {
					// wildcard
					potentialNamesOfAtomOnGroup1 = null; // search all atoms
					break;
				}
			}

			List<String> potentialNamesOfAtomOnGroup2 = linkage.getPDBNameOfPotentialAtomsOnComponent2();
			for (String name : potentialNamesOfAtomOnGroup2) {
				if (name.equals("*")) {
					// wildcard
					potentialNamesOfAtomOnGroup2 = null; // search all atoms
					break;
				}
			}

			for (Group g1 : groups1) {
				for (Group g2 : groups2) {
					if (g1.equals(g2)) {
						continue;
					}

					// only for wildcard match of two residues
					boolean ignoreNCLinkage =
						potentialNamesOfAtomOnGroup1 == null &&
						potentialNamesOfAtomOnGroup2 == null &&
						residues.contains(g1) &&
						residues.contains(g2);

					Atom[] atoms = StructureUtil.findNearestAtomLinkage(
							g1, g2,
							potentialNamesOfAtomOnGroup1,
							potentialNamesOfAtomOnGroup2,
							ignoreNCLinkage,
							bondLengthTolerance);
					if (atoms!=null) {
						list.add(atoms);
					}
				}
			}

			if (list.isEmpty()) {
				// broken linkage
				break;
			}

			matchedAtomsOfLinkages.add(list);
		}

		return matchedAtomsOfLinkages;
	}
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)
800061
def get_covalent_bonds(self, tol=0.2):
        
        bonds = []
        for site1, site2 in itertools.combinations(self._sites, 2):
            if CovalentBond.is_bonded(site1, site2, tol):
                bonds.append(CovalentBond(site1, site2))
        return bonds
151523
def molecules2symbols(molecules, add_hydrogen=True):
    
    symbols = sorted(
        list(set(
            ase.symbols.string2symbols(''.join(
                map(
                    lambda _x:
                    ''.join(ase.symbols.string2symbols(_x)), molecules)
            ))
        )),
        key=lambda _y: ase.data.atomic_numbers[_y])

    if add_hydrogen and 'H' not in symbols:
        symbols.insert(0, 'H')

    return symbols
333972
private final double[] getAprxBounds(IAtomContainer mol) {
        double maxX = -Double.MAX_VALUE;
        double maxY = -Double.MAX_VALUE;
        double minX = Double.MAX_VALUE;
        double minY = Double.MAX_VALUE;
        IAtom[] boundedAtoms = new IAtom[4];
        for (int i = 0; i < mol.getAtomCount(); i++) {
            IAtom atom = mol.getAtom(i);
            if (atom.getPoint2d() != null) {
                if (atom.getPoint2d().x < minX) {
                    minX = atom.getPoint2d().x;
                    boundedAtoms[0] = atom;
                }
                if (atom.getPoint2d().y < minY) {
                    minY = atom.getPoint2d().y;
                    boundedAtoms[1] = atom;
                }
                if (atom.getPoint2d().x > maxX) {
                    maxX = atom.getPoint2d().x;
                    boundedAtoms[2] = atom;
                }
                if (atom.getPoint2d().y > maxY) {
                    maxY = atom.getPoint2d().y;
                    boundedAtoms[3] = atom;
                }
            }
        }
        double[] minmax = new double[4];
        minmax[0] = minX;
        minmax[1] = minY;
        minmax[2] = maxX;
        minmax[3] = maxY;
        double minXAdjust = adjustForHydrogen(boundedAtoms[0], mol);
        double maxXAdjust = adjustForHydrogen(boundedAtoms[1], mol);
        if (minXAdjust < 0) minmax[0] += minXAdjust;
        if (maxXAdjust > 0) minmax[1] += maxXAdjust;
        return minmax;
    }
1249841
def show_hbonds(self):
        """"""
        grp = self.getPseudoBondGroup("Hydrogen Bonds-%i" % self.tid, associateWith=[self.model])
        grp.lineWidth = 3
        for i in self.plcomplex.hbonds.ldon_id:
            b = grp.newPseudoBond(self.atoms[i[0]], self.atoms[i[1]])
            b.color = self.colorbyname('blue')
            self.bs_res_ids.append(i[0])
        for i in self.plcomplex.hbonds.pdon_id:
            b = grp.newPseudoBond(self.atoms[i[0]], self.atoms[i[1]])
            b.color = self.colorbyname('blue')
            self.bs_res_ids.append(i[1])
872895
def assign_charge(mol, force_recalc=False):
    """"""
    # TODO: not implemented yet
    mol.require("Aromatic")
    for i, nbrs in mol.neighbors_iter():
        atom = mol.atom(i)
        nbrcnt = len(nbrs)
        if atom.symbol == "N":
            if not atom.pi:
                # non-conjugated amines are anion
                mol.atom(i).charge_phys = 1
            elif nbrcnt == 1 and atom.pi == 2:
                # amidine, guanidine are conjugated cation
                ni = list(nbrs.keys())[0]
                conj = False
                sp2n = None
                for nni, nnb in mol.neighbors(ni).items():
                    if mol.atom(nni).symbol == "N" and nnb.order == 2 \
                            and not mol.atom(nni).aromatic:
                        mol.atom(nni).charge_conj = 1
                        conj = True
                    elif mol.atom(nni).symbol == "N" and nni != i:
                        sp2n = nni
                if conj:
                    mol.atom(i).charge_phys = 1
                    if sp2n is not None:
                        mol.atom(sp2n).charge_conj = 1
        elif atom.symbol == "O" and nbrcnt == 1 and atom.pi == 2:
            # oxoacid are conjugated anion
            ni = list(nbrs.keys())[0]
            conj = False
            if mol.atom(ni).symbol == "N":
                mol.atom(i).n_oxide = True
                mol.atom(ni).n_oxide = True
            for nni, nnb in mol.neighbors(ni).items():
                if mol.atom(nni).symbol in ("O", "S") \
                        and nnb.order == 2 and not mol.atom(ni).n_oxide:
                    mol.atom(nni).charge_conj = -1
                    conj = True
            if conj:
                mol.atom(i).charge_phys = -1
        elif atom.symbol == "S" and nbrcnt == 1:
            # thiophenols are anion
            ni = list(nbrs.keys())[0]
            if mol.atom(ni).aromatic:
                mol.atom(i).charge_phys = -1
    mol.charge_assigned = True
    mol.descriptors.add("Phys_charge")
455606
private List<Integer> getAtomNosForRingGroup(IAtomContainer molecule, List<Integer> ringGroup, IRingSet ringSet) {
        List<Integer> atc = new ArrayList<Integer>();
        for (Integer i : ringGroup) {
            for (IAtom atom : ringSet.getAtomContainer(i).atoms()) {
                if (atc.size() > 0) {
                    if (!atc.contains(molecule.indexOf(atom))) {
                        atc.add(molecule.indexOf(atom));
                    }
                } else {
                    atc.add(molecule.indexOf(atom));
                }
            }
        }
        return atc;
    }
1108134
def draw_salt_bridges(self,color="blue"):
		
		self.draw_saltbridges=""
		if self.saltbridges!=None:
			for bond in self.saltbridges.saltbridges_for_drawing:
				self.draw_saltbridges ="<g class='SaltBridges' transform='translate("+str((self.molecule.x_dim-self.molecule.molsize1)/2)+","+str((self.molecule.y_dim-self.molecule.molsize2)/2)+")'>'"
				atom = self.topology_data.universe.atoms[bond[0]-1] #zero-based index vs one-based index
				residue = (atom.resname, str(atom.resid), atom.segid)
				self.draw_saltbridges=self.draw_saltbridges+"<line x1='"+str(int(self.molecule.nearest_points_coords[residue][0]))+"' y1='"+str(int(self.molecule.nearest_points_coords[residue][1]))+"' x2='"+str(float(self.molecule.ligand_atom_coords_from_diagr[bond[1]][0]))+"' y2='"+str(float(self.molecule.ligand_atom_coords_from_diagr[bond[1]][1]))+"' style='stroke:white;stroke-width:15' />"
				self.draw_saltbridges=self.draw_saltbridges+"<line x1='"+str(int(self.molecule.nearest_points_coords[residue][0]))+"' y1='"+str(int(self.molecule.nearest_points_coords[residue][1]))+"' x2='"+str(float(self.molecule.ligand_atom_coords_from_diagr[bond[1]][0]))+"' y2='"+str(float(self.molecule.ligand_atom_coords_from_diagr[bond[1]][1]))+"' style='stroke:"+color+";stroke-width:4' />"
				self.draw_saltbridges= self.draw_saltbridges+"</g>"
794504
def makeAsp(segID, N, CA, C, O, geo):
    ''''''
    ##R-Group
    CA_CB_length=geo.CA_CB_length
    C_CA_CB_angle=geo.C_CA_CB_angle
    N_C_CA_CB_diangle=geo.N_C_CA_CB_diangle

    CB_CG_length=geo.CB_CG_length
    CA_CB_CG_angle=geo.CA_CB_CG_angle
    N_CA_CB_CG_diangle=geo.N_CA_CB_CG_diangle

    CG_OD1_length=geo.CG_OD1_length
    CB_CG_OD1_angle=geo.CB_CG_OD1_angle
    CA_CB_CG_OD1_diangle=geo.CA_CB_CG_OD1_diangle

    CG_OD2_length=geo.CG_OD2_length
    CB_CG_OD2_angle=geo.CB_CG_OD2_angle
    CA_CB_CG_OD2_diangle=geo.CA_CB_CG_OD2_diangle

    carbon_b= calculateCoordinates(N, C, CA, CA_CB_length, C_CA_CB_angle, N_C_CA_CB_diangle)
    CB= Atom("CB", carbon_b, 0.0 , 1.0, " "," CB", 0,"C")
    carbon_g= calculateCoordinates(N, CA, CB, CB_CG_length, CA_CB_CG_angle, N_CA_CB_CG_diangle)
    CG= Atom("CG", carbon_g, 0.0, 1.0, " ", " CG", 0, "C")
    oxygen_d1= calculateCoordinates(CA, CB, CG, CG_OD1_length, CB_CG_OD1_angle, CA_CB_CG_OD1_diangle)
    OD1= Atom("OD1", oxygen_d1, 0.0, 1.0, " ", " OD1", 0, "O")
    oxygen_d2= calculateCoordinates(CA, CB, CG, CG_OD2_length, CB_CG_OD2_angle, CA_CB_CG_OD2_diangle)
    OD2= Atom("OD2", oxygen_d2, 0.0, 1.0, " ", " OD2", 0, "O")

    ##Create Residue Data Structure
    res= Residue((' ', segID, ' '), "ASP", '    ')
    res.add(N)
    res.add(CA)
    res.add(C)
    res.add(O)
    res.add(CB)
    res.add(CG) 
    res.add(OD1)
    res.add(OD2)
    return res
394538
private static int countAlignedBonds(IAtomContainer mol) {
        final double ref = Math.toRadians(30);
        final double diff = Math.toRadians(1);
        int count = 0;
        for (IBond bond : mol.bonds()) {
            Point2d beg = bond.getBegin().getPoint2d();
            Point2d end = bond.getEnd().getPoint2d();
            if (beg.x > end.x) {
                Point2d tmp = beg;
                beg = end;
                end = tmp;
            }
            Vector2d vec = new Vector2d(end.x - beg.x, end.y - beg.y);
            double angle = Math.atan2(vec.y, vec.x);

            if (Math.abs(angle) - ref < diff) {
                count++;
            }
        }
        return count;
    }
1179936
public double calculateBondPolarizability(IAtomContainer atomContainer, IBond bond) {
        double polarizabilitiy = 0;
        IAtomContainer acH = atomContainer.getBuilder().newInstance(IAtomContainer.class, atomContainer);
        addExplicitHydrogens(acH);
        if (bond.getAtomCount() == 2) {
            polarizabilitiy += getKJPolarizabilityFactor(acH, bond.getBegin());
            polarizabilitiy += getKJPolarizabilityFactor(acH, bond.getEnd());
        }
        return (polarizabilitiy / 2);
    }
1320921
def is_anion_cation_bond(valences, ii, jj):
    
    if valences == 'undefined':
        return True
    if valences[ii] == 0 or valences[jj] == 0:
        return True
    return (valences[ii] > 0 > valences[jj]) or (valences[jj] > 0 > valences[ii])
912680
def makeHis(segID, N, CA, C, O, geo):
    ''''''
    ##R-Group
    CA_CB_length=geo.CA_CB_length
    C_CA_CB_angle=geo.C_CA_CB_angle
    N_C_CA_CB_diangle=geo.N_C_CA_CB_diangle

    CB_CG_length=geo.CB_CG_length
    CA_CB_CG_angle=geo.CA_CB_CG_angle
    N_CA_CB_CG_diangle=geo.N_CA_CB_CG_diangle
    
    CG_ND1_length=geo.CG_ND1_length
    CB_CG_ND1_angle=geo.CB_CG_ND1_angle
    CA_CB_CG_ND1_diangle=geo.CA_CB_CG_ND1_diangle
    
    CG_CD2_length=geo.CG_CD2_length
    CB_CG_CD2_angle=geo.CB_CG_CD2_angle
    CA_CB_CG_CD2_diangle=geo.CA_CB_CG_CD2_diangle
    
    ND1_CE1_length=geo.ND1_CE1_length
    CG_ND1_CE1_angle=geo.CG_ND1_CE1_angle
    CB_CG_ND1_CE1_diangle=geo.CB_CG_ND1_CE1_diangle
    
    CD2_NE2_length=geo.CD2_NE2_length
    CG_CD2_NE2_angle=geo.CG_CD2_NE2_angle
    CB_CG_CD2_NE2_diangle=geo.CB_CG_CD2_NE2_diangle
    
    carbon_b= calculateCoordinates(N, C, CA, CA_CB_length, C_CA_CB_angle, N_C_CA_CB_diangle)
    CB= Atom("CB", carbon_b, 0.0 , 1.0, " "," CB", 0,"C")
    carbon_g= calculateCoordinates(N, CA, CB, CB_CG_length, CA_CB_CG_angle, N_CA_CB_CG_diangle)
    CG= Atom("CG", carbon_g, 0.0, 1.0, " ", " CG", 0, "C")
    nitrogen_d1= calculateCoordinates(CA, CB, CG, CG_ND1_length, CB_CG_ND1_angle, CA_CB_CG_ND1_diangle)
    ND1= Atom("ND1", nitrogen_d1, 0.0, 1.0, " ", " ND1", 0, "N")
    carbon_d2= calculateCoordinates(CA, CB, CG, CG_CD2_length, CB_CG_CD2_angle, CA_CB_CG_CD2_diangle)
    CD2= Atom("CD2", carbon_d2, 0.0, 1.0, " ", " CD2", 0, "C")
    carbon_e1= calculateCoordinates(CB, CG, ND1, ND1_CE1_length, CG_ND1_CE1_angle, CB_CG_ND1_CE1_diangle)
    CE1= Atom("CE1", carbon_e1, 0.0, 1.0, " ", " CE1", 0, "C")
    nitrogen_e2= calculateCoordinates(CB, CG, CD2, CD2_NE2_length, CG_CD2_NE2_angle, CB_CG_CD2_NE2_diangle)
    NE2= Atom("NE2", nitrogen_e2, 0.0, 1.0, " ", " NE2", 0, "N")

    ##Create Residue Data Structure
    res= Residue((' ', segID, ' '), "HIS", '    ')
    res.add(N)
    res.add(CA)
    res.add(C)
    res.add(O)
    res.add(CB)
    res.add(CG)
    res.add(ND1)
    res.add(CD2)
    res.add(CE1)
    res.add(NE2)
    return res
1105619
def _jit_give_bond_array(pos, bond_radii, self_bonding_allowed=False):
        
        n = pos.shape[0]
        bond_array = np.empty((n, n), dtype=nb.boolean)

        for i in range(n):
            for j in range(i, n):
                D = 0
                for h in range(3):
                    D += (pos[i, h] - pos[j, h])**2
                B = (bond_radii[i] + bond_radii[j])**2
                bond_array[i, j] = (B - D) >= 0
                bond_array[j, i] = bond_array[i, j]
        if not self_bonding_allowed:
            for i in range(n):
                bond_array[i, i] = False
        return bond_array