4075885785

Build Type

push

github

Committed by Shyue Ping Ong

Commit Message

Merge branch 'master' of github.com:materialsproject/pymatgen

Run Details

96 of 96 new or added lines in 27 files covered. (100.0%)

81013 of 102710 relevant lines covered (78.88%)

0.79 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

97.85

/pymatgen/analysis/tests/test_surface_analysis.py

from __future__ import annotations

import json
import os
import unittest
import warnings

from pytest import approx
from sympy import Number, Symbol

from pymatgen.analysis.surface_analysis import (
    NanoscaleStability,
    SlabEntry,
    SurfaceEnergyPlotter,
    WorkFunctionAnalyzer,
)
from pymatgen.entries.computed_entries import ComputedStructureEntry
from pymatgen.util.testing import PymatgenTest

__author__ = "Richard Tran"
__copyright__ = "Copyright 2012, The Materials Project"
__version__ = "0.1"
__maintainer__ = "Richard Tran"
__email__ = "rit001@eng.ucsd.edu"
__date__ = "Aug 24, 2017"


def get_path(path_str):
    path = os.path.join(PymatgenTest.TEST_FILES_DIR, "surface_tests", path_str)
    return path


class SlabEntryTest(PymatgenTest):
    def setUp(self):
        with warnings.catch_warnings():
            warnings.simplefilter("ignore")

        with open(os.path.join(get_path(""), "ucell_entries.txt")) as ucell_entries:
            ucell_entries = json.loads(ucell_entries.read())
        self.ucell_entries = ucell_entries

        # Load objects for O adsorption tests
        self.metals_O_entry_dict = load_O_adsorption()

        # Load objects for Cu test
        self.Cu_entry_dict = get_entry_dict(os.path.join(get_path(""), "Cu_entries.txt"))
        assert len(self.Cu_entry_dict) == 13
        self.Cu_ucell_entry = ComputedStructureEntry.from_dict(self.ucell_entries["Cu"])

        # Load dummy MgO slab entries
        self.MgO_ucell_entry = ComputedStructureEntry.from_dict(self.ucell_entries["MgO"])
        self.Mg_ucell_entry = ComputedStructureEntry.from_dict(self.ucell_entries["Mg"])
        self.MgO_slab_entry_dict = get_entry_dict(os.path.join(get_path(""), "MgO_slab_entries.txt"))

    def test_properties(self):
        # Test cases for getting adsorption related quantities for a 1/4
        # monolalyer adsorption of O on the low MMI surfaces of Pt, Ni and Rh

        for el, val in self.metals_O_entry_dict.items():
            el_ucell = ComputedStructureEntry.from_dict(self.ucell_entries[el])
            for hkl in val:
                for clean in self.metals_O_entry_dict[el][hkl]:
                    for ads in self.metals_O_entry_dict[el][hkl][clean]:
                        ml = ads.get_unit_primitive_area
                        assert ml == approx(4, abs=1e-2)
                        assert ads.get_monolayer == approx(1 / 4, abs=1e-2)
                        Nads = ads.Nads_in_slab
                        assert Nads == 1
                        assert ads.Nsurfs_ads_in_slab == 1

                        # Determine the correct binding energy
                        with open(os.path.join(get_path(""), "isolated_O_entry.txt")) as isolated_O_entry:
                            isolated_O_entry = json.loads(isolated_O_entry.read())
                        O = ComputedStructureEntry.from_dict(isolated_O_entry)
                        gbind = (ads.energy - ml * clean.energy) / Nads - O.energy_per_atom
                        assert gbind == ads.gibbs_binding_energy()
                        # Determine the correction Gibbs adsorption energy
                        eads = Nads * gbind
                        assert eads == ads.gibbs_binding_energy(eads=True)
                        se = ads.surface_energy(el_ucell)
                        assert se.as_coefficients_dict()[Symbol("delu_O")] == approx(
                            (-1 / 2) * ads.surface_area ** (-1)
                        )

    def test_create_slab_label(self):
        for el, val in self.metals_O_entry_dict.items():
            for hkl in val:
                # Test WulffShape for adsorbed surfaces
                for clean in self.metals_O_entry_dict[el][hkl]:
                    label = clean.create_slab_label
                    comp = str(clean.composition.reduced_composition)
                    assert f"{hkl} {comp}" == label

                    for ads in self.metals_O_entry_dict[el][hkl][clean]:
                        label = ads.create_slab_label
                        assert label == f"{hkl} {comp}+O, 0.250 ML"

    def test_surface_energy(self):
        # For a non-stoichiometric case, the chemical potentials do not
        # cancel out, they serve as a reservoir for any missing atoms
        for slab_entry in self.MgO_slab_entry_dict[(1, 1, 1)]:
            se = slab_entry.surface_energy(self.MgO_ucell_entry, ref_entries=[self.Mg_ucell_entry])
            assert tuple(se.as_coefficients_dict()) == (Number(1), Symbol("delu_Mg"))

        # For the case of a clean, stoichiometric slab, the surface energy
        # should be constant (i.e. surface energy is a constant).
        all_se = []
        ECu = self.Cu_ucell_entry.energy_per_atom
        for val in self.Cu_entry_dict.values():
            slab_entry = list(val)[0]
            se = slab_entry.surface_energy(self.Cu_ucell_entry)
            all_se.append(se)
            # Manually calculate surface energy
            manual_se = (slab_entry.energy - ECu * len(slab_entry.structure)) / (2 * slab_entry.surface_area)
            self.assertArrayAlmostEqual(float(se), manual_se, 10)

        # The (111) facet should be the most stable
        clean111_entry = list(self.Cu_entry_dict[(1, 1, 1)])[0]
        se_Cu111 = clean111_entry.surface_energy(self.Cu_ucell_entry)
        assert min(all_se) == se_Cu111

    def test_cleaned_up_slab(self):
        # The cleaned up slab should have the same reduced formula as a clean slab
        for el, val in self.metals_O_entry_dict.items():
            for hkl in val:
                for clean in self.metals_O_entry_dict[el][hkl]:
                    for ads in self.metals_O_entry_dict[el][hkl][clean]:
                        s = ads.cleaned_up_slab
                        assert s.composition.reduced_composition == clean.composition.reduced_composition


