Assign aromatic ring
    sp2 atom:
    pi=1 -> +1
    N, O, S, C- -> +2
    >C=O, B, C+ -> 0
    sp3 atom -> not aromatic
    sum of the score satisfies 4n+2 -> aromatic
505072
def assign_aromatic(mol):
    
    mol.require("Valence")
    mol.require("MinifiedRing")
    for ring in mol.rings:
        pi_cnt = 0
        for r in ring:
            if mol.atom(r).pi == 0:
                if mol.atom(r).symbol == "C":
                    if mol.atom(r).charge == 1:
                        pass
                    elif mol.atom(r).charge == -1:
                        pi_cnt += 2
                    else:
                        break
                elif mol.atom(r).charge == 0:
                    if mol.atom(r).symbol in ("N", "O", "S"):
                        pi_cnt += 2
                    elif mol.atom(r).symbol == "B":
                        pass
                    else:
                        break
                else:
                    break
            elif mol.atom(r).pi == 1:
                if mol.atom(r).carbonyl_C:
                    pass
                else:
                    pi_cnt += 1
            else:
                break
        else:
            if pi_cnt % 4 == 2:
                for u, v in iterator.consecutive(ring + [ring[0]], 2):
                    mol.atom(u).aromatic = True
                    mol.bond(u, v).aromatic = True
    mol.descriptors.add("Aromatic")
522718
private int prepare(IAtom atom, IBond prev) {
        int count = 0;
        amap[atom.getIndex()] = 1;
        for (IBond bond : atom.bonds()) {
            if (bond == prev)
                continue;
            IAtom nbr = bond.getOther(atom);
            if (amap[nbr.getIndex()] == 0) {
                qbonds[numBonds++] = (IQueryBond) bond;
                count += prepare(nbr, bond) + 1;
            } else if (nbr.getIndex() < atom.getIndex()) {
                ++count; // ring closure
                qbonds[numBonds++] = (IQueryBond) bond;
            }
        }
        return count;
    }
537524
def aromatize(molecule, usedPyroles=None):    
    """"""
    
    pyroleLike = getPyroleLikeAtoms(molecule)

    if usedPyroles is None:
        usedPyroles = {}
        
    cyclesToCheck = []
    # determine which cycles came in marked as aromatic
    # and which need to be checked form the kekular form
    #
    # if a cycle came in as aromatic, convert it
    # before going on.
    for cycle in molecule.cycles:
        for atom in cycle.atoms:
            if not atom.aromatic:
                cyclesToCheck.append(cycle)
                break
        else:
            if not convert(cycle, pyroleLike, usedPyroles):
                # XXX FIX ME
                # oops, an aromatic ring came in but
                # we can't convert it.  This is an error
                # daylight would conjugate the ring
                raise PinkyError("Bad initial aromaticity")

    # keep checking rings until something happens
    while 1:
        # assume nothing happened
        needToCheckAgain = 0

        _cyclesToCheck = []
        for cycle in cyclesToCheck:
            canAromatic = canBeAromatic(cycle, pyroleLike)
            if canAromatic == NEVER:
                # the ring can NEVER EVER be aromatic, so remove it for good
                pass
            elif canAromatic and convert(cycle, pyroleLike, usedPyroles):
                needToCheckAgain = 1
            else:
                _cyclesToCheck.append(cycle)

        cyclesToCheck = _cyclesToCheck
        if not needToCheckAgain:
            break

    # fix bonds that have no bondorder if necessary
    molecule = fixBonds(molecule, pyroleLike)
    # add implicit hydrogens
    return addHydrogens(molecule, usedPyroles)
1087558
public static void markAromaticRings(IRing ring) {
        // return as soon as the conditions are not met:
        // 1. all atoms are labeled aromatic
        for (IAtom atom : ring.atoms())
            if (!atom.getFlag(CDKConstants.ISAROMATIC)) return;

        // 2. all bonds are labeled aromatic
        for (IBond bond : ring.bonds())
            if (!bond.getFlag(CDKConstants.ISAROMATIC)) return;

        // OK, all conditions are met, so ring is aromatic
        ring.setFlag(CDKConstants.ISAROMATIC, true);
    }
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;
    }
