def test_attributes(self):
        is_true = {("Xe", "Kr") : "is_noble_gas",
                   ("Fe", "Ni") : 'is_transition_metal',
                   ('Li', 'Cs') : 'is_alkali',
                   ('Ca', 'Mg') : 'is_alkaline',
                   ('F', 'Br', 'I') : 'is_halogen',
                   ('La',) : 'is_lanthanoid',
                   ('U', 'Pu') : 'is_actinoid',
                   ('Si', 'Ge') : 'is_metalloid'
                   }

        for k, v in is_true.items():
            for sym in k:
                self.assertTrue(getattr(Element(sym), v), sym + ' is false')

        keys = ["name", "Z", "mendeleev_no", "atomic_mass", "electronic_structure", "X", "atomic_radius", "min_oxidation_state",
                        "max_oxidation_state", "electrical_resistivity", "velocity_of_sound",
                        "reflectivity", "refractive_index", "poissons_ratio", "molar_volume", "thermal_conductivity", "melting_point", "boiling_point",
                        "liquid_range", "critical_temperature", "superconduction_temperature",
                        "bulk_modulus", "youngs_modulus", "brinell_hardness", "rigidity_modulus", "mineral_hardness",
                        "vickers_hardness", "density_of_solid", "coefficient_of_linear_thermal_expansion", "oxidation_states", "common_oxidation_states", 'average_ionic_radius', 'ionic_radii']

        #Test all elements up to Uranium
        for i in range(1, 93):
            for k in keys:
                self.assertIsNotNone(getattr(Element.from_Z(i), k))
            el = Element.from_Z(i)
            if len(el.oxidation_states) > 0:
                self.assertEqual(max(el.oxidation_states), el.max_oxidation_state)
                self.assertEqual(min(el.oxidation_states), el.min_oxidation_state)

    def test_attributes(self):
        is_true = {("Xe", "Kr"): "is_noble_gas",
                   ("Fe", "Ni"): "is_transition_metal",
                   ("Li", "Cs"): "is_alkali",
                   ("Ca", "Mg"): "is_alkaline",
                   ("F", "Br", "I"): "is_halogen",
                   ("La",): "is_lanthanoid",
                   ("U", "Pu"): "is_actinoid",
                   ("Si", "Ge"): "is_metalloid",
                   ("O", "Te"): "is_chalcogen"}

        for k, v in is_true.items():
            for sym in k:
                self.assertTrue(getattr(Element(sym), v), sym + " is false")

        keys = ["mendeleev_no", "atomic_mass",
                "electronic_structure", "atomic_radius",
                "min_oxidation_state", "max_oxidation_state",
                "electrical_resistivity", "velocity_of_sound", "reflectivity",
                "refractive_index", "poissons_ratio", "molar_volume",
                "thermal_conductivity", "melting_point", "boiling_point",
                "liquid_range", "critical_temperature",
                "superconduction_temperature",
                "bulk_modulus", "youngs_modulus", "brinell_hardness",
                "rigidity_modulus", "mineral_hardness",
                "vickers_hardness", "density_of_solid",
                "atomic_orbitals", "coefficient_of_linear_thermal_expansion",
                "oxidation_states",
                "common_oxidation_states", "average_ionic_radius",
                "average_cationic_radius", "average_anionic_radius",
                "ionic_radii", "long_name", "metallic_radius", "iupac_ordering"]

        # Test all elements up to Uranium
        for i in range(1, 104):
            el = Element.from_Z(i)
            d = el.data
            for k in keys:
                k_str = k.capitalize().replace("_", " ")
                if k_str in d and (not str(d[k_str]).startswith("no data")):
                    self.assertIsNotNone(getattr(el, k))
                elif k == "long_name":
                    self.assertEqual(getattr(el, "long_name"), d["Name"])
                elif k == "iupac_ordering":
                    self.assertTrue("IUPAC ordering" in d)
                    self.assertIsNotNone(getattr(el, k))
            el = Element.from_Z(i)
            if len(el.oxidation_states) > 0:
                self.assertEqual(max(el.oxidation_states),
                                 el.max_oxidation_state)
                self.assertEqual(min(el.oxidation_states),
                                 el.min_oxidation_state)

            if el.symbol not in ["He", "Ne", "Ar"]:
                self.assertTrue(el.X > 0, "No electroneg for %s" % el)

        self.assertRaises(ValueError, Element.from_Z, 1000)

 def test_equals(self):
     random_z = random.randint(1, 92)
     fixed_el = Element.from_Z(random_z)
     other_z = random.randint(1, 92)
     while other_z == random_z:
         other_z = random.randint(1, 92)
     comp1 = Composition({fixed_el:1, Element.from_Z(other_z) :0})
     other_z = random.randint(1, 92)
     while other_z == random_z:
         other_z = random.randint(1, 92)
     comp2 = Composition({fixed_el:1, Element.from_Z(other_z) :0})
     self.assertEqual(comp1, comp2, "Composition equality test failed. %s should be equal to %s" % (comp1.formula, comp2.formula))
     self.assertEqual(comp1.__hash__(), comp2.__hash__(), "Hashcode equality test failed!")

    def from_dict(cls, d, lattice=None):
        """
        Create PeriodicSite from dict representation.

        Args:
            d (dict): dict representation of PeriodicSite
            lattice: Optional lattice to override lattice specified in d.
                Useful for ensuring all sites in a structure share the same
                lattice.

        Returns:
            PeriodicSite
        """
        atoms_n_occu = {}
        for sp_occu in d["species"]:
            if "oxidation_state" in sp_occu and Element.is_valid_symbol(
                    sp_occu["element"]):
                sp = Specie.from_dict(sp_occu)
            elif "oxidation_state" in sp_occu:
                sp = DummySpecie.from_dict(sp_occu)
            else:
                sp = Element(sp_occu["element"])
            atoms_n_occu[sp] = sp_occu["occu"]
        props = d.get("properties", None)
        lattice = lattice if lattice else Lattice.from_dict(d["lattice"])
        return cls(atoms_n_occu, d["abc"], lattice, properties=props)

    def test_element(self):
        symbols = list()
        for i in range(1, 102):
            el = Element.from_Z(i)
            self.assertGreater(el.atomic_mass, 0,
                               "Atomic mass cannot be negative!")
            self.assertNotIn(el.symbol, symbols,
                             "Duplicate symbol for " + el.symbol)
            symbols.append(""" + el.symbol + """)
            self.assertIsNotNone(el.group,
                                 "Group cannot be none for Z=" + str(i))
            self.assertIsNotNone(el.row, "Row cannot be none for Z=" + str(i))

            #Test all properties
            all_attr = ["Z", "symbol", "X", "name", "atomic_mass",
                        "atomic_radius", "max_oxidation_state",
                        "min_oxidation_state", "mendeleev_no",
                        "electrical_resistivity", "velocity_of_sound",
                        "reflectivity", "refractive_index", "poissons_ratio",
                        "molar_volume", "electronic_structure",
                        "thermal_conductivity", "boiling_point",
                        "melting_point", "critical_temperature",
                        "superconduction_temperature", "liquid_range",
                        "bulk_modulus", "youngs_modulus", "brinell_hardness",
                        "rigidity_modulus", "mineral_hardness",
                        "vickers_hardness", "density_of_solid",
                        "coefficient_of_linear_thermal_expansion"]

            for a in all_attr:
                self.assertIsNotNone(el, a)

        def _parse_chomp_and_rank(m, f, m_dict, m_points):
            """
            A helper method for formula parsing that helps in interpreting and
            ranking indeterminate formulas
            Author: Anubhav Jain

            Args:
                m: A regex match, with the first group being the element and
                    the second group being the amount
                f: The formula part containing the match
                m_dict: A symbol:amt dictionary from the previously parsed
                    formula
                m_points: Number of points gained from the previously parsed
                    formula

            Returns:
                A tuple of (f, m_dict, points) where m_dict now contains data
                from the match and the match has been removed (chomped) from
                the formula f. The "goodness" of the match determines the
                number of points returned for chomping. Returns
                (None, None, None) if no element could be found...
            """

            points = 0
            # Points awarded if the first element of the element is correctly
            # specified as a capital
            points_first_capital = 100
            # Points awarded if the second letter of the element is correctly
            # specified as lowercase
            points_second_lowercase = 100

            #get element and amount from regex match
            el = m.group(1)
            if len(el) > 2 or len(el) < 1:
                raise CompositionError("Invalid element symbol entered!")
            amt = float(m.group(2)) if m.group(2).strip() != "" else 1

            #convert the element string to proper [uppercase,lowercase] format
            #and award points if it is already in that format
            char1 = el[0]
            char2 = el[1] if len(el) > 1 else ""

            if char1 == char1.upper():
                points += points_first_capital
            if char2 and char2 == char2.lower():
                points += points_second_lowercase

            el = char1.upper() + char2.lower()

            #if it's a valid element, chomp and add to the points
            if Element.is_valid_symbol(el):
                if el in m_dict:
                    m_dict[el] += amt * factor
                else:
                    m_dict[el] = amt * factor
                return f.replace(m.group(), "", 1), m_dict, m_points + points

            #else return None
            return None, None, None

    def test_pickle(self):
        el1 = Element.Fe
        o = pickle.dumps(el1)
        self.assertEqual(el1, pickle.loads(o))

        #Test all elements up to Uranium
        for i in range(1, 93):
            self.serialize_with_pickle(Element.from_Z(i), test_eq=True)