class SurfaceEnergyPlotterTest(PymatgenTest):
    def setUp(self):
        entry_dict = get_entry_dict(os.path.join(get_path(""), "Cu_entries.txt"))
        self.Cu_entry_dict = entry_dict
        with open(os.path.join(get_path(""), "ucell_entries.txt")) as ucell_entries:
            ucell_entries = json.loads(ucell_entries.read())

        self.Cu_ucell_entry = ComputedStructureEntry.from_dict(ucell_entries["Cu"])
        self.Cu_analyzer = SurfaceEnergyPlotter(entry_dict, self.Cu_ucell_entry)

        self.metals_O_entry_dict = load_O_adsorption()
        ucell_entry = ComputedStructureEntry.from_dict(ucell_entries["Pt"])
        self.Pt_analyzer = SurfaceEnergyPlotter(self.metals_O_entry_dict["Pt"], ucell_entry)
        ucell_entry = ComputedStructureEntry.from_dict(ucell_entries["Ni"])
        self.Ni_analyzer = SurfaceEnergyPlotter(self.metals_O_entry_dict["Ni"], ucell_entry)
        ucell_entry = ComputedStructureEntry.from_dict(ucell_entries["Rh"])
        self.Rh_analyzer = SurfaceEnergyPlotter(self.metals_O_entry_dict["Rh"], ucell_entry)
        self.Oads_analyzer_dict = {
            "Pt": self.Pt_analyzer,
            "Ni": self.Ni_analyzer,
            "Rh": self.Rh_analyzer,
        }

    def test_get_stable_entry_at_u(self):
        for plotter in self.Oads_analyzer_dict.values():
            for hkl in plotter.all_slab_entries:
                # Test that the surface energy is clean for specific range of chempot
                entry1, gamma1 = plotter.get_stable_entry_at_u(hkl, delu_dict={Symbol("delu_O"): -7})
                entry2, gamma2 = plotter.get_stable_entry_at_u(hkl, delu_dict={Symbol("delu_O"): -6})
                assert gamma1 == gamma2
                assert entry1.label == entry2.label

                # Now test that for a high chempot, adsorption
                # occurs and gamma is not equal to clean gamma
                entry3, gamma3 = plotter.get_stable_entry_at_u(hkl, delu_dict={Symbol("delu_O"): -1})
                assert entry3.label != entry2.label
                assert gamma3 != gamma2

                # For any chempot greater than -6, surface energy should vary
                # but the configuration should remain the same
                entry4, gamma4 = plotter.get_stable_entry_at_u(hkl, delu_dict={Symbol("delu_O"): 0})
                assert entry3.label == entry4.label
                assert gamma3 != gamma4

    def test_wulff_from_chempot(self):
        # Test if it generates a Wulff shape, test if
        # all the facets for Cu wulff shape are inside.
        Cu_wulff = self.Cu_analyzer.wulff_from_chempot()
        area_frac_dict = Cu_wulff.area_fraction_dict
        facets_hkl = [
            (1, 1, 1),
            (3, 3, 1),
            (3, 1, 0),
            (1, 0, 0),
            (3, 1, 1),
            (2, 1, 0),
            (2, 2, 1),
        ]
        for hkl in area_frac_dict:
            if hkl in facets_hkl:
                assert area_frac_dict[hkl] != 0
            else:
                assert area_frac_dict[hkl] == 0

        for analyzer in self.Oads_analyzer_dict.values():
            # Test WulffShape for adsorbed surfaces
            # chempot = analyzer.max_adsorption_chempot_range(0)
            wulff = analyzer.wulff_from_chempot(delu_default=-6)
            wulff.weighted_surface_energy

        # Test if a different Wulff shape is generated
        # for Ni when adsorption comes into play
        wulff_neg7 = self.Oads_analyzer_dict["Ni"].wulff_from_chempot(delu_default=-7)
        wulff_neg6 = self.Oads_analyzer_dict["Ni"].wulff_from_chempot(delu_default=-6)
        assert wulff_neg7.weighted_surface_energy == wulff_neg6.weighted_surface_energy
        wulff_neg55 = self.Oads_analyzer_dict["Ni"].wulff_from_chempot(delu_default=-5.5)
        assert wulff_neg55.weighted_surface_energy != wulff_neg6.weighted_surface_energy
        wulff_neg525 = self.Oads_analyzer_dict["Ni"].wulff_from_chempot(delu_default=-5.25)
        assert wulff_neg55.weighted_surface_energy != wulff_neg525.weighted_surface_energy

    def test_color_palette_dict(self):
        for el, val in self.metals_O_entry_dict.items():
            analyzer = self.Oads_analyzer_dict[el]
            color_dict = analyzer.color_palette_dict()
            for hkl in val:
                for clean in self.metals_O_entry_dict[el][hkl]:
                    _ = color_dict[clean]
                    for ads in self.metals_O_entry_dict[el][hkl][clean]:
                        _ = color_dict[ads]

    def test_get_surface_equilibrium(self):
        # For clean stoichiometric system, the two equations should
        # be parallel because the surface energy is a constant. Then
        # get_surface_equilibrium should return None
        clean111_entry = list(self.Cu_entry_dict[(1, 1, 1)])[0]
        clean100_entry = list(self.Cu_entry_dict[(1, 0, 0)])[0]
        soln = self.Cu_analyzer.get_surface_equilibrium([clean111_entry, clean100_entry])
        assert not soln

        # For adsorbed system, we should find one intercept
        Pt_entries = self.metals_O_entry_dict["Pt"]
        clean = list(Pt_entries[(1, 1, 1)])[0]
        ads = Pt_entries[(1, 1, 1)][clean][0]
        Pt_analyzer = self.Oads_analyzer_dict["Pt"]
        soln = Pt_analyzer.get_surface_equilibrium([clean, ads])

        assert list(soln.values())[0] != list(soln.values())[1]

        # Check if the number of parameters for adsorption are correct
        assert (Symbol("delu_O"), Symbol("gamma")) == tuple(soln)
        # Adsorbed systems have a b2=(-1*Nads) / (Nsurfs * Aads)
        se = ads.surface_energy(Pt_analyzer.ucell_entry, Pt_analyzer.ref_entries)
        assert se.as_coefficients_dict()[Symbol("delu_O")] == approx(-1 / (2 * ads.surface_area))

    def test_stable_u_range_dict(self):
        analyzer = list(self.Oads_analyzer_dict.values())[-1]

        stable_u_range = analyzer.stable_u_range_dict([-1, 0], Symbol("delu_O"), no_doped=False)
        all_u = []
        for entry in stable_u_range:
            all_u.extend(stable_u_range[entry])
        assert len(all_u) > 1

    def test_entry_dict_from_list(self):
        # Plug in a list of entries to see if it works
        all_Pt_slab_entries = []
        Pt_entries = self.Pt_analyzer.all_slab_entries
        for hkl in Pt_entries:
            for clean in Pt_entries[hkl]:
                all_Pt_slab_entries.append(clean)
                all_Pt_slab_entries.extend(Pt_entries[hkl][clean])

        a = SurfaceEnergyPlotter(all_Pt_slab_entries, self.Pt_analyzer.ucell_entry)
        assert type(a).__name__ == "SurfaceEnergyPlotter"

    # def test_monolayer_vs_BE(self):
    #     for el in self.Oads_analyzer_dict:
    #         # Test WulffShape for adsorbed surfaces
    #         analyzer = self.Oads_analyzer_dict[el]
    #         plt = analyzer.monolayer_vs_BE()
    #
    # def test_area_frac_vs_chempot_plot(self):
    #
    #     for el in self.Oads_analyzer_dict:
    #         # Test WulffShape for adsorbed surfaces
    #         analyzer = self.Oads_analyzer_dict[el]
    #         plt = analyzer.area_frac_vs_chempot_plot(x_is_u_ads=True)
    #
    # def test_chempot_vs_gamma_clean(self):
    #
    #     plt = self.Cu_analyzer.chempot_vs_gamma_clean()
    #     for el in self.Oads_analyzer_dict:
    #         # Test WulffShape for adsorbed surfaces
    #         analyzer = self.Oads_analyzer_dict[el]
    #         plt = analyzer.chempot_vs_gamma_clean(x_is_u_ads=True)
    #
    # def test_chempot_vs_gamma_facet(self):
    #
    #     for el, val in self.metals_O_entry_dict.items():
    #         for hkl in val:
    #             # Test WulffShape for adsorbed surfaces
    #             analyzer = self.Oads_analyzer_dict[el]
    #             plt = analyzer.chempot_vs_gamma_facet(hkl)
    # def test_surface_chempot_range_map(self):
    #
    #     for el, val in self.metals_O_entry_dict.items():
    #         for hkl in val:
    #             # Test WulffShape for adsorbed surfaces
    #             analyzer = self.Oads_analyzer_dict[el]
    #             plt = analyzer.chempot_vs_gamma_facet(hkl)


