From 0e35962ddf75e7cf84c848d793c3e778fed7b667 Mon Sep 17 00:00:00 2001 From: Pierre Beaujean Date: Wed, 23 Oct 2024 15:01:10 +0200 Subject: [PATCH] add reference --- DOCUMENTATION.md | 3 ++- ec_interface/ec_results.py | 12 ++++++------ ec_interface/scripts/compute_fee.py | 12 +++++++----- tests/test_ec_results.py | 6 +++--- 4 files changed, 18 insertions(+), 15 deletions(-) diff --git a/DOCUMENTATION.md b/DOCUMENTATION.md index 4b8ba66..e7e6a5a 100644 --- a/DOCUMENTATION.md +++ b/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). diff --git a/ec_interface/ec_results.py b/ec_interface/ec_results.py index 0a3e069..6e39260 100644 --- a/ec_interface/ec_results.py +++ b/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 diff --git a/ec_interface/scripts/compute_fee.py b/ec_interface/scripts/compute_fee.py index 739bf11..b4a1c53 100644 --- a/ec_interface/scripts/compute_fee.py +++ b/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 diff --git a/tests/test_ec_results.py b/tests/test_ec_results.py index 6b305ba..42f3fc5 100644 --- a/tests/test_ec_results.py +++ b/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)