Contents

Note This tutorial should be run on Google Colab and not on JupyterHub/noto.epfl.ch.

This is due to the fact that Colab offers a free GPU for the simulation whereas on JupyterHub the simulation runs on the CPU only which is much slower.

On Colab: Change the runtime to use a GPU by selecting Runtime -> Change Runtime type -> select GPU

6.3. Practical: TRP Cage using OpenMM

We will use a modern and open source simulation program to run our molecular dynamics simulation of the TRP cage protein: OpenMM

OpenMM is different from most other simulations tools in that it is not controlled by executing different subprograms in the command line but instead one uses a python program to control the flow of our simulation. This allows us to easily rerun parts or the whole simulation.

Below is a detailed description of our workflow. The respective parts you need to change during the exercise are indicated.

6.3.1. Imports

First we need to install the necessary modules and import them. First we install the conda package manager which is not available by default on Colab. The runtime will display a notification that it crashed, this is normal and it will automatically restart. Then, we can install our simulation package Openmm and parmed, which we use to load our Amber trajectory files.

# Uncomment below code on Colab

#!pip install -q condacolab
#import condacolab
#condacolab.install()

# Uncomment on google colab 
#!pip install py3Dmol
#!conda install openmm ambertools=20 &> /dev/null  
# The last part simply makes the command silent. Remove &> /dev/null  if there are problems

import py3Dmol

6.3.2. Creating our system

Creating the system will be done using tleap which is part of the Amber software package. Linear peptides can be easily created using this tool by specifying the sequence of our TRP cage miniprotein which is:

NLYIQWLKDGGPSSGRPPPS

Below we create this file using python and save it as leap.in to the filesystem.

If you open the sidebar on Colab you can check whether it was correctly created and view its contents.

Then we execute tleap using our input file using the command !tleap -f leap.in.

  leapin = """source leaprc.protein.ff19SB
  set default PBradii mbondi3
  TRP = sequence {  NASN LEU TYR ILE GLN TRP LEU LYS ASP GLY GLY PRO SER SER GLY ARG PRO PRO PRO CSER}
  savepdb TRP trp_cage_linear.pdb
  saveamberparm TRP trp_cage.prmtop trp_cage_linear.nc
  quit
  """

  with open("leap.in","w+") as f:
    f.writelines(leapin)

6.3.3. Execute tleap

This will create our input files for the MD simulation.

!tleap -f leap.in

-I: Adding /home/lcbc/miniconda3/envs/jupyterbook/dat/leap/prep to search path.
-I: Adding /home/lcbc/miniconda3/envs/jupyterbook/dat/leap/lib to search path.
-I: Adding /home/lcbc/miniconda3/envs/jupyterbook/dat/leap/parm to search path.
-I: Adding /home/lcbc/miniconda3/envs/jupyterbook/dat/leap/cmd to search path.
-f: Source leap.in.

Welcome to LEaP!
(no leaprc in search path)
Sourcing: ./leap.in
----- Source: /home/lcbc/miniconda3/envs/jupyterbook/dat/leap/cmd/leaprc.protein.ff19SB
----- Source of /home/lcbc/miniconda3/envs/jupyterbook/dat/leap/cmd/leaprc.protein.ff19SB done
Log file: ./leap.log

Loading parameters: /home/lcbc/miniconda3/envs/jupyterbook/dat/leap/parm/parm19.dat
Reading title:
PARM99 + frcmod.ff99SB + frcmod.parmbsc0 + OL3 for RNA + ff19SB
Loading parameters: /home/lcbc/miniconda3/envs/jupyterbook/dat/leap/parm/frcmod.ff19SB
Reading force field modification type file (frcmod)
Reading title:
ff19SB AA-specific backbone CMAPs for protein 07/25/2019
Loading library: /home/lcbc/miniconda3/envs/jupyterbook/dat/leap/lib/amino19.lib
Loading library: /home/lcbc/miniconda3/envs/jupyterbook/dat/leap/lib/aminoct12.lib
Loading library: /home/lcbc/miniconda3/envs/jupyterbook/dat/leap/lib/aminont12.lib
Using ArgH and AspGluO modified Bondi2 radii
Writing pdb file: trp_cage_linear.pdb

/home/lcbc/miniconda3/envs/jupyterbook/bin/teLeap: Warning!
 Converting N-terminal residue name to PDB format: NASN -> ASN

/home/lcbc/miniconda3/envs/jupyterbook/bin/teLeap: Warning!
 Converting C-terminal residue name to PDB format: CSER -> SER
Checking Unit.

/home/lcbc/miniconda3/envs/jupyterbook/bin/teLeap: Warning!
The unperturbed charge of the unit (1.000000) is not zero.

/home/lcbc/miniconda3/envs/jupyterbook/bin/teLeap: Note.
Ignoring the warning from Unit Checking.

Building topology.
Building atom parameters.
Building bond parameters.
Building angle parameters.
Building proper torsion parameters.
Building improper torsion parameters.
 total 63 improper torsions applied
Building H-Bond parameters.
Incorporating Non-Bonded adjustments.
Not Marking per-residue atom chain types.
Marking per-residue atom chain types.
  (Residues lacking connect0/connect1 - 
   these don't have chain types marked:

	res	total affected

	CSER	1
	NASN	1
  )
 (no restraints)
	Quit

Exiting LEaP: Errors = 0; Warnings = 3; Notes = 1.

You will be greeted with plenty output, lets go through it step by step.

>$ tleap
-I: Adding /home/packages/amber20/dat/leap/prep to search path.
-I: Adding /home/packages/amber20/dat/leap/lib to search path.
-I: Adding /home/packages/amber20/dat/leap/parm to search path.
-I: Adding /home/packages/amber20/dat/leap/cmd to search path.
Welcome to LEaP!
>(no leaprc in search path)

This is leap starting up telling us that it has added several directories to its search path. These directories contain the parameters of the amber force field that we now will load.

Next we need to load a specific force field. We do this using the following command:

source leaprc.protein.ff19SB

If you want you can take a look at this file by creating a new cell and using the following command:

!cat  /usr/local/dat/leap/parm/parm19.dat

Here you will find the masses, bond distance and force constants etc. that make up our forcefield.

This produces the following output which are all force field parameters that are loaded.

----- Source: /usr/local/dat/leap/cmd/leaprc.protein.ff19SB
----- Source of /usr/local/dat/leap/cmd/leaprc.protein.ff19SB done
Log file: ./leap.log
Loading parameters: /usr/local/dat/leap/parm/parm19.dat
Reading title:
PARM99 + frcmod.ff99SB + frcmod.parmbsc0 + OL3 for RNA + ff19SB
Loading parameters: /usr/local/dat/leap/parm/frcmod.ff19SB
Reading force field modification type file (frcmod)
Reading title:
ff19SB AA-specific backbone CMAPs for protein 07/25/2019
Loading library: /usr/local/dat/leap/lib/amino19.lib
Loading library: /usr/local/dat/leap/lib/aminoct12.lib
Loading library: /usr/local/dat/leap/lib/aminont12.lib

In the next line, we set the radii of the atoms to be used for the implicit solvation.

set default PBradii mbondi3

This produces the following line in the leap output:

Using ArgH and AspGluO modified Bondi2 radii

Next, we create a new UNIT in AMBER using the sequence command.

TRP = sequence {  NASN LEU TYR ILE GLN TRP LEU LYS ASP GLY GLY PRO SER SER GLY ARG PRO PRO PRO CSER}

