[gmx-users] Free energy calculation on a frozen fragment

[Ignacio Fernández Galván]

Hi all,

I'm interested in calculating the free-energy change when a frozen part of the system changes. The change, in general, implies a change in coordinates and in non-bonded parameters.

As far as I know, there's no efficient way to calculate this currently in gromacs. The only way I see is adding the two endpoints of the frozen fragment, and progressively turn on and off the electrostatic and vdw interactions (that is, use a dual-topology approach). This can be calculated either with TI or with BAR.

But this is quite cumbersome, and it would be more efficient if one could change the frozen coordinates with lambda. Something similar can be done with position restraints, but these don't give a frozen enough geometry, and I'm not sure the results would be quite reliable.

When calculating dH/dlambda with a frozen fragment, the contributions from the changing charges and non-bonded parameters are correctly included, what is missing is the contribution from the changing coordinates, and I believe this can be obtained if I can get an output of the average force (in vector form) on each affected atom. Is there some way to get this information (without storing the forces for the whole system)? Assuming that forces on frozen atoms are calculated at all.

With the average force I can get the contribution to dH/dlambda, but this would only work with TI, for BAR I would need to be able to calculate the energy with different coordinates, and I'm afraid that's harder.

Any other ideas on how to calculate the free energy when the coordinates of a frozen fragment change?

Thanks.
Ignacio