Dear LAMMPS users and developers,
In order to estimate the free hydration energy of an ion (Na+) through the Thermodynamic Integration Method, I made use of the trajectories obtained at the full interaction potential (lambda=1) as the reference state, and successively changing the coupling parameter to 0.0, traced the ion-water interaction potential. The original system was made charge-neutral by insertion of a chloride anion into bulk water. The perturbation was accomplished by applying the soft pair potential and using the rerun command. It seems that the same approach is adopted by the “compute fep” command where the atom coordinates remain unchanged during such perturbation.
However, the reviewers of my paper have suggested that different series of simulations need to be conducted at each lambda value. I’m not sure if that is the right approach, because when the snapshots are not identical, it is not meaningful to compare the changes in the potential energy. On the other hand, at lambda values other than 1.0, my system won’t be charge-neutral any more. Could you please let me know about your opinion on this matter?
Many thanks in advance.
Best regards,
Monir Hosseini
1 Like
Did you just do runs at lambda = 0 and 1? Or also at in-between values?
The idea for using (that is, simulating at) lambda at in between values is so that you can more accurately integrate from the reference to the target state. However, if your trajectories between lambda=0 and lambda=1 overlap significantly, you can use Bennett’s acceptance ratio (BAR) and still get good results. However, just mentioning that to the reviewers might trigger a “then please show me that it’s the same” kind of response.
Technically it is fine to not have a neutral system when doing TI but it will make electrostatic calculations harder because you cannot use FFTs. You could technically scale down the anion charge or concentration with lambda to always remain neutral, but you have to ask yourself if that path in phase space still gives you the free energy difference that you’re interested in.
Hope this helps.
Hi Stefan,
I really appreciate your response.
Yes, lambda was lowered from 1 to 0 with a step size of 0.1, but I did not conduct a dynamic run at those values. The reference snapshots generated at lambda=1 were used, the coupling parameter was consecutively decreased and the time-average changes in the potential energy were monitored (by means of the rerun command).
So according to your comment, I can subject the system to a dynamic run at other lambda values, where the system is not really neutral? I’d rather not change the interaction strength of the charge-balancing anion, as it will definitly affect the free energy (of the cation) which I’m trying to calculate.
Thank you so much for your help.
Best regards,
Monir
Yes, you technically can do runs at where your charge is not neutral. However, I think you cannot use FFTs in that case which might make the computational time a lot longer. Someone with more experience with electrostatic might be able to help with that.
I would definitely recommend to at least do a run for lambda = 0 and then rerun those trajectories for lambda=1, since that will give you an estimate for the error in your free energy. This should give you the negative of the free energy difference you got from your original approach, if it converges reasonably.
1 Like
Sure, I will definitly try that.
Thanks a lot for your help.
Best regards,
Monir