385247
def to_rdmol(mol):
    """"""
    rwmol = Chem.RWMol(Chem.MolFromSmiles(''))
    key_to_idx = {}
    bond_type = {1: Chem.BondType.SINGLE,
                 2: Chem.BondType.DOUBLE,
                 3: Chem.BondType.TRIPLE}
    conf = Chem.Conformer(rwmol.GetNumAtoms())
    for k, a in mol.atoms_iter():
        i = rwmol.AddAtom(Chem.Atom(atom_number(a.symbol)))
        key_to_idx[k] = i
        conf.SetAtomPosition(i, a.coords)
    rwmol.AddConformer(conf)
    for u, v, b in mol.bonds_iter():
        ui = key_to_idx[u]
        vi = key_to_idx[v]
        rwmol.AddBond(ui, vi, bond_type[b.order])
    Chem.GetSSSR(rwmol)  # Ring recognition is required for fingerprint
    rwmol.UpdatePropertyCache(strict=False)
    return rwmol.GetMol()
160961
@Override
    public DescriptorValue calculate(IAtomContainer atomContainer) {

        int n = 0; // count all heavy atom - heavy atom bonds
        for (IBond bond : atomContainer.bonds()) {
            if (bond.getBegin().getAtomicNumber() != 1 && bond.getEnd().getAtomicNumber() != 1) {
                n++;
            }
        }

        double vadjMa = 0;
        if (n > 0) {
            vadjMa += (Math.log(n) / Math.log(2)) + 1;
        }
        return new DescriptorValue(getSpecification(), getParameterNames(), getParameters(), new DoubleResult(vadjMa),
                getDescriptorNames());
    }
953326
def pymatgen_mol(self):
        
        sp = []
        coords = []
        for atom in ob.OBMolAtomIter(self._obmol):
            sp.append(atom.GetAtomicNum())
            coords.append([atom.GetX(), atom.GetY(), atom.GetZ()])
        return Molecule(sp, coords)
999
private IRingSet recoverRingSystem(IAtomContainer mol) {
        IRingSet ringSet = mol.getBuilder().newInstance(IRingSet.class);
        for (Integer[] bondNumbers : listOfRings) {
            IRing ring = mol.getBuilder().newInstance(IRing.class, bondNumbers.length);
            for (int bondNumber : bondNumbers) {
                IBond bond = mol.getBond(bondNumber);
                ring.addBond(bond);
                if (!ring.contains(bond.getBegin())) ring.addAtom(bond.getBegin());
                if (!ring.contains(bond.getEnd())) ring.addAtom(bond.getEnd());
            }
            ringSet.addAtomContainer(ring);
        }
        return ringSet;
    }
1282083
def detect_aromatic_rings_in_ligand(self):
        """"""
        self.ligrings = {}
        try:
            ring_info = self.topology_data.mol.GetRingInfo()
            self.ligand_ring_num = ring_info.NumRings()
        except Exception as e:
            m = Chem.MolFromPDBFile("lig.pdb")
            ring_info = m.GetRingInfo()
            self.ligand_ring_num = ring_info.NumRings()
        i=0
        for ring in range(self.ligand_ring_num):
            if 4 < len(ring_info.AtomRings()[ring]) <= 6 and  False not in [self.topology_data.mol.GetAtomWithIdx(x).GetIsAromatic() for x in ring_info.AtomRings()[ring]]: #narrow ring definition
                atom_ids_in_ring = []
                for atom in ring_info.AtomRings()[ring]:
                    atom_ids_in_ring.append(self.topology_data.universe.ligand.atoms[atom].name)
                self.ligrings[i]=atom_ids_in_ring
                i+=1
1113868
@Override
    public DescriptorValue calculate(IAtomContainer container) {
        int c1sp1 = 0;
        int c2sp1 = 0;
        int c1sp2 = 0;
        int c2sp2 = 0;
        int c3sp2 = 0;
        int c1sp3 = 0;
        int c2sp3 = 0;
        int c3sp3 = 0;
        int c4sp3 = 0;

        for (IAtom atom : container.atoms()) {
            if (!atom.getSymbol().equals("C") && !atom.getSymbol().equals("c")) continue;
            List<IAtom> connectedAtoms = container.getConnectedAtomsList(atom);

            int cc = 0;
            for (IAtom connectedAtom : connectedAtoms) {
                if (connectedAtom.getSymbol().equals("C") || connectedAtom.getSymbol().equals("c")) cc++;
            }

            IBond.Order maxBondOrder = getHighestBondOrder(container, atom);

            if (maxBondOrder == IBond.Order.TRIPLE && cc == 1)
                c1sp1++;
            else if (maxBondOrder == IBond.Order.TRIPLE && cc == 2)
                c2sp1++;
            else if (maxBondOrder == IBond.Order.DOUBLE && cc == 1)
                c1sp2++;
            else if (maxBondOrder == IBond.Order.DOUBLE && cc == 2)
                c2sp2++;
            else if (maxBondOrder == IBond.Order.DOUBLE && cc == 3)
                c3sp2++;
            else if (maxBondOrder == IBond.Order.SINGLE && cc == 1)
                c1sp3++;
            else if (maxBondOrder == IBond.Order.SINGLE && cc == 2)
                c2sp3++;
            else if (maxBondOrder == IBond.Order.SINGLE && cc == 3)
                c3sp3++;
            else if (maxBondOrder == IBond.Order.SINGLE && cc == 4) c4sp3++;
        }

        IntegerArrayResult retval = new IntegerArrayResult(9);
        retval.add(c1sp1);
        retval.add(c2sp1);
        retval.add(c1sp2);
        retval.add(c2sp2);
        retval.add(c3sp2);
        retval.add(c1sp3);
        retval.add(c2sp3);
        retval.add(c3sp3);
        retval.add(c4sp3);

        return new DescriptorValue(getSpecification(), getParameterNames(), getParameters(), retval,
                getDescriptorNames());
    }
