What value to use for rmax when integrating over a solute?

[jeff231li]

Hi,

I played around with the benzene example and I'm a bit confused with the choice of rmax. It seems that the solvation free energy is sensitive to this variable.

• 3.0 = 3.21
• 4.0 = 5.29
• 5.0 = 2.74
• 6.0 = 2.70
• 7.0 = 3.31
• 8.0 = 1.90

Is there a general rule for choosing an appropriate value for rmax?

Thank you

[fwaibl]

Hi,

yes, you are right, the solvation free energy is very sensitive to the integration cutoff. This is a general difficulty in GIST analyses, but it is also not optimally shown in the benzene example at the moment.

All GIST quantities are computed relative to the corresponding value in bulk water. For the water-water energy, this is done via the reference energy eww_ref. For the entropy, the cpptraj output should already be relative to bulk, but there is often a tiny offset which can sum to large numbers when the integration region is large. This can be avoided by choosing a small rmax, but then we miss parts of the long-range solute-water interaction. What you see is the combination of both effects: if rmax is too small, we don't see the full interaction, and if it is too large, we add increasing amounts of noise from the bulk solution.

To answer your question:

• If you are only interested in relative values, or in mapping thermodynamic values to the solvent surface, you should use a small rmax to be more robust against noise - maybe 4.5 Angstrom in this example, to include the first solvation shell.
• If you want to obtain a total solvation free energy, you need to get rid of the bias in bulk, and use a large cutoff to include all interactions. You should try to have well converged GIST calculations (maybe use more frames than in this example). If that is not enough, you can subtract a reference entropy from your GIST output (in the same way as the Eww reference). Then you can use a large integration cutoff rmax without adding much noise from the bulk.

The error introduced by subtracting a reference entropy is surprisingly small. You can test this by running GIST calculations with slightly different reference densities (this should lead to non-zero bulk entropies), and then integrate the solvation entropy after referencing each calculation to the respective bulk value - you should get almost the same result!

There is another thing to consider when computing absolute solvation free energies: In https://doi.org/10.1021/acs.jctc.0c01185 it was shown that the higher-order entropy contributions, which are omitted in GIST, scale roughly with -0.4 times the first-order entropy. This means you can improve the agreement between GIST and TI (or other "exact" methods) by scaling the entropy with a factor of 0.6. (In that paper, they also use an entropy reference as described above).

I will update the benzene example to show how to obtain more robust free energy values, and how to validate the convergence with rmax.

Best regards, Franz

[jeff231li]

@fwaibl thank you very much for your reply, this is very informative. I'm interested in estimating the change in solvation free energy upon binding. Will a short cutoff work in this case?

Also, can gisttools process files generated with PME-GIST?

Thank you, Jeff

[fwaibl]

Hi.

For the change in solvation free energy upon binding, I suggest that you try to obtain well converged free energies (using a larger cutoff).

I updated the benzene example to show how to subtract a reference entropy, and how this improves the radial convergence.

We also plan to implement some improvements to the entropy calculation in cpptraj, that might allow us to obtain converged entropy integrals without subtracting a reference.

The old gisttools version did not work with PME-GIST. I uploaded a patch that is able to work with PME-GIST output files today. However, be aware that the meaning of the energy columns is slightly different. The PME_norm and PME_dens columns (which will be used as Eww_norm by gisttools) correspond to the old Eww plus HALF the old Esw columns. Therefore, you need to add the potential energy of the solute to the integral to get the solvation free energy (corresponding to half the old Esw integral). This is printed out by PME-GIST, but I am not sure whether they include solute self-interactions in this number. I would have to dig deeper to find that out.

Best regards, Franz

[jeff231li]

I updated the benzene example to show how to subtract a reference entropy, and how this improves the radial convergence. We also plan to implement some improvements to the entropy calculation in cpptraj, that might allow us to obtain converged entropy integrals without subtracting a reference.

Thank you very much for your great work and help!

The old gisttools version did not work with PME-GIST. I uploaded a patch that is able to work with PME-GIST output files today. However, be aware that the meaning of the energy columns is slightly different. The PME_norm and PME_dens columns (which will be used as Eww_norm by gisttools) correspond to the old Eww plus HALF the old Esw columns. Therefore, you need to add the potential energy of the solute to the integral to get the solvation free energy (corresponding to half the old Esw integral). This is printed out by PME-GIST, but I am not sure whether they include solute self-interactions in this number. I would have to dig deeper to find that out.

I see, so the change in enthalpy will be something like total_energy = gistfile['PME_dens'] + 0.5*gistfile['Esw_dens'], integrate the voxels up to a given radius, and then subtract this value by the solute's self-interaction term, right?

[fwaibl]

Esw is set to zero in PME_GIST, because it is more difficult to split the energy terms obtained by PME. Eww is also set to zero, but I substitute PME_dens for Eww_unref_dens in gisttools, so that I can use it to detect rho0, reference energy etc.

The total change in enthalpy would be something like gistfile.integrate_around('Eww', rmax=...) + solute_energy - solute_self_energy). Here, solute_energy is the "Ensemble total solute energy on the grid" reported by cpptraj with PME-GIST, but you would have to calculate the solute self-energy yourself. For example, you could do a short MD run of only the solute and compute the electrostatic + VdW energy. In my testing, I read it out of an old benzene TI that I had.

[jeff231li]

Thank you very much for your help. I am able to reproduce the hydration free energies from https://doi.org/10.1021/acs.jctc.0c01185 using gisttools. I obtained the solute_self_energy term using the energy command in CPPTRAJ.

[fwaibl]

@jeff231li it turns out I was using an old GitHub version of cpptraj when updating gisttools for PME. The current version properly splits Esw and Eww, and can be used the same way as non-PME GIST. Therefore, I will probably undo some of the changes to gisttools that I implemented back then, such that the Eww columns are not replaced for PME any more. I am just posting this here because I don't want to mess up your calculations the next time you update gisttools...