normal mode analysis
A normal mode analysis was made with GROMACS 3.3.1. grompp,mdrun and g_nmeig
were used in double precision.
Two systems were tested: BPTI in SPC and Protein/DNA in vacuum
Therefore, energy was minimized twice with cg/cg and steep/cg respectively.
Analyzation of the Hessian ALWAYS yielded negative eigenvalues, which was
Both systems were then tested with the same mdp-parameters in GROMACS 3.2.1 and
yielded good results (compared to BPTI with 3.1.4) and not so good results
(first 6 eigenvalues not zero in the DNA-system, due to bad energy minimization(?)).
Then, a nm-tpr-file was generated with nm.mdp and grompp_d from 3.2.1, which was
then processed by mdrun_d and g_nmeig_d from 3.3.1. This also just yielded
Conclusion: Bug either in mdrun or g_nmeig. I wasn't able to go deeper, cause
g_nmeig_3.3.1 can't read the Hessian generated by 3.2.1 and vice versa.
Attached is a bugreport.tar.gz, which contains:
root: topol.top, em.mdp, nm.mdp and confout.gro (top and gro from the BPTI-System).
321: em: topol.tpr(3.2.1); nma: topol.tpr(3.2.1), eigenval.xvg(3.2.1)
331: em: topol.tpr(3.3.1); nma: topol.tpr(3.3.1), eigenval.xvg(3.3.1)
where em contains the system for minimization and nma that for normal mode analysis.
I didn't want to include the Hessians because of the file-size.
#2 Updated by Erik Lindahl about 14 years ago
This should be enough for us to reproduce it. If you don't have any problems with making a small hack
and recompiling you could help us fix it faster by trying to force the full hessian representation if the
system is reasonably small:
In file src/mdlib/minimize.c, around line 1850, we set the variable bSparse to TRUE if the system
contains more than 1000 atoms and cutoffs are used. If you set bSparse to FALSE after these
statements we will use a different representation and standard LAPACK diagonalizer instead of ARPACK
sparse matrix diagonalization. After editing, just go to the top gromacs directory and rerun "make".
Alternatively I'll try to have a look next week!
#3 Updated by David van der Spoel over 13 years ago
I tried what you suggested and it did not help. However the energy minimization
was not good:
Maximum force: 2.34618e+02
Maximum force probably not small enough to ensure that you are in a
energy well. Be aware that negative eigenvalues may occur when the
resulting matrix is diagonalized.
Nevertheless it gives more reasonable values in 3.2.1.
#4 Updated by Erik Lindahl over 13 years ago
OK, that force is so high that you're not even close to any minimium, so we can't blame the NMA.
Try using the L-BFGS minimizer (not parallelized yet, unfortunately), preferrably with switch interactions
and neighborlist that is 2-3 angstrom longer than the cutoff. In double precision you should be able to
achieve forces in the order of 1e-7.
#5 Updated by David van der Spoel about 13 years ago
It seems that this problem is due to the minimization. Your em.mdp tells us that you only did 500 steps of minimization which may not be enough. Indeed minimization can go on for more than 10000 steps before converging. It may be that 3.2.1 was better at minimizing that is another problem in my opinion. And as a matter of fact this does not seem to be the case.