How to Write the Methods Section for a Molecular Dynamics Study (GROMACS Example)

A molecular dynamics methods section is reproducible when a stranger could rerun your simulation from it alone. Report the software and version, the starting structure and how you prepared it, the force field and water model, the box and ions, energy minimization, equilibration, and the production run: integrator, timestep, thermostat, barostat, constraints, electrostatics, cutoffs, length, and number of replicas.

Most rejected methods sections do not fail because the simulation was wrong. They fail because the write-up leaves out the numbers that decide the result. “The protein was simulated for 100 ns in GROMACS” tells a reviewer almost nothing they can check or reproduce. This guide from the StemSkills Lab team (10+ years in sequence and structural bioinformatics, drug discovery and design, and multiscale molecular modeling) lists exactly what a GROMACS methods section must report, in the order you ran it, and ends with a worked example paragraph you can adapt. It sits inside our pillar guide on learning molecular dynamics with GROMACS.

What must a molecular dynamics methods section actually contain?

It must contain every choice that would change the trajectory if someone set it differently. A molecular dynamics run is a chain of decisions, and each one is a parameter a reader needs to reproduce your work or judge whether you did it correctly. The test is simple: hand your methods section to a labmate who has never seen your project. If they cannot rebuild the same system and launch the same run without asking you a question, the section is incomplete.

The single most useful habit is to write the methods as you go, not at the end. Every value you type into an .mdp file or a gmx command is a line your methods section will need. Keep those files, because they are the primary record of what you did.

Here is the difference between a line that gets flagged and a line that survives review.

The rest of this guide walks each block in the order you ran it.

How do you document system preparation and the force field?

Start with the structure and be specific about what you changed. State the source (a PDB accession code, or a homology or AlphaFold model with its details), then every edit you made before simulating: chains kept or removed, ligands and crystallographic waters kept or deleted, missing residues or side chains rebuilt and with which tool, and the protonation states you assigned. Protonation is the step students skip most often, and it changes hydrogen-bond networks and salt bridges. If you set histidine tautomers or assigned states with a tool such as PROPKA or H++, name it and give the pH.

Next, the force field and water model, reported together because they are validated together. GROMACS asks you to pick both when you build the topology with gmx pdb2gmx, so both belong in the methods. Common pairings are CHARMM36m with CHARMM-modified TIP3P, AMBER ff99SB-ILDN with TIP3P, and a GROMOS force field with SPC. Do not write “a standard force field”. Name the exact version and cite its paper, because two versions of the same family can give measurably different behaviour.

Then the box and the solvent. Report the box shape (a rhombic dodecahedron is common because it holds fewer water molecules than a cube for the same minimum image distance), the minimum distance you kept between the solute and the box edge, the water model used to fill it, and how you neutralised charge. If you added salt to a physiological concentration, give the number, for example 0.15 M NaCl, and say which ions.

How do you report equilibration and the production run?

This is where the extractable parameters live, and where reviewers look first. Report energy minimization as an algorithm plus a stopping criterion: steepest descent until the maximum force drops below your tolerance (the GROMACS mdp options manual describes the steep integrator as “A steepest descent algorithm for energy minimization”). Then report equilibration in its two standard phases and name the coupling schemes.

For temperature coupling, do not just write “equilibrated at 300 K”. Name the thermostat and its time constant. The velocity-rescale thermostat is a frequent choice for the production ensemble because, as the same manual notes, it is “similar to Berendsen coupling” but adds a stochastic term so that “the stochastic term ensures that a proper canonical ensemble is generated”. A pure Berendsen thermostat is fine for early equilibration but does not produce a correct canonical ensemble, which is exactly the kind of detail a reviewer checks.

For pressure coupling in the NPT phase and production, name the barostat (Parrinello-Rahman and C-rescale are common in current GROMACS), the reference pressure, and the coupling constant. State whether you applied position restraints during equilibration and to which atoms, because that is what lets the solvent relax around a fixed solute before you release it.

The production run needs the tightest reporting. List the ensemble (usually NPT), the integrator (the manual describes md as “A leap-frog algorithm for integrating Newton’s equations of motion”), the timestep, the constraint algorithm, the electrostatics method, the cutoffs, the total simulated time, and how many independent replicas you ran. The default GROMACS timestep is 0.001 ps, that is 1 fs, but constraining bonds to hydrogen atoms with LINCS lets most protein simulations use 2 fs safely, which halves the cost for the same simulated time. State the constraint choice (h-bonds or all-bonds) because it sets which timestep is valid. For long-range electrostatics, report Particle Mesh Ewald, described in the manual as “Fast smooth Particle-Mesh Ewald (SPME) electrostatics”, and give the van der Waals and Coulomb cutoffs. Use the cutoff scheme your force field specifies (CHARMM36, for instance, expects a force-switched van der Waals cutoff), and report the exact values rather than assuming the defaults were right for your force field.

Two things separate a strong methods section from an average one. First, replicas: one trajectory is an anecdote, so state how many independent runs you did and how you seeded them (usually different initial velocities). Second, honesty about analysis windows: say how much of the trajectory you discarded as equilibration before computing averages, and give that number.

What does a finished GROMACS methods paragraph look like?