def from_string(string, symbols=None):
    structures = []
    timesteps = []
    steps = string.split("ITEM: TIMESTEP")
    steps.pop(0)
    for step in steps:
        lines = tuple(step.split("\n"))
        mdstep = int(lines[1])
        natoms = int(lines[3])
        xbox   = tuple((lines[5] + " 0").split())
        ybox   = tuple((lines[6] + " 0").split())
        zbox   = tuple((lines[7] + " 0").split())
        xlo = float(xbox[0])
        xhi = float(xbox[1])
        xy  = float(xbox[2])
        ylo = float(ybox[0])
        yhi = float(ybox[1])
        xz  = float(ybox[2])
        zlo = float(zbox[0])
        zhi = float(zbox[1])
        yz  = float(zbox[2])
        xlo -= np.min([0, xy, xz, xy+xz])
        xhi -= np.max([0, xy, xz, xy+xz])
        ylo -= np.min([0, yz])
        yhi -= np.max([0, yz])
      
        lattice = [xhi-xlo, 0, 0, xy, yhi-ylo, 0, xz, yz, zhi-zlo]
        atoms = [[float(j) for j in line.split()] for line in lines[9:-1]]
        atoms = sorted(atoms, key=lambda a_entry: a_entry[0], reverse=True) 

        coords = []
        atomic_sym = []
        while (len(atoms)):
           one = atoms.pop()
           typ = one[1]
           coo = one[2:5]
           sym = symbols[typ] if symbols else Element.from_Z(typ).symbol
           atomic_sym.append(sym)
           coords.append(coo)
        struct = Structure(lattice, atomic_sym, coords,
                   to_unit_cell=False, validate_proximity=False,
                   coords_are_cartesian=False)
        structures.append(struct)
        timesteps.append(mdstep)
        
    return lmpdump(structures, timesteps)

    def test_element(self):
        symbols = list()
        for i in range(1, 102):
            el = Element.from_Z(i)
            self.assertGreater(el.atomic_mass, 0, "Atomic mass cannot be negative!")
            self.assertNotIn(el.symbol, symbols, "Duplicate symbol for " + el.symbol)
            symbols.append('"' + el.symbol + '"')
            self.assertIsNotNone(el.group, "Group cannot be none for Z=" + str(i))
            self.assertIsNotNone(el.row, "Row cannot be none for Z=" + str(i))

            #Test all properties
            all_attr = ['Z', 'symbol', 'X', 'name', 'atomic_mass', 'atomic_radius', 'max_oxidation_state', 'min_oxidation_state', 'mendeleev_no',
            'electrical_resistivity', 'velocity_of_sound', 'reflectivity', 'refractive_index', 'poissons_ratio', 'molar_volume' , 'electronic_structure',
            'thermal_conductivity', 'boiling_point', 'melting_point', 'critical_temperature', 'superconduction_temperature', 'liquid_range', 'bulk_modulus',
            'youngs_modulus', 'brinell_hardness', 'rigidity_modulus', 'mineral_hardness', 'vickers_hardness', 'density_of_solid', 'coefficient_of_linear_thermal_expansion']

            for a in all_attr:
                self.assertIsNotNone(el, a)

