Constraint failure in RATTLE algorithm occurs during NAMD simulation with amber force field

From: sunyeping (sunyeping_at_aliyun.com)
Date: Thu Aug 29 2013 - 21:32:22 CDT

Dear all,
 
I am now trying do NAMD simulation of proteins with amber field field. The systems crash rapidly after several dozens of minimization steps. Three chains of the proteins break down apart and run out of the water box. The program stops just when the minimization and initiation of heating, and "ERROR: Constraint failure in RATTLE algorithm for atom......" appears. The discussion about RATTLE algorithm is very common on web, I have tried many suggestions put up in these disscusions, such as desceasing the timestep, increasing margin, and increasing the temperature gradurally during heating, but none of them works.
    I prepared the amberparm files, filename.crd and filename.top by tleap using the follow leap.in script:
 
source leaprc.ff03.r1complex = loadpdb complex_mod.pdbbond complex.101.SG complex.164.SGbond complex.203.SG complex.259.SGbond complex.300.SG complex.355.SGsavepdb complex complex1.pdbsolvateBox complex TIP3PBOX 12.0 0.75 addIons2 complex Na+ 0 saveamberparm complex complex_wb1.top complex_wb1.crd
savepdb complex complex_wb.pdb quit
 
And the NAMD configure file is as follow:
 
############################################################# ## JOB DESCRIPTION                                         ## #############################################################  # Minimization and Equilibration of  # H2PA in a Water Box   ############################################################# ## ADJUSTABLE PARAMETERS                                   ## #############################################################  amber              onparmfile          complex_wb.top ambercoor         complex_wb.crd  set temperature    310 set outputname     eq1_o  firsttimestep      0   ############################################################# ## SIMULATION PARAMETERS                                   ## #############################################################  # Input #paraTypeCharmm        on #parameters          par_all27_prot_lipid.inp  temperature         $temperature   # Force-Field Parameters exclude             scaled1-4 1-4scaling          1.0 cutoff              12.0 switching           on switchdist          10.0 pairlistdist        14.0   # Integrator Parameters timestep            2  ;# 2fs/step rigidBonds          all  ;# needed for 2fs steps nonbondedFreq       1 fullElectFrequency  2   stepspercycle       10   # Constant Temperature Control langevin            on    ;# do langevin dynamics langevinDamping     1     ;# damping coefficient (gamma) of 1/ps langevinTemp        $temperature langevinHydrogen    off    ;# don't couple langevin bath to hydrogens   # Periodic Boundary Conditions cellBasisVector1    86.1     0.0    0.0 cellBasisVector2     0.0     80.1   0.0 cellBasisVector3     0.0     0.0    81.3 cellOrigin          0.08574733138084412 0.287835955619812 0.16694538295269012 wrapAll             onmargin              3   # PME (for full-system periodic electrostatics) PME                 yes PMEGridSpacing      1.0  #manual grid definition PMEGridSizeX        90 PMEGridSizeY        90 PMEGridSizeZ        90   # Constant Pressure Control (variable volume) useGroupPressure      yes ;# needed for rigidBonds useFlexibleCell       no useConstantArea       no  langevinPiston        on langevinPistonTarget  1.01325 ;#  in bar -> 1 atm langevinPistonPeriod  100 langevinPistonDecay   50 langevinPistonTemp    $temperature  #Harmonic constrain if {1} { constraints           on consref               constr.pdb conskfile             constr.pdb conskcol              B constraintscaling     10.0 }   # Output outputName          $outputname  restartfreq         500     ;# 500steps = every 1ps dcdfreq             1000 xstFreq             1000 outputEnergies      1000 outputPressure      1000   ############################################################# ## EXTRA PARAMETERS                                        ## #############################################################   ############################################################# ## EXECUTION SCRIPT                                        ## #############################################################  # Minimization if {1} { minimize            5000reinitvels          $temperature }#Raising temperaturefor { set TEMP 10 } { $TEMP <= 300} { incr TEMP 5 } {langevinTemp $TEMPoutput md.$TEMPrun 500 ; # 2ps } run 500000 ;# 1000ps
 
The last energy output in namd log files:
 
 
ETITLE:      TS           BOND          ANGLE          DIHED          IMPRP               ELECT            VDW       BOUNDARY           MISC        KINETIC               TOTAL           TEMP      POTENTIAL         TOTAL3        TEMPAVG            PRESSURE      GPRESSURE         VOLUME       PRESSAVG      GPRESSAVGENERGY:     5000    132692.3181     88075.3164      9460.6170         0.0000        -158435.0636 -9999999999.9999  17717558.9656         0.0000         0.0000      -9999999999.9999         0.0000 -9999999999.9999 -9999999999.9999         0.0000      -3390108488.2031 -3390100133.8221    560694.3930 -3390108488.2031 -3390100133.8221
 
Could you help me with these? Thanks in advance!
 
 
 
Institute of Microbiology, Chinese Academy of Sciences

Institute of Microbiology, Chinese Academy of Sciences

This archive was generated by hypermail 2.1.6 : Wed Dec 31 2014 - 23:21:35 CST