Dear Yicheng,
Regarding your initial question
However, I am uncertain whether this value includes the ideal gas energy difference between these two systems. In other words, if I wish to calculate the free energy difference between the NaCl+63H2O and KCl+63H2O systems, should the energy difference between the ideal gas corresponding to these two systems (which is equal to the free energy difference between gaseous Na and K in these two systems) be added to the thermodynamic integration result?
The method of thermodynamic integration with VCAIMAGES will give you the free energy difference between your two states aqueous NaCl and KCl. There is no need to add the free energy difference between the ideal gas states of NaCl and KCl. For your purpose the VCAIMAGES method in VASP is the appropriate method because it directly determines the free energy difference for your states of interest.
Another possibility would be to first determine the free energy for the aqueous NaCl system with respect to some reference system as the ideal gas or the Einstein crystal (harmonic solid). Both are valid reference systems. But in principle you could also use any other reference potential for which you can write down the free energy analytically. Next you would have to obtain the free energy for your aqueous KCl system with respect to the same reference potential and determine its free energy as well. Like this you would have obtained the free energies of NaCl and of KCl with respect to some artificial reference potential. The next step is now to compute the free energy difference between NaCl and KCl.
So to summarize, in the first method you directly obtain the free energy difference between your states. And in the second case you obtain the free energy difference by doing a detour by first computing the free energies of NaCl and KCl wrt to an artificial reference system and then compute the free energy difference. In theory these methods give the same answer because thermodynamic free energy differences only depend on the initial and the final state. They are independent of the path which was taken two reach from the initial to the final state. For more information I would recommend to take a look in the book Thermodynamics and an introduction to thermostatistics by Herbert Callen.
Because both methods are only exact in the adiabatic limit, you would have to make infinitely many images with the VCAIMAGES method. And you would have to varySCALEE for the reference system method infinitely slow. In practice this is not possible and therefore I would suspect the methods to give slightly different answers. But qualitatively the two approaches will yield the same result.
There is also a nice chapter in Daan Frenkels book Understanding Molecular Simulation about thermodynamic integration which I highly recommend. Also the references my colleague sent are very helpful. In your case, especially the second paper will be interesting.
If you are interested in how strong the two approaches deviate you are very welcome to try. It is definitely interesting to look into this. But there is no simple relation between the two methods because in principle they should give the same results.
I hope this is of help, all the best
Jonathan