Added Site.saturation_parameter.

Added the saturation parameter as described in https://doi.org/10.1016/S0167-2738(02)00484-8 and https://doi.org/10.1038/s43588-021-00041-y

pyscses/site.py

@@ -22,9 +22,21 @@ class Site:
         defects (list): List of DefectAtSite objects, containing the properties of all individual defects at the site.
         scaling (float): A scaling factor that can be applied in the charge calculation.
         valence (float): The charge of the defect present at the site (in atomic units).
+        saturation_parameter (float): Optional saturation parameter as described in
+            `Hendricks et al. Sol. Stat. Ionics (2002)`_
+            and `Swift et al. Nature Comp. Sci. (2021)`_.
+            Setting `saturation_parameter` < `1.0` sets some proportion of excluded sites that are
+            unavailable for occupation by any explicit defects.
+            Default value is `1.0`, i.e., 100% of sites may be occupied.
         defects (list): List of Defect_Species objects for all defects present at the site.
         sites (list): List containing all x coordinates and corresponding  defect segregation energies.
 
+    .. _Hendricks et al. Sol. Stat. Ionics (2002):
+       https://doi.org/10.1016/S0167-2738(02)00484-8
+
+    .. _Swift et al. Nature Comp. Sci. (2021):
+       https://doi.org/10.1038/s43588-021-00041-y
+
     """
 
     def __init__(self,
@@ -33,7 +45,8 @@ def __init__(self,
                  defect_species: List[DefectSpecies],
                  defect_energies: List[float],
                  scaling: Optional[np.ndarray] = None,
-                 valence: float = 0.0) -> None:
+                 valence: float = 0.0,
+                 saturation_parameter: float = 1.0) -> None:
         """Initialise a Site object.
 
         Args:
