Bug fixes.

pyscses/set_up_calculation.py

@@ -9,8 +9,8 @@
 def site_from_input_file(site,
                          defect_species,
                          site_charge,
-                         core,
-                         temperature):
+                         temperature,
+                         zero_energies_less_than_kT=False):
     """
     Takes the data from the input file and converts it into a site.
     The input data file is a .txt file where each line in the file corresponds to a site. The values in each line are formatted and separated into the corresponding properties before creating a Site object for each site.
@@ -19,36 +19,32 @@ def site_from_input_file(site,
         site (str): A line in the input file.
         defect_species (object): Class object containing information about the defect species present in the system.
         site_charge (bool): The site charge refers to the contribution to the overall charge of a site given by the original, non-defective species present at that site. True if the site charge contribution is to be included in the calculation, False if it is not to be included.
+        temperature (float): Temperature in Kelvin.
+        zero_energies_less_than_kT (optional(bool)): If True segregation energies with absolute values < kT will be set to zero.
+            Default is False.
 
     Returns:
         :obj:`Site`
 
     """
-
     label = site[0]
     if site_charge == True:
         valence = float(site[1])
     if site_charge == False:
         valence = 0.0
     x = float(site[2])
     defect_labels = site[3::2]
-    defect_energies = [ float(e) for e in site[4::2] ]
-    min_energy = min(defect_energies)
-    # TODO: **Not well-defined what this is *supposed* to do.**
-    # TODO: Intended to define the core as being a single plane versus multiple planes.
-    # TODO: For multiple planes consider everything with |E_seg| < kT.
-    # TODO: For single plane, what is the intended behaviour if we have >1 site with the
-    # TODO: **same** lowest segregation energy?
-    if core == 'single':
-        for d_e in defect_energies:
-            if d_e > min_energy:
-                d_e = 0.0
-    if core == 'multi-site':
+    defect_energies = [float(e) for e in site[4::2]]
+    if zero_energies_less_than_kT:
+        kT = boltzmann_eV * temperature
         for d_e in defect_energies:
-            if ( -boltzmann_eV * temperature) <= d_e <= ( boltzmann_eV * temperature ):
+            if abs(d_e) < kT:
                 d_e = 0.0
-    #defect_energies = [ 0.0 for e in site[4::2] ]
-    return Site( label, x, [ defect_species[l] for l in defect_labels ], defect_energies, valence=valence )
+    return Site(label,
+                x, 
+                [defect_species[l] for l in defect_labels],\
+                defect_energies,
+                valence=valence)
 
 def format_line(line,
                 site_charge,

pyscses/site.py

@@ -71,6 +71,9 @@ 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))
 
     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.
@@ -141,10 +144,37 @@ def probabilities(self,
             if defect.fixed:
                 probabilities.append(defect.mole_fraction)
             else:
-                probabilities.append(defect.mole_fraction * defect.boltzmann_factor(phi, temp) /
-                                     (1.0 + sum([d.mole_fraction * (d.boltzmann_factor(phi, temp) - 1.0)
-                                                for d in self.defects])))
+                numerator = defect.mole_fraction * defect.boltzmann_factor(phi, temp)
+                denominator = (1.0 + sum([d.mole_fraction * (d.boltzmann_factor(phi, temp) - 1.0)
+                                          for d in self.defects]))
+                probabilities.append(numerator / denominator)
         return probabilities
+
+    def probabilities_fixed(self:
+                            phi: float,
+                            temp: float) -> List[float]:
+        """Calculates the probabilities of this site being occupied by each defect species.
+
+        Args:
+            phi (float): Electrostatic potential at this site in Volts.
+            temp (float): Temperature in Kelvin.
+
+        Returns:
+            list(float): Probabilities of site occupation on a 1D grid.
+
+        """
+        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) 
+                            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]   
 
     def defect_valences(self) -> np.ndarray:
         """Returns an array of valences for each defect from self.defects """