[gmx-users] Free Energy Calculations

[David Mobley]

Justin,

Following up on this, I looked back through my old e-mails and found
that, indeed, you do appear to be running into a known bug. It's one I
helped a guy named Bharat find, so the relevant bugzilla is from him
rather than from me: http://bugzilla.gromacs.org/show_bug.cgi?id=175.

It appears you will need to use the CVS version in order to get around
this. Or just give up (for the time being) the idea of doing the
calculations in reverse also.

TO THE DEVELOPERS: I've said this before, and I'll say it again: There
REALLY should be a list of KNOWN BUGS in each particular version of
gromacs somewhere obvious (aside from buried in Bugzilla). That is,
someone running a certain calculation in a particular version of
gromacs should be able to check some specific page and find out if any
of the things they're doing will be affected by any known bugs. Better
yet, one could input a topology file and mdp file to a web page which
would run a script to check it to see if it would use any routines
which are known to be problematic. Right now, I still don't know how
folks like Justin would be expected to know that they're doing
something that is known to be affected by a bug. It seems like he
would have needed to (a) read the entire history of the mailing list,
and (b) wade through all of the old bugzillas.


Thanks,
David


On Feb 16, 2008 3:19 PM, David Mobley <dmobley at gmail.com> wrote:
> 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
>