class WorkfunctionAnalyzerTest(PymatgenTest):
    def setUp(self):
        self.kwargs = {
            "poscar_filename": get_path("CONTCAR.relax1.gz"),
            "locpot_filename": get_path("LOCPOT.gz"),
            "outcar_filename": get_path("OUTCAR.relax1.gz"),
        }
        self.wf_analyzer = WorkFunctionAnalyzer.from_files(**self.kwargs)

    def test_shift(self):
        wf_analyzer_shift = WorkFunctionAnalyzer.from_files(shift=-0.25, blength=3.7, **self.kwargs)
        assert self.wf_analyzer.ave_bulk_p == approx(wf_analyzer_shift.ave_bulk_p, abs=1e-1)

    def test_is_converged(self):
        assert self.wf_analyzer.is_converged()


class NanoscaleStabilityTest(PymatgenTest):
    def setUp(self):
        # Load all entries
        La_hcp_entry_dict = get_entry_dict(os.path.join(get_path(""), "La_hcp_entries.txt"))
        La_fcc_entry_dict = get_entry_dict(os.path.join(get_path(""), "La_fcc_entries.txt"))
        with open(os.path.join(get_path(""), "ucell_entries.txt")) as ucell_entries:
            ucell_entries = json.loads(ucell_entries.read())
        La_hcp_ucell_entry = ComputedStructureEntry.from_dict(ucell_entries["La_hcp"])
        La_fcc_ucell_entry = ComputedStructureEntry.from_dict(ucell_entries["La_fcc"])

        # Set up the NanoscaleStabilityClass
        self.La_hcp_analyzer = SurfaceEnergyPlotter(La_hcp_entry_dict, La_hcp_ucell_entry)
        self.La_fcc_analyzer = SurfaceEnergyPlotter(La_fcc_entry_dict, La_fcc_ucell_entry)
        self.nanoscale_stability = NanoscaleStability([self.La_fcc_analyzer, self.La_hcp_analyzer])

    def test_stability_at_r(self):
        # Check that we have a different polymorph that is
        # stable below or above the equilibrium particle size
        r = self.nanoscale_stability.solve_equilibrium_point(self.La_hcp_analyzer, self.La_fcc_analyzer) * 10

        # hcp phase of La particle should be the stable
        # polymorph above the equilibrium radius
        hcp_wulff = self.La_hcp_analyzer.wulff_from_chempot()
        bulk = self.La_hcp_analyzer.ucell_entry
        ghcp, rhcp = self.nanoscale_stability.wulff_gform_and_r(hcp_wulff, bulk, r + 10, from_sphere_area=True)
        fcc_wulff = self.La_fcc_analyzer.wulff_from_chempot()
        bulk = self.La_fcc_analyzer.ucell_entry
        gfcc, rfcc = self.nanoscale_stability.wulff_gform_and_r(fcc_wulff, bulk, r + 10, from_sphere_area=True)
        assert gfcc > ghcp

        # fcc phase of La particle should be the stable
        # polymorph below the equilibrium radius
        hcp_wulff = self.La_hcp_analyzer.wulff_from_chempot()
        bulk = self.La_hcp_analyzer.ucell_entry
        ghcp, rhcp = self.nanoscale_stability.wulff_gform_and_r(hcp_wulff, bulk, r - 10, from_sphere_area=True)
        fcc_wulff = self.La_fcc_analyzer.wulff_from_chempot()
        bulk = self.La_fcc_analyzer.ucell_entry
        gfcc, rfcc = self.nanoscale_stability.wulff_gform_and_r(fcc_wulff, bulk, r - 10, from_sphere_area=True)
        assert gfcc < ghcp

    def test_scaled_wulff(self):
        # Ensure for a given radius, the effective radius
        # of the Wulff shape is the same (correctly scaled)
        hcp_wulff = self.La_hcp_analyzer.wulff_from_chempot()
        fcc_wulff = self.La_fcc_analyzer.wulff_from_chempot()
        w1 = self.nanoscale_stability.scaled_wulff(hcp_wulff, 10)
        w2 = self.nanoscale_stability.scaled_wulff(fcc_wulff, 10)
        assert w1.effective_radius == approx(w2.effective_radius)
        assert w1.effective_radius == approx(10)
        assert 10 == approx(w2.effective_radius)


