Coverage for dynasor/post_processing/neutron_scattering_lengths.py: 100%

1import json 

2from importlib.resources import files 

3from typing import Optional 

4 

5import numpy as np 

6from pandas import DataFrame 

7from .weights import Weights 

8 

9 

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. 

15 

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>`_. 

20 

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'] = {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 """ 

36 

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() 

43 

44 # Sub select only the relevant species 

45 scat_lengths = scat_lengths[scat_lengths.species.isin(atom_types)].reset_index() 

46 

47 for species in atom_types: 

48 if not np.any(scat_lengths.species == species): 

49 raise ValueError('Missing tabulated values ' 

50 f'for requested species {species}.') 

51 

52 # Update the abundances if another weighting is desired 

53 if abundances is not None: 

54 for species in abundances: 

55 scat_lengths.loc[scat_lengths.species == species, 'abundance'] = 0 

56 for Z, frac in abundances[species].items(): 

57 if not np.isfinite(frac) or not (0.0 <= frac <= 1.0): 

58 raise ValueError(f'Abundance fraction for {species} isotope {Z} ' 

59 f'must be finite and between 0 and 1 (got {frac}).') 

60 match = (scat_lengths.species == species) & (scat_lengths.isotope == Z) 

61 if not np.any(match): # Check if any row+column matches 

62 raise ValueError(f'No match in database for {species} and isotope {Z}') 

63 scat_lengths.loc[match, 'abundance'] = frac 

64 

65 self._scattering_lengths = scat_lengths 

66 

67 # Check if any of the fetched scattering lengths is None, 

68 # indicating that it is missing in the experimental database. 

69 # Only raise an error if the desired abundance is greater than 0. 

70 nan_rows = scat_lengths[scat_lengths.isnull().any(axis=1) & (scat_lengths.abundance > 0.0)] 

71 if not nan_rows.empty: 

72 # Grab first offending entry 

73 row = nan_rows.iloc[0] 

74 raise ValueError(f'{row.isotope}{row.species} is missing tabulated values for either' 

75 ' the coherent or incoherent scattering length.' 

76 ' Adjust the abundance parameter to set the fraction of' 

77 f' {row.isotope}{row.species} to zero.') 

78 

79 # Make sure abundances add up to 100% 

80 by_species = self._scattering_lengths.groupby('species') 

81 for species, species_df in by_species: 

82 if not np.isclose(species_df.abundance.sum(), 1): 

83 raise ValueError(f'Abundance values for {species} do not sum up to 1.0') 

84 

85 # Compute scattering lengths weighted by abundance 

86 weights_coh = by_species.apply( 

87 lambda s: (s.b_coh * s.abundance).sum(), 

88 include_groups=False 

89 ).to_dict() 

90 # First compute the average scattering length, then take the square 

91 # since the incoherent scattering lengths enter as b_incoh**2, but 

92 # dynasor only applies a single weighting factor w_incoh. 

93 weights_inc = by_species.apply( 

94 lambda s: (s.b_inc * s.abundance).sum(), 

95 include_groups=False 

96 ).to_dict() 

97 

98 supports_currents = False 

99 super().__init__(weights_coh, weights_inc, supports_currents) 

100 

101 @property 

102 def abundances(self) -> dict[str, dict[int, float]]: 

103 """Abundances used for calculating scattering lengths.""" 

104 abundance_dict = {} 

105 for (species), species_df in self._scattering_lengths.groupby('species'): 

106 abundance_dict[species] = {} 

107 for (isotope, abundance), _ in species_df.groupby(['isotope', 'abundance']): 

108 abundance_dict[species][isotope] = abundance 

109 return abundance_dict 

110 

111 @property 

112 def parameters(self) -> DataFrame: 

113 """Scattering lengths used to compute the coherent and 

114 incoherent weights for the selected isotopes. 

115 """ 

116 return self._scattering_lengths 

117 

118 

119def _read_scattering_lengths() -> DataFrame: 

120 """ 

121 Extracts the scattering lengths from the file `neutron-scattering-lengths-nist-1992.json` 

122 for each of the supplied species. Scattering lengths are in units of fm. 

123 

124 The scattering lengths have been extracted from the following NIST 

125 database: https://www.ncnr.nist.gov/resources/n-lengths/list.html, 

126 which in turn have been extracted from 

127 Neutron News **3**, No. 3***, 26 (1992). 

128 """ 

129 data_file = files(__package__) / 'form-factors/neutron-scattering-lengths-nist-1992.json' 

130 with open(data_file) as fp: 

131 scattering_lengths = json.load(fp) 

132 

133 data = [] 

134 for species in scattering_lengths: 

135 for isotope in scattering_lengths[species]: 

136 for fld in 'b_coh b_inc'.split(): 

137 val = scattering_lengths[species][isotope][fld] 

138 if 'None' in val: 

139 scattering_lengths[species][isotope][fld] = np.nan 

140 elif 'j' in val: 

141 scattering_lengths[species][isotope][fld] = complex(val) 

142 else: 

143 scattering_lengths[species][isotope][fld] = float(val) 

144 data.append( 

145 dict( 

146 species=species, 

147 isotope=int(isotope), 

148 abundance=float(scattering_lengths[species][isotope]['abundance']) 

149 / 100, # % -> fraction 

150 b_coh=complex(scattering_lengths[species][isotope]['b_coh']), 

151 b_inc=complex(scattering_lengths[species][isotope]['b_inc']), 

152 ) 

153 ) 

154 scattering_lengths = DataFrame.from_dict(data) 

155 return scattering_lengths