Issue with GROMACS Simulation : MDP file parameters and NVT equilibration

User discussions
, ,
1

Dear GROMACS Community,

I am encountering an issue with my GROMACS simulation, specifically when analyzing the output from the gmx energy command. I am running an NVT equilibration simulation using GROMACS version 2022, but the energy output suggests that the simulation might not have completed successfully or that there’s a problem with the energy file. Below, I’ve provided a detailed description of the problem, simulation settings, and troubleshooting steps I’ve taken so far. I’d greatly appreciate your insights!

This is a simulation of Docked structure of two proteins**

Simulation Details

  • GROMACS Version: 2022

  • Simulation Type: NVT equilibration

  • Command Used for Simulation:
    gmx mdrun -deffnm nvt

  • MDP File Settings:
    itle = FAS Charmm36 NVT equilibration
    define = -DPOSRES ; position restrain the protein
    ; Run parameters
    integrator = sd ; Langevin dynamics integrator
    nsteps = 50000 ; 2 * 50000 = 100 ps
    dt = 0.002 ; 2 fs
    ; Output control
    nstxout = 500 ; save coordinates every 1.0 ps
    nstvout = 500 ; save velocities every 1.0 ps
    nstenergy = 500 ; save energies every 1.0 ps
    nstlog = 500 ; update log file every 1.0 ps
    ; Bond parameters
    continuation = no ; first dynamics run
    constraint_algorithm = shake ; shake algorithm
    constraints = h-bonds ; bonds involving H are constrained
    ; Nonbonded settings
    cutoff-scheme = Verlet ; Buffered neighbor searching
    ns_type = grid ; search neighboring grid cells
    nstlist = 10 ; 20 fs, largely irrelevant with Verlet
    rcoulomb = 1.4 ; short-range electrostatic cutoff (in nm)
    rvdw = 1.4 ; short-range van der Waals cutoff (in nm)
    rvdw-switch = 1.2 ; Start switching at 12 Å
    vdw-modifier = Force-switch ; Enable switching function
    DispCorr = EnerPres ; account for cut-off vdW scheme
    ; Electrostatics
    coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics
    pme_order = 4 ; cubic interpolation
    fourierspacing = 0.16 ; grid spacing for FFT
    ; Temperature coupling is on
    tcoupl = no ; Temperature coupling is handled by the sd integrator
    tc-grps = Protein Non-Protein ; two coupling groups - more accurate
    tau_t = 1.0 1.0 ; time constant, in ps
    ref_t = 300 300 ; reference temperature, one for each group, in K
    ; Pressure coupling is off
    pcoupl = no ; no pressure coupling in NVT
    ; Periodic boundary conditions
    pbc = xyz ; 3-D PBC
    ; Velocity generation
    gen_vel = yes ; assign velocities from Maxwell distribution
    gen_temp = 300 ; temperature for Maxwell distribution
    gen_seed = -1 ; generate a random seed

Issue Description

When I run the following command to extract the temperature:
gmx energy -f nvt.edr -o temperature.xvg
and select option 16 (Temperature), the output is:
Temperature 301.807 – 0 0 (K)
This shows that only one energy frame—at time 0.000 ps—was read, instead of the expected 101 frames over 100 ps. This suggests the simulation either didn’t progress beyond the initial step or the nvt.edr file wasn’t written correctly.

Troubleshooting Steps Taken

  1. Checked Simulation Completion:
    The nvt.log file ends at step 0, indicating the simulation didn’t proceed beyond the initial step.

  2. Inspected Energy File:
    The gmx energy output confirms that nvt.edr contains only one frame, consistent with the log.

  3. Reviewed MDP Settings:
    Verified that nstenergy = 500 should save energies every 1 ps, which seems correct for a 100 ps run.

  4. Attempted to Visualize Trajectory:
    I tried opening nvt.trr in Chimera and PyMOL, but got an error:

I’d like your help with the following:

  1. What could cause the simulation to stop at step 0?
  2. Is the extremely negative pressure (-5,101.99 bar) a sign of an unstable initial configuration, and how can I fix this?
  3. Are there specific settings in the MDP file I should check or adjust to ensure the simulation runs properly?

Thank you so much for your time and assistance—I’m eager to get this simulation running smoothly!

2

hi!

try with this option: tcoupl=berendsen (for NVT equilibration) or tcoupl=nose-hoover (for NVT production).

3

As I am using sd integrator, is it necessary to add a separate thermostat?

4

yes. The sd integrator is the way it is solving the equations of motion. The tcoupl is the way the ensemble (in this case NVT) is coupling the temperature.

5

Hi @pritam_biotech

Are you sure the simulation produced any output? How large are the .trr and .edr files? It might be that the simulation never started and you only printed the first frame. Given that the log file also reports only step 0, I am pretty sure this was the case.
Also regarding the .trr I suggest you to open it with VMD or to convert the .trr in an .xtc which should be handled correctly by Pymol.

The mdp seems okay to me. The only thing is that you are really printing a ton, and you are printing positions/velocities/forces which will end up in a huge chunky file. Do you really need all this data?

The manual states that tcoupl entry is ignored when using the sd integrator, so the thermostat entry will not be considered anyway.

1 Like