288801
def assign_wctype(mol):
    
    mol.require("Aromatic")
    p_block = ["Al", "B", "Si", "Ga", "Ge", "As", "Se", "Sn", "Te", "Pb",
               "Ne", "Ar", "Kr", "Xe", "Rn"]
    d_block = ["Fe", "Co", "Cu", "Zn", "Tc", "Cd", "Pt", "Au", "Hg", "Gd"]
    for i, atom in mol.atoms_iter():
        nbrs = mol.neighbors(i)
        nbrcnt = len(nbrs)
        nbratms = [mol.atom(i) for i, _ in nbrs.items()]
        isnbrarm = any([mol.atom(i).aromatic for i, _ in nbrs.items()])

        if atom.symbol == "C":
            if atom.aromatic:  # Aromatic carbon
                ex = [(mol.atom(i), b) for i, b in nbrs.items()
                      if not b.aromatic]
                if ex:
                    sbst, sbstb = ex[0]
                    cnosx = {"C": "C21", "N": "C22", "O": "C23", "S": "C24",
                             "F": "C14", "Cl": "C15", "Br": "C16", "I": "C17"}
                    if sbst.symbol not in cnosx:
                        atom.wctype = "C13"
                    elif sbst.aromatic:
                        atom.wctype = "C20"
                    elif sbst.symbol in ("C", "N", "O") and sbstb.order == 2:
                        atom.wctype = "C25"
                    else:
                        atom.wctype = cnosx[sbst.symbol]
                else:
                    if nbrcnt == 2:  # No substituent
                        atom.wctype = "C18"
                    elif nbrcnt == 3:  # Bridgehead
                        atom.wctype = "C19"
            elif atom.pi == 0:  # Aliphatic carbon
                hcnopsx = ["H", "C", "N", "O", "P", "S", "F", "Cl", "Br", "I"]
                if not all([a.symbol in hcnopsx for a in nbratms]):
                    atom.wctype = "C27"
                elif isnbrarm:  # Adjacent to aromatic
                    if atom.H_count == 3:
                        if nbratms[0].symbol == "C":
                            atom.wctype = "C8"
                        else:
                            atom.wctype = "C9"
                    else:
                        numh = {0: "C12", 1: "C11", 2: "C10"}
                        atom.wctype = numh[atom.H_count]
                elif all([a.symbol == "C" for a in nbratms]):
                    if atom.H_count >= 2:
                        atom.wctype = "C1"
                    else:
                        atom.wctype = "C2"
                else:
                    if atom.H_count >= 2:
                        atom.wctype = "C3"
                    else:
                        atom.wctype = "C4"
            elif atom.pi == 1:
                dbi = [i for i, b in nbrs.items() if b.order == 2][0]
                dba = mol.atom(dbi)
                if dba.symbol == "C":
                    if isnbrarm:  # C=C adjacent to aromatic
                        atom.wctype = "C26"
                    else:  # C=C
                        atom.wctype = "C6"
                else:  # C=(!C)
                    atom.wctype = "C5"
                    # Overwrite carbonyl oxygens
                    if dba.symbol == "O":
                        if all(a.symbol != "C" for a in nbratms):
                            mol.atom(dbi).wctype = "O11"
                        elif any(a.aromatic for a in nbratms):
                            mol.atom(dbi).wctype = "O10"
                        else:
                            mol.atom(dbi).wctype = "O9"
            elif atom.pi == 2:
                od = max([b.order for _, b in nbrs.items()])
                if od == 3:  # Alkyne
                    atom.wctype = "C7"
                elif od == 2:  # Allene
                    atom.wctype = "C6"
            if atom.wctype is None:
                atom.wctype = "CS"

        elif atom.symbol == "N":
            if atom.charge == 0:
                if atom.aromatic:
                    atom.wctype = "N11"
                elif atom.pi == 0:  # amine
                    amines = {
                        0: {0: "N7", 1: "N2", 2: "N1"},
                        1: {0: "N8", 1: "N4", 2: "N3"}
                    }
                    atom.wctype = amines[isnbrarm][atom.H_count]
                elif atom.pi == 1:  # imine
                    if atom.H_count:
                        atom.wctype = "N5"
                    else:
                        atom.wctype = "N6"
                elif atom.pi == 2:  # nitrile
                    atom.wctype = "N9"
            elif atom.charge > 0:
                if atom.aromatic:
                    atom.wctype = "N12"
                elif atom.pi == 0:  # ammonium
                    if atom.H_count:
                        atom.wctype = "N10"
                    else:
                        atom.wctype = "N13"
                elif atom.pi == 1:  # iminium
                    atom.wctype = "N13"
                elif atom.pi == 2:
                    if any(a.charge < 1 for a in nbratms):  # 1,3-dipole
                        if all(a.symbol == "N" for a in nbratms):  # azide
                            atom.wctype = "N14"
                        else:  # diazo
                            atom.wctype = "N13"
                    else:  # nitrilium
                        atom.wctype = "N14"
            else:  # nitride
                atom.wctype = "N14"
            if atom.wctype is None:
                atom.wctype = "NS"

        elif atom.symbol == "O":
            if atom.wctype:  # skip if already assigned by carbonyl C
                continue
            if atom.aromatic:
                atom.wctype = "O1"
            elif atom.pi == 0 and atom.charge == 0:
                if atom.H_count == 1:
                    if nbratms[0].symbol in ("O", "S"):  # ROOH, RSOH
                        atom.wctype = "O2a"
                    elif nbratms[0].carbonyl_C:  # Carboxylic acid
                        atom.wctype = "O2a"
                    else:  # Alcohol
                        atom.wctype = "O2"
                elif atom.H_count == 2:  # Water
                    atom.wctype = "O2"
                elif isnbrarm:  # Ether adjacent to aromatic
                    atom.wctype = "O4"
                else:  # Ether
                    atom.wctype = "O3"
            elif nbrcnt == 1:
                nbr = nbratms[0]
                if nbr.aromatic:  # O=Aromatic
                    atom.wctype = "O8"
                elif nbr.symbol in ("N", "O"):  # N-oxide, O2, NO
                    atom.wctype = "O5"
                elif nbr.symbol == "S":  # S-oxide
                    atom.wctype = "O6"
                elif nbr.carbonyl_C and atom.charge < 0:  # Carboxylate anion
                    atom.wctype = "O12"
                else:  # ex. Phosphine oxide
                    atom.wctype = "O7"
            if atom.wctype is None:
                atom.wctype = "OS"

        elif atom.symbol in ("F", "Cl", "Br", "I"):
            if atom.charge:
                atom.wctype = "Hal"
            else:
                atom.wctype = atom.symbol
        elif atom.symbol == "P":
            atom.wctype = "P"
        elif atom.symbol == "S":
            if atom.aromatic:
                atom.wctype = "S3"
            elif atom.charge:
                atom.wctype = "S2"
            else:
                atom.wctype = "S1"
        elif atom.symbol in p_block:
            atom.wctype = "Me1"
        elif atom.symbol in d_block:
            atom.wctype = "Me2"
        elif atom.symbol == "H":  # Explicit hydrogen
            nbr = nbratms[0]
            if nbr.symbol == "C":
                atom.wctype = "H1"
            elif nbr.symbol == "N":
                atom.wctype = "H3"
            else:
                atom.wctype = "H2"
        else:
            raise ValueError("Undefined atom symbol: {}".format(atom.symbol))
    mol.descriptors.add("Wildman-Crippen")
