Exercise: Molecular Dynamics with OpenMM

Inspect the system

Go to the website of the RCSB protein data bank and download the solution NMR structure of a small peptide with the PDB-ID 6a5j as .pdb-textfile. Open the .pdb-file with the text editor of your choice and identify the section in which the molecular structure is described.

  • a) How many residues and atoms does the molecule have?

  • b) Copy the atom definitions for the second amino acid residue and include it in your report.

Open the .pdb-file in VMD.

  • c) How many structures are contained in the file?

  • d) Set the background colour to white and remove the axis icon. Visualise the peptide structure in a representation that emphasises the secondary structure (e.g. NewCartoon, NewRibbons, …) and choose Secondary Structure as colouring method. Visualise the side chain of lysines in the structure in a representation that shows bonds and atoms (e.g. Licorice, CPK, …) and color them by chemical element. Render an image of the first structure frame and include it in your report.

  • e) The structure file contains hydrogens. What is the total charge of the peptide?

Prepare the structure

Save the first structure included in the .pdb-file to a new separate .pdb file, e.g. conf.pdb. We want to use this frame as starting structure for our simulation. Use OpenMM to load the structure into memory.

  • a) Investigate the topology attribute of the loaded molecule. How many covalent bonds are present?

  • b) We want to simulate the molecule using the AMBER 14SB force field in TIP3P water. Create a forcefield object suitable for this purpose.

  • c) Our structure file does not include solvent and has a non-zero total charge. Put the molecule in a cubic water box with a minimum distance of 1 nm to the box borders. Make sure that also counter ions are added to neutralise the system.

  • d) Save the modified structure (peptide plus solvent and ions) to a .pdb-file. Load the structure into VMD and visualise the peptide as before. Additionally show ions as van-der-Waals spheres and water molecules within a distance of 2 Å in the Licorice representation. Optional Add a transparent surface representation (e.g. QuickSurf) for the peptide.

Simulation

Setup

  • a) Create a system object from the topology of the solvated peptide with the forcefield you setup earlier. Use H-bonds constraints, and PME for nonbonded interactions (cutoff 1 nm).

  • b) We want to simulate in the NVT ensemble. Add an Andersen thermostat to the system to keep the temperature at 300 K (coupling rate 1/ps).

  • c) We want to propagate the simulation using a Verlet integrator in timesteps of 2 fs (we need a contraint tolerance of 0.00001). Setup a simulation object with the solvated topology, the system and the integrator. Use "CPU" as the simulation platform.

  • d) Add the current atomic positions of the solvated peptide to the context of the simulation.

Energy minimisation

  • a) Minimise the energy of the system. Save the new positions to a .pdb-file.

  • b) Open the structure in VMD and render an image the same way as in the last step in which the original and the minimised atoms are overlaid. Choose different colors for the two states. Can you see the difference?

Equilibration

  • a) We now want to equilibrate our system in the NVT ensemble for 100 ps. Add a state reporter to the simulation object that reports the total and potential energy and the temperature at every picosecond and write it to a file (e.g. equilibration.log).

  • b) Run the simulation. It should take about 10 minutes to complete on a processor with 4 cores. Reduce the number of steps, if it takes much longer on your hardware.

  • c) Plot the energies against the simulation time step.

  • d) Plot the temperature against the simulation time step.

  • e) Save the current state of the simulation to an .xml-file.

Production

  • a) Load the state of the equilibrated system into the simulation object.

  • b) Reset the simulation reporters. Let’s now assume we want to simulate the system for 100 ns. Add a state reporter to the simulation object that reports the potential energy every 10 ps and writes it to a file (e.g. simulation.log). Add a structure reporter that saves the current atom positions to a .pdb-file every 10 ps.

  • c) The simulation can easily take several days to finish, depending on the used hardware, so you do not have to run it yourself completely. Just start the simulation to see if everything works and run it until the first structure has been saved by the structure reporter (after 10 ps).

  • d) Open the saved .pdb-file in VMD. Render an image that shows only the peptide in a representation which emphasises the secondary structure and compares the original to the minimised and the current conformation. Align the structures (e.g. using the RMSD Trajectory Tool).