Hello all,
I am attempting to follow an explicit solvent NMR refinement protocol
described by Xia, et al. (J Biomol NMR, 22: 317–311, 2002).
The protocol involves a slow heating phase, followed by two constant
pressure simulations with increasing pressure relaxation time (0.2 ps
in the first and 1.0 in the second). I have found that my solvated
protein consistently blows up at restart, following successful
completion of the first constant pressure simulation (TAUP=0.2). The
target temperature is 300 K, but in the first steps following restart,
the temperature soars to well over 800 K, then sometimes recovers due
to weak harmonic restraints on the protein to its starting structure.
I have seen this behavior regardless of which NMR model I use (have
tried 4 out of 10 models previously refined in GB solvent), regardless
of the pressure relaxation time used in the second NPT simulation
(have tried TAUP=1.0, 0.4, 0.2 and 5.0), regardless of timestep length
(2 or 1 fs), and regardless of the thermostat used. Following the
Xia, et al. protocol, I have used Berendsen temperature control
mainly, but I also tried a Langevin thermostat while troubleshooting.
I have checked my starting models for close contacts (using
Swiss-PdbViewer) and have found none. I checked also for water
molecules inserted by Leap into gaps in the folded protein but found
none. The only hints I have are as follows: 1) In cases where SHAKE
failure occurred, the atoms enumerated in the SHAKE error were in
residues near two different amino acids in non-default protonation
states (based on H++ predictions prior to the GB refinement); 2) the
density of the system during the first NPT simulation, with TAUP=0.2,
is lower than expected (~0.98 atm, where previous simulations with
this protein gave an equilibrated density of ~1.02 atm). I tried
recording every step for 1000 steps following the restart (1000 steps
into the second NPT simulation), but I could not see anything obvious
that could be causing the blow-up, as my harmonic restraints keep the
protein intact. Perhaps I will try repeating this process without the
restraints.
Does anyone have insight as to why the blow-up would consistently
occur at restart, even when using the same pressure coupling time
constant (TAUP=0.2) and when using the Berendsen thermostat? Is there
definitely a problem with my structure, or have I made an error in
parameterization? Below are my input files for the two NPT runs and
excerpts from the output files. I would be happy to supply further
information, if needed.
Thank you for your time.
Sally Pias
md1_NPT.in :
====
40 ps constant pressure MD at 300 K with weak restraints on protein
- Timestep is 1 fs
- Using default nonbonded cutoff (CUT = 8 Angstroms)
- Write to output file every 0.5 ps and to mdcrd file every 2 ps
- Pressure relaxation time (TAUP) is 0.2 ps (first cycle, following Xia 2002)
- Heat bath coupling constant is 1.0 ps (TAUTP)
&cntrl
imin = 0,
irest = 1, ntx = 5,
ntb = 2, pres0 = 1.0, ntp = 1, taup = 0.2,
ntr = 1, restraint_wt = 5.0, restraintmask = ':1-130',
ntc = 2, ntf = 2,
tempi = 300.0, temp0 = 300.0,
ntt = 1, tautp = 1.0,
nstlim = 40000, dt = 0.001,
ntpr = 500, ntwx = 2000, ntwr = 20000,
/
====
md2_NPT.in :
====
40 ps constant pressure MD at 300 K with weak restraints on protein
- Timestep is 1 fs
- Using default nonbonded cutoff (CUT = 8 Angstroms)
- Write to output file every 0.5 ps and to mdcrd file every 2 ps
- Pressure relaxation time (TAUP) is 1.0 ps (second cycle, following Xia 2002)
- Heat bath coupling constant is 1.0 ps (TAUTP)
&cntrl
imin = 0,
irest = 1, ntx = 5,
ntb = 2, pres0 = 1.0, ntp = 1, taup = 1.0,
ntr = 1, restraint_wt = 5.0, restraintmask = ':1-130',
ntc = 2, ntf = 2,
tempi = 300.0, temp0 = 300.0,
ntt = 1, tautp = 1.0,
nstlim = 40000, dt = 0.001,
ntpr = 500, ntwx = 2000, ntwr = 20000,
/
====
md1_NPT.out, last part of file:
====
------------------------------------------------------------------------------
check COM velocity, temp: 0.001720 0.00(Removed)
NSTEP = 40000 TIME(PS) = 160.000 TEMP(K) = 298.53 PRESS = -16.0
Etot = -62826.2928 EKtot = 15307.2402 EPtot = -78133.5330
BOND = 397.9615 ANGLE = 1129.7260 DIHED = 1422.9532
1-4 NB = 450.2222 1-4 EEL = 3405.7787 VDWAALS = 9825.9915
EELEC = -95069.2046 EHBOND = 0.0000 RESTRAINT = 303.0384
EAMBER (non-restraint) = -78436.5715
EKCMT = 6924.0838 VIRIAL = 7012.4230 VOLUME = 256073.5477
Density = 1.0001
Ewald error estimate: 0.9295E-04
------------------------------------------------------------------------------
A V E R A G E S O V E R 40000 S T E P S
NSTEP = 40000 TIME(PS) = 160.000 TEMP(K) = 300.24 PRESS = -13.7
Etot = -62660.2834 EKtot = 15395.0727 EPtot = -78055.3560
BOND = 395.9384 ANGLE = 1106.0447 DIHED = 1417.1549
1-4 NB = 444.9742 1-4 EEL = 3428.4524 VDWAALS = 9794.6675
EELEC = -94961.7665 EHBOND = 0.0000 RESTRAINT = 319.1783
EAMBER (non-restraint) = -78374.5343
EKCMT = 6934.5099 VIRIAL = 7015.8141 VOLUME = 258937.1979
Density = 0.9896
Ewald error estimate: 0.6336E-04
------------------------------------------------------------------------------
R M S F L U C T U A T I O N S
NSTEP = 40000 TIME(PS) = 160.000 TEMP(K) = 1.59 PRESS = 99.0
Etot = 274.2779 EKtot = 81.7028 EPtot = 251.8786
BOND = 15.8161 ANGLE = 21.9277 DIHED = 10.5646
1-4 NB = 7.3627 1-4 EEL = 21.0463 VDWAALS = 94.7523
EELEC = 254.3587 EHBOND = 0.0000 RESTRAINT = 10.8263
EAMBER (non-restraint) = 241.0523
EKCMT = 59.6722 VIRIAL = 575.3798 VOLUME = 6365.1963
Density = 0.0228
Ewald error estimate: 0.4788E-04
------------------------------------------------------------------------------
====
md2_NPT.out, first few hundred steps:
====
NSTEP = 500 TIME(PS) = 160.500 TEMP(K) = 849.32 PRESS = 2379.2
Etot = -9410.0858 EKtot = 43549.1015 EPtot = -52959.1872
BOND = 2031.7633 ANGLE = 5084.1456 DIHED = 2651.9693
1-4 NB = 890.0949 1-4 EEL = 3289.0518 VDWAALS = 10006.9733
EELEC = -86475.7602 EHBOND = 0.0000 RESTRAINT = 9562.5747
EAMBER (non-restraint) = -62521.7619
EKCMT = 19517.3019 VIRIAL = 4857.9887 VOLUME = 285373.0589
Density = 0.8974
Ewald error estimate: 0.1484E-04
------------------------------------------------------------------------------
check COM velocity, temp: 0.002222 0.01(Removed)
NSTEP = 1000 TIME(PS) = 161.000 TEMP(K) = 684.45 PRESS = 1223.4
Etot = -20517.6091 EKtot = 35095.5673 EPtot = -55613.1764
BOND = 2634.9845 ANGLE = 5954.0205 DIHED = 2542.3932
1-4 NB = 898.1212 1-4 EEL = 3244.1036 VDWAALS = 10108.4924
EELEC = -86862.8428 EHBOND = 0.0000 RESTRAINT = 5867.5509
EAMBER (non-restraint) = -61480.7273
EKCMT = 11148.2243 VIRIAL = 3361.0106 VOLUME = 294812.4374
Density = 0.8687
Ewald error estimate: 0.5358E-04
------------------------------------------------------------------------------
NSTEP = 1500 TIME(PS) = 161.500 TEMP(K) = 593.98 PRESS = 598.0
Etot = -29061.6347 EKtot = 30456.8660 EPtot = -59518.5007
BOND = 2512.6410 ANGLE = 5829.3987 DIHED = 2570.7454
1-4 NB = 820.6017 1-4 EEL = 3377.9469 VDWAALS = 10135.5203
EELEC = -88893.0697 EHBOND = 0.0000 RESTRAINT = 4127.7151
EAMBER (non-restraint) = -63646.2158
EKCMT = 8550.4287 VIRIAL = 4669.1921 VOLUME = 300617.6259
Density = 0.8519
Ewald error estimate: 0.3996E-04
------------------------------------------------------------------------------
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Fri Jul 10 2009 - 01:09:30 PDT
