From 97970648d3c64df28ff127b3781c24937bf627ae Mon Sep 17 00:00:00 2001 From: henrikjacobsenfys Date: Sat, 17 Oct 2026 13:04:37 +0200 Subject: [PATCH 1/2] PhononSimple_process, Phonon_simple: scale the one-phonon cross-section to Squires 3.120 Both components traced 1/(4 pi x 1.00866) of the one-phonon cross-section they implement (Squires Eq. 3.120), see #2792: - 1/4pi: PhononSimple_process returned mu = rho b^2 X where Union expects mu = 4 pi rho dsigma/dOmega; Phonon_simple divided its weight by 4 pi. - m_n/u: the recoil hbar^2 kappa^2/2m_n was divided by M as if M were in neutron masses, while both headers document M as an atomic mass. The test targets that recorded the old scale are updated: the Phonon_simple based ones (Samples_Phonon, Test_Dispersion_relation_TAS examples and scan) scale by exactly 4 pi x 1.00866 at a fixed seed, and the Union one (Test_inhomogenous_process, sample=phonon) moves by +4.6%. The TAS test's dispersion generator mirrors Phonon_simple's cross-section, so it gets the same two factors. Fixes #2792 Co-Authored-By: Claude Opus 5.5 --- .../examples/Tests_samples/Samples_Phonon/README.md | 4 ++-- .../Samples_Phonon/Samples_Phonon.instr | 4 ++-- .../Test_Dispersion_relation_TAS/README.md | 6 +++--- .../Test_Dispersion_relation_TAS.instr | 12 ++++++------ .../Test_Dispersion_relation_TAS.md | 4 ++-- .../generate_phonon_dispersion.py | 9 +++++---- .../Tests_union/Test_inhomogenous_process/README.md | 2 +- .../Test_inhomogenous_process.instr | 2 +- mcstas-comps/samples/Phonon_simple.comp | 6 +++--- mcstas-comps/union/PhononSimple_process.comp | 4 ++-- 10 files changed, 27 insertions(+), 26 deletions(-) diff --git a/mcstas-comps/examples/Tests_samples/Samples_Phonon/README.md b/mcstas-comps/examples/Tests_samples/Samples_Phonon/README.md index 490327f6df..20439ba665 100644 --- a/mcstas-comps/examples/Tests_samples/Samples_Phonon/README.md +++ b/mcstas-comps/examples/Tests_samples/Samples_Phonon/README.md @@ -18,8 +18,8 @@ Refer to the component documentation for further instructions. ## Examples -- **Test: E=10 -n 1e5 focus_r=0.001 Detector: mon1_I=2.86265e-25** -- **Test: E=10 -n 1e5 focus_a=0.1 Detector: mon1_I=2.86265e-25** +- **Test: E=10 -n 1e5 focus_r=0.001 Detector: mon1_I=3.62846e-24** +- **Test: E=10 -n 1e5 focus_a=0.1 Detector: mon1_I=3.62846e-24** ## Input parameters diff --git a/mcstas-comps/examples/Tests_samples/Samples_Phonon/Samples_Phonon.instr b/mcstas-comps/examples/Tests_samples/Samples_Phonon/Samples_Phonon.instr index 9ba9e83fbd..62c9c9a1bf 100644 --- a/mcstas-comps/examples/Tests_samples/Samples_Phonon/Samples_Phonon.instr +++ b/mcstas-comps/examples/Tests_samples/Samples_Phonon/Samples_Phonon.instr @@ -17,8 +17,8 @@ * Simple test instrument for the Phonon_simple component. * Refer to the component documentation for further instructions. * -* %Example: E=10 -n 1e5 focus_r=0.001 Detector: mon1_I=2.86265e-25 -* %Example: E=10 -n 1e5 focus_a=0.1 Detector: mon1_I=2.86265e-25 +* %Example: E=10 -n 1e5 focus_r=0.001 Detector: mon1_I=3.62846e-24 +* %Example: E=10 -n 1e5 focus_a=0.1 Detector: mon1_I=3.62846e-24 * * %Parameters * E: [meV] Mean energy at source diff --git a/mcstas-comps/examples/Tests_samples/Test_Dispersion_relation_TAS/README.md b/mcstas-comps/examples/Tests_samples/Test_Dispersion_relation_TAS/README.md index 75964742af..d0d68b0933 100644 --- a/mcstas-comps/examples/Tests_samples/Test_Dispersion_relation_TAS/README.md +++ b/mcstas-comps/examples/Tests_samples/Test_Dispersion_relation_TAS/README.md @@ -55,9 +55,9 @@ generate_dispersion_files.py (Python with numpy). ## Examples -- **Test: h=1.2 l=1 dE=5 sample_choice=0 Detector: detector_I=1.46e-11** -- **Test: h=1.2 l=1 dE=5 sample_choice=1 Detector: detector_I=1.41e-11** -- **Test: h=1.2 l=1 dE=5 sample_choice=0 primitive_cell=1 Detector: detector_I=1.46e-11** +- **Test: h=1.2 l=1 dE=5 sample_choice=0 Detector: detector_I=1.85e-10** +- **Test: h=1.2 l=1 dE=5 sample_choice=1 Detector: detector_I=1.79e-10** +- **Test: h=1.2 l=1 dE=5 sample_choice=0 primitive_cell=1 Detector: detector_I=1.85e-10** - **Test: h=1.15 l=0 dE=3 T=2 sample_choice=2 Detector: detector_I=8.91e-10** - **Test: h=1.15 l=0 dE=3 T=2 sample_choice=3 Detector: detector_I=9.41e-10** diff --git a/mcstas-comps/examples/Tests_samples/Test_Dispersion_relation_TAS/Test_Dispersion_relation_TAS.instr b/mcstas-comps/examples/Tests_samples/Test_Dispersion_relation_TAS/Test_Dispersion_relation_TAS.instr index 46b4fe57ad..1948234897 100644 --- a/mcstas-comps/examples/Tests_samples/Test_Dispersion_relation_TAS/Test_Dispersion_relation_TAS.instr +++ b/mcstas-comps/examples/Tests_samples/Test_Dispersion_relation_TAS/Test_Dispersion_relation_TAS.instr @@ -47,15 +47,15 @@ * instrument is compiled, if they do not exist yet, by * generate_dispersion_files.py (Python with numpy). * -* %Example: h=1.2 l=1 dE=5 sample_choice=0 Detector: detector_I=1.46e-11 -* %Example: h=1.2 l=1 dE=5 sample_choice=1 Detector: detector_I=1.41e-11 -* %Example: h=1.2 l=1 dE=5 sample_choice=0 primitive_cell=1 Detector: detector_I=1.46e-11 +* %Example: h=1.2 l=1 dE=5 sample_choice=0 Detector: detector_I=1.85e-10 +* %Example: h=1.2 l=1 dE=5 sample_choice=1 Detector: detector_I=1.79e-10 +* %Example: h=1.2 l=1 dE=5 sample_choice=0 primitive_cell=1 Detector: detector_I=1.85e-10 * %Example: h=1.15 l=0 dE=3 T=2 sample_choice=2 Detector: detector_I=8.91e-10 * %Example: h=1.15 l=0 dE=3 T=2 sample_choice=3 Detector: detector_I=9.41e-10 * %Scan: mcrun Test_Dispersion_relation_TAS.instr -M -N 5,5 -n 2e5 h=0.6,1.4 l=0.6,1.4 dE=5 sample_choice=0 Detector: detector_I={ -* 0.0,0.0,0.0,0.0,0.0,0.0,1.2167e-12,3.52579e-11, -* 2.26406e-12,0.0,0.0,1.17383e-11,0.0,1.39525e-11,0.0,0.0, -* 6.10596e-12,1.00579e-11,8.66923e-12,0.0,0.0,0.0,0.0,0.0, +* 0.0,0.0,0.0,0.0,0.0,0.0,1.5422e-11,4.46901e-10, +* 2.86974e-11,0.0,0.0,1.48785e-10,0.0,1.76851e-10,0.0,0.0, +* 7.73942e-11,1.27486e-10,1.09884e-10,0.0,0.0,0.0,0.0,0.0, * 0.0 } * * %Parameters diff --git a/mcstas-comps/examples/Tests_samples/Test_Dispersion_relation_TAS/Test_Dispersion_relation_TAS.md b/mcstas-comps/examples/Tests_samples/Test_Dispersion_relation_TAS/Test_Dispersion_relation_TAS.md index 6581a5cb40..344fdd616e 100644 --- a/mcstas-comps/examples/Tests_samples/Test_Dispersion_relation_TAS/Test_Dispersion_relation_TAS.md +++ b/mcstas-comps/examples/Tests_samples/Test_Dispersion_relation_TAS/Test_Dispersion_relation_TAS.md @@ -75,7 +75,7 @@ of the two scanned parameters. `generate_phonon_dispersion.py` writes the `Phonon_simple` dispersion of Pb to a grid. The last column is the `Phonon_simple` cross section, -b²/(4π) · (ħ²κ²/2m_n)/ω / M per unit cell, in the form `Dispersion_relation` +b² · (ħ²κ²/2M)/ω per unit cell, M in amu, in the form `Dispersion_relation` reads (see the component header). It depends on the full scattering vector κ, while the file is read at q folded into one unit cell, so the file is written for the Brillouin zone around one reciprocal lattice point (`--zone`, default @@ -95,7 +95,7 @@ broadened by the resolution of the spectrometer: | | `Dispersion_relation` | `Phonon_simple` | |---|---|---| -| Total intensity over the scan (n/s) | 1.572e-9 | 1.585e-9 | +| Total intensity over the scan (n/s) | 1.993e-8 | 2.009e-8 | | Ring radius from intensity (r.l.u.), contour 0.215 | 0.215 | 0.215 | The intensity maps have a correlation of 0.999, and the total intensities diff --git a/mcstas-comps/examples/Tests_samples/Test_Dispersion_relation_TAS/generate_phonon_dispersion.py b/mcstas-comps/examples/Tests_samples/Test_Dispersion_relation_TAS/generate_phonon_dispersion.py index 35cfc1cd73..4c63de2f14 100644 --- a/mcstas-comps/examples/Tests_samples/Test_Dispersion_relation_TAS/generate_phonon_dispersion.py +++ b/mcstas-comps/examples/Tests_samples/Test_Dispersion_relation_TAS/generate_phonon_dispersion.py @@ -10,7 +10,7 @@ The last column is the Phonon_simple cross section written in the form that Dispersion_relation reads, in barn/sr per unit cell of side a/2: - I = (rho_atom/rho_cell) * b^2/(4 pi) * (hbar^2 kappa^2/2m_n)/omega / M. + I = (rho_atom/rho_cell) * b^2 * (hbar^2 kappa^2/2M)/omega, M in amu. It depends on the full scattering vector kappa, which the folded grid cannot hold, so it is written for the Brillouin zone around one reciprocal lattice point given by --zone: each folded coordinate is unfolded to lie within half @@ -32,6 +32,7 @@ import numpy as np HBAR2_OVER_2MN = 2.0721 # hbar^2/(2 m_n) [meV AA^2], K2V^2*VS2E in McStas +NEUTRON_MASS_U = 1.00866 # m_n/u, so that M is in amu def phonon_simple_omega(qx, qy, qz, a, c): @@ -46,7 +47,7 @@ def phonon_simple_omega(qx, qy, qz, a, c): def phonon_simple_cross_section(kappa2, omega, a, b, M, cell_volume=None): """Phonon_simple cross section in barn/sr per unit cell of the given volume. - Phonon_simple weights a ray by rho_atom/(4 pi) * (hbar^2 kappa^2/2m_n)/omega * b^2/M, + Phonon_simple weights a ray by rho_atom * (hbar^2 kappa^2/2m_n)/omega * b^2 * (m_n/u)/M, with b in fm, while Dispersion_relation uses rho_cell * I with I in barn/sr. The default cell is the cube of side a/2. """ @@ -54,7 +55,7 @@ def phonon_simple_cross_section(kappa2, omega, a, b, M, cell_volume=None): rho_cell = 1/(cell_volume if cell_volume is not None else (a/2)**3) b2_barn = b**2*0.01 with np.errstate(divide="ignore", invalid="ignore"): - cross_section = (rho_atom/rho_cell)*b2_barn/(4*np.pi)*HBAR2_OVER_2MN*kappa2/omega/M + cross_section = (rho_atom/rho_cell)*b2_barn*HBAR2_OVER_2MN*kappa2/omega*NEUTRON_MASS_U/M # omega = 0 only at the reciprocal lattice points, where no phonon is excited return np.where(omega > 0, cross_section, 0) @@ -88,7 +89,7 @@ def main(): parser.add_argument("--a", type=float, default=4.95, help="fcc lattice constant [AA] (Pb)") parser.add_argument("--c", type=float, default=10.0, help="Velocity of sound [meV AA] (Pb)") parser.add_argument("--b", type=float, default=9.405, help="Scattering length [fm] (Pb)") - parser.add_argument("--M", type=float, default=207.2, help="Atomic mass [a.u.] (Pb)") + parser.add_argument("--M", type=float, default=207.2, help="Atomic mass [amu] (Pb)") parser.add_argument("--zone", type=float, nargs=3, default=[1, 0, 1], metavar=("H", "K", "L"), help="Reciprocal lattice point, in units of 4 pi/a, whose zone the cross section is written for") parser.add_argument("--cell", choices=["cubic", "primitive"], default="cubic", diff --git a/mcstas-comps/examples/Tests_union/Test_inhomogenous_process/README.md b/mcstas-comps/examples/Tests_union/Test_inhomogenous_process/README.md index 1dc47007eb..4fb6eeb925 100644 --- a/mcstas-comps/examples/Tests_union/Test_inhomogenous_process/README.md +++ b/mcstas-comps/examples/Tests_union/Test_inhomogenous_process/README.md @@ -54,7 +54,7 @@ gaussian function and the heavyside function. - **Test: sample=inc_linear Detector: lin_det_I=0.017457** - **Test: sample=inho_linear Detector: lin_det_I=0.017457** - **Test: sample=inc_and_inho Detector: lin_det_I=0.017457** -- **Test: sample=phonon thick=0.001 Detector: lin_det_I=0.00505839** +- **Test: sample=phonon thick=0.001 Detector: lin_det_I=0.00529178** - **Test: sample=inc_trans Detector: lin_det_I=0.081029** - **Test: sample=inho_trans Detector: lin_det_I=0.081040** - **Test: sample=gauss_heavy Detector: lin_det_I=0.0985658** diff --git a/mcstas-comps/examples/Tests_union/Test_inhomogenous_process/Test_inhomogenous_process.instr b/mcstas-comps/examples/Tests_union/Test_inhomogenous_process/Test_inhomogenous_process.instr index 0a6725f33f..66d063554b 100755 --- a/mcstas-comps/examples/Tests_union/Test_inhomogenous_process/Test_inhomogenous_process.instr +++ b/mcstas-comps/examples/Tests_union/Test_inhomogenous_process/Test_inhomogenous_process.instr @@ -53,7 +53,7 @@ * %Example: sample=inc_linear Detector: lin_det_I=0.017457 * %Example: sample=inho_linear Detector: lin_det_I=0.017457 * %Example: sample=inc_and_inho Detector: lin_det_I=0.017457 -* %Example: sample=phonon thick=0.001 Detector: lin_det_I=0.00505839 +* %Example: sample=phonon thick=0.001 Detector: lin_det_I=0.00529178 * %Example: sample=inc_trans Detector: lin_det_I=0.081029 * %Example: sample=inho_trans Detector: lin_det_I=0.081040 * %Example: sample=gauss_heavy Detector: lin_det_I=0.0985658 diff --git a/mcstas-comps/samples/Phonon_simple.comp b/mcstas-comps/samples/Phonon_simple.comp index 1931584ced..116a40b675 100644 --- a/mcstas-comps/samples/Phonon_simple.comp +++ b/mcstas-comps/samples/Phonon_simple.comp @@ -37,7 +37,7 @@ * sigma_inc: [barns] Incoherent scattering cross section per atom * a: [AA] fcc Lattice constant * b: [fm] Scattering length -* M: [a.u.] Atomic mass +* M: [amu] Atomic mass * c: [meV/AA^(-1)] Velocity of sound * DW: [1] Debye-Waller factor * T: [K] Temperature @@ -484,10 +484,10 @@ TRACE my_a_f = V_my_a_v / v_f; bose_factor = nbose (omega, T); p1 = exp (-(V_my_s * (l_i + l_o) + my_a_i * l_i + my_a_f * l_o)); /* Absorption factor */ - p2 = nf * solid_angle * l_full * V_rho / (4 * PI); /* Focusing factors; assume random choice of n_f possibilities */ + p2 = nf * solid_angle * l_full * V_rho; /* Focusing factors; assume random choice of n_f possibilities */ p3 = (v_f / v_i) * DW * (kappa2 * K2V * K2V * VS2E) / fabs (omega) * bose_factor; /* Cross section factor 1 */ p4 = 2 * VS2E * v_f / J_factor; /* Jacobian of delta functions in cross section */ - p5 = b * b / M; /* Cross section factor 2 */ + p5 = b * b * 1.00866 / M; /* Cross section factor 2; 1.00866 = m_n/u, M in amu */ p *= p1 * p2 * p3 * p4 * p5; } /* else transmit: Neutron did not hit the sample */ free (vf_list); diff --git a/mcstas-comps/union/PhononSimple_process.comp b/mcstas-comps/union/PhononSimple_process.comp index adb3a67886..063548eb2c 100644 --- a/mcstas-comps/union/PhononSimple_process.comp +++ b/mcstas-comps/union/PhononSimple_process.comp @@ -39,7 +39,7 @@ * unit_cell_volume: [AA^3] Unit cell volume * a: [AA] fcc lattice constant * c: [meV*AA] Velocity of sound -* M: [units] Nucleus atomic mass in units +* M: [amu] Atomic mass * b: [fm] Scattring length * T: [K] Temperature * longitudinal: [0/1] Simulate longitudinal branches @@ -505,7 +505,7 @@ INITIALIZE exit (1); } - PhononSimple_storage.my_constant = packing_factor * 4.0 / (a * a * a) * DW * K2V * K2V * VS2E * 2 * VS2E * K2V * b * b / M; + PhononSimple_storage.my_constant = packing_factor * 4.0 / (a * a * a) * DW * K2V * K2V * VS2E * 2 * VS2E * K2V * b * b * 4 * PI * 1.00866 / M; /* 4 PI: mu = 4 PI rho dsigma/dOmega; 1.00866 = m_n/u, M in amu */ // my_constant = V_rho*DW*K2V*K2V*VS2E*2*VS"E*K2V*b*b/M // done in initialize // First initialise This_process with default values: From 2fdabb52cea2ca22b66f23d71a67484d42f485ed Mon Sep 17 00:00:00 2001 From: Peter Willendrup Date: Sun, 18 Oct 2026 22:09:07 +0200 Subject: [PATCH 2/2] =?UTF-8?q?Test=5FDispersion=5Frelation=5FTAS:=20resca?= =?UTF-8?q?le=20the=20moved=20%Scan=20targets=20by=204=CF=80=20=C3=97=201.?= =?UTF-8?q?00866?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The merge of main brought in the new scan mesh (h=0.8..1.2, l=1.1..1.3) with references recorded before the Phonon_simple fix, so test 6 failed at ~13x. The references (mean of 16 seeds) are scaled by the exact factor 12.6752 the fix applies at fixed seed. mctest with this branch (macOS arm64, seed 1000): 25/25 points within 20% without MPI (worst 106%) and with 2 MPI ranks (worst 105%). Co-Authored-By: Claude Opus 5.5 --- .../Test_Dispersion_relation_TAS.instr | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/mcstas-comps/examples/Tests_samples/Test_Dispersion_relation_TAS/Test_Dispersion_relation_TAS.instr b/mcstas-comps/examples/Tests_samples/Test_Dispersion_relation_TAS/Test_Dispersion_relation_TAS.instr index b3f2ef3d62..9dfbdf56b6 100644 --- a/mcstas-comps/examples/Tests_samples/Test_Dispersion_relation_TAS/Test_Dispersion_relation_TAS.instr +++ b/mcstas-comps/examples/Tests_samples/Test_Dispersion_relation_TAS/Test_Dispersion_relation_TAS.instr @@ -53,10 +53,10 @@ * %Example: h=1.15 l=0 dE=3 T=2 sample_choice=2 Detector: detector_I=8.91e-10 * %Example: h=1.15 l=0 dE=3 T=2 sample_choice=3 Detector: detector_I=9.41e-10 * %Scan: mcrun Test_Dispersion_relation_TAS.instr -M -N 5,5 -n 2e5 h=0.8,1.2 l=1.1,1.3 dE=5 sample_choice=0 Detector: detector_I={ -* 8.964e-12,1.086e-11,8.771e-12,1.261e-11,1.151e-11,7.748e-12,1.057e-11,1.137e-11, -* 1.199e-11,1.014e-11,6.126e-12,8.926e-12,1.006e-11,1.007e-11,8.424e-12,3.55e-12, -* 7.184e-12,8.291e-12,8.243e-12,5.35e-12,9.55e-13,4.109e-12,5.612e-12,4.924e-12, -* 1.74e-12 } +* 1.136e-10,1.377e-10,1.112e-10,1.598e-10,1.459e-10,9.821e-11,1.34e-10,1.441e-10, +* 1.52e-10,1.285e-10,7.765e-11,1.131e-10,1.275e-10,1.276e-10,1.068e-10,4.5e-11, +* 9.106e-11,1.051e-10,1.045e-10,6.781e-11,1.21e-11,5.208e-11,7.113e-11,6.241e-11, +* 2.205e-11 } * * %Parameters * h: [r.l.u.] Scattering vector component along a* (sample x axis)