Re: Free-energy perturbation calcs of zinc<->peptide binding
I have parameterised the zinc-ion in AMBER 3.0A as an 'H-bonding' type
of atom ie. it has 10-12 well-depth and equilibrium bond distances with
all relevent donor atoms on the peptide (peptide C=O, carboxylate-O,
amine H2N's etc).
Why do you see zinc as hbonding? If it were negatively charged there might
be some basis. 6-12 would be the way to go, I believe.
These calculations are sensitive to the nonbonded potential you choose,
so I suggest serious thought on how you determine 6-12 parameters (or
10-12 if there's a reason for it). See e.g. J. Aqvist, J. Phys. Chem.
1990, 94, 8021-8024, and remember that ion-water parameters might be
different from ion-peptide ones when water is modelled as a big sphere
with a few charges in it.
In addition, it has a point-charge calculated by fitting
the electrostatic potential with MOPAC.
What charge results? Was this for the ion alone? Why not use 2+ ?
Am I right in thinking that I need to do 3 molecular dynamics
simulations:
(i) using MD to equilibrate the previously minimised structure to 298K -
this may take 10ps or so with 1-2fs steps. Save the 'equilibrated'
velocities and coords.
You might want to warm carefully and measure the geometry around the
ion especially to get an idea of the time period of the motions and
as a guide to when equilibration is satisfactory.
(ii) start the GIBBS slow-growth using the equilibrated velocities and
coords with IELPER=1 and vary LAMBDA from 1->0 over say, 20ps.
Then run backward to check your hysteresis before continuing. You
could even run forward to say .8 & then back to see how long your
run should be (longer run => less hysteresis). Since the electrostatic
effect is strong and long range, you may need 50-200 ps each way. You
may want to break it up into several runs so that finer granularity
is used at higher charge density.
A safety measure might be to equilibrate some more at 0 charge before
continuing.. (requires a special parm topology file for the 0 charge
or running a window perturbation with only a long equilibration 1st window).
(iii) start GIBBS again using the equilibrated velocities and coords
with IELPER=-1 and again vary LAMBDA from 1->0 over the same interval.
Then run the gibbs(iii) back from small, non-0 lambda to check hysteresis.
(Running backward from 0 lambda gives a ghost that can 'appear' on top of
another atom.) My guess is that this perturbation can be much shorter
than the charge one.
How do I decide if 20ps was long enough?
Do it again for 50ps from scratch?
This practice is catching on; if you have the time.. My experience with
changing the size of monovalent cations in DNA is that owing to the
simplicity of the perturbation, getting, say, <5% hysteresis is sufficient.
But for big charge perturbations, this may not apply.
Another note: don't use your Gibbs parm topology file with the md program.
Another paper: Straatsma & Berendsen, J. Chem. Phys. 89:9, 5876-5886,
1988 (ion charge perturbations in water).
Bill Ross