Date: Thu, 25 Apr 2002 17:13:23 +0200
From: David Smith
Subject: Re: PMF calculations

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.