Returns the multiplicity of a defect site within the structure (needed for concentration analysis)
969452
def multiplicity(self):
        
        sga = SpacegroupAnalyzer(self.bulk_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(self.bulk_structure[defindex])
        return len(equivalent_sites)
1097478
def determine_min_cell(disordered_structure):
        

        def lcm(n1, n2):
            """
            Find least common multiple of two numbers
            """
            return n1 * n2 / gcd(n1, n2)

        # assumes all order parameters for a given species are the same
        mag_species_order_parameter = {}
        mag_species_occurrences = {}
        for idx, site in enumerate(disordered_structure):
            if not site.is_ordered:
                op = max(site.species.values())
                # this very hacky bit of code only works because we know
                # that on disordered sites in this class, all species are the same
                # but have different spins, and this is comma-delimited
                sp = str(list(site.species.keys())[0]).split(",")[0]
                if sp in mag_species_order_parameter:
                    mag_species_occurrences[sp] += 1
                else:
                    mag_species_order_parameter[sp] = op
                    mag_species_occurrences[sp] = 1

        smallest_n = []

        for sp, order_parameter in mag_species_order_parameter.items():
            denom = Fraction(order_parameter).limit_denominator(100).denominator
            num_atom_per_specie = mag_species_occurrences[sp]
            n_gcd = gcd(denom, num_atom_per_specie)
            smallest_n.append(lcm(int(n_gcd), denom) / n_gcd)

        return max(smallest_n)
1351268
def _get_disorder_mappings(composition, partitions):
        

        def _get_replacement_dict_from_partition(partition):
            d = {}  # to be passed to Structure.replace_species()
            for sp_list in partition:
                if len(sp_list) > 1:
                    total_occ = sum([composition[sp] for sp in sp_list])
                    merged_comp = {sp: composition[sp] / total_occ for sp in sp_list}
                    for sp in sp_list:
                        d[sp] = merged_comp
            return d

        disorder_mapping = [_get_replacement_dict_from_partition(p)
                            for p in partitions]

        return disorder_mapping
1134973
def get_refined_structure(self):
        
        # Atomic positions have to be specified by scaled positions for spglib.
        lattice, scaled_positions, numbers \
            = spglib.refine_cell(self._cell, self._symprec, self._angle_tol)

        species = [self._unique_species[i - 1] for i in numbers]
        s = Structure(lattice, species, scaled_positions)
        return s.get_sorted_structure()
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
789939
def _read_charge_and_multiplicity(self):
        
        temp_charge = read_pattern(
            self.text, {
                "key": r"\$molecule\s+([\-\d]+)\s+\d"
            },
            terminate_on_match=True).get('key')
        if temp_charge != None:
            self.data["charge"] = int(temp_charge[0][0])
        else:
            temp_charge = read_pattern(
                self.text, {
                    "key": r"Sum of atomic charges \=\s+([\d\-\.\+]+)"
                },
                terminate_on_match=True).get('key')
            if temp_charge == None:
                self.data["charge"] = None
            else:
                self.data["charge"] = int(float(temp_charge[0][0]))

        temp_multiplicity = read_pattern(
            self.text, {
                "key": r"\$molecule\s+[\-\d]+\s+(\d)"
            },
            terminate_on_match=True).get('key')
        if temp_multiplicity != None:
            self.data["multiplicity"] = int(temp_multiplicity[0][0])
        else:
            temp_multiplicity = read_pattern(
                self.text, {
                    "key": r"Sum of spin\s+charges \=\s+([\d\-\.\+]+)"
                },
                terminate_on_match=True).get('key')
            if temp_multiplicity == None:
                self.data["multiplicity"] = 1
            else:
                self.data["multiplicity"] = int(
                    float(temp_multiplicity[0][0])) + 1
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]
910048
public static CeSymmResult analyzeLevel(Atom[] atoms,
			CESymmParameters params) throws StructureException {

		if (atoms.length < 1)
			throw new IllegalArgumentException("Empty Atom array given.");

		CeSymmResult result = align(atoms, params);

		if (result.isRefined()) {
			// STEP 5: symmetry alignment optimization
			if (result.getParams().getOptimization()) {
				try {
					MultipleAlignment msa = result.getMultipleAlignment();
					SymmOptimizer optimizer = new SymmOptimizer(result);
					msa = optimizer.optimize();
					result.setMultipleAlignment(msa);
				} catch (RefinerFailedException e) {
					logger.debug("Optimization failed:" + e.getMessage());
				}
			}
		}
		return result;
	}
