[gmx-users] FEP trajectory errors

[David Mobley]

Robert,

> I'm computing the free energy of binding of a DNA base on a carbon
> nanotube. I think it's a pretty simple calculation and I'm proceeding
> in a very standard way. This is what I'm doing:

An absolute free energy? This isn't necessarily straightforward --
there are a lot of wrinkles. Some of my recent work has some
discussions if this if it's helpful.

> I have the optimal orientation of the base on the nanotube. I'm
> constraining the positions of the base atoms with a soft harmonic
> potential.

OK, that's a good idea for standard state reasons (see the 2003
Boresch paper) among other things.

> I then am running two different FEP calculations: one where I turn off
> the charges on the base atoms and a second where I then turn off all
> the lennard-jones parameters for the base atoms. For each of these I
> use the following:
> delta_lambda = 0
> sc_alpha = 0.7
> sc_power = 1
> sc_sigma = 0.3
>
> I run a series of trajectories at constant lambda values from 0 to 1.

OK. I have some information posted on what I think the best settings
are for these on alchemistry.org, especially at
http://www.alchemistry.org/wiki/index.php/Best_Practices. It is still
in progress but may be somewhat helpful. I doubt this is the source of
your problems though, although I prefer a sc_alpha = 0.5 which may
improve the situation somewhat.

> However, I notice some problems with the trajectories when I turn off
> the LJ parameters. As lambda is varied from 0 to 1, it seems that the
> position restraints no longer are being applied. Additionally, the DNA
> base geometry starts to become severely distorted at lambda values
> greater than about 0.6. This happens despite the fact that I am not
> perturbing the internal bonded interactions of the base. Here is a
> sample of my topology (included just the first line of each section):

Hm, on the position restraints, depending on your version of gromacs
you may need to additionally specifiy the B state position restraints,
else these may be assumed to be zero. I believe you can do this by
adding additional columns in your posre.itp that give the B state
parameters.

In terms of distortion, it's not clear to me why you think the
molecule should *not* be distorted once you turn off the LJ
interactions (which includes internal LJ interactions). The bonded
parameters (especially torsions) are of course derived after LJ
parameters are already fitted, in most cases, so they only can be
relied on to give meaningful conformations when the LJ interactions
are on. The beauty of the thermodynamic cycle approach, though, is
that you can (at least in principle) get correct binding free energies
even if your molecule of interest samples wacky, artificial
conformations when it's noninteracting.

> Also, I am experiencing additional problems when lambda=1. After about
> 2ns, all the motion in the system begins to freeze and all the atoms
> simply vibrate about a fixed position. Soon after that the simulation
> crashes.

My experience has been that sc-alpha = 0.7 can lead to a fairly steep
dV/dlambda and some large forces near lambda=1, which is why I
recommend sc-alpha = 0.5. This parameter is actually fairly sensitive.
This *could* lead to crashes, but I don't know about
freezing/vibrating.

Let me know if this helps.

Best wishes,
David Mobley
UCSF
http://www.dillgroup.ucsf.edu/~dmobley

> Can anyone comment on what's going on here?
> Thanks,
> Bob Johnson
> _______________________________________________
> 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
>