estimate capacitance

ec_interface/ec_results.py

@@ -167,12 +167,10 @@ class ECResults:
     def __init__(
         self,
         ec_parameters: ECParameters,
-        data: NDArray,
-        verbose: bool = False
+        data: NDArray
     ):
         assert data.shape[1] == 5
         self.ec_parameters = ec_parameters
-        self.verbose = verbose
 
         # gather data from directories
         self.data = data
@@ -184,15 +182,11 @@ def __init__(
 
     @classmethod
     def from_calculations(cls, ec_parameters: ECParameters, directory: pathlib.Path, verbose: bool = False):
-        return cls(ec_parameters, _extract_data_from_directories(ec_parameters, directory, verbose), verbose=verbose)
+        return cls(ec_parameters, _extract_data_from_directories(ec_parameters, directory, verbose))
 
     def __len__(self):
         return self.nelects.shape[0]
 
-    def _outverb(self, *args, **kwargs):
-        if self.verbose:
-            print(*args, **kwargs)
-
     def to_hdf5(self, path: pathlib.Path):
         """Save data in a HDF5 file"""
 
@@ -201,7 +195,7 @@ def to_hdf5(self, path: pathlib.Path):
             dset.attrs['version'] = 1
 
     @classmethod
-    def from_hdf5(cls, ec_parameters: ECParameters, path: pathlib.Path, verbose: bool = False):
+    def from_hdf5(cls, ec_parameters: ECParameters, path: pathlib.Path):
         with h5py.File(path) as f:
             if 'ec_results' not in f:
                 raise Exception('invalid h5 file: no `ec_results` dataset')
@@ -210,7 +204,7 @@ def from_hdf5(cls, ec_parameters: ECParameters, path: pathlib.Path, verbose: boo
             if 'version' not in dset.attrs or dset.attrs['version'] > 1:
                 raise Exception('unknown version for dataset, use a more recent version of this package!')
 
-            return cls(ec_parameters, dset[:], verbose)
+            return cls(ec_parameters, dset[:])
 
     def estimate_active_fraction(self, shift: float = .0) -> float:
         """Estimate the active fraction from estimates of the surface capacitance.
@@ -225,10 +219,6 @@ def estimate_active_fraction(self, shift: float = .0) -> float:
         fit_2 = numpy.polyfit(work_function, fee, 2)  # grand pot vs work function
         cap_2 = -fit_2[0] * 2
 
-        self._outverb('Capacitance [e/V] = {:.5f} (charge), {:.5f} (grand potential), active fraction: {:.4f}'.format(
-            cap_1, cap_2, cap_2 / cap_1
-        ))
-
         return cap_2 / cap_1
 
     def compute_fee_hbm(self, alpha: float, shift: float = .0):

ec_interface/scripts/compute_fee.py

@@ -15,7 +15,6 @@ def main():
     parser = argparse.ArgumentParser(description=__doc__)
     parser.add_argument('-i', '--h5', default='ec_results.h5', help='H5 file')
     parser.add_argument('-o', '--output', default=sys.stdout, type=argparse.FileType('w'))
-    parser.add_argument('-v', '--verbose', action='store_true', help='Verbose')
 
     g_analysis = parser.add_mutually_exclusive_group()
     g_analysis.add_argument('--pbm', action='store_true', help='Assume PBM approach')
@@ -35,7 +34,7 @@ def main():
     ec_parameters = get_ec_parameters(this_directory)
 
     # extract data
-    ec_results = ECResults.from_hdf5(ec_parameters, args.h5, verbose=args.verbose)
+    ec_results = ECResults.from_hdf5(ec_parameters, args.h5)
 
     # write results
     args.output.write(
@@ -56,17 +55,29 @@ def main():
 
     if args.pbm:
         args.output.write('Grand potential (PBM) [V]\n')
-        numpy.savetxt(args.output, ec_results.compute_fee_pbm(shift=args.shift), delimiter='\t')
+        results = ec_results.compute_fee_pbm(shift=args.shift)
+        numpy.savetxt(args.output, results, delimiter='\t')
     elif args.hbm:
         args.output.write('Grand potential (HBM, alpha={:.4f}) [V]\n'.format(args.hbm))
-        numpy.savetxt(args.output, ec_results.compute_fee_hbm(alpha=args.hbm, shift=args.shift), delimiter='\t')
+        results = ec_results.compute_fee_hbm(alpha=args.hbm, shift=args.shift)
+        numpy.savetxt(args.output, results, delimiter='\t')
     elif args.hbm_ideal:
         alpha = ec_results.estimate_active_fraction(shift=args.shift)
         args.output.write('Grand potential (HBM, alpha={:.4f}) [V]\n'.format(alpha))
-        numpy.savetxt(args.output, ec_results.compute_fee_hbm(alpha=alpha, shift=args.shift), delimiter='\t')
+        results = ec_results.compute_fee_hbm(alpha=alpha, shift=args.shift)
+        numpy.savetxt(args.output, results, delimiter='\t')
     else:
         args.output.write('Grand potential (HBM, WF=Fermi) [V]\n')
-        numpy.savetxt(args.output, ec_results.compute_fee_hbm_fermi(shift=args.shift), delimiter='\t')
+        results = ec_results.compute_fee_hbm_fermi(shift=args.shift)
+        numpy.savetxt(args.output, results, delimiter='\t')
+
+    # Estimate differential capacitance
+    fit_2 = numpy.polyfit(results[:, 1], results[:, 2], 2)  # grand pot vs work function
+    args.output.write(
+        '\n\n'
+        'Capacitance [e/V]\n'
+        '{:.3f}'.format(-fit_2[0] * 2)
+    )
 
     args.output.close()