I’m currently using GROMACS to perform normal mode analysis (NMA) on a protein system, and I’m encountering a few issues I’d appreciate help with.
Setup Summary:
I placed the protein in a simulation box and added a 5 Å layer of water around it — the box is not fully solvated.
I neutralized the system by adding positive ions.
I performed energy minimization using the steepest descent algorithm in single precision.
Following it, I did NVT, NPT equilibration and ran 5 ns md simulation.
For more accurate minimization, I performed energy minimization using steepest descent in double precision.
I followed it with conjugate gradient minimization in double precision with em step size=0.001 and emtol =0.001
My target emtol was 0.001, but the minimization stopped early, reporting:
Energy minimization has stopped, but the forces have not converged to the
requested precision Fmax < 0.001 (which may not be possible for your system).
It stopped because the algorithm tried to make a new step whose size was too
small, or there was no change in the energy since last step.
Either way, we regard the minimization as converged to within the available
machine precision, given your starting configuration and EM parameters.
After this, I proceeded with the normal mode analysis, but the results are troubling:
Eigenfrequencies are all zero.
Eigenvalues are negative.
The output message suggests that this could be due to improper minimization.
My Questions:
Is this issue due to the simulation box not being fully solvated? Would full solvation make a difference in NMA accuracy?
Is it valid to perform normal mode analysis in a partially solvated or solvated system at all?
What are other possible reasons behind these unphysical eigenvalues and zero frequencies?
Would a different minimization strategy or a better starting configuration help?
Any insights or suggestions from those who’ve done similar work with GROMACS and NMA would be greatly appreciated!
That should be small enough. How did you transfer the minimized coordinates to grompp for preparing NM? You would need to use a trr or cpt file to retain sufficient precision.
I have used gmx_d grompp -f nm.mdp -c em.gro -p topol.top -o nm.tpr to get the .tpr file and gmx_d mdrun -s nm.tpr -mtx hessian.mtx for creating hessian matrix. I haven’t used checkpoint or trajectory files here. probably this may be the reason for the inaccuracy. I will check to see if this works.
Also, one more question, I have done the energy minimization after the 5 ns md simulation, for preparing the system for the NMA. Could this be the reason for the zero eigen frequencies as doing energy minimization close to my final step, likely means that the protein is at 0K (no thermal agitation).
You were right, everything was because of the rounding off error. Thank you for the insight. While this worked for the small system, I tried to see what will happen for the large one.
After completing a 5 ns MD simulation, I performed a three-step energy minimization protocol to prepare my system for normal mode analysis (NMA). For a smaller system (~4400 atoms with a 5 Å water shell), I used steepest descent minimization with emtol = 1000, followed by conjugate gradient (CG) with emtol = 0.001. While the CG run did not initially reach the force convergence criterion, a subsequent CG minimization did, ultimately yielding physically meaningful eigenfrequencies. I used minimized conformation from the full precision trajectory file to get .tpr file while preparing the system for NMA. The first six frequencies were close to zero (as expected for rigid body motions), and the remaining were non-zero. However, when applying the same protocol to a larger system (~15,000 atoms), I consistently obtain exactly zero values for the first 16 eigenfrequencies, which I suspect is due to the system not reaching the desired emtol during minimization. Given this, I would like advice on how best to proceed. Should I revert to steepest descent minimization for better force relaxation before switching to CG? Or would switching to the L-BFGS minimizer be more appropriate for a system of this size and complexity? I would appreciate any suggestions on how to reliably reach convergence and obtain accurate normal modes for larger systems.