of effective electronegativity
482253
public double getAtomicSoftnessCore(IAtomContainer ac, int atomPosition) throws CDKException {
        if (factory == null) {
            factory = AtomTypeFactory
                    .getInstance("org/openscience/cdk/config/data/jmol_atomtypes.txt", ac.getBuilder());
        }
        IAtom target = null;
        double core = 0;
        double radiusTarget = 0;
        target = ac.getAtom(atomPosition);
        double partial = 0;
        double radius = 0;
        String symbol = null;
        IAtomType type = null;
        try {
            symbol = ac.getAtom(atomPosition).getSymbol();
            type = factory.getAtomType(symbol);
            if (getCovalentRadius(symbol, ac.getMaximumBondOrder(target)) > 0) {
                radiusTarget = getCovalentRadius(symbol, ac.getMaximumBondOrder(target));
            } else {
                radiusTarget = type.getCovalentRadius();
            }

        } catch (Exception ex1) {
            logger.debug(ex1);
            throw new CDKException("Problems with AtomTypeFactory due to " + ex1.getMessage(), ex1);
        }

        Iterator<IAtom> atoms = ac.atoms().iterator();
        while (atoms.hasNext()) {
            IAtom atom = atoms.next();
            if (!target.equals(atom)) {
                symbol = atom.getSymbol();
                partial = 0;
                try {
                    type = factory.getAtomType(symbol);
                } catch (Exception ex1) {
                    logger.debug(ex1);
                    throw new CDKException("Problems with AtomTypeFactory due to " + ex1.getMessage(), ex1);
                }
                if (getCovalentRadius(symbol, ac.getMaximumBondOrder(atom)) > 0) {
                    radius = getCovalentRadius(symbol, ac.getMaximumBondOrder(atom));
                } else {
                    radius = type.getCovalentRadius();
                }
                partial += radius * radius;
                partial += (radiusTarget * radiusTarget);
                partial = partial / (calculateSquaredDistanceBetweenTwoAtoms(target, atom));
                core += partial;
            }
        }
        core = 2 * core;
        core = 0.172 * core;
        return core;
    }
5182
def basis_comparison_report(bs1, bs2, uncontract_general=False):
    '''
    
    '''

    all_bs1 = list(bs1['elements'].keys())

    if uncontract_general:
        bs1 = manip.uncontract_general(bs1)
        bs2 = manip.uncontract_general(bs2)

    not_in_bs1 = []  # Found in bs2, not in bs1
    not_in_bs2 = all_bs1.copy()  # Found in bs1, not in bs2
    no_diff = []  # Elements for which there is no difference
    some_diff = []  # Elements that are different
    big_diff = []  # Elements that are substantially different

    for k, v in bs2['elements'].items():
        if k not in all_bs1:
            not_in_bs1.append(k)
            continue

        print()
        print("-------------------------------------")
        print(" Element ", k)
        bs1_el = bs1['elements'][k]

        max_rdiff_el = 0.0
        max_rdiff_ecp = 0.0

        # Check to make sure that neither or both have ecp/electron shells
        if 'electron_shells' in v and 'electron_shells' not in bs1_el:
            print("bs2 has electron_shells, but bs1 does not")
            max_rdiff_el = float('inf')
        if 'electron_shells' in bs1_el and 'electron_shells' not in v:
            print("bs1 has electron_shells, but bs2 does not")
            max_rdiff_el = float('inf')
        if 'ecp_potentials' in v and 'ecp_potentials' not in bs1_el:
            print("bs2 has ecp_potentials, but bs1 does not")
            max_rdiff_ecp = float('inf')
        if 'ecp_potentials' in bs1_el and 'ecp_potentials' not in v:
            print("bs1 has ecp_potentials, but bs2 does not")
            max_rdiff_ecp = float('inf')

        if 'electron_shells' in v and 'electron_shells' in bs1_el:
            max_rdiff_el = max(max_rdiff_el, shells_difference(v['electron_shells'], bs1_el['electron_shells']))
        if 'ecp_potentials' in v and 'ecp_potentials' in bs1_el:
            nel1 = v['ecp_electrons']
            nel2 = bs1_el['ecp_electrons']
            if int(nel1) != int(nel2):
                print('Different number of electrons replaced by ECP ({} vs {})'.format(nel1, nel2))
                max_rdiff_ecp = float('inf')
            else:
                max_rdiff_ecp = max(max_rdiff_ecp, potentials_difference(v['ecp_potentials'],
                                                                         bs1_el['ecp_potentials']))

        max_rdiff = max(max_rdiff_el, max_rdiff_ecp)

        # Handle some differences
        if max_rdiff == float('inf'):
            big_diff.append(k)
        elif max_rdiff == 0.0:
            no_diff.append(k)
        else:
            some_diff.append(k)

        not_in_bs2.remove(k)

    print()
    print("     Not in bs1: ", _print_list(not_in_bs1))
    print("     Not in bs2: ", _print_list(not_in_bs2))
    print("  No difference: ", _print_list(no_diff))
    print("Some difference: ", _print_list(some_diff))
    print(" BIG difference: ", _print_list(big_diff))
    print()

    return (len(not_in_bs1) == 0 and len(not_in_bs2) == 0 and len(some_diff) == 0 and len(big_diff) == 0)
