Hands-On Session 2: Introduction to Molecular Dynamics Simulations with Amber

Erika McCarthy¹, Julie Puyo-Fourtine¹, 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

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.

Relevant Literature

The following references provide background and context for the methods and best practices used in this tutorial:

McCarthy, E.; Ekesan, Ş.; Giese, T. J.; Wilson, T. J.; Deng, J.; Huang, L.; Lilley, D. J.; York, D. M. Catalytic mechanism and pH dependence of a methyltransferase ribozyme (MTR1) from computational enzymology. Nucleic Acids Research 2023, 51, 4508–4518. (https://doi.org/10.1093/nar/gkad260)
Deng, J.; Wilson, T. J.; Wang, J.; Peng, X.; Li, M.; Lin, X.; Liao, W.; Lilley, D. M. J.; Huang, L. Structure and mechanism of a methyltransferase ribozyme. Nature Chemical Biology 2022, 18, 556–564. (https://doi.org/10.1038/s41589-022-00982-z)
Zgarbová, M.; Otyepka, M.; Šponer, J.; Mládek, A.; Banáš, P.; Cheatham, T. E., III; Jurečka, P. Refinement of the Cornell et al. nucleic acids force field based on reference quantum chemical calculations of glycosidic torsion profiles. Journal of Chemical Theory and Computation 2011, 7, 2886–2902. (https://doi.org/10.1021/ct200162x)
Horn, H. W.; Swope, W. C.; Pitera, J. W.; Madura, J. D.; Dick, T. J.; Hura, G. L.; Head-Gordon, T. Development of an improved four-site water model for biomolecular simulations: TIP4P-Ew. Journal of Chemical Physics 2004, 120, 9665–9678. (https://doi.org/10.1063/1.1683075)

Tutorial

Starting structure and input files 

This tutorial starts from pre-generated AMBER topology and coordinate files that define the initial structure of the fully solvated RNA system with counterions. These files serve as the starting point for all subsequent equilibration and production molecular dynamics simulations performed in this tutorial.

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 tutorial are:

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

Accessing the tutorial files 

All input files and reference output files required for this tutorial are provided in two ways.

Option 1: On the cluster

For convenience during the workshop, the full directory associated with Hands-On Session 2 is already available on the cluster and can be copied directly to your working directory from:

/path/Day1/HS2/Files

This directory contains all input files as well as reference output files corresponding to the equilibration and production steps described in this tutorial.

Option 2: Downloadable archive

All input files and reference output files can also be downloaded as a single compressed archive:

Download Files.tar.gz

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

1. Equilibrating the Box 

Copy the inputs to your working directory:

mkdir MTR1_Equil
cd MTR1_Equil
cp -r /path/Day1/HS2/Files/input/Equil/* ./

Now you will learn how to equilibrate the simulation box. Electrostatic interactions are extremely important for RNA structure given the negatively charged backbone. Therefore, a rigorous equilibration process allows the solvent to properly organize around the molecule in a controlled way while avoiding large fluctuations in density that could cause errors. This is especially important when your system contains divalent metal ions. The electrostatic interaction energy is proportional to the product of the two charges, so increasing the charge will lead to a larger interaction energy at a given distance. Nucleic acids are highly charged systems and must be equilibrated slowly.

You have been provided the Equil directory containing eqFromStratch.slurm and a directory called mdins.

List the contents of the mdins directory:

ls mdins
000_Min.mdin  003_NPT.mdin  03_Equil_NPT_solvent.mdin  06_Equil_solute_Rwt25.mdin  09_Equil_solute_Rwt2.mdin
001_NPT.mdin  01_Min.mdin   04_Min_solute.mdin         07_Equil_solute_Rwt10.mdin
002_NPT.mdin  02_Heat.mdin  05_Heat_solute_Rwt25.mdin  08_Equil_solute_Rwt5.mdin

These are the input files for the stepwise equilibration procedure and subsequent production simulation.

Note

The general approach is to slowly equilibrate the solvent around the solute, then gradually allow the solute to move. We cycle through minimizations to remove clashes and find a local minimum, heating to explore new conformations, then simulations in the NPT ensemble to approach the global minimum. There are some specific points to note about how the provided mdin files have been written.

1. 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.

You have been provided the production MD input file:

11_Production_0.5ns.mdin

This input file 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.

You have also been provided the SLURM submission script used to run the production simulation:

prod.slurm

The production SLURM script 

The following SLURM script is used to run the production simulation. It takes as input the prefix of the AMBER topology and restart files and automatically increments the output file numbering.

#!/bin/bash
#SBATCH -o slurm-%j.out
#SBATCH -e slurm-%j.err
#SBATCH --nodes=1
#SBATCH --gpus=1
#SBATCH --ntasks-per-node=1
#SBATCH --partition=gpu-shared
#SBATCH -t 0-00:20:00
#SBATCH --no-requeue
#SBATCH --reservation=amber24gpu
#SBATCH --account=gue998

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

mdin=mdins/11_Production_0.5ns.mdin

module load workshop/amber24/default

newname=""

if [[ $Base == *-run-* ]]; then
  prefix=${Base%-*}
  suffix=${Base##*-}
  if [[ $suffix =~ [0-9]+ ]]; then
    number=${suffix#*-}
    new_number=$(printf "%02d" $((number + 1)))
    newname="${prefix}-${new_number}"
  else
    newname="${Base}-01"
  fi
else
  newname="${Base}-01"
fi

exe=${AMBERHOME}/bin/pmemd.cuda_SPFP

$exe -O -p ${Parm} -c ${Base}.rst7 -i ${mdin} \
     -o ${newname}.mdout -inf ${newname}.mdinfo \
     -x ${newname}.nc -r ${newname}.rst7 -ref ${Base}.rst7

This script appends -01 to the input prefix when generating output files. If the input prefix already contains a numeric suffix (e.g. when restarting a simulation), the suffix is incremented to -02, -03, and so on. This allows longer production simulations to be split into multiple consecutive runs to accommodate HPC wall-time limits.

Preparing the production input files 

Copy the equilibrated restart file and rename it for production:

cp 09_Equil_solute_Rwt2-run.rst7 MTR1_counterIon_water_bulk-run.rst7

Production MD will be run with a 4 fs timestep. To do this stably, hydrogen mass repartitioning (HMR) is applied. You have been provided the cpptraj input file used for this purpose:

hmr.in

The cpptraj input file 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.

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.

Wrapping and visualization of the production trajectory 

After the production simulation has completed, the resulting trajectory must be wrapped to account for periodic boundary conditions before visualization.

You have been provided a cpptraj input file for wrapping the production trajectory:

wrap_RNA_prod.in

Wrap the structure:

cpptraj -i wrap_RNA_prod.in

This command produces the wrapped trajectory file centered.MTR1_counterIon_water_bulk-run-01.nc. This trajectory can now be visualized in VMD.

Open the trajectory in VMD:

vmd centered.MTR1_counterIon_water_bulk-run-01.nc \
    MTR1_counterIon_water_bulk-run.parm7

To use the same visualization representations as in the equilibration section, you can source the VMD representation file provided for this tutorial. (/path/Day1/HS2/Files/input/build/MTR_waterreps.txt)

MTR_waterreps.txt

Alternatively, you can load the trajectory and apply the representations directly from the command line:

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.

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.

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

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.

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.

You have been provided the original equilibration SLURM script:

eqFromStratch.slurm

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

[user@cluster] 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 -o slurm-%j.out
#SBATCH -e slurm-%j.err
#SBATCH --nodes=1
#SBATCH --ntasks-per-node=1
#SBATCH --partition=shared
#SBATCH --account=gue998
#SBATCH --reservation=amber24
#SBATCH -t 0-00:30:00
#SBATCH --no-requeue

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.

Before submitting the job, verify that the Equil directory contains the necessary files. If you are also running production simulations in the same directory, that is fine. However, if you previously copied 000_Min-run.mdout from the reference output, you may need to delete it, otherwise the script will detect that the step has already completed and exit immediately.

The directory should contain the following:

[user@cluster] ls
MTR1_counterIon_water_bulk.parm7
MTR1_counterIon_water_bulk.rst7
eqFromStratch.slurm
eqFromStratch_000.slurm
mdins
restrnt.disang

Submit the job:

[user@cluster] 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:

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

Plot the energy using xmgrace. Make sure the appropriate modules are loaded first. It may be more convenient to run xmgrace in a separate terminal:

[user@cluster] 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.