[gmx-users] Free energy of discharging and then recharging not zero

David Mobley dmobley at gmail.com
Thu May 8 23:45:27 CEST 2008 

Dear Bob,

> Yes...sorry...I'm doing a 7.5 ns for each lambda value (lambda = 0,
> 0.25, 0.5, 0.75, 1.0). Every single simulation takes place from the
> exact some starting configuration. No restraints are being used. The
> calculation I'm planning on doing (once everything is worked out) is a
> relative binding free energy calculation where I morph G to A. But for
> now I'm just taking G and turning the charges off and then turning
> them back on. I do not perturb the LJ parameters at all so I'm not
> using soft-core potentials or any of that stuff.
>
> The trajectories are all qualitatively the same - irrespective of
> lambda value. The base adsorbs very strongly to the nanotube. For each
> trajectory, the base simply diffuses around on the surface of the
> nanotube. To the base, the nanotube appears almost as a uniform
> cylinder - the base feels very little corrugation. Plus, there is no
> electrostatic interaction between the base and the nanotube. It is
> completely a van der Waals effect. Thus, the base doesn't have the
> chance to desorb because I'm not perturbing the LJ parameters.

Hmm, that's interesting. What's your equilibration protocol? I assume
you're allowing sufficiently long for equilibration? If not, you could
see what appear to be artificial fluctuations due to that... I can
imagine equilibration could be rather slow here as these are charged,
right?

Also -- you are changing the net charge of the system, now that I
think about it. There are some complications associated with this that
I haven't entirely worked out for myself yet, but I believe there may
be a correction you need to apply to the computed free energy that has
to due with long range effects and/or boundary conditions. I believe
Gerhard Hummer worked out this correction in the 1990s so you may want
to look there. I believe the correction is system size dependent, and
it's possible the correction might go in different directions
depending on whether you're turning on or turning off charges. It
shouldn't cause fluctuations, I wouldn't think, but it's probably
another thing worth thinking about.

> It seems that in all the trials I have done, I always obtain a larger
> error that doesn't converge whenever I turn ON charges. I don't really
> see anything drastically different occuring in any of the
> trajectories, so it's hard to believe that it could be a sampling
> problem.

On the other hand, it's hard for me to believe it could *not* be a
sampling problem. Directionality *can't* matter here if the
simulations are sufficiently long. That is, if you're linearly scaling
lambda (sc-alpha=0), then, <dV/dlambda> at 1-lambda going in one
direction has to be exactly equivalent to <dV/dlambda> at lambda in
the other direction. Fluctuations should be the same too, right?

You could actually do a numerical test of this, I suppose. In
particular, suppose you save a full precision trajectory consisting of
a few hundred frames for your system. For the purposes of the test,
the particular lambda value you use to generate the trajectory
shouldn't matter. Then, take that stored trajectory use mdrun -rerun
to reevaluate the energies and dv/dlambda for your trajectory in two
different cases. Case 1: Use a topology for turning off the charges,
and run at a lambda value you pick (say, lambda = 0.3). Case 2: Use a
topology for turning on the charges, and run at the matching lambda
value (if lambda =0.3 for case 1, use lambda = 1-0.3 = 0.7 here).
Check that (a) the <dV/dlambda> you get from the two cases is equal,
and (b) the fluctuations in <dV/dlambda> that you get are equal.

If, indeed, they *are* equal, then you've confirmed that there's no
problem with the numerics, and that forward transformations are
equivalent to reverse transformations (as it they be). I suspect
that's what you'll find.

If, indeed, the numerics are OK here, then I think one of the
following things must be going on:
(a) You are not equilibrated
(b) You have some sort of sampling problem
(c) Your lambda values going in one direction are not equivalent to
the (1-lambda) values going in the other direction so you have
different integration errors.
(d) You need a correction because of changing the net charge of the system

I really don't see any reason (still) that you'd see bigger
fluctuations going in one direction than in the other, except if you
aren't equilibrated or you have a sampling problem.

Let me know what you figure out.

For what it's worth, reading ahead, I'm not happy with Maik's
solution. First of all, you *should* be able to do this (whether it's
the ideal way to do it or not is another matter). Second, the behavior
you're describing strikes me as really odd and I think you should get
to the bottom of it to make sure there isn't a bug.

Hope this helps,
David

> Bob
>
>
> On Wed, May 7, 2008 at 4:12 PM, David Mobley <dmobley at gmail.com> wrote:
>> Bob,
>>
>>
>>  > Hello everyone,
>>  >  I'm trying to calculate the free energy of binding of DNA bases on a
>>  >  carbon nanotube. I'm running some tests to make sure that I'm doing
>>  >  everything correctly. One thing I tried was turning off all the atom
>>  >  charges of the DNA base and then turning them back on again.
>>  >  Theoretically, the free energy changes of these two processes should
>>  >  be equal and opposite and thus sum to zero. However, this is not what
>>  >  I'm finding.
>>  >
>>  >  For guanine, I get a free energy change of 648 kJ/mol and -618 kJ/mol
>>  >  for turning off and turning on the charges, respectively. Obviously,
>>  >  they are not equal by 30 kJ/mol, which seems pretty big. I have done
>>  >  some error estimation using the g_analyze -ee program. One thing I
>>  >  find strange is that the error estimates in dV/dl for TURNING ON the
>>  >  charges is large (over 2) and do not even converge for a 7.5 ns
>>  >  simulation. In contrast, the error in dV/dl for TURNING OFF the
>>  >  charges converges extremely quickly (using small block sizes of 50 or
>>  >  less) and is smaller at 0.3. So it seems like I have some sampling
>>  >  problems with the TURNING ON portion. Is there some reason why you
>>  >  must sample a longer trajectory when turning on the charges?
>>
>>  Er, can you provide a few more details here? i.e., when you say "a 7.5
>>  ns simulation", do you mean 7.5 ns at each lambda value? (I appreciate
>>  that you're doing the error analysis separately at each lambda,
>>  though).
>>
>>  Are you restraining the base to be near the nanotube in a particular
>>  orientation? For example, I could imagine that if the base had drifted
>>  away from the nanotube after turning off the charges, then it might be
>>  hard to converge the calculation for turning on the charges (and the
>>  calculation for turning off the charges would probably not be
>>  converged either, though it might appear converged).
>>
>>  The only reason I can see why you might have different convergence
>>  properties for turning on the charges than turning off the charges is
>>  if the starting state in the two sets of simulations is different (or,
>>  if you are starting the simulation at each lambda value from the end
>>  point from the preceding lambda value). Are you doing one of those
>>  things?
>>
>>  David
>>
>>
>>  >  I'm following the procedures of
>>  >  http://www.dillgroup.ucsf.edu/group/wiki/index.php/Free_Energy:_Tutorial
>>  >
>>  >  Does anyone know the reason for the discrepancy between these two
>>  >  (seemingly identical) processes?
>>  >
>>  >  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
>>  >
>>  _______________________________________________
>>  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
>>
> _______________________________________________
> 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
>

More information about the gromacs.org_gmx-users mailing list