2. Hands-On Session 2: Perform molecular dynamics (MD) simulations using AMBER

Julie Puyo-Fourtine¹, Erika McCarthy¹, Lauren Lerew¹, Solen Ekesan¹, and Darrin M. York¹

¹Laboratory for Biomolecular Simulation Research, Institute for Quantitative Biomedicine and Department of Chemistry and Chemical Biology, Rutgers University, Piscataway, NJ 08854, USA

2.1. Learning Objectives

Understand the basic workflow of molecular dynamics (MD) simulations using AMBER, from equilibration to production.
Perform stepwise equilibration, apply and gradually release positional restraints during equilibration.
Read and modify AMBER input files (.mdin) for common MD stages (minimization, heating, equilibration, production).
Use cpptraj to process MD outputs and prepare structures for visualization.
Visualize MD results with VMD and compare structures before and after equilibration.
Understand the special considerations made to stably equilibrating a nucleic acid system.

2.2. Activities

        flowchart LR

%% =========================
%% Equilibration
%% =========================

A2["Input files<br/>MTR1_counterIon_water_bulk.parm7"]

A1["Input files<br/>MTR1_counterIon_water_bulk.rst7<br/>mdins/*.mdin<br/>restrnt.disang"]

P1{{"<b>[§2.3.4.2]</b><br/>Run equilibration<br/>sander / pmemd.cuda"}}

O1_1["<b>[from §2.3.4.2]</b><br/>Equilibration outputs<br/>mdout<br/>equil.nc<br/>"]

O1_2["<b>[from §2.3.4.2]</b><br/>Final restart<br/>09_Equil_solute_Rwt2-run.rst7"]

A1 --> P1
A2 --> P1
P1 --> O1_1
P1 --> O1_2


%% =========================
%% Wrap equilibrated restart and trajectory
%% =========================

P2{{"<b>[§2.3.4.3]</b><br/>Wrap equilibration<br/>cpptraj"}}

O2["<b>[from §2.3.4.3]</b><br/>Wrapped equilibration files<br/>centered.09_Equil_solute_Rwt2-run.rst7<br/>centered_equil.nc"]

O1_1 --> P2
O1_2 --> P2
P2 --> O2


%% =========================
%% Inspect equilibration
%% =========================

P3{{"<b>[§2.3.4.3]</b><br/>Inspect equilibration<br/>VMD"}}

O3["Equilibration inspection<br/>RNA RMSD<br/>ion positions"]

O2 --> P3
P3 --> O3


%% =========================
%% Prepare and Run Production
%% =========================

P4{{"<b>[§2.3.5.2]</b><br/>Apply HMR<br/>cpptraj"}}
O4["<b>[from §2.3.5.2]</b><br/>HMR topology<br/>MTR1_counterIon_water_bulk-run.parm7"]

A2 --> P4
P4 --> O4


P5{{"<b>[§2.3.5.3]</b><br/>Run production MD<br/>pmemd.cuda"}}
O5_1["<b>[from §2.3.5.3]</b><br/>Production logs<br/>MTR1_counterIon_water_bulk-run-01.mdout<br/>MTR1_counterIon_water_bulk-run-01.mdinfo"]

O5_2["<b>[from §2.3.5.3]</b><br/>Production trajectory/restart<br/>MTR1_counterIon_water_bulk-run-01.nc<br/>MTR1_counterIon_water_bulk-run-01.rst7"]


O1_2 --> P5
O4 --> P5
P5 --> O5_1
P5 --> O5_2

P6{{"<b>[§2.3.5.4]</b><br/>Wrap production<br/>cpptraj"}}
O6["<b>[from §2.3.5.4]</b><br/>Wrapped production files<br/>centered.MTR1_counterIon_water_bulk-run-01.nc<br/>centered.MTR1_counterIon_water_bulk-run-01.rst7"]

O5_2 --> P6
P6 --> O6

O7["Production inspection<br/>RNA RMSD<br/>ion positions"]

