Protocol for solvation free energy calculations in LAMMPS

Hello,

From what I can see, your use of the soft TIP4P pair styles is correct.
And yes, only solute-solvent interactions should be scaled.
I’ve known of other uses of this package for solutes like ethanol and, upon comparison with other codes, users didn’t report discrepancies.

One thing you can try to do is a run 0 at a given lambda_i, using the compute fep to calculate the energy difference upon perturbation with dlambda. Atoms won’t move. Output the c_cFEP[1].
Then compare the “perturbed” energy with that of a run 0 at lambda_{i+1} = lambda_i+dlambda, using the same configuration.
Do they match? If so, then the compute fep is doing it’s job.

If you want I can do this test for you, just send me the input files.
Please keep me updated. I’d like to know if something needs fixing.

Agilio