def get_transformed_entries(entries, elements):
    comp_matrix = get_comp_matrix(entries, elements)
    newmat = []
    energies = []
    for i in xrange(len(elements)):
        col = comp_matrix[:, i]
        maxval = max(col)
        maxind = list(col).index(maxval)
        newmat.append(comp_matrix[maxind])
        energies.append(entries[i].energy_per_atom)
    invm = np.linalg.inv(np.array(newmat).transpose())
    newentries = []
    for i in xrange(len(entries)):
        entry = entries[i]
        lincomp = np.dot(invm, comp_matrix[i])
        lincomp = np.around(lincomp, 5)
        comp = Composition({Element.from_Z(j + 1):lincomp[j] for j in xrange(len(elements))})
        scaled_energy = entry.energy_per_atom - sum(lincomp * energies)
        newentries.append(TransformedPDEntry(comp, scaled_energy, entry))
    return newentries

    def write_xdatcar(self, filepath="XDATCAR", groupby_type=True, overwrite=False, to_unit_cell=False):
        """
        Write Xdatcar file with unit cell and atomic positions to file ``filepath``.

        Args:
            filepath: Xdatcar filename. If None, a temporary file is created.
            groupby_type: If True, atoms are grouped by type. Note that this option
                may change the order of the atoms. This option is needed because
                there are post-processing tools (e.g. ovito) that do not work as expected
                if the atoms in the structure are not grouped by type.
            overwrite: raise RuntimeError, if False and filepath exists.
            to_unit_cell (bool): Whether to translate sites into the unit cell.

        Return:
            path to Xdatcar file.
        """
        if filepath is not None and os.path.exists(filepath) and not overwrite:
            raise RuntimeError("Cannot overwrite pre-existing file `%s`" % filepath)
        if filepath is None:
            import tempfile
            fd, filepath = tempfile.mkstemp(text=True, suffix="_XDATCAR")

        # int typat[natom], double znucl[npsp]
        # NB: typat is double in the HIST.nc file
        typat = self.reader.read_value("typat").astype(int)
        znucl = self.reader.read_value("znucl")
        ntypat = self.reader.read_dimvalue("ntypat")
        num_pseudos = self.reader.read_dimvalue("npsp")
        if num_pseudos != ntypat:
            raise NotImplementedError("Alchemical mixing is not supported, num_pseudos != ntypat")
        #print("znucl:", znucl, "\ntypat:", typat)

        symb2pos = OrderedDict()
        symbols_atom = []
        for iatom, itype in enumerate(typat):
            itype = itype - 1
            symbol = Element.from_Z(int(znucl[itype])).symbol
            if symbol not in symb2pos: symb2pos[symbol] = []
            symb2pos[symbol].append(iatom)
            symbols_atom.append(symbol)

        if not groupby_type:
            group_ids = np.arange(self.reader.natom)
        else:
            group_ids = []
            for pos_list in symb2pos.values():
                group_ids.extend(pos_list)
            group_ids = np.array(group_ids, dtype=np.int)

        comment = " %s\n" % self.initial_structure.formula
        with open(filepath, "wt") as fh:
            # comment line  + scaling factor set to 1.0
            fh.write(comment)
            fh.write("1.0\n")
            for vec in self.initial_structure.lattice.matrix:
                fh.write("%.12f %.12f %.12f\n" % (vec[0], vec[1], vec[2]))
            if not groupby_type:
                fh.write(" ".join(symbols_atom) + "\n")
                fh.write("1 " * len(symbols_atom) + "\n")
            else:
                fh.write(" ".join(symb2pos.keys()) + "\n")
                fh.write(" ".join(str(len(p)) for p in symb2pos.values()) + "\n")

            # Write atomic positions in reduced coordinates.
            xred_list = self.reader.read_value("xred")
            if to_unit_cell:
                xred_list = xred_list % 1

            for step in range(self.num_steps):
                fh.write("Direct configuration= %d\n" % (step + 1))
                frac_coords = xred_list[step, group_ids]
                for fs in frac_coords:
                    fh.write("%.12f %.12f %.12f\n" % (fs[0], fs[1], fs[2]))

        return filepath

 def test_print_periodic_table(self):
     Element.print_periodic_table()

 def validate(self, value):
     if not all(Element.is_valid_symbol(k) for k in value.keys()):
         self.error('Keys should be element symbols')
     if not all(isinstance(v, (float, int)) for v in value.values()):
         self.error('Values should be numbers')
     super(DictField, self).validate(value)