@@ -43,10 +56,21 @@ def __init__(self,
             defect_energies (list(float)): List of defect segregation energies for each defect species at this site.
             scaling (optional, list(float): Optional list of scaling factors for the net charge at this site. Default scaling for each defect species is 1.0.
             valence (optional, float): Optional formal valence for this site in the absence of any defects. Default is 0.0.
+            saturation_parameter (optional, float): Optional saturation parameter as described in
+                Hendricks et al. Sol. Stat. Ionics (2002) [#HendricksEtAl_SolStatIonics2002]_
+                and Swift et al. Nature Comp. Sci. (2021). [#SwiftEtAl_NatureCompSci2021]_.
+                A saturation parameter < 1.0 introduces some proportion of excluded sites that are
+                unavailable for occupation by any explicit defects. Default is 1.0.
 
         Raises:
             ValueError if the number of DefectSpecies != the number of defect segregation energies != the number of scaling factors (if passed).
 
+        .. [HendricksEtAl_SolStatIonics2002]:
+           https://doi.org/10.1016/S0167-2738(02)00484-8
+
+        .. [SwiftEtAl_NatureCompSci2021]:
+           https://doi.org/10.1038/s43588-021-00041-y
+
         """
         if len(defect_species) != len(defect_energies):
             raise ValueError("len(defect_species) must be equal to len(defect_energies)")
@@ -71,9 +95,10 @@ def __init__(self,
             self.scaling = np.ones_like(defect_energies, dtype=float)
         self.grid_point: Optional[GridPoint] = None
         self.valence = valence
-        self.fixed_defects = tuple(d for d in self.defects if d.fixed=True)
-        self.mobile_defects = tuple(d for d in self.defects if d.fixed=False)
-        self.alpha = 1.0 - sum((d.mole_fraction for d in fixed_defects))
+        self.saturation_parameter = saturation_parameter
+        self.fixed_defects = tuple(d for d in self.defects if d.fixed)
+        self.mobile_defects = tuple(d for d in self.defects if not d.fixed)
+        self.alpha = self.saturation_parameter - sum((d.mole_fraction for d in self.fixed_defects))
 
     def competing_defect_species(self) -> Dict[str, int]:
         """Returns a dictionary reporting the number of fixed and / or mobile defect species that can occupy this site.
@@ -125,8 +150,8 @@ def average_local_energy(self,
             return self.grid_point.average_site_energy(method)
         else:
             raise ValueError("TODO")
-        
-    def probabilities(self:
+
+    def probabilities(self,
                       phi: float,
                       temp: float) -> List[float]:
         """Calculates the probabilities of this site being occupied by each defect species.
@@ -141,16 +166,16 @@ def probabilities(self:
         """
         probabilities_dict = {}
         boltzmann_factors = {d.label: d.boltzmann_factor(phi, temp) for d in self.mobile_defects}
-        denominator = (self.alpha + 
-                       sum([d.mole_fraction * (boltzmann_factors[d.label] - 1.0) 
+        denominator = (self.alpha +
+                       sum([d.mole_fraction * (boltzmann_factors[d.label] - 1.0)
                             for d in self.mobile_defects]))
         for defect in self.defects:
             if defect.fixed:
                 probabilities_dict[defect.label] = defect.mole_fraction
             else:
                 numerator = self.alpha * defect.mole_fraction * boltzmann_factors[defect.label]
                 probabilities_dict[defect.label] = numerator / denominator
-        return [probabilities_dict[d.label] for d in self.defects]   
+        return [probabilities_dict[d.label] for d in self.defects]
 
     def defect_valences(self) -> np.ndarray:
         """Returns an array of valences for each defect from self.defects """

tests/test_site.py

@@ -22,11 +22,31 @@ def create_mock_defect_species(n):
         mock_defect_species.append(m)
     return mock_defect_species
 
+def create_mock_defects_at_site(n):
+    labels = ['A', 'B', 'C', 'D', 'E']
+    valence = [-2.0, -1.0, 0.0, 1.0, 2.0]
+    mole_fraction = [0.15, 0.25, 0.35, 0.45, 0.55]
+    mobility = [0.1, 0.2, 0.3, 0.4, 0.5]
+    energies = [-0.1, -0.2, -0.3, -0.4, -0.5]
+    mock_defects_at_site = []
+    for i in range(n):
+        m = Mock(spec=DefectAtSite)
+        m.label = labels.pop()
+        m.valence = valence.pop()
+        m.mole_fraction = mole_fraction.pop()
+        m.mobility = mobility.pop()
+        m.energy = energies.pop()
+        m.fixed = False
+        mock_defects_at_site.append(m)
+    return mock_defects_at_site
+
 class TestSiteInit(unittest.TestCase):
 
     def test_site_is_initialised(self):
         mock_defect_species = create_mock_defect_species(2)
+        mock_defects_at_site = create_mock_defects_at_site(2)
         with patch('pyscses.site.DefectAtSite', autospec=True) as mock_DefectAtSite:
+            mock_DefectAtSite.side_effect = mock_defects_at_site
             site = Site(label='A',
                         x=1.5,
                         defect_species=mock_defect_species,
@@ -37,27 +57,58 @@ def test_site_is_initialised(self):
         self.assertEqual(site.defect_energies, [-0.2, +0.2])
         np.testing.assert_equal(site.scaling, np.array([1.0, 1.0]))
         self.assertEqual(site.valence, 0.0)
+        self.assertEqual(site.saturation_parameter, 1.0)
+        self.assertEqual(site.fixed_defects, ())
+        self.assertEqual(site.mobile_defects, tuple(mock_defects_at_site))
+        self.assertEqual(site.alpha, 1.0)
 
     def test_site_is_initialised_with_optional_args(self):
         mock_defect_species = create_mock_defect_species(2)
         with patch('pyscses.site.DefectAtSite', autospec=True) as mock_DefectAtSite:
+            mock_DefectAtSite.side_effect = create_mock_defects_at_site(2)
             site = Site(label='B',
                         x=1.5,
                         defect_species=mock_defect_species,
                         defect_energies=[-0.2, +0.2],
                         scaling=[0.5, 0.4],
-                        valence=-2.0)
+                        valence=-2.0,
+                        saturation_parameter=0.1)
         self.assertEqual(site.label, 'B')
         self.assertEqual(site.x, 1.5)
         self.assertEqual(site.defect_species, mock_defect_species)
         self.assertEqual(site.defect_energies, [-0.2, +0.2])
         np.testing.assert_equal(site.scaling, np.array([0.5, 0.4]))
         self.assertEqual(site.valence, -2.0)
+        self.assertEqual(site.saturation_parameter, 0.1)
+        self.assertEqual(site.alpha, 0.1)
+
+    def test_site_init_with_mixed_mobile_and_fixed_defects(self):
+        mock_defect_species = create_mock_defect_species(3)
+        mock_defects_at_site = create_mock_defects_at_site(3)
+        mock_defects_at_site[0].fixed = False
+        mock_defects_at_site[0].mole_fraction = 0.4
+        mock_defects_at_site[1].fixed = True
+        mock_defects_at_site[1].mole_fraction = 0.3
+        mock_defects_at_site[2].fixed = True
+        mock_defects_at_site[2].mole_fraction = 0.2
+        with patch('pyscses.site.DefectAtSite', autospec=True) as mock_DefectAtSite:
+            mock_DefectAtSite.side_effect = mock_defects_at_site
+            site = Site(label='C',
+            x=1.5,
+            defect_species=mock_defect_species,
+            defect_energies=[-0.2, +0.2, 0.0])
+        self.assertEqual(site.fixed_defects, (mock_defects_at_site[1], mock_defects_at_site[2]))
+        self.assertEqual(site.mobile_defects[0], mock_defects_at_site[0])
+        self.assertEqual(site.alpha, 0.5)
+    # TODO: test that fixed defects are added to the site.fixed_defects tuple and mobile defects are added to the site.mobiile_defects tuple
+
+    # TODO: test that site.alpha is computed correctly for various combinations of saturation_parameter and fixed defect species.
 
     def test_site_init_data_check_1(self):
         """Checks that initialising a Site object raises a ValueError if n(defect_species) != n(defect_energies)"""
         mock_defect_species = create_mock_defect_species(1)
         with patch('pyscses.site.DefectAtSite', autospec=True) as mock_DefectAtSite:
+            # mock_DefectAtSite.side_effect = create_mock_defects_at_site(1)
             with self.assertRaises(ValueError):
                 site = Site(label='A',
                             x=1.5,
@@ -79,12 +130,13 @@ class TestSite(unittest.TestCase):
 
     def setUp(self):
         mock_defect_species = create_mock_defect_species(2)
+        mock_defects_at_site = create_mock_defects_at_site(2)
         with patch('pyscses.site.DefectAtSite', autospec=True) as mock_DefectAtSite:
+            mock_DefectAtSite.side_effect = mock_defects_at_site
             self.site = Site(label='A',
                         x=1.5,
                         defect_species=mock_defect_species,
                         defect_energies=[-0.2, +0.2])
-        self.site.defects = [Mock(spec=DefectAtSite), Mock(spec=DefectAtSite)]
 
 
     def test_defect_with_label(self):