From 575047c2787e457297813f54c5bdd43ad9322cc4 Mon Sep 17 00:00:00 2001
From: Evangelos Voyiatzis <evoyiatzis@gmail.com>
Date: Fri, 7 Jun 2024 14:40:21 +0200
Subject: [PATCH 1/2] Update fep.py

---
 tools/fep/fep.py | 47 ++++++++++++++++++++++++++++++++---------------
 1 file changed, 32 insertions(+), 15 deletions(-)

diff --git a/tools/fep/fep.py b/tools/fep/fep.py
index c7656c5ef5..0ce1927d73 100755
--- a/tools/fep/fep.py
+++ b/tools/fep/fep.py
@@ -2,23 +2,40 @@
 # fep.py - calculate free energy from compute fep results
 
 import sys
+from argparse import ArgumentParser
 import math
 
-if len(sys.argv) < 2:
-    print("Free Energy Perturbation")
-    print("usage: fep.py temperature < out.fep")
-    sys.exit()
+def compute_fep():
 
-rt = 0.008314 / 4.184 * float(sys.argv[1]) # in kcal/mol
+    parser = ArgumentParser(description='A python script to calculate free energy from compute fep results')
 
-v = 1.0
-sum = 0.0
-for line in sys.stdin:
-    if line.startswith("#"):
-        continue
-    tok = line.split()
-    if len(tok) == 4:
-        v = float(tok[3])
-    sum += math.log(float(tok[2]) / v)
+    parser.add_argument("units", help="unit system can be lj, real or si")
+    parser.add_argument("Temperature", type=float, help="The temperature of the system")
+    parser.add_argument("InputFile", help="The name of the file with the data from compute fep.")
+
+    args = parser.parse_args()
+
+    r_value = {'lj': 1.0, 'real': 0.0019872036, 'si': 8.314}
+
+    if args.units in r_value:
+        rt = r_value[args.units] * args.Temperature
+    else:
+        sys.exit("The provided units keyword is not valid")
+
+    v = 1.0
+    mysum = 0.0
+
+    with open(args.InputFile, "r") as input_file:
+        for line in input_file:
+            if line.startswith("#"):
+                continue
+            tok = line.split()
+            if len(tok) == 4:
+                v = float(tok[3])
+            mysum += math.log(float(tok[2]) / v)
+
+    print(-rt * mysum)
+
+if __name__ == "__main__":
+    compute_fep()
 
-print(-rt * sum)

From 567ba1f43752006c35e1f05828a3eacfc2ba1881 Mon Sep 17 00:00:00 2001
From: Evangelos Voyiatzis <evoyiatzis@gmail.com>
Date: Mon, 10 Jun 2024 19:37:35 +0200
Subject: [PATCH 2/2] improve R value for SI units

---
 tools/fep/fep.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/tools/fep/fep.py b/tools/fep/fep.py
index 0ce1927d73..cdcf3f10c6 100755
--- a/tools/fep/fep.py
+++ b/tools/fep/fep.py
@@ -15,7 +15,7 @@ def compute_fep():
 
     args = parser.parse_args()
 
-    r_value = {'lj': 1.0, 'real': 0.0019872036, 'si': 8.314}
+    r_value = {'lj': 1.0, 'real': 0.0019872036, 'si': 8.31446}
 
     if args.units in r_value:
         rt = r_value[args.units] * args.Temperature