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

@@ -3,9 +3,9 @@ steps                     NSTEPS2
 stepsize                  STEPSIZE
 temperature               50
 bath_coupling             STEPSIZE
+shake_solvent             on
 shake_hydrogens           STEPTOGGLE
 shake_solute              off
-shake_solvent             on
 lrf                       on
 separate_scaling          on
 

src/QligFEP/INPUTS/eq3.inp

@@ -3,9 +3,9 @@ steps                     NSTEPS2
 stepsize                  STEPSIZE
 temperature               150
 bath_coupling             STEPSIZE
+shake_solvent             on
 shake_hydrogens           STEPTOGGLE
 shake_solute              off
-shake_solvent             on
 lrf                       on
 separate_scaling          on
 

src/QligFEP/INPUTS/eq4.inp

@@ -3,9 +3,9 @@ steps                     NSTEPS2
 stepsize                  STEPSIZE
 temperature               275
 bath_coupling             STEPSIZE
+shake_solvent             on
 shake_hydrogens           STEPTOGGLE
 shake_solute              off
-shake_solvent             on
 lrf                       on
 separate_scaling          on
 

src/QligFEP/INPUTS/eq5.inp

@@ -3,9 +3,9 @@ steps                     NSTEPS1
 stepsize                  STEPSIZE
 temperature               T_VAR
 bath_coupling             10.0
-shake_hydrogens           STEPTOGGLE
-shake_solute              STEPTOGGLE
 shake_solvent             on
+shake_hydrogens           STEPTOGGLE
+shake_solute              off
 lrf                       on
 separate_scaling          on
 
@@ -26,7 +26,7 @@ polarisation              on
 polarisation_force        20.0
 
 [intervals]
-output                    5
+output                    25
 trajectory                100
 non_bond                  25
 

src/QligFEP/INPUTS/md_0000_1000.inp

@@ -3,9 +3,9 @@ steps                     NSTEPS2
 stepsize                  STEPSIZE
 temperature               T_VAR
 bath_coupling             10
-shake_hydrogens           STEPTOGGLE
-shake_solute              STEPTOGGLE
 shake_solvent             on
+shake_hydrogens           STEPTOGGLE
+shake_solute              off
 lrf                       on
 separate_scaling          on
 
@@ -26,7 +26,7 @@ polarisation              on
 polarisation_force        20.0
 
 [intervals]
-output                    5
+output                    25
 energy                    10
 trajectory                100
 non_bond                  25

src/QligFEP/INPUTS/md_0500_0500.inp

@@ -3,9 +3,9 @@ steps                     NSTEPS2
 stepsize                  STEPSIZE
 temperature               T_VAR
 bath_coupling             10
-shake_hydrogens           STEPTOGGLE
-shake_solute              STEPTOGGLE
 shake_solvent             on
+shake_hydrogens           STEPTOGGLE
+shake_solute              off
 lrf                       on
 separate_scaling          on
 
@@ -26,7 +26,7 @@ polarisation              on
 polarisation_force        20.0
 
 [intervals]
-output                    5
+output                    25
 energy                    10
 trajectory                100
 non_bond                  25

src/QligFEP/INPUTS/md_1000_0000.inp

@@ -3,9 +3,9 @@ steps                     NSTEPS2
 stepsize                  STEPSIZE
 temperature               T_VAR
 bath_coupling             10
-shake_hydrogens           STEPTOGGLE
-shake_solute              STEPTOGGLE
 shake_solvent             on
+shake_hydrogens           STEPTOGGLE
+shake_solute              off
 lrf                       on
 separate_scaling          on
 
@@ -26,7 +26,7 @@ polarisation              on
 polarisation_force        20.0
 
 [intervals]
-output                    5
+output                    25
 energy                    10
 trajectory                100
 non_bond                  25

src/QligFEP/INPUTS/md_XXXX_XXXX.inp

@@ -3,9 +3,9 @@ steps                     NSTEPS2
 stepsize                  STEPSIZE
 temperature               T_VAR
 bath_coupling             10
