Coverage for local_installation/dynasor/post_processing/neutron_scattering_lengths.py: 100%
56 statements
« prev ^ index » next coverage.py v7.3.2, created at 2024-09-27 15:43 +0000
« prev ^ index » next coverage.py v7.3.2, created at 2024-09-27 15:43 +0000
1import json
2from importlib.resources import files
3from typing import Dict, List
5import numpy as np
6from pandas import DataFrame
7from .weights import Weights
10class NeutronScatteringLengths(Weights):
11 """A class for generating 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, however, be overwritten by the :attr:`abundances` argument.
16 The scattering lengths have been extracted from the following NIST
17 database: https://www.ncnr.nist.gov/resources/n-lengths/list.html,
18 which in turn have been extracted from
19 Neutron News **3**, 29 (1992).
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: 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'Non-zero abundance of {row.isotope}{row.species}'
67 ' with missing tabulated scattering length.')
69 # Make sure abundances add up to 100%
70 by_species = self._scattering_lengths.groupby('species')
71 for species, species_df in by_species:
72 if not np.isclose(species_df.abundance.sum(), 1):
73 raise ValueError(f'Abundance values for {species} do not sum up to 1.0')
75 # Compute scattering lengths weighted by abundance
76 weights_coh = by_species.apply(
77 lambda s: (s.b_coh * s.abundance)
78 .sum().real # weights in dynasor can only be real atm
79 ).to_dict()
80 # First compute the average scattering length, then take the square
81 # since the incoherent scattering lengths enter as b_incoh**2, but
82 # dynasor only applies a single weighting factor w_incoh.
83 weights_inc = by_species.apply(
84 lambda s: ((s.b_inc * s.abundance).sum()**2).real
85 ).to_dict()
87 supports_currents = False
88 super().__init__(weights_coh, weights_inc, supports_currents)
90 @property
91 def abundances(self) -> Dict[str, Dict[int, float]]:
92 abundance_dict = {}
93 for (species), species_df in self._scattering_lengths.groupby('species'):
94 abundance_dict[species] = {}
95 for (isotope, abundance), _ in species_df.groupby(['isotope', 'abundance']):
96 abundance_dict[species][isotope] = abundance
97 return abundance_dict
100def _read_scattering_lengths() -> DataFrame:
101 """
102 Extracts the scattering lengths from the file `neutron_scattering_lengths.json`
103 for each of the supplied species. Scattering lengths are in units of fm.
105 The scattering lengths have been extracted from the following NIST
106 database: https://www.ncnr.nist.gov/resources/n-lengths/list.html,
107 which in turn have been extracted from
108 Neutron News **3**, No. 3***, 29 (1992).
109 """
110 data_file = files(__package__) / 'neutron_scattering_lengths.json'
111 with open(data_file) as fp:
112 scattering_lengths = json.load(fp)
114 data = []
115 for species in scattering_lengths:
116 for isotope in scattering_lengths[species]:
117 for fld in 'b_c b_i'.split():
118 val = scattering_lengths[species][isotope][fld]
119 if 'None' in val:
120 scattering_lengths[species][isotope][fld] = np.nan
121 elif 'j' in val:
122 scattering_lengths[species][isotope][fld] = complex(val)
123 else:
124 scattering_lengths[species][isotope][fld] = float(val)
125 data.append(
126 dict(
127 species=species,
128 isotope=int(isotope),
129 abundance=scattering_lengths[species][isotope]['abundance']
130 / 100, # % -> fraction
131 b_coh=scattering_lengths[species][isotope]['b_c'],
132 b_inc=scattering_lengths[species][isotope]['b_i'],
133 )
134 )
135 scattering_lengths = DataFrame.from_dict(data)
136 return scattering_lengths