[gmx-users] Free Energy Calculations

[David Mobley]

Dear Justin,

> I'll preface this message by saying that this is my first foray into free energy
> calculations, and I am hoping that someone more well-versed in the matter may be
> able to help me sort out some strange results.  I have done a good bit of
> reading in terms of tutorials and in the archives, so I decided to give my
> system a try.  I am attempting to determine DeltaG of solvation for a
> polyphenolic compound.  My approach is to turn off charges in one step, and to
> use soft-core to turn off Lennard-Jones interactions in a second step, using
> thermodynamic integration to determine my DeltaG values.  Bonded parameters
> remain untouched.  I follow closely the procedure posted by David Mobley at the
> DillGroup Wiki site, changing only a few parameters for my simulation, as I am
> using Gromos96, not OPLS-AA.

Thanks. This is a good intro -- helps to know what level of help you need.

Which gromacs version are you using?

I vaguely remember, that at one point I may have found a bug relating
to atoms that initially started out with no charges and then had
charges that were turned on; I think it might have been a problem with
PME or something similar. This would have been a bug with 3.3.1 or
CVS, prior to the release of 3.3.2, I think; you can probably find it
on Bugzilla. I would think it would have been fixed for 3.3.2 if
that's what you're using, but it might be worth digging up the
bugzilla to see for sure. Let me know if you have any trouble finding
it.


> Relevant snapshots of my topologies and .mdp files follow this message.  The
> problem I am having (well, at least I think it's a problem), is that I get
> different results for the forward and reverse simulations.  I define "forward"
> as turning off the charges/LJ parameters, and "reverse" as turning said
> parameters back on.

OK. Can you explain where your starting configurations come from?
Ordinarily, one would expect that "forward" and "reverse" only are
meaningful when you're (a) doing nonequilibrium simulations, or (b)
beginning simulations at lambda values from the ending configurations
at preceding lambda values. You're not doing (a); are you doing (b)?
That is, does the simulation at lambda number i+1 begin with the
ending coordinates from the simulation at lambda number i?

Also, it might be useful to know how many (and which) lambda values
you are using for each transformation. Your equilibration protocol
might be useful to know as well.

I am not seeing anything wrong with your topologies, at least at the
level of detail you've sent.


I think if I were you, I'd begin by:
(a) Checking bugzilla
(b) Looking at the energy components of dV/dlambda at different lambda
values for your forward and reverse coulomb transformations. In
particular, the dV/dlambda components at lambda=0.2 in the forward
direction, say, should be equal (within uncertainty) to those at
lambda=0.8 in the reverse direction. Looking at the contributions
component-by-component (LJ, PME, etc) will help you nail down whether
there's a specific energy term which is causing most of the
difference. Figure out where the big differences you're seeing are
coming from. The differences are so big, even short test simulations
should be able to figure this out. This may ultimately point to either
(i) a bug, or (ii) some kind of problem with your topology that we've
missed.

Good luck!

David