[gmx-users] non zero free energy values in mass transformation

[deepak bapat]

Dear Gmx Users

I am trying to find out free energy of transforming one neutral LI to one
CL in aqueous solution.

Water model is SPC/E. Gromacs version: 5.1.2

I am doing something like this

Insert 1 Li  olecule in water. (Li parameters were taken from Kasper P.
Jensen and William L. Jorgensen 2006.)

step 1: increase mass of atom from 6.94100 to 35.45300.

Step 2: change LJ parameters of LI to LJ of CL


Now my doubt is in step 1:  i get non zero free energy values.

*required nonbonded.itp file content is*

;LI+
opls_406   LI  3   6.94100     0.000       A    2.8700e-01 2.0929e-03
; Li+ parameters from jensen and jorgensen JCTC 2006 2 1499-1509

;LI+ grown
opls_407   LI  17   35.45300     0.000       A    2.8700e-01 2.0929e-03
; Li+ parameters from jensen and jorgensen JCTC 2006 2 1499-1509 ;; Mass of
chlorine


*LI.itp file content of alchemical transformation for mass is*
[ atoms ]
; id    at type     res nr  residu name at name  cg nr  charge   mass
1       opls_406    1       LI          LI       1      0.0        6.941
opls_407 0.0 35.45300


*.mdp file used is as follows*


integrator             = sd
tinit                      = 0
dt                         = 0.001
nsteps                  = 1000000
nstcomm              = 100
nstxout                 = 0
nstvout                 = 0
nstfout                  = 0
nstlog                   = 5000
nstenergy             = 5000
nstxtcout              = 500
cutoff-scheme            = verlet
nstlist                  = 10
ns_type                  = grid
pbc                      = xyz
rlist                    = 1.0
coulombtype              = PME
rcoulomb                 = 1.0
vdwtype                  = cutoff
rvdw                     = 1.0
DispCorr                  = EnerPres
fourierspacing           = 0.12
pme_order                = 6
tc_grps                  = system
tau_t                    = 1.0
ref_t                    = 300
Pcoupl                   = Parrinello-Rahman
tau_p                    = 1.0
compressibility          = 4.5e-05
ref_p                    = 1.0
free_energy              = yes
init_lambda  = 0.05 ;;; depending on state, lambda value to be changed from
0 to 1 ;;;;;
delta_lambda             = 0
nstdhdl                  = 10
gen_vel                  = no
constraints              = h-bonds
constraint-algorithm     = lincs
continuation             = yes
lincs-order              = 12


Free enegy of step 1: i.e. mass transfomation is around +6.2 kJ/mol.

In Gromacs 5.1.2 manual section 4.5 p.no 99 dependance of free energy on
lambda equation is given.

Now my doubt is if I am using the NPT ensemble then the total kinetic
energy of the system should be constant i.e. summation(p^2/2m) which is
being controlled by the thermostat.

But i am getting non zero value in this case? Is it normal to get something
like this or i have blundered somewhere?

Will be happy to share input files if it helps.

Thanks and Regards

-- 
Deepak U. Bapat