P7{{"<b>[§2.3.5.4]</b><br/>Inspect production<br/>VMD"}}
O6 --> P7
P7 --> O7

%% =========================
%% Styling
%% =========================

classDef file fill:#fff7e6,stroke:#d98c00,stroke-width:1.5px,color:#111;
classDef program fill:#e8f1ff,stroke:#1f77b4,stroke-width:1.8px,color:#111;
classDef result fill:#eaf7ea,stroke:#2ca02c,stroke-width:1.5px,color:#111;

class A1,A2 file;
class P1,P2,P3,P4,P5,P6,P7 program;
class O1_1,O1_2,O2,O3,O4,O5_1,O5_2,O6,O7 result;

2.2.1. Starting structure and input files 

This activity starts from pre-generated AMBER topology and coordinate files that define the initial structure of the fully solvated MTR1 ribozyme system [1][2] with counterions, described using the OL3 RNA force field [3] and the TIP4P-Ew water model [4]. These files serve as the starting point for all subsequent equilibration and production molecular dynamics simulations performed in this activity.

Starting from these files, the system will be equilibrated stepwise and then used to run production molecular dynamics simulations, generating new restart and trajectory files at each stage.

The specific starting files used in this activity are MTR1_counterIon_water_bulk.parm7 and MTR1_counterIon_water_bulk.rst7.

These files are used to demonstrate the complete equilibration and production molecular dynamics workflow.

2.2.2. Accessing the activity files 

All input files and reference output files required for this activity are available on the DIPC cluster and as a downloadable archive.

