AIQM-PBSA

AIQM-PBSA is a hybrid endpoint binding free-energy framework for protein-ligand systems. It combines machine-learning interatomic-potential (MLIP) energies with classical MM/PBSA terms, using a full-system or ligand-centred hybrid treatment.

The framework is designed for AIQM/MLIP-MM calculations. The current initial release provides the MACE-OFF backend and provides the mace_S and mace_M models described below.

The repository provides two fitted models:

Method	MACE calculation	MM/PBSA inputs
mace_S	Complete complex, receptor, and ligand trajectories	Full-system MM/PBSA
mace_M	Hydrogenated ligand-centred QM-region trajectories	Full-system and QM-region MM/PBSA

Key ideas

• MLIP/MM hybrid treatment: an AIQM/MLIP backend supplies molecular energies while PBSA and MM/PBSA retain the classical solvation and energy terms.
• Initial MACE support: mace_S is the simpler full-structure workflow; mace_M replaces the local MM gas contribution with a QM-region MACE calculation (it is recommended to perform this procedure with neutral ligands, as the force field has not yet been trained in the charged systems).

Why AIQM-PBSA?

Conventional MM/PBSA is efficient but relies entirely on a molecular-mechanics description. High-level quantum calculations can improve local interactions but are often impractical for many trajectory frames. AIQM-PBSA is intended to use AIQM/MLIP energies where they are needed while retaining the efficiency of endpoint PBSA calculations. This initial implementation uses MACE-OFF24 for that energy component.

Method	Accuracy	Speed	Practical for screening
QM methods	⭐⭐⭐⭐	⭐	❌
MMPBSA	⭐⭐	⭐⭐⭐⭐	✅
AIQM-PBSA	⭐⭐⭐⭐	⭐⭐⭐	✅

Applications

• Protein-ligand binding-affinity estimation
• Lead optimization
• Retrospective virtual-screening analysis
• Comparative studies of full-structure and local-region MLIP corrections

Conceptual workflow

The figure illustrates the AIQM-PBSA hybrid endpoint-energy concept. The executable workflow below is the initial MACE-OFF24 implementation; molecular-dynamics preparation and production simulation are external prerequisites.

Requirements

• AmberTools or Amber with cpptraj, ante-MMPBSA.py, and MMPBSA.py (or MMPBSA.py.MPI).
• Python 3.10+ with the packages in requirements.txt.
• A local MACE-OFF checkpoint for the current MACE implementation.
• PyMOL with its Python module available when running mace_M hydrogen capping.

Install the Python dependencies in the environment used for the calculation:

pip install -r requirements.txt

MACE-OFF is a short-range organic force field. Check that all selected atoms use supported elements and that the protonation and charge state are appropriate for the model before a production calculation.

Inputs after MD

Start with a production trajectory and compatible Amber topology files for one receptor-ligand system.

• Full-system MM/PBSA requires solvated, complex, receptor, and ligand topologies, plus one or more trajectories with compatible atom ordering.
• The ligand mask must unambiguously select the ligand, for example :MOL.
• mace_S requires a solvent-free complex topology and a trajectory with exactly the same atoms.
• mace_M requires the solvated topology and trajectory used to define the ligand-centred QM region.
• The available fitted coefficient sets support only indi=1.0 and indi=2.1. Use the same value in the MMPBSA input file and in the final affinity command.

Quick start: mace_S

1. Run full-system MM/PBSA

bash scripts/run_mmpbsa.sh \
  --solvated-top <solvated.prmtop> --complex-top <complex.prmtop> \
  --receptor-top <receptor.prmtop> --ligand-top <ligand.prmtop> \
  --trajectory <production.mdcrd> --input config/mmpbsa_full.in \
  --output <results/full_mmpbsa.dat> --mpi-procs <N>

2. Export full-structure PDB trajectories

