MD Simulation at Different Temperatures in GROMACS

To run MD at different temperatures in GROMACS, build and minimise one system, then give each temperature its own NVT and NPT equilibration and production run. In every .mdp file set ref_t to the target (one value per tc-grps group) and, in NVT only, gen_vel = yes with gen_temp equal to that target. Keep everything else identical.

A thermal stability comparison, typically 300 K against 350 K and 400 K, is one of the most common MSc molecular dynamics projects. It looks simple: change a number and rerun. In practice most failed projects come from three places: a forgotten gen_temp that silently starts every run at 300 K, runs that differ in more than temperature, and plots that are read as if MD temperature were an experimental melting point. This guide from the StemSkills Lab team (10+ years in structural bioinformatics, drug design and molecular modeling) fixes all three. It assumes you can already run the standard protein-in-water workflow from our GROMACS molecular dynamics learning path.

Which .mdp parameters control temperature in GROMACS?

Five parameters matter, and the GROMACS mdp options reference defines each one. If you want a line-by-line walk through the rest of the file, see our GROMACS .mdp file explained guide.

• tcoupl: the thermostat. Use V-rescale, the velocity rescaling thermostat with a stochastic term published by Bussi, Donadio and Parrinello (J. Chem. Phys. 126, 014101, 2007). It samples the correct canonical ensemble. Avoid Berendsen: the GROMACS manual calls it a historical thermostat and says “we strongly recommend not to use it for new production runs.”
• tc-grps: the groups coupled to separate baths. For a protein in water, Protein Non-Protein is the usual choice.
• tau_t: the coupling time constant in ps, one value per group (for example 0.1 0.1).
• ref_t: the reference temperature in K, one value per group. This is the number you change, and it must be changed for every group (350 350, not 350).
• gen_vel, gen_temp, gen_seed: with gen_vel = yes, grompp draws starting velocities from a Maxwell distribution at gen_temp using gen_seed. The default gen_temp is 300 K, so if you only edit ref_t, your 400 K run starts from 300 K velocities and spends its first picoseconds heating up.

Pressure coupling stays at ref_p = 1.0 bar for every temperature. Use pcoupl = C-rescale, which the manual describes as exponential relaxation pressure coupling with a stochastic term “to enforce correct volume fluctuations.”

Why does each temperature need its own NVT and NPT equilibration?

Because the equilibrium state of the water box depends on temperature. An NPT box equilibrated at 300 K has the density and volume of 300 K water. If you start a 400 K production run from it, the first part of production is really equilibration: the box expands, the pressure fluctuates hard, and your RMSD and Rg curves carry that transient.

The size of the change is real. Liquid water at 1 bar has a density of about 997 kg/m³ at 25 °C and about 958 kg/m³ at its 100 °C boiling point (NIST Chemistry WebBook, thermophysical properties of fluids), a drop of roughly 4 percent. Water models in a simulation box show the same kind of trend, so each temperature needs its own NPT phase to find its own box volume. The pattern is:

1. Shared steps (run once): pdb2gmx, editconf, solvate, genion and energy minimisation. One em.gro is the common starting point for every temperature.
2. NVT at T: gen_vel = yes, gen_temp = T, ref_t = T T, continuation = no.
3. NPT at T: gen_vel = no, continuation = yes, ref_t = T T, read the NVT checkpoint.
4. Production at T: same as NPT without position restraints, read the NPT checkpoint.

Our NVT and NPT equilibration guide covers how to check that each phase has converged (temperature and density plateaus from gmx energy). Do that check at every temperature, not just at 300 K.

How do I set up 300 K, 350 K and 400 K runs without editing files by hand?

Write each .mdp once as a template with a placeholder, then generate the three versions with sed. In nvt_template.mdp, the temperature lines look like this:

tcoupl      = V-rescale
tc-grps     = Protein Non-Protein
tau_t       = 0.1     0.1
ref_t       = TEMP    TEMP
gen_vel     = yes
gen_temp    = TEMP
gen_seed    = -1
continuation = no

In npt_template.mdp and md_template.mdp, keep the same ref_t = TEMP TEMP line but set gen_vel = no and continuation = yes. Then loop:

for T in 300 350 400; do
  for stage in nvt npt md; do
    sed "s/TEMP/$T/g" ${stage}_template.mdp > ${stage}_$T.mdp
  done

  gmx grompp -f nvt_$T.mdp -c em.gro -r em.gro -p topol.top -o nvt_$T.tpr
  gmx mdrun -deffnm nvt_$T

  gmx grompp -f npt_$T.mdp -c nvt_$T.gro -r nvt_$T.gro -t nvt_$T.cpt -p topol.top -o npt_$T.tpr
  gmx mdrun -deffnm npt_$T

  gmx grompp -f md_$T.mdp -c npt_$T.gro -t npt_$T.cpt -p topol.top -o md_$T.tpr
  gmx mdrun -deffnm md_$T
done

Before launching, check the generated files with grep -E "ref_t|gen_temp" nvt_*.mdp. It takes five seconds and catches the most common mistake in this whole project. On a shared cluster, submit each temperature as its own job instead of running the loop serially.

How do I keep runs at different temperatures comparable?

Temperature must be the only variable. A reviewer or examiner will check these, so fix them before you run anything:

Run at least three replicas per temperature, each with a different seed (gen_seed = -1 picks a pseudo-random seed each time). One trajectory per temperature cannot separate a temperature effect from chance. How long each replica should be depends on what you are measuring; our guide on how long to run an MD simulation explains how to decide and justify it.

