Breaking: fix SubstrateAnalyzer film + substrate vectors not using or…

…iginal crystal coordinates (#3572)

* Fix the problem that the film vectors and substrate vectors were not expressed in the original crystal coordinates when making the intreface

This change is to help the users to be capable to track the real film vectors and substrate vectors expressed in their own original crystal coordinates. This is important when they make interface structures because they need to figure out the exact crystal vectors in along the a, b, c vectors in the supercell

Signed-off-by: Jason Xie <[email protected]>

* pre-commit auto-fixes

* breaking: make film + substrate args of generate_surface_vectors()

before, method failed if called before calling calculate() and passing in film + substrate

* add test_generate_surface_vectors

* fix missing doc str white space

---------

Signed-off-by: Jason Xie <[email protected]>
Co-authored-by: Janosh Riebesell <[email protected]>

pymatgen/analysis/interfaces/substrate_analyzer.py

@@ -10,6 +10,8 @@
 from pymatgen.core.surface import SlabGenerator, get_symmetrically_distinct_miller_indices
 
 if TYPE_CHECKING:
+    from numpy.typing import ArrayLike
+
     from pymatgen.core import Structure
 
 
@@ -95,50 +97,59 @@ def __init__(self, film_max_miller=1, substrate_max_miller=1, **kwargs):
         Initializes the substrate analyzer
 
         Args:
-            zslgen(ZSLGenerator): Defaults to a ZSLGenerator with standard
+            zslgen (ZSLGenerator): Defaults to a ZSLGenerator with standard
                 tolerances, but can be fed one with custom tolerances
-            film_max_miller(int): maximum miller index to generate for film
+            film_max_miller (int): maximum miller index to generate for film
                 surfaces
-            substrate_max_miller(int): maximum miller index to generate for
+            substrate_max_miller (int): maximum miller index to generate for
                 substrate surfaces.
         """
         self.film_max_miller = film_max_miller
         self.substrate_max_miller = substrate_max_miller
         self.kwargs = kwargs
         super().__init__(**kwargs)
 
-    def generate_surface_vectors(self, film_millers, substrate_millers):
+    def generate_surface_vectors(
+        self, film: Structure, substrate: Structure, film_millers: ArrayLike, substrate_millers: ArrayLike
+    ):
         """
         Generates the film/substrate slab combinations for a set of given
         miller indices.
 
         Args:
-            film_millers(array): all miller indices to generate slabs for
+            film (Structure): film structure
+            substrate (Structure): substrate structure
+            film_millers (array): all miller indices to generate slabs for
                 film
-            substrate_millers(array): all miller indices to generate slabs
+            substrate_millers (array): all miller indices to generate slabs
                 for substrate
         """
         vector_sets = []
 
-        for f in film_millers:
-            film_slab = SlabGenerator(self.film, f, 20, 15, primitive=False).get_slab()
-            film_vectors = reduce_vectors(film_slab.lattice.matrix[0], film_slab.lattice.matrix[1])
+        for f_miller in film_millers:
+            film_slab = SlabGenerator(film, f_miller, 20, 15, primitive=False).get_slab()
+            film_vectors = reduce_vectors(
+                film_slab.oriented_unit_cell.lattice.matrix[0], film_slab.oriented_unit_cell.lattice.matrix[1]
+            )
 
-            for s in substrate_millers:
-                substrate_slab = SlabGenerator(self.substrate, s, 20, 15, primitive=False).get_slab()
-                substrate_vectors = reduce_vectors(substrate_slab.lattice.matrix[0], substrate_slab.lattice.matrix[1])
+            for s_miller in substrate_millers:
+                substrate_slab = SlabGenerator(substrate, s_miller, 20, 15, primitive=False).get_slab()
+                substrate_vectors = reduce_vectors(
+                    substrate_slab.oriented_unit_cell.lattice.matrix[0],
+                    substrate_slab.oriented_unit_cell.lattice.matrix[1],
+                )
 
-                vector_sets.append((film_vectors, substrate_vectors, f, s))
+                vector_sets.append((film_vectors, substrate_vectors, f_miller, s_miller))
 
         return vector_sets
 
     def calculate(
         self,
-        film,
-        substrate,
+        film: Structure,
+        substrate: Structure,
         elasticity_tensor=None,
-        film_millers=None,
-        substrate_millers=None,
+        film_millers: ArrayLike = None,
+        substrate_millers: ArrayLike = None,
         ground_state_energy=0,
         lowest=False,
     ):
@@ -148,33 +159,28 @@ def calculate(
         ground state energy are provided:
 
         Args:
-            film(Structure): conventional standard structure for the film
-            substrate(Structure): conventional standard structure for the
+            film (Structure): conventional standard structure for the film
+            substrate (Structure): conventional standard structure for the
                 substrate
-            elasticity_tensor(ElasticTensor): elasticity tensor for the film
+            elasticity_tensor (ElasticTensor): elasticity tensor for the film
                 in the IEEE orientation
-            film_millers(array): film facets to consider in search as defined by
+            film_millers (array): film facets to consider in search as defined by
                 miller indices
-            substrate_millers(array): substrate facets to consider in search as
+            substrate_millers (array): substrate facets to consider in search as
                 defined by miller indices
-            ground_state_energy(float): ground state energy for the film
-            lowest(bool): only consider lowest matching area for each surface
+            ground_state_energy (float): ground state energy for the film
+            lowest (bool): only consider lowest matching area for each surface
         """
-        self.film = film
-        self.substrate = substrate
-
         # Generate miller indices if none specified for film
         if film_millers is None:
-            film_millers = sorted(get_symmetrically_distinct_miller_indices(self.film, self.film_max_miller))
+            film_millers = sorted(get_symmetrically_distinct_miller_indices(film, self.film_max_miller))
 
         # Generate miller indices if none specified for substrate
         if substrate_millers is None:
-            substrate_millers = sorted(
-                get_symmetrically_distinct_miller_indices(self.substrate, self.substrate_max_miller)
-            )
+            substrate_millers = sorted(get_symmetrically_distinct_miller_indices(substrate, self.substrate_max_miller))
 
         # Check each miller index combination
-        surface_vector_sets = self.generate_surface_vectors(film_millers, substrate_millers)
+        surface_vector_sets = self.generate_surface_vectors(film, substrate, film_millers, substrate_millers)
         for [
             film_vectors,
             substrate_vectors,

pymatgen/analysis/interfaces/zsl.py

@@ -116,8 +116,8 @@ def generate_sl_transformation_sets(self, film_area, substrate_area):
         lattices with a maximum area
 
         Args:
-            film_area(int): the unit cell area for the film
-            substrate_area(int): the unit cell area for the substrate
+            film_area (int): the unit cell area for the film
+            substrate_area (int): the unit cell area for the substrate
 
         Returns:
             transformation_sets: a set of transformation_sets defined as:

pymatgen/analysis/piezo_sensitivity.py

@@ -229,7 +229,7 @@ def get_rand_IST(self, max_force=1):
         symmetry and the acoustic sum rule.
 
         Args:
-            max_force(float): maximum born effective charge value
+            max_force (float): maximum born effective charge value
 
         Return:
             InternalStrainTensor object

pymatgen/analysis/xas/spectrum.py

@@ -96,7 +96,7 @@ def stitch(self, other: XAS, num_samples: int = 500, mode: Literal["XAFS", "L23"
 
         Args:
             other: Another XAS object.
-            num_samples(int): Number of samples for interpolation.
+            num_samples (int): Number of samples for interpolation.
             mode("XAFS" | "L23"): Either XAFS mode for stitching XANES and EXAFS
                 or L23 mode for stitching L2 and L3.
 

pymatgen/io/aims/inputs.py

@@ -476,7 +476,7 @@ def get_aims_control_parameter_str(self, key: str, value: Any, fmt: str) -> str:
             fmt (str): The format string to apply to the value
 
         Returns:
-            The line to add to the control.in file
+            str: The line to add to the control.in file
         """
         return f"{key:35s}{fmt % value}\n"
 

pymatgen/io/lobster/inputs.py

@@ -619,7 +619,7 @@ def _get_potcar_symbols(POTCAR_input: str) -> list:
         Will return the name of the species in the POTCAR.
 
         Args:
-            POTCAR_input(str): string to potcar file
+            POTCAR_input (str): string to potcar file
 
         Returns:
             list of the names of the species in string format
@@ -662,8 +662,8 @@ def standard_calculations_from_vasp_files(
         Will generate Lobsterin with standard settings.
 
         Args:
-            POSCAR_input(str): path to POSCAR
-            INCAR_input(str): path to INCAR
+            POSCAR_input (str): path to POSCAR
+            INCAR_input (str): path to INCAR
             POTCAR_input (str): path to POTCAR
             dict_for_basis (dict): can be provided: it should look the following:
                 dict_for_basis={"Fe":'3p 3d 4s 4f', "C": '2s 2p'} and will overwrite all settings from POTCAR_input

pymatgen/io/vasp/sets.py

@@ -2162,7 +2162,7 @@ class MVLGBSet(DictSet):
     Class for writing a vasp input files for grain boundary calculations, slab or bulk.
 
     Args:
-        structure(Structure): provide the structure
+        structure (Structure): provide the structure
         k_product: Kpoint number * length for a & b directions, also for c direction in
             bulk calculations. Default to 40.
         slab_mode (bool): Defaults to False. Use default (False) for a bulk supercell.

tests/analysis/interfaces/test_substrate_analyzer.py

@@ -1,20 +1,22 @@
 from __future__ import annotations
 
+from numpy.testing import assert_allclose
+
 from pymatgen.analysis.elasticity.elastic import ElasticTensor
 from pymatgen.analysis.interfaces.substrate_analyzer import SubstrateAnalyzer
 from pymatgen.symmetry.analyzer import SpacegroupAnalyzer
 from pymatgen.util.testing import PymatgenTest
 
+VO2 = PymatgenTest.get_structure("VO2")
+TiO2 = PymatgenTest.get_structure("TiO2")
 
-def test_substrate_analyzer_init():
-    # Film VO2
-    film = SpacegroupAnalyzer(PymatgenTest.get_structure("VO2"), symprec=0.1).get_conventional_standard_structure()
+# Film VO2
+film = SpacegroupAnalyzer(VO2, symprec=0.1).get_conventional_standard_structure()
+# Substrate TiO2
+substrate = SpacegroupAnalyzer(TiO2, symprec=0.1).get_conventional_standard_structure()
 
-    # Substrate TiO2
-    substrate = SpacegroupAnalyzer(
-        PymatgenTest.get_structure("TiO2"), symprec=0.1
-    ).get_conventional_standard_structure()
 
+def test_substrate_analyzer_init():
     film_elastic_tensor = ElasticTensor.from_voigt(
         [
             [324.32, 187.3, 170.92, 0, 0, 0],
@@ -32,3 +34,19 @@ def test_substrate_analyzer_init():
     assert len(matches) == 296
     for match in matches:
         assert isinstance(match.match_area, float)
+
+
+def test_generate_surface_vectors():
+    film_miller_indices = [(1, 0, 0)]
+    substrate_miller_indices = [(1, 1, 1)]
+
+    vector_sets = SubstrateAnalyzer().generate_surface_vectors(
+        film, substrate, film_miller_indices, substrate_miller_indices
+    )
+    assert len(vector_sets) == 1
+    film_vectors, substrate_vectors, film_millers, substrate_millers = vector_sets[0]
+
+    assert [film_millers] == film_miller_indices
+    assert [substrate_millers] == substrate_miller_indices
+    assert_allclose(film_vectors, [[0, 0, 3.035429], [-2.764654e-16, 4.515023, 2.764654e-16]], atol=1e-6)
+    assert_allclose(substrate_vectors, [[-3.766937, -1.928326, -6.328967], [3.766937, -12.307154, 0.0]], atol=1e-6)