[gmx-users] non zero free energy values in mass transformation
Dan Gil
dan.gil9973 at gmail.com
Mon Aug 14 17:35:51 CEST 2017
Hi Deepak,
I'd like to share a reference with you. http://www.tandfonline.
com/doi/abs/10.1080/00268970600893060
If you could take a look at results (14) versus (15), and results (28)
versus (29)... You see that when only the potential energy is considered,
the mass doesn't contribute to the free energy. However this is based on
the assumption that the mass between the two states are equal. The mass
between Li+ and Cl- are rather different.
Please correct me if I am wrong, but I believe TPI compares between two
states, a solute in the ideal gas phase and the same solute in water (in
your case). E.g. [Li]_idealGas -> [Li]_aqueous. The mass is not changing
here.
On the other hand TI considers two states with different species. In your
case, it is [Li]_aqueous -> [Cl]_aqueous. The two states being compared is
different than TPI. In TI, the masses of the solutes are different. I think
the idea is that when you do a TI in water and do a TI in gas, the mass
contributions to the free energy cancel out. Can you also try doing
[Li]_idealGas -> [Cl]_idealGas? Hopefully, the free energy change due to
the mass perturbation for the two transformations cancel out to zero.
A disclaimer though, I am also working on a very similar problem without
success (yet). I am very curious to see if this works out for you.
Dan
On Sun, Aug 13, 2017 at 11:08 PM, deepak bapat <dubapat at gmail.com> wrote:
> Hi
>
> To study the efect of particle mass on excess chemical potential i did test
> particle insertion (tpi) in SPC/E water. I created a lennard jones
> particle with very small LJ parameters so that the probability of gettining
> it inside the cavity increases.
>
> In one case I kept mass of tpi particle equal to that of Lithium. In second
> simulation the mass was equal to that of Chlorine atome.
>
> LJ parameters for tpi particles in both cases were 2.0e-02 7.64793e-02.
>
> In both cases the excess chemical potential was in the range of 0.0798
> kJ/mol.
>
> tpi only considers energy as a function of coordinates only so i assume
> that it would depend only on LJ parameters as long as we are delaing with
> neutral systems.
>
> So looking at this I again come back to the mass perturbation which started
> all this.
>
> Since I was doing NPT calculations at 300 K I expected that average total
> kinetic energy of the system should not be different. The themostat would
> just adjust velocity distribution so that it brings the temperature of the
> system to the set temperature. Any increase or decrease in mass should be
> accounted for by the thermostat.
>
> As I said, in free energy calculations I got non zero value. (~ +6.2
> kJ/mol). So I am still worried about this non zero behavior.
>
> Thanks and regards
>
> Deepak
>
>
>
>
>
>
>
>
>
>
>
> On Sat, Aug 12, 2017 at 10:11 AM, deepak bapat <dubapat at gmail.com> wrote:
>
> > Thank you Andre and Dan.
> >
> > I'll read more about.
> >
> > Deepak
> >
> > On Fri, Aug 11, 2017 at 7:35 PM, Dan Gil <dan.gil9973 at gmail.com> wrote:
> >
> >> Hi Deepak,
> >>
> >> I am battling the same problem too. I am doing free energy calculation
> for
> >> heptane -> perfluoro-heptane in water. Also SPC/E water, Gromacs 5.1.
> >>
> >> Like Andre said, there is expected to be a free energy change when you
> >> change the mass of a particle. But if you do the same transformation
> from
> >> Li+ to Cl- in the gas phase, hopefully the free energy change will be
> >> opposite and (approximately) equal. Have you tried doing this second
> >> calculation [Li+ -> Cl-] in the gas phase? I think this website:
> >> http://www.alchemistry.org/wiki/Thermodynamic_Cycle explains why this
> >> second step might be needed better.
> >>
> >> Dan
> >>
> >> On Fri, Aug 11, 2017 at 9:37 AM, André Farias de Moura <moura at ufscar.br
> >
> >> wrote:
> >>
> >> > Deepak,
> >> >
> >> > I never did a mass perturbation myself, but I would expect some free
> >> energy
> >> > change, since the translational partition function depends on the mass
> >> of
> >> > the molecules.
> >> >
> >> > Also, having the same kinetic energy is not meaningful at all, the
> >> correct
> >> > quantity you should always use in statistical mechanics reasoning is
> the
> >> > momentum of the particle, p=mv, which clearly changes during your
> >> > alchemical transformation. Put in another way: the phase space of the
> >> > system is a 6N dimensional space of positions and momenta, and you
> have
> >> > just perturbed the momenta part.
> >> >
> >> > And by the way: that is a typical case in which Monte Carlo methods
> >> would
> >> > fail to reproduce an MD result, since only the configurational phase
> >> space
> >> > is sampled in MC (masses effects would have to be added afterwards,
> >> > assuming that some analytical result could be computed)
> >> >
> >> > I hope it helps
> >> >
> >> > Andre
> >> >
> >> > On Fri, Aug 11, 2017 at 6:36 AM, deepak bapat <dubapat at gmail.com>
> >> wrote:
> >> >
> >> > > 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
> >> > > --
> >> > > Gromacs Users mailing list
> >> > >
> >> > > * Please search the archive at http://www.gromacs.org/
> >> > > Support/Mailing_Lists/GMX-Users_List before posting!
> >> > >
> >> > > * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
> >> > >
> >> > > * For (un)subscribe requests visit
> >> > > https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users
> or
> >> > > send a mail to gmx-users-request at gromacs.org.
> >> > >
> >> >
> >> >
> >> >
> >> > --
> >> > _____________
> >> >
> >> > Prof. Dr. André Farias de Moura
> >> > Department of Chemistry
> >> > Federal University of São Carlos
> >> > São Carlos - Brazil
> >> > phone: +55-16-3351-8090
> >> > --
> >> > Gromacs Users mailing list
> >> >
> >> > * Please search the archive at http://www.gromacs.org/
> >> > Support/Mailing_Lists/GMX-Users_List before posting!
> >> >
> >> > * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
> >> >
> >> > * For (un)subscribe requests visit
> >> > https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
> >> > send a mail to gmx-users-request at gromacs.org.
> >> >
> >> --
> >> Gromacs Users mailing list
> >>
> >> * Please search the archive at http://www.gromacs.org/Support
> >> /Mailing_Lists/GMX-Users_List before posting!
> >>
> >> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
> >>
> >> * For (un)subscribe requests visit
> >> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
> >> send a mail to gmx-users-request at gromacs.org.
> >>
> >
> >
> >
> > --
> > Deepak U. Bapat
> >
> >
>
>
> --
> Deepak U. Bapat
> --
> Gromacs Users mailing list
>
> * Please search the archive at http://www.gromacs.org/
> Support/Mailing_Lists/GMX-Users_List before posting!
>
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
> * For (un)subscribe requests visit
> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
> send a mail to gmx-users-request at gromacs.org.
>
More information about the gromacs.org_gmx-users
mailing list