Using Agents to Parameterize Molecular Dynamics Biomolecules

Coarse-Grained Molecular Dynamics (CG-MD)¹ lets scientists understand the behavior of complex molecular systems through the relationship between the behavior of the molecular structure of the system and its end properties, adding a higher-resolution level of understanding that might otherwise be expensive to obtain in a wet lab setting. This is a promising alternative especially in fields that deal with large molecular structures, such as drug design and materials science.

CG-MD, however, still requires human expertise to set up new simulations at multiple stages of the pipeline: (1) parameterizing the individual molecules, (2) setting up the simulation, (3) running the simulation, and (4) interpreting the results.

This post investigates the utility of using agents to do (1) directly from experimental measurements. The core idea being that for unknown molecules, teams should be able to take a measurement of the new compound and use an agent to parameterize the molecules in CG-MD grounded in experimental data, ready for downstream use.

Optimization goal

To keep things computationally feasible, I chose the common biomolecule glycerol and decided to parameterize only the three Lennard-Jones terms² using two experimentally measured properties found in the scientific literature:

1. The density of glycerol at atmospheric pressure and a set temperature³: 1258.4 ± 0.86 kg/m³.
2. The enthalpy of vaporization⁴ of glycerol: 91 ± 1 kJ/mol.

Note: All other parameters (bond length, angle, bonded force constants, etc.) were left at their reported values⁵.

Since the intention here is to minimally demonstrate parameterization from experimentally validated data points, the optimization problem is fitting the starting molecule to match the experimental density and vaporization enthalpy.

This can be formally expressed by the following loss function:

where:

• $\rho$ is the simulated density of the system, and $\rho^*$ is the target (experimental) density.
• $\Delta H_{\mathrm{vap}}$ is the simulated enthalpy of vaporization, and $\Delta H_{\mathrm{vap}}^*$ is the target (experimental) enthalpy of vaporization.

Both terms are weighted equally which means that the agent is penalized equally for missing the correct density or the correct enthalpy of vaporization.

The simulator

For simplicity, I use a purpose-built three-bead engine which maps each glycerol onto three beads, one per C–OH group at 30.7 u apiece, and simulates 125 of them. It runs at constant temperature and pressure (NPT) with a Berendsen thermostat and barostat to derive density and vaporization enthalpy.

Note: A simple extension here is to swap out my custom simulator for GROMACS or LAMMPS.

The agentic framework

The agent framework should enable having multiple agents that can be compared against each other. I want to compare (1) whether they are able to converge during minimization of the loss, and (2) how efficiently and reliably they are able to do this.

This requires a shared surface from which I can fork various implementations:

Now, as you can already see, the agentic surface above does not imply LLMs. This is intentional, making sure simpler paradigms such as Bayesian Optimization (BO) can be benchmarked, building in sophistication up to using LLMs.

This is the agentic spectrum I decided to test:

• BO: A basic Gaussian Process (GP) optimizing over the scalar loss by suggesting the next parameterization.
• BO plus a feasibility model: The same GP as above, this time I'm fitting it only on runs that did not crash (i.e. blow up), and delegate to a second GP to predict when a run may or may not crash. The latter should improve BO's ability to navigate the feasibility threshold.
• BO with an LLM warm-start: I hypothesize that the LLM's prior can save significant exploration runs. So in this case, the LLM seeds BO before switching to standard Bayesian optimization.
• LLM (Claude Opus 5): Finally, we hand the driving wheel to an LLM that can see the target experimental data, and can propose new parameters based on the scalar loss.

The four agents can branch out from the above abstraction by swapping out the policy where needed:

Parameterizing the biomolecule

Running the portfolio of agents over many simulations and many random seeds⁶ gives me the behavior shown below, which shows the representative behavior across the settings. The plot below shows a few things in one animation to understand (1) the optimization landscape across bead size and well-depth, (2) the trajectory taken by each agent, and (3) overall performance of each agent as more simulations are accrued.

Notice that the simulation can blow up, meaning that the parameterization becomes unstable; this feasibility region is denoted by a gray background with a grid of "x" markers.

Four agents searching the loss landscape over well depth and bead size on seed 7. Lighter shading is closer to the optimum, marked with a star; the hatched region is where simulations blow up. Improving moves are joined by lines, other proposals are faint, crashes are crosses. The shared convergence trace below plots best loss against simulations spent, with the achievable floor as a dashed line.

BO versus LLM-based approaches differ clearly in how they traverse the optimization landscape. While BO needs to spend a large portion of its simulation budget building a prior, an LLM already implicitly carries a prior from its pretraining phase, resulting in a more intentional optimization process.

This is not to say that BO methods don't earn their keep! BO in tandem with the feasibility model—though needing ~30 simulations—eventually builds a fine prior over the optimization landscape, allowing it to close much of the gap to the LLM-based approach, albeit less effectively.

Visualizing the results

Since this is a CG-MD parameterization exercise, let's look at the three-bead system being parameterized by the different agents, and how the final parameterization behaves in an actual CG-MD simulation.

Looking at the three-bead system at its first, middle, and final proposals, it's actually the level of cohesion among beads that changes: bead size moves by under 20%, while well depth moves by up to 75% (i.e., the visible change is the hue rather than the radius of the bead).

One glycerol bead as each of the four agents proposed it, at the first, middle and final proposal. Radius is bead size sigma, hue is well depth epsilon, and the dashed ring is the interaction cutoff. Gauges below each bead show the measured density and enthalpy of vaporization against their experimental targets represented as a red vertical line.

Finally, simulating the fitted system over 0, 30, and 60 picoseconds shows a disordered, voided slice with beads exchanging neighbors; this is the behavior of a liquid, not a lattice. It is also visibly too mobile: beads travel 0.29–0.43 nm rms over those 60 ps, about one bead diameter, where a real glycerol molecule would cover some 0.02 nm in the same window⁷ — the fitted models diffuse one to two hundred times too fast. This is a limitation of this study, where I'm focusing on the agentic setup to parameterize a glycerol stub three-bead system.

A 0.50 nm slice through each of the four fitted liquids at 0, 30 and 60 picoseconds. This animation shows the behavior of the parameterized molecules when interacting with each other.

Conclusions

In this blogpost I went through a demonstration of using an agentic framework for parameterizing a stub glycerol molecule to match its experimental data. The aim of this experiment is to understand whether agentic frameworks can let us more speedily configure new molecules in CG-MD when all we have is an experimental target, in this case density and vaporization enthalpy of glycerol.

I found out that indeed both Bayesian Optimization and, more directly, LLMs are able to parameterize the glycerol stub when slotted into an agentic framework next to a CG-MD simulator. With Claude Opus 5 being the clear winner, though using BO + a feasibility model is a competitive approach, considering that it doesn't involve any API spending.

Future work

Finally, as I've noted along this write-up, two things need to be swapped out in order for the output of this agentic framework to be physically usable in CG-MD pipelines:

1. Swap the custom-built simulator for GROMACS or LAMMPS: this will let the agentic framework interface with standard simulation pipelines.
2. Swap the stub bead glycerol model for one resulting from an all-atom simulation.

Disclaimer: I used Claude Opus 5 to collect glycerol's experimental data from the academic literature as well as to build the purpose-built glycerol simulator.