I am considering trying a Bennett Acceptance Ratio (BAR) set of calculations using Gromacs. 

One of the complexed proteins (e.g., protein A interacting with protein B) one of them has a lysine and an arginine pair. Everything else is neutral. As the coul_lamdas and vdw_lambdas evolve, I will be left with a non-neutral system (2-). I am trying to deal with the charges in this system. 

I thought of coupling the two Cl ions to the vdw_lambdas and coul_lambdas of the protein, but surely this would bias the free energy estimation and when I compare this particular system to the wild-type (which has a neutral charge) the results will be inaccurate.  

However, if I run a separate series of lambda simulations on the two ions in bulk water, and also couple the ions as mentioned above, then I should be able to be the DeltaG_bind providing the following is correct.,

DG_bind = DG_complex(inc_ions) + DG_desolv + DG_ion_desolv

DG_solv = -DG_desolv

DG_ion_solv = -DG_ion_solv

Therefore, DG_bind = DG_complex(inc_ions) - DG_solv - DG_ion_solv 

Appreciate any feed back. Thanks