Merge pull request #129 from timcallow/MIS_calc

make a MIS property for the Density object

atoMEC/staticKS.py

@@ -620,6 +620,8 @@ def __init__(self, orbs):
             "N": np.zeros((config.spindims)),
         }
 
+        self._MIS = 0.0
+
         self._orbs = orbs
 
     @property
@@ -656,6 +658,13 @@ def unbound(self):
                 self._unbound = self.construct_rho_unbound()
         return self._unbound
 
+    @property
+    def MIS(self):
+        """ndarray: the mean ionization state."""
+        occs_pos = np.where(self._orbs.eigvals > 0, self._orbs.occnums_w, 0)
+        self._MIS = np.sum(occs_pos, axis=(0, 2, 3)) + self.unbound["N"]
+        return self._MIS
+
     @staticmethod
     def construct_rho_orbs(eigfuncs, occnums, xgrid):
         """

atoMEC/writeoutput.py

@@ -282,8 +282,7 @@ def write_final(self, energy, orbitals, density, conv_vals):
         output_str += self.write_final_energies(energy) + spc
 
         # write the chemical potential and mean ionization state
-        occs_pos = np.where(orbitals.eigvals > 0, orbitals.occnums_w, 0)
-        N_ub = np.sum(occs_pos, axis=(0, 2, 3)) + density.unbound["N"]
+        N_ub = density.MIS
 
         if config.spindims == 2:
             mu_str = "Chemical potential (u/d)"