1268436
def get_electron_number(self, charge=0):
        
        atomic_number = constants.elements['atomic_number'].to_dict()
        return sum([atomic_number[atom] for atom in self['atom']]) - charge
1111957
public double calculateKJMeanMolecularPolarizability(IAtomContainer atomContainer) {
        double polarizabilitiy = 0;
        IAtomContainer acH = atomContainer.getBuilder().newInstance(IAtomContainer.class, atomContainer);
        addExplicitHydrogens(acH);
        for (int i = 0; i < acH.getAtomCount(); i++) {
            polarizabilitiy += getKJPolarizabilityFactor(acH, acH.getAtom(i));
        }
        return polarizabilitiy;
    }
1372505
def write_nwchem(basis):
    '''
    '''

    # Uncontract all but SP
    basis = manip.uncontract_spdf(basis, 1, True)
    basis = sort.sort_basis(basis, True)

    s = ''

    # Elements for which we have electron basis
    electron_elements = [k for k, v in basis['elements'].items() if 'electron_shells' in v]

    # Elements for which we have ECP
    ecp_elements = [k for k, v in basis['elements'].items() if 'ecp_potentials' in v]

    if len(electron_elements) > 0:
        # basis set starts with a string
        s += 'BASIS "ao basis" PRINT\n'

        # Electron Basis
        for z in electron_elements:
            data = basis['elements'][z]
            sym = lut.element_sym_from_Z(z, True)
            s += '#BASIS SET: {}\n'.format(misc.contraction_string(data))

            for shell in data['electron_shells']:
                exponents = shell['exponents']
                coefficients = shell['coefficients']
                ncol = len(coefficients) + 1

                am = shell['angular_momentum']
                amchar = lut.amint_to_char(am).upper()
                s += '{}    {}\n'.format(sym, amchar)

                point_places = [8 * i + 15 * (i - 1) for i in range(1, ncol + 1)]
                s += printing.write_matrix([exponents, *coefficients], point_places)

        s += 'END\n'

    # Write out ECP
    if len(ecp_elements) > 0:
        s += '\n\nECP\n'

        for z in ecp_elements:
            data = basis['elements'][z]
            sym = lut.element_sym_from_Z(z, True)
            max_ecp_am = max([x['angular_momentum'][0] for x in data['ecp_potentials']])

            # Sort lowest->highest, then put the highest at the beginning
            ecp_list = sorted(data['ecp_potentials'], key=lambda x: x['angular_momentum'])
            ecp_list.insert(0, ecp_list.pop())

            s += '{} nelec {}\n'.format(sym, data['ecp_electrons'])

            for pot in ecp_list:
                rexponents = pot['r_exponents']
                gexponents = pot['gaussian_exponents']
                coefficients = pot['coefficients']

                am = pot['angular_momentum']
                amchar = lut.amint_to_char(am).upper()

                if am[0] == max_ecp_am:
                    s += '{} ul\n'.format(sym)
                else:
                    s += '{} {}\n'.format(sym, amchar)

                point_places = [0, 10, 33]
                s += printing.write_matrix([rexponents, gexponents, *coefficients], point_places)

        s += 'END\n'

    return s
537825
def ξ_bounds(self):
        
        λ, Δμ = self._dichalcogenide.λ, self.Δμ
        return (abs(Δμ) - 2 * λ, abs(Δμ))