Then, we save a PDB file of our system.

savepdb TRP trp_cage_linear.pdb

This creates the pdb file on disk and produces the following output:

Writing pdb file: trp_cage_linear.pdb

/usr/local/bin/teLeap: Warning!
 Converting N-terminal residue name to PDB format: NASN -> ASN

/usr/local/bin/teLeap: Warning!
 Converting C-terminal residue name to PDB format: CSER -> SER
Checking Unit.

As you can see, the N and C-terminal residues have a different residue namer in AMBER compared e.g to a normal serine. This is due to the fact that here the force field parameters slightly deviate from the parameters of a serine involved in two peptide bonds versus one single peptide bond for a terminal uncapped serine.

In the last step, we now generate the parameters in AMBER format and the input coordinates in the NetCDF format (which is a binary format).

saveamberparm TRP trp_cage.prmtop trp_cage_linear.nc

This produces plenty of output and a warning:

/usr/local/bin/teLeap: Warning!
The unperturbed charge of the unit (1.000000) is not zero.

/usr/local/bin/teLeap: Note.
Ignoring the warning from Unit Checking.

Building topology.
Building atom parameters.
Building bond parameters.
Building angle parameters.
Building proper torsion parameters.
Building improper torsion parameters.
 total 63 improper torsions applied
Building H-Bond parameters.
Incorporating Non-Bonded adjustments.
Not Marking per-residue atom chain types.
Marking per-residue atom chain types.
  (Residues lacking connect0/connect1 - 
   these don't have chain types marked:

	res	total affected

	CSER	1
	NASN	1
  )
 (no restraints)

The first warning we can safely ignore as in the implicit solvent treatment we will use, we cannot use counter ions to neutralize our system. Instead we will use a screened charge model to emulate the presence of ions in the solvent. Next we see terms, that align with the functional form of the AMBER forcefield:

Again, we can look at the contents of the parameter file !cat trp_cage.prmtop and we see that it does not contain the coordinates anymore in contrast to the pdb file.

6.3.3.1. The starting structure visualized

with open("trp_cage_linear.pdb") as ifile:
    system = "".join([x for x in ifile])

view = py3Dmol.view(width=400, height=300)
view.addModelsAsFrames(system)
view.setStyle({'model': -1}, {"stick": {}})
view.zoomTo()
view.show()

You appear to be running in JupyterLab (or JavaScript failed to load for some other reason). You need to install the 3dmol extension:
jupyter labextension install jupyterlab_3dmol

6.3.4. Simulating the system

Now we can start the simulation by writing our simulation script in python.

First, lets import the required modules:

  • openmm for running the simulation

  • parmed to read/write the input and output files

# OpenMM Imports
import simtk.openmm as mm
import simtk.openmm.app as app

# ParmEd Imports
from parmed import load_file, unit as u
from parmed.openmm import StateDataReporter, NetCDFReporter

# System modules
import sys

6.3.4.1. Loading Amber files

In this stage, we simply instantiate the AmberParm object from the input topology and coordinate files. After this command, trp_cage_gas will contain a full description of every particle, the parameters defining their interactions, and their position.

# Load the Amber files
print('Loading AMBER files...')
trp_cage_gas = load_file('trp_cage.prmtop', 'trp_cage_linear.nc')
print('Loaded positions')
print(trp_cage_gas.positions)
print(trp_cage_gas.angles)
print(trp_cage_gas.dihedrals)

Loading AMBER files...
Loaded positions
[Vec3(x=3.32577, y=1.547909, z=-1.6e-06), Vec3(x=4.046154, y=0.839991, z=-2.9e-06), Vec3(x=2.823094, y=1.499508, z=-0.874687), Vec3(x=2.823097, y=1.499507, z=0.874685), Vec3(x=3.970048, y=2.845795, z=-1e-07), Vec3(x=3.671663, y=3.400129, z=-0.88982), Vec3(x=3.576965, y=3.653838, z=1.232143), Vec3(x=2.496995, y=3.801075, z=1.241379), Vec3(x=3.877484, y=3.115795, z=2.131197), Vec3(x=4.2537, y=5.017112, z=1.232144), Vec3(x=5.005299, y=5.340406, z=0.315072), Vec3(x=3.984885, y=5.817909, z=2.265917), Vec3(x=4.408015, y=6.733702, z=2.314743), Vec3(x=3.359611, y=5.504297, z=2.994464), Vec3(x=5.485541, y=2.705207, z=-4.4e-06), Vec3(x=6.008824, y=1.593175, z=-8.4e-06), Vec3(x=6.1910071, y=3.8385836, z=-5.1e-06), Vec3(x=5.7152868, y=4.7295323, z=-2.8e-06), Vec3(x=7.6400076, y=3.8385837, z=-7.4e-06), Vec3(x=8.0038585, y=3.3248391, z=0.8898118), Vec3(x=8.1890012, y=3.1272109, z=-1.2321515), Vec3(x=7.8406894, y=2.0944027, z=-1.2413883), Vec3(x=7.8406911, y=3.6356234, z=-2.1312048), Vec3(x=9.713384, y=3.148898, z=-1.1945863), Vec3(x=10.0624752, y=4.1814477, z=-1.1857965), Vec3(x=10.201034, y=2.4370787, z=0.0628467), Vec3(x=9.8527007, y=1.4042699, z=0.0545459), Vec3(x=11.290585, y=2.4522959, z=0.0901921), Vec3(x=9.8088607, y=2.9451708, z=0.9438364), Vec3(x=10.2623786, y=2.437527, z=-2.4267313), Vec3(x=9.9140694, y=2.9459412, z=-3.3257836), Vec3(x=11.3519542, y=2.452742, z=-2.4003739), Vec3(x=9.9140658, y=1.4047192, z=-2.4359701), Vec3(x=8.1879242, y=5.2585384, z=-2e-06), Vec3(x=7.4245351, y=6.2216985, z=3.4e-06), Vec3(x=9.5167782, y=5.3864922, z=-2e-06), Vec3(x=10.1032885, y=4.5642362, z=-5.3e-06), Vec3(x=10.1610562, y=6.6843782, z=1.2e-06), Vec3(x=9.8626701, y=7.2387139, z=-0.8898172), Vec3(x=9.7679746, y=7.4924189, z=1.2321464), Vec3(x=8.6880046, y=7.6396559, z=1.2413839), Vec3(x=10.0684947, y=6.9543742, z=2.131199), Vec3(x=10.4583376, y=8.8348479, z=1.1949511), Vec3(x=10.2509378, y=9.75749, z=2.2272871), Vec3(x=9.5911198, y=9.5067954, z=3.0578844), Vec3(x=10.8910087, y=11.0021261, z=2.1927997), Vec3(x=10.7295337, y=11.7204685, z=2.9965482), Vec3(x=11.7384805, y=11.3241231, z=1.1259783), Vec3(x=12.3602644, y=12.5331981, z=1.0924779), Vec3(x=12.1473373, y=13.0855637, z=1.8482252), Vec3(x=11.9458833, y=10.401483, z=0.0936413), Vec3(x=12.6057013, y=10.6521796, z=-0.736956), Vec3(x=11.3058143, y=9.156846, z=0.1281237), Vec3(x=11.4672934, y=8.4385045, z=-0.6756239), Vec3(x=11.6765492, y=6.5437902, z=-5.1e-06), Vec3(x=12.1998322, y=5.4317582, z=-1.19e-05), Vec3(x=12.3820153, y=7.6771668, z=-4.5e-06), Vec3(x=11.906295, y=8.5681155, z=0.0), Vec3(x=13.8310157, y=7.6771669, z=-8.4e-06), Vec3(x=14.1948677, y=7.1634206, z=0.8898093), Vec3(x=14.3393155, y=6.9582766, z=-1.2451795), Vec3(x=13.9762283, y=5.9305276, z=-1.2456538), Vec3(x=13.8309501, y=7.6772011, z=-2.4903046), Vec3(x=14.1940364, y=8.7049505, z=-2.4907693), Vec3(x=14.1939901, y=7.163751, z=-3.3806257), Vec3(x=12.7409506, y=7.6774893, z=-2.4907437), Vec3(x=15.8643157, y=6.9582741, z=-1.2451831), Vec3(x=16.2281675, y=7.9857529, z=-1.2451758), Vec3(x=16.2281673, y=6.4445292, z=-0.3553642), Vec3(x=16.3726164, y=6.2393856, z=-2.4903561), Vec3(x=16.0095299, y=6.7528612, z=-3.380643), Vec3(x=17.4626159, y=6.2390974, z=-2.490852), Vec3(x=16.0095275, y=5.2116375, z=-2.4908324), Vec3(x=14.3789324, y=9.0971216, z=-1.1e-06), Vec3(x=13.6155433, y=10.0602817, z=7.1e-06), Vec3(x=15.7077863, y=9.2250754, z=-2.3e-06), Vec3(x=16.2942967, y=8.4028194, z=-7.9e-06), Vec3(x=16.3520643, y=10.5229614, z=2.6e-06), Vec3(x=16.0536772, y=11.0772988, z=-0.8898145), Vec3(x=15.9589842, y=11.3309998, z=1.2321497), Vec3(x=14.8790142, y=11.4782367, z=1.2413888), Vec3(x=16.2595053, y=10.7929534, z=2.131201), Vec3(x=16.6562051, y=12.6867639, z=1.1945872), Vec3(x=17.7362901, y=12.5403399, z=1.1857951), Vec3(x=16.355798, y=13.2256253, z=0.295985), Vec3(x=16.288998, y=13.5256003, z=2.4103162), Vec3(x=15.520986, y=13.0855031, z=3.2629074), Vec3(x=16.8408782, y=14.73868, z=2.4883505), Vec3(x=16.630458, y=15.336199, z=3.2749863), Vec3(x=17.4663015, y=15.0533698, z=1.7603963), Vec3(x=17.8675573, y=10.3823734, z=-5.8e-06), Vec3(x=18.3908403, y=9.2703414, z=-1.53e-05), Vec3(x=18.5730235, y=11.51575, z=-3.8e-06), Vec3(x=18.0973032, y=12.4066987, z=2.9e-06), Vec3(x=20.0220239, y=11.5157501, z=-9.5e-06), Vec3(x=20.3858769, y=11.0020022, z=0.8898069), Vec3(x=20.5710147, y=10.804382, z=-1.2321576), Vec3(x=20.2227028, y=9.7715738, z=-1.2413975), Vec3(x=20.2227025, y=11.3127979, z=-2.1312082), Vec3(x=22.0775181, y=10.7530289, z=-1.3211062), Vec3(x=22.8092912, y=10.1917656, z=-2.293264), Vec3(x=22.3161361, y=9.7057441, z=-3.1350959), Vec3(x=24.1832075, y=10.3900395, z=-1.9498351), Vec3(x=24.9747361, y=10.0763641, z=-2.4931489), Vec3(x=24.249489, y=11.0441895, z=-0.8167882), Vec3(x=25.394902, y=11.4466851, z=-0.1196306), Vec3(x=26.3835782, y=11.217224, z=-0.5170795), Vec3(x=25.1750674, y=12.1329251, z=1.0690013), Vec3(x=26.0209601, y=12.4766372, z=1.6643395), Vec3(x=23.9224911, y=12.3846886, z=1.5050842), Vec3(x=23.7689132, y=12.9242424, z=2.4396402), Vec3(x=22.7832204, y=11.9693148, z=0.7856205), Vec3(x=21.7866856, y=12.190115, z=1.1680684), Vec3(x=23.0046322, y=11.2781377, z=-0.4115644), Vec3(x=20.5699405, y=12.9357048, z=-1e-07), Vec3(x=19.8065514, y=13.8988649, z=1.08e-05), Vec3(x=21.8987945, y=13.0636587, z=-2.6e-06), Vec3(x=22.4853049, y=12.2414027, z=-1.05e-05), Vec3(x=22.5430725, y=14.3615447, z=3.9e-06), Vec3(x=22.2446844, y=14.9158837, z=-0.8898117), Vec3(x=22.1499938, y=15.1695807, z=1.2321531), Vec3(x=21.0700239, y=15.3168176, z=1.2413937), Vec3(x=22.450516, y=14.6315325, z=2.1312029), Vec3(x=22.8472147, y=16.5253449, z=1.1945923), Vec3(x=23.9272997, y=16.3789209, z=1.1857987), Vec3(x=22.4264533, y=17.278645, z=-0.0628359), Vec3(x=21.3464733, y=17.425863, z=-0.0545313), Vec3(x=22.9245372, y=18.2478022, z=-0.0901782), Vec3(x=22.7071802, y=16.701459, z=-0.9438292), Vec3(x=22.454138, y=17.3333809, z=2.4267423), Vec3(x=22.7546622, y=16.7953328, z=3.3257912), Vec3(x=22.9522309, y=18.302561, z=2.400388), Vec3(x=21.3741681, y=17.4806168, z=2.4359849), Vec3(x=24.0585655, y=14.2209567, z=-6.5e-06), Vec3(x=24.5818485, y=13.1089247, z=-1.87e-05), Vec3(x=24.7640316, y=15.3543333, z=-3.2e-06), Vec3(x=24.2883113, y=16.2452819, z=5.8e-06), Vec3(x=26.213032, y=15.3543334, z=-1.06e-05), Vec3(x=26.5768862, y=14.8405837, z=0.8898044), Vec3(x=26.7620214, y=14.6429676, z=-1.2321607), Vec3(x=26.4137095, y=13.6101594, z=-1.2414021), Vec3(x=26.4137081, y=15.1513852, z=-2.1312099), Vec3(x=28.2864043, y=14.6646544, z=-1.1946007), Vec3(x=28.6354955, y=15.6972041, z=-1.1858063), Vec3(x=28.6354982, y=14.1559766, z=-0.2960006), Vec3(x=28.8353945, y=13.9532904, z=-2.4267517), Vec3(x=28.4870817, y=12.9204827, z=-2.4359952), Vec3(x=28.4870821, y=14.4617098, z=-3.3257999), Vec3(x=30.3597779, y=13.9749759, z=-2.3891898), Vec3(x=30.70887, y=15.0075251, z=-2.3803943), Vec3(x=30.7088696, y=13.4662958, z=-1.4905916), Vec3(x=30.8889677, y=13.2892695, z=-3.5769035), Vec3(x=30.5659591, y=13.7607247, z=-4.4096671), Vec3(x=31.8985592, y=13.3036302, z=-3.5520265), Vec3(x=30.5659568, y=12.3323488, z=-3.5850196), Vec3(x=26.7609487, y=16.7742881, z=9e-07), Vec3(x=25.9975596, y=17.7374481, z=1.45e-05), Vec3(x=28.0898027, y=16.9022419, z=-3e-06), Vec3(x=28.676313, y=16.0799859, z=-1.31e-05), Vec3(x=28.7340807, y=18.2001279, z=5.3e-06), Vec3(x=28.4356915, y=18.7544686, z=-0.889809), Vec3(x=28.3410034, y=19.0081615, z=1.2321564), Vec3(x=27.2610335, y=19.1553985, z=1.2413986), Vec3(x=28.6415267, y=18.4701117, z=2.1312049), Vec3(x=29.0391393, y=20.3657038, z=1.1945474), Vec3(x=28.4331436, y=21.3092682, z=0.6200414), Vec3(x=30.1717718, y=20.4461987, z=1.7406827), Vec3(x=30.2495737, y=18.0595399, z=-7.1e-06), Vec3(x=30.7728567, y=16.9475079, z=-2.21e-05), Vec3(x=30.9550398, y=19.1929165, z=-2.6e-06), Vec3(x=30.4793195, y=20.0838652, z=8.7e-06), Vec3(x=32.4040402, y=19.1929166, z=-1.16e-05), Vec3(x=32.7678954, y=18.6791653, z=0.8898019), Vec3(x=32.7678847, y=18.6791867, z=-0.8898421), Vec3(x=32.9345676, y=20.6194592, z=2.1e-06), Vec3(x=32.1594688, y=21.5732218, z=1.83e-05), Vec3(x=34.2617591, y=20.7636388, z=-2.9e-06), Vec3(x=34.8582716, y=19.9486098, z=-1.53e-05), Vec3(x=34.8901321, y=22.0692994, z=7.2e-06), Vec3(x=34.5849916, y=22.6199549, z=-0.8898056), Vec3(x=34.5850071, y=22.6199341, z=0.8898384), Vec3(x=36.4056248, y=21.9287114, z=-7.3e-06), Vec3(x=36.9289091, y=20.8166796, z=-2.49e-05), Vec3(x=37.111091, y=23.062088, z=-1.4e-06), Vec3(x=36.4673778, y=24.3702918, z=1.79e-05), Vec3(x=36.1600381, y=24.6245613, z=1.0144092), Vec3(x=35.5918554, y=24.3450062, z=-0.6487657), Vec3(x=37.5399055, y=25.2797951, z=-0.521985), Vec3(x=37.4958194, y=26.2264683, z=0.0164993), Vec3(x=37.3897655, y=25.4604948, z=-1.586366), Vec3(x=38.8474038, y=24.5657281, z=-0.2756996), Vec3(x=39.343812, y=25.0126395, z=0.5856645), Vec3(x=39.4836402, y=24.6684176, z=-1.154766), Vec3(x=38.5620908, y=23.0620881, z=-1.22e-05), Vec3(x=38.9602182, y=22.5667425, z=-0.8855773), Vec3(x=39.1086404, y=22.3843484, z=1.2483615), Vec3(x=38.346055, y=21.8704502, z=2.0637218), Vec3(x=40.4354217, y=22.3843549, z=1.396268), Vec3(x=41.0213467, y=22.8199592, z=0.6983853), Vec3(x=41.0782314, y=21.7717283, z=2.5412948), Vec3(x=40.7309927, y=22.2543195, z=3.4548758), Vec3(x=40.7534343, y=20.283924, z=2.6223488), Vec3(x=39.6755642, y=20.151372, z=2.7157583), Vec3(x=41.1033072, y=19.7853514, z=1.718405), Vec3(x=41.4038056, y=19.7154753, z=3.761991), Vec3(x=41.1993442, y=18.77889, z=3.8130142), Vec3(x=42.5914678, y=21.9137696, z=2.4611454), Vec3(x=43.1143703, y=22.4900044, z=1.509848), Vec3(x=43.295434, y=21.3856156, z=3.464993), Vec3(x=42.8201046, y=20.9212472, z=4.2256002), Vec3(x=44.7422265, y=21.4547893, z=3.5051009), Vec3(x=45.1548162, y=20.9552032, z=2.6285835), Vec3(x=45.2226604, y=22.9021164, z=3.5123194), Vec3(x=44.8747569, y=23.4054773, z=2.6102716), Vec3(x=44.8250804, y=23.4130952, z=4.3892064), Vec3(x=46.6518404, y=22.9297262, z=3.5521115), Vec3(x=46.9542773, y=23.8408306, z=3.556657), Vec3(x=45.2887762, y=20.777048, z=4.7534735), Vec3(x=44.5261893, y=20.263148, z=5.5688313), Vec3(x=46.6155574, y=20.7770556, z=4.9013806), Vec3(x=47.2014823, y=21.2126618, z=4.2034991), Vec3(x=47.2583672, y=20.1644275, z=6.0464066), Vec3(x=46.9111276, y=20.6470166, z=6.9599883), Vec3(x=47.0093358, y=19.1037627, z=6.0791618), Vec3(x=48.7697822, y=20.3155565, z=5.9501363), Vec3(x=49.2790837, y=20.8942712, z=4.992983), Vec3(x=49.4875251, y=19.7923744, z=6.9468061), Vec3(x=49.0230977, y=19.3256772, z=7.7127044), Vec3(x=50.9342162, y=19.8703203, z=6.9715369), Vec3(x=51.3405143, y=19.3726199, z=6.091017), Vec3(x=51.4059322, y=21.3205296, z=6.9727376), Vec3(x=51.0454538, y=21.8211429, z=6.0741048), Vec3(x=51.0145776, y=21.8297134, z=7.8534608), Vec3(x=52.9302149, y=21.3592158, z=6.9990019), Vec3(x=53.2914437, y=20.85916, z=7.8976446), Vec3(x=53.3223211, y=20.8505922, z=6.1182881), Vec3(x=53.4019318, y=22.809425, z=7.0002046), Vec3(x=53.0414524, y=23.3100398, z=6.1015732), Vec3(x=53.0105782, y=23.3186072, z=7.8809289), Vec3(x=54.87971, y=22.8864583, z=7.0254756), Vec3(x=55.407037, y=22.0251518, z=7.0387715), Vec3(x=55.5450355, y=24.0380686, z=7.0314766), Vec3(x=54.9211681, y=25.2125551, z=7.0149989), Vec3(x=53.9118302, y=25.2445701, z=6.9973309), Vec3(x=55.4596084, y=26.0670458, z=7.0202276), Vec3(x=56.8740084, y=23.9912734, z=7.0547596), Vec3(x=57.344659, y=23.0977228, z=7.0672264), Vec3(x=57.4064711, y=24.8495019, z=7.0598688), Vec3(x=51.4980495, y=19.1967744, z=8.2144818), Vec3(x=50.7472657, y=18.6788487, z=9.0382032), Vec3(x=52.8262986, y=19.2049059, z=8.3483217), Vec3(x=53.6967467, y=19.8259937, z=7.357194), Vec3(x=53.7271373, y=20.9028145, z=7.5234282), Vec3(x=53.3121237, y=19.6239176, z=6.3575294), Vec3(x=55.0216051, y=19.1673865, z=7.6041471), Vec3(x=55.8134735, y=19.9068467, z=7.4848077), Vec3(x=55.1727391, y=18.3558549, z=6.8923351), Vec3(x=54.9674728, y=18.6408169, z=9.0183233), Vec3(x=55.5428887, y=19.2994189, z=9.6688914), Vec3(x=55.3964511, y=17.6391372, z=9.0450648), Vec3(x=53.485805, y=18.5961422, z=9.4884351), Vec3(x=53.2752477, y=17.5282636, z=9.5467368), Vec3(x=53.0628124, y=19.2675559, z=10.7871904), Vec3(x=52.2143457, y=20.1566675, z=10.7821979), Vec3(x=53.6577034, y=18.839179, z=11.9029068), Vec3(x=54.6619712, y=17.7823743, z=11.8835731), Vec3(x=55.6257575, y=18.2015548, z=11.5946025), Vec3(x=54.3698016, y=17.0158817, z=11.1657813), Vec3(x=54.6539477, y=17.2763554, z=13.2956211), Vec3(x=55.6757336, y=17.0431305, z=13.5950519), Vec3(x=54.0411805, y=16.3776465, z=13.3659072), Vec3(x=54.0744783, y=18.390702, z=14.1338017), Vec3(x=54.8804554, y=18.8960038, z=14.665939), Vec3(x=53.3720465, y=17.9704601, z=14.8535836), Vec3(x=53.3416167, y=19.3988294, z=13.2037838), Vec3(x=52.291399, y=19.2519112, z=13.4558871), Vec3(x=53.6719609, y=20.8836686, z=13.2548696), Vec3(x=54.0845481, y=21.4656847, z=12.2541352), Vec3(x=53.4892248, y=21.496729, z=14.4266159), Vec3(x=52.9984553, y=20.7846628, z=15.6004434), Vec3(x=53.8220222, y=20.2474837, z=16.0708438), Vec3(x=52.2273823, y=20.0753592, z=15.299721), Vec3(x=52.4581208, y=21.8836348, z=16.466641), Vec3(x=52.6903075, y=21.6615088, z=17.5082016), Vec3(x=51.3776307, y=21.9579132, z=16.343657), Vec3(x=53.1411452, y=23.1495776, z=16.00733), Vec3(x=53.9210041, y=23.4158632, z=16.7207801), Vec3(x=52.4077657, y=23.9542185, z=15.9543456), Vec3(x=53.7677284, y=22.9096057, z=14.6044233), Vec3(x=53.1504095, y=23.5227891, z=13.9479008), Vec3(x=55.2368605, y=23.2135681, z=14.3480211), Vec3(x=56.0015099, y=22.3265996, z=13.9751424), Vec3(x=55.6306209, y=24.4732139, z=14.5491932), Vec3(x=54.9645167, y=25.1676327, z=14.8560923), Vec3(x=57.002764, y=24.8892139, z=14.339982), Vec3(x=57.6597774, y=24.3203003, z=14.9978345), Vec3(x=57.4342894, y=24.6522328, z=12.8966351), Vec3(x=57.352448, y=23.590941, z=12.6619845), Vec3(x=56.7913307, y=25.2226501, z=12.2263156), Vec3(x=58.7904189, y=25.0738914, z=12.7292498), Vec3(x=59.062067, y=24.9247111, z=11.8206506), Vec3(x=57.1779785, y=26.3709208, z=14.6405292), Vec3(x=56.2219886, y=27.0501032, z=15.0082724), Vec3(x=58.2775843, y=26.9062427, z=14.5191329)] A
TrackedList([
	<Angle <Atom O [15]; In ASN 0>--<Atom C [14]; In ASN 0>--<Atom N [16]; In LEU 1>; type=<AngleType; k=80.000, theteq=122.900>>
	<Angle <Atom C [14]; In ASN 0>--<Atom N [16]; In LEU 1>--<Atom CA [18]; In LEU 1>; type=<AngleType; k=50.000, theteq=121.900>>
	<Angle <Atom OD1 [10]; In ASN 0>--<Atom CG [9]; In ASN 0>--<Atom ND2 [11]; In ASN 0>; type=<AngleType; k=80.000, theteq=122.900>>
	<Angle <Atom CB [6]; In ASN 0>--<Atom CA [4]; In ASN 0>--<Atom C [14]; In ASN 0>; type=<AngleType; k=63.000, theteq=111.100>>
	<Angle <Atom CB [6]; In ASN 0>--<Atom CG [9]; In ASN 0>--<Atom OD1 [10]; In ASN 0>; type=<AngleType; k=80.000, theteq=120.400>>
	<Angle <Atom CB [6]; In ASN 0>--<Atom CG [9]; In ASN 0>--<Atom ND2 [11]; In ASN 0>; type=<AngleType; k=70.000, theteq=116.600>>
	<Angle <Atom CA [4]; In ASN 0>--<Atom CB [6]; In ASN 0>--<Atom CG [9]; In ASN 0>; type=<AngleType; k=63.000, theteq=111.100>>
	<Angle <Atom CA [4]; In ASN 0>--<Atom C [14]; In ASN 0>--<Atom O [15]; In ASN 0>; type=<AngleType; k=80.000, theteq=120.400>>
	<Angle <Atom CA [4]; In ASN 0>--<Atom C [14]; In ASN 0>--<Atom N [16]; In LEU 1>; type=<AngleType; k=70.000, theteq=116.600>>
	<Angle <Atom N [0]; In ASN 0>--<Atom CA [4]; In ASN 0>--<Atom CB [6]; In ASN 0>; type=<AngleType; k=80.000, theteq=111.200>>
	<Angle <Atom N [0]; In ASN 0>--<Atom CA [4]; In ASN 0>--<Atom C [14]; In ASN 0>; type=<AngleType; k=80.000, theteq=111.200>>
	<Angle <Atom O [34]; In LEU 1>--<Atom C [33]; In LEU 1>--<Atom N [35]; In TYR 2>; type=<AngleType; k=80.000, theteq=122.900>>
	<Angle <Atom C [33]; In LEU 1>--<Atom N [35]; In TYR 2>--<Atom CA [37]; In TYR 2>; type=<AngleType; k=50.000, theteq=121.900>>
	<Angle <Atom CD1 [25]; In LEU 1>--<Atom CG [23]; In LEU 1>--<Atom CD2 [29]; In LEU 1>; type=<AngleType; k=40.000, theteq=109.500>>
	<Angle <Atom CB [20]; In LEU 1>--<Atom CA [18]; In LEU 1>--<Atom C [33]; In LEU 1>; type=<AngleType; k=63.000, theteq=111.100>>
	<Angle <Atom CB [20]; In LEU 1>--<Atom CG [23]; In LEU 1>--<Atom CD1 [25]; In LEU 1>; type=<AngleType; k=40.000, theteq=109.500>>
	<Angle <Atom CB [20]; In LEU 1>--<Atom CG [23]; In LEU 1>--<Atom CD2 [29]; In LEU 1>; type=<AngleType; k=40.000, theteq=109.500>>
	<Angle <Atom CA [18]; In LEU 1>--<Atom CB [20]; In LEU 1>--<Atom CG [23]; In LEU 1>; type=<AngleType; k=40.000, theteq=109.500>>
	<Angle <Atom CA [18]; In LEU 1>--<Atom C [33]; In LEU 1>--<Atom O [34]; In LEU 1>; type=<AngleType; k=80.000, theteq=120.400>>
	<Angle <Atom CA [18]; In LEU 1>--<Atom C [33]; In LEU 1>--<Atom N [35]; In TYR 2>; type=<AngleType; k=70.000, theteq=116.600>>
	<Angle <Atom N [16]; In LEU 1>--<Atom CA [18]; In LEU 1>--<Atom CB [20]; In LEU 1>; type=<AngleType; k=80.000, theteq=109.700>>
	<Angle <Atom N [16]; In LEU 1>--<Atom CA [18]; In LEU 1>--<Atom C [33]; In LEU 1>; type=<AngleType; k=63.000, theteq=110.100>>
	<Angle <Atom O [55]; In TYR 2>--<Atom C [54]; In TYR 2>--<Atom N [56]; In ILE 3>; type=<AngleType; k=80.000, theteq=122.900>>
	<Angle <Atom C [54]; In TYR 2>--<Atom N [56]; In ILE 3>--<Atom CA [58]; In ILE 3>; type=<AngleType; k=50.000, theteq=121.900>>
	...
	<Angle <Atom HA [295]; In SER 19>--<Atom CA [294]; In SER 19>--<Atom C [301]; In SER 19>; type=<AngleType; k=50.000, theteq=109.500>>
	<Angle <Atom CA [294]; In SER 19>--<Atom CB [296]; In SER 19>--<Atom HB2 [297]; In SER 19>; type=<AngleType; k=50.000, theteq=109.500>>
	<Angle <Atom CA [294]; In SER 19>--<Atom CB [296]; In SER 19>--<Atom HB3 [298]; In SER 19>; type=<AngleType; k=50.000, theteq=109.500>>
	<Angle <Atom H [293]; In SER 19>--<Atom N [292]; In SER 19>--<Atom CA [294]; In SER 19>; type=<AngleType; k=50.000, theteq=118.040>>
	<Angle <Atom N [292]; In SER 19>--<Atom CA [294]; In SER 19>--<Atom HA [295]; In SER 19>; type=<AngleType; k=50.000, theteq=109.500>>
])
TrackedList([
	<Dihedral; <Atom O [15]; In ASN 0>--<Atom C [14]; In ASN 0>--<Atom N [16]; In LEU 1>--<Atom CA [18]; In LEU 1>; type=<DihedralType; phi_k=2.500, per=2, phase=180.000,  scee=1.200, scnb=2.000>>
	<Dihedral; <Atom C [14]; In ASN 0>--<Atom N [16]; In LEU 1>--<Atom CA [18]; In LEU 1>--<Atom CB [20]; In LEU 1>; type=<DihedralType; phi_k=0.000, per=1, phase=0.000,  scee=1.200, scnb=2.000>>
	<Dihedral [ign]; <Atom C [14]; In ASN 0>--<Atom N [16]; In LEU 1>--<Atom CA [18]; In LEU 1>--<Atom CB [20]; In LEU 1>; type=<DihedralType; phi_k=0.000, per=2, phase=0.000,  scee=1.200, scnb=2.000>>
	<Dihedral [ign]; <Atom C [14]; In ASN 0>--<Atom N [16]; In LEU 1>--<Atom CA [18]; In LEU 1>--<Atom CB [20]; In LEU 1>; type=<DihedralType; phi_k=0.000, per=3, phase=0.000,  scee=1.200, scnb=2.000>>
	<Dihedral [ign]; <Atom C [14]; In ASN 0>--<Atom N [16]; In LEU 1>--<Atom CA [18]; In LEU 1>--<Atom CB [20]; In LEU 1>; type=<DihedralType; phi_k=0.000, per=4, phase=0.000,  scee=1.200, scnb=2.000>>
	<Dihedral; <Atom C [14]; In ASN 0>--<Atom N [16]; In LEU 1>--<Atom CA [18]; In LEU 1>--<Atom C [33]; In LEU 1>; type=<DihedralType; phi_k=0.000, per=1, phase=0.000,  scee=1.200, scnb=2.000>>
	<Dihedral [ign]; <Atom C [14]; In ASN 0>--<Atom N [16]; In LEU 1>--<Atom CA [18]; In LEU 1>--<Atom C [33]; In LEU 1>; type=<DihedralType; phi_k=0.000, per=2, phase=0.000,  scee=1.200, scnb=2.000>>
	<Dihedral [ign]; <Atom C [14]; In ASN 0>--<Atom N [16]; In LEU 1>--<Atom CA [18]; In LEU 1>--<Atom C [33]; In LEU 1>; type=<DihedralType; phi_k=0.000, per=3, phase=0.000,  scee=1.200, scnb=2.000>>
	<Dihedral [ign]; <Atom C [14]; In ASN 0>--<Atom N [16]; In LEU 1>--<Atom CA [18]; In LEU 1>--<Atom C [33]; In LEU 1>; type=<DihedralType; phi_k=0.000, per=4, phase=0.000,  scee=1.200, scnb=2.000>>
	<Dihedral; <Atom CG [9]; In ASN 0>--<Atom CB [6]; In ASN 0>--<Atom CA [4]; In ASN 0>--<Atom C [14]; In ASN 0>; type=<DihedralType; phi_k=1.046, per=1, phase=0.000,  scee=1.200, scnb=2.000>>
	<Dihedral [ign]; <Atom CG [9]; In ASN 0>--<Atom CB [6]; In ASN 0>--<Atom CA [4]; In ASN 0>--<Atom C [14]; In ASN 0>; type=<DihedralType; phi_k=0.303, per=2, phase=180.000,  scee=1.200, scnb=2.000>>
	<Dihedral [ign]; <Atom CG [9]; In ASN 0>--<Atom CB [6]; In ASN 0>--<Atom CA [4]; In ASN 0>--<Atom C [14]; In ASN 0>; type=<DihedralType; phi_k=0.033, per=3, phase=0.000,  scee=1.200, scnb=2.000>>
	<Dihedral [ign]; <Atom CG [9]; In ASN 0>--<Atom CB [6]; In ASN 0>--<Atom CA [4]; In ASN 0>--<Atom C [14]; In ASN 0>; type=<DihedralType; phi_k=0.107, per=4, phase=0.000,  scee=1.200, scnb=2.000>>
	<Dihedral; <Atom CB [6]; In ASN 0>--<Atom CA [4]; In ASN 0>--<Atom C [14]; In ASN 0>--<Atom O [15]; In ASN 0>; type=<DihedralType; phi_k=0.000, per=2, phase=0.000,  scee=1.200, scnb=2.000>>
	<Dihedral; <Atom CB [6]; In ASN 0>--<Atom CA [4]; In ASN 0>--<Atom C [14]; In ASN 0>--<Atom N [16]; In LEU 1>; type=<DihedralType; phi_k=0.200, per=1, phase=0.000,  scee=1.200, scnb=2.000>>
	<Dihedral [ign]; <Atom CB [6]; In ASN 0>--<Atom CA [4]; In ASN 0>--<Atom C [14]; In ASN 0>--<Atom N [16]; In LEU 1>; type=<DihedralType; phi_k=0.200, per=2, phase=0.000,  scee=1.200, scnb=2.000>>
	<Dihedral [ign]; <Atom CB [6]; In ASN 0>--<Atom CA [4]; In ASN 0>--<Atom C [14]; In ASN 0>--<Atom N [16]; In LEU 1>; type=<DihedralType; phi_k=0.400, per=3, phase=0.000,  scee=1.200, scnb=2.000>>
	<Dihedral [ign]; <Atom CB [6]; In ASN 0>--<Atom CA [4]; In ASN 0>--<Atom C [14]; In ASN 0>--<Atom N [16]; In LEU 1>; type=<DihedralType; phi_k=0.000, per=4, phase=0.000,  scee=1.200, scnb=2.000>>
	<Dihedral; <Atom CA [4]; In ASN 0>--<Atom CB [6]; In ASN 0>--<Atom CG [9]; In ASN 0>--<Atom OD1 [10]; In ASN 0>; type=<DihedralType; phi_k=0.000, per=1, phase=0.000,  scee=1.200, scnb=2.000>>
	<Dihedral; <Atom CA [4]; In ASN 0>--<Atom CB [6]; In ASN 0>--<Atom CG [9]; In ASN 0>--<Atom ND2 [11]; In ASN 0>; type=<DihedralType; phi_k=0.828, per=1, phase=180.000,  scee=1.200, scnb=2.000>>
	<Dihedral [ign]; <Atom CA [4]; In ASN 0>--<Atom CB [6]; In ASN 0>--<Atom CG [9]; In ASN 0>--<Atom ND2 [11]; In ASN 0>; type=<DihedralType; phi_k=0.485, per=2, phase=180.000,  scee=1.200, scnb=2.000>>
	<Dihedral [ign]; <Atom CA [4]; In ASN 0>--<Atom CB [6]; In ASN 0>--<Atom CG [9]; In ASN 0>--<Atom ND2 [11]; In ASN 0>; type=<DihedralType; phi_k=0.301, per=3, phase=180.000,  scee=1.200, scnb=2.000>>
	<Dihedral [ign]; <Atom CA [4]; In ASN 0>--<Atom CB [6]; In ASN 0>--<Atom CG [9]; In ASN 0>--<Atom ND2 [11]; In ASN 0>; type=<DihedralType; phi_k=0.008, per=4, phase=0.000,  scee=1.200, scnb=2.000>>
	<Dihedral; <Atom CA [4]; In ASN 0>--<Atom C [14]; In ASN 0>--<Atom N [16]; In LEU 1>--<Atom CA [18]; In LEU 1>; type=<DihedralType; phi_k=2.500, per=2, phase=180.000,  scee=1.200, scnb=2.000>>
	...
	<Dihedral [imp]; <Atom CZ [241]; In ARG 15>--<Atom HH21 [246]; In ARG 15>--<Atom NH2 [245]; In ARG 15>--<Atom HH22 [247]; In ARG 15>; type=<DihedralType; phi_k=1.000, per=2, phase=180.000,  scee=0.000, scnb=0.000>>
	<Dihedral [imp]; <Atom CZ [241]; In ARG 15>--<Atom HH11 [243]; In ARG 15>--<Atom NH1 [242]; In ARG 15>--<Atom HH12 [244]; In ARG 15>; type=<DihedralType; phi_k=1.000, per=2, phase=180.000,  scee=0.000, scnb=0.000>>
	<Dihedral [imp]; <Atom CD [236]; In ARG 15>--<Atom CZ [241]; In ARG 15>--<Atom NE [239]; In ARG 15>--<Atom HE [240]; In ARG 15>; type=<DihedralType; phi_k=1.000, per=2, phase=180.000,  scee=0.000, scnb=0.000>>
	<Dihedral [imp]; <Atom C [224]; In GLY 14>--<Atom H [227]; In ARG 15>--<Atom N [226]; In ARG 15>--<Atom CA [228]; In ARG 15>; type=<DihedralType; phi_k=1.100, per=2, phase=180.000,  scee=0.000, scnb=0.000>>
	<Dihedral [imp]; <Atom C [290]; In PRO 18>--<Atom CA [294]; In SER 19>--<Atom N [292]; In SER 19>--<Atom H [293]; In SER 19>; type=<DihedralType; phi_k=1.100, per=2, phase=180.000,  scee=0.000, scnb=0.000>>
])

