AMBER Archive (2009)

Subject: Re: [AMBER] unexplained blow-up at restart

From: Carlos Simmerling (carlos.simmerling_at_gmail.com)
Date: Fri Jul 10 2009 - 05:55:51 CDT

can you send the earlier part of run 2? it looks here like there are already
serious problems, maybe if we see what happens earlier it will help. maybe
setting ntpr to something smaller (even ntpr=1, and set nstlim to 1000 or
so) will help debug it, you might get more info on what goes wrong first.
for example, is the restraint energy really high already at step 1? compare
the energies at the end of run1 with those at the very start of run 2, see
if any terms differ.

On Fri, Jul 10, 2009 at 4:02 AM, Sally Pias <sallypias_at_gmail.com> wrote:

> 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
>
_______________________________________________
AMBER mailing list
AMBER_at_ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber