[gmx-users] FEP trajectory errors

[Maik Goette]

Robert

A value for alpha of 0.5 should work, although I tested some different 
alphas with different test systems and I have to say, that a proper 
value can be between 0.2 and 0.7 depending on the system itself.
The distortion, you describe is really strange. I calculated DNA and 
single bases a lot and never saw something like that. Even a fully 
"dummy-base" is stable, if your topology is correct. The bonded terms 
are sufficient to keep it in the right conformation.
Are you sure, your position restrains in the topology are used with the 
correct structure file?

Maybe, I should also mention, that, especially with such large 
perturbations as a full DNA-base, you should split your topology into 2, 
where you switch off the QQ hardcore and the VdW softcore in two steps.

Regards

Maik Goette, Dipl. Biol.
Max Planck Institute for Biophysical Chemistry
Theoretical & computational biophysics department
Am Fassberg 11
37077 Goettingen
Germany
Tel.  : ++49 551 201 2310
Fax   : ++49 551 201 2302
Email : mgoette[at]mpi-bpc.mpg.de
         mgoette2[at]gwdg.de
WWW   : http://www.mpibpc.gwdg.de/groups/grubmueller/


Robert Johnson wrote:
> Thank you David and Maik for your detailed replies. Yes, I am trying
> to obtain an absolute free energy of binding. My thermodynamic cycle
> is:
> 
> NT+DNA(adsorbed) -> NT+ DNA(desorbed)  (obviously this is not the one
> that I'm obtaining with the alchemic method)
> 
> Thus, the one I'm using in the simulation is:
> 
> Water + NT + DNA(adsorbed) -> Water + NT + nothing (DNA dissapears)
> Water + DNA -> Water + nothing
> 
> Any suggestions about this? I will check out your recent publications.
> 
> I believe my topology is properly formatted. I wasn't expecting the
> large distortions because I was using position restraints. Thus, I was
> expecting the base to remain in about the same geometry. Perhaps I
> will increase the force constant for the restraining potential.
> 
> I am going to repeat some of these computations using sc_alpha=0.5.
> 
> Thanks for your help.
> Bob
> 
> On Jan 17, 2008 2:24 PM, David Mobley <dmobley at gmail.com> wrote:
>> 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
>>>
>> _______________________________________________
>> 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
> 
> .
>