Normal Mode Analysis – Minimization Doesn’t Reach Target emtol & Zero Eigenfrequencies

User discussions
1

Hi everyone,

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:

    1. Is this issue due to the simulation box not being fully solvated? Would full solvation make a difference in NMA accuracy?

    2. Is it valid to perform normal mode analysis in a partially solvated or solvated system at all?

    3. What are other possible reasons behind these unphysical eigenvalues and zero frequencies?

    4. 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!

2

What is your max force after the last minimization?

3

Hi. This is the maximum force at the end.

MStep 892008, Epot=-6.790620e+04, Fnorm=6.433e-04, Fmax=9.834e-03 (atom 1345)

4

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.

5

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).

6

Likely everything is caused by moving away from the minimum by rounding the coordinates.