Improve equilibration and md configuration, improve analysis step

src/QligFEP/INPUTS/eq1.inp

@@ -6,7 +6,7 @@ bath_coupling             0.2
 random_seed               SEED_VAR
 initial_temperature       1
 shake_solvent             on
-shake_hydrogens           off
+shake_hydrogens           on
 shake_solute              off
 lrf                       on
 separate_scaling          on

src/QligFEP/INPUTS/eq2.inp

@@ -5,7 +5,7 @@ temperature               50
 bath_coupling             STEPSIZE
 shake_hydrogens           STEPTOGGLE
 shake_solute              off
-shake_solvent             on
+shake_solvent             off
 lrf                       on
 separate_scaling          on
 

src/QligFEP/INPUTS/eq3.inp

@@ -5,7 +5,7 @@ temperature               150
 bath_coupling             STEPSIZE
 shake_hydrogens           STEPTOGGLE
 shake_solute              off
-shake_solvent             on
+shake_solvent             off
 lrf                       on
 separate_scaling          on
 

src/QligFEP/INPUTS/eq5.inp

@@ -4,7 +4,7 @@ stepsize                  STEPSIZE
 temperature               T_VAR
 bath_coupling             10.0
 shake_hydrogens           STEPTOGGLE
-shake_solute              STEPTOGGLE
+shake_solute              off
 shake_solvent             on
 lrf                       on
 separate_scaling          on

src/QligFEP/INPUTS/md_0000_1000.inp

@@ -4,7 +4,7 @@ stepsize                  STEPSIZE
 temperature               T_VAR
 bath_coupling             10
 shake_hydrogens           STEPTOGGLE
-shake_solute              STEPTOGGLE
+shake_solute              off
 shake_solvent             on
 lrf                       on
 separate_scaling          on

src/QligFEP/INPUTS/md_0500_0500.inp

@@ -4,7 +4,7 @@ stepsize                  STEPSIZE
 temperature               T_VAR
 bath_coupling             10
 shake_hydrogens           STEPTOGGLE
-shake_solute              STEPTOGGLE
+shake_solute              off
 shake_solvent             on
 lrf                       on
 separate_scaling          on

src/QligFEP/INPUTS/md_1000_0000.inp

@@ -4,7 +4,7 @@ stepsize                  STEPSIZE
 temperature               T_VAR
 bath_coupling             10
 shake_hydrogens           STEPTOGGLE
-shake_solute              STEPTOGGLE
+shake_solute              off
 shake_solvent             on
 lrf                       on
 separate_scaling          on

src/QligFEP/INPUTS/md_XXXX_XXXX.inp

@@ -4,7 +4,7 @@ stepsize                  STEPSIZE
 temperature               T_VAR
 bath_coupling             10
 shake_hydrogens           STEPTOGGLE
-shake_solute              STEPTOGGLE
+shake_solute              off
 shake_solvent             on
 lrf                       on
 separate_scaling          on

src/QligFEP/analyze_FEP.py

@@ -442,8 +442,18 @@ def load_experimental_data(self, exp_key: str):
             _to = self.data[self.system][fep]["to"]
 
             # search ligands in the edges
+            ddG = None
             for edge in self.mapping_json["edges"]:
                 if edge["from"] == _from and edge["to"] == _to:
+                    if exp_key not in edge:
+                        logger.error(
+                            f"Experimental key '{exp_key}' not found in edge {_from} -> {_to}. "
+                            f"Available keys in this edge: {', '.join(edge.keys())}. "
+                            f"Please check your mapping JSON file and ensure the experimental key is correct."
+                        )
+                        raise KeyError(
+                            f"Key '{exp_key}' not found in edge data. Check your mapping JSON file."
+                        )
                     ddG = edge[exp_key]
                     break
 
@@ -489,15 +499,39 @@ def populate_mapping_dictionary(self, method, output_file: str | Path | None = N
             with output_file.open("w") as f:
                 json.dump(self.mapping_json, f, indent=4)
 
-    @staticmethod
-    def prepare_df(json_dict, experimental_data: bool = True):
-        pref = "dg" if "dg_error" in json_dict["edges"][0] else "ddg"
+    def prepare_df(self, json_dict, experimental_data: bool = True):
         df = pd.DataFrame(json_dict["edges"])
         if experimental_data:
+            if self.exp_key is None:
+                logger.error(
+                    "No experimental key has been set. Call load_experimental_data() first "
+                    "or set experimental_data=False."
+                )
+                raise ValueError("exp_key is None - cannot prepare dataframe with experimental data")
+
+            # For custom keys, edges have "delta_{node key}", which needs to be passed by the user
+            # The keys we use for experimental values are "ddg_value" (edge) or "dg_value" (node)
+            expected_col = f"{self.exp_key}"
+            if expected_col not in df.columns:
+                if "dg_value" in df.columns:
+                    expected_col = "dg_value"
+                elif "ddg_value" in df.columns:
+                    expected_col = "ddg_value"
+                else:
+                    available_cols = [col for col in df.columns if "delta_" in col or "_value" in col or "dg" in col.lower()]
+                    logger.error(
+                        f"Expected experimental data column '{self.exp_key}' not found in edges data. "
+                        f"Available columns that might contain experimental data: {', '.join(available_cols) if available_cols else 'none'}. "
+                        f"Please check for the correct key with experimental value on your mapping JSON file."
+                    )
+                    raise KeyError(
+                        f"Column '{self.exp_key}' not found. Check your mapping JSON file has the correct experimental data."
+                    )
+
             df = (
                 df.assign(
-                    ddg_value=lambda x: x[pref + "_value"],
-                    residual=lambda x: x[pref + "_value"] - x["Q_ddG_avg"],
+                    ddg_value=lambda x: x[expected_col],
+                    residual=lambda x: x[expected_col] - x["Q_ddG_avg"],
                     residual_abs=lambda x: x["residual"].abs(),
                 )
                 .sort_values("residual_abs", ascending=False)
@@ -646,7 +680,7 @@ def result_to_latex(res, latexify_each=False):  # TODO: move this out of this me
         # set labels, make it square and add legend
         plt.title(
             f"{(target_name + ' ' if target_name is not None else '')}"
-            r"$\Delta\Delta \text{G}_{\text{BAR}}$ ($\mathrm{N}="
+            r"$\Delta\Delta \mathrm{G}_{\mathrm{BAR}}$ ($\mathrm{N}="
             f"{len(exp_values)}$)"
         )
         plt.xlabel("$\Delta\Delta G_{exp} (kcal/mol)$")  # noqa: W605
@@ -887,7 +921,7 @@ def main(args: argparse.Namespace):
         results_df = fep_reader.prepare_df(results_json)
         if fep_reader.ignored_edges:
             results_df = results_df.query("~fep_name.isin(@fep_reader.ignored_edges)").reset_index(drop=True)
-        fig, ax = fep_reader.create_ddG_plot(results_df=results_df)
+        fig, _ = fep_reader.create_ddG_plot(results_df=results_df)
         fig.savefig(f"{args.target}_ddG_plot.png", dpi=300, bbox_inches="tight")
     else:
         results_json = json.loads((Path.cwd() / results_file).read_text())