Hello ambers,
I have been struggling with the best (most correct?) way to do
temperature control for some coupled potential runs. Let me set
things up:
The system is a DNA-enzyme complex embedded in a 35A radius sphere
of waters, fully neutralized with Na+ ions. The solute structures
were taken from the last frame of a sequential series of fairly
long regular MD runs (ntt=0,temp0=300.0) using FF99 (approximately
5ns, total). I am using amber8, hence sander.QMMM is the coupled
potential program. There are 81 atoms in the QM set, needing 4
link atoms at the sites where the DNA and two residues of the
enzyme join to the rest of the system. For those who have not
tried coupled potential runs, sander.QMMM cannot do PBC or use any
of the implicit solvent treatments. In particular, we cannot use
the Poisson-Boltzmann model. So the boundary of the sphere is
vacuum with no treatment of the long range electrostatic effects.
I have done simulations two basic ways: with fcap=10.0 and no
restraints on the system, and with fcap=0.0 and restraints on all
the residues > 15A from the center of the active site and
RESTRAINT_WT=10.0. The center of the active site is near, but not
identical to the center of the sphere; next time I will make them
identical when I cut out the sphere of waters.
The problem is that the temperature will not stabilize at 300K
using ntt=1. ntt=2 or 3 gave similar behaviors. If I use a 0.5fs
time step (longer gives vlimit problems), a linear ramp of 10000
or 15000 steps, tempi=1, temp0=300, and tautp=0.5, the temperature
comes to a steady state at approximately 330K. If I relax tautp
to 5.0, the temp continues to rise. This behavior is illustrated
in the attached graph. (Since attachments are problematic on
mailing lists, you can go to
http://siegfried.utmb.edu/bdodson/temp1.pdf to pick it up.) I
have put the sander.QMMM input which generated this graph at the
end of the email. The system had been very briefly minimized (5
steps) before this MD run, really just to get an appropriate
restart file for this run.
I am assuming that this behavior is due to the electrostatics at
the vacuum-sphere boundary. Does this sound right? I am less
interested in the dynamical behavior of the system than in the
ensemble of structures represented in the trajectory. Would it be
appropriate to just set temp0 at 270, tautp at 0.5 and get on with
my simulations? Or is something else needed? Any suggestions or
critiques would be appreciated.
Thanks in advance,
Bud Dodson
-----------------------cut-8<-here------------------------
cap, equil-1 QMMM dynamics without D&C qm calculation
&cntrl
nmropt = 1,
ntx = 1, irest = 0,
ntpr = 10, ntwx = 100, ntwr = 100,
ntb = 0,
cut = 8.1,
ntc = 2, ntf = 2,
nstlim = 45000,
dt = 0.0005,
ig = 71277,
ntt = 1, gamma_ln = 0, tautp = 0.5,
igb = 0, fcap = 0.0,
tempi = 1.0, temp0 = 270.0,
ntr = 1, RESTRAINT_WT = 10.0, RESTRAINTMASK = ".2472,2477 > :15.0",
ifqt = 1, nqt = 81, idc = 0,
/
1 2 3 4 5 6 7 8 9 10 11 12 13 14
356 357 358 359 360 361
2465 2466 2467 2468 2469 2470 2471 2472 2473 2474 2475
2476 2477 2478 2479 2480 2481 2482 2483 2484 2485
2486 2487 2488 2489 2490 2491 2492 2493 2494 2495 2496
2497 2498 2499 2500 2501 2502 2503 2504 2505 2506 2507
2508 2509 2510 2511 2512 2513 2514 2515 2516 2517
2518 2519 2520 2521 2522 2523 2524 2525
&wt
type='TEMP0', istep1=0, istep2=15000,
value1=1.0, value2=270.0,
/
&wt
type='TEMP0', istep1=15000,istep2=45000,
value1=270.0, value2=270.0,
/
&wt
type='TAUTP', istep1=0, istep2=30000,
value1=0.5, value2=0.5,
/
&wt
type='TAUTP', istep1=30000,istep2=45000,
value1=0.5, value2=5.0,
/
&wt
type='END',
/
&rst
iat=0,
/
END
END
-----------------------cut-8<-here------------------------
--
M. L. Dodson bdodson.scms.utmb.edu
409-772-2178 FAX: 409-772-1790
-----------------------------------------------------------------------
The AMBER Mail Reflector
To post, send mail to amber.scripps.edu
To unsubscribe, send "unsubscribe amber" to majordomo.scripps.edu
Received on Mon Aug 16 2004 - 20:53:00 PDT
