[gmx-users] Large dV/dlambda value at lambda = 0

[David Mobley]

Bob,
> I'm computing the relative binding free energies of DNA bases on a
> nanotube using the protocol listed here:
> http://www.dillgroup.ucsf.edu/group/wiki/index.php/Free_Energy:_Tutorial
> I'm using the following steps:
> 1. Turn off all charges on base
> 2. Transform LJ parameters
> 3. Turn on all charges
>
> Because the transformation (from G to A for example) involves changing
> every single atom charge, I decided to first reduce all charges to
> zero. The problem I'm experiencing occurs at small lambda values
> during step 2 (transforming LJ parameters) of the above routine. I
> obtain large dV/dlambda values of about 2E3 for lambda 0.05 or less.
> The rest of the curve contains dV/dlambda of about 6E2. I don't think
> this is due to any equilibration issues. Each lambda point is run for
> 15 ns and I incur a small error of ~1E-1 for each dV/dl value. The
> error is a little bit larger (~1) for the lambda=0 point. I believe
> the topology and simulation protocol is correct. I have already run
> calculations for the absolute binding free energy of this system and I
> get good agreement for steps 1 and 3. However, it seems that step 2 is
> problematic. If I can't get rid of the large dV/dlambda value, how
> should I handle that? Is this large value an artifact? Or should I
> assume that the "real" dV/dl curve can be extrapolated to a value that
> is more reasonable and in better agreement with the rest of the curve?
> If anyone has some suggestions on how to alleviate this problem, I'd
> appreciate reading them.

This doesn't sound that good, though I don't see any reason in
principle why dV/dlambda can't be particularly large for some region
of the curve. I assume you're using soft core potentials with the
parameters I recommend? Also, you doublechecked your topology and you
haven't forgotten to turn off some of the charges?

If this doesn't go away with longer simulations, it is probably real
and hence you shouldn't attempt to "fit" it out.

What exactly do you mean by error? Integration error?

Generally with TI you want to avoid any assumptions about smoothness,
especially if your lambda spacing is not that close. If one point is
substantially different from the others it may mean that there is a
feature in dV/dlambda there that you are mostly missing (and only
seeing at one lambda value) and hence that you may need to space your
lambda values more closely. You could run some additional simulations
nearby to see if it is real. Do NOT assume the curve is smooth and the
feature is an artifact unless you have evidence to the contrary. (I
can send you a couple references with bad things that could happen if
you do this, if you don't see why it's a problem).

This is, of course, one of the reasons TI is less than ideal -- there
is no way error analysis in TI can tell you that your lambda points
aren't spaced closely enough, and if there are significant features
these may lead to large integration error, or worse you may be
*missing* significant features. That's why we typically use BAR, but
this requires evaluating energiesf rom each trajectory in other lambda
values. This can be done using mdrun -rerun, as I've mentioned before,
but requires storing lots of snapshots.

David

>
> Thanks,
> Bob
> _______________________________________________
> gmx-users mailing list    gmx-users at gromacs.org
> http://www.gromacs.org/mailman/listinfo/gmx-users
> Please search the archive at http://www.gromacs.org/search before posting!
> Please don't post (un)subscribe requests to the list. Use the
> www interface or send it to gmx-users-request at gromacs.org.
> Can't post? Read http://www.gromacs.org/mailing_lists/users.php
>