def periodic_table_heatmap(elemental_data, cbar_label="",
                           show_plot=False, cmap="YlOrRd", blank_color="grey",
                           value_format=None, max_row=9):
    """
    A static method that generates a heat map overlapped on a periodic table.

    Args:
         elemental_data (dict): A dictionary with the element as a key and a
            value assigned to it, e.g. surface energy and frequency, etc.
            Elements missing in the elemental_data will be grey by default
            in the final table elemental_data={"Fe": 4.2, "O": 5.0}.
         cbar_label (string): Label of the colorbar. Default is "".
         figure_name (string): Name of the plot (absolute path) being saved
            if not None.
         show_plot (bool): Whether to show the heatmap. Default is False.
         value_format (str): Formatting string to show values. If None, no value
            is shown. Example: "%.4f" shows float to four decimals.
         cmap (string): Color scheme of the heatmap. Default is 'coolwarm'.
         blank_color (string): Color assigned for the missing elements in
            elemental_data. Default is "grey".
         max_row (integer): Maximum number of rows of the periodic table to be
            shown. Default is 9, which means the periodic table heat map covers
            the first 9 rows of elements.
    """

    # Convert elemental data in the form of numpy array for plotting.
    max_val = max(elemental_data.values())
    min_val = min(elemental_data.values())
    max_row = min(max_row, 9)

    if max_row <= 0:
        raise ValueError("The input argument 'max_row' must be positive!")

    value_table = np.empty((max_row, 18)) * np.nan
    blank_value = min_val - 0.01

    for el in Element:
        if el.row > max_row: continue
        value = elemental_data.get(el.symbol, blank_value)
        value_table[el.row - 1, el.group - 1] = value

    # Initialize the plt object
    import matplotlib.pyplot as plt
    fig, ax = plt.subplots()
    plt.gcf().set_size_inches(12, 8)

    # We set nan type values to masked values (ie blank spaces)
    data_mask = np.ma.masked_invalid(value_table.tolist())
    heatmap = ax.pcolor(data_mask, cmap=cmap, edgecolors='w', linewidths=1,
                        vmin=min_val-0.001, vmax=max_val+0.001)
    cbar = fig.colorbar(heatmap)

    # Grey out missing elements in input data
    cbar.cmap.set_under(blank_color)
    cbar.set_label(cbar_label, rotation=270, labelpad=15)
    cbar.ax.tick_params(labelsize=14)

    # Refine and make the table look nice
    ax.axis('off')
    ax.invert_yaxis()

    # Label each block with corresponding element and value
    for i, row in enumerate(value_table):
        for j, el in enumerate(row):
            if not np.isnan(el):
                symbol = Element.from_row_and_group(i+1, j+1).symbol
                plt.text(j + 0.5, i + 0.25, symbol,
                         horizontalalignment='center',
                         verticalalignment='center', fontsize=14)
                if el != blank_value and value_format is not None:
                    plt.text(j + 0.5, i + 0.5, value_format % el,
                             horizontalalignment='center',
                             verticalalignment='center', fontsize=10)

    plt.tight_layout()

    if show_plot:
        plt.show()

    return plt

    def from_string(data, default_names=None):
        """
        Reads a Poscar from a string.

        The code will try its best to determine the elements in the POSCAR in
        the following order:
        1. If default_names are supplied and valid, it will use those. Usually,
        default names comes from an external source, such as a POTCAR in the
        same directory.
        2. If there are no valid default names but the input file is Vasp5-like
        and contains element symbols in the 6th line, the code will use that.
        3. Failing (2), the code will check if a symbol is provided at the end
        of each coordinate.

        If all else fails, the code will just assign the first n elements in
        increasing atomic number, where n is the number of species, to the
        Poscar. For example, H, He, Li, ....  This will ensure at least a
        unique element is assigned to each site and any analysis that does not
        require specific elemental properties should work fine.

        Args:
            data:
                string containing Poscar data.
            default_names:
                default symbols for the POSCAR file, usually coming from a
                POTCAR in the same directory.

        Returns:
            Poscar object.
        """

        chunks = re.split("^\s*$", data.strip(), flags=re.MULTILINE)

        #Parse positions
        lines = tuple(clean_lines(chunks[0].split("\n"), False))
        comment = lines[0]
        scale = float(lines[1])
        lattice = np.array([map(float, line.split())
                            for line in lines[2:5]])
        if scale < 0:
            # In vasp, a negative scale factor is treated as a volume. We need
            # to translate this to a proper lattice vector scaling.
            vol = abs(det(lattice))
            lattice *= (-scale / vol) ** (1 / 3)
        else:
            lattice *= scale

        vasp5_symbols = False
        try:
            natoms = map(int, lines[5].split())
            ipos = 6
        except ValueError:
            vasp5_symbols = True
            symbols = lines[5].split()
            natoms = map(int, lines[6].split())
            atomic_symbols = list()
            for i in xrange(len(natoms)):
                atomic_symbols.extend([symbols[i]] * natoms[i])
            ipos = 7

        postype = lines[ipos].split()[0]

        sdynamics = False
        # Selective dynamics
        if postype[0] in "sS":
            sdynamics = True
            ipos += 1
            postype = lines[ipos].split()[0]

        cart = postype[0] in "cCkK"
        nsites = sum(natoms)

        # If default_names is specified (usually coming from a POTCAR), use
        # them. This is in line with Vasp"s parsing order that the POTCAR
        # specified is the default used.
        if default_names:
            try:
                atomic_symbols = []
                for i in xrange(len(natoms)):
                    atomic_symbols.extend([default_names[i]] * natoms[i])
                vasp5_symbols = True
            except IndexError:
                pass

        if not vasp5_symbols:
            ind = 3 if not sdynamics else 6
            try:
                #check if names are appended at the end of the coordinates.
                atomic_symbols = [l.split()[ind]
                                  for l in lines[ipos + 1:ipos + 1 + nsites]]
                #Ensure symbols are valid elements
                if not all([Element.is_valid_symbol(sym)
                            for sym in atomic_symbols]):
                    raise ValueError("Non-valid symbols detected.")
                vasp5_symbols = True
            except (ValueError, IndexError):
                #Defaulting to false names.
                atomic_symbols = []
                for i in xrange(len(natoms)):
                    sym = Element.from_Z(i + 1).symbol