6.3.4.2. Create the OpenMM System

This command creates an OpenMM System object from the information stored in trp_cage_gas. It contains multiple Force instances for the bonds, angles, periodic torsions, and nonbonded (electrostatic and van der Waals) interactions. It is in this function that we define the potential parameters we want to use. In this example, we have chosen the default values for each parameter except the ones specified. In particular:

  • nonbondedMethod=app.NoCutoff indicates we do not want to use a cutoff for nonbonded interactions.

  • constraints=app.HBonds indicates we want to constrain all bonds in which at least one atom is a Hydrogen (i.e., SHAKE or SETTLE for water). Other options are None (no constraints), app.AllBonds, or app.HAngles.

  • implicitSolvent=app.GBn2 indicates we want to use the second GBneck model described in Nguyen et al., J. Chem. Theory Comput., 2014 9(4) p.2020-2034. This a Generalized Born solvation model

  • implicitSolventSaltConc=0.1*u.liters/u.mole indicates we want to model a ca. 0.1 molar solution of monovalent ions using a Debye screening model.

# Create the OpenMM system
print('Creating OpenMM System')
system = trp_cage_gas.createSystem(nonbondedMethod=app.NoCutoff,
                               constraints=app.HBonds, implicitSolvent=app.GBn2,
                               implicitSolventSaltConc=0.1*u.moles/u.liter,
)