Which analyses show thermal stability, and how do I read them?

First remove periodic boundary artefacts, then run the same four analyses on every trajectory. Use one common reference for RMSD (the minimised structure, em.tpr), so all temperatures are measured against the same starting point.

for T in 300 350 400; do
  gmx trjconv -s md_$T.tpr -f md_$T.xtc -o md_${T}_noPBC.xtc -pbc mol -center
  gmx rms    -s em.tpr   -f md_${T}_noPBC.xtc -o rmsd_$T.xvg -tu ns
  gmx gyrate -s md_$T.tpr -f md_${T}_noPBC.xtc -o gyrate_$T.xvg -tu ns
  gmx sasa   -s md_$T.tpr -f md_${T}_noPBC.xtc -o sasa_$T.xvg -tu ns
  gmx dssp   -s md_$T.tpr -f md_${T}_noPBC.xtc -o dssp_$T.dat -num dssp_num_$T.xvg -tu ns
done

Each tool asks for groups interactively: choose Protein to centre and output in trjconv, Backbone for RMSD fitting and calculation, and Protein for gyrate and sasa. gmx dssp exists in GROMACS 2023 and later; older versions used gmx do_dssp with an external DSSP binary. Step-by-step reading of each plot is in our guides on RMSD and RMSF, radius of gyration and SASA, and secondary structure with gmx dssp.

Plot the three temperatures on the same axes (with replica mean and spread), and report averages over the same equilibrated window of each run. A clean result reads like: “Rg and DSSP are unchanged between 300 and 350 K; at 400 K two of three replicas lose the C-terminal helix after about 60 ns.” That sentence states the effect, the threshold and the replica support.

Why did my protein not unfold at 400 K?

Often because it should not, on your timescale. Two reasons apply to most student systems.

First, unfolding is slow. Even at high temperature, a compact globular protein can stay folded for the whole of a 100 ns run. That is a result, not a failure: report it as “no loss of fold within X ns at 400 K”.

Second, covalent cross-links hold the fold together. The classic tutorial protein, hen egg-white lysozyme (PDB 1AKI, used in Justin Lemkul’s GROMACS lysozyme tutorial, published as Lemkul, LiveCoMS 2018), has four disulfide bonds: Cys6 to Cys127, Cys30 to Cys115, Cys64 to Cys80 and Cys76 to Cys94. pdb2gmx builds them as bonds in the topology, and classical MD cannot break covalent bonds at any temperature. Expect local loosening, not a full unfolding. When you choose a protein for your own project, check for disulfides first and mention them when you interpret the results.

Tap to explore in 3D

The folded proteinThe four disulfide bonds

Can I compare MD temperatures to experimental melting temperatures?

Not directly. Force fields and water models are parameterised mainly near room temperature, and simulations are far shorter than experiments, so a protein that unfolds at 400 K in 100 ns has not been shown to melt at 400 K. Use high temperature as an accelerant that ranks stability. The defensible claim compares variants or conditions at the same temperature, for example a mutant that loses its fold at 400 K while the wild type holds. That is the design in our wild-type vs mutant MD simulation guide, and temperature combines well with it. Keep both axes in the methods section so the examiner sees exactly what was varied.

What goes wrong in high-temperature runs, and how do I fix it?

• grompp error “Invalid T coupling input: 2 groups, 1 ref-t values”. You wrote one ref_t value for two tc-grps. Give one value per group: ref_t = 400 400. The same rule applies to tau_t.
• NVT temperature starts low and climbs. gen_temp was left at its 300 K default. Set it equal to ref_t and regenerate the .tpr.
• LINCS warnings or a crash early in the 400 K run. Atoms move faster, so constrained bonds rotate further per step. Confirm NVT and NPT at 400 K both converged first. If it persists, reduce dt from 0.002 to 0.001 ps, and apply the same dt to all temperatures to keep them comparable.
• “Pressure scaling more than 1%” warnings. The box is still far from its equilibrium volume at that temperature. Lengthen NPT at that temperature, or check you started it from the NVT output of the same temperature.
• RMSD jumps of several nanometres. The protein crossed the periodic boundary. Rerun gmx trjconv with -pbc mol -center before any analysis.
• Every temperature looks identical. Check the log files: grep "ref-t" md_*.log prints the parameters mdrun actually used. If they all say 300, the template substitution did not happen.

Once the study is running cleanly, place it in the bigger picture of skills employers and PhD panels look for with our computational biology skills roadmap.

Frequently asked questions

What temperatures should I use for a thermal stability MD study?

A common design is 300 K as the physiological baseline plus two higher temperatures such as 350 K and 400 K. Choose values that bracket the behaviour you want to compare, and justify them in your methods rather than copying them from another paper.

Do I need to re-run energy minimisation for each temperature?

No. Minimisation has no temperature, so one minimised structure serves every run. Only NVT, NPT and production are repeated per temperature.

Should I change pressure when I raise the temperature?

No. Keep ref_p = 1.0 bar so temperature is the only variable. The box volume adjusts itself during the NPT phase at each temperature.

Can I continue my 300 K run at 400 K with a checkpoint?

You can, but it is a heating protocol, not an independent 400 K run, and the early frames are not equilibrated. For a comparison study, start each temperature from the same minimised structure with fresh velocities.

How many replicas do I need per temperature?

Three per temperature, with different gen_seed values, is the minimum most supervisors accept. It lets you report a mean and spread and shows that an unfolding event is reproducible.