344324
void primaryCharges(IAtomContainer mol, int[][] adjList, GraphUtil.EdgeToBondMap edgeMap) {


        for (int v = 0; v < mol.getAtomCount(); v++) {
            IAtom atom = mol.getAtom(v);
            String symbType = atom.getAtomTypeName();
            BigDecimal fc = mmffParamSet.getFormalCharge(symbType);
            
            atom.setCharge(0d);

            if (fc != null) {
                atom.setCharge(fc.doubleValue());
            }
            // charge sharing between equivalent terminal oxygens
            else if (symbType.equals("O2S") || symbType.equals("O3S") || symbType.equals("O2P") || symbType.equals("O3P") || symbType.equals("O4P")) {

                // already handled
                if (atom.getCharge() != 0)
                    continue;

                // find the central atom (S or P)
                int focus = -1;
                for (int w : adjList[v]) {
                    int elem = mol.getAtom(w).getAtomicNumber();
                    if (elem == Elements.Sulfur.number() || elem == Elements.Phosphorus.number()) {
                        if (focus >= 0) {
                            focus = -2;
                            break;
                        }
                        focus = w;
                    }
                }

                // log - multiple or unfound focus
                if (focus < 0)
                    continue;

                // ensure [P+]-[O-] vs P=O are same by including the charge from
                // the focus
                double qSum = fCharge(mol.getAtom(focus));
                int nTerm = 0;

                for (int w : adjList[focus]) {
                    if (mol.getAtom(w).getAtomTypeName().equals(symbType)) {
                        qSum += fCharge(mol.getAtom(w));
                        nTerm++;
                    }
                }
                double qSplt = qSum / nTerm;

                for (int w : adjList[focus]) {
                    if (mol.getAtom(w).getAtomTypeName().equals(symbType)) {
                        atom.setCharge(qSplt);
                    }
                }

            }
            // charge sharing between nitrogen anions 
            else if (symbType.equals("N5M")) {

                if (atom.getCharge() != 0)
                    continue;

                Set<IAtom> eqiv = new HashSet<>();
                Set<Integer> visit = new HashSet<>();
                Deque<Integer> queue = new ArrayDeque<>();

                queue.add(v);

                while (!queue.isEmpty()) {
                    int w = queue.poll();
                    visit.add(w);

                    if (mol.getAtom(w).getAtomTypeName().equals("N5M"))
                        eqiv.add(mol.getAtom(w));

                    for (int u : adjList[w]) {
                        IBond bond = edgeMap.get(w, u);
                        if (bond.getProperty(MMFF_AROM) != null && !visit.contains(u)) {
                            queue.add(u);
                        }
                    }
                }

                double q = 0;
                for (IAtom eqivAtom : eqiv) {
                    q += fCharge(eqivAtom);
                }
                q /= eqiv.size();
                for (IAtom eqivAtom : eqiv) {
                    eqivAtom.setCharge(q);
                }
            }
        }
    }
809629
private void bondAtom(IAtomContainer container, IAtom atom) {
        double myCovalentRadius = atom.getCovalentRadius();
        double searchRadius = myCovalentRadius + maxCovalentRadius + bondTolerance;
        Point tupleAtom = new Point(atom.getPoint3d().x, atom.getPoint3d().y, atom.getPoint3d().z);
        for (Bspt.EnumerateSphere e = bspt.enumHemiSphere(tupleAtom, searchRadius); e.hasMoreElements();) {
            IAtom atomNear = ((TupleAtom) e.nextElement()).getAtom();
            if (!atomNear.equals(atom) && container.getBond(atom, atomNear) == null) {
                boolean bonded = isBonded(myCovalentRadius, atomNear.getCovalentRadius(), e.foundDistance2());
                if (bonded) {
                    IBond bond = atom.getBuilder().newInstance(IBond.class, atom, atomNear, IBond.Order.SINGLE);
                    container.addBond(bond);
                }
            }
        }
    }
292482
def get_edot(F, mc, e):
    

    # chirp mass
    mc *= SOLAR2S

    dedt = -304/(15*mc) * (2*np.pi*mc*F)**(8/3) * e * \
        (1 + 121/304*e**2) / ((1-e**2)**(5/2))

    return dedt
1363852
def Eg(self, **kwargs):
        '''
        
        '''
        return min(self.Eg_Gamma(**kwargs),
                   self.Eg_L(**kwargs),
                   self.Eg_X(**kwargs))
