add reference

DOCUMENTATION.md

@@ -173,7 +173,8 @@ At the end of the procedure, a `results.csv` file (containing tab-separated data
 
 And then another dataset, with:
 + The charge added or removed to the system;
-+ The (absolute) work function, as `vacuum_potential - fermi_energy` (you might want to shift those value w.r.t. a reference such as the SHE);
++ The (absolute) work function, as `vacuum_potential - fermi_energy`;
++ The relative potential versus reference, `work_function - ref` (`ref=4.5V` by default, your can change this value with `--ref`);
 + The corresponding grand potential value. By default, it is computed as `free_energy - dn * fermi_energy`.
 
   **Note:** this might not be the correct free electrochemical energy you are looking for, and other methods are available: see [this document](white_paper/potential.pdf) for more information (TL;DR: use either `--pbm` or `--hbm xxx`, where `xxx` is the fraction of active electrons).

ec_interface/ec_results.py

@@ -212,7 +212,7 @@ def estimate_active_fraction(self, shift_with_avg: bool = False) -> float:
 
         return cap_2 / cap_1
 
-    def compute_fee_hbm(self, alpha: float, shift_with_avg: bool = False):
+    def compute_fee_hbm(self, alpha: float, shift_with_avg: bool = False, ref: float = 4.5):
         """Compute the Free electrochemical energy (grand potential) assuming a homogeneous background method
         calculation. `alpha` is the vacuum fraction.
         """
@@ -241,9 +241,9 @@ def compute_fee_hbm(self, alpha: float, shift_with_avg: bool = False):
         fee = fe0 + alpha * (self.free_energies - fe0 + dnelect * work_function - integ_average_pot)
 
         # get fee:
-        return numpy.array([dnelect, work_function, fee - shift_fee]).T
+        return numpy.array([dnelect, work_function, work_function - ref, fee - shift_fee]).T
 
-    def compute_fee_hbm_fermi(self, shift_with_avg: bool = False):
+    def compute_fee_hbm_fermi(self, shift_with_avg: bool = False, ref: float = 4.5):
         """Compute the Free electrochemical energy (grand potential) assuming a homogeneous background method
         calculation, and use the Fermi energy as the work function.
         """
@@ -257,9 +257,9 @@ def compute_fee_hbm_fermi(self, shift_with_avg: bool = False):
             index_0 = numpy.where(dnelect == .0)[0][0]
             shift_fee = self.average_potentials[index_0]
 
-        return numpy.array([dnelect, work_function, fee - shift_fee]).T
+        return numpy.array([dnelect, work_function, work_function - ref, fee - shift_fee]).T
 
-    def compute_fee_pbm(self, shift_with_avg: bool = False) -> NDArray:
+    def compute_fee_pbm(self, shift_with_avg: bool = False, ref: float = 4.5) -> NDArray:
         """Compute the Free electrochemical energy (grand potential) assuming a Poisson-Boltzmann method
         calculation"""
 
@@ -272,4 +272,4 @@ def compute_fee_pbm(self, shift_with_avg: bool = False) -> NDArray:
             index_0 = numpy.where(dnelect == .0)[0][0]
             shift_fee = self.average_potentials[index_0]
 
-        return numpy.array([dnelect, work_function, fee - shift_fee]).T
+        return numpy.array([dnelect, work_function, work_function - ref, fee - shift_fee]).T

ec_interface/scripts/compute_fee.py

@@ -15,6 +15,7 @@ def main():
     parser.add_argument('-p', '--parameters', default=INPUT_NAME, type=get_ec_parameters)
     parser.add_argument('-i', '--h5', default=H5_NAME, help='H5 file')
     parser.add_argument('-o', '--output', default=sys.stdout, type=argparse.FileType('w'))
+    parser.add_argument('-r', '--ref', default=4.5, type=float, help='Reference value')
 
     g_analysis = parser.add_mutually_exclusive_group()
     g_analysis.add_argument('--pbm', action='store_true', help='Assume PBM approach')
@@ -47,25 +48,26 @@ def main():
     args.output.write(
         '\n\n'  # just skip a few lines so that it is another dataset
         'Charge [e]\t'
-        'Work function [V]{}\t'.format(' (shifted with vacuum)' if args.shift_avg else '')
+        'Work function{} [V]\t'.format(' (shifted with vacuum)' if args.shift_avg else '') + \
+        'Potential vs ref [V]\t'
     )
 
     if args.pbm:
         args.output.write('Grand potential (PBM) [V]\n')
-        results = ec_results.compute_fee_pbm(shift_with_avg=args.shift_avg)
+        results = ec_results.compute_fee_pbm(shift_with_avg=args.shift_avg, ref=args.ref)
         numpy.savetxt(args.output, results, delimiter='\t')
     elif args.hbm:
         args.output.write('Grand potential (HBM, alpha={:.4f}) [V]\n'.format(args.hbm))
-        results = ec_results.compute_fee_hbm(alpha=args.hbm, shift_with_avg=args.shift_avg)
+        results = ec_results.compute_fee_hbm(alpha=args.hbm, shift_with_avg=args.shift_avg, ref=args.ref)
         numpy.savetxt(args.output, results, delimiter='\t')
     elif args.hbm_ideal:
         alpha = ec_results.estimate_active_fraction(shift_with_avg=args.shift_avg)
         args.output.write('Grand potential (HBM, alpha={:.4f}) [V]\n'.format(alpha))
-        results = ec_results.compute_fee_hbm(alpha=alpha, shift_with_avg=args.shift_avg)
+        results = ec_results.compute_fee_hbm(alpha=alpha, shift_with_avg=args.shift_avg, ref=args.ref)
         numpy.savetxt(args.output, results, delimiter='\t')
     else:
         args.output.write('Grand potential (HBM, WF=Fermi) [V]\n')
-        results = ec_results.compute_fee_hbm_fermi(shift_with_avg=args.shift_avg)
+        results = ec_results.compute_fee_hbm_fermi(shift_with_avg=args.shift_avg, ref=args.ref)
         numpy.savetxt(args.output, results, delimiter='\t')
 
     # Estimate differential capacitance

tests/test_ec_results.py

@@ -79,19 +79,19 @@ def test_compute_fee():
 
     # get similar results with fermi and "ideal" approach
     data_hbm_fermi = ec_results.compute_fee_hbm_fermi()
-    capacitance_hbm_fermi = -numpy.polyfit(data_hbm_fermi[:, 1], data_hbm_fermi[:, 2], 2)[0] * 2
+    capacitance_hbm_fermi = -numpy.polyfit(data_hbm_fermi[:, 1], data_hbm_fermi[:, 3], 2)[0] * 2
 
     assert capacitance_hbm_fermi == pytest.approx(0.05, abs=0.01)
 
     data_hbm_ideal = ec_results.compute_fee_hbm(alpha=0.334)
     assert numpy.allclose(data_hbm_ideal[:, 2], data_hbm_fermi[:, 2], rtol=1e-3)
 
-    capacitance_hbm_ideal = -numpy.polyfit(data_hbm_ideal[:, 1], data_hbm_ideal[:, 2], 2)[0] * 2
+    capacitance_hbm_ideal = -numpy.polyfit(data_hbm_ideal[:, 1], data_hbm_ideal[:, 3], 2)[0] * 2
     assert capacitance_hbm_ideal == pytest.approx(capacitance_hbm_fermi, abs=0.001)
 
     # get something different (and larger!) with the estimated vacuum fraction
     data_hbm_vac = ec_results.compute_fee_hbm(alpha=0.59)
-    capacitance_hbm_vac = -numpy.polyfit(data_hbm_vac[:, 1], data_hbm_vac[:, 2], 2)[0] * 2
+    capacitance_hbm_vac = -numpy.polyfit(data_hbm_vac[:, 1], data_hbm_vac[:, 3], 2)[0] * 2
 
     assert capacitance_hbm_vac == pytest.approx(0.09, abs=0.01)
     assert capacitance_hbm_vac != pytest.approx(capacitance_hbm_fermi, abs=0.001)