AMBER Archive (2008)

Subject: AMBER: QM/MM Heating

From: Steven Winfield (saw44_at_cam.ac.uk)
Date: Mon Jan 14 2008 - 09:22:27 CST

Dear all,

I'm having a look at AMBER 9's energy conservation when using its internal
QM/MM code. I took a simple system of pure water (1222 molecules in total
   in a periodic box of length 34.419A), with one water molecule
flagged as the QM region (using the PM3 Hamiltonian), turned off SHAKE for
this molecule and used a 0.2fs time step. The other molecules were TIP3P.
The input restart file had been equilibrated at 300K
using 100ps of NVT dynamics with the Langevin thermostat. The
QM/MM simulation lasted 1ns in total.

I saved the position and velocity data so that I could examine the
temperature of the MM and QM regions afterwards. I have attached a bzipped
tar with my input files, an output text file, and some temperature plots.
Measuring the temperature of a single molecule is obviously very noisy, so
I also calculated a weighted average of temperature values at each point
with weights exponentially decaying with a time constant of 50ps.

As can be seen from the plot (pme_temp_50.ps) and the text output file, I
observed a steady temperature increase from 300K to around 500K, and the
average temperature of the quantum water becomes quite large. I tried the
same simulation but using a full Ewald sum (ewald_temp_50.ps) and observed
similar results.

I ran a pure MM simulation with the same input files for comparison (with
SHAKE on all molecules),
because page 143 of the manual says that QM/MM with default
parameters should "generally conserve energy about as well as one would
find for a corresponding pure MM simulation". I increased the time step by
a factor of 10 to 2.0fs, which I thought would be on the very edge of
becoming unstable, but the temperature plot (water_mm.ps) shows no signs
of heating.

So my questions are:

What non-default parameters have I used, if any?
I've had others look at my input files and the only comment I've had is
that my density may be slightly low. I've kept my timestep small and my
cutoffs large (9.0A)
Is there a problem with AMBER's QM/MM code? I'm using a fully patched
version.

Any comments would be most helpful.

Regards,
Steven Winfield.

P.S. The input restart file is called 'zeroed' because I wanted initially
to see if energy was being successfully transferred from the rigid MM
molecules to the non-rigid bonds of the QM molecule, so I manually zeroed
the velocities of the QM molecule before starting. This shouldn't affect
the results though.

-----------------------------------------------------------------------
The AMBER Mail Reflector
To post, send mail to amber_at_scripps.edu
To unsubscribe, send "unsubscribe amber" to majordomo_at_scripps.edu