1240068
private double getTopologicalFactors(IAtomContainer atomContainer, IAtomContainer ac) {
        /* factor for separation of charge */
        int totalNCharge1 = AtomContainerManipulator.getTotalNegativeFormalCharge(atomContainer);
        int totalPCharge1 = AtomContainerManipulator.getTotalPositiveFormalCharge(atomContainer);

        double fQ = 1.0;
        if (totalNCharge1 != 0.0) {
            fQ = 0.5;
            for (int i = 0; i < atomContainer.getBondCount(); i++) {
                IBond bond = atomContainer.getBond(i);
                if (bond.getBegin().getFormalCharge() != 0.0 && bond.getEnd().getFormalCharge() != 0.0) {
                    fQ = 0.25;
                    break;
                }
            }
        }
        /* factor, if the number of covalents bonds is decreased */
        double fB = 1.0;

        int numBond1 = 0;
        int numBond2 = 0;
        for (int i = 0; i < atomContainer.getBondCount(); i++) {
            if (atomContainer.getBond(i).getOrder() == IBond.Order.DOUBLE) numBond1 += 1;
            if (ac.getBond(i).getOrder() == IBond.Order.DOUBLE) numBond2 += 1;
        }

        if (numBond1 < /* > */numBond2) fB = 0.8;

        double fPlus = 1.0;
        if (totalNCharge1 == 0.0 && totalPCharge1 == 0.0) fPlus = 0.1;

        /* aromatic */
        double fA = 1.0;
        try {
            if (Aromaticity.cdkLegacy().apply(ac)) if (!Aromaticity.cdkLegacy().apply(atomContainer)) fA = 0.3;
        } catch (CDKException e) {
            e.printStackTrace();
        }
        logger.debug("return " + fQ * fB * fPlus * fA + "= sp:" + fQ + ", dc:" + fB + ", fPlus:" + fPlus + ", fA:" + fA);

        return fQ * fB * fPlus * fA;
    }
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")
1277879
def _get_valences(self):
        
        try:
            bv = BVAnalyzer()
            self._structure = bv.get_oxi_state_decorated_structure(
                self._structure)
            valences = bv.get_valences(self._structure)
        except:
            try:
                bv = BVAnalyzer(symm_tol=0.0)
                self._structure = bv.get_oxi_state_decorated_structure(
                    self._structure)
                valences = bv.get_valences(self._structure)
            except:
                valences = []
                for site in self._structure.sites:
                    if len(site.specie.common_oxidation_states) > 0:
                        valences.append(site.specie.common_oxidation_states[0])
                    # Handle noble gas species
                    # which have no entries in common_oxidation_states.
                    else:
                        valences.append(0)
                if sum(valences):
                    valences = [0] * self._structure.num_sites
                else:
                    self._structure.add_oxidation_state_by_site(valences)
                # raise

        # el = [site.specie.symbol for site in self._structure.sites]
        # el = [site.species_string for site in self._structure.sites]
        # el = [site.specie for site in self._structure.sites]
        # valence_dict = dict(zip(el, valences))
        # print valence_dict
        return valences
70614
def invariants(mol):
    """"""
    atom_ids = {}
    for a in mol.atoms:
        components = []
        components.append(a.number)
        components.append(len(a.oatoms))
        components.append(a.hcount)
        components.append(a.charge)
        components.append(a.mass)
        if len(a.rings) > 0:
            components.append(1)

        atom_ids[a.index] = gen_hash(components)

    return atom_ids
605362
def compare_electron_shells(shell1, shell2, compare_meta=False, rel_tol=0.0):
    '''
     
    '''

    if shell1['angular_momentum'] != shell2['angular_momentum']:
        return False

    # Zip together exponents and coeffs for sorting
    # This basically creates the typical matrix with exponents
    # being in the first column
    tmp1 = list(zip(shell1['exponents'], *shell1['coefficients']))
    tmp2 = list(zip(shell2['exponents'], *shell2['coefficients']))

    # Now sort by first non-zero coefficient
    tmp1 = sorted(tmp1)
    tmp2 = sorted(tmp2)

    if not _compare_matrix(tmp1, tmp2, rel_tol):
        return False
    if compare_meta:
        if shell1['region'] != shell2['region']:
            return False
        if shell1['function_type'] != shell2['function_type']:
            return False
        return True
    else:
        return True
806680
def ycbcr2rgb(y__, cb_, cr_):
    

    kb_ = 0.114
    kr_ = 0.299

    r__ = 2 * cr_ / (1 - kr_) + y__
    b__ = 2 * cb_ / (1 - kb_) + y__
    g__ = (y__ - kr_ * r__ - kb_ * b__) / (1 - kr_ - kb_)

    return r__, g__, b__
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;
    }