Hello,
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
solute-solute
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
water-water
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
solute-water
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.
Thanks,
Michelle