Creating OpenMM System

6.3.4.3. Create the integrator to do Langevin Dynamics

In this stage we specify an integrator. A common choices are for example LangevinIntegrator (as we’ve chosen here) to do simulations in the NVT ensemble or VerletIntegrator to do simulations in the NVE ensemble. In this example, we’ve chosen the Langevin integrator with a target temperature of 300 K, a friction coefficient of 1/ps and a time step of 0.5 fs.

# Create the integrator to do Langevin dynamics
integrator = mm.LangevinIntegrator(
                        0*u.kelvin,       # Temperature of heat bath
                        1.0/u.picoseconds,  # Friction coefficient
                        0.5*u.femtoseconds, # Time step
)

6.3.4.4. Create the Simulation object

This step creates a Simulation object that will be used to run the actual simulations. OpenMM picks the fastest platform for us if we do not specify a platform.

# Create the Simulation object
sim = app.Simulation(trp_cage_gas.topology, system, integrator)

This stage is very important. In this step, we set the particle positions stored in the trp_cage_gas object to our object. If you omit this step, you can get strange results or other errors like segmentation violations. These particle positions have been parsed from the input coordinate file.

# Set the particle positions
sim.context.setPositions(trp_cage_gas.positions)

6.3.4.5. Minimize the energy

This stage performs a basic energy minimization to relax particle positions. This particular invocation will perform at most 500 iterations.