1062426
private static boolean multiPointCoverageMultiPoint_(
			MultiPoint _multipointA, MultiPoint _multipointB, double tolerance,
			boolean bPerformWithin, boolean bPerformEquals,
			boolean bPerformOverlaps, ProgressTracker progress_tracker) {
		boolean bPerformContains = false;
		MultiPoint multipoint_a;
		MultiPoint multipoint_b;

		if (_multipointA.getPointCount() > _multipointB.getPointCount()) {
			if (bPerformWithin) {
				bPerformWithin = false;
				bPerformContains = true;
			}

			multipoint_a = _multipointB;
			multipoint_b = _multipointA;
		} else {
			multipoint_a = _multipointA;
			multipoint_b = _multipointB;
		}

		AttributeStreamOfInt8 contained = null;

		if (bPerformEquals || bPerformOverlaps || bPerformContains) {
			contained = new AttributeStreamOfInt8(multipoint_b.getPointCount());

			for (int i = 0; i < multipoint_b.getPointCount(); i++)
				contained.write(i, (byte) 0);
		}

		Envelope2D env_a = new Envelope2D();
		Envelope2D env_b = new Envelope2D();
		Envelope2D envInter = new Envelope2D();
		multipoint_a.queryEnvelope2D(env_a);
		multipoint_b.queryEnvelope2D(env_b);
		env_a.inflate(tolerance, tolerance);
		env_b.inflate(tolerance, tolerance);
		envInter.setCoords(env_a);
		envInter.intersect(env_b);

		Point2D ptA = new Point2D();
		Point2D ptB = new Point2D();

		boolean bWithin = true; // starts off true by default

		QuadTreeImpl quadTreeB = InternalUtils.buildQuadTree(
				(MultiPointImpl) (multipoint_b._getImpl()), envInter);
		QuadTreeImpl.QuadTreeIteratorImpl qtIterB = quadTreeB.getIterator();
		double tolerance_sq = tolerance * tolerance;

		for (int vertex_a = 0; vertex_a < multipoint_a.getPointCount(); vertex_a++) {
			multipoint_a.getXY(vertex_a, ptA);

			if (!envInter.contains(ptA)) {
				if (bPerformEquals || bPerformWithin)
					return false;
				else {
					bWithin = false;
					continue;
				}
			}

			boolean bPtACovered = false;
			env_a.setCoords(ptA.x, ptA.y, ptA.x, ptA.y);
			qtIterB.resetIterator(env_a, tolerance);
			for (int elementHandleB = qtIterB.next(); elementHandleB != -1; elementHandleB = qtIterB
					.next()) {
				int vertex_b = quadTreeB.getElement(elementHandleB);
				multipoint_b.getXY(vertex_b, ptB);

				if (Point2D.sqrDistance(ptA, ptB) <= tolerance_sq) {
					if (bPerformEquals || bPerformOverlaps || bPerformContains)
						contained.write(vertex_b, (byte) 1);

					bPtACovered = true;

					if (bPerformWithin)
						break;
				}
			}

			if (!bPtACovered) {
				bWithin = false;

				if (bPerformEquals || bPerformWithin)
					return false;
			}
		}

		if (bPerformOverlaps && bWithin)
			return false;

		if (bPerformWithin)
			return true;

		for (int i = 0; i < multipoint_b.getPointCount(); i++) {
			if (contained.read(i) == 1) {
				if (bPerformOverlaps)
					return true;
			} else {
				if (bPerformEquals || bPerformContains)
					return false;
			}
		}

		if (bPerformOverlaps)
			return false;

		return true;
	}
72637
def recognize_nx(mol):
    
    mol.rings = nx.cycle_basis(mol.graph)
    mol.isolated = sorted(nx.connected_components(mol.graph),
                          key=len, reverse=True)[1:]
    mol.scaffolds = [x for x in nx.biconnected_components(mol.graph)
                     if len(x) > 2]
631194
def get_oxidation_state_decorated_structure(self):
        
        structure = self.chgcar.structure
        charges = [-self.get_charge_transfer(i)
                   for i in range(len(structure))]
        structure.add_oxidation_state_by_site(charges)
        return structure
