[gmx-users] elastic network model - normal mode analysis

[Mark Abraham]

On 23/09/2012 5:22 PM, mohan maruthi sena wrote:
> Hi all,
>             I want to  do normal mode  analysis for a protein using elastic
> network model containing 691 residues. I consider only c-alpha atoms
> connected by a spring constant of 81600 kj/nm2.  When I do normal mode
> analysis I  get the following error :
>
> Maximum force: 6.30230e+02
> Maximum force probably not small enough to ensure that you are in an
> energy well. Be aware that negative eigenvalues may occur when the
> resulting matrix is diagonalized.
> Finished step 691 out of 691
>
> I have minimised the structure , still i get the same error. I use  in
> gromacs 4.5.5 double precision  for  these simulations. When I load the
> trajectory eigen.trr  in vmd(to see normal modes) , I could see the protein
> far away and all atoms overlapped. I use the following .mdp file parameters
> as input
>
> ; RUN CONTROL PARAMETERS
> integrator               = nm
> dt                       = 0.002 ; time step (ps)
> nsteps                   = 1000000 ; number of steps
>
>
> ; OUTPUT CONTROL OPTIONS
> nstenergy                = 500
> nstxout                  = 500
> nstvout                  = 500
> nstfout                  = 500
> energygrps               = System
>
>
> ; NEIGHBORSEARCHING PARAMETERS
> nstlist                   = 0      #Frequency to update neighbourlist
>
> pbc                       = no
>
> comm_mode                 = ANGULAR
>
> ; OPTIONS FOR ELECTROSTATICS AND VDW
>
> ; Temperature coupling
> tcoupl                   = v-rescale           ; Couple temperature to
> external heat bath according to Berendsen method
> tc-grps                  = System  ;Non-Protein ; Use separate heat baths
> for Protein and Non-Protein groups
> tau_t                    = 0.2      ;0.2         ; Coupling time constant,
> controlling strength of coupling
> ref_t                    = 300        ; Temperature of heat bath
>
>
> ; GENERATE VELOCITIES FOR STARTUP RUN
> gen_vel                  =  no   ; Assign velocities to particles by taking
> them randomly from a Maxwell distribution
> gen_temp                 = 300.0  ; Temperature to generate corresponding
> Maxwell distribution
> gen_seed                 = 9999   ; Seed for (semi) random number
> generation. Different numbers give different sets of velocities
>
> ; OPTIONS FOR BONDS
> constraints              = none ; No constraints except for those defined
> explicitly in the topology
>
> refcoord_scaling         = com
>
> Please suggest me a way to rectify this error.

You're doing something wrong in the construction of your model. Find a 
toy example in the literature that you can attempt to copy so that you 
know you know how to walk before you try to run. Failing that, construct 
some geometric solid that should be stable, like a tetragonal pyramid, 
or something. Until you can do that, worrying about a protein is a waste 
of time. The strength of your forces, temperature and size of your time 
step are all closely interlinked if you want the integration to do 
something sensible, and probably you're inventing all of them at the moment.

Mark