fix(diffraction): correct reciprocal lattice cross terms in get_interplanar_angle (#4685)

get_interplanar_angle computes each plane's reciprocal vector norm and
their dot product using the standard triclinic metric formulas. The
k*l cross term in both r1_norm and r2_norm used cos_gamma_star instead
of cos_alpha_star, and the h*l cross term in the dot product indexed
p1[1] instead of p1[2]. Since a, b, c* worth of terms cancel out
whenever two of the lattice angles are 90 degrees, this was invisible
on the existing cubic/tetragonal/hexagonal test fixtures and only
shows up for genuinely triclinic or low-symmetry monoclinic
structures, where it silently returns the wrong angle.

Added a regression test against an independently computed reference
using the crystallographic reciprocal lattice directly on a triclinic
cell with three distinct angles.

Signed-off-by: Mike German <mike@stepsventures.com>
Co-authored-by: Shyue Ping Ong <shyuep@users.noreply.github.com>

Author: steps-re

src/pymatgen/analysis/diffraction/tem.py

@@ -408,22 +408,22 @@ def get_interplanar_angle(structure: Structure, p1: tuple[int, int, int], p2: tu
             + p1[2] ** 2 * c_star**2
             + 2 * p1[0] * p1[1] * a_star * b_star * cos_gamma_star
             + 2 * p1[0] * p1[2] * a_star * c_star * cos_beta_star
-            + 2 * p1[1] * p1[2] * b_star * c_star * cos_gamma_star
+            + 2 * p1[1] * p1[2] * b_star * c_star * cos_alpha_star
         )
         r2_norm = np.sqrt(
             p2[0] ** 2 * a_star**2
             + p2[1] ** 2 * b_star**2
             + p2[2] ** 2 * c_star**2
             + 2 * p2[0] * p2[1] * a_star * b_star * cos_gamma_star
             + 2 * p2[0] * p2[2] * a_star * c_star * cos_beta_star
-            + 2 * p2[1] * p2[2] * b_star * c_star * cos_gamma_star
+            + 2 * p2[1] * p2[2] * b_star * c_star * cos_alpha_star
         )
         r1_dot_r2 = (
             p1[0] * p2[0] * a_star**2
             + p1[1] * p2[1] * b_star**2
             + p1[2] * p2[2] * c_star**2
             + (p1[0] * p2[1] + p2[0] * p1[1]) * a_star * b_star * cos_gamma_star
-            + (p1[0] * p2[2] + p2[0] * p1[1]) * a_star * c_star * cos_beta_star
+            + (p1[0] * p2[2] + p2[0] * p1[2]) * a_star * c_star * cos_beta_star
             + (p1[1] * p2[2] + p2[1] * p1[2]) * b_star * c_star * cos_alpha_star
         )
         phi = np.arccos(r1_dot_r2 / (r1_norm * r2_norm))

tests/analysis/diffraction/test_tem.py

@@ -217,6 +217,30 @@ def test_interplanar_angle(self):
         phi = tem_calc.get_interplanar_angle(hexagonal, (0, 0, 1), (1, 0, 6))
         assert phi == approx(21.0517, rel=1e-4)
 
+    def test_interplanar_angle_triclinic(self):
+        # Regression test: get_interplanar_angle's reciprocal-lattice-vector
+        # norm and dot-product formulas had copy-paste errors (cos_gamma_star
+        # used instead of cos_alpha_star in the k*l cross term of each norm,
+        # and p1[1] used instead of p1[2] in the h*l cross term of the dot
+        # product). These only show up when alpha/beta/gamma are all
+        # different from 90 degrees and from each other, which none of the
+        # cubic/tetragonal/hexagonal cases above exercise (their cos_alpha_star
+        # / cos_beta_star terms are zero, hiding the bug). Compare against an
+        # independent calculation using the crystallographic reciprocal
+        # lattice directly.
+        tem_calc = TEMCalculator()
+        lattice = Lattice.from_parameters(4.0, 5.0, 6.0, 80, 95, 70)
+        struct = Structure(lattice, ["Si"], [[0, 0, 0]])
+        p1, p2 = (1, 1, 1), (1, 0, 2)
+
+        recip = lattice.reciprocal_lattice_crystallographic
+        v1, v2 = recip.get_cartesian_coords(p1), recip.get_cartesian_coords(p2)
+        cos_phi_expected = np.dot(v1, v2) / (np.linalg.norm(v1) * np.linalg.norm(v2))
+        phi_expected = np.rad2deg(np.arccos(cos_phi_expected))
+
+        phi = tem_calc.get_interplanar_angle(struct, p1, p2)
+        assert phi == approx(phi_expected, rel=1e-6)
+
     def test_get_plot_coeffs(self):
         # Test if x * p1 + y * p2 yields p3.
         tem_calc = TEMCalculator()