346405
private static boolean multiPointIntersectsMultiPoint_(
			MultiPoint _multipointA, MultiPoint _multipointB, double tolerance,
			ProgressTracker progress_tracker) {
		MultiPoint multipoint_a;
		MultiPoint multipoint_b;

		if (_multipointA.getPointCount() > _multipointB.getPointCount()) {
			multipoint_a = _multipointB;
			multipoint_b = _multipointA;
		} else {
			multipoint_a = _multipointA;
			multipoint_b = _multipointB;
		}

		Envelope2D env_a = new Envelope2D();
		Envelope2D env_b = new Envelope2D();
		Envelope2D envInter = new Envelope2D();
		multipoint_a.queryEnvelope2D(env_a);
		multipoint_b.queryEnvelope2D(env_b);
		env_a.inflate(tolerance, tolerance);
		env_b.inflate(tolerance, tolerance);
		envInter.setCoords(env_a);
		envInter.intersect(env_b);

		Point2D ptA = new Point2D();
		Point2D ptB = new Point2D();
		double tolerance_sq = tolerance * tolerance;

		QuadTreeImpl quadTreeB = InternalUtils.buildQuadTree(
				(MultiPointImpl) (multipoint_b._getImpl()), envInter);
		QuadTreeImpl.QuadTreeIteratorImpl qtIterB = quadTreeB.getIterator();

		for (int vertex_a = 0; vertex_a < multipoint_a.getPointCount(); vertex_a++) {
			multipoint_a.getXY(vertex_a, ptA);

			if (!envInter.contains(ptA))
				continue;

			env_a.setCoords(ptA.x, ptA.y, ptA.x, ptA.y);
			qtIterB.resetIterator(env_a, tolerance);

			for (int elementHandleB = qtIterB.next(); elementHandleB != -1; elementHandleB = qtIterB
					.next()) {
				int vertex_b = quadTreeB.getElement(elementHandleB);
				multipoint_b.getXY(vertex_b, ptB);

				if (Point2D.sqrDistance(ptA, ptB) <= tolerance_sq)
					return true;
			}
		}

		return false;
	}
796743
def add_oxidation_state_by_site_fraction(structure, oxidation_states):
        
        try:
            for i, site in enumerate(structure):
                new_sp = collections.defaultdict(float)
                for j, (el, occu) in enumerate(get_z_ordered_elmap(site
                        .species)):
                    specie = Specie(el.symbol, oxidation_states[i][j])
                    new_sp[specie] += occu
                structure[i] = new_sp
            return structure
        except IndexError:
            raise ValueError("Oxidation state of all sites must be "
                             "specified in the list.")
385019
def _process_structures(structures):
        
        input_structures = structures
        structures = [input_structures[0]]
        for s in input_structures[1:]:
            prev = structures[-1]
            for i in range(len(s)):
                t = np.round(prev[i].frac_coords - s[i].frac_coords)
                if np.any(np.abs(t) > 0.5):
                    s.translate_sites([i], t, to_unit_cell=False)
            structures.append(s)
        return structures
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
1374150
def get_zmatrix(self):
        
        output = []
        outputvar = []
        for i, site in enumerate(self._mol):
            if i == 0:
                output.append("{}".format(site.specie))
            elif i == 1:
                nn = self._find_nn_pos_before_site(i)
                bondlength = self._mol.get_distance(i, nn[0])
                output.append("{} {} B{}".format(self._mol[i].specie,
                                                 nn[0] + 1, i))
                outputvar.append("B{}={:.6f}".format(i, bondlength))
            elif i == 2:
                nn = self._find_nn_pos_before_site(i)
                bondlength = self._mol.get_distance(i, nn[0])
                angle = self._mol.get_angle(i, nn[0], nn[1])
                output.append("{} {} B{} {} A{}".format(self._mol[i].specie,
                                                        nn[0] + 1, i,
                                                        nn[1] + 1, i))
                outputvar.append("B{}={:.6f}".format(i, bondlength))
                outputvar.append("A{}={:.6f}".format(i, angle))
            else:
                nn = self._find_nn_pos_before_site(i)
                bondlength = self._mol.get_distance(i, nn[0])
                angle = self._mol.get_angle(i, nn[0], nn[1])
                dih = self._mol.get_dihedral(i, nn[0], nn[1], nn[2])
                output.append("{} {} B{} {} A{} {} D{}"
                              .format(self._mol[i].specie, nn[0] + 1, i,
                                      nn[1] + 1, i, nn[2] + 1, i))
                outputvar.append("B{}={:.6f}".format(i, bondlength))
                outputvar.append("A{}={:.6f}".format(i, angle))
                outputvar.append("D{}={:.6f}".format(i, dih))
        return "\n".join(output) + "\n\n" + "\n".join(outputvar)
2388
private int nAdjacentCentres(int i) {
        int n = 0;
        for (IAtom atom : tetrahedralElements[i].getLigands())
            if (tetrahedralElements[atomToIndex.get(atom)] != null) n++;
        return n;
    }