Re: AMBER: Thermodynamic Integration for conformational changes?

From: Francesco Pietra <chiendarret.yahoo.com>
Date: Tue, 21 Aug 2007 23:53:15 -0700 (PDT)

--- Gustavo Seabra <gustavo.seabra.gmail.com> wrote:

> On 8/21/07, Francesco Pietra <chiendarret.yahoo.com> wrote:
> > I have carried out QM-MM (AM1) simulations for different conformations of a
> > 98-atoms non-polymeric neutral molecule (solute) in explicit solvents.
>
> Unfortunately, that doesn't tell us much. What kind of calculation did
> you do? Did you use Amber? Was it a minimization or molecular
> dynamics? Or was it just a single point calculation for each
> structure?

Started from EQE and EQA conformers (conformational inversion about chirality
axes changes EQE into EQA, and viceversa; the two conformers are in
diastereomeric relationship and are isolatable). Energy difference measured
from NMR populations (EQE/EQA 90:10 at r.t. in CHCl3). Very extensive MMFF94
molecular mechanics and simulated annealing by Amber, both in vacuum, agree
for EQE and EQA being the absolute minima for each diastereomer. Simulated
annealing at 1000K interchanges the conformers - and gives also minor ones -
though it gives wrong relative energies. Nonetheless, the molecular shapes
grossly agree with what can be drawn from NMR.

After minimization as from min.in below, I carried out sets of MD (for ca.
160ps up to now) with the following command, where md.in is specified below.
That for the couple of conformers in either TIP3PBOX or CHCL3BOX. The gross
integrity of both EQE and EQA is conserved during MD.

~$ mpirun -np 4 $AMBERHOME/exe/sander.MPI -O -i md_20000_qmmm.in -o EQE_md3.out
-p EQE.prmtop -c EQE_md2.rst -r EQE_md3.rst -x EQE_md3.mdcrd

min.in:

Initial min of our structure QMMM nano written
&cntrl
imin=1, maxcyc=500, ncyc=200,
cut=8.0, ntb=1, ntc=2, ntf=2,
ifqnt=1
/
&qmmm
qmmask=':1',
qmcharge=0,
qmtheory=2,
qmshake=1,
qm_ewald=1, qm_pme=1
/

md.in

300K constant temp QMMMMD
 &cntrl
  imin=0, ntb=1
  cut=8.0, ntc=2, ntf=2,
  tempi=300.0, temp0=300.0,
  ntt=3, gamma_ln=1.0,
  nstlim=20000, dt=0.002,
  ntpr=100, ntwx=100,ifqnt=1
 /
 &qmmm
  qmmask=':1',
  qmcharge=0,
  qmtheory=2,
  qmshake=1,
  qm_ewald=1, qm_pme=1
 /



After the 160ps, averaged AM1ESCF (EPtot) for EQE in H2O: -225.7 (-29983.3);
for EQA in H2O: -222.1 (-30480.4).

To my shame, I must confess that your indications below about ESCF/EPtot are
not clear to me. In particular, I do not understand the difference between "the
energy of the QM region, which is influenced by the presence of the MM atoms."
and "The energy of the MM atoms, however, is *not* included there." I.e, I do
not understand what is meant for "QM region" and "MM atoms". Has the ESCF
difference (3.6 kcal/mol) any meaning?


>
> > Would like to compare with DFTB, if I get a license for that.
>
> About that, Marcus Elstner is out of the country (his country) right
> now, so it's unlikely you'll get an answer from him soon. But you
> could also try to download the parameters from www.dftb.org. You just
> have to fill a registration form, then fax it to the number in the
> form. Then you'll get an username and password, which you can use to
> download the parameters. The specific set you want is 'mio-0-1', which
> is the 'organic' set. The others contain atoms not supported in
> Amber9.
>
> When you get the parameter set from them, you will need to rename the
> files. Their naming scheme is, for example, C-C.skf (Upper case,
> 'dash' and .skf extension). Then, all you have to do is to copy those
> files to the $AMBERHOME/dat/slko directory, and rename them to
> something like 'cc.spl' (lowercase, no dash, .spl extension), and that
> should do.
>
> > It is unclear to me how to relate above results to conformational energies,
> in
> > particular how to deal with ESCF. The results seem to suggest that
> > perturbations by the solvent cancel out by comparing ESCF for the
> conformers,
> > though the theory behind has not been examined.
>
>
> ESCF only gives you the energy of the QM region, which is influenced
> by the presence of the MM atoms. The energy of the MM atoms, however,
> is *not* included there. If you want to look at the energy of the
> *system*, what you really want is the total energy instead, which does
> include the solvent molecules.
>
> > At any event, I would like to calculate free energy changes for this
> system.
> > The question whether thermodynamic integration is appropriate for
> complicated
> > conformational changes received authoritative NO at Amber8's time
> (archives)
> > and is probably to be extended to Amber9. Are umbrella sampling or MM_PBSA
> > potential solutions? The first one, from the manual, also seems unable to
> deal
> > with multiple dihedral changes from one conformation to the other one. What
> > about MM_PBSA? Is any example to which to refer for the evaluation of
> absolute
> > free energies in explicit solvent by MM_PBSA?
>
> I'll not pretend to understand much here, and I'd advise you to take
> whatever I say with a whole pound of salt, but notice that you have a
> really big system and, for what I can imagine, some big conformational
> changes from one point to the other. Have you tried to imagine what
> would you use as a reaction coordinate? The big problem for you (I
> believe) is that for most free energy methods you need to define not
> only initial and final states, but you also need some idea about the
> reaction path you want to probe. AFAIK, there's no 'magical' method
> that will get initial and final point and locate the minimum free
> energy path between them. And that would basically be the reason those
> methods won't work for you: with that many dihedrals, how many
> different reaction paths can you think of? When can you be satisfied
> that you have finally found the lowest energy one? (But I can be wrong
> here... anybody?)
>

If I understand, the initial and final states are what I am interested in, at
least now. I would be happy enough for now to calculate relative potential
energies, and detailed conformational differences, for EQE vs. EQA under the
influence of explicit water and explicit chloroform (and other explicit
solvents). This would be of immediate use in letaion to CD spectra of these
diastereomers, and of future use for their interaction with pharmacologically
relevant receptors.

> In principle, you can also try for example a real long REMD
> simulation, until you see enough transitions beteween the two states
> that you can consider the populations to be equilibrated. (Follow the
> populations of the conformers until they stabilize) Then, get the free
> energy difference from the population of each conformer. But that
> would probably take longer than what you (or anyone I know) would be
> willing to wait if you are going to do that with QM.
>
> But I don't mean to be negative here. Hopefully, someone with more
> experience will chime in here with a better solution.


Thanks
francesco
>
> Good luck,
> Gustavo.
> -----------------------------------------------------------------------
> The AMBER Mail Reflector
> To post, send mail to amber.scripps.edu
> To unsubscribe, send "unsubscribe amber" to majordomo.scripps.edu
>



      ____________________________________________________________________________________
Shape Yahoo! in your own image. Join our Network Research Panel today! http://surveylink.yahoo.com/gmrs/yahoo_panel_invite.asp?a=7


-----------------------------------------------------------------------
The AMBER Mail Reflector
To post, send mail to amber.scripps.edu
To unsubscribe, send "unsubscribe amber" to majordomo.scripps.edu
Received on Sun Aug 26 2007 - 06:07:04 PDT
Custom Search