Dear all:
I’m using Lammps-5jun19 version to calculate the hydration free energy of a water molecule with SPC/E model. I used the fix adapt/fep + compute fep to calculate the free energy. My input file is similar to the inputs in the example files. My simulation system consists of 2000 water, as for the water molecule which is selected to be eliminated, the oxygen and hydrogen atom types are defined differently from the other 1999 water molecules (type 1 and 2 for the solvent and type 3 and 4 for the one selected). I analysed the result with FEP or BAR methods, finally, i got a value about -56.8 kcal/mol, which is obviously higher than the value reported in other literatures (about -41.2 kj/mol). But i can’t find out the problems in my simulation, my input file is as below, does anyone knows how to calculate the hydration free energy of water molecule?
variable lambda equal ramp(1.0,0.0)
variable coeff_wat equal v_lambda
variable qeow equal -0.8476*(v_coeff_wat)
variable qehw equal 0.4238*(v_coeff_wat)
fix ADAPT all adapt/fep 400000 &
pair lj/cut/coul/long/soft lambda 12 34 v_coeff_wat &
atom charge 3 v_qeow &
atom charge 4 v_qehw
variable dlam_wat equal -0.02
variable dqeow equal -0.8476v_dlam_wat
variable dqehw equal 0.4238v_dlam_wat
compute FEP all fep ${TK} &
pair lj/cut/coul/long/soft lambda 12 34 v_dlam_wat &
atom charge 3 v_dqeow &
atom charge 4 v_dqehw &
volume yes tail yes
fix FEP all ave/time 10 20000 400000 c_FEP[1] c_FEP[2] c_FEP[3] file fep_10.lmp
fix bar1 all ave/time 1 1 10 c_FEP[1] c_FEP[2] file bar_10.lmp
variable dwat2 equal 0.02
variable deow2 equal -0.8476v_dwat2
variable dehw2 equal 0.4238v_dwat2
compute FEP2 all fep ${TK} &
pair lj/cut/coul/long/soft lambda 12 34 v_dwat2 &
atom charge 3 v_deow2 &
atom charge 4 v_dehw2 &
volume yes tail yes
fix FEP2 all ave/time 10 20000 400000 c_FEP2[1] c_FEP2[2] c_FEP2[3] file fep_01.lmp
fix bar2 all ave/time 1 1 10 c_FEP2[1] c_FEP2[2] file bar_01.lmp
run 20000000 start 0 stop 20000000
unfix ADAPT
run 400000