Date: Fri, 19 Mar 2004 15:38:40 -0600
From: Chris Moth
Subject: AMBER: About RESP Calculation

>At 01:20 AM 3/19/2004, Ilyas Yildirim wrote:
>Hi Amber Members,
>
>I have a couple of question about RESP calculation. I found the .pdb
>file of a Cytosine molecule in Cambridge Crystalographic Database,
>CYTIDI10. I want to find the RESP charges of this particular molecule and
>have a couple of questions to ask:

I have been performing RESP parameterizations of small molecule enzyme
inhibitors. As a relative newcomer (lowly grad student) to RESP
parameterization myself, I thought it might help if I mention a few of my
personal experiences, observations, and conclusions.

I have found there to be a gulf between the (excellent) treatment of RESP
in the literature, and specific steps needed to parameterize a molecule
prior to starting an MD simulation. My comments reflect _my_ navigation of
this gulf, and share some lessons I have learned the hard way.

If any of these ideas are misguided, I encourage correction from the more
experienced folks on the list!

>1. I have tried to use the program RED written by Mr. Dupradeau. I entered
>the .pdb file and it took 3 days and still didnt end up with the
>calculation. I killed the job. Are the gaussian calculations taking
>that long?

I have seen g98 geometry optimization calculations take this long when the
start structures are far from a reasonable local minimum in energy.

BTW, here is the third line of my g98 geometry optimization input file:

#RHF/3-21G* OPT Test

Then, for MEP calculation prior to RESP we use:

#RHF/6-31G* SP Pop=ESP iop(6/33=2) TEST

The other folks on the list who have already responded know far more about
this however...


By first coarsely minimizing the molecular structure using a
quick-and-dirty electrostatic assignments in a package like MOE first, I've
found that g98 runs to completion within 12 hours or so. (50 to 60 atom
molecules - 300MHZ R12000 SGI Uniprocessor machine) You can probably get
the same kind of results with INSIGHTII. If there is a MMF94 force field
option, you might give that a try.

> I did not change the ATOM NAMES, and later I realized that in
>the sugar molecule, there is a CH2 term. Is it possible that this is the
>reason for this problem? (If it is a problem for the gaussian to calculate
>a particular job more than 3 days)

Gaussian does not care about atom types or names or bonds - just atomic
numbers (element symbols) and their coordinates.


>2. Do we have to optimize the structure of a particular molecule in order
>to calculate the RESP charges? For instance, can I just use the crystal
>structure I got from the CCD and just do the MEP calculation and then the
>RESP calculations? Is this legitimate?

Here is exactly the rub. My read of the literature is that the ideal
electrostatic parameterization would involve Bolzmann weighted averaging of
charges derived from RESP fits of all reasonable conformers. However, that
is not computationally tractable.

So - there is a real human factor here. Things are going to depend on your
particular project, and your chemical intuition.

Charges do vary considerably with conformation, and slightly with
orientation, to be sure. (I've played around with this myself - and have
seen variations of up to 3/10ths of a q through conformational change) I
get through the day by taking the Cornell et. al. charges as they come -
since they have been extensively put to the test, and have been shown to
exquisitely mimic natural phenomena (see papers predicting A vs B form DNA
duplex transitions, oligo/ligand binding, and on and on etc).

For parameterization of my own molecules, I have had to factor in the
downstream analysis needs of my work. (I.e. - I goofed things up the first
time around, and now am doing a better job). For example - a couple of
rules of thumb I work with now are:

- Ligand coformers with intramolecular hydrogen bonds should not considered
(this was mentioned in the original Cornell et al paper). RESP fitting of
these conformers places far too much electronegativity on the hydrogen
bonding heavy atoms. And, this conformation is not going to be "right" for
ligand bound to protein

- Stereoisomers should be parameterized identically. This seems obvious on
the surface... and, indeed, in my experience, RESP fitting mirror image
molecules did result in differences in only the 1/100ths of a q. By
bypassing RESP fitting of the sceond enantiomer, and manually setting
identical charges, you honor the MM model. And, at the end of the day, this
is the only way you can make clean comparisons of binding energetics
between R vs S enantiomers at the MM level of theory, anyway.

So, I have concluded that is more important for electrostatics to be
consistently derived throughout a project than be "correct". But yes, at
some level, it must be the case that including more conformers in a RESP
fit will improve "correctness" So, for your work, perhaps averaging a
putative global minimum (use coarse MOE systematic conformational search)
conformation, and a crystal structure conformation might be a way to go.

And of course, if you see a positive charge on an oxygen, you've goofed
somewhere :)

>3. In the RED program, it is talking about re-orientation. I followed the
>manual, and did not understand exactly what it is meant by reorientation.
>In the tutorials there are examples for re-orientation, but after getting
>an optimized structure, why do we have to reorient the structure? The
>structure has a clear structure when viewed in InsightII or Molden. So,
>isnt it legitimate to do the MEP calculation in Gaussian and then the RESP
>calculations without reorient the structure?

I think others have addressed this far better than I can. I have not yet
used RED - but Gaussian does build a rectangular box of electrostatic
points around the ligand, and the orientation of the box vs the ligand will
affect final charges assignments a bit.

Perhaps in a perfect world, g98 would have an option to build a perfect
unvarying sphere of points around the ligand center-of-mass :) But, I'm
sure that would defeat other requirements of the MEP which involve building
ESP points at fixed multiples of VdW radii......

I hope some of this has provided useful "food for thought" without being
too informal or too much in error.

Best of luck with your AMBER simulation work!


Sincerely

Chris
For trajectory
visualization: http://www.structbio.vanderbilt.edu/~cmoth/mddisplay