How to Run an MD Simulation of DNA or RNA in GROMACS

The GROMACS workflow for DNA or RNA is the same chain of commands you use for a protein, but four steps change. You must install a contributed AMBER nucleic acid force field, neutralise a much larger backbone charge with matched ion parameters, pad the box for an elongated molecule, and replace protein-only analyses such as gmx dssp with base pair and helical parameters.

This guide is part of our molecular dynamics with GROMACS series, written by the StemSkills Lab team, who have spent more than 10 years in structural bioinformatics, drug design and multiscale molecular modeling. If you have just finished docking a ligand to DNA or RNA and want to simulate the complex you built, this is the next step.

What actually changes when you simulate a nucleic acid instead of a protein?

Less than you fear, but the parts that change will stop you cold if you do not know about them in advance. The command chain is unchanged: download the structure, build the topology with pdb2gmx, set the box with editconf, solvate, add ions, minimise, equilibrate, run production, analyse.

Four of those steps need different decisions:

1. Force field. The force fields GROMACS installs by default are not the ones the nucleic acid field actually uses. You have to install one.
2. Ions. Every phosphate carries a formal charge of minus one, so counterions are a structural part of the system rather than a rounding correction.
3. Box. A B-DNA duplex is elongated rather than globular, so padding that is fine for a small protein can let the molecule meet its own periodic image.
4. Analysis. Secondary structure tools are written for amino acids. You need base pair, base step and helical parameters instead.

Which force field should you use for DNA or RNA in GROMACS?

For double-stranded DNA, use an AMBER OL package (OL24 is the current release) or parmbsc1. For RNA, use the chi-OL3 parameters, which the OL packages ship alongside the DNA set. None of these comes with GROMACS, so all of them have to be downloaded and installed into your working directory.

The GROMACS reference manual lists native support for AMBER94, AMBER96, AMBER99, AMBER99SB, AMBER99SB-ILDN, AMBER03, AMBERGS, AMBER14SB and AMBER19SB. All of those are protein parameter sets. On the CHARMM side the manual is explicit that “The CHARMM27 force field has been ported to GROMACS and is officially supported”, and CHARMM27 does cover nucleic acids, but it predates two decades of nucleic acid refinement. The nucleic acid community has moved on, and so should your thesis.

Where the field has moved is documented in the long-running Cheatham group assessment series. In Assessing the Current State of Amber Force Field Modifications for DNA (Journal of Chemical Theory and Computation, 2016, volume 12, pages 4114 to 4127), Galindo-Murillo and co-workers compared parmbsc1 and OL15 on the Drew-Dickerson dodecamer, the twelve base pair B-DNA duplex deposited as PDB 1BNA, using 100 independent simulations of at least 10 microseconds each, concatenated into a one millisecond trajectory per force field and water model combination. The 2023 edition of that assessment, by Love and colleagues, aggregated roughly 7.75 milliseconds of sampling and concluded that the improved parameters in OL21 make it the optimal double-stranded DNA force field of those then available, used with the OPC water model.

Parmbsc1 remains a defensible and very widely used choice. Its primary reference is Parmbsc1: a refined force field for DNA simulations (Nature Methods, 2016, volume 13, pages 55 to 58), whose abstract states that it “has been parameterized from high-level quantum mechanical data and tested for nearly 100 systems (representing a total simulation time of ~140 micro-seconds) covering most of DNA structural space”.

One warning that the OL force field group at Palacky University Olomouc prints on its own front page, and that students regularly miss: “Note that currently different parameter sets are needed for DNA and RNA simulations!” A package that is excellent for your duplex is not automatically the right answer for an RNA hairpin.

A caveat worth planning around: the OPC water model recommended in the 2023 assessment is not in the watermodels.dat list of any of the three GROMACS packages above. Each of them offers TIP3P, TIP4P, TIP4P-Ew, SPC and SPC/E, with TIP3P flagged as recommended. If you want OPC you will have to add it yourself, which is a real piece of extra work rather than a flag. For an MSc project, running OL24 or parmbsc1 with the TIP3P that ships in the package is the sane decision, and you state the water model in your methods.

How do you install a contributed force field so pdb2gmx can see it?

You do not install it into the GROMACS tree. You unpack the .ff directory next to your structure file. The gmx pdb2gmx documentation defines the search path precisely:

“gmx pdb2gmx will search for force fields by looking for a forcefield.itp file in subdirectories <forcefield>.ff of the current working directory and of the GROMACS library directory as inferred from the path of the binary or the GMXLIB environment variable.”

So for OL24:

wget https://fch.upol.cz/ff_ol/amber14sb_OL24.ff.tar.gz
tar -xzf amber14sb_OL24.ff.tar.gz
ls amber14sb_OL24.ff/forcefield.itp

