[gmx-users] non zero free energy values in mass transformation
deepak bapat
dubapat at gmail.com
Wed Aug 16 10:51:20 CEST 2017
Hi Dan
Thanks a lot. I haven't tried a simulation in vacuum yet. I'll give it a
try.
Deepak
On Mon, Aug 14, 2017 at 9:05 PM, Dan Gil <dan.gil9973 at gmail.com> wrote:
> 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.
> >
> --
> 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
More information about the gromacs.org_gmx-users
mailing list