Looks for R5, R4, R3 and R2 axes in spherical top molecules.  Point
        group T molecules have only one unique 3-fold and one unique 2-fold
        axis. O molecules have one unique 4, 3 and 2-fold axes. I molecules
        have a unique 5-fold axis.
168758
def _find_spherical_axes(self):
        
        rot_present = defaultdict(bool)
        origin_site, dist_el_sites = cluster_sites(self.centered_mol, self.tol)
        test_set = min(dist_el_sites.values(), key=lambda s: len(s))
        coords = [s.coords for s in test_set]
        for c1, c2, c3 in itertools.combinations(coords, 3):
            for cc1, cc2 in itertools.combinations([c1, c2, c3], 2):
                if not rot_present[2]:
                    test_axis = cc1 + cc2
                    if np.linalg.norm(test_axis) > self.tol:
                        op = SymmOp.from_axis_angle_and_translation(test_axis,
                                                                    180)
                        rot_present[2] = self.is_valid_op(op)
                        if rot_present[2]:
                            self.symmops.append(op)
                            self.rot_sym.append((test_axis, 2))

            test_axis = np.cross(c2 - c1, c3 - c1)
            if np.linalg.norm(test_axis) > self.tol:
                for r in (3, 4, 5):
                    if not rot_present[r]:
                        op = SymmOp.from_axis_angle_and_translation(
                            test_axis, 360 / r)
                        rot_present[r] = self.is_valid_op(op)
                        if rot_present[r]:
                            self.symmops.append(op)
                            self.rot_sym.append((test_axis, r))
                            break
            if rot_present[2] and rot_present[3] and (
                        rot_present[4] or rot_present[5]):
                break
330897
def spin2_a(self):
        """"""
        return coordinates.cartesian_to_spherical_rho(
                                    self.spin1x, self.spin1y, self.spin1z)
871138
def mass_3d(self, r, rho0, Rs):
        
        mass_3d = 2*np.pi*Rs**3*rho0 * r**2/(r + Rs)**2
        return mass_3d
887616
def get_molecule(self, index=0):
        
        return Molecule(self.numbers, self.geometries[index], self.titles[index], symbols=self.symbols)
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
581528
def multiplicity(self):
        
        if self._multiplicity is None:
            # generate multiplicity based on space group symmetry operations performed on defect coordinates
            try:
                d_structure = create_saturated_interstitial_structure(self)
            except ValueError:
                logger.debug('WARNING! Multiplicity was not able to be calculated adequately '
                             'for interstitials...setting this to 1 and skipping for now...')
                return 1

            sga = SpacegroupAnalyzer(d_structure)
            periodic_struc = sga.get_symmetrized_structure()
            poss_deflist = sorted(
                periodic_struc.get_sites_in_sphere(self.site.coords, 2, include_index=True),
                key=lambda x: x[1])
            defindex = poss_deflist[0][2]

            equivalent_sites = periodic_struc.find_equivalent_sites(periodic_struc[defindex])
            return len(equivalent_sites)

        else:
            return self._multiplicity
