Allow Structure.interpolate to extrapolate (#3467)

* added parameter to allow Structure.interpolate to extrapolate beyond the end_structure

* generalized extrapolation to end_amplitude, a cofactor of the distortion vector connecting structure to end_structure. This allows for reverse inter/extrapolation as well as partial inter/extrapolation. With this argument, 0 implies no distortion, 1 implies full distortion to end_structure (default), 0.5 implies distortion to a point halfway between structure and end_structure, and -1 implies full distortion in the opposite direction to end_structure.

* added tests for new end_amplitude argument in Structure.interpolate()

* added note about consistency between start and end structures to docstring -- ensuring consistency is required for obtaining useful results

pymatgen/core/structure.py

@@ -2143,13 +2143,16 @@ def interpolate(
         interpolate_lattices: bool = False,
         pbc: bool = True,
         autosort_tol: float = 0,
+        end_amplitude: float = 1,
     ) -> list[IStructure | Structure]:
         """Interpolate between this structure and end_structure. Useful for
-        construction of NEB inputs.
+        construction of NEB inputs. To obtain useful results, the cell setting
+        and order of sites must consistent across the start and end structures.
 
         Args:
             end_structure (Structure): structure to interpolate between this
-                structure and end.
+                structure and end. Must be in the same setting and have the
+                same site ordering to yield useful results.
             nimages (int,list): No. of interpolation images or a list of
                 interpolation images. Defaults to 10 images.
             interpolate_lattices (bool): Whether to interpolate the lattices.
@@ -2162,6 +2165,13 @@ def interpolate(
                 closest points in this particular structure. This is usually
                 what you want in a NEB calculation. 0 implies no sorting.
                 Otherwise, a 0.5 value usually works pretty well.
+            end_amplitude (float): The fractional amplitude of the endpoint
+                of the interpolation, or a cofactor of the distortion vector
+                connecting structure to end_structure. Thus, 0 implies no
+                distortion, 1 implies full distortion to end_structure
+                (default), 0.5 implies distortion to a point halfway
+                between structure and end_structure, and -1 implies full
+                distortion in the opposite direction to end_structure.
 
         Returns:
             List of interpolated structures. The starting and ending
@@ -2217,7 +2227,7 @@ def interpolate(
 
             end_coords = sorted_end_coords
 
-        vec = end_coords - start_coords
+        vec = end_amplitude * (end_coords - start_coords)
         if pbc:
             vec[:, self.pbc] -= np.round(vec[:, self.pbc])
         sp = self.species_and_occu
@@ -2227,7 +2237,7 @@ def interpolate(
             # interpolate lattice matrices using polar decomposition
             # u is a unitary rotation, p is stretch
             u, p = polar(np.dot(end_structure.lattice.matrix.T, np.linalg.inv(self.lattice.matrix.T)))
-            lvec = p - np.identity(3)
+            lvec = end_amplitude * (p - np.identity(3))
             lstart = self.lattice.matrix.T
 
         for x in images:

tests/core/test_structure.py

@@ -392,6 +392,42 @@ def test_interpolate(self):
         int_s_pbc = struct_pbc.interpolate(struct2_pbc, nimages=2)
         assert_allclose(int_s_pbc[1][0].frac_coords, [1.05, 1.05, 0.55])
 
+        # Test end_amplitude =/= 1
+        coords = [[0, 0, 0], [0.75, 0.5, 0.75]]
+        struct = IStructure(self.lattice, ["Si"] * 2, coords)
+        coords2 = []
+        coords2.extend(([0, 0, 0], [0.5, 0.5, 0.5]))
+        struct2 = IStructure(self.struct.lattice, ["Si"] * 2, coords2)
+        # testing large positive values
+        interpolated_structs = struct.interpolate(struct2, 20, end_amplitude=2)
+        for inter_struct in interpolated_structs:
+            assert inter_struct is not None, "Interpolation Failed!"
+            assert interpolated_structs[0].lattice == inter_struct.lattice
+        assert_array_equal(interpolated_structs[0][1].frac_coords, [0.75, 0.5, 0.75])
+        assert_array_equal(interpolated_structs[10][1].frac_coords, [0.5, 0.5, 0.5])
+        assert_array_equal(interpolated_structs[20][1].frac_coords, [0.25, 0.5, 0.25])
+        # testing large negative values
+        interpolated_structs = struct.interpolate(struct2, 20, end_amplitude=-2)
+        for inter_struct in interpolated_structs:
+            assert inter_struct is not None, "Interpolation Failed!"
+            assert interpolated_structs[0].lattice == inter_struct.lattice
+        assert_array_equal(interpolated_structs[0][1].frac_coords, [0.75, 0.5, 0.75])
+        assert_array_equal(interpolated_structs[10][1].frac_coords, [1.0, 0.5, 1.0])
+        assert_array_equal(interpolated_structs[20][1].frac_coords, [1.25, 0.5, 1.25])
+        # testing partial interpolation
+        interpolated_structs = struct.interpolate(struct2, 5, end_amplitude=-0.5)
+        for inter_struct in interpolated_structs:
+            assert inter_struct is not None, "Interpolation Failed!"
+            assert interpolated_structs[0].lattice == inter_struct.lattice
+        assert_array_equal(interpolated_structs[0][1].frac_coords, [0.75, 0.5, 0.75])
+        assert_array_equal(interpolated_structs[5][1].frac_coords, [0.875, 0.5, 0.875])
+        # testing end_amplitude=0
+        interpolated_structs = struct.interpolate(struct2, 5, end_amplitude=0)
+        for inter_struct in interpolated_structs:
+            assert inter_struct is not None, "Interpolation Failed!"
+            assert interpolated_structs[0].lattice == inter_struct.lattice
+            assert_array_equal(inter_struct[1].frac_coords, [0.75, 0.5, 0.75])
+
     def test_interpolate_lattice(self):
         coords = [[0, 0, 0], [0.75, 0.5, 0.75]]
         struct = IStructure(self.lattice, ["Si"] * 2, coords)
@@ -405,22 +441,60 @@ def test_interpolate_lattice(self):
         assert_allclose(struct2.lattice.angles, int_s[2].lattice.angles)
         int_angles = [110.3976469, 94.5359731, 64.5165856]
         assert_allclose(int_angles, int_s[1].lattice.angles)
+        # Assert that volume is monotonic
+        assert struct2.volume >= int_s[1].volume
+        assert int_s[1].volume >= struct.volume
+
+        # Repeat for end_amplitude = 0.5
+        int_s = struct.interpolate(struct2, 2, interpolate_lattices=True, end_amplitude=0.5)
+        assert_allclose(struct.lattice.abc, int_s[0].lattice.abc)
+        assert_allclose(struct.lattice.angles, int_s[0].lattice.angles)
+        assert_allclose(int_angles, int_s[2].lattice.angles)
+        # Assert that volume is monotonic
+        assert struct2.volume >= int_s[1].volume
+        assert int_s[1].volume >= struct.volume
 
+        # Repeat for end_amplitude = 2
+        int_s = struct.interpolate(struct2, 4, interpolate_lattices=True, end_amplitude=2)
+        assert_allclose(struct.lattice.abc, int_s[0].lattice.abc)
+        assert_allclose(struct.lattice.angles, int_s[0].lattice.angles)
+        assert_allclose(int_angles, int_s[1].lattice.angles)
         # Assert that volume is monotonic
         assert struct2.volume >= int_s[1].volume
         assert int_s[1].volume >= struct.volume
 
+        # Repeat for end_amplitude = -1
+        int_s = struct.interpolate(struct2, 2, interpolate_lattices=True, end_amplitude=-1)
+        assert_allclose(struct.lattice.abc, int_s[0].lattice.abc)
+        assert_allclose(struct.lattice.angles, int_s[0].lattice.angles)
+        int_angles = [127.72010946461334, 86.27613506707404, 56.52554566317311]
+        assert_allclose(int_angles, int_s[1].lattice.angles)
+        # Assert that volume is monotonic (should be shrinking for negative end_amplitude)
+        assert int_s[1].volume <= struct.volume
+        # Assert that coordinate shift is reversed
+        assert_array_equal(int_s[1][1].frac_coords, [0.875, 0.5, 0.875])
+        assert_array_equal(int_s[2][1].frac_coords, [1.0, 0.5, 1.0])
+
     def test_interpolate_lattice_rotation(self):
         l1 = Lattice([[1, 0, 0], [0, 1, 0], [0, 0, 1]])
         l2 = Lattice([[-1.01, 0, 0], [0, -1.01, 0], [0, 0, 1]])
         coords = [[0, 0, 0], [0.75, 0.5, 0.75]]
         struct1 = IStructure(l1, ["Si"] * 2, coords)
         struct2 = IStructure(l2, ["Si"] * 2, coords)
-        int_s = struct1.interpolate(struct2, 2, interpolate_lattices=True)
 
-        # Assert that volume is monotonic
-        assert struct2.volume >= int_s[1].volume
-        assert int_s[1].volume >= struct1.volume
+        # Test positive end_amplitudes
+        for end_amplitude in [0, 0.5, 1, 2]:
+            int_s = struct1.interpolate(struct2, 2, interpolate_lattices=True, end_amplitude=end_amplitude)
+            # Assert that volume is monotonic
+            assert struct2.volume >= int_s[1].volume
+            assert int_s[1].volume >= struct1.volume
+
+        # Test negative end_amplitudes
+        for end_amplitude in [-2, -0.5, 0]:
+            int_s = struct1.interpolate(struct2, 2, interpolate_lattices=True, end_amplitude=end_amplitude)
+            # Assert that volume is monotonic
+            assert struct2.volume >= int_s[1].volume
+            assert int_s[1].volume <= struct1.volume
 
     def test_get_primitive_structure(self):
         coords = [[0, 0, 0], [0.5, 0.5, 0], [0, 0.5, 0.5], [0.5, 0, 0.5]]