Estimating Solvation Free Energies with Boltzmann Generators

Accurate calculations of solvation free energies remain a central challenge in molecular simulations, often requiring extensive sampling and numerous alchemical intermediates to ensure sufficient overlap between phase-space distributions of a solute in the gas phase and in solution. Here, we introduce a computational framework based on normalizing flows that directly maps solvent configurations between solutes of different sizes, and compare the accuracy and efficiency to conventional free energy estimates. For a Lennard-Jones solvent, we demonstrate that this approach yields acceptable accuracy in estimating free energy differences for challenging transformations, such as solute growth or increased solute-solute separation, which typically demand multiple intermediate simulation steps along the transformation. Analysis of radial distribution functions indicates that the flow generates physically meaningful solvent rearrangements, substantially enhancing configurational overlap between states in configuration space. These results suggest flow-based models as a promising alternative to traditional free energy estimation methods.

💡 Research Summary

Accurate estimation of solvation free energies is a long‑standing challenge in molecular simulation because the transition from the gas phase to solution involves large, collective rearrangements of solvent molecules. Traditional alchemical approaches such as free‑energy perturbation (FEP) or multistate Bennett acceptance ratio (MBAR) mitigate the poor phase‑space overlap between the end states by introducing a series of intermediate Hamiltonians. While effective, this strategy incurs a substantial computational cost and requires careful selection of an order parameter.

In this work, Schebék et al. propose a fundamentally different framework based on normalizing flows, specifically Boltzmann Generators (BGs), to learn an invertible mapping that directly transforms solvent configurations from one solute state to another. The authors focus on a simple Lennard‑Jones (LJ) fluid containing two fixed solute particles (A and B) and 98 solvent particles in a cubic box (N = 100). Two test cases are examined: (i) systematic variation of the LJ radius σ_B of solute B (2.0 Å → 4.5 Å) over a temperature range of 100 K–140 K, and (ii) variation of the distance between the two solutes while keeping σ_B fixed.

The BG architecture employed is a coupling flow, which partitions the coordinates of each solvent particle into two channels and updates one channel conditioned on the other. This yields a triangular Jacobian with analytically tractable determinant, enabling exact likelihood evaluation. Training minimizes the reverse Kullback–Leibler divergence D_KL(p_θ‖p_B) (energy‑based training) using configurations sampled from a single base state (σ_0 = 3.0 Å, T_0 = 100 K). After training, the learned mapping f_θ is used to generate transformed configurations for any target state. The free‑energy difference Δf_AB is then obtained via the targeted free‑energy perturbation (TFEP) estimator
Δf_AB = −log ⟨w(x)⟩_{x∼p_A},
with importance weight w(x) = exp