From c19f486670b042603ead9bd49849cfa6ce472d1c Mon Sep 17 00:00:00 2001
From: Anthony Scemama <scemama@irsamc.ups-tlse.fr>
Date: Tue, 21 Mar 2023 17:31:53 +0100
Subject: [PATCH] Added qp_extract_cipsi_data.py

---
 scripts/qp_extract_cipsi_data.py | 54 ++++++++++++++++++++++++++++++++
 src/tools/save_natorb.irp.f      |  1 +
 2 files changed, 55 insertions(+)
 create mode 100755 scripts/qp_extract_cipsi_data.py

diff --git a/scripts/qp_extract_cipsi_data.py b/scripts/qp_extract_cipsi_data.py
new file mode 100755
index 00000000..8f0b1f3c
--- /dev/null
+++ b/scripts/qp_extract_cipsi_data.py
@@ -0,0 +1,54 @@
+#!/usr/bin/env python3
+
+import re
+import sys
+
+# Read output file
+with open(sys.argv[1], 'r') as file:
+  output = file.read()
+
+
+def extract_data(output):
+    lines = output.split("\n")
+    data = []
+
+    n_det = None
+    e = None
+    pt2 = None
+    err_pt2 = None
+    rpt2 = None
+    err_rpt2 = None
+    e_ex = None
+
+
+    reading = False
+    for iline, line in enumerate(lines):
+        if not reading and line.startswith(" N_det "):
+            n_det = int(re.search(r"N_det\s+=\s+(\d+)", line).group(1))
+            reading = True
+
+        if reading:
+            if line.startswith(" E "):
+                e = float(re.search(r"E\s+=\s+(-?\d+\.\d+)", line).group(1))
+            elif line.startswith(" PT2 "):
+                pt2 = float(re.search(r"PT2\s+=\s+(-?\d+\.\d+E?.\d*)", line).group(1))
+                err_pt2 = float(re.search(r"\+/-\s+(-?\d+\.\d+E?.\d*)", line).group(1))
+            elif line.startswith(" rPT2 "):
+                rpt2 = float(re.search(r"rPT2\s+=\s+(-?\d+\.\d+E?.\d*)", line).group(1))
+                err_rpt2 = float(re.search(r"\+/-\s+(-?\d+\.\d+E?.\d*)", line).group(1))
+            elif "minimum PT2 Extrapolated energy" in line:
+                e_ex_line = lines[iline+2]
+                e_ex = float(e_ex_line.split()[1])
+                reading = False
+
+                data.append((n_det, e, pt2, err_pt2, rpt2, err_rpt2, e_ex))
+                n_det = e = pt2 = err_pt2 = rpt2 = err_rpt2 = e_ex = None
+
+    return data
+
+data = extract_data(output)
+
+for item in data:
+    print(" ".join(str(x) for x in item))
+
+
diff --git a/src/tools/save_natorb.irp.f b/src/tools/save_natorb.irp.f
index f6331d13..b4e9f5dc 100644
--- a/src/tools/save_natorb.irp.f
+++ b/src/tools/save_natorb.irp.f
@@ -12,6 +12,7 @@ program save_natorb
 ! matrices of each state with the corresponding
 ! :option:`determinants state_average_weight`
   END_DOC
+  PROVIDE nucl_coord
   read_wf = .True.
   touch read_wf
   call save_natural_mos