Coverage for dynasor / post_processing / neutron_scattering_lengths.py: 100%
59 statements
« prev ^ index » next coverage.py v7.13.4, created at 2026-03-16 12:31 +0000
1import json
2from importlib.resources import files
3from typing import Optional
5import numpy as np
6from pandas import DataFrame
7from .weights import Weights
10class NeutronScatteringLengths(Weights):
11 """This class provides sample weights corresponding to neutron scattering lengths.
12 By default, the coherent and incoherent scattering lengths are weighted by the natural
13 abundance of each isotope of the considered atomic species.
14 This weighting can be overwritten using the :attr:`abundances` argument.
16 The scattering lengths have been extracted from `this NIST
17 database <https://www.ncnr.nist.gov/resources/n-lengths/list.html>`__,
18 which in turn have been taken from Table 1 of Neutron News **3**, 26 (1992);
19 `doi: 10.1080/10448639208218770 <https://doi.org/10.1080/10448639208218770>`_.
21 Parameters
22 ----------
23 atom_types
24 List of atomic species for which to retrieve scattering lengths.
25 abundances
26 Dict of the desired fractional abundance of each isotope for
27 each species in the sample. For example, to use an equal
28 weighting of all isotopes of oxygen, one can write
29 ``abundances['O'] = dict(16=1/3, 17=1/3, 18=1/3)``. Note that
30 the abundance for any isotopes that are *not* included in this
31 dict is automatically set to zero. In other words, you need to
32 ensure that the abundances provided sum up to 1. By default
33 the neutron scattering lengths are weighted proportionally to
34 the natural abundance of each isotope.
35 """
37 def __init__(
38 self,
39 atom_types: list[str],
40 abundances: Optional[dict[str, dict[int, float]]] = None,
41 ):
42 scat_lengths = _read_scattering_lengths()
44 # Sub select only the relevant species
45 scat_lengths = scat_lengths[scat_lengths.species.isin(atom_types)].reset_index()
47 # Update the abundances if another weighting is desired
48 if abundances is not None:
49 for species in abundances:
50 scat_lengths.loc[scat_lengths.species == species, 'abundance'] = 0
51 for Z, frac in abundances[species].items():
52 match = (scat_lengths.species == species) & (scat_lengths.isotope == Z)
53 if not np.any(match): # Check if any row+column matches
54 raise ValueError(f'No match in database for {species} and isotope {Z}')
55 scat_lengths.loc[match, 'abundance'] = frac
57 self._scattering_lengths = scat_lengths
59 # Check if any of the fetched scattering lengths is None,
60 # indicating that it is missing in the experimental database.
61 # Only raise an error if the desired abundance is greater than 0.
62 nan_rows = scat_lengths[scat_lengths.isnull().any(axis=1) & (scat_lengths.abundance > 0.0)]
63 if not nan_rows.empty:
64 # Grab first offending entry
65 row = nan_rows.iloc[0]
66 raise ValueError(f'{row.isotope}{row.species} is missing tabulated values for either'
67 ' the coherent or incoherent scattering length.'
68 ' Adjust the abundance parameter to set the fraction of'
69 f' {row.isotope}{row.species} to zero.')
71 # Make sure abundances add up to 100%
72 by_species = self._scattering_lengths.groupby('species')
73 for species, species_df in by_species:
74 if not np.isclose(species_df.abundance.sum(), 1):
75 raise ValueError(f'Abundance values for {species} do not sum up to 1.0')
77 # Compute scattering lengths weighted by abundance
78 weights_coh = by_species.apply(
79 lambda s: (s.b_coh * s.abundance).sum(),
80 include_groups=False
81 ).to_dict()
82 # First compute the average scattering length, then take the square
83 # since the incoherent scattering lengths enter as b_incoh**2, but
84 # dynasor only applies a single weighting factor w_incoh.
85 weights_inc = by_species.apply(
86 lambda s: (s.b_inc * s.abundance).sum(),
87 include_groups=False
88 ).to_dict()
90 supports_currents = False
91 super().__init__(weights_coh, weights_inc, supports_currents)
93 @property
94 def abundances(self) -> dict[str, dict[int, float]]:
95 """Abundances used for calculating scattering lengths."""
96 abundance_dict = {}
97 for (species), species_df in self._scattering_lengths.groupby('species'):
98 abundance_dict[species] = {}
99 for (isotope, abundance), _ in species_df.groupby(['isotope', 'abundance']):
100 abundance_dict[species][isotope] = abundance
101 return abundance_dict
103 @property
104 def parameters(self) -> DataFrame:
105 """Scattering lengths used to compute the coherent and
106 incoherent weights for the selected isotopes.
107 """
108 return self._scattering_lengths
111def _read_scattering_lengths() -> DataFrame:
112 """
113 Extracts the scattering lengths from the file `neutron_scattering_lengths.json`
114 for each of the supplied species. Scattering lengths are in units of fm.
116 The scattering lengths have been extracted from the following NIST
117 database: https://www.ncnr.nist.gov/resources/n-lengths/list.html,
118 which in turn have been extracted from
119 Neutron News **3**, No. 3***, 26 (1992).
120 """
121 data_file = files(__package__) / 'form-factors/neutron_scattering_lengths.json'
122 with open(data_file) as fp:
123 scattering_lengths = json.load(fp)
125 data = []
126 for species in scattering_lengths:
127 for isotope in scattering_lengths[species]:
128 for fld in 'b_coh b_inc'.split():
129 val = scattering_lengths[species][isotope][fld]
130 if 'None' in val:
131 scattering_lengths[species][isotope][fld] = np.nan
132 elif 'j' in val:
133 scattering_lengths[species][isotope][fld] = complex(val)
134 else:
135 scattering_lengths[species][isotope][fld] = float(val)
136 data.append(
137 dict(
138 species=species,
139 isotope=int(isotope),
140 abundance=float(scattering_lengths[species][isotope]['abundance'])
141 / 100, # % -> fraction
142 b_coh=complex(scattering_lengths[species][isotope]['b_coh']),
143 b_inc=complex(scattering_lengths[species][isotope]['b_inc']),
144 )
145 )
146 scattering_lengths = DataFrame.from_dict(data)
147 return scattering_lengths