GROMACS version: 5.0.7
GROMACS modification: No
Dear GROMACS users!
We’ve try to modify OPLS model of THF by changing its charges and wanted to estimate its vapor pressure. To do this we perform free energy calculation according well-known tutorial by Justin Lemkul. We’ve decoupled one THF molecule (solute) from THF bath (solvent). If I’m right a solvation free energy should be equal to RT ln(P_vap/cRT), where P_vap is a vapor pressure of our liquid at temperature T, and c is a particle concentration of decoupled molecules in a reference state which is ideal gas. That means that a result of free energy calculation should be size-dependent, by changing system size we change c in reference state. We’ve try to calculate it for 512 and for 4096 molecules and obtained the same results.
So, my questions are:
Why results are independent of system size while theoretically they should be?
What is the reference state we should consider?
How do pressure coupling work for the system, where one molecule is decoupled from other system? Do decoupled molecule and rest of liquid have same pressure (and hence different box volumes) or they have common box volume (what means we have in fact NVT ensemble for decoupled molecule)?
Is calculation of solvation free energy for THF molecule in liquid THF the best way to estimate vapor pressure?
First of all, you should really consider using a recent version of Gromacs. There have been performance upgrades and added features.
System size effects depend on the size and complexity of the molecule. THF is small and doesn’t have much conformational freedom so it doesn’t take a lot of molecules to generate a system that isn’t significantly perturbed by periodic boundary effects.
The reference state will be the ideal gas. You can treat the molecule as a single particle and the thermodynamics will still check out.
The pressure arises from intermolecular forces. By removing the forces of a molecule, you slightly decrease the pressure and the system will slightly shrink to compensate but this effect is negligible in normal conditions (i.e. the change in pV contribution to the free energy will be much smaller than the standard error you are probably getting for the whole free energy change).
Thank you for your reply!
Let me clarify my question. A reference state should be ideal gas, but what is the pressure of that ideal gas? Should we consider the pressure given by barostat (i.e. 1bar in my case), or we should recalculate it from the liquid volume in the state with lambda=1?
The first case should be true if pressure is coupled independently for removing molecule and the rest of the system. The second case occurs if solvent and removing molecule has common box and barostat hold the pressure for liquid. How do pressure coupling work in free energy calculation?
We use sd integrator.