Repository: bjmorgan/py-sc-fermi

Issue:
add a warning when defect concentrations are "beyond dilute limit"
It would be nice to warn the user if the code predicts defect concentrations beyond the dilute limit. This would be simple to implement, we would just need a definition of when the dilute limit approximation is becoming questionable. This could easily be heuristic.

Discussion:

Patch:
diff --git a/README.md b/README.md
index 483d32c..db3d6cd 100644
--- a/README.md
+++ b/README.md
@@ -5,28 +5,12 @@
 [![Documentation Status](https://readthedocs.org/projects/py-sc-fermi/badge/?version=latest)](https://py-sc-fermi.readthedocs.io/en/latest/?badge=latest)
 [![status](https://joss.theoj.org/papers/2fa8b3d80281c36371472485adbc83e7/status.svg)](https://joss.theoj.org/papers/2fa8b3d80281c36371472485adbc83e7)
 
-`py-sc-fermi` is a materials modelling code for calculating self-consistent Fermi energies and defect concentrations under thermodynamic equilibrium (or quasi-equilibrium) in the dilute limit given pre-calculated defect formation energies. For the theory, see [this paper](https://doi.org/10.1016/j.cpc.2019.06.017).   
+`py-sc-fermi` is a materials modelling code for calculating self-consistent Fermi energies and defect concentrations under thermodynamic equilibrium (or quasi-equilibrium) given defect formation energies. For the theory, see [this paper](https://doi.org/10.1016/j.cpc.2019.06.017).   
 
 The necessary inputs are (charged) defect formation energies, an (electronic) density of states, and the volume of the unit cell. Having this data, a `DefectSystem` object can be inititalised, from which the self consistent Fermi energy, defect concentrations, defect transition levels and electron and hole concentrations can be calculated. 
 
 Documentation and usage guides can be found [here](https://py-sc-fermi.readthedocs.io/en/latest/).
 
-## Installation
-
-py-sc-fermi can be installed with pip:
-
-```
-pip install py-sc-fermi
-```
-
-Alternatively, the latest development build can be installed by cloning this repository and installing with pip:
-
-```
-git clone https://github.com/bjmorgan/py-sc-fermi.git
-cd py-sc-fermi
-pip install .
-```
-
 ## Contributing
 
 ### Bugs reports and feature requests
diff --git a/docs/source/usage_notes.rst b/docs/source/usage_notes.rst
index 1af7889..12f12ab 100644
--- a/docs/source/usage_notes.rst
+++ b/docs/source/usage_notes.rst
@@ -4,8 +4,8 @@ Usage notes
   of states data was calculated, and any volumes, degeneracies and numbers of electrons provided should be 
   consitent with this structure for the reported defect concentrations to be accurate. In most cases, inconsistent
   specification of these values will lead to incorrect Fermi energies.
-- The reported Fermi energy and transition levels reported are referenced to 0 eV; the code expects that the input
-  density of states data is zeroed on the valence band maximum.
+- All Fermi energies reported by the code, and provided in inputs are referenced to the VBM, which in turn
+  must be referenced to 0 in the input density of states data.
 - The code operates using the following units,
   and all user-input must be consistent:
   
@@ -18,8 +18,3 @@ Usage notes
 - Concentrations are a special case, internally, the :py:mod:`py-sc-fermi` operates in the concentration of sites in the unit cell 
   which are defective, but will report in cm :superscript:`-3`. The documentation should always specify what kind of concentration
   is expected by a particular function, if it does not, this is a bug! Please report it. 
- 
-- Although in principle the code could be patched to account for temperature-dependent defect formation energy data
-  (if this is something you would be interested in, please let us know!), as this is a rare practice in the field, 
-  we currently only calculate defect concentrations with respect to the "0K" formation energies, following the model 
-  of `Zhang and Northrup`_https://doi.org/10.1103/PhysRevLett.67.2339
diff --git a/py_sc_fermi/defect_system.py b/py_sc_fermi/defect_system.py
index 8206caf..a0f2e86 100644
--- a/py_sc_fermi/defect_system.py
+++ b/py_sc_fermi/defect_system.py
@@ -149,7 +149,6 @@ class DefectSystem(object):
         direction = +1.0
         e_fermi = (emin + emax) / 2.0
         step = 1.0
-        converged = False
         reached_e_min = False
         reached_e_max = False
 
@@ -167,7 +166,6 @@ class DefectSystem(object):
                 reached_e_min = True
                 direction = +1.0
             if abs(q_tot) < self.convergence_tolerance:
-                converged = True
                 break
             if q_tot > 0.0:
                 if direction == +1.0:
@@ -192,7 +190,6 @@ class DefectSystem(object):
         """generate string to facilitate self.report()"""
         string = ""
         e_fermi = self.get_sc_fermi()[0]
-        string += f"Temperature :      {self.temperature}  (K)\n"
         string += f"SC Fermi level :      {e_fermi}  (eV)\n"
         p0, n0 = self.dos.carrier_concentrations(e_fermi, self.temperature)
         string += "Concentrations:\n"
@@ -201,7 +198,7 @@ class DefectSystem(object):
         for ds in self.defect_species:
             concall = ds.get_concentration(e_fermi, self.temperature)
             if ds.fixed_concentration == None:
-                string += f"{ds.name:9}      : {concall * 1e24 / self.volume} cm^-3\n"
+                string += f"{ds.name:9}      : {concall * 1e24 / self.volume} cm^-3, (percentage of defective sites: {(concall / ds.nsites) * 100:.3} %)\n"
             else:
                 string += (
                     f"{ds.name:9}      : {concall * 1e24 / self.volume} cm^-3 [fixed]\n"
@@ -336,3 +333,25 @@ class DefectSystem(object):
                 for ds in self.defect_species
             }
             return {**run_stats, **decomp_concs}
+
+    def site_percentages(
+        self, 
+    ) -> Dict[str, float]:
+        """Returns a dictionary of the DefectSpecies in the DefectSystem which
+        giving the percentage of the sites in the structure that will host that 
+        defect.
+
+        Returns:
+            Dict[str, Any]: dictionary specifying the per-DefectSpecies site
+            concentrations.
+        """
+
+        e_fermi = self.get_sc_fermi()[0]
+
+        sum_concs = {
+                str(ds.name): float(
+                    (ds.get_concentration(e_fermi, self.temperature) / ds.nsites) * 100
+                )
+                for ds in self.defect_species
+            }
+        return sum_concs


Test patch:
diff --git a/tests/dummy_inputs/report_string.txt b/tests/dummy_inputs/report_string.txt
index f25ed41..f6254db 100644
--- a/tests/dummy_inputs/report_string.txt
+++ b/tests/dummy_inputs/report_string.txt
@@ -1,4 +1,3 @@
-Temperature :      298  (K)
 SC Fermi level :      0.5  (eV)
 Concentrations:
 n (electrons)  : 1e+24 cm^-3
@@ -12,4 +11,4 @@ v_O        : Charge Concentration(cm^-3) Total
            :  1  1.000000e+25          100.00  [fixed]
 ---------------------------------------------------------
 O_i        : Charge Concentration(cm^-3) Total
-           : -1  1.000000e+25          100.00  [fixed]
+           : -1  1.000000e+25          100.00  [fixed]
\ No newline at end of file
diff --git a/tests/test_defect_system.py b/tests/test_defect_system.py
index 1903787..d68afb0 100644
--- a/tests/test_defect_system.py
+++ b/tests/test_defect_system.py
@@ -121,6 +121,19 @@ class TestDefectSystem(unittest.TestCase):
             {"Fermi Energy": 1, "p0": 1, "n0": 1, "O_i": 1, "v_O": 1},
         )
 
+    def test_site_percentages(self):
+        self.defect_system.get_sc_fermi = Mock(return_value=[1, {}])
+        self.defect_system.dos.carrier_concentrations = Mock(return_value=(1, 1))
+        self.defect_system.defect_species[0].get_concentration = Mock(return_value=1)
+        self.defect_system.defect_species[1].get_concentration = Mock(return_value=1)
+        self.defect_system.defect_species[0].nsites = 1
+        self.defect_system.defect_species[1].nsites = 1
+        self.defect_system.defect_species[0].name = "v_O"
+        self.defect_system.defect_species[1].name = "O_i"
+        self.assertEqual(
+            self.defect_system.site_percentages(), {"v_O": 100, "O_i": 100}
+        )
+
     def test__get_report_string(self):
         self.defect_system.get_sc_fermi = Mock(return_value=[0.5, {}])
         self.defect_system.dos.carrier_concentrations = Mock(return_value=(100, 100))