# Minimize the energy
sim.minimizeEnergy(maxIterations=500)

state = sim.context.getState(getEnergy=True)
print(state.getPotentialEnergy())

-3089.34254023513 kJ/mol

6.3.4.6. Heating phase

After having obtained the minimized coordinates of our system we can begin slowly heating it. We define reporters that will “report” on the status of the simulation periodically throughout the heating phase. The first is a StateDataReporter which will print out a summary of energies and temperatures every 500 steps. This reporter directs the printout to standard output (the screen), sys.stdout or to a file name (if set).

The second reporter is a NetCDF trajectory reporter, which is written in the Amber NetCDF format.

6.3.4.7. CHANGE NUMBER OF STEPS HERE

Use the sim.step(nsteps) command to adjust the numer of steps such that 50 ps of heating is run in total. Check the integrator to see what time step we are using.

print('Heating phase')
sim.reporters.append(StateDataReporter('heating.csv', 500, step=True, potentialEnergy=True,  kineticEnergy=True, temperature=True)
 )
sim.reporters.append(
         NetCDFReporter('trp_cage_heating.nc', 500, crds=True)
)

for i in range(5):
    integrator.setTemperature(65*(i+1)*u.kelvin)
    sim.step(500)  ## Change this in order to run 50 ps of heating in 5 steps

