In explicitly solvated simulations, the time spent in _compute_replica_energies() can still be significant.

When cycling over different thermodynamic states to compute energies where the only the alchemical parameters are changing but the positions are fixed, it should be possible to achieve a significant speedup by assigning the alchemically-modified Force objects to a different force group and only computing the energies for the alchemical force group. This could be added to the static component (independent of alchemical parameters) to compute the total reduced potential.

More generally, this same optimization can be used frequently, especially in multidimensional umbrella sampling with SAMS where only a fast-to-compute umbrella potential is changing during a lot of energy evaluations. I wonder if we could incorporate some sort of automatic caching mechanism into ThermodynamicState.reduced_potential() that handles this under the hood, maybe with some hints that the alchemical parameters all appear in one force group.