Chemical potential calculations for water

User discussions
1

GROMACS version: GROMACS/2023.3
GROMACS modification: No
Here post your question

Hi,
I have been performing chemical potential calculations for water using GROMACS by partitioning the solvent in a simulation of water + monovalent salts (ex: KCl and NaCl). Hence the simulation box contains SOL, 1 SWM (special water molecule but with same itp file with tip4pew–changing SOL to SWM) and corresponding ions. I have followed the decoupling of the van der Waals and electrostatic interactions in two ways (i) first decoupling the electrostatics, then van der Waals interactions (using coul_lambdas and vdw_lambdas sets of lambda parameters), and (ii) uncoupling both at the same time (using one common set of lambda parameters controlled by fep_lambdas). However the two methods, provide drastically different results (with error bars far away from overlapping) for chemical potential as a function of salt concentration. The trend is similar but the potential goes down much faster using fep_lambdas (second method). The first method provides more reasonable results (compared to literature). I am using the same lambda spacing in both methods and have repeated the calculation multiple times but the outcome persists.

I was wondering if anyone has noted this before for water and whether there is a reasonable explanation to such conflicting results. Any feedback would be appreciated.

2

You don’t write with method you use to compute the free-energy differences. If you use thermodynamic integration, the quadrature errors can be extremely large with such large decoupling effects. You should use BAR.

3

Thank you for your reply. I am using gmx bar to compute the free energy difference from the dhdl.xvg files for each lambda milestone (discarding first 100 ps at each milestone).

4

But have you done multiple calculations for each method? I would think that the results can differ massively, depending on how close, and how long, to an ion the water molecule gets.