Protocol for solvation free energy calculations in LAMMPS


I am trying to get solvation free energies of small molecules in LAMMPS to match up to those calculated in AMBER. Both programs are using thermodynamic integration with softcore potentials for both vdW and Coulombic interactions, and integrated via 12-point Gaussian quadrature. The Lennard-Jones, bond, angle, and dihedral parameters I am using are exactly the same. I am getting the right ∆Gsolv for methane, but slightly larger solutes (isobutane, methanol) are getting incorrect values (the errors are on the order of 1-2 kcal/mol off for values around -5.5 kcal/mol).

I want to clarify what settings I should be using for the free energy calculation. I currently have these inputs:

units real

boundary p p p

force field parameters

atom_style full

bond_style harmonic

angle_style harmonic

dihedral_style charmm

special_bonds lj 0.0 0.0 0.5 coul 0.0 0.0 0.8333

read_data data/data_fftool_mod.lmp

pair_style hybrid &

lj/cut/coul/long 9.0 9.0 &

lj/cut/tip4p/long/soft 3 4 3 4 0.125 1 0.5 10.0 9.0 9.0

pair_modify tail no

kspace_style pppm/tip4p 1.0e-4


pair_coeff 1 1 lj/cut/coul/long 0.10939999991572773 3.399669508450741

pair_coeff 2 2 lj/cut/coul/long 0.01570000002623629 2.649532787260222

pair_coeff 1 2 lj/cut/coul/long 0.04144369681313082 3.0246011476505674


pair_coeff 3 3 lj/cut/tip4p/long/soft 0.16275000000964407 3.164349302954738 1.0

pair_coeff 4 4 lj/cut/tip4p/long/soft 0.0 1.0 1.0

pair_coeff 3 4 lj/cut/tip4p/long/soft 0.0 1.0 1.0


pair_coeff 1 3 lj/cut/tip4p/long/soft 0.14918465883419652 3.1545217994348804 ${lambda}

pair_coeff 1 4 lj/cut/tip4p/long/soft 0.0 1.0 ${lambda}

pair_coeff 2 3 lj/cut/tip4p/long/soft 0.056515208119110004 2.8075835657262154 ${lambda}

pair_coeff 2 4 lj/cut/tip4p/long/soft 0.0 1.0 ${lambda}

… and I am running 12 simulations with the appropriate values of ${lambda}. The commands which compute dU/d(lambda) are:

compute cFEP all fep ${temp} pair lj/cut/tip4p/long/soft lambda 12 34 v_dlambda

fix fFEP all ave/time 25 ${nave} ${nsteps} c_cFEP[1] c_cFEP[2] v_dlambda v_lambda v_weight file ti_out.l${lambda}.lmp

Are these settings doing what I think they are? Am I using the tip4p & soft pair styles correctly? Any other suggestions for why my calculated ∆Gs are not matching those in AMBER? My understanding is that this should ONLY scale the solute-water vdW and Coulombic interactions, which should be all that is necessary to compute ∆Gsolv. I believe that my parameters exactly match those I am using in AMBER, because I ran a straight potential energy calculation of the same XYZ configuration, and the breakdowns of the energies are matching exactly.