[gmx-users] about Normal Mode Analysis

[Mark Abraham]

Anna Marabotti wrote:
> Dear All,
> I'm trying to perform my first Normal Mode Analysis on a protein of 158 residues, therefore I read the Gromacs
> manual searching for the proper way to do it.
> I understand that the first step is a deep minimization of my protein, but when I tried to apply it to my
> system, it failed to reach the required convergence. I tried several times without success. I'm using Gromacs
> v. 3.3.3 which is compiled in double precision (I also checked in the .tpr file and it gives to me confidence
> that the precision is really double, since coordinates are written with six decimals).

That doesn't strike me as enough for double precision. I'd be expecting 
more like 12 digits of precision for double. However this is not 
necessarily the problem.

> I first made a steepest
> descent EM, then I used the output .trr file as an input for grompp (the command was grompp_d -f emdeep.mdp -c
> em1st.gro -t em1st.trr -o emdeep.tpr). My emdeep.mdp file was:
> 
> title               =  emdeep 
> cpp                 =  /usr/bin/cpp
> define              =  -DFLEXIBLE (however, I minimized the protein only, in vacuum)
> constraints         =  none
> integrator          =  cg
> nsteps              =  20000
> emtol               =  0.00001 
> emstep              =  0.01, then 0.05, then 0.1
> nstlist             =  5 
> ns_type             =  grid
> rlist               =  1.0
> coulombtype         =  PME
> rcoulomb            =  1.0
> rvdw                =  1.0
> fourierspacing	=  0.12
> fourier_nx		=  0
> fourier_ny		=  0
> fourier_nz		=  0
> pme_order		=  4
> ewald_rtol		=  1e-5
> optimize_fft		=  yes
> 
> The run did not reach convergence to the requested Emtol, the Fmax was oscillating around a value of 100.
> Since the program claimed (as usually) that the stepsize was too short, I enlarged it, but nothing changed
> with cg. I changed the integrator with l-bfgs, and it continued until 20,000 steps, without success (I noticed
> that the Fmax was more or less the same during the whole minimization, so I don't think that it would be
> useful to prolonge further the simulation). I also read the messages in the gmx-users list, but did not find
> useful infos to solve my problem.

The possibility exists that you are in a local minimum that lacks a 
point for which the gradients are as low as you require for the step 
size available to the algorithm. emtol (see manual) is some kind of 
maximum gradient (i.e. force), and the value you've chosen is extremely 
small. I've no idea if that tolerance is consistent with what people do 
before NMA - check literature. Increasing the step size eventually lets 
you escape the current local minimum, but there's no guarantee of 
reaching a suitably flat minimum, since the stepsize may then be too 
large to stay within some potentially suitable potential well.

In extremis, run a short equilibration and try minimizing from a few 
points taken from the run.

> I would like to ask you two things. First of all, obviously, if you could help me to find what's wrong in my
> settings and/or in my system. Then, my main doubt is how to set properly the .mdp file for minimization. Could
> you inform me if there are somewhere tutorials related to NMA using GROMACS? I think that such kind of
> exercise would be very useful: in fact, I think that it is not easy to translate the theory of NMA in practice
> when you have to deal with all the parameters present in a .mdp file, and an example would serve as a starting
> guideline.

Other than doing a literature search for reported methodology for 
GROMACS NMA, I can only suggest 
http://wiki.gromacs.org/index.php/normal_mode_analysis, which was based 
on a post on this list, IIRC.

Mark