1157570
private IRingSet getRingSetOfAtom(List ringSystems, IAtom atom) {
        IRingSet ringSetOfAtom = null;
        for (int i = 0; i < ringSystems.size(); i++) {
            if (((IRingSet) ringSystems.get(i)).contains(atom)) {
                return (IRingSet) ringSystems.get(i);
            }
        }
        return ringSetOfAtom;
    }
295785
private IAtom getNextAtomWithAliphaticUnplacedNeigbors() {
        IBond bond;
        for (int f = 0; f < molecule.getBondCount(); f++) {
            bond = molecule.getBond(f);

            if (bond.getEnd().getFlag(CDKConstants.ISPLACED) && !bond.getBegin().getFlag(CDKConstants.ISPLACED)) {
                return bond.getEnd();
            }

            if (bond.getBegin().getFlag(CDKConstants.ISPLACED) && !bond.getEnd().getFlag(CDKConstants.ISPLACED)) {
                return bond.getBegin();
            }
        }
        return null;
    }
1098914
def pication(rings, pos_charged, protcharged):
    
    data = namedtuple(
        'pication', 'ring charge distance offset type restype resnr reschain restype_l resnr_l reschain_l protcharged')
    pairings = []
    if len(rings) == 0 or len(pos_charged) == 0:
        return pairings
    for ring in rings:
        c = ring.center
        for p in pos_charged:
            d = euclidean3d(c, p.center)
            # Project the center of charge into the ring and measure distance to ring center
            proj = projection(ring.normal, ring.center, p.center)
            offset = euclidean3d(proj, ring.center)
            if not config.MIN_DIST < d < config.PICATION_DIST_MAX or not offset < config.PISTACK_OFFSET_MAX:
                continue
            if type(p).__name__ == 'lcharge' and p.fgroup == 'tertamine':
                # Special case here if the ligand has a tertiary amine, check an additional angle
                # Otherwise, we might have have a pi-cation interaction 'through' the ligand
                n_atoms = [a_neighbor for a_neighbor in OBAtomAtomIter(p.atoms[0].OBAtom)]
                n_atoms_coords = [(a.x(), a.y(), a.z()) for a in n_atoms]
                amine_normal = np.cross(vector(n_atoms_coords[0], n_atoms_coords[1]),
                                        vector(n_atoms_coords[2], n_atoms_coords[0]))
                b = vecangle(ring.normal, amine_normal)
                # Smallest of two angles, depending on direction of normal
                a = min(b, 180 - b if not 180 - b < 0 else b)
                if not a > 30.0:
                    resnr, restype = whichresnumber(ring.atoms[0]), whichrestype(ring.atoms[0])
                    reschain = whichchain(ring.atoms[0])
                    resnr_l, restype_l = whichresnumber(p.orig_atoms[0]), whichrestype(p.orig_atoms[0])
                    reschain_l = whichchain(p.orig_atoms[0])
                    contact = data(ring=ring, charge=p, distance=d, offset=offset, type='regular',
                                   restype=restype, resnr=resnr, reschain=reschain,
                                   restype_l=restype_l, resnr_l=resnr_l, reschain_l=reschain_l,
                                   protcharged=protcharged)
                    pairings.append(contact)
                break
            resnr = whichresnumber(p.atoms[0]) if protcharged else whichresnumber(ring.atoms[0])
            resnr_l = whichresnumber(ring.orig_atoms[0]) if protcharged else whichresnumber(p.orig_atoms[0])
            restype = whichrestype(p.atoms[0]) if protcharged else whichrestype(ring.atoms[0])
            restype_l = whichrestype(ring.orig_atoms[0]) if protcharged else whichrestype(p.orig_atoms[0])
            reschain = whichchain(p.atoms[0]) if protcharged else whichchain(ring.atoms[0])
            reschain_l = whichchain(ring.orig_atoms[0]) if protcharged else whichchain(p.orig_atoms[0])
            contact = data(ring=ring, charge=p, distance=d, offset=offset, type='regular', restype=restype,
                           resnr=resnr, reschain=reschain, restype_l=restype_l, resnr_l=resnr_l,
                           reschain_l=reschain_l, protcharged=protcharged)
            pairings.append(contact)
    return filter_contacts(pairings)
317633
def assign_descriptors(mol):
    
    topology.recognize(mol)
    descriptor.assign_valence(mol)
    descriptor.assign_rotatable(mol)
    topology.minify_ring(mol)
    descriptor.assign_aromatic(mol)
1370701
@Override
    public DescriptorValue calculate(IAtomContainer atomContainer) {
        IAtomContainer cloneContainer = AtomContainerManipulator.removeHydrogens(atomContainer);
        double petitjeanNumber; //weinerPath
        int diameter = PathTools.getMolecularGraphDiameter(cloneContainer);
        int radius = PathTools.getMolecularGraphRadius(cloneContainer);

        if (diameter == 0)
            petitjeanNumber = 0;
        else
            petitjeanNumber = (diameter - radius) / (double) diameter;
        return new DescriptorValue(getSpecification(), getParameterNames(), getParameters(), new DoubleResult(
                petitjeanNumber), getDescriptorNames());
    }