From f2aa65c8a54e6f94dbff8d654160e23064130c7f Mon Sep 17 00:00:00 2001 From: "T. Fuchs" Date: Wed, 29 Jul 2026 20:23:37 -0400 Subject: [PATCH] fix(phys): accumulate epitaxy strain per layer --- CHANGELOG.md | 5 + doc/source/ctr_structure_factors.rst | 15 +++ doc/source/release_notes.rst | 11 +++ orgui/datautils/xrayutils/CTRfilm.py | 58 ++++++++--- .../datautils/xrayutils/test/test_CTRcalc.py | 99 ++++++++++++++++--- 5 files changed, 158 insertions(+), 30 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 980a086..cbf02fa 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -37,6 +37,11 @@ Scientific and analysis additions: Scientific correctness and performance fixes: +- Corrected ``EpitaxyInterface`` to accumulate the strain-coupled height of + every structural layer. The previous calculation advanced a complete unit + cell using only its final layer's strain, creating a layer-cycle-periodic + displacement and artificial quasi-Bragg intensity at systematic extinctions + for broad interfaces. - Optical profiles now preserve areal optical content under surface-normal strain, keep ionic forward scattering factors for charged species, avoid merging layers beyond the requested z tolerance, and remain finite at the diff --git a/doc/source/ctr_structure_factors.rst b/doc/source/ctr_structure_factors.rst index b38e90a..118a3b0 100644 --- a/doc/source/ctr_structure_factors.rst +++ b/doc/source/ctr_structure_factors.rst @@ -240,6 +240,21 @@ the linear interpolation z_m(\kappa,0) = z_{0,m} + \kappa\left(z_1-z_{0,m}\right). +For the fully strain-coupled discrete profile, every structural layer supplies +its own physical interval. If :math:`d_{\mathrm{top},i}` and +:math:`d_{\mathrm{bottom},i}` are the upper- and lower-lattice intervals for +layer :math:`i`, the coupled interval is + +.. math:: + + d_{1,i} + = P_i d_{\mathrm{top},i} + + (1-P_i)d_{\mathrm{bottom},i}, + +and the layer origins are the cumulative sum of these intervals. The +calculation must therefore accumulate all structural layers independently; +one layer's strain cannot represent a complete multi-layer unit-cell cycle. + The ``offset`` parameter :math:`o` is expressed in fractional lower-bulk :math:`c` coordinates. Its physical displacement is diff --git a/doc/source/release_notes.rst b/doc/source/release_notes.rst index c7eb1c2..e414878 100644 --- a/doc/source/release_notes.rst +++ b/doc/source/release_notes.rst @@ -9,11 +9,22 @@ Unreleased (2026-07-19) Scientific and analysis additions: +- Added interface-based optical ``delta``/``beta`` profiles and layered wavefield calculations for s- and p-polarized X-rays. Scalar and multidimensional angle inputs are supported, and specular reflectivity can be calculated for s, p, or unpolarized incidence. This provides the optical wavefield foundation for future distorted-wave Born approximation (DWBA) calculations; the full distorted-wave scattering amplitude is not included yet. - Added optional CTR intensity-resolution modeling with constant or gamma-dependent box and Gaussian functions. Calculations can convolve irregular existing L points or sample the crystal structure factor with deterministic quadrature. +- Added Poisson-distributed surface occupancies and coherent out-of-plane epitaxy/strain coupling for film-interface models. - Added py3Dmol atom-sphere rendering for Jupyter notebooks. ``plot3d`` now selects py3Dmol automatically in a notebook, can be directed to either py3Dmol or Mayavi explicitly, and can incrementally add unit cells to a shared viewer. Covalent radii are interpreted consistently by both backends and can be adjusted with the dimensionless ``radius_scale`` parameter. - ``PoissonSurface``'s basis now requires three parameters (``W``, ``alpha``, ``offset``) instead of two. The previous two-parameter form was only ever used in test code, so no migration path is provided; a saved ``.xtal``/ ``.xpr`` file with a two-parameter ``PoissonSurface`` basis will fail to load. - Wyckoff-parameterized fit values are now interpreted as absolute atomic positions instead of deltas relative to the Wyckoff position. Wyckoff parameter fitting was development-only and never part of a release, so no migration is provided. +Scientific correctness and performance fixes: + +- Corrected ``EpitaxyInterface`` to accumulate the strain-coupled height of every structural layer. The previous calculation advanced a complete unit cell using only its final layer's strain, creating a layer-cycle-periodic displacement and artificial quasi-Bragg intensity at systematic extinctions for broad interfaces. +- Optical profiles now preserve areal optical content under surface-normal strain, keep ionic forward scattering factors for charged species, avoid merging layers beyond the requested z tolerance, and remain finite at the exact p-polarized critical-angle limit. +- Optical reflectivity and full wavefield calculations now evaluate angle arrays with NumPy-vectorized kernels while preserving the caller's input shape. +- Corrected film and epitaxy-interface anchoring, stacking, and support ownership across strained interfaces. +- Fixed missing support in one-dimensional electron-density calculations and corrected atomic-coordinate stacking when splitting unit cells into layers. +- Added a C++ electron-density backend and bounded caches for atomic form factors, anomalous scattering factors, and accelerated form-factor lookup. + A ***critical bug*** was fixed that affects bulk CTR calculations: - ``UnitCell.F_bulk``'s semi-infinite geometric lattice sum used the raw, untransformed ``l`` index instead of the index converted by ``refHKLTransform`` when computing the out-of-plane attenuation phase. This was only correct for the default case where a component uses its own bulk cell as the reference (no ``reference_uc`` set); any explicit ``reference_uc`` whose out-of-plane reciprocal axis differs from the bulk's — including a plain scale difference between the reference and bulk out-of-plane axis length, not only a rotated or reindexed reference — gave incorrect bulk structure-factor amplitudes. This bug was present in both the accelerated (numba/C++) and plain-Python code paths in all previous released versions, up to and including v1.5.0. See the CTR structure-factor documentation for details. diff --git a/orgui/datautils/xrayutils/CTRfilm.py b/orgui/datautils/xrayutils/CTRfilm.py index 9d80cb4..b380bbe 100644 --- a/orgui/datautils/xrayutils/CTRfilm.py +++ b/orgui/datautils/xrayutils/CTRfilm.py @@ -763,23 +763,46 @@ def createInterfaceCells(self): mat_0 = np.vstack((np.identity(3).T, np.array([0, 0, 0]))).T - ratio_top = a3_bottom / a3_top - ratio_bottom = 1 / ratio_top coupled_top = [[] for _ in self.top_layers] coupled_bottom = [[] for _ in self.bottom_layers] ideal_top = [[] for _ in self.top_layers] ideal_bottom = [[] for _ in self.bottom_layers] coupled_height = 0.0 + top_layer_space = np.diff( + self.layerpos, + append=self.layerpos[0] + 1.0, + ) + bottom_layerpos = _unwrapped_layer_positions( + np.array( + [ + self.uc_bottom.layerpos[layer_id] + for layer_id in self._layer_ids + ] + ), + self._layer_order_indices, + ) + bottom_layer_space = np.diff( + bottom_layerpos, + append=bottom_layerpos[0] + 1.0, + ) for cycle_index, (p_t, p_b) in enumerate( zip(probability_top, probability_bottom) ): - top_strains = p_t + ratio_top * p_b - bottom_strains = ratio_bottom * p_t + p_b - physical_top_index = np.flatnonzero( - self._layer_order_indices == n_layers - 1 - )[0] - cell_height = a3_top * top_strains[physical_top_index] + # Accumulate the coherent height one structural layer at a + # time. Using a single layer's strain for the whole unit-cell + # cycle creates a displacement modulation with the layer-cycle + # period and leaks intensity into systematic extinctions. + physical_layer_heights = ( + p_t * a3_top * top_layer_space + + p_b * a3_bottom * bottom_layer_space + ) + top_strains = physical_layer_heights / ( + a3_top * top_layer_space + ) + bottom_strains = physical_layer_heights / ( + a3_bottom * bottom_layer_space + ) for i, (uc_t, uc_b) in enumerate( zip(self.top_layers, self.bottom_layers) @@ -787,27 +810,34 @@ def createInterfaceCells(self): mat_top_i = np.copy(mat_0) top_strain_and_h = top_strains[i] mat_top_i[2, 2] = top_strain_and_h - relative_layer_position = self.layerpos[i] - self.layerpos[0] + relative_top_position = ( + self.layerpos[i] - self.layerpos[0] + ) top_layer_offset = ( - relative_layer_position + relative_top_position - (self.uc_top.layerpos[self.layer_order[i]]) ) mat_top_i[2, 3] = ( coupled_height / a3_top - + top_strain_and_h * top_layer_offset + - top_strain_and_h + * self.uc_top.layerpos[self.layer_order[i]] ) coupled_top[i].append(mat_top_i) mat_bottom_i = np.copy(mat_0) bottom_strain_and_h = bottom_strains[i] mat_bottom_i[2, 2] = bottom_strain_and_h + relative_bottom_position = ( + bottom_layerpos[i] - bottom_layerpos[0] + ) bottom_layer_offset = ( - relative_layer_position + relative_bottom_position - (self.uc_bottom.layerpos[self.layer_order[i]]) ) mat_bottom_i[2, 3] = ( coupled_height / a3_bottom - + bottom_strain_and_h * bottom_layer_offset + - bottom_strain_and_h + * self.uc_bottom.layerpos[self.layer_order[i]] ) coupled_bottom[i].append(mat_bottom_i) @@ -818,7 +848,7 @@ def createInterfaceCells(self): ideal_bottom_i = np.copy(mat_0) ideal_bottom_i[2, 3] = cycle_index + bottom_layer_offset ideal_bottom[i].append(ideal_bottom_i) - coupled_height += cell_height + coupled_height += physical_layer_heights[i] loc_rescaled = loc - unitcells[0] uc_no_loc = int(np.floor(loc_rescaled)) // n_layers diff --git a/orgui/datautils/xrayutils/test/test_CTRcalc.py b/orgui/datautils/xrayutils/test/test_CTRcalc.py index 4ec59c1..b968bee 100644 --- a/orgui/datautils/xrayutils/test/test_CTRcalc.py +++ b/orgui/datautils/xrayutils/test/test_CTRcalc.py @@ -2028,6 +2028,42 @@ def test_epitaxy_strain_coupled_field_is_anchored_to_deep_bulk(self): interface.height_absolute - independent.height_absolute, ) + def test_epitaxy_accumulates_each_structural_layer_height(self): + interface = self.make_epitaxy_interface(strain_coupling=1.0) + n_layers = len(interface.layers) + unitcells = np.arange(-interface.fixed_ucs, interface.fixed_ucs) + probability_top = interface.profile.occupancy( + unitcells, + n_layers, + ).reshape((-1, n_layers)) + probability_bottom = 1.0 - probability_top + layer_space = np.diff( + interface.layerpos, + append=interface.layerpos[0] + 1.0, + ) + expected_layer_heights = ( + probability_top * interface.uc_top.a[2] * layer_space + + probability_bottom * interface.uc_bottom.a[2] * layer_space + ) + + generated_origins = [] + for cycle_index in range(probability_top.shape[0]): + for layer_index, layer in enumerate(interface.top_layers): + matrix = layer.coherentDomainMatrix[cycle_index] + layer_id = interface.layer_order[layer_index] + generated_origins.append( + ( + matrix[2, 2] * interface.uc_top.layerpos[layer_id] + + matrix[2, 3] + ) + * interface.uc_top.a[2] + ) + + np.testing.assert_allclose( + np.diff(generated_origins), + expected_layer_heights.ravel()[:-1], + ) + def test_epitaxy_strain_coupling_interpolates_generated_positions(self): independent = self.make_epitaxy_interface(strain_coupling=0.0) halfway = self.make_epitaxy_interface(strain_coupling=0.5) @@ -2373,8 +2409,8 @@ def test_epitaxy_split_bulk_domains_feed_all_profile_paths(self): np.testing.assert_allclose( structure_factor, [ - 12.290590935844 + 23.878872280467j, - -1.767926251260 - 0.197899416448j, + 9.783621364160 + 25.055258758458j, + -0.401679480851 - 1.415519651389j, ], rtol=1e-8, ) @@ -2382,16 +2418,16 @@ def test_epitaxy_split_bulk_domains_feed_all_profile_paths(self): density[[0, 100, 200, 300, 400]], [ 0.0, - -3.006212046217e-03 - 1.665286072977e-05j, - 9.802746914877e-01 + 3.098533417278e-03j, - 1.082893246501e-05, + -3.093100733396e-03 - 1.727080407634e-05j, + 7.183801385126e-01 + 2.427528867636e-03j, + 8.622439679978e-05, 0.0, ], atol=1e-9, ) np.testing.assert_allclose( np.sum(density), - 40.033115045982 + 0.034943462837j, + 40.028295180092 + 0.034223926465j, rtol=1e-8, ) mid_row = optical_profile.shape[0] // 2 @@ -2399,7 +2435,7 @@ def test_epitaxy_split_bulk_domains_feed_all_profile_paths(self): optical_profile[[0, mid_row, -1]], [ [-12.0, -7.72508074e-10, -7.18561204e-13], - [-0.209418853, 2.94939896e-06, 2.74343239e-09], + [0.0, 2.94939896e-06, 2.74343239e-09], [10.0, 0.0, 0.0], ], atol=1e-9, @@ -2407,6 +2443,36 @@ def test_epitaxy_split_bulk_domains_feed_all_profile_paths(self): class TestLegacyLayeredCTR(unittest.TestCase): + def test_broad_interface_does_not_create_rutile_quasi_bragg_peak(self): + fixture_root = os.path.join(os.path.dirname(__file__), "testdata") + xtal = CTRcalc.SXRDCrystal.fromFile( + os.path.join(fixture_root, "RuO2_TiO2_Poisson_etching.xtal") + ) + surface = xtal["RuO2surface"] + film = xtal["RuO2"] + interface = xtal["TiO2toRuO2"] + surface_width = 10.0 / ( + film.unitcell.a[2] / len(film.layers) + ) + interface_width = 10.0 / interface.uc_bottom.a[2] + surface.basis[:] = [surface_width, 1.0, 0.0] + film.basis[0] = 105.0 + interface.basis[:] = [interface_width, 0.0, 1.0, 0.0] + + l_values = np.array([3.9, 4.0, 4.1]) + amplitudes = np.abs( + xtal.F( + np.ones(3, dtype=np.float64), + np.zeros(3, dtype=np.float64), + l_values, + ) + ) + + self.assertLess( + amplitudes[1], + np.mean(amplitudes[[0, 2]]), + ) + def test_legacy_xtal_uses_corrected_interface_support(self): fixture_root = os.path.join(os.path.dirname(__file__), "testdata") xtal_path = os.path.join( @@ -2440,8 +2506,9 @@ def test_legacy_xtal_uses_corrected_interface_support(self): # UNTESTED reference values: computed with this module's own # CTRcalc.SXRDCrystal.F() after the interfacial-strain fix in # EpitaxyInterface (commit beee738), the F_bulk lattice-sum - # denominator fix, and the deep-bulk strain-field anchoring fix, but - # not independently verified against real data or an analytic result + # denominator fix, the deep-bulk strain-field anchoring fix, and the + # cumulative per-structural-layer interface-height fix, but not + # independently verified against real data or an analytic result # -- unlike test_simple_unitcell_bulk_structure_factor_matches_ # reference below, which does compare against real fit-derived # intensities for the un-layered unit-cell + bulk case. @@ -2464,13 +2531,13 @@ def test_legacy_xtal_uses_corrected_interface_support(self): expected = np.array( [ 17238.5440045948, - 7001.3340402313, - 3852.2553000918, - 2726.9468772925, - 2175.6251991885, - 1855.7582644482, - 1648.1329478054, - 1501.1768061163, + 7001.3031826720, + 3852.1926461220, + 2726.8587419552, + 2175.5183472486, + 1855.6392611730, + 1648.0076716502, + 1501.0502790582, ] ) np.testing.assert_allclose(calculated, expected, rtol=2e-5, atol=2e-5)