-shake_hydrogens           STEPTOGGLE
-shake_solute              STEPTOGGLE
 shake_solvent             on
+shake_hydrogens           STEPTOGGLE
+shake_solute              off
 lrf                       on
 separate_scaling          on
 
@@ -26,7 +26,7 @@ polarisation              on
 polarisation_force        20.0
 
 [intervals]
-output                    5
+output                    25
 energy                    10
 trajectory                100
 non_bond                  25

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
 
@@ -472,9 +482,9 @@ def populate_mapping_dictionary(self, method, output_file: str | Path | None = N
             for fep in self.feps:
                 fep_dict = self.data["result"][method][fep]
                 if fep_dict["from"] == _from and fep_dict["to"] == _to:
-                    avg_val = fep_dict[f"{method}_avg"] if not np.isnan(fep_dict[f"{method}_avg"]) else None
-                    sem_val = fep_dict[f"{method}_sem"] if not np.isnan(fep_dict[f"{method}_sem"]) else None
-                    std_val = fep_dict[f"{method}_std"] if not np.isnan(fep_dict[f"{method}_std"]) else None
+                    avg_val = fep_dict[f"{method}_avg"] if not pd.isna(fep_dict[f"{method}_avg"]) else None
+                    sem_val = fep_dict[f"{method}_sem"] if not pd.isna(fep_dict[f"{method}_sem"]) else None
+                    std_val = fep_dict[f"{method}_std"] if not pd.isna(fep_dict[f"{method}_std"]) else None
 
                     edge.update(
                         {
@@ -490,14 +500,34 @@ def populate_mapping_dictionary(self, method, output_file: str | Path | None = N
                 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(json_dict, experimental_data: bool = True, exp_key: Optional[str] = None) -> pd.DataFrame:
         df = pd.DataFrame(json_dict["edges"])
         if 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)
+            if exp_key is None:
+                if "dg_value" in df.columns:
+                    expected_col = "dg_value"
+                elif "ddg_value" in df.columns:
+                    expected_col = "ddg_value"
+            else:
+                expected_col = exp_key
+
+            if expected_col not in df.columns:
+                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 '{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 '{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 +676,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 +917,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())

src/QligFEP/pdb_utils.py

@@ -442,6 +442,17 @@ def write_dataframe_to_pdb(df, output_file, header: Optional[str] = None):
             file.write(pdb_line)
 
 
+def reindex_pdb_residues(pdb_path: Path, out_pdb_path: Path):
+    pdb_df = read_pdb_to_dataframe(pdb_path)
+    uniq_indexes = pdb_df.set_index(
+        ["residue_seq_number", "residue_name", "chain_id", "insertion_code"]
+    ).index
+    resn_mapping = {resn: idx for idx, resn in enumerate(uniq_indexes.unique(), 1)}
+    pdb_df["residue_seq_number"] = uniq_indexes.map(resn_mapping)
+    pdb_df["insertion_code"] = ""
+    write_dataframe_to_pdb(pdb_df, str(out_pdb_path.resolve().absolute()))
+
+
 def sdf_to_pdb(in_sdf_file, out_pdb_file):
     """Converts an SDF file to a PDB file.
 

src/QligFEP/qligfep.py

@@ -20,6 +20,7 @@
     pdb_parse_in,
     pdb_parse_out,
     read_pdb_to_dataframe,
+    reindex_pdb_residues,
     rm_HOH_clash_NN,
 )
 from .restraints.restraint_setter import RestraintSetter
@@ -1006,6 +1007,9 @@ def write_qprep(self, writedir):
         replacements["SOLUTEDENS"] = f"{density:.5f}"
 
         with open(qprep_in) as infile, open(qprep_out, "w") as outfile:
+            # We reindex the residues prior to defining the cysbonds because Q considers
+            # the first residue to be always 1, regardless of the numbering in the PDB file.
+            reindex_pdb_residues(Path(writedir) / self.pdb_fname, Path(writedir) / self.pdb_fname)
             cysbond_str = handle_cysbonds(
                 self.cysbond, Path(writedir) / self.pdb_fname, comment_out=(self.system != "protein")
             )