I have been trying to finish the polymer tutorial for MARTINI, I was able to reproduce the bonded distribution, and I was also able to get the density, end to end distance and gyration radius within 5%. But I am not able to reproduce the solvation free energy in case of PEG9 in a solution. According to the tutorial I should get 7kJ/mol, but I am getting a value of 150kJ/mol. 

md.mdp file:

integrator = sd

; start time and timestep in ps =

tinit = 0

dt = 0.01

nsteps = 2500000

cutoff-scheme = group

; We remove center of mass motion. In periodic boundary conditions, the center of mass motion is spurious; the periodic system is the same in all translational directions.

comm-mode = Linear

; number of steps for center of mass motion removal =

nstcomm = 10;5 ; ORIGINALLY 10

ns_type = grid ; USER ADDED NOT IN TUTORIAL

pbc = xyz ; USER ADDED NOT IN TUTORIAL

table-extension = 2

; Output frequency for energies to log file and energy file =

nstxout = 0

nstvout = 0

nstfout = 0

nstlog = 10000

nstenergy = 10000

nstcalcenergy = 10

nstxtcout = 10000

nstlist = 10;10 ;5 ; ORIGINALLY 10

rlist = 1.4

coulombtype = Shift

;rcoulomb_switch = 0.0

;rcoulomb = 1.2

;epsilon_r = 15

vdw_type = Shift

rvdw_switch = 0.9

rvdw = 1.2

tc-grps = System

tcoupl = v-rescale

tau_t = 1.0

ref_t = 300

; Pressure coupling =

Pcoupl = Parrinello-Rahman

tau_p = 12.0 ; PREVIOUSLY 4

compressibility = 3e-4

ref_p = 1.0

;--------------------

; Free energy parameters

free-energy = yes

sc-power = 1

sc-alpha = 0.5sc-r-power = 6

sc-r-power = 6

; Which intermediate state do we start with?

-------

init-lambda-state = sedstate

; What are the values of lambda at the intermediate states?

;-------

vdw-lambdas = 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0

; This makes sure we print out the differences in Hamiltonians between all states, and not just the neighboring states

;--------

calc-lambda-neighbors = -1

; the frequency the free energy information is calculated. This

; frequency (every 0.2 ps) is pretty good for small molecule solvation.

;-------

nstdhdl = 10

; not required, but useful if you are doing any temperature reweighting. Without

; temperature reweighting, you don't need the total energy -- differences are enough

dhdl-print-energy = yes

couple-moltype = PEG

; we are changing both the vdw and the charge. In the initial state, both are on

couple-lambda0 = vdw

; in the final state, both are off.

couple-lambda1 = none

; we are keeping the intramolecular interactions ON in all the interactions from state 0 to state 8

couple-intramol = no

Please help, I have been stuck at this tutorial for far too long.