python python/prepare_full_mace_trajectories.py \
  --topology <solvent_free_complex.prmtop> --trajectory <dry_production.mdcrd> \
  --ligand-mask ':MOL' --start-frame <1> --end-frame <100> --stride <1> \
  --output-dir <work/full>

3. Calculate MACE energies

python python/calculate_mace_energies.py \
  --complex-pdb <work/full/full_complex.pdb> \
  --receptor-pdb <work/full/full_receptor.pdb> \
  --ligand-pdb <work/full/full_ligand.pdb> \
  --model <MACE-OFF24_medium.model> --device <cuda:0-or-cpu> \
  --output <results/mace_s_energy.json>

4. Calculate the affinity

python python/calculate_binding_affinity.py --method mace_S \
  --mmpbsa-result <results/full_mmpbsa.dat> \
  --mace-energy <results/mace_s_energy.json> --indi <1.0-or-2.1> \
  --output <results/mace_s_affinity.json>

Quick start: mace_M

First complete the full-system MM/PBSA calculation from the mace_S workflow.

1. Extract and hydrogenate the QM region

python python/prepare_qm_region.py \
  --topology <solvated.prmtop> --trajectory <production.mdcrd> \
  --ligand-mask ':MOL' --cutoff <6.0> \
  --start-frame <1> --end-frame <100> --stride <1> \
  --add-hydrogens --output-dir <work/qm>

This creates a stripped QM topology and trajectory. The receptor QM frames are hydrogenated and combined with the all-atom ligand frames to produce QM_complex_h.pdb.

2. Generate QM-region MMPBSA topologies

bash scripts/prepare_qm_mmpbsa_topologies.sh \
  --solvated-top <work/qm/QM_solvated_complex.prmtop> \
  --ligand-mask ':MOL' --output-dir <work/qm>

ante-MMPBSA.py generates matching complex, receptor, and ligand topologies for the QM-region trajectory. It prepares topology files; it does not calculate an energy.

3. Run QM-region MM/PBSA

bash scripts/run_qm_mmpbsa.sh \
  --solvated-top <work/qm/QM_solvated_complex.prmtop> \
  --complex-top <work/qm/QM_complex.prmtop> \
  --receptor-top <work/qm/QM_receptor.prmtop> \
  --ligand-top <work/qm/QM_ligand.prmtop> \
  --trajectory <work/qm/QM_trajectory.mdcrd> \
  --input config/mmpbsa_qm.in --output <results/qm_mmpbsa.dat> \
  --mpi-procs <N>

4. Calculate QM-region MACE energies and the final affinity

python python/calculate_mace_energies.py \
  --complex-pdb <work/qm/QM_complex_h.pdb> \
  --receptor-pdb <work/qm/QM_receptor_h.pdb> \
  --ligand-pdb <work/qm/QM_ligand.pdb> \
  --model <MACE-OFF24_medium.model> --device <cuda:0-or-cpu> \
  --output <results/mace_m_energy.json>

python python/calculate_binding_affinity.py --method mace_M \
  --mmpbsa-result <results/full_mmpbsa.dat> \
  --qm-mmpbsa-result <results/qm_mmpbsa.dat> \
  --mace-energy <results/mace_m_energy.json> --indi <1.0-or-2.1> \
  --output <results/mace_m_affinity.json>

Affinity model

For mace_M, the final affinity is:

c_mm * (Ggas_full_complex - Ggas_full_receptor - Ggas_QM_complex + Ggas_QM_receptor)
+ c_MACE * DeltaE_MACE
+ c_PB * EPB + c_SASA * ESURF + c_SAV * ENPOLAR

For mace_S, the MM replacement term is zero. The output JSON reports the final affinity in kcal/mol and every component of the formula.

Repository layout

config/   MMPBSA input templates
python/   QM extraction, MACE energy, and affinity programs
scripts/  AmberTools and MMPBSA command wrappers
images/   Documentation figures

Philosophy

Use molecular-learning potentials where they add value, and retain efficient endpoint solvation terms where they are sufficient.

Issues and pull requests are welcome.