Run gmx pdb2gmx from that directory and the package appears at the top of the interactive menu, above the shipped force fields. You can skip the menu with -ff, which the same page describes as a way “to specify one of the short names in the list on the command line instead”, in which case pdb2gmx “just looks for the corresponding <forcefield>.ff directory”:

gmx pdb2gmx -f dna.pdb -o dna_processed.gro -p topol.top -ff amber14sb_OL24 -water tip3p

The parmbsc1 equivalent lives on the GROMACS user contributions page as amber14sb_parmbsc1.ff.tar.gz. Checked on 1 October 2026, that file downloads correctly, and so does amber14sb_OL15.ff_corrected-Na-cation-params.tar.gz, but the plain amber14sb_OL15.ff.tar.gz link on the same page returns a 404. If you want OL15 rather than OL21 or OL24, use the corrected-cation archive or take it from the OL group site instead of assuming the download is broken on your end. Our wider guide to choosing a force field in GROMACS covers how to justify the choice in writing.

Which residue and atom names does the force field expect?

This is the question that generates most pdb2gmx failures on nucleic acids, and the answer is readable directly from the package. Inside amber14sb_OL24.ff the file dna.r2b maps the names DA, DG, DC and DT, while rna.r2b maps the single letters A, U, C and G. Those are exactly the residue names the RCSB PDB uses, so a structure downloaded unmodified will normally be recognised.

Each line in those tables lists four building blocks, not one: the internal nucleotide, the 5′ terminal form, the 3′ terminal form and the form used for an isolated nucleotide. For adenine in DNA those are DA, DA5, DA3 and DAN. GROMACS selects between them from the position in the chain, so you should not rename terminal residues yourself. Doing so is a common and self-inflicted error.

Atom names are handled the same way. The dna.arn and rna.arn files translate the modern PDB version 3 names into the internal AMBER names, including OP1 to O1P, OP2 to O2P, H2' to H2'1 and HO5' to H5T. That translation is why a current PDB file works without editing, and it is also why hand-editing atom names usually makes things worse.

Why do ions matter more for a nucleic acid than for a protein?

Because the charge is large, systematic and structural. Each phosphate group in the backbone carries a formal charge of minus one, so the net charge of a duplex scales with its length instead of hovering near zero the way a globular protein’s often does. Counterions are not there to tidy up a residual charge; they condense around the backbone and are part of the physics you are trying to reproduce.

Do not calculate the charge by hand. Whether a 5′ terminus carries a phosphate differs between deposited structures, and pdb2gmx reports the total charge of the topology it has just written. Read that number and neutralise it.

The gmx genion documentation describes the tool as a utility that “randomly replaces solvent molecules with monoatomic ions”, and its defaults are easy to trip over: -pname defaults to NA, -nname to CL, -np and -nn to 0, -conc to 0, and -neutral is off. Nothing happens unless you ask. For a nucleic acid you normally want both neutralisation and background salt:

gmx genion -s ions.tpr -o solv_ions.gro -p topol.top -pname NA -nname CL -neutral -conc 0.15

The manual states that -neutral “will add enough ions to neutralize the system” and that “These ions are added on top of those specified with -np/-nn or -conc”, so the two options combine rather than compete. The -conc value you choose is a reported methods parameter, not a default to leave unexamined.

Match the ion parameters to the force field rather than accepting whatever is nearest. The OL21 package documents Joung-Cheatham monovalent ion parameters and Allner magnesium parameters, which matters because the Joung-Cheatham set, published as Determination of Alkali and Halide Monovalent Ion Parameters for Use in Explicitly Solvated Biomolecular Simulations (Journal of Physical Chemistry B, 2008, volume 112, pages 9020 to 9041), was parameterised against specific water models including TIP3P, TIP4P-Ew and SPC/E. Pairing ions from one water model with another is a quiet way to introduce error. If your project involves a magnesium-dependent RNA, check that the package actually contains magnesium parameters before you start.

How should you set the box for a DNA duplex?

Think about the longest dimension, not the average one. A globular protein tumbles inside roughly the same envelope whichever way it turns. An elongated duplex does not, so padding chosen from its width will be too small once it rotates, and the molecule can then interact with its own periodic image.

In gmx editconf the relevant options are -d, documented as the “Distance between the solute and the box” and defaulting to 0, and -bt, the “Box type for -box and -d: triclinic, cubic, dodecahedron, octahedron”, defaulting to triclinic. A rhombic dodecahedron holds the same minimum image distance in fewer water molecules than a cube, which is why it is the usual choice:

gmx editconf -f dna_processed.gro -o dna_newbox.gro -c -d 1.5 -bt dodecahedron

Treat the -d value as something to check rather than copy. Pick a value, build the box, and confirm that the duplex cannot reach its image even when rotated end over end. If it can, increase -d and accept the extra solvent. There is no universal number here, which is precisely why you should verify it for your own system instead of inheriting a figure from a tutorial written for a small protein.

Which analyses replace gmx dssp for DNA and RNA?