state = sim.context.getState(getPositions=True)
minimized_and_heated=state.getPositions()

Heating phase

/home/lcbc/miniconda3/envs/jupyterbook/lib/python3.9/site-packages/parmed/amber/netcdffiles.py:409: UserWarning: Could not find netCDF4 module. Falling back on scipy implementation, which can significantly slow down simulations if used as a reporter
  warnings.warn('Could not find netCDF4 module. Falling back on '

#### Plotting the results from the trajectory run

import matplotlib.pyplot as plt
import pandas as pd
import numpy as np

df=pd.read_csv('heating.csv')

fig, ax = plt.subplots(1)

ax.plot(df['Time (ps)'],df['Potential Energy (kilocalorie/mole)'], label='Epot')
ax.plot(df['Time (ps)'],df['Kinetic Energy (kilocalorie/mole)'], label='Ekin')
ax.plot(df['Time (ps)'],df['Total Energy (kilocalorie/mole)'], label='Etot')
ax.set_xlabel('t (ps)')
ax.set_ylabel('Energy (kcal/mol)')
ax.legend()
plt.show()
../_images/TRP_Cage_OpenMM_colab_32_0.png 
fig, ax = plt.subplots(1)

ax.plot(df['Time (ps)'],df['Temperature (K)'], label='Temperature')
ax.set_xlabel('t (ps)')
ax.set_ylabel('Temp (K)')
ax.legend()
plt.show()
../_images/TRP_Cage_OpenMM_colab_33_0.png

6.3.5. Running dynamics

We first reset our Integrator to update the timestep, create the simulation object and then we run the MD of our heated system to see if our linear peptide folds up into its 3D form. We use again the NetCDF and StateDataReporter to retrieve information about our simulation.

6.3.5.1. CHANGE STEPS HERE

Below you again need to adjust the number of steps to run 20 ns of simulation in total. Again check with the integrator that you are using.

# Create the integrator to do Langevin dynamics
integrator = mm.LangevinIntegrator(
                        325*u.kelvin,       # Temperature of heat bath
                        1.0/u.picoseconds,  # Friction coefficient
                        2.0*u.femtoseconds, # Time step
)

# Create the Simulation object
sim = app.Simulation(trp_cage_gas.topology, system, integrator)

# Set the particle positions
sim.context.setPositions(minimized_and_heated)

# Set up the reporters to report energies and coordinates every 50000 steps
sim.reporters.append(
        StateDataReporter('production.csv', 50000, step=True, potentialEnergy=True,
                               kineticEnergy=True, temperature=True)
)
# Additionally print to the output of the notebook, so that we can track the simulation.
sim.reporters.append(
        StateDataReporter(sys.stdout, 50000, step=True, potentialEnergy=True,
                               kineticEnergy=True, temperature=True)
)
sim.reporters.append(
        NetCDFReporter('trp_cage_gb.nc', 50000, crds=True)
)

# Run dynamics
print('Running dynamics')
sim.step(1) # change this to run 20 ns of simulation

Running dynamics

# Plotting energy from production simulation

df=pd.read_csv('production.csv')

fig, ax = plt.subplots(1)

ax.plot(df['Time (ps)'],df['Potential Energy (kilocalorie/mole)'], label='Epot')
ax.plot(df['Time (ps)'],df['Kinetic Energy (kilocalorie/mole)'], label='Ekin')
ax.plot(df['Time (ps)'],df['Total Energy (kilocalorie/mole)'], label='Etot')
ax.set_xlabel('t (ps)')
ax.set_ylabel('Energy (kcal/mol)')
ax.legend()
plt.show()

---------------------------------------------------------------------------
EmptyDataError                            Traceback (most recent call last)
<ipython-input-21-2830dfba72c6> in <module>
      1 # Plotting energy from production simulation
      2 
----> 3 df=pd.read_csv('production.csv')
      4 
      5 fig, ax = plt.subplots(1)

~/miniconda3/envs/jupyterbook/lib/python3.9/site-packages/pandas/io/parsers.py in read_csv(filepath_or_buffer, sep, delimiter, header, names, index_col, usecols, squeeze, prefix, mangle_dupe_cols, dtype, engine, converters, true_values, false_values, skipinitialspace, skiprows, skipfooter, nrows, na_values, keep_default_na, na_filter, verbose, skip_blank_lines, parse_dates, infer_datetime_format, keep_date_col, date_parser, dayfirst, cache_dates, iterator, chunksize, compression, thousands, decimal, lineterminator, quotechar, quoting, doublequote, escapechar, comment, encoding, dialect, error_bad_lines, warn_bad_lines, delim_whitespace, low_memory, memory_map, float_precision, storage_options)
    603     kwds.update(kwds_defaults)
    604 
--> 605     return _read(filepath_or_buffer, kwds)
    606 
    607 

~/miniconda3/envs/jupyterbook/lib/python3.9/site-packages/pandas/io/parsers.py in _read(filepath_or_buffer, kwds)
    455 
    456     # Create the parser.
--> 457     parser = TextFileReader(filepath_or_buffer, **kwds)
    458 
    459     if chunksize or iterator:

~/miniconda3/envs/jupyterbook/lib/python3.9/site-packages/pandas/io/parsers.py in __init__(self, f, engine, **kwds)
    812             self.options["has_index_names"] = kwds["has_index_names"]
    813 
--> 814         self._engine = self._make_engine(self.engine)
    815 
    816     def close(self):

~/miniconda3/envs/jupyterbook/lib/python3.9/site-packages/pandas/io/parsers.py in _make_engine(self, engine)
   1043             )
   1044         # error: Too many arguments for "ParserBase"