DataDir="/data/amber26/TutorialData/HandsOn2_MD"
TutorialDir="/scratch/$USER/TutorialData/HandsOn2_MD"
mkdir -p $TutorialDir
cp -r $DataDir/* $TutorialDir
cd $TutorialDir
cp -r input working_dir
cd working_dir

Option B — Downloading the archive

Hands-On Session 2 files archive (tar.gz)

tar -xzf HandsOn2_MD.tar.gz
cd HandsOn2_MD

2.2.3. Equilibration and Production Workflow 

        flowchart LR

A["Initial system<br/>Inputs: parm7, rst7 (ions and solvent)"]

subgraph EQ["Equilibration"]
   direction LR
   B["Solvent equilibration<br/>Restrained solute"]
   C["Heating to 300 K"]
   D["Gradual solute release"]
   E["Wrap and inspect equilibration<br/>Centered equilibration structure"]
   B --> C --> D --> E
end

subgraph PR["Production"]
   direction LR
   F["Production MD<br/>Unrestrained system"]
   G["Wrap and inspect production<br/>Centered production trajectory"]
   F --> G
end

subgraph AN["Analysis"]
   direction LR
   H["Final inspection and analysis<br/>RMSD, visualization"]
end

A --> B
E --> F
G --> H

2.2.5. Running Production Simulations 

In this section, you will run a very short 0.5 ns production molecular dynamics simulation without positional restraints on the RNA. In previous studies of this RNA system, production simulations converged on the order of tens of nanoseconds. Although longer time scales are desirable, such simulations cannot be performed within the time constraints of this workshop.

2.2.5.1. Preparing the production input files 

Create a new working directory for production and copy all necessary files:

cd ../Prod
cp ../Equil/MTR_waterreps.txt ./

The production MD input file 11_Production_0.5ns.mdin performs 0.5 ns of unrestrained dynamics in the isothermal–isobaric ensemble (NPT) and is intended to provide a brief introduction to running production simulations in AMBER.

Production MD will be run with a 4 fs timestep. To do this stably, hydrogen mass repartitioning (HMR) is applied. The cpptraj input file hmr.in used to apply hydrogen mass repartitioning is shown below:

parm MTR1_counterIon_water_bulk.parm7
hmassrepartition
parmwrite out MTR1_counterIon_water_bulk-run.parm7
go
quit

Run cpptraj to generate the HMR-modified topology file:

cpptraj -i hmr.in

This command generates MTR1_counterIon_water_bulk-run.parm7, which will be used for the production simulation.

Create an empty restraint file in the production directory:

touch restrnt.disang

2.2.5.2. The production SLURM script 

The following SLURM script prod.slurm is used to run the production simulation. It takes as input the prefix of the AMBER topology and restart files and appends -01 to generate the output file names.

#!/bin/bash
#SBATCH --qos=regular
#SBATCH --time=00:20:00
#SBATCH --job-name=prod
#SBATCH --reservation=amber26
#SBATCH --nodes=1
#SBATCH --gpus=1
#SBATCH --ntasks-per-node=1
#SBATCH --mem=10G
#SBATCH --output=prod-%j.out
#SBATCH --error=prod-%j.err

module load Amber/25.0-foss-2022a-AmberTools-25-CUDA-12.1.1-DeePMD

Base=${1}
Parm=${1}.parm7

mdin=mdins/11_Production_0.5ns.mdin
newname=${Base}-01

${AMBERHOME}/bin/pmemd.cuda_SPFP -O -p ${Parm} -c ${Base}.rst7 -i ${mdin} -o ${newname}.mdout -inf ${newname}.mdinfo -x ${newname}.nc -r ${newname}.rst7 -ref ${Base}.rst7

The core AMBER command is pmemd.cuda_SPFP, which runs GPU-accelerated molecular dynamics. Its key flags are:

-O: overwrite existing output files
-i: (input) control data for the MD run
-o: (output) user readable state info and diagnostics
-inf: (output) latest mdout-format energy info
-p: (input) molecular topology, force field, periodic box type, atom and residue names
-c: (input) initial coordinates and (optionally) velocities and periodic box size
-ref: (input, optional) reference coords for position restraints
-x: (output) coordinate sets saved over trajectory
-r: (output) final coordinates, velocity, and box dimensions — used for restarting the run

2.2.5.3. Running the production simulation 

Submit the production job:

sbatch prod.slurm MTR1_counterIon_water_bulk-run

The job should take approximately 10 minutes to complete. Output files will be written with the prefix MTR1_counterIon_water_bulk-run-01.

Note

If you do not want to wait for the production job to finish, reference output files are available. Copy them to a subdirectory:

cp -r ../../output/Prod ./output

2.2.5.4. Wrapping and visualization of the production trajectory 

Wrap the production trajectory:

cpptraj -i wrap_RNA_prod.in

Open the production trajectory in VMD locally. First, download the necessary files:

rsync -avz "username@hyperion.sw.ehu.es:/scratch/amber26-XX/TutorialData/HandsOn2_MD/working_dir/Prod/{centered.MTR1_counterIon_water_bulk-run-01.nc,MTR1_counterIon_water_bulk-run.parm7,MTR_waterreps.txt}" .

Locally: Open the trajectory in VMD:

vmd -e MTR_waterreps.txt \
    centered.MTR1_counterIon_water_bulk-run-01.nc \
    MTR1_counterIon_water_bulk-run.parm7

Inspect the dynamics of the different components. Below is a movie demonstration:

Warning

This was a very short production simulation. Realistically, it would be desirable to perform multiple independent long production simulations. In previous studies of MTR1, five independent production simulations of 150 ns each were performed. Analysis of the root mean square deviation (RMSD) indicated that the first 50 ns corresponded to equilibration, while the remaining 100 ns of each simulation were used for production analysis, resulting in a total of 500 ns of production data.

The figure below shows the RMSD of the RNA heavy atoms with respect to the first trajectory frame over four independent 150 ns production simulations.

RMSD relative to the first trajectory frame — **Figure 3.** Root mean square deviation of the positions of heavy atoms with respect to the first trajectory frame over five independent 150 ns production simulations.

Hint

You can compute the RMSD of your own trajectory using the following cpptraj input file. The rms command calculates the RMSD of all RNA heavy atoms (hydrogens excluded, !@H=) with respect to the first frame (first). The mass keyword enables mass-weighting:

parm MTR1_counterIon_water_bulk-run.parm7
trajin centered.MTR1_counterIon_water_bulk-run-01.nc
rms ToFirst :1-69&!@H= first out rmsd1.dat mass
run
quit

cpptraj -i rmsd.in

2.2.6. (Optional) Additional Exercise: Performing a brief energy minimization 

For this optional exercise, change directories back to the equilibration directory:

cd ../Equil

The following exercise should only be performed if you have time left over or are waiting for other simulations to run. In this exercise, you will run the 000_Min step of the equilibration procedure independently using 000_Min.mdin.

You have already been introduced to the equilibration script eqFromStratch.slurm. In this optional exercise, you will modify this script to run only the initial minimization step on the CPU.

Copy the equilibration script and create a new script for the minimization only run:

cp eqFromStratch.slurm eqFromStratch_000.slurm

Open eqFromStratch_000.slurm in a text editor (e.g. vim) and modify the SLURM header so that the job runs on the CPU rather than on a GPU:

#!/bin/bash
#SBATCH --qos=regular
#SBATCH --time=00:30:00
#SBATCH --job-name=min
#SBATCH --reservation=amber26
#SBATCH --nodes=1
#SBATCH --ntasks-per-node=1
#SBATCH --mem=10G
#SBATCH --output=equil-%j.out
#SBATCH --error=equil-%j.err

Next, edit the script so that only the 000_Min step is executed. Delete all other steps and retain only the following definition:

steps="\
000_Min \
"

In this case, only sander will be invoked, so no GPU resources are required. The job is expected to finish in approximately 10 minutes.

The Equil directory should then contain the following:

ls
MTR1_counterIon_water_bulk.parm7
MTR1_counterIon_water_bulk.rst7
eqFromStratch.slurm
eqFromStratch_000.slurm
mdins
restrnt.disang
wrap_RNA_equil.in

Submit the job:

sbatch eqFromStratch_000.slurm MTR1_counterIon_water_bulk

When the job is complete, monitor how the total energy (EAMBER) decreases over the course of the minimization:

grep -r EAMBER 000_Min-run.mdout | awk '{print $3}' | head -n10 > ene.dat

Plot the energy using xmgrace. Load the module first:

module load Grace
xmgrace ene.dat

Total energy during minimization — **Figure 4.** Total energy output every 50 optimization steps during the `000_Min` minimization.

You will observe that the total energy (in kcal/mol) decreases rapidly at the beginning of the minimization and then begins to plateau. In this run, the energy was printed every 50 minimization steps.

Open 000_Min-run.mdout in a text editor and inspect the individual energy components:

  NSTEP       ENERGY          RMS            GMAX         NAME    NUMBER
     1       1.7197E+06     6.7827E+04     2.2501E+07     N3       1730

BOND    =     5038.2053  ANGLE   =      397.2193  DIHED      =     1511.0059
VDWAALS =  1926187.1298  EEL     =  -206390.2786  HBOND      =        0.0000
1-4 VDW =      697.4115  1-4 EEL =    -7778.5268  RESTRAINT  =        0.0000


  NSTEP       ENERGY          RMS            GMAX         NAME    NUMBER
    50      -2.3201E+05     4.6787E+00     6.6797E+02     C6       1723

BOND    =    22954.6607  ANGLE   =      746.4316  DIHED      =     1649.0907
VDWAALS =    27936.2925  EEL     =  -278521.9124  HBOND      =        0.0000
1-4 VDW =      650.9327  1-4 EEL =    -7954.8135  RESTRAINT  =      533.4714
EAMBER  =  -232539.3177

You will notice that the van der Waals energy decreases rapidly at the beginning of the minimization as bad contacts are relieved. This illustrates the importance of a gradual and well-controlled equilibration procedure.

2.3. Relevant literature

References cited in this tutorial are listed below.