def get_entry_dict(filename):
    # helper to generate an entry_dict

    entry_dict = {}
    with open(filename) as entries:
        entries = json.loads(entries.read())
    for k in entries:
        n = k[25:]
        miller_index = []
        for i, s in enumerate(n):
            if s == "_":
                break
            if s == "-":
                continue
            t = int(s)
            if n[i - 1] == "-":
                t *= -1
            miller_index.append(t)
        hkl = tuple(miller_index)
        if hkl not in entry_dict:
            entry_dict[hkl] = {}
        entry = ComputedStructureEntry.from_dict(entries[k])
        entry_dict[hkl][SlabEntry(entry.structure, entry.energy, hkl, label=k)] = []

    return entry_dict


def load_O_adsorption():
    # Loads the dictionary for clean and O adsorbed Rh, Pt, and Ni entries

    # Load the adsorbate as an entry
    with open(os.path.join(get_path(""), "isolated_O_entry.txt")) as isolated_O_entry:
        isolated_O_entry = json.loads(isolated_O_entry.read())
    O = ComputedStructureEntry.from_dict(isolated_O_entry)

    # entry_dict for the adsorption case, O adsorption on Ni, Rh and Pt
    metals_O_entry_dict = {
        "Ni": {(1, 1, 1): {}, (1, 0, 0): {}},
        "Pt": {(1, 1, 1): {}},
        "Rh": {(1, 0, 0): {}},
    }

    with open(os.path.join(get_path(""), "csentries_slabs.json")) as entries:
        entries = json.loads(entries.read())
    for k in entries:
        entry = ComputedStructureEntry.from_dict(entries[k])
        for el in metals_O_entry_dict:  # pylint: disable=C0206
            if el in k:
                if "111" in k:
                    clean = SlabEntry(entry.structure, entry.energy, (1, 1, 1), label=k + "_clean")
                    metals_O_entry_dict[el][(1, 1, 1)][clean] = []
                if "110" in k:
                    clean = SlabEntry(entry.structure, entry.energy, (1, 1, 0), label=k + "_clean")
                    metals_O_entry_dict[el][(1, 1, 0)][clean] = []
                if "100" in k:
                    clean = SlabEntry(entry.structure, entry.energy, (1, 0, 0), label=k + "_clean")
                    metals_O_entry_dict[el][(1, 0, 0)][clean] = []

    with open(os.path.join(get_path(""), "csentries_o_ads.json")) as entries:
        entries = json.loads(entries.read())
    for k in entries:
        entry = ComputedStructureEntry.from_dict(entries[k])
        for el, val in metals_O_entry_dict.items():
            if el in k:
                if "111" in k:
                    clean = list(val[(1, 1, 1)])[0]
                    ads = SlabEntry(
                        entry.structure,
                        entry.energy,
                        (1, 1, 1),
                        label=k + "_O",
                        adsorbates=[O],
                        clean_entry=clean,
                    )
                    metals_O_entry_dict[el][(1, 1, 1)][clean] = [ads]
                if "110" in k:
                    clean = list(val[(1, 1, 0)])[0]
                    ads = SlabEntry(
                        entry.structure,
                        entry.energy,
                        (1, 1, 0),
                        label=k + "_O",
                        adsorbates=[O],
                        clean_entry=clean,
                    )
                    metals_O_entry_dict[el][(1, 1, 0)][clean] = [ads]
                if "100" in k:
                    clean = list(val[(1, 0, 0)])[0]
                    ads = SlabEntry(
                        entry.structure,
                        entry.energy,
                        (1, 0, 0),
                        label=k + "_O",
                        adsorbates=[O],
                        clean_entry=clean,
                    )
                    metals_O_entry_dict[el][(1, 0, 0)][clean] = [ads]

    return metals_O_entry_dict


if __name__ == "__main__":
    unittest.main()