-> 1045         return mapping[engine](self.f, **self.options)  # type: ignore[call-arg]
   1046 
   1047     def _failover_to_python(self):

~/miniconda3/envs/jupyterbook/lib/python3.9/site-packages/pandas/io/parsers.py in __init__(self, src, **kwds)
   1891 
   1892         try:
-> 1893             self._reader = parsers.TextReader(self.handles.handle, **kwds)
   1894         except Exception:
   1895             self.handles.close()

pandas/_libs/parsers.pyx in pandas._libs.parsers.TextReader.__cinit__()

EmptyDataError: No columns to parse from file

print('Average temperature')
print(f"{np.mean(df['Temperature (K)']):.3f} +- {np.std(df['Temperature (K)']):.3f} ")

Average temperature
83.891 +- 55.240 

fig, ax = plt.subplots(1)

ax.plot(df['Time (ps)'],df['Temperature (K)'], label='Temperature')
ax.set_xlabel('t (ps)')
ax.set_ylabel('Temp (K)')
ax.legend()
plt.show()
../_images/TRP_Cage_OpenMM_colab_39_0.svg

6.3.6. Save trajectory for later processing

Now you need to download the parameter file and trajectory file to view them later in VMD. Those can be found in the left sidebar in the files tab on the left (Bottom icon ). Your files will be deleted from Colab so it is important that you save them.

6.2. Protein Modelling 6.4. Hbond moving average and RMSD