Date: Thu, 25 Apr 2002 13:36:41 -0500
From: David Pearlman
Subject: Re: PMF calculations


NCORC is for determining the contributions of any constrained
lambda-dependent INTERNAL coordinate (bonds, angles, torsions) to the free energy.
If you turn off shake, you don't remove the contributions of the internals, you
just make it harder to calculate them--you still need to include the contributions
of the changing internal coordinates.

In general, for ANY Gibbs free energy simulation you should turn on SHAKE for
all bonds, set NCORC=1 and set INTPRT = 5. These should be the defaults for
a free energy simulation.

dap

David Smith wrote:
>
> Sophia Kondratova wrote:
> >
> > Hello,
> >
> > I have a question concerning "PMF bond contribution method", which is used to
> > calculate the intra-perturbed group interaction contributions to the free
> > energy while doing FEP calculations (NCORC and INTPRT flag).
> >
> > What other parameters besides NCORC=1 and INTPRT=1,2,3,4 or 5 have to changed
> > in order to calculate these intra-pert group contributions? I am trying to do
> > an FEP of guanine to adenine and when I set NCORC=1, I get values of ~600
> > kcal/mol, which is completely wrong (I do get the right answer with NCORC=0
> > and INTPRT=3)
> >
> > Any help would be appreciated
> >
>
> I have another answer and related question.
>
> ANSWER:
>
> As far as I undertand NCORC is really for calculating the contributions
> of (lambda dependent) constraints to the free energy. Thus you really
> need to have some constraints for it to make sense. You can impose these
> individually with INTR but the more general method is to use holonomic
> constraints like SHAKE. Thus if you want calculate the bond length
> contribution for your G->A mutation you could set NTC=3 and NTF=3. If I
> understand correctly this should constrain all bond lengths for which
> the equilibrium distance is a function of lambda, and then calculate the
> energy required to impose these constraints.
>
> QUESTION:
>
> Is this answer correct ? My above answer comes from what I think I
> understand from the literature. However, the manual says that NTC=3 in
> Gibbs constrains all bonds. Is this really all bonds in the system, or
> all bonds in the perturbed group, or all bonds in the perturbed group
> with a lambda dependent equilibrium distance.
>
> The second part of my question relates more directly to my own work.
>
> I have a C-C single bond for which I developed my own parameters.
> Because of the 1-4 electrostatics I needed to set the equilibrium
> bondlength for this bond (in the lambda=1 state) to something like 1.486
> with a force constant of 317 to reproduce the X-Ray distance of 1.532.
> For the lambda=0 state, the equilibrium bond for the corresponding C-C
> bond is 1.469 (K=317) to reproduce the X-Ray distance of 1.563. As a
> relative newbie in the field, I developed these parameters before I had
> heard of the overlooked bond-stretching contribution in FEP
> calculations. In the gas phase and water simulations they end up with an
> average value similar to the X-Ray.
>
> Now I would like to run the perturbation. It seems to me that using
> SHAKE to constrain these bonds to their equilibrium distances would be
> disasterous.
>
> Should I:
>
> (1) redo the parameters by playing with the force constant instead of
> r(e) ?
> (2) overlook the overlooked bond stretching constribution ? or
> (3) define specific restraints with INTR to keep the bonds at their
> "natural" lengths ?
>
> Thanks once again for your all of your interest, patience, and
> wonderfully enlightening replies...
>
> Dave.
>