#.........这里部分代码省略.........

    def _parse(self, filename):

        start_patt = re.compile(" \(Enter \S+l101\.exe\)")
        route_patt = re.compile(" #[pPnNtT]*.*")
        charge_mul_patt = re.compile("Charge\s+=\s*([-\\d]+)\s+" "Multiplicity\s+=\s*(\d+)")
        num_basis_func_patt = re.compile("([0-9]+)\s+basis functions")
        pcm_patt = re.compile("Polarizable Continuum Model")
        stat_type_patt = re.compile("imaginary frequencies")
        scf_patt = re.compile("E\(.*\)\s*=\s*([-\.\d]+)\s+")
        mp2_patt = re.compile("EUMP2\s*=\s*(.*)")
        oniom_patt = re.compile("ONIOM:\s+extrapolated energy\s*=\s*(.*)")
        termination_patt = re.compile("(Normal|Error) termination of Gaussian")
        std_orientation_patt = re.compile("Standard orientation")
        end_patt = re.compile("--+")
        orbital_patt = re.compile("Alpha\s*\S+\s*eigenvalues --(.*)")
        thermo_patt = re.compile("(Zero-point|Thermal) correction(.*)=" "\s+([\d\.-]+)")

        self.properly_terminated = False
        self.is_pcm = False
        self.stationary_type = "Minimum"
        self.structures = []
        self.corrections = {}
        self.energies = []

        coord_txt = []
        read_coord = 0
        orbitals_txt = []
        parse_stage = 0
        num_basis_found = False
        terminated = False

        with zopen(filename, "r") as f:
            for line in f:
                if parse_stage == 0:
                    if start_patt.search(line):
                        parse_stage = 1
                    elif route_patt.search(line):
                        self.route = {}
                        for tok in line.split():
                            sub_tok = tok.strip().split("=")
                            key = sub_tok[0].upper()
                            self.route[key] = sub_tok[1].upper() if len(sub_tok) > 1 else ""
                            m = re.match("(\w+)/([^/]+)", key)
                            if m:
                                self.functional = m.group(1)
                                self.basis_set = m.group(2)
                elif parse_stage == 1:
                    if charge_mul_patt.search(line):
                        m = charge_mul_patt.search(line)
                        self.charge = int(m.group(1))
                        self.spin_mult = int(m.group(2))
                        parse_stage = 2
                elif parse_stage == 2:

                    if self.is_pcm:
                        self._check_pcm(line)

                    if "FREQ" in self.route and thermo_patt.search(line):
                        m = thermo_patt.search(line)
                        if m.group(1) == "Zero-point":
                            self.corrections["Zero-point"] = float(m.group(3))
                        else:
                            key = m.group(2).strip(" to ")
                            self.corrections[key] = float(m.group(3))

                    if read_coord:
                        if not end_patt.search(line):
                            coord_txt.append(line)
                        else:
                            read_coord = (read_coord + 1) % 4
                            if not read_coord:
                                sp = []
                                coords = []
                                for l in coord_txt[2:]:
                                    toks = l.split()
                                    sp.append(Element.from_Z(int(toks[1])))
                                    coords.append(map(float, toks[3:6]))
                                self.structures.append(Molecule(sp, coords))
                    elif termination_patt.search(line):
                        m = termination_patt.search(line)
                        if m.group(1) == "Normal":
                            self.properly_terminated = True
                        terminated = True
                    elif (not num_basis_found) and num_basis_func_patt.search(line):
                        m = num_basis_func_patt.search(line)
                        self.num_basis_func = int(m.group(1))
                        num_basis_found = True
                    elif (not self.is_pcm) and pcm_patt.search(line):
                        self.is_pcm = True
                        self.pcm = {}
                    elif "FREQ" in self.route and "OPT" in self.route and stat_type_patt.search(line):
                        self.stationary_type = "Saddle"
                    elif mp2_patt.search(line):
                        m = mp2_patt.search(line)
                        self.energies.append(float(m.group(1).replace("D", "E")))
                    elif oniom_patt.search(line):
                        m = oniom_patt.matcher(line)
                        self.energies.append(float(m.group(1)))
                    elif scf_patt.search(line):
                        m = scf_patt.search(line)
