Pymatgen Structure Advanced

from mp_api.client import MPRester

MP_API_KEY = "api-key"

from pymatgen.core import Structure
Si_unitcell = Structure.from_id(id_="mp-149",api_key=MP_API_KEY)
Si_unitcell.to(filename="Si_unitcell.cif")

Si_supercell = Si_unitcell.make_supercell([2,2,2], in_place=False)

print(f"Si unitcell has {len(Si_unitcell.sites)} sites")
print(f"Si supercell has {len(Si_supercell.sites)} sites")

LiFePO4_structure = Structure.from_id(id_="mp-19017",api_key=MP_API_KEY)
LiFePO4_structure.to(fmt="cif", filename="LiFePO4.cif")
LiFePO4_structure.add_oxidation_state_by_element({"Li":1,"Fe":2,"P":5,"O":-2})

print(LiFePO4_structure)

from pymatgen.transformations.standard_transformations import OrderDisorderedStructureTransformation, SubstitutionTransformation, AutoOxiStateDecorationTransformation

transformation = SubstitutionTransformation({"Li+":{"Li+":0.5}})
transformation.apply_transformation(LiFePO4_structure).to(fmt="cif", filename="LiFePO4_substituted.cif")

from pymatgen.alchemy.transmuters import StandardTransmuter

transmuter = StandardTransmuter.from_structures([LiFePO4_structure])
transmuter.append_transformation(AutoOxiStateDecorationTransformation())
transmuter.append_transformation(SubstitutionTransformation({"Li+":{"Li+":0.5}}))
transmuter.append_transformation(OrderDisorderedStructureTransformation(), extend_collection=500)
print ("Total ordering:", len(transmuter))

from pymatgen.alchemy.filters import RemoveDuplicatesFilter

transmuter.apply_filter(RemoveDuplicatesFilter())
print(f"Total unique orderings: {len(transmuter)}")

# Initialize the MP Rester
with MPRester(MP_API_KEY) as mpr:
    docs = mpr.materials.surface_properties.search(
        material_ids=["mp-149", "mp-13"],
    )

# dump results to a json file
import json
with open('materials_data.json', 'w') as f:
    json.dump([doc.dict() for doc in docs], f, indent=4)

for doc in docs:
    print(doc.structure.composition.reduced_formula)
    for surface in doc.surfaces:
        print(f"Miller index {surface.miller_index} {surface.surface_energy}")
        surface_structure = Structure.from_str(surface.structure, fmt="cif")
        surface_structure.to(fmt="cif", filename=f"surface_{doc.structure.composition.reduced_formula}_{surface.miller_index}.cif")

from pymatgen.analysis.wulff import WulffShape

# Create a Wulff shape from the surface energies
for doc in docs:
    wulffshape = WulffShape(
        lattice=doc.structure.lattice,
        miller_list=[surface.miller_index for surface in doc.surfaces],
        e_surf_list=[surface.surface_energy for surface in doc.surfaces],
    )
    print(f"{doc.structure.composition.reduced_formula}: shape factor: {wulffshape.shape_factor:.3f}, anisotropy: {wulffshape.anisotropy:.3f}, weighted surface energy: {wulffshape.weighted_surface_energy:.3f} J/m^2")
    # Plot the Wulff shape
    wulffshape.show()

!pip install wulffpack

from wulffpack import SingleCrystal

for doc in docs:
    surface_energies = {tuple(surface.miller_index): surface.surface_energy for surface in doc.surfaces}
    particle = SingleCrystal(surface_energies)
    particle.view()

from pymatgen.core.surface import SlabGenerator

LiFePO4_structure.add_oxidation_state_by_element({"Fe": 2, "Li": 1, "P": 5, "O": -2})

slabgen = SlabGenerator(initial_structure=LiFePO4_structure, 
                        miller_index=(1, 0, 0), 
                        min_slab_size=2, 
                        min_vacuum_size=10)

all_slabs = slabgen.get_slabs(bonds={("P", "O"): 2})
print(f"For the (001) slab of LiFePO4, there are {len(all_slabs)} terminations.")

for slab in all_slabs:
    print(f"Is polar? {slab.is_polar()}, is symmetric? {slab.is_symmetric()}")

from itertools import product
from pymatgen.core.surface import generate_all_slabs
# Generate all possible miller indices
miller_indices = list(product(range(0,4), repeat=3))[1:] # Skip (0, 0, 0)
    
all_slabs = generate_all_slabs(structure=LiFePO4_structure,
                               max_index=2,
                               min_slab_size=10,
                               min_vacuum_size=10,
                               bonds={("P", "O"): 3})
    
for slab in all_slabs:
    print(f"{slab.miller_index} slabs of LiFePO4")
    print(f"Is polar? {slab.is_polar()}, is symmetric? {slab.is_symmetric()}")
    if slab.is_symmetric() and not slab.is_polar():
        print(f"Is polar? {slab.is_polar()}, is symmetric? {slab.is_symmetric()}")
        slab.to(fmt="cif", filename=f"LiFePO4_{slab.miller_index}.cif")
    

from pymatgen.analysis.interfaces.substrate_analyzer import SubstrateAnalyzer

sub_analyzer = SubstrateAnalyzer(max_area_ratio_tol=0.09,
                                 max_area=800,
                                 max_length_tol=0.03,
                                 max_angle_tol=0.01,
                                 bidirectional=False
                                )

matches = list(sub_analyzer.calculate(film=LiFePO4_structure,
                                      substrate=Si_supercell,
                                      film_millers=[(1,0,0)],
                                      substrate_millers=[(1,1,1)]))
print(f"Found {len(matches)} matches")
for match in matches:
    print(f"Match area: {match.match_area:.2f}, strain {match.von_mises_strain:.3f}")

from pymatgen.analysis.interfaces.coherent_interfaces import CoherentInterfaceBuilder
from pymatgen.analysis.interfaces.zsl import ZSLGenerator


zsl = ZSLGenerator(max_area=800)
cib = CoherentInterfaceBuilder(film_structure=LiFePO4_structure,
                               substrate_structure=Si_unitcell,
                               film_miller=(1,0,0),
                               substrate_miller=(1,1,1),
                               zslgen=zsl)

print(f"Terminations: {cib.terminations}")

interfaces=list(cib.get_interfaces(termination= ('O2_Pmmm_1', 'Si_R-3m_1'),
                                   film_thickness=2,
                                   substrate_thickness=5))
len(interfaces)
interface = interfaces[0]

interface.to(fmt="cif", filename="LiFePO4_Si_interface.cif")