979174
private void calculateRotationalAxis(Matrix rotation, Atom translation,
			double c) {
		// Calculate magnitude of rotationAxis components, but not signs
		double sum=0;
		double[] rotAx = new double[3];
		for(int i=0;i<3;i++) {
			rotAx[i] = Math.sqrt(rotation.get(i, i)-c);
			sum+=rotAx[i]*rotAx[i];
		}
		for(int i=0;i<3;i++) {
			rotAx[i] /= Math.sqrt(sum);
		}

		// Now determine signs
		double d0 = rotation.get(2,1)-rotation.get(1,2); //=2u[0]*sin(theta)
		double d1 = rotation.get(0,2)-rotation.get(2,0); //=2u[1]*sin(theta)
		double d2 = rotation.get(1,0)-rotation.get(0,1); //=2u[2]*sin(theta)

		double s12 = rotation.get(2,1)+rotation.get(1,2); //=2*u[1]*u[2]*(1-cos(theta))
		double s02 = rotation.get(0,2)+rotation.get(2,0); //=2*u[0]*u[2]*(1-cos(theta))
		double s01 = rotation.get(1,0)+rotation.get(0,1); //=2*u[0]*u[1]*(1-cos(theta))

		//Only use biggest d for the sign directly, for numerical stability around 180deg
		if( Math.abs(d0) < Math.abs(d1) ) { // not d0
			if( Math.abs(d1) < Math.abs(d2) ) { //d2
				if(d2>=0){ //u[2] positive
					if( s02 < 0 ) rotAx[0] = -rotAx[0];
					if( s12 < 0 ) rotAx[1] = -rotAx[1];
				} else { //u[2] negative
					rotAx[2] = -rotAx[2];
					if( s02 >= 0 ) rotAx[0] = -rotAx[0];
					if( s12 >= 0 ) rotAx[1] = -rotAx[1];
				}
			} else { //d1
				if(d1>=0) {//u[1] positive
					if( s01 < 0) rotAx[0] = -rotAx[0];
					if( s12 < 0) rotAx[2] = -rotAx[2];
				} else { //u[1] negative
					rotAx[1] = -rotAx[1];
					if( s01 >= 0) rotAx[0] = -rotAx[0];
					if( s12 >= 0) rotAx[2] = -rotAx[2];
				}
			}
		} else { // not d1
			if( Math.abs(d0) < Math.abs(d2) ) { //d2
				if(d2>=0){ //u[2] positive
					if( s02 < 0 ) rotAx[0] = -rotAx[0];
					if( s12 < 0 ) rotAx[1] = -rotAx[1];
				} else { //u[2] negative
					rotAx[2] = -rotAx[2];
					if( s02 >= 0 ) rotAx[0] = -rotAx[0];
					if( s12 >= 0 ) rotAx[1] = -rotAx[1];
				}
			} else { //d0
				if(d0>=0) { //u[0] positive
					if( s01 < 0 ) rotAx[1] = -rotAx[1];
					if( s02 < 0 ) rotAx[2] = -rotAx[2];
				} else { //u[0] negative
					rotAx[0] = -rotAx[0];
					if( s01 >= 0 ) rotAx[1] = -rotAx[1];
					if( s02 >= 0 ) rotAx[2] = -rotAx[2];
				}
			}
		}

		rotationAxis = new AtomImpl();
		rotationAxis.setCoords(rotAx);

		// Calculate screw = (rotationAxis dot translation)*u
		double dotProduct = Calc.scalarProduct(rotationAxis, translation);

		screwTranslation = Calc.scale(rotationAxis, dotProduct);
		otherTranslation = Calc.subtract(translation, screwTranslation);

		Atom hypot = Calc.vectorProduct(otherTranslation,rotationAxis);
		Calc.scaleEquals(hypot,.5/Math.tan(theta/2.0));

		// Calculate rotation axis position
		rotationPos = Calc.scaleAdd(.5,otherTranslation, hypot);
	}
852942
def rotatable_count(mol):
    
    mol.require("Rotatable")
    return sum(1 for _, _, b in mol.bonds_iter() if b.rotatable)
703466
function() {
    var center = this.center();
    var radiusSquare = 0.0;
    this.eachAtom(function(atom) {
      radiusSquare = Math.max(radiusSquare, vec3.sqrDist(center, atom.pos()));
    });
    return new geom.Sphere(center, Math.sqrt(radiusSquare));
  }
36748
def axes(self):
        """"""
        axes = self.element.findall(GML_NS + 'axis')
        return [Axis(axis.find(GML_NS + 'CoordinateSystemAxis')) for
                axis in axes]
65347
def get_least_orbits(atom_index, cell, site_symmetry, symprec=1e-5):
    """"""
    orbits = _get_orbits(atom_index, cell, site_symmetry, symprec)
    mapping = np.arange(cell.get_number_of_atoms())

    for i, orb in enumerate(orbits):
        for num in np.unique(orb):
            if mapping[num] > mapping[i]:
                mapping[num] = mapping[i]

    return np.unique(mapping)
819269
def calc2dcoords(mol):
    
    topology.recognize(mol)
    g = set(i for i, _ in mol.atoms_iter())
    # 1: get nodes in scaffolds
    scaffolds = []
    belongs = {}
    for i, rkeys in enumerate(sorted(mol.scaffolds, key=len)):
        scf = []
        for rkey in rkeys:
            ring = mol.rings[rkey]
            for r in ring:
                belongs[r] = i
            scf.append(ring)
            g -= set(ring)
        scaffolds.append(scf)
    # 2: traverse nodes and scaffolds
    # the node and scaffold graph should be a tree (no cycles)
    f = True
    # print(scaffolds)
    coords = {}
    while g:
        if f and scaffolds:  # largest scaffold is first
            stack = [scaffolds[-1][0][0]]
            f = False
        else:
            stack = [g.pop()]
        pred = {}
        branch = {}
        while stack:
            # print("stack: {}".format(stack))
            tail = stack.pop()
            # print("tail: {}".format(tail))
            if tail in belongs:  # scaffolds
                scf = scaffold_coords(scaffolds[belongs[tail]])
                # print(scf.keys())
                # rotate and translate
                if not coords:
                    coords = scf
                else:
                    u = coords[pred[tail]]
                    v = scf[tail]
                    op = [u[0] + math.cos(u[2]), u[1] + math.sin(u[2])]
                    translate(scf, gm.vector(v[:2], op))
                    rotate(scf, op, gm.rad(u[2] + math.pi - v[2]))
                    coords.update(scf)
                # stack nbrs of scaffold
                for k in scf.keys():
                    pred[k] = None
                    for nbr in mol.neighbors(k):
                        if nbr not in scf.keys():
                            stack.append(nbr)
                            pred[nbr] = k
            else:  # append linker
                if tail not in pred:  # isolated
                    coords[tail] = [0, 0, 0, 1]
                    continue
                p = pred[tail]
                x, y, ang, d = coords[p]
                # TODO: ring configuration
                coords[tail] = [x + math.cos(ang), y + math.sin(ang),
                                ang + d * math.pi / 3, d * -1]
                if p not in branch:
                    coords[p][2] = gm.rad(coords[p][2] + math.pi * 2 / 3 * d)
                    branch[p] = 1
                elif branch[p] == 1:
                    coords[p][2] = gm.rad(coords[p][2] + math.pi * d)
                    branch[p] += 1
                elif branch[p] == 2:
                    coords[p][2] = gm.rad(coords[p][2] + math.pi * 2 / 3 * d)
                    branch[p] += 1
                for nbr in mol.neighbors(tail):
                    if nbr not in pred:
                        stack.append(nbr)
                        pred[nbr] = tail
        g -= set(pred)
    resolve_overlap(coords)
    for i, a in mol.atoms_iter():
        mol.atom(i).coords = coords[i][:2]
658172
function markDiastereotopicAtoms(molecule) {
  // changed from markDiastereo(); TLS 9.Nov.2015
  let ids = getAtomIDs(molecule);
  let analyzed = {};
  let group = 0;
  for (let id of ids) {
    console.log(`${id} - ${group}`);
    if (!analyzed.contains(id)) {
      analyzed[id] = true;
      for (let iAtom = 0; iAtom < ids.length; iAtom++) {
        if (id.equals(ids[iAtom])) {
          molecule.setAtomCustomLabel(iAtom, group);
        }
      }
      group++;
    }
  }
}
709214
def _get_bonds(self, mol):
        
        num_atoms = len(mol)
        # index starting from 0
        if self.ignore_ionic_bond:
            covalent_atoms = [i for i in range(num_atoms) if mol.species[i].symbol not in self.ionic_element_list]
        else:
            covalent_atoms = list(range(num_atoms))
        all_pairs = list(itertools.combinations(covalent_atoms, 2))
        pair_dists = [mol.get_distance(*p) for p in all_pairs]
        elements = mol.composition.as_dict().keys()
        unavailable_elements = list(set(elements) -
                                    set(self.covalent_radius.keys()))
        if len(unavailable_elements) > 0:
            raise ValueError("The covalent radius for element {} is not "
                             "available".format(unavailable_elements))
        bond_13 = self.get_13_bonds(self.priority_bonds)
        max_length = [(self.covalent_radius[mol.sites[p[0]].specie.symbol] +
                       self.covalent_radius[mol.sites[p[1]].specie.symbol]) *
                      (1 + (self.priority_cap if p in self.priority_bonds
                            else (self.bond_length_cap if p not in bond_13
                                  else self.bond_13_cap))) *
                      (0.1 if (self.ignore_halogen_self_bond and p not in self.priority_bonds and
                               mol.sites[p[0]].specie.symbol in self.halogen_list and
                               mol.sites[p[1]].specie.symbol in self.halogen_list)
                       else 1.0)
                      for p in all_pairs]

        bonds = [bond
                 for bond, dist, cap in zip(all_pairs, pair_dists, max_length)
                 if dist <= cap]
        return bonds
333383
def spin2_polar(self):
        """"""
        return coordinates.cartesian_to_spherical_polar(
                                     self.spin2x, self.spin2y, self.spin2z)
612755
def sssr(molecule):
    """"""
    results = {}
    lookup = {}
    fullSet = {}
    oatoms = {}
    # XXX FIX ME
    # copy atom.oatoms to atom._oatoms
    # atom._oatoms will be modified my the routine

    for atom in molecule.atoms:
        atom.rings = []
	fullSet[atom.handle] = 1
	lookup[atom.handle] = atom
        oatoms[atom.handle] = atom.oatoms[:]

    for bond in molecule.bonds:
        bond.rings = []
        
    trimSet = []

    while fullSet:
	nodesN2 = []
        minimum, minimum_degree = None, 100000


	# find the N2 atoms and remove atoms with degree 0
	for atomID in fullSet.keys():
	    atom = lookup[atomID]
	    degree = len(oatoms[atom.handle])
	    if degree == 0:
		del fullSet[atomID]
		#fullSet.remove(atomID)
	    elif degree == 2:
		nodesN2.append(atom)

	    # keep track of the minimum degree
	    if (degree > 0) and ( (not minimum) or 
		                (degree < minimum_degree)):
		minimum, minimum_degree = atom, degree

	if not minimum:
	    # nothing to do!  (i.e. can't have a ring)
	    break

	if minimum_degree == 1:
	    # these cannot be in rings so trim and remove
	    # my version of trimming
            for oatom in oatoms[minimum.handle]:
		oatoms[oatom.handle].remove(minimum)
            oatoms[minimum.handle] = []
	    del fullSet[minimum.handle]

	elif minimum_degree == 2:
	    # find the rings!
	    startNodes = []
	    for atom in nodesN2:
		ring, bonds = getRing(atom, fullSet, lookup, oatoms)

		if ring:
                    rlookup = ring[:]
                    rlookup.sort()
                    rlookup = tuple(rlookup)
		    if (not results.has_key(rlookup)):# not in results):
			results[rlookup] = ring, bonds
			startNodes.append(atom)

	    # in case we didn't get a ring remove the head of the nodesN2
	    startNodes = startNodes or [nodesN2[0]]
	    for atom in startNodes:
		# again, my version of trimming
                if oatoms[atom.handle]:
		    oatom = oatoms[atom.handle].pop()
		    oatoms[oatom.handle].remove(atom)

	elif minimum_degree > 2:
	    # no N2 nodes so remove the "optimum" edge to create
	    # N2 nodes in the next go-around.
	    ring, bonds = getRing(minimum, fullSet, lookup, oatoms)
            if ring:
                key = ring[:]
                key.sort()
                key = tuple(key)
                if not results.has_key(key):
                    results[key] = ring, bonds
                    atoms = map(lookup.get, ring)
                    atoms, bonds = toposort(atoms, bonds)
                    checkEdges(atoms, lookup, oatoms)
            else:
                del fullSet[minimum.handle]
 	else:
	    raise ShouldntGetHereError

    # assign the ring index to the atom
    rings = []
    index = 0

    # transform the handles back to atoms
    for result, bonds in results.values():
	ring = []
        for atomID in result:
	    atom = lookup[atomID]
            assert atom.handle == atomID
	    ring.append(atom)
	rings.append((ring, bonds))
        index = index + 1

    molecule.rings = rings
    potentialCycles = []
    index = 0
    for atoms, bonds in rings:
        # due to the dictionaries used in getRing
        # the atoms are not in the order found
        # we need to topologically sort these
        # for the cycle
        atoms, bonds = toposort(atoms, bonds)
        potentialCycles.append((atoms, bonds))

    rings = potentialCycles#checkRings(potentialCycles)
    molecule.rings = rings
    molecule.cycles = [Cycle(atoms, bonds) for atoms, bonds in rings]
    return molecule
1271874
def _set_cluster(self):
        
        center = self.struct[self.center_index].coords
        sphere = self.struct.get_neighbors(self.struct[self.center_index], self.radius)

        symbols = [self.absorbing_atom]
        coords = [[0, 0, 0]]
        for i, site_dist in enumerate(sphere):
            site_symbol = re.sub(r"[^aA-zZ]+", "", site_dist[0].species_string)
            symbols.append(site_symbol)
            coords.append(site_dist[0].coords - center)
        return Molecule(symbols, coords)