#.........这里部分代码省略.........

def get_excluded_list():
    
    filename = "excluded_compounds.p"
    if os.path.exists(filename):
        with open(filename) as f:
            return pickle.load(f)
    
    from ga_optimization_ternary.fitness_evaluators import FitnessEvaluator, eval_fitness_simple
    all_AB = FitnessEvaluator(eval_fitness_simple, 10)._reverse_dict.keys()
    print 'generating exclusion ranks...'
    exclusions = []
    
    for a in all_AB:
        for b in all_AB:
            for x in range(7):
                
                #nelectrons must be even
                ne_a = Element.from_Z(a).Z
                ne_b = Element.from_Z(b).Z
                ne_x = None
                if x == 0:
                    ne_x = Element("O").Z * 3
                elif x == 1:
                    ne_x = Element("O").Z * 2 + Element("N").Z
                elif x == 2:
                    ne_x = Element("O").Z + Element("N").Z * 2
                elif x == 3:
                    ne_x = Element("N").Z
                elif x == 4:
                    ne_x = Element("O").Z * 2 + Element("F").Z
                elif x == 5:
                    ne_x = Element("O").Z + Element("F").Z + Element("N").Z
                elif x == 6:
                    ne_x = Element("O").Z * 2 + Element("S").Z
                
                #modify the score based on charge-balance
                even_found = False
                el_a = Element.from_Z(a)
                el_b = Element.from_Z(b)
                val_x = 0
                if x == 0:
                    val_x = -2 * 3
                elif x == 1:
                    val_x = (-2 * 2) + (-3 * 1)
                elif x == 2:
                    val_x = (-2 * 1) + (-3 * 2)
                elif x == 3:
                    val_x = -3 * 3
                elif x == 4:
                    val_x = (-2 * 2) + (-1 * 1)
                elif x == 5:
                    val_x = (-2 * 1) + (-1 * 1) + (-3 * 1)
                elif x == 6:
                    val_x = (-2 * 2) + (-2 * 1)
                
                for a_oxi in el_a.oxidation_states:
                    for b_oxi in el_b.oxidation_states:
                        if (ne_a + ne_b + ne_x) % 2 == 0 and (a_oxi + b_oxi + val_x) == 0:
                            even_found = True
                
                if not even_found:
                    exclusions.append((a, b, x))
                
    with open(filename, "wb") as f:
        pickle.dump(exclusions, f)
    return exclusions

