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.

5

Yes, I have used averages from multiple independent calculations for both methods. however; the error bars are far from overlapping so averaging over more trajectories would probably decrease the error bars but not the discrepancy.

6

This process is highly problematic. First removing the charges leads to hydrophobic LJ spheres, which behave opposite from hydrophilic water molecules. Removing both charge and LJ at the same time could lead to intermediate states that bind extremely strongly to ions, depending on the soft-core parameters. So my guess is that both processes have severe, but opposite, sampling issues.