Start by accepting that the protein tool does not transfer. The gmx dssp documentation defines it as applying the DSSP algorithm “by detecting specific patterns of hydrogen bonds between amino acid residues” to determine “the secondary structure of a protein”. Our post on secondary structure analysis with gmx dssp therefore does not apply to your duplex.

The nucleic acid equivalents are helical and base pair descriptors, and the standard route from a GROMACS trajectory is do_x3dna, published by Kumar and Grubmuller as do_x3dna: a tool to analyze structural fluctuations of dsDNA or dsRNA from molecular dynamics simulations (Bioinformatics, 2015, volume 31, issue 15, pages 2583 to 2585). It wraps the 3DNA package so it can run frame by frame over a trajectory, taking a .tpr, a .trr or .xtc, and optionally an index file.

What it gives you is the vocabulary a nucleic acid paper is written in:

• Base pair parameters: shear, stretch, stagger, buckle, propeller, opening.
• Base step parameters: shift, slide, rise, tilt, roll, twist.
• Helical parameters: X-displacement, Y-displacement, helical rise, inclination, tip, helical twist.
• Groove geometry: major and minor groove widths, which is the observable that matters if you are studying a groove binder.
• Backbone and sugar torsions for both strands, plus helical radii and the local helical axis.

The standard measures still work, with one change in how you apply them. Terminal base pairs fray: they open and close at the ends of a duplex, and this is expected behaviour rather than a failed simulation. If you fit and measure across the whole molecule, that fraying dominates your numbers. In the Drew-Dickerson dodecamer of PDB 1BNA, for instance, the two strands are numbered 1 to 12 and 13 to 24, so the fraying pairs are residues 1 with 24 and 12 with 13, and the core you want is residues 2 to 11 with 14 to 23. Build an index group that excludes the terminal base pairs and use it as both the fit and the measurement group, then run RMSD and RMSF and radius of gyration, SASA and hydrogen bonds on the core as usual. As with any MD result, report means across independent replicas rather than a single trajectory, and write the parameters up using our MD methods section guide.

Tap to explore in 3D

The duplex as depositedBoth strands as sticksThe core you fit onTerminal base pairs that fray

What about a protein-DNA complex?

This is the easy case, and it is the reason the combined packages exist. amber14sb_OL24.ff, amber14sb_OL21.ff and amber14sb_parmbsc1.ff all pair their nucleic acid parameters with the ff14SB protein force field in a single directory, so pdb2gmx can build a transcription factor and its operator from one PDB file without any manual merging. The OL group additionally notes that OL24 is recommended for protein-DNA work because those interfaces often involve A-form DNA, whose stability older AMBER variants underestimated.

If a small molecule is also present, it still needs separate ligand parameterisation, exactly as in a protein-ligand simulation.

Troubleshooting: real errors and what fixes them

Frequently asked questions

Can I use the same force field for DNA and RNA?

Only if the package explicitly contains both sets, and even then they are different parameters. The OL packages ship OL24 or OL21 for DNA and chi-OL3 with bsc0 for RNA in one directory, which is convenient, but the OL group states plainly that different parameter sets are needed for DNA and RNA. Name the one you used for the molecule you simulated.

Does CHARMM36 work for nucleic acids in GROMACS?

Yes. The GROMACS reference manual states that GROMACS supports the CHARMM force field for proteins, lipids and nucleic acids, and CHARMM36 files in GROMACS format are maintained by the MacKerell lab. The manual also lists the .mdp settings CHARMM36 requires, including constraints = h-bonds, cutoff-scheme = Verlet, vdwtype = cutoff, vdw-modifier = force-switch, rlist = 1.2, rvdw = 1.2, rvdw-switch = 1.0, coulombtype = PME, rcoulomb = 1.2 and DispCorr = no. Our .mdp file guide explains what each of those does.

How long should I run a nucleic acid simulation?

Longer than you would like, and the assessment literature is the honest guide. The Cheatham group studies behind the current recommendations use microsecond-scale sampling per replica and aggregate into the millisecond range to reach convergence on a twelve base pair duplex. An MSc project will not match that, so choose a length you can justify, run several independent replicas rather than one long run, and state the limitation instead of overclaiming.

Do I need a special water model?

TIP3P is flagged as recommended inside the OL24, OL21 and parmbsc1 GROMACS packages, and it is a defensible choice. The 2023 Amber DNA assessment prefers OPC, but OPC is not in those packages’ water model lists, so using it means adding the model yourself. Decide before you start and report what you used.

Can I simulate a G-quadruplex or an aptamer this way?

The setup chain is the same, but both bring extra requirements. G-quadruplexes depend on specific coordinated cations in the central channel, so check that your ion parameters and starting coordinates include them rather than letting genion place ions at random. Aptamers are RNA, so use the RNA parameters. Our MD project ideas post has scoping advice, and the computational biology skills roadmap shows where this sits in a full training path.