I am currently studying the interaction of methane in water solvents, specifically using the TIP4P and SPC models. My objective is to investigate the effect of each solvent on methane by calculating the free energy difference between a reference state (methane in the TIP4P water model) and a perturbed state (methane in the SPC water model). I aim to evaluate the energy difference as follows:
Energy(Methane + TIP4P interacting) - Energy(TIP4P with Frozen Methane)
Energy(Methane + SPC interacting) - Energy(SPC with Frozen Methane)
I successfully ran the simulation of methane in TIP4P and calculated the energy of each resulting snapshot using the SPC water model. However, I am encountering issues when attempting to simulate a system of frozen methane in the TIP4P water solvent for the same snapshots. My NPT simulations are failing with the following error message:
Step 1 Warning: pressure scaling more than 1%, mu: 4.95406e+23 4.95406e+23 4.95406e+23
Program gmx mdrun, VERSION 5.1.4
Source code file: /gromacs-5.1.4/src/gromacs/ewald/pme-redistribute.cpp, line: 276
Fatal error:
4 particles communicated to PME rank 4 are more than 2/3 times the cut-off out of the domain decomposition cell of their charge group in dimension x.
This usually means that your system is not well equilibrated.
For more information and tips for troubleshooting, please check the GROMACS website at http://www.gromacs.org/Documentation/Errors
Could anyone provide guidance on how to resolve this issue?
Thank you for your questions. Let me clarify my goal.
I intend to evaluate the solvation free energy difference using the Zwanzig equation (Free Energy Perturbation method), with methane in TIP4P water as the reference state and methane in SPC water as the perturbed state.
Initially, I calculated the free energies of methane in both solvents using thermodynamic integration (TI) and found the difference between these two free energies, which I consider my ground truth. I expect to obtain a similar (if not nearly identical) value using the free energy perturbation method.
However, using the free energy perturbation method yields a significantly different value. I suspect that this discrepancy arises from the contributions of each water model, since the energy of a system, such as potential energy, is composed of E(Solute-Solvent interaction), E(Solvent-Solvent interaction), and E(Solute-Solute interaction). Therefore, I aim to subtract E(Solvent-Solvent interaction) for each simulation (methane in TIP4P and methane in SPC) to isolate E(Solute-Solvent interaction), as E(Solute-Solute interaction) will cancel out in the Zwanzig formula. By doing this, I will focus solely on methane interacting with the TIP4P and SPC water models, thereby eliminating the discrepancy.
By “frozen methane,” I mean that the methane is static and not interacting with the solvent, allowing me to remove E(Solute-Solvent interaction) from the total energy.
You can not use free energy perturbation to move between water models. The geometry of the water models is different.
And considering energies only is meaningless. Solvation free-energy in water have very large enthalpy-entropy compensation and you can not decompose the entropy.
Thank you for your comment. I was thinking that I can use the free energy perturbation (FEP) method to move between water models. However, I am curious to know if it is possible to use the FEP method along with any other techniques of adjustments to move between water models. Thank you for your patience and for caring about my concerne.
You can add more steps, e.g. a TI step for the geometry alone, but I don’t see the point of that. You already have the difference you want from TI. Or am I misunderstanding something.
Note that TI has been superseded by the Bennett acceptance ratio method. This should always given more accurate results. In GROMACS you can use the gmx bar tool on the same output as used for TI. The AWH method is another option. There are tutorials for both at: https://tutorials.gromacs.org/