AMBER Archive (2008)

Subject: RE: AMBER: RESP GAUSSIAN AMBER

From: FyD (fyd_at_q4md-forcefieldtools.org)
Date: Fri Dec 12 2008 - 08:25:06 CST

Hi Ross,

> Which reference is this specified in? I could never find a definitive

> document that specified exactly what options were used for Gaussian fits.

The link is @
http://www.teokem.lu.se/~ulf/Methods/resp.html#How_to_determine_RESP

This means the tutorial @
http://ambermd.org/tutorials/advanced/tutorial1/section1.htm
demonstrates the way used by Pr Ryde's lab to derive RESP charge, &
not that used in the building of the AMBER force field topology
database. In your tutorial (dedicated to charge derivation), it is
important to explain this choice, but also the reason of this choice.
Thus, new users will understand which options to use, & will have the
opportunity to tests the difference cases. I think this would be
really less confusing...

There is obviously not a unique answer: It is why the different
options should be clearly established.

> I think these options were found by the students that wrote this tutorial

> from the attached webpage. I have no idea what the history of this page is.

"these options were found _by_ the students" & "I have no idea".

  Waouh ! Here, I prefer not answering...

> However, do you know if these options actually make a measurable difference?

Sure, there is a difference. Why do you think those options are used ?
& indeed, IOp(6/33=2,6/41=10,6/42=17) might be particularly
interesting in some cases.

> Sure they might give slightly different charges but then my experience with

> RESP is that it is only marginally reproducible at best anyway.

I am not sure I agree here (see below).

> It seems to

> me that, for example, upping the number of points per layer and the number

> of layers should simply result in the charges generated converging to a

> specific value as a function of the number of points used. If one gets

> different charges at the default settings then the default options are

> surely not large enough to obtain convergence in the results and will just

> lead to large error bars.

May be you should do some tests here, & explain this point in your
tutorial as the use of a fine grid for organic molecules is reported.

> I guess people should read the disclaimer at the top of all the tutorials:

>

> "These tutorials are meant to provide illustrative examples of how to use

> the AMBER software suite to carry out simulations that can be run on a

> simple workstation in a reasonable period of time. They do not necessarily

> provide the optimal choice of parameters or methods for the particular

> application area."

>

> The main purpose of the tutorials here is to take people through the steps

> involved in a way that lets them appreciate (and hopefully think!!!) about

> what is going on and avoid use of any black box approaches as much as

> possible. They don't attempt to discuss the complexities of MD simulations

> which make MD such a dynamic area of research.

Yes I agree, although convenient ;-)

A new user generally follows what is written (if there is no clear
explanation) because she/he gets lost between all the possible
options. This tutorial is about "Charge derivation", and as it is a
complex topic it is crucial to well define the conditions in which the
computations are done. This is why in a R.E.DD.B. project all the
computation conditions are stored. See for instance the W-46 project @
http://q4md-forcefieldtools.org/REDDB/up/W-46/ & the "Information
regarding Quantum Calculations" section.

I have a last comment about this tutorial; the worst one:
Usually, one uses (an) intra-molecular charge constraint(s) for charge
value(s) for the atom(s) you _keep_ in the target fragment; NOT for
atoms that are removed. The tutorial reports the opposite (I think):
This arbitrary increases the RRMS for nothing. Usually, one also uses
a global charge constraint for the chemical group to be removed, and
not individual constraints otherwise one might observe artifacts in
the fit and/or an increase of the RRMS value.
http://pubs.acs.org/doi/abs/10.1021/ja964372c is one interesting
reference illustrating this point...

> Note that Antechamber (from the Antechamber source) uses:

>

> #HF/6-31G* SCF=tight Test Pop=MK iop(6/33=2) iop(6/42=6) opt

>

> Which includes the 6/42 option although with a value of 6 instead of 17.

Whatever the choices made - as I already said there is not a unique
answer, explaining those choices in a tutorial dedicated to charge
derivation is suitable - in particular if one does not follow the
standard procedure.

> I'm not sure the SCF Tight comment is strictly true here.

SCF(Conver=8) in geometry optimization & SCF(Conver=6) in single point
computation: Yes.

This is indeed "strictly" true. Using SCF(Conver=6) one generates MEP
with identical number of MEP points than with SCF(Conver=8). We have
done many tests to set the SCF keyword to SCF(Conver=6) in Gaussian
[but also to set the corresponding keywords in GAMESS-US &
PC-GAMESS/Firefly; more than 200 calculations on organic molecules
using different molecules/conformations/orientations].

> Note Gaussian (at

> least in Gaussian 98) uses a very loose SCF convergence criteria for single

> point SCF calculations(something like 10^-4 if I remember, I don't have the

> manual to hand) which will certainly affect the results. I have not checked

> to what extent it effects the final RESP charges though. It certainly can't

> hurt to include SCF=Tight - it will just up the computational requirements.

For Water-3 atoms- & HF/6-31G* ? Sure.

But for complex force field topology database building using
B3LYP/cc-pVTZ or aug-cc-pVTZ, the gain is several hours.

>> I would also suggest you to use HF/6-31G* (or HF/6-31G**, Duan et al.

>> FF) in the geometry optimization step - if you do want to rigorously

>> follow what has been done (& not MP2/6-31G* or B3LYP/6-31G* as it is

>> described in this tutorial).

>

> Are you sure HF/6-31G* was used for geometry optimization in FF94 and

> subsequently FF99, 99SB etc? The original FF94 JACS paper is awful when it

> comes to describing exactly what was done but I find it hard to believe that

> HF/6-31G* was used for the optimization.

Yes, you are right this is difficult to not be lost. However, once
again, in the context of writing a tutorial, choices made in that
tutorial should be clearly established to guide new users in all the
possible options.

> Note Cornell et al. JACS 1993, 115,

> 9620 talks about MP3/6-31+G**//HF/6-31G* and MP2/6-31G*//HF/6-31G* and MM2

> but never mentions HF/6-31G*//HF/6-31G*. I don't have the FF03 paper to hand

> to check what Yong did (no internet connection at the moment).

Duan et al. use HF/6-31G** in the geometry optimization step before
MEP computation (not HF/6-31G*); you can trust me ;-)

> I am

> intrigued to know where you get the recommendation for HF/6-31G* geometry

> optimization from.

If one has as objective to follow a standard procedure or to develop a
tutorial, I would differentiate the choices done in some papers, and
that that has been done in the AMBER force field topology database
building for instance, whatever the fact that one theory level is
better than another one (in some cases).

>> > To get a feeling of whether this procedure is correct I had a look

>> > at /usr/local/amber9/examples/resp_charge_fit/water. I created a

>> > water molecule in molden and applied the same procedure (this time

>> > with and without option iop(6/41=10)). The obtained charge are

>> > basically identical:

>> > O1 -0.81327

>> > H1 0.40664

>> > H2 0.40664

>

> As I thought - upping the Gaussian precision does not impact the results.

Please, consider this water molecule case (-3 atoms-) quite basic as
water should have a single "Standard orientation" in Gaussian.

>> I looked at your problem & computed RESP or ESP charges for water using

>> HF/6-31G*//HF/6-31G* (Cornell at al. FF, 1994-...)

>> or the olds:

>> HF/STO-3G//HF/6-31G* (Weiner at al. FF, 1984/1986)

>> HF/STO-3G//HF/STO-3G (Weiner at al. FF, 1984/1986)

>

>> This means the data available @

>> /usr/local/amber10/examples/resp_charge_fit/water/ have been generated

>> using HF/STO-3G//HF/STO-3G !

>> I guess they come from Amber... 4/3 ?. Do not use that as a reference !

>

> So if you use HF/6-31G*//HF-6-31G* can you reproduce the charges in the FF94

> paper?

What does it mean "reproduce" ? with what error ? None of the charge
values in the AMBER force fields are "reproducible": The maximum
differences for a conformation and a fitting approach is 0.06-0.07 e.
If the fitting approach or the conformation is different the
differences are really larger (how much ? impossible to say, but
terribly disturbing...).

If you use R.E.D. the charge reproducibility achieved is +/-.0001
whatever the method used is, whatever the QM program used, & whatever
the initial structure is. Please, use R.E.D. & see. The charge values
generated using R.E.D. & available in R.E.DD.B. are highly
reproducible: max diff = 0.0001!

> All the force field parameter derivations, like bond lengths and angles were

> done using MP2/6-31G* (same for GAFF) so it makes sense that they would use

> the same for charges, although there is the discussion that HF/6-31G*

> overestimates the gas phase dipole moment which is beneficial for solution

> phase simulations but from what I can see the discussion was never very

> clear what geometry this referred to.

See above.

> Perhaps someone should volunteer to provide a tutorial on the AMBER website

> that goes through a step by step process for reproducing the FF94 / FF99 /

> FF99SB and FF03 parameters (charges, VDW and valence terms) for say Alanine.

> I think this would really help clear up lots and lots of confusion that

> exists in the "black art" of force field derivation.

Concerning charge derivation, I think you have more than that @
http://q4md-forcefieldtools.org/Tutorial/ - taking the
dimethyl-alanine dipeptide as an example. You will find answers about
the impact of the molecular orientation/conformation/fitting option on
charge values with full R.E.D. jobs/data to be able to re-do the jobs.
Charge derivations are available for whole molecules but also for
different molecular fragments. The choices of using specific charge
equivalencing approaches and specific charge constraints are also
described, the goal being not just "how to get charges" but "how to
get charges with the lowest RRMS value".

I hope this helps.

regards, Francois

-----------------------------------------------------------------------
The AMBER Mail Reflector
To post, send mail to amber_at_scripps.edu
To unsubscribe, send "unsubscribe amber" (in the *body* of the email)
      to majordomo_at_scripps.edu