Here is a model paragraph built from the blocks above. The values are an example, so replace every one with the setting you actually used. Read it as a template for structure and level of detail, not as a set of numbers to copy.

“All molecular dynamics simulations were performed in GROMACS 2023.3. The starting coordinates were taken from PDB entry 1AKI. Crystallographic waters and ions were removed, missing side-chain atoms were rebuilt, and protonation states were assigned at pH 7.0, with all histidines modelled as the epsilon tautomer. The protein was described with the CHARMM36m force field and solvated with CHARMM-modified TIP3P water in a rhombic dodecahedron box with a minimum solute-to-edge distance of 1.0 nm. The system was neutralised and brought to 0.15 M with Na+ and Cl- ions. Energy was minimised by steepest descent until the maximum force fell below 1000 kJ mol-1 nm-1. The system was then equilibrated for 100 ps in the NVT ensemble followed by 100 ps in the NPT ensemble, with position restraints on all protein heavy atoms. Temperature was maintained at 300 K with the velocity-rescale thermostat (tau_t = 0.1 ps) and pressure at 1 bar with the Parrinello-Rahman barostat (tau_p = 2.0 ps). Production simulations used a leap-frog integrator with a 2 fs timestep; bonds to hydrogen were constrained with LINCS. Long-range electrostatics were treated with Particle Mesh Ewald, and van der Waals and Coulomb interactions used the cutoffs recommended for CHARMM36. Three independent 100 ns replicas were run from different initial velocities. The first 10 ns of each replica was discarded before analysis. RMSD and RMSF were computed with gmx rms and gmx rmsf on backbone atoms, and hydrogen bonds with gmx hbond.”

Notice what the paragraph never says: it never says “standard settings” or “default parameters”. Every choice is named and, where a choice has a paper behind it, that paper is cited in the real version.

How should you cite GROMACS and your force field?

Cite three things: the software, the force field, and any tutorial or protocol you followed closely. For the software, cite the GROMACS reference papers listed on the official GROMACS website and state the exact version you ran, because features and defaults change between releases. The widely used software reference is Abraham et al., “GROMACS: High performance molecular simulations through multi-level parallelism from laptops to supercomputers”, SoftwareX (2015), doi:10.1016/j.softx.2015.06.001.

For the force field, cite its own publication, not the GROMACS paper. The MD engine and the force field are separate scientific objects, and a reviewer will notice if you conflate them. If you learned the workflow from a published tutorial, cite it: a common one for beginners is Justin A. Lemkul, “From Proteins to Perturbed Hamiltonians: A Suite of Tutorials for the GROMACS-2018 Molecular Simulation Package”, Living Journal of Computational Molecular Science (2018), 1(1), 5068, doi:10.33011/livecoms.1.1.5068. Citing the protocol you followed is a mark of good practice, not a weakness.

Troubleshooting: the reporting gaps reviewers flag most

• “Which water model?” You named the force field but not the water. Report both together, because they were parameterised together.
• “How were protonation states assigned?” You wrote “the structure was prepared”. Name the pH, the tool, and any manually set histidine or cysteine states.
• “Was this a single run?” You reported one trajectory. State the number of independent replicas and how they were seeded; if it was genuinely one run, say so and justify it.
• “What thermostat and barostat?” You wrote “300 K and 1 bar”. Give the coupling method and time constant for each, and say which phase used which.
• “What timestep and constraints?” These belong together: a 2 fs timestep is only valid with bonds to hydrogen constrained, so report both or neither makes sense.
• “How much did you discard as equilibration?” You averaged over “the trajectory”. State the analysis window and the discarded portion explicitly.
• “What are your cutoffs?” You relied on defaults. Report the actual van der Waals and Coulomb cutoffs, and match them to what your force field expects.

Frequently asked questions

How long should a molecular dynamics methods section be?

Long enough to reproduce the run and no longer. For a single-system study that is usually one to two paragraphs. Depth comes from naming every parameter, not from prose. If a value would change the trajectory, it belongs in the section.

Do I need to report the random seed?

Report how replicas were made independent, which usually means different initial velocities generated from different seeds. Exact bitwise reproducibility is not guaranteed across different hardware or processor counts anyway, so the useful thing to report is the number of replicas and how they differed, not a single seed presented as if it fixed everything.

Should the force field citation go in the methods or the references?

Both. Name the force field and its version in the methods sentence, with the citation marker there, and let the full reference sit in your bibliography. Citing the GROMACS paper does not cover the force field; they are separate works.

Can I just write “simulated with default parameters”?

No. Defaults change between GROMACS versions, and many force fields need non-default cutoffs to be used correctly. “Default parameters” is unverifiable, which is the opposite of what a methods section is for. Write the values out.

What is the fastest way to keep my methods accurate?

Keep your .mdp files and the exact gmx commands you ran, and write the methods section from them rather than from memory. Those files are the ground truth; anything you reconstruct later is a guess.

A methods section is not paperwork bolted on at the end. It is the part of your study that lets the work be trusted, reused, and cited, and it is often the first thing a reviewer reads closely. If you want to see where documenting a run sits alongside preparing structures, running docking, and analysing trajectories, our computational biology skills roadmap lays out the order we teach these skills in and why.