1	from __future__ import annotations	1✔
2
3	import json	1✔
4	import os	1✔
5	import unittest	1✔
6	import warnings	1✔
7
8	from pytest import approx	1✔
9	from sympy import Number, Symbol	1✔
10
11	from pymatgen.analysis.surface_analysis import (	1✔
12	NanoscaleStability,
13	SlabEntry,
14	SurfaceEnergyPlotter,
15	WorkFunctionAnalyzer,
16	)
17	from pymatgen.entries.computed_entries import ComputedStructureEntry	1✔
18	from pymatgen.util.testing import PymatgenTest	1✔
19
20	__author__ = "Richard Tran"	1✔
21	__copyright__ = "Copyright 2012, The Materials Project"	1✔
22	__version__ = "0.1"	1✔
23	__maintainer__ = "Richard Tran"	1✔
24	__email__ = "rit001@eng.ucsd.edu"	1✔
25	__date__ = "Aug 24, 2017"	1✔
26
27
28	def get_path(path_str):	1✔
29	path = os.path.join(PymatgenTest.TEST_FILES_DIR, "surface_tests", path_str)	1✔
30	return path	1✔
31
32
33	class SlabEntryTest(PymatgenTest):	1✔
34	def setUp(self):	1✔
35	with warnings.catch_warnings():	1✔
36	warnings.simplefilter("ignore")	1✔
37
38	with open(os.path.join(get_path(""), "ucell_entries.txt")) as ucell_entries:	1✔
39	ucell_entries = json.loads(ucell_entries.read())	1✔
40	self.ucell_entries = ucell_entries	1✔
41
42	# Load objects for O adsorption tests
43	self.metals_O_entry_dict = load_O_adsorption()	1✔
44
45	# Load objects for Cu test
46	self.Cu_entry_dict = get_entry_dict(os.path.join(get_path(""), "Cu_entries.txt"))	1✔
47	assert len(self.Cu_entry_dict) == 13	1✔
48	self.Cu_ucell_entry = ComputedStructureEntry.from_dict(self.ucell_entries["Cu"])	1✔
49
50	# Load dummy MgO slab entries
51	self.MgO_ucell_entry = ComputedStructureEntry.from_dict(self.ucell_entries["MgO"])	1✔
52	self.Mg_ucell_entry = ComputedStructureEntry.from_dict(self.ucell_entries["Mg"])	1✔
53	self.MgO_slab_entry_dict = get_entry_dict(os.path.join(get_path(""), "MgO_slab_entries.txt"))	1✔
54
55	def test_properties(self):	1✔
56	# Test cases for getting adsorption related quantities for a 1/4
57	# monolalyer adsorption of O on the low MMI surfaces of Pt, Ni and Rh
58
59	for el, val in self.metals_O_entry_dict.items():	1✔
60	el_ucell = ComputedStructureEntry.from_dict(self.ucell_entries[el])	1✔
61	for hkl in val:	1✔
62	for clean in self.metals_O_entry_dict[el][hkl]:	1✔
63	for ads in self.metals_O_entry_dict[el][hkl][clean]:	1✔
64	ml = ads.get_unit_primitive_area	1✔
65	assert ml == approx(4, abs=1e-2)	1✔
66	assert ads.get_monolayer == approx(1 / 4, abs=1e-2)	1✔
67	Nads = ads.Nads_in_slab	1✔
68	assert Nads == 1	1✔
69	assert ads.Nsurfs_ads_in_slab == 1	1✔
70
71	# Determine the correct binding energy
72	with open(os.path.join(get_path(""), "isolated_O_entry.txt")) as isolated_O_entry:	1✔
73	isolated_O_entry = json.loads(isolated_O_entry.read())	1✔
74	O = ComputedStructureEntry.from_dict(isolated_O_entry)	1✔
75	gbind = (ads.energy - ml * clean.energy) / Nads - O.energy_per_atom	1✔
76	assert gbind == ads.gibbs_binding_energy()	1✔
77	# Determine the correction Gibbs adsorption energy
78	eads = Nads * gbind	1✔
79	assert eads == ads.gibbs_binding_energy(eads=True)	1✔
80	se = ads.surface_energy(el_ucell)	1✔
81	assert se.as_coefficients_dict()[Symbol("delu_O")] == approx(	1✔
82	(-1 / 2) * ads.surface_area ** (-1)
83	)
84
85	def test_create_slab_label(self):	1✔
86	for el, val in self.metals_O_entry_dict.items():	1✔
87	for hkl in val:	1✔
88	# Test WulffShape for adsorbed surfaces
89	for clean in self.metals_O_entry_dict[el][hkl]:	1✔
90	label = clean.create_slab_label	1✔
91	comp = str(clean.composition.reduced_composition)	1✔
92	assert f"{hkl} {comp}" == label	1✔
93
94	for ads in self.metals_O_entry_dict[el][hkl][clean]:	1✔
95	label = ads.create_slab_label	1✔
96	assert label == f"{hkl} {comp}+O, 0.250 ML"	1✔
97
98	def test_surface_energy(self):	1✔
99	# For a non-stoichiometric case, the chemical potentials do not
100	# cancel out, they serve as a reservoir for any missing atoms
101	for slab_entry in self.MgO_slab_entry_dict[(1, 1, 1)]:	1✔
102	se = slab_entry.surface_energy(self.MgO_ucell_entry, ref_entries=[self.Mg_ucell_entry])	1✔
103	assert tuple(se.as_coefficients_dict()) == (Number(1), Symbol("delu_Mg"))	1✔
104
105	# For the case of a clean, stoichiometric slab, the surface energy
106	# should be constant (i.e. surface energy is a constant).
107	all_se = []	1✔
108	ECu = self.Cu_ucell_entry.energy_per_atom	1✔
109	for val in self.Cu_entry_dict.values():	1✔
110	slab_entry = list(val)[0]	1✔
111	se = slab_entry.surface_energy(self.Cu_ucell_entry)	1✔
112	all_se.append(se)	1✔
113	# Manually calculate surface energy
114	manual_se = (slab_entry.energy - ECu * len(slab_entry.structure)) / (2 * slab_entry.surface_area)	1✔
115	self.assertArrayAlmostEqual(float(se), manual_se, 10)	1✔
116
117	# The (111) facet should be the most stable
118	clean111_entry = list(self.Cu_entry_dict[(1, 1, 1)])[0]	1✔
119	se_Cu111 = clean111_entry.surface_energy(self.Cu_ucell_entry)	1✔
120	assert min(all_se) == se_Cu111	1✔
121
122	def test_cleaned_up_slab(self):	1✔
123	# The cleaned up slab should have the same reduced formula as a clean slab
124	for el, val in self.metals_O_entry_dict.items():	1✔
125	for hkl in val:	1✔
126	for clean in self.metals_O_entry_dict[el][hkl]:	1✔
127	for ads in self.metals_O_entry_dict[el][hkl][clean]:	1✔
128	s = ads.cleaned_up_slab	1✔
129	assert s.composition.reduced_composition == clean.composition.reduced_composition	1✔
130
131
132	class SurfaceEnergyPlotterTest(PymatgenTest):	1✔
133	def setUp(self):	1✔
134	entry_dict = get_entry_dict(os.path.join(get_path(""), "Cu_entries.txt"))	1✔
135	self.Cu_entry_dict = entry_dict	1✔
136	with open(os.path.join(get_path(""), "ucell_entries.txt")) as ucell_entries:	1✔
137	ucell_entries = json.loads(ucell_entries.read())	1✔
138
139	self.Cu_ucell_entry = ComputedStructureEntry.from_dict(ucell_entries["Cu"])	1✔
140	self.Cu_analyzer = SurfaceEnergyPlotter(entry_dict, self.Cu_ucell_entry)	1✔
141
142	self.metals_O_entry_dict = load_O_adsorption()	1✔
143	ucell_entry = ComputedStructureEntry.from_dict(ucell_entries["Pt"])	1✔
144	self.Pt_analyzer = SurfaceEnergyPlotter(self.metals_O_entry_dict["Pt"], ucell_entry)	1✔
145	ucell_entry = ComputedStructureEntry.from_dict(ucell_entries["Ni"])	1✔
146	self.Ni_analyzer = SurfaceEnergyPlotter(self.metals_O_entry_dict["Ni"], ucell_entry)	1✔
147	ucell_entry = ComputedStructureEntry.from_dict(ucell_entries["Rh"])	1✔
148	self.Rh_analyzer = SurfaceEnergyPlotter(self.metals_O_entry_dict["Rh"], ucell_entry)	1✔
149	self.Oads_analyzer_dict = {	1✔
150	"Pt": self.Pt_analyzer,
151	"Ni": self.Ni_analyzer,
152	"Rh": self.Rh_analyzer,
153	}
154
155	def test_get_stable_entry_at_u(self):	1✔
156	for plotter in self.Oads_analyzer_dict.values():	1✔
157	for hkl in plotter.all_slab_entries:	1✔
158	# Test that the surface energy is clean for specific range of chempot
159	entry1, gamma1 = plotter.get_stable_entry_at_u(hkl, delu_dict={Symbol("delu_O"): -7})	1✔
160	entry2, gamma2 = plotter.get_stable_entry_at_u(hkl, delu_dict={Symbol("delu_O"): -6})	1✔
161	assert gamma1 == gamma2	1✔
162	assert entry1.label == entry2.label	1✔
163
164	# Now test that for a high chempot, adsorption
165	# occurs and gamma is not equal to clean gamma
166	entry3, gamma3 = plotter.get_stable_entry_at_u(hkl, delu_dict={Symbol("delu_O"): -1})	1✔
167	assert entry3.label != entry2.label	1✔
168	assert gamma3 != gamma2	1✔
169
170	# For any chempot greater than -6, surface energy should vary
171	# but the configuration should remain the same
172	entry4, gamma4 = plotter.get_stable_entry_at_u(hkl, delu_dict={Symbol("delu_O"): 0})	1✔
173	assert entry3.label == entry4.label	1✔
174	assert gamma3 != gamma4	1✔
175
176	def test_wulff_from_chempot(self):	1✔
177	# Test if it generates a Wulff shape, test if
178	# all the facets for Cu wulff shape are inside.
179	Cu_wulff = self.Cu_analyzer.wulff_from_chempot()	1✔
180	area_frac_dict = Cu_wulff.area_fraction_dict	1✔
181	facets_hkl = [	1✔
182	(1, 1, 1),
183	(3, 3, 1),
184	(3, 1, 0),
185	(1, 0, 0),
186	(3, 1, 1),
187	(2, 1, 0),
188	(2, 2, 1),
189	]
190	for hkl in area_frac_dict:	1✔
191	if hkl in facets_hkl:	1✔
192	assert area_frac_dict[hkl] != 0	1✔
193	else:
194	assert area_frac_dict[hkl] == 0	1✔
195
196	for analyzer in self.Oads_analyzer_dict.values():	1✔
197	# Test WulffShape for adsorbed surfaces
198	# chempot = analyzer.max_adsorption_chempot_range(0)
199	wulff = analyzer.wulff_from_chempot(delu_default=-6)	1✔
200	wulff.weighted_surface_energy	1✔
201
202	# Test if a different Wulff shape is generated
203	# for Ni when adsorption comes into play
204	wulff_neg7 = self.Oads_analyzer_dict["Ni"].wulff_from_chempot(delu_default=-7)	1✔
205	wulff_neg6 = self.Oads_analyzer_dict["Ni"].wulff_from_chempot(delu_default=-6)	1✔
206	assert wulff_neg7.weighted_surface_energy == wulff_neg6.weighted_surface_energy	1✔
207	wulff_neg55 = self.Oads_analyzer_dict["Ni"].wulff_from_chempot(delu_default=-5.5)	1✔
208	assert wulff_neg55.weighted_surface_energy != wulff_neg6.weighted_surface_energy	1✔
209	wulff_neg525 = self.Oads_analyzer_dict["Ni"].wulff_from_chempot(delu_default=-5.25)	1✔
210	assert wulff_neg55.weighted_surface_energy != wulff_neg525.weighted_surface_energy	1✔
211
212	def test_color_palette_dict(self):	1✔
213	for el, val in self.metals_O_entry_dict.items():	1✔
214	analyzer = self.Oads_analyzer_dict[el]	1✔
215	color_dict = analyzer.color_palette_dict()	1✔
216	for hkl in val:	1✔
217	for clean in self.metals_O_entry_dict[el][hkl]:	1✔
218	_ = color_dict[clean]	1✔
219	for ads in self.metals_O_entry_dict[el][hkl][clean]:	1✔
220	_ = color_dict[ads]	1✔
221
222	def test_get_surface_equilibrium(self):	1✔
223	# For clean stoichiometric system, the two equations should
224	# be parallel because the surface energy is a constant. Then
225	# get_surface_equilibrium should return None
226	clean111_entry = list(self.Cu_entry_dict[(1, 1, 1)])[0]	1✔
227	clean100_entry = list(self.Cu_entry_dict[(1, 0, 0)])[0]	1✔
228	soln = self.Cu_analyzer.get_surface_equilibrium([clean111_entry, clean100_entry])	1✔
229	assert not soln	1✔
230
231	# For adsorbed system, we should find one intercept
232	Pt_entries = self.metals_O_entry_dict["Pt"]	1✔
233	clean = list(Pt_entries[(1, 1, 1)])[0]	1✔
234	ads = Pt_entries[(1, 1, 1)][clean][0]	1✔
235	Pt_analyzer = self.Oads_analyzer_dict["Pt"]	1✔
236	soln = Pt_analyzer.get_surface_equilibrium([clean, ads])	1✔
237
238	assert list(soln.values())[0] != list(soln.values())[1]	1✔
239
240	# Check if the number of parameters for adsorption are correct
241	assert (Symbol("delu_O"), Symbol("gamma")) == tuple(soln)	1✔
242	# Adsorbed systems have a b2=(-1Nads) / (Nsurfs Aads)
243	se = ads.surface_energy(Pt_analyzer.ucell_entry, Pt_analyzer.ref_entries)	1✔
244	assert se.as_coefficients_dict()[Symbol("delu_O")] == approx(-1 / (2 * ads.surface_area))	1✔
245
246	def test_stable_u_range_dict(self):	1✔
247	analyzer = list(self.Oads_analyzer_dict.values())[-1]	1✔
248
249	stable_u_range = analyzer.stable_u_range_dict([-1, 0], Symbol("delu_O"), no_doped=False)	1✔
250	all_u = []	1✔
251	for entry in stable_u_range:	1✔
252	all_u.extend(stable_u_range[entry])	1✔
253	assert len(all_u) > 1	1✔
254
255	def test_entry_dict_from_list(self):	1✔
256	# Plug in a list of entries to see if it works
257	all_Pt_slab_entries = []	1✔
258	Pt_entries = self.Pt_analyzer.all_slab_entries	1✔
259	for hkl in Pt_entries:	1✔
260	for clean in Pt_entries[hkl]:	1✔
261	all_Pt_slab_entries.append(clean)	1✔
262	all_Pt_slab_entries.extend(Pt_entries[hkl][clean])	1✔
263
264	a = SurfaceEnergyPlotter(all_Pt_slab_entries, self.Pt_analyzer.ucell_entry)	1✔
265	assert type(a).__name__ == "SurfaceEnergyPlotter"	1✔
266
267	# def test_monolayer_vs_BE(self):
268	# for el in self.Oads_analyzer_dict:
269	# # Test WulffShape for adsorbed surfaces
270	# analyzer = self.Oads_analyzer_dict[el]
271	# plt = analyzer.monolayer_vs_BE()
272	#
273	# def test_area_frac_vs_chempot_plot(self):
274	#
275	# for el in self.Oads_analyzer_dict:
276	# # Test WulffShape for adsorbed surfaces
277	# analyzer = self.Oads_analyzer_dict[el]
278	# plt = analyzer.area_frac_vs_chempot_plot(x_is_u_ads=True)
279	#
280	# def test_chempot_vs_gamma_clean(self):
281	#
282	# plt = self.Cu_analyzer.chempot_vs_gamma_clean()
283	# for el in self.Oads_analyzer_dict:
284	# # Test WulffShape for adsorbed surfaces
285	# analyzer = self.Oads_analyzer_dict[el]
286	# plt = analyzer.chempot_vs_gamma_clean(x_is_u_ads=True)
287	#
288	# def test_chempot_vs_gamma_facet(self):
289	#
290	# for el, val in self.metals_O_entry_dict.items():
291	# for hkl in val:
292	# # Test WulffShape for adsorbed surfaces
293	# analyzer = self.Oads_analyzer_dict[el]
294	# plt = analyzer.chempot_vs_gamma_facet(hkl)
295	# def test_surface_chempot_range_map(self):
296	#
297	# for el, val in self.metals_O_entry_dict.items():
298	# for hkl in val:
299	# # Test WulffShape for adsorbed surfaces
300	# analyzer = self.Oads_analyzer_dict[el]
301	# plt = analyzer.chempot_vs_gamma_facet(hkl)
302
303
304	class WorkfunctionAnalyzerTest(PymatgenTest):	1✔
305	def setUp(self):	1✔
306	self.kwargs = {	1✔
307	"poscar_filename": get_path("CONTCAR.relax1.gz"),
308	"locpot_filename": get_path("LOCPOT.gz"),
309	"outcar_filename": get_path("OUTCAR.relax1.gz"),
310	}
311	self.wf_analyzer = WorkFunctionAnalyzer.from_files(**self.kwargs)	1✔
312
313	def test_shift(self):	1✔
314	wf_analyzer_shift = WorkFunctionAnalyzer.from_files(shift=-0.25, blength=3.7, **self.kwargs)	1✔
315	assert self.wf_analyzer.ave_bulk_p == approx(wf_analyzer_shift.ave_bulk_p, abs=1e-1)	1✔
316
317	def test_is_converged(self):	1✔
318	assert self.wf_analyzer.is_converged()	1✔
319
320
321	class NanoscaleStabilityTest(PymatgenTest):	1✔
322	def setUp(self):	1✔
323	# Load all entries
324	La_hcp_entry_dict = get_entry_dict(os.path.join(get_path(""), "La_hcp_entries.txt"))	1✔
325	La_fcc_entry_dict = get_entry_dict(os.path.join(get_path(""), "La_fcc_entries.txt"))	1✔
326	with open(os.path.join(get_path(""), "ucell_entries.txt")) as ucell_entries:	1✔
327	ucell_entries = json.loads(ucell_entries.read())	1✔
328	La_hcp_ucell_entry = ComputedStructureEntry.from_dict(ucell_entries["La_hcp"])	1✔
329	La_fcc_ucell_entry = ComputedStructureEntry.from_dict(ucell_entries["La_fcc"])	1✔
330
331	# Set up the NanoscaleStabilityClass
332	self.La_hcp_analyzer = SurfaceEnergyPlotter(La_hcp_entry_dict, La_hcp_ucell_entry)	1✔
333	self.La_fcc_analyzer = SurfaceEnergyPlotter(La_fcc_entry_dict, La_fcc_ucell_entry)	1✔
334	self.nanoscale_stability = NanoscaleStability([self.La_fcc_analyzer, self.La_hcp_analyzer])	1✔
335
336	def test_stability_at_r(self):	1✔
337	# Check that we have a different polymorph that is
338	# stable below or above the equilibrium particle size
339	r = self.nanoscale_stability.solve_equilibrium_point(self.La_hcp_analyzer, self.La_fcc_analyzer) * 10	1✔
340
341	# hcp phase of La particle should be the stable
342	# polymorph above the equilibrium radius
343	hcp_wulff = self.La_hcp_analyzer.wulff_from_chempot()	1✔
344	bulk = self.La_hcp_analyzer.ucell_entry	1✔
345	ghcp, rhcp = self.nanoscale_stability.wulff_gform_and_r(hcp_wulff, bulk, r + 10, from_sphere_area=True)	1✔
346	fcc_wulff = self.La_fcc_analyzer.wulff_from_chempot()	1✔
347	bulk = self.La_fcc_analyzer.ucell_entry	1✔
348	gfcc, rfcc = self.nanoscale_stability.wulff_gform_and_r(fcc_wulff, bulk, r + 10, from_sphere_area=True)	1✔
349	assert gfcc > ghcp	1✔
350
351	# fcc phase of La particle should be the stable
352	# polymorph below the equilibrium radius
353	hcp_wulff = self.La_hcp_analyzer.wulff_from_chempot()	1✔
354	bulk = self.La_hcp_analyzer.ucell_entry	1✔
355	ghcp, rhcp = self.nanoscale_stability.wulff_gform_and_r(hcp_wulff, bulk, r - 10, from_sphere_area=True)	1✔
356	fcc_wulff = self.La_fcc_analyzer.wulff_from_chempot()	1✔
357	bulk = self.La_fcc_analyzer.ucell_entry	1✔
358	gfcc, rfcc = self.nanoscale_stability.wulff_gform_and_r(fcc_wulff, bulk, r - 10, from_sphere_area=True)	1✔
359	assert gfcc < ghcp	1✔
360
361	def test_scaled_wulff(self):	1✔
362	# Ensure for a given radius, the effective radius
363	# of the Wulff shape is the same (correctly scaled)
364	hcp_wulff = self.La_hcp_analyzer.wulff_from_chempot()	1✔
365	fcc_wulff = self.La_fcc_analyzer.wulff_from_chempot()	1✔
366	w1 = self.nanoscale_stability.scaled_wulff(hcp_wulff, 10)	1✔
367	w2 = self.nanoscale_stability.scaled_wulff(fcc_wulff, 10)	1✔
368	assert w1.effective_radius == approx(w2.effective_radius)	1✔
369	assert w1.effective_radius == approx(10)	1✔
370	assert 10 == approx(w2.effective_radius)	1✔
371
372
373	def get_entry_dict(filename):	1✔
374	# helper to generate an entry_dict
375
376	entry_dict = {}	1✔
377	with open(filename) as entries:	1✔
378	entries = json.loads(entries.read())	1✔
379	for k in entries:	1✔
380	n = k[25:]	1✔
381	miller_index = []	1✔
382	for i, s in enumerate(n):	1✔
383	if s == "_":	1✔
384	break	1✔
385	if s == "-":	1✔
386	continue	1✔
387	t = int(s)	1✔
388	if n[i - 1] == "-":	1✔
389	t *= -1	1✔
390	miller_index.append(t)	1✔
391	hkl = tuple(miller_index)	1✔
392	if hkl not in entry_dict:	1✔
393	entry_dict[hkl] = {}	1✔
394	entry = ComputedStructureEntry.from_dict(entries[k])	1✔
395	entry_dict[hkl][SlabEntry(entry.structure, entry.energy, hkl, label=k)] = []	1✔
396
397	return entry_dict	1✔
398
399
400	def load_O_adsorption():	1✔
401	# Loads the dictionary for clean and O adsorbed Rh, Pt, and Ni entries
402
403	# Load the adsorbate as an entry
404	with open(os.path.join(get_path(""), "isolated_O_entry.txt")) as isolated_O_entry:	1✔
405	isolated_O_entry = json.loads(isolated_O_entry.read())	1✔
406	O = ComputedStructureEntry.from_dict(isolated_O_entry)	1✔
407
408	# entry_dict for the adsorption case, O adsorption on Ni, Rh and Pt
409	metals_O_entry_dict = {	1✔
410	"Ni": {(1, 1, 1): {}, (1, 0, 0): {}},
411	"Pt": {(1, 1, 1): {}},
412	"Rh": {(1, 0, 0): {}},
413	}
414
415	with open(os.path.join(get_path(""), "csentries_slabs.json")) as entries:	1✔
416	entries = json.loads(entries.read())	1✔
417	for k in entries:	1✔
418	entry = ComputedStructureEntry.from_dict(entries[k])	1✔
419	for el in metals_O_entry_dict: # pylint: disable=C0206	1✔
420	if el in k:	1✔
421	if "111" in k:	1✔
422	clean = SlabEntry(entry.structure, entry.energy, (1, 1, 1), label=k + "_clean")	1✔
423	metals_O_entry_dict[el][(1, 1, 1)][clean] = []	1✔
424	if "110" in k:	1✔
425	clean = SlabEntry(entry.structure, entry.energy, (1, 1, 0), label=k + "_clean")	×
426	metals_O_entry_dict[el][(1, 1, 0)][clean] = []	×
427	if "100" in k:	1✔
428	clean = SlabEntry(entry.structure, entry.energy, (1, 0, 0), label=k + "_clean")	1✔
429	metals_O_entry_dict[el][(1, 0, 0)][clean] = []	1✔
430
431	with open(os.path.join(get_path(""), "csentries_o_ads.json")) as entries:	1✔
432	entries = json.loads(entries.read())	1✔
433	for k in entries:	1✔
434	entry = ComputedStructureEntry.from_dict(entries[k])	1✔
435	for el, val in metals_O_entry_dict.items():	1✔
436	if el in k:	1✔
437	if "111" in k:	1✔
438	clean = list(val[(1, 1, 1)])[0]	1✔
439	ads = SlabEntry(	1✔
440	entry.structure,
441	entry.energy,
442	(1, 1, 1),
443	label=k + "_O",
444	adsorbates=[O],
445	clean_entry=clean,
446	)
447	metals_O_entry_dict[el][(1, 1, 1)][clean] = [ads]	1✔
448	if "110" in k:	1✔
449	clean = list(val[(1, 1, 0)])[0]	×
450	ads = SlabEntry(	×
451	entry.structure,
452	entry.energy,
453	(1, 1, 0),
454	label=k + "_O",
455	adsorbates=[O],
456	clean_entry=clean,
457	)
458	metals_O_entry_dict[el][(1, 1, 0)][clean] = [ads]	×
459	if "100" in k:	1✔
460	clean = list(val[(1, 0, 0)])[0]	1✔
461	ads = SlabEntry(	1✔
462	entry.structure,
463	entry.energy,
464	(1, 0, 0),
465	label=k + "_O",
466	adsorbates=[O],
467	clean_entry=clean,
468	)
469	metals_O_entry_dict[el][(1, 0, 0)][clean] = [ads]	1✔
470
471	return metals_O_entry_dict	1✔
472
473
474	if __name__ == "__main__":	1✔
475	unittest.main()	×

materialsproject / pymatgen / 4075885785

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous