# Re: FEP question

From: Jerome Henin (jhenin_at_cmm.chem.upenn.edu)
Date: Fri Sep 07 2007 - 10:33:47 CDT

Hi Luca,
Fixed atoms do constitute a bias. I'm afraid the argument about
infinite masses, though formally interesting, does not change that
state of affairs.
One can adopt either point of view:
1) Constraints are applied by adding constraint forces (deriving from
a constrained Hamiltonian). Hence a bias in configurational sampling.
2) The Hamiltonian is unconstrained, but constraints are applied by
setting infinite masses on a set of atoms. As a consequence, that part
of the system has an infinite relaxation time. Obviously, relaxation
does not occur during the simulation, which is therefore not sampling
the equilibrium ensemble, but another ensemble dictated by the initial
conditions, aka the constrained ensemble.

So the "no force, no bias" idea only goes so far...
Jerome

On 9/7/07, luca <bellucci14_at_unisi.it> wrote:
> Hi all,
> I think you are right....but I have one observation.
> Can you test use of fixed atom only in initial step?
> In this way you do not have any additional force therefore no bias.
> If you set fixed atom you have put infinite mass to some atoms but you do not
> change potential form of hamiltonian value.
> In this way you Infuence only the sampling at initial step.
> However you should minimize the number of restraint for a good sampling.
>
> Best,
> Luca
>
>
> > Matthew,
> > I think your ideas are essentially the right things to do.
> > In principle, you can use a custom Tcl script that enforces soft
> > restraints to, as you said, keep the two rings roughly in the same
> > place.
> >
> > When dealing with isolated molecules, this is not a problem at all,
> > since at the end-points (lambda=0 or 1), one of the molecules is
> > virtually non-existent, so tethering it to the other does not bias the
> > energies. Only points at 0 < lambda < 1 are affected, which just means
> > that the free energy is computed along a different pathway - a
> > better-behaved pathway, actually.
> > Now, the issue with your DNA residues is that even at the end-points,
> > the ghost residue still has bonded terms, so tethering it may induce a
> > bias.
> >
> > Still, I see two possible courses of action:
> > 1) add such restraints coupling the two residues, but in a way that
> > does not add a significant conformational bias - this may be possible
> > if "native" conformational fluctuations are small.
> > 2) as you also suggest, add torsional restraints (same kind of custom
> > Tcl script) to the ghost residue until its nonbonded interactions are
> > strong enough. The issue there is that in theory, it would require an
> > FEP step to determine the FE difference between the restrained residue
> > and the unrestrained one... this may be done by growing the restraints
> > using the "conformational free energy" module of NAMD, but it is still
> > adding a nontrivial step to a nontrivial calculation.
> >
> > You see, there is no cheap way out, but I hope it helps nonetheless.
> > Best,
> > Jerome
> >
> >
> > On 9/6/07, Matthew WIlce <matthew.wilce_at_med.monash.edu.au> wrote:
> > > Dear NAMDers,
> > > I have a couple of questions regarding FEP.
> > >
> > > I am using FEP and mutating a Cytosine into a Thymine, ADE and GUA. I
> > > have set up the dual topologies so that it is the purine or
> > > pyrimidine that is mutating while the rest of residue in not
> > > changing. I have noticed that at the beginning of the FEP run the
> > > purine or pyrimidine that is starting to grow moves around too much
> > > and sometimes the ring flips ~180 degrees. When the ring flips 180
> > > degrees it prevents a sensible outcome from the FEP. This also occurs
> > > for the Cytosine ring towards the end of the FEP. Once lambda has
> > > become greater than about 0.01 and the growing ring starts to
> > > interact with the rest of the molecules it behaves much better.
> > >
> > > For a purine to purine mutation eg ADE to GUA I could of course
> > > modify the dual topology so that the all of the common atoms are
> > > maintained throughout the FEP, and the same goes for the pyrimidines.
> > > But I cannot see how a work around when the rings are quite different
> > > eg. pyrimidine vs purine.
> > >
> > > Is there a way to sensibly restrain/constrain the movement of the
> > > growing component of the FEP so that it cannot 'float' around until
> > > lambda1 passes say 0.01, and also apply similar restraints to the
> > > shrinking component when lambda 1 is greater then ~0.99.
> > >
> > > Has anyone done this? What impact could this have on the free energies?
> > >
> > > I have thought of writing a tcl script that constrains the atoms at
> > > the beginning and end of the FEP depending on lambda ?? but I am not
> > > sure how. I want the growing and shrinking rings of the purine or
> > > pyrimidine to stay roughly in the same place relative to each other
> > > until their interactions with the rest of the molecule may cause them
> > > to move if necessary and not have them moving because they are not
> > > interacting.
> > >
> > > I have tried using the namd parameters as suggested in the most
> > > recent FEP tutorial as well as varying them but I have not been able
> > > to resolve this. Even using a very low number of steps in each FEP
> > > window (as low as 1000) does not stop the ring from sometimes flipping.
> > >
> > > Thanks in advance,
> > > Matthew
> > >
> > >
> >
> >
>
>
>

This archive was generated by hypermail 2.1.6 : Wed Feb 29 2012 - 15:46:56 CST