Re: [AMBER] Big potential difference between conformation in MD simulation and optimized conformation in Gaussian.

From: Jason Swails <>
Date: Mon, 17 Aug 2015 16:18:36 -0400

On Mon, Aug 17, 2015 at 3:08 PM, Eg eh <> wrote:

> Dear all,
> I am doing simulation with a polymer in TIP3P water now. The monomer has a
> 7membered ring. I used PARAMFIT to get a set of force field for it. I used
> TLEAP to combine monomers to a polymer and pack it with TIP3P water.
> After running for 20ns (NPT), I extracted some fragment from the
> trajectory file. These fragments were all monomer with missing terminals.

​You need to make both terminal and internal residue templates. This is
also done for amino acids, where N-terminal residues have a leading "N" in
front of them and C-terminal ones have a leading "C". For example, if you
wanted to make sure that you were simulating a zwitterionic tri-alanine
molecule in tleap, you would need to build the sequence as such:

triala = sequence {NALA ALA CALA}

If you just used ALA three times, the N-terminus would be missing 2 H's and
the C-terminus missing an oxygen. The leaprc files in
$AMBERHOME/dat/leap/cmd (e.g., leaprc.ff14SB) show the commands used to
tell tleap to map the N-terminal amino acids (and nucleic acids) to the
N-terminated (or 5'-terminated)​

​residue template.

For your monomer, you need to define charges for all termini (with
different residue names) as well as the internal residue (and set the HEAD
and TAIL atoms correctly to make sure tleap forms the correct intra-residue

> I added hydrogen atoms to makeup the terminals. And then I used Gaussian09
> to do SP and OPT with MP2/6-31G* to see the potential difference.
> The results showed the conformation extracted from MD simulation had much
> higher potential than its optimized conformation. I had 24 conformations.
> The biggest potential difference was 66kcal/mol. The smallest potential
> difference was 35kcal/mol.
> I know AMBER could overcome small energy barrier, but the potential
> difference was too big. Is this condition normal?

​No, crossing 35 kcal/mol barriers is not normal, but that is also probably
not what is happening here. MD will cross barriers based on the full
(temperature-dependent) free energy surface. The 35 kcal/mol barrier seems
to correspond to a gas-phase QM scan, whereas the MM simulation is using
the parameters you assigned *in explicit solvent*. It is *these* barriers
that influence how quickly you should expect conformations to interconvert,
not what the QM says.

​Note that gas-phase and solution-phase molecules often behave *very*
differently, so there's no reason to expect, a-priori, that a gas phase QM
calculation should be representative of what should happen in a
solution-phase MM simulation.


Jason M. Swails
Rutgers University
Postdoctoral Researcher
AMBER mailing list
Received on Mon Aug 17 2015 - 13:30:07 PDT
Custom Search