def get_ranked_list_goldschmidt_halffill():
    #TODO: get the oxidation state right!!!
    
    filename = "goldschmidt_rank_halffill.p"
    if os.path.exists(filename):
        with open(filename) as f:
            return pickle.load(f)
    
    from ga_optimization_ternary.fitness_evaluators import FitnessEvaluator, eval_fitness_simple
    all_AB = FitnessEvaluator(eval_fitness_simple, 10)._reverse_dict.keys()
    print 'generating goldschmidt ranks...'
    cand_score = {}  # dictionary of cand_tuple:score. a high score is BAD
    
    for a in all_AB:
        for b in all_AB:
            for x in range(7):
                r_a = Element.from_Z(a).average_ionic_radius  # TODO: get the correct oxidation state!
                r_b = Element.from_Z(b).average_ionic_radius
                r_x = None
                if x == 0:
                    r_x = Element("O").ionic_radii[-2]
                elif x == 1:
                    r_x = Element("O").ionic_radii[-2] * 2/3 + Element("N").ionic_radii[-3] * 1/3
                elif x == 2:
                    r_x = Element("O").ionic_radii[-2] * 1/3 + Element("N").ionic_radii[-3] * 2/3
                elif x == 3:
                    r_x = Element("N").ionic_radii[-3]
                elif x == 4:
                    r_x = Element("O").ionic_radii[-2] * 2/3 + Element("F").ionic_radii[-1] * 1/3
                elif x == 5:
                    r_x = Element("O").ionic_radii[-2] * 1/3 + Element("F").ionic_radii[-1] * 1/3 + Element("N").ionic_radii[-3] * 1/3
                elif x == 6:
                    r_x = Element("O").ionic_radii[-2] * 2/3 + Element("S").ionic_radii[-2] * 1/3
                
                goldschmidt = (r_a + r_x)/(math.sqrt(2) *(r_b+r_x))
                score = abs(goldschmidt - 1)  # a high score is bad, like golf
                
                #nelectrons must be even
                ne_a = Element.from_Z(a).Z
                ne_b = Element.from_Z(b).Z
                ne_x = None
                if x == 0:
                    ne_x = Element("O").Z * 3
                elif x == 1:
                    ne_x = Element("O").Z * 2 + Element("N").Z
                elif x == 2:
                    ne_x = Element("O").Z + Element("N").Z * 2
                elif x == 3:
                    ne_x = Element("N").Z
                elif x == 4:
                    ne_x = Element("O").Z * 2 + Element("F").Z
                elif x == 5:
                    ne_x = Element("O").Z + Element("F").Z + Element("N").Z
                elif x == 6:
                    ne_x = Element("O").Z * 2 + Element("S").Z
                
                #modify the score based on charge-balance
                even_found = False
                el_a = Element.from_Z(a)
                el_b = Element.from_Z(b)
                val_x = 0
                if x == 0:
                    val_x = -2 * 3
                elif x == 1:
                    val_x = (-2 * 2) + (-3 * 1)
                elif x == 2:
                    val_x = (-2 * 1) + (-3 * 2)
                elif x == 3:
                    val_x = -3 * 3
                elif x == 4:
                    val_x = (-2 * 2) + (-1 * 1)
                elif x == 5:
                    val_x = (-2 * 1) + (-1 * 1) + (-3 * 1)
                elif x == 6:
                    val_x = (-2 * 2) + (-2 * 1)
                
                for a_oxi in el_a.oxidation_states:
                    for b_oxi in el_b.oxidation_states: