Michel Becker wrote:
> During my search in the archive, I found that it's strongly recommanded
> to use as well for the geometry optimisation 6-31G* (NOT the ESP
> calculation), but actually I am quite sure, that I've read once here
> that some people used higher ones for the optimisation.
Dear Michel,
I've often wondered this as well, I think I've asked this question of
people before and not got a particularly fulfilling answer. I think it
is a complex area. The following is what I think and have learnt so far
in my PhD. It may not be fully correct and I beg any onlist AMBER gurus
to chime in and correct anything that I may have wrong in my reasoning.
My own conclusion is that one should use the "best" QM method for the
optimisation of a target structure, but then at this optimised
structure, one MUST follow the charge fitting protocol for that
particular force field parameter set. For example, in FF94 one would use
a HF/6-31G* wavefunction (at this optimised structure) to generate the
Merz Kollman(MK) point charge grid for use in the sequential RESP fit.
Bayly et al's RESP paper:
J. Phys. Chem.;1993; 97(40); 10269-10280
http://pubs.acs.org/cgi-bin/searchRedirect.cgi/jpchax/1993/97/i40/pdf/j100142a004.pdf
stipulates this method and basis for the wavefunction. Deviation from
this would lead to a charge inconsistency with the other molecular
"units" within that force field . I think the HF/6-31G* wavefunction is
used with the RESP protocol since there is a fortunate over estimation
(~20%) of the in-vacuo dipole moment. This is also touched upon bottom
right on page 100 of:
http://dx.doi.org/10.1002/jcc.20157
This over estimation corresponds to the value of the dipole moment one
would get if carrying out the the calculation in condensed phase. Hence
exploiting this over estimation, the cheaper in-vacuo HF/6-31G* RESP
charges can be used to obtain charges that would be suitable for a
condensed phase biological simulation. Remember, during the development
of the RESP protocol in the early 90's, such a QM calculation would have
been quite time consuming in itself; the addition of Polarizable
Continuum models or even explicit solvent to the QM calculation would
have been infeasible at that point in time. Also see bottom of page 21
of AMBER 9 manual:
http://amber.scripps.edu/doc9/amber9.pdf
Remember, as Francois points out, some of the newer force field
parameters use a different wavefunction for their RESP fitting.
The whole point of the charges from the RESP protocol is to reproduce
the "correct" long range electrostatic behavior of a molecule in
condensed phase. This is perhaps the most important factor in a
biological MD simulation. The father of AMBER, the late Peter Kollman,
is described as impressing the importance of this point:
http://www.biophysj.org/cgi/content/full/81/4/2422 (beginning of 2nd
paragraph)
Ok, now I'm (slowly) getting to the point here, looking at the
evaluation of the AMBER force field function, the charges from the RESP
fit are only used in the non-bonded electrostatic term and in the
dihedral term.
Therefore the interaction between two atoms that are classified as
bonded in the model's topology is dominated by the bond parameters for a
bond of those two atom types. Ideally, covalent parameters such as the
equilibrium bond length and force constant for such a term should be
taken from the best QM optimisation or experiment.
But this is where my argument falls apart and the email starts to go a
bit offtopic. :( Each of the atoms defined as being bonded also
experience energy contributions from some or all other terms in the
AMBER function. If one was to optimise to an energy minimum using solely
the AMBER function for a parameterized system, one may end up with a
structure that is not the same as experiment/high level QM due to
additions from other terms. One may even need to tweak either of the
bond parameters to be incorrect, so that correct structure is lowest
AMBER energy conformation.
The situation is more complex with dihedrals since there is a scaled
RESP charge term in it. (via a factor of 1.2 (read: multiplied by
0.833333) )
There seems to be a dependency chain with each line depending on the
line above it:
Method of Optimisation
^
|
Optimised structure
^
|
RESP Charges
^
|
Dihedral terms <-> bonded terms <-> (+other terms?)
Hence all the terms within the parameter site are dependant on each
other which in turn depend on the RESP charges.
The only way I see of developing missing parameters in a non-standard
residue for a specific force field parameter set is to do the following:
1) Start from an experimental conformation of the non-standard residue.
2) Fit RESP charges using the charge fitting protocol of that force
field parameter set.
3) Assign all parameters that can be assigned within the "vocabulary" of
the force field parameter set and consider these as FIXED. VDW
parameters will also be considered FIXED.
4) Take reasonable guesses at any missing parameters and consider these
as VARIABLE.
5) Generate an ensemble of non-standard residue structures by running
some GB MD at a temperature that is relevant to biological processes
(i.e. ~298K). In this way, the ensemble is exploring the region of phase
space where it will be used, hence the eventual fitted parameters will
be most relevant.
6) Post process the trajectory, evaluating each structure's QM energy
using the "best QM method available"
7) Evaluate each structure's MM energy using the AMBER parameters
8) Start to adjust the VARIABLE parameters so that the majority of the
ensembles AMBER energies match the QM energies for each structure.
9) Goto 5, running the MD with the new generation of VARIABLE parameters.
10) Stop if a pre-defined convergence in VARIABLE parameters in reached.
The above recipe is based upon an off list informal chat ages ago with
Ross Walker.
Returning to the topic, one could actually do the experiment (which I
should actually do myself) whereby one carries the out default protocol
RESP fit charges on two target structures of the same system; one having
been optimised at say, MP4/{some very large basis set} and the other at
HF/6-31G*. This will give an indication of the RESP charge's sensitivity
to method of optimisation.
In addition I recall a recent thread along these lines as well:
http://amber.ch.ic.ac.uk/archive/200611/0003.html
I started writing this email last night, and ended up getting quite
confused by the end of it all. Then this morning, I saw Francois'
response to this thread and it is much more concise than mine. I'm going
to post what I've written anyway; and it will display my ongoing
confusion of the topic. Anyway.... I'm not sure if that's helped or just
clouded the water, but it's starting point for discussion.
regards,
Mark Williamson
http://dumb.ch.ic.ac.uk/~mjw99/
-----------------------------------------------------------------------
The AMBER Mail Reflector
To post, send mail to amber.scripps.edu
To unsubscribe, send "unsubscribe amber" to majordomo.scripps.edu
Received on Sun Mar 25 2007 - 06:07:14 PDT
