AMBER Archive (2009)

Subject: [AMBER] unexplained blow-up at restart

From: Sally Pias (sallypias_at_gmail.com)
Date: Fri Jul 10 2009 - 03:02:42 CDT

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_at_ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber