[gmx-users] FEP trajectory errors
Robert Johnson
bobjohnson1981 at gmail.com
Tue Jan 22 16:29:38 CET 2008
Yep, all atoms begin to freeze....even the ones not being perturbed
(i.e. the water molecules). It's pretty strange. It only happens for
lambda=1.0. Below is my .mdp file.
; Run Control
integrator = md ; calculation type
dt = 0.0015 ; time step (ps)
emstep = 0.1 ; step size
for steepest descent (nm)
nsteps = 3000000 ; number of steps
gen_vel = no ; velocity
generation on startup
comm_mode = linear ; center of
mass restrictions
nstcomm = 100 ; center of
mass motion removal frequency
; Free Energy Perturbation
free_energy = yes
init_lambda = 1.0
delta_lambda = 0
sc_alpha = 0.5
sc_power = 1
sc_sigma = 0.3
; Data Output Options
nstxout = 5000 ; coordinate
output frequency
nstvout = 5000 ; velocity
output frequency
nstlog = 5000 ; energy
output frequency to log file
nstenergy = 5000 ; energy
output frequency to energy file
; Neighbor-list
nstlist = 15 ;
neighbor-list update frequency
ns_type = grid ; type of neighbor-list
pbc = full ; periodic
boundary conditions
; Cut-off distances
rlist = 1.2 ; cut-off
distance for neighbor-list (nm)
rvdw = 1.2 ; cut-off
distance for vdw (nm)
rcoulomb = 1.2 ; cut-off
distance for coulomb (nm)
dispcorr = EnerPres
; Electrostatics
coulombtype = PME ; electrostatics method
pme_order = 4 ;
interpolation order of PME
optimize_fft = yes ; optimization
of FFT grid at startup
fourierspacing = 0.12 ; grid spacing (nm)
On Jan 22, 2008 3:45 AM, Maik Goette <mgoette at mpi-bpc.mpg.de> wrote:
> Do really all atoms "freeze". So also the atoms, which are not perturbed?
> Can you post your mdp-file?
>
>
> 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:
> > The distortion was occuring because my position restraint potential
> > was too weak. Increasing the force constant keeps the base in the
> > right conformation.
> >
> > I believe my topology is correct, but it seems that turning off the
> > VdW parameters causes the base to distort. I will look into this more.
> >
> > Also, even with sc_alpha=0.5, I'm still experiencing this weird
> > "freezing" when lambda=1.0. The trajectory proceeds normally for about
> > 2.1 ns, then all of a sudden all motion freezes and the atoms simply
> > vibrate about their "equilibrium positions". Any other ideas what
> > might cause this?
> > Thanks,
> > Bob
> >
> > On Jan 21, 2008 3:39 AM, Maik Goette <mgoette at mpi-bpc.mpg.de> wrote:
> >> 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
> >>>
> >>> .
> >>>
> >> _______________________________________________
> >> 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