Add non-equilibrium (NEQ²) alchemical free-energy workflow

README.md

@@ -17,6 +17,8 @@ This repository is devoted to **QligFEP**, an automated workflow for small molec
   - [Compiling Q for local use (non-MPI)](#compiling-q-for-local-use-non-mpi)
   - [Setting up HPC configurations](#setting-up-hpc-configurations)
 - [⌨️ Command line interface (CLI)](#️-command-line-interface-cli)
+- [Tutorials](#tutorials)
+    - [Non-equilibrium FEP (NEQ²)](#non-equilibrium-fep-neq2)
 - [📊 Benchmarking](#-benchmarking)
 - [📚 Citations](#-citations)
 - [⏩ Q-GPU](#-q-gpu)
@@ -142,30 +144,48 @@ Now you're set with the qligfep package. This includes the command-linde-interfa
 5. `qlomap`: wraps `Lomap` to generate the `.json` perturbation mapping;
 6. `qmapfep`: in-house developed method to generate the `.json` perturbation mapping, interactively visualize and add or remove edges.
 7. `qligfep`: main CLI for running QligFEP simulations.
-8. `setupFEP`: sets up all the the QligFEP files for a simulation, including protein and water systems.
+8. `setupFEP`: sets up all the QligFEP files for a simulation, including protein and water systems. Pass `--neq` to set up the non-equilibrium (NEQ²) workflow instead of the windowed one.
 9. `qligfep_analyze`: CLI to analyze the results of a QligFEP simulation.
 10. `ligalign`: aligns a set of ligands to a reference ligand based on their maximum common substructure (MCS).
+11. `qligfep_neq_analyze`: CLI to analyze the results of a non-equilibrium (NEQ²) QligFEP simulation.
 
 ## Tutorials
 
-We are working on the documentation and tutorials for QligFEP. In the meantime, please refer to the Tyk2 case study available in the [tutorials directory](/tutorials/Tyk2/README.md). In addition to that, you can check the [benchmarking section](#-benchmarking) below, which contains the link to our benchmarking repository with scripts to reproduce the results.
+We are working on the documentation and tutorials for QligFEP. In the meantime, please refer to the Tyk2 case study available in the [tutorials directory](/tutorials/Tyk2/README.md). A dedicated [non-equilibrium (NEQ²) tutorial](/tutorials/Tyk2/neq2/README.md) walks through the NEQ² workflow end to end. In addition to that, you can check the [benchmarking section](#-benchmarking) below, which contains the link to our benchmarking repository with scripts to reproduce the results.
 
-# 📊 Benchmarking
+<a id="non-equilibrium-fep-neq2"></a>
 
-To check and reproduce QligFEP performance results, please refer to our [benchmarking repository](https://github.com/qusers/qligfepv2-BenchmarkExperiments).
+### Non-equilibrium FEP (NEQ²)
 
-For the preprint describing the benchmarking results, see:
+Alongside the standard windowed (equilibrium) protocol, QligFEP supports a **non-equilibrium** alchemical workflow, referred to as **NEQ²**. Rather than sampling many fixed-λ windows, NEQ² drives λ continuously from one end state to the other over many short, independent switching trajectories, and recovers ΔΔG from the Bennett Acceptance Ratio (BAR) over the resulting forward and reverse work distributions. Because the switching trajectories are independent, they parallelize trivially across a cluster.
 
-> Alencar Araripe D, Díaz Holguín A, Poso A, van Westen GJP, Åqvist J, Gutiérrez-de-Terán H, et al. Doing More with Less: Accurate and Scalable Ligand Free Energy Calculations by Focusing on the Binding Site. ChemRxiv. 2025; [doi:10.26434/chemrxiv-2025-x3r3z](https://doi.org/10.26434/chemrxiv-2025-x3r3z-v3)
+Set up a non-equilibrium calculation by passing `--neq` to `setupFEP`, and analyze the accumulated switching work with the `qligfep_neq_analyze` CLI. The non-equilibrium engine (`qdyn_neq`) is built together with the other Q binaries by `make all` in `src/q6`. See the [Tyk2 NEQ² tutorial](/tutorials/Tyk2/neq2/README.md) for an end-to-end walkthrough.
 
-# 📚 Citations
-Q6:       https://doi.org/10.1016/j.softx.2017.12.001
+# 📊 Benchmarking
 
-Q         https://doi.org/10.1016/S1093-3263(98)80006-5
+To check and reproduce QligFEP performance results, please refer to our [benchmarking repository](https://github.com/qusers/qligfepv2-BenchmarkExperiments).
 
-QligFEP:  https://doi.org/10.1186/s13321-019-0348-5
+# 📚 Citations
+
+To cite the latest version of QligFEP, cite:
+```bibtex
+@article{araripe2026qligfepv2,
+  author  = {Alencar Araripe, David and Díaz-Holguín, Alejandro and Poso, Antti and van Westen, Gerard J. P. and Åqvist, Johan and Gutiérrez-de-Terán, Hugo and Jespers, Willem},
+  title   = {Doing More with Less: Accurate and Scalable Ligand Free Energy Calculations by Focusing on the Binding Site},
+  journal = {Journal of Chemical Information and Modeling},
+  year    = {2026},
+  volume  = {66},
+  number  = {6},
+  pages   = {3164--3172},
+  doi     = {10.1021/acs.jcim.5c02932},
+  url     = {https://doi.org/10.1021/acs.jcim.5c02932},
+}
+```
+**Other relevant references:**
 
-QresFEP:  https://doi.org/10.1021/acs.jctc.9b00538
+- Q:        https://doi.org/10.1016/S1093-3263(98)80006-5
+- QligFEP:  https://doi.org/10.1186/s13321-019-0348-5
+- QresFEP:  https://doi.org/10.1021/acs.jctc.9b00538
 
 # ⏩ Q-GPU
 

pyproject.toml

@@ -32,14 +32,24 @@ dependencies = [
     "tqdm",
     "py3Dmol",
     "numpy",
+    "pandas",
     "matplotlib",
     "rdkit",
     "scikit-learn",
     "loguru",
-    "cinnabar",
     "MolClusterkit @ git+https://github.com/David-Araripe/MolClusterkit.git#master",
 ]
 
+[project.optional-dependencies]
+test = ["pytest"]
+
+[tool.pytest.ini_options]
+testpaths = ["test"]
+filterwarnings = [
+    # third-party (gufe/pydantic) deprecation surfaced when importing the package
+    "ignore:.*pydantic.config.Extra.*:DeprecationWarning",
+]
+
 [tool.setuptools]
 include-package-data = true
 
@@ -106,4 +116,5 @@ qmapfep = "QligFEP.CLI.qmapfep_cli:main_exe"
 qligfep = "QligFEP.CLI.qligfep_cli:main_exe"
 setupFEP = "QligFEP.CLI.setupFEP:main_exe"
 qligfep_analyze = "QligFEP.analyze_FEP:main_exe"
+qligfep_neq_analyze = "QligFEP.analyze_neq:main_exe"
 ligalign = "QligFEP.CLI.lig_align_cli:main_exe"

src/QligFEP/CLI/lomap_wrap_cli.py

@@ -19,9 +19,18 @@
 class LomapWrap:
     """Class to wrap the lomap package for the QligFEP CLI."""
 
-    def __init__(self, inp: str, out: Optional[str] = None, time=30, verbose="info", **kwargs):
+    def __init__(
+        self,
+        inp: str,
+        out: Optional[str] = None,
+        time=30,
+        verbose="info",
+        exp_key: Optional[str] = None,
+        **kwargs,
+    ):
         self.nodes = {}
         self.inp = inp
+        self.exp_key = exp_key
         self.out = self._parse_output(out)
         self.cores = self._setup_cores()
         self.lomap_args = {
@@ -158,17 +167,24 @@ def run_lomap(self) -> None:
         db_mol = lomap.DBMolecules(**self.lomap_args)
         rundir = db_mol.options["directory"]
 
-        extra_data_numerical = []
         for lomap_mol in db_mol._list:
             mol = lomap_mol.getMolecule()
             molpath = Path(rundir) / lomap_mol.getName()
             extra_data = self.extract_user_defined_properties(molpath)
-            extra_data_numerical.extend([k for k, v in extra_data.items() if isinstance(v, (int, float))])
             name = molpath.stem  # remove the file extension
             smiles = Chem.MolToSmiles(mol)
             formal_charge = Chem.GetFormalCharge(mol)
             self.nodes.update({name: {"smiles": smiles, "formal_charge": formal_charge, **extra_data}})
-        extra_data_numerical = list(set(extra_data_numerical))  # keep unique
+
+        # Rename experimental key to standardized dg_value on nodes
+        if self.exp_key:
+            for name, node in self.nodes.items():
+                if self.exp_key not in node:
+                    raise KeyError(
+                        f"exp_key '{self.exp_key}' not found in node '{name}'. "
+                        f"Available keys: {', '.join(k for k in node if k not in ('smiles', 'formal_charge'))}"
+                    )
+                node["dg_value"] = node.pop(self.exp_key)
 
         # Calculate the similarity matrices
         strict, loose = db_mol.build_matrices()
@@ -180,25 +196,11 @@ def run_lomap(self) -> None:
         for edge in result_dict["edges"]:
             _from = edge["from"]
             _to = edge["to"]
-            # check for the same formal charge
             same = self.nodes[_to]["formal_charge"] == self.nodes[_from]["formal_charge"]
-            for key in extra_data_numerical:
-                try:
-                    delta = self.nodes[_from][key] - self.nodes[_to][key]
-                    edge.update({f"delta_{key}": delta})
-                except TypeError:
-                    logger.warning(
-                        f"The {key} property is not numerical for one the ligands: {_from} | {_to}"
-                    )
-                    edge.update({f"delta_{key}": None})
-                except KeyError:
-                    logger.info(
-                        f"The {key} property is not present for one the ligands: {_from} | {_to} "
-                        "Information won't be kept in edge..."
-                    )
+            if self.exp_key:
+                edge["ddg_value"] = self.nodes[_to]["dg_value"] - self.nodes[_from]["dg_value"]
             edge.update({"same_charge": same})
             same_charges.append(same)
-            # update the potential ddG value
         if all(same_charges):
             logger.info("All ligands have the same formal charge.")
         else:
@@ -234,12 +236,19 @@ def parse_arguments() -> argparse.Namespace:
         "--time", "-t", type=int, default=30, help="Maximum time in seconds used to perform the MCS search."
     )
     parser.add_argument("--verbose", "-v", type=str, default="info", help="Verbosity level.")
+    parser.add_argument(
+        "--exp-key",
+        "-exp",
+        type=str,
+        default=None,
+        help="SDF property containing experimental dG. Stored as dg_value on nodes, ddg_value on edges.",
+    )
     return parser.parse_args()
 
 
 def main(args):
     # TODO: could implement other lomap arguemnts here
-    lomap = LomapWrap(args.input, args.output, args.time, args.verbose)
+    lomap = LomapWrap(args.input, args.output, args.time, args.verbose, exp_key=args.exp_key)
     lomap.run_lomap()
 
 

src/QligFEP/CLI/parser_base.py

@@ -227,6 +227,88 @@ def parse_arguments(program: str) -> argparse.Namespace:
             "Reproducible random state for the random FEP seed generator. Defaults to None (random FEP seeds)."
         ),
     )
+    parser.add_argument(
+        "-neq",
+        "--neq",
+        dest="neq",
+        action="store_true",
+        help=(
+            "Set up a non-equilibrium (NEQ) FEP instead of the windowed equilibrium approach. "
+            "Instead of many fixed-lambda windows, NEQ runs the `qdyn_neq` engine to drive lambda "
+            "from one endpoint to the other over a single simulation, accumulating the switching "
+            "work. Free energies are obtained from BAR over the forward/reverse work distributions "
+            "(see `qligfep_neq_analyze`). When set, the windowed parameters `--windows` and "
+            "`--sampling` are not used."
+        ),
+    )
+    parser.add_argument(
+        "-neqr",
+        "--neq-reps",
+        dest="neq_reps",
+        type=int,
+        default=5,
+        help=(
+            "NEQ only: number of forward/reverse switching pairs run per replicate. Each replicate "
+            "keeps a continuous endpoint equilibration and fires a forward and reverse switch from "
+            "successive snapshots to decorrelate the work samples. Defaults to 5."
+        ),
+    )
+    parser.add_argument(
+        "-neqs",
+        "--neq-steps",
+        dest="neq_steps",
+        type=int,
+        default=50000,
+        help=(
+            "NEQ only: length of each lambda-switching simulation in MD steps. Recommended >16000. "
+            "Defaults to 50000."
+        ),
+    )
+    parser.add_argument(
+        "-neqes",
+        "--neq-eq-steps",
+        dest="neq_eq_steps",
+        type=int,
+        default=1000,
+        help=(
+            "NEQ only: number of endpoint equilibration steps between successive switches "
+            "(spacing). Recommended >250. Defaults to 1000."
+        ),
+    )
+    parser.add_argument(
+        "-neqrs",
+        "--neq-relax-steps",
+        dest="neq_relax_steps",
+        type=int,
+        default=5000,
+        help=(
+            "NEQ only: length (MD steps) of the one-time endpoint relaxation run at lambda=0 "
+            "and lambda=1 before the first switch, to settle the nearly-decoupled ligand at each "
+            "endpoint. Applied to the first switching iteration; later iterations use the shorter "
+            "--neq-eq-steps (tEQ) spacing. ~10 ps at 2 fs. Defaults to 5000."
+        ),
+    )
+    parser.add_argument(
+        "-L",
+        "--neq-steepness",
+        dest="neq_L",
+        type=float,
+        default=8.0,
+        help=(
+            "NEQ only: steepness `L` of the sigmoidal lambda schedule l(t) = 1/[1+e^(L(t-0.5))]. "
+            "Higher L spends more time near lambda=0 and lambda=1; lower L approaches a linear "
+            "schedule. Recommended between 4 and 16. Defaults to 8. Ignored when "
+            "`--neq-schedule linear` is used."
+        ),
+    )
+    parser.add_argument(
+        "-neqsched",
+        "--neq-schedule",
+        dest="neq_schedule",
+        default="sigmoidal",
+        choices=["sigmoidal", "linear"],
+        help="NEQ only: lambda switching schedule. Defaults to `sigmoidal`.",
+    )
     parser.add_argument(
         "-log",
         "--log-level",

src/QligFEP/CLI/qligfep_cli.py

@@ -41,6 +41,13 @@ def main(args: Optional[argparse.Namespace] = None, **kwargs) -> None:
             "dr_force": args.dr_force,
             "random_state": args.random_state,
             "wath_ligand_only": args.wath_ligand_only,
+            "neq": args.neq,
+            "neq_reps": args.neq_reps,
+            "neq_steps": args.neq_steps,
+            "neq_eq_steps": args.neq_eq_steps,
+            "neq_relax_steps": args.neq_relax_steps,
+            "neq_L": args.neq_L,
+            "neq_schedule": args.neq_schedule,
         }
     else:
         param_dict = {}
@@ -65,6 +72,13 @@ def main(args: Optional[argparse.Namespace] = None, **kwargs) -> None:
                 command_str += f" --{k}".replace("_", "-")
         elif k == "dr_force":
             command_str += f" --{k} {v}".replace("dr_force", "distance_restraint_force")
+        elif k == "neq":
+            if v:
+                command_str += " --neq"
+        elif k == "neq_L":
+            command_str += f" --neq-steepness {v}"
+        elif k in ("neq_reps", "neq_steps", "neq_eq_steps", "neq_relax_steps", "neq_schedule"):
+            command_str += f" --{k.replace('_', '-')} {v}"
         else:
             command_str += f" --{k} {v}"
     command_str += f" --restraint_method {args.restraint_method}"
@@ -94,8 +108,9 @@ def main(args: Optional[argparse.Namespace] = None, **kwargs) -> None:
 
     run.write_water_pdb(inputdir)
 
-    logger.debug("Getting the lambdas")
-    lambdas = run.get_lambdas(args.windows, args.sampling)
+    if not run.neq:
+        logger.debug("Getting the lambdas")
+        lambdas = run.get_lambdas(args.windows, args.sampling)
     logger.debug("Writing atom mapping for distance restraints")
 
     run.avoid_water_protein_clashes(
@@ -109,6 +124,15 @@ def main(args: Optional[argparse.Namespace] = None, **kwargs) -> None:
     run.write_FEP_file(change_charges, change_vdw, FEP_vdw, inputdir, lig_size1, lig_size2)
     overlapping_atoms = run.set_restraints(writedir, args.restraint_method, strict_check=True)
 
+    if run.neq:
+        logger.debug("Writing the non-equilibrium MD files")
+        file_list = run.write_MD_neq(inputdir, lig_size1, lig_size2, overlapping_atoms)
+        run.write_neq_runfile(inputdir, file_list)
+        logger.debug(f"Generated files: {file_list}")
+        logger.debug("Writing the submit files")
+        run.write_submitfile(writedir)
+        return
+
     # Handling the correct offset here
     logger.debug("Writing the MD files")
     if args.start == "0.5":

src/QligFEP/CLI/setupFEP.py

@@ -43,6 +43,12 @@ def create_call(**kwargs):
         template += " -wath {water_thresh}"
     if "wath_ligand_only" in kwargs and kwargs["wath_ligand_only"]:
         template += " -wath-ligo"
+    if kwargs.get("neq"):
+        template += (
+            " --neq --neq-reps {neq_reps} --neq-steps {neq_steps} "
+            "--neq-eq-steps {neq_eq_steps} --neq-relax-steps {neq_relax_steps} "
+            "-L {neq_L} --neq-schedule {neq_schedule}"
+        )
     return template.format(**kwargs)
 
 
@@ -120,6 +126,13 @@ def main(args: Optional[argparse.Namespace] = None, **kwargs) -> None:
                 water_thresh=args.water_thresh,
                 log=args.log,
                 wath_ligand_only=args.wath_ligand_only,
+                neq=args.neq,
+                neq_reps=args.neq_reps,
+                neq_steps=args.neq_steps,
+                neq_eq_steps=args.neq_eq_steps,
+                neq_relax_steps=args.neq_relax_steps,
+                neq_L=args.neq_L,
+                neq_schedule=args.neq_schedule,
             )
             logger.info(f"Submitting the command:\n{command}")
             dst = sys_dir / f"FEP_{lig1}_{lig2}"