Hello everyone.
I hope you are all doing great.
Recently, I was doing a test check on FEP in LAMMPS. In doing so, I started from solvation of methane in TIP4P-ew, an available example on Github from Agilio Padua (compute_fep/examples/USER/fep/CH4hyd/fdti01/fdti01.lmp at 7adab86f0ec19539d6a526666e1452b9bceb3c63 · agiliopadua/compute_fep · GitHub)
The first thing that got my attention was how to handle the dummy atom in TIP4P when there is a need to use soft-core potential. In other words, if I want to have efficiency when it comes to simulating 4-point water models, it is better to use the implicit approach, so basically I use long-range corrections with the style of tip4p (lj/cut/tip4p/long) with a Kspace as pppm/tip4p.
Thus, when I use this pair style and combine it with a different pair style for the solute, let’s say lj/cut/coul/long, I get the WARNING as it says the computation of coulomb interactions could be inconsistent. I read several topics here and studied the manual as well. Axel a couple of times mentioned if tip4p style is used, it needs to be used for all, otherwise in the case of soft/non-soft mix pair styles, the explicit approach should be taken.
I did two simulations for methane (OPLS-AA) in water (TIP4P-ew) using LAMMPS and GROMACS, with identical system geometry and particle numbers, tho the time step was bigger in GROMACS (2 fs) than LAMMPS (1 fs).
when I do finite difference thermodynamic integration on the original output from github I get this value 2.377 kcal/mol. keep in mind, they are combining tip4 style with a non tip4p style as below:
pair_style hybrid &
lj/cut/coul/long 10.0 10.0 &
lj/cut/tip4p/long/soft 3 4 2 2 0.125 1 0.5 10.0 10.0 10.0
Like I said, this brings on the warning.
So, I did tip4p style for all in this job script, I did the same simulation with the same time spent as the original example from github, and I got around 2.67 kcal/mol. To have a solid result, I increased the equlibration time and FEP time so that I equilibrated it for 10 ns, got FEP averages for 10 ns. By doing so I made sure, there is no poor sampling, and I got this value 2.303 kcal/mol.
Also, in the case of the GROMACS, I got the following results from TI, BAR, and MBAR respectively as follows
2.308 kcal/mol, 2.305 kcal/mol, 2.309 kcal/mol
I did the simulation in GROMACS for 10 ns as equilibration, and 10 ns in each window state to get those derivatives.
Now, it seems reasonable results to me. But I am not sure if I did correct in LAMMPS when I was doing FEP, I mean tip4p style. Here are my LAMMPS input commands. Could you please give me a hint if that is the way or not because I did not understand quite well when Axel said it is gotta be all tip4p style.
data.lmp (139.1 KB)
fdti01.lmp (2.8 KB)
in.sim (1.9 KB)
Other than this, if even this is the way, I do want to do it using the explicit approach for TIP4P as well. I did search a lot around here, and did not figure out how one person can model TIP4P/2005f using the explicit approach. I found an example here (WaterModelComparison/NPT simulations at main · VisualizationAndModelProduction/WaterModelComparison · GitHub) where they were doing TIP4P/2005f as a rigid model in the NPT / NVT ensembles. What is the point in having TIP4P/2005f if I am going to treat it as a rigid model? Could you please direct me to the right resource?