Skip to content

Commit 8820ef0

Browse files
pavoljuhaspavoljuhas
committed
Add class SFAverage for scattering factors stats.
Calculate compositional average of scattering factors from a structure object or a dictionary-passed composition.
1 parent 09fa456 commit 8820ef0

3 files changed

Lines changed: 267 additions & 3 deletions

File tree

‎diffpy/srreal/scatteringfactortable.py‎

Lines changed: 5 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -14,21 +14,23 @@
1414
##############################################################################
1515

1616

17-
"""class ScatteringFactorTable -- scattering factors for atoms, ions and
18-
isotopes.
17+
"""\
18+
class ScatteringFactorTable -- scattering factors for atoms, ions and isotopes.
1919
"""
2020

2121

2222
# exported items, these also makes them show in pydoc.
2323
__all__ = ['ScatteringFactorTable',
24-
'SFTXray', 'SFTElectron', 'SFTNeutron', 'SFTElectronNumber']
24+
'SFTXray', 'SFTElectron', 'SFTNeutron', 'SFTElectronNumber',
25+
'SFAverage']
2526

2627
from diffpy.srreal.srreal_ext import ScatteringFactorTable
2728
from diffpy.srreal.srreal_ext import SFTXray
2829
from diffpy.srreal.srreal_ext import SFTElectron
2930
from diffpy.srreal.srreal_ext import SFTNeutron
3031
from diffpy.srreal.srreal_ext import SFTElectronNumber
3132
from diffpy.srreal.wraputils import _pickle_getstate, _pickle_setstate
33+
from diffpy.srreal.sfaverage import SFAverage
3234

3335
# Pickling Support -----------------------------------------------------------
3436

‎diffpy/srreal/sfaverage.py‎

Lines changed: 161 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,161 @@
1+
#!/usr/bin/env python
2+
##############################################################################
3+
#
4+
# diffpy.srreal Complex Modeling Initiative
5+
# (c) 2016 Brookhaven Science Associates,
6+
# Brookhaven National Laboratory.
7+
# All rights reserved.
8+
#
9+
# File coded by: Pavol Juhas
10+
#
11+
# See AUTHORS.txt for a list of people who contributed.
12+
# See LICENSE.txt for license information.
13+
#
14+
##############################################################################
15+
16+
17+
"""\
18+
class SFAverage -- calculate average scattering factors and other statistics.
19+
20+
This class calculates compositional average of scattering factors and squared
21+
factors from a structure object or from atom concentrations.
22+
23+
24+
Examples
25+
--------
26+
::
27+
28+
import numpy as np
29+
from diffpy.srreal.scatteringfactortable import SFTXray
30+
xtb = SFTXray()
31+
qa = np.linspace(0, 25)
32+
sfavg1 = SFAverage.fromComposition(xtb, {'Na' : 1, 'Cl' : 1}, qa)
33+
sfavg2 = SFAverage.fromComposition(xtb, [('Mn', 1), ('O', 2)], qa)
34+
35+
from diffpy.Structure import loadStructure
36+
from diffpy.srreal.pdfcalculator import DebyePDFCalculator
37+
dpdfc = DebyePDFCalculator()
38+
dpdfc(loadStructure('rutile.cif'))
39+
sfavg3 = SFAverage.fromStructure(dpdfc.scatteringfactortable,
40+
dpdfc.getStructure(), dpdfc.qgrid)
41+
"""
42+
43+
44+
# class SFAverage ------------------------------------------------------------
45+
46+
class SFAverage(object):
47+
"""\
48+
Calculate compositional statistics of atom scattering factors.
49+
50+
Compositional averages can be calculated for an array of Q-values.
51+
Results are stored in the class attributes.
52+
53+
Attributes
54+
----------
55+
f1sum : float or NumPy array
56+
sum of scattering factors from all atoms.
57+
f2sum : float or NumPy array
58+
sum of squared scattering factors from all atoms.
59+
count : float
60+
total number of atoms. Can be non-integer in case of
61+
fractional occupancies
62+
f1avg : float or NumPy array
63+
compositional average of scattering factors.
64+
f2avg : float or NumPy array
65+
compositional average of squared scattering factors.
66+
composition : dict
67+
dictionary of atom symbols and their total abundancies.
68+
"""
69+
70+
f1sum = 0
71+
f2sum = 0
72+
count = 0
73+
f1avg = 0
74+
f2avg = 0
75+
composition = None
76+
77+
@classmethod
78+
def fromStructure(cls, sftb, stru, q=0):
79+
"""\
80+
Calculate average scattering factors from a structure object.
81+
82+
Parameters
83+
----------
84+
sftb : ScatteringFactorTable
85+
The ScatteringFactorTable object for looking up the values.
86+
stru : diffpy Structure or pyobjcryst Crystal or StructureAdapter
87+
The structure object that stores the atom species and their
88+
occupancies. Can be any type with a registered conversion
89+
to the StructureAdapter class.
90+
q : float or NumPy array (optional)
91+
The Q value in inverse Angstroms for which to lookup
92+
the scattering factor values.
93+
94+
See also
95+
--------
96+
RegisterStructureAdapter : to add support for more structure types.
97+
98+
Returns
99+
-------
100+
SFAverage
101+
The calculated scattering factor averages.
102+
"""
103+
# a bit of duck-typing for faster handling of diffpy.Structure
104+
if hasattr(type(stru), 'composition'):
105+
composition = stru.composition
106+
if isinstance(composition, dict):
107+
return cls.fromComposition(sftb, composition, q)
108+
# otherwise let's convert to a known structure type
109+
from diffpy.srreal.structureadapter import createStructureAdapter
110+
adpt = createStructureAdapter(stru)
111+
composition = {}
112+
for i in range(adpt.countSites()):
113+
smbl = adpt.siteAtomType(i)
114+
cnt = adpt.siteOccupancy(i) * adpt.siteMultiplicity(i)
115+
composition[smbl] = composition.get(smbl, 0) + cnt
116+
return cls.fromComposition(sftb, composition, q)
117+
118+
119+
@classmethod
120+
def fromComposition(cls, sftb, composition, q=0):
121+
"""\
122+
Calculate average scattering factors from atom concentrations.
123+
124+
Parameters
125+
----------
126+
sftb : ScatteringFactorTable
127+
The ScatteringFactorTable object for looking up the values.
128+
composition : dictionary or a list of (symbol, amount) pairs.
129+
The chemical composition for evaluating the average. Atom
130+
symbols may repeat when it is a list of (symbol, amount) pairs.
131+
q : float or NumPy array (optional)
132+
The Q value in inverse Angstroms for which to lookup
133+
the scattering factor values.
134+
135+
Returns
136+
-------
137+
SFAverage
138+
The calculated scattering factor averages.
139+
"""
140+
sfa = cls()
141+
sfa.composition = {}
142+
if isinstance(composition, dict):
143+
sfa.composition.update(composition)
144+
else:
145+
for smbl, cnt in composition:
146+
if not smbl in sfa.composition:
147+
sfa.composition[smbl] = 0
148+
sfa.composition[smbl] += cnt
149+
sfa.f1sum = 0.0 * q
150+
sfa.f2sum = 0.0 * q
151+
for smbl, cnt in sfa.composition.items():
152+
sfq = sftb.lookup(smbl, q)
153+
sfa.f1sum += cnt * sfq
154+
sfa.f2sum += cnt * sfq**2
155+
sfa.count += cnt
156+
denom = sfa.count if sfa.count > 0 else 1
157+
sfa.f1avg = sfa.f1sum / denom
158+
sfa.f2avg = sfa.f2sum / denom
159+
return sfa
160+
161+
# End of class SFAverage

‎diffpy/srreal/tests/testsfaverage.py‎

Lines changed: 101 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,101 @@
1+
#!/usr/bin/env python
2+
3+
"""Unit tests for diffpy.srreal.sfaverage
4+
"""
5+
6+
7+
import unittest
8+
import numpy
9+
10+
from diffpy.srreal.sfaverage import SFAverage
11+
from diffpy.srreal.scatteringfactortable import ScatteringFactorTable
12+
from diffpy.srreal.tests.testutils import loadDiffPyStructure
13+
from diffpy.srreal.tests.testutils import TestCaseObjCrystOptional
14+
from diffpy.srreal.tests.testutils import loadObjCrystCrystal
15+
16+
17+
##############################################################################
18+
class TestSFAverage(unittest.TestCase):
19+
20+
def setUp(self):
21+
self.sftx = ScatteringFactorTable.createByType('X')
22+
self.sftn = ScatteringFactorTable.createByType('N')
23+
return
24+
25+
def tearDown(self):
26+
return
27+
28+
29+
def test_fromStructure_CdSe(self):
30+
"""check SFAverage.fromStructure() for CdSe
31+
"""
32+
cdse = loadDiffPyStructure('CdSe_cadmoselite.cif')
33+
sfavg = SFAverage.fromStructure(self.sftx, cdse)
34+
fcd = self.sftx.lookup('Cd')
35+
fse = self.sftx.lookup('Se')
36+
self.assertTrue(isinstance(sfavg.f1sum, float))
37+
self.assertAlmostEqual(0.5 * (fcd + fse), sfavg.f1avg)
38+
self.assertAlmostEqual(0.5 * (fcd**2 + fse**2), sfavg.f2avg)
39+
self.assertEqual(4, sfavg.count)
40+
self.assertEqual(cdse.composition, sfavg.composition)
41+
qa = numpy.arange(0, 25, 0.1)
42+
sfavg2 = SFAverage.fromStructure(self.sftx, cdse, qa)
43+
self.assertTrue(isinstance(sfavg2.f1sum, numpy.ndarray))
44+
self.assertNotEqual(sfavg2.f1sum[0], sfavg2.f1sum[-1])
45+
self.assertEqual(sfavg.f1sum, sfavg2.f1sum[0])
46+
self.assertEqual(sfavg.f2sum, sfavg2.f2sum[0])
47+
sfavg3 = SFAverage.fromStructure(self.sftn, cdse, qa)
48+
self.assertEqual(sfavg3.f1sum[0], sfavg3.f1sum[-1])
49+
self.assertRaises(TypeError, SFAverage.fromStructure,
50+
self.sftx, 'notastructure')
51+
return
52+
53+
54+
def test_fromComposition(self):
55+
"""check SFAverage.fromComposition()
56+
"""
57+
sfavg1 = SFAverage.fromComposition(self.sftx, {'Na' : 1, 'Cl' : 1})
58+
fm = ['Na', 0.25, 'Cl', 0.75, 'Cl', 0.25, 'Na', 0.5, 'Na', 0.25]
59+
smblcnts = zip(fm[0::2], fm[1::2])
60+
sfavg2 = SFAverage.fromComposition(self.sftx, smblcnts)
61+
self.assertEqual(sfavg1.f1sum, sfavg2.f1sum)
62+
self.assertEqual(sfavg1.f2sum, sfavg2.f2sum)
63+
self.assertEqual(sfavg1.count, sfavg2.count)
64+
sfempty = SFAverage.fromComposition(self.sftx, [])
65+
self.assertEqual(0, sfempty.f1avg)
66+
self.assertEqual(0, sfempty.f2avg)
67+
self.assertEqual(0, sfempty.count)
68+
return
69+
70+
# End of class TestSFAverage
71+
72+
##############################################################################
73+
class TestSFAverageObjCryst(TestCaseObjCrystOptional):
74+
75+
def setUp(self):
76+
self.sftx = ScatteringFactorTable.createByType('X')
77+
return
78+
79+
def tearDown(self):
80+
return
81+
82+
83+
def test_from_rutile(self):
84+
"""check SFAverage.fromStructure for pyobjcryst Crystal of rutile.
85+
"""
86+
rutile = loadObjCrystCrystal('rutile.cif')
87+
qa = numpy.arange(0, 25, 0.1)
88+
sfavg = SFAverage.fromStructure(self.sftx, rutile, qa)
89+
fti = self.sftx.lookup('Ti', qa)
90+
fo = self.sftx.lookup('O', qa)
91+
self.assertTrue(numpy.allclose((fti + 2 * fo) / 3.0, sfavg.f1avg))
92+
fti2, fo2 = fti**2, fo**2
93+
self.assertTrue(numpy.allclose((fti2 + 2 * fo2) / 3.0, sfavg.f2avg))
94+
self.assertEqual(6, sfavg.count)
95+
96+
# End of class TestSFAverageObjCryst
97+
98+
if __name__ == '__main__':
99+
unittest.main()
100+
101+
# End of file

0 commit comments

Comments
 (0)