RATTLE errors in equilibration

From: regafan_at_usc.es
Date: Wed Nov 29 2006 - 09:46:32 CST

Hello,
I am trying to simulate a protein using NAMD, but I am having a lot of
problems with the RATTLE error. I have looked at the NAMD-L archive,
but it did not help. I have tried a lot of things, and at this moment,
I don't know what to do...Please, couls somebody help me?

What I did was the following:

-First, from a XR protein structure, I introduced the hydrogens with
xLeaP and minimized the whole system with AMBER (I had tried first
introducing the hydrogens with "guesscoord" and VMD, but the resultant
structure involved deformations related to the hydrogens added). I
create a psf and a pdb for using with NAMD and I solvated and ionized
this structure with the corresponding modules of VMD.

-The resultant structure is minimized with NAMD in 3 steps:
      -> First, only water and ions minimization (5000 steps)
      -> Second, minimization with the backbone of the proteine fixed
(5000 steps)
      -> Third, minimization of the whole system (10000 steps)
The final structure has a stable energy, I show you one of the inputs
of the minimization:
----------------------------------------------------------------------------
                            EXAMPLE MINIMIZATION
---------------------------------------------------------------------------
#protocol parameters
numsteps 5000
cutoff 9
switching off
binaryoutput off
binaryrestart no

#minimization parameters
minimization on
minTinyStep 1.0e-6
minBabyStep 1.0e-2
minLineGoal 1.0e-4

#fixed atoms parameters
fixedAtoms on
fixedAtomsFile proteina_fijada.pdb
fixedAtomsCol B

# files de entrada
  structure ionized.psf
  paraTypeCharmm on
  parameters par_all22_prot.inp
  coordinates ionized.pdb
  outputname proteina_min1
  restartname proteina_min1
  exclude scaled1-4
  1-4scaling 1.0
--------------------------------------------------------------------------------

Then, in a try to equilibrate the system, I heat it gradually
mantaining the skeleton of the protein fixed and constraints in the
rest of the atoms. Here you have my input file:

--------------------------------------------------------------------------------
                      INPUT FILE FOR EQUILIBRATION (1)
--------------------------------------------------------------------------------
rigidBonds all
rigidTolerance 0.0005

numsteps 40000
stepspercycle 20
firsttimestep 0
timestep 1
restartfreq 1000
DCDfreq 1000
DCDfile proteina_md1_dcd
outputEnergies 1000
temperature 0
reassignFreq 1000
reassignTemp 10
reassignIncr 10
reassignHold 310

cutoff 10
pairlistdist 11
switching on
switchdist 9
binaryoutput no
binaryrestart no

wrapAll on
wrapNearest on

#Temperature coupling
langevin off
#langevinDamping 10
#langevinTemp 310
#langevinHydrogen no

  PMEGridSizeX 120
  PMEGridSizeY 180
  PMEGridSizeZ 180
  margin 5

#molecular system
  structure ionized.psf
  #forcefield
  paraTypeCharmm on
  parameters par_all22_prot.inp

  # files de entrada
coordinates proteina_min3.coor
outputname proteina_md1
restartname proteina_md1
exclude scaled1-4
1-4scaling 1.0

#pressure control:
useGroupPressure yes
useFlexibleCell yes
useConstantRatio yes
langevinPiston on
langevinPistonTarget 1.01325 # in bar -> 1 atm
langevinPistonPeriod 200
langevinPistonDecay 200
langevinPistonTemp 310

#fixed atoms parameters
fixedAtoms on
fixedAtomsFile proteina_min3_backbone.pdb
fixedAtomsCol B

#constraints
constraints on
consref proteina_min3_99.pdb
conskfile proteina_min3_99.pdb
conskcol B

--------------------------------------------------------------------------------

After this, I turned on langevin, and decreased the constraints in the
protein, here you have an example of one of this inputs:

-------------------------------------------------------------------------------
                           INPUT FILE FOR EQUILIBRATION (2)
-------------------------------------------------------------------------------

rigidBonds all
rigidTolerance 0.0005

numsteps 50000
stepspercycle 20
firsttimestep 40000
timestep 1
restartfreq 1000
DCDfreq 1000
DCDfile proteina_md2_dcd
outputEnergies 1000
cutoff 10
pairlistdist 11
switching on
switchdist 9
binaryoutput no
binaryrestart no

wrapAll on
wrapNearest on

#Temperature coupling
langevin on
langevinDamping 10
langevinTemp 310
langevinHydrogen no

  PME on
  cellOrigin 3.8 -20.6 7.6
  cellBasisVector1 90 0 0
  cellBasisVector2 0 140 0
  cellBasisVector3 0 0 125
  PMEGridSizeX 120
  PMEGridSizeY 180
  PMEGridSizeZ 180
  margin 5
#molecular system
  structure ionized.psf
  #forcefield
  paraTypeCharmm on
  parameters par_all22_prot.inp

  # files de entrada
coordinates proteina_md1.coor
velocities proteina_md1.vel
outputname proteina_md2
restartname proteina_md2
1-4scaling 1.0

#pressure control:
useGroupPressure yes
useFlexibleCell yes
useConstantRatio yes
langevinPiston on
langevinPistonTarget 1.01325
langevinPistonPeriod 200
langevinPistonDecay 200
langevinPistonTemp 310

#fixed atoms parameters
fixedAtoms on
fixedAtomsFile proteina_md1_backbone.pdb
fixedAtomsCol B

#constraints
constraints on
consref proteina_md1_50.pdb
conskfile proteina_md1_50.pdb
conskcol B

-------------------------------------------------------------------------------

The problem appears in the next step, where the constraints in the
protein are "25". I get error of the type:

ERROR: Constraint failure in RATTLE algorithm for atom 10493!
ERROR: Exiting prematurely.

Changing the dimensions of the box, does not help, only changes the
atom involved in the error. Besides, in some of the cases, I can see
DEFORMATIONS in the residues of the protein, perhaps the source of the
error (but why are these atoms deformed). In another cases, the atoms
affected by the error do not present any kind of deformation.

Please, does anybody idea of what can be happening?

Thanks a lot in advance.

Rebeca
Postdoctoral student
Parc Cientific de Barcelona
regafan_at_usc.es

This archive was generated by hypermail 2.1.6 : Wed Feb 29 2012 - 15:44:13 CST