DUE: 11:59 PM TUESDAY, SEPTEMBER 15
In this assignment you will be using OpenMM to calculate the energies, entropies, and free energies of the C-C$\alpha$ dihedral bond of selected amino acid residues.
For the purposes of our analysis, we will assume the bond has 360 discrete states, one for each degree of rotation. We will use the Amber14 molecular dynamics force field to calculate potential energies of molecular conformers.
For simplicity, we are going to ignore kinetic energy (that is, we are only considering the configurational partition function).
Since we are using a molecular dynamics force field calibrated to proteins, which are chains of amino acids connected to one another, we will be analyzing individual residues, not an actual independent amino acid. Specifically, the chemically astute will notice there is a hydrogen missing from the nitrogen. For our pedagogical purposes, we aren't going to worry about this, but be aware if you find yourself needing to simulate a single amino acid (e.g. glycine in a glycine receptor), you will need an additional set of parameters to what is in a typical force field (which assumes each amino acid is bonded to at least one other).
The amino acid structures we will be evaluating are here:
We provide starter code that takes care of rotating atoms around the desired dihedral. Note that the code requires the alpha carbon be at the origin, and this is true for all our provided files.
Note that OpenMM defines the constants BOLTZMANN_CONSTANT_kB, AVOGADRO_CONSTANT_NA, and MOLAR_GAS_CONSTANT_R with appropriate units. Quantity objects can be converted to different units with in_units_of().
Read through this code — there are additional instructions embedded within the code.
You need to modify this code in the designated place to produce the desired output. For example:
./assign2.py ala.pdb --out out.png
T=100 K, Z=1.14e-51, <U>=101.40 kJ/mol, S=38.69 J/(K mol), F=97.53 kJ/mol
T=300 K, Z=3.08e-16, <U>=102.50 kJ/mol, S=44.69 J/(K mol), F=89.09 kJ/mol
T=500 K, Z=4.5e-09, <U>=103.35 kJ/mol, S=46.90 J/(K mol), F=79.90 kJ/mol
T=1000 K, Z=0.00119, <U>=104.33 kJ/mol, S=48.34 J/(K mol), F=55.98 kJ/mol
At 300 K:
Probability of state A is 0.4725
Probability of state B is 0.5275
State A: <U>=102.70 kJ/mol S=39.14 J/(K mol) F=90.96 kJ/mol
State B: <U>=102.31 kJ/mol S=38.75 J/(K mol) F=90.69 kJ/mol
ΔG is 0.275 kJ/mol
Your code will also produce a fancy graph that you will upload in the next part.
Submission
You will submit your single .py script to the GradeScope autograder (Part 1). You must match the reference output exactly (do not change the provided print statements). You may submit as many times as you want, but only your last submission will count.
You also need to complete Part 2, which calls for written answers.
Your code is evaluated on a Ubuntu 22.04 base image with OpenMM installed using conda:
wget https://repo.anaconda.com/miniconda/Miniconda3-py312_26.7.1-1-Linux-x86_64.sh
chmod +x Miniconda3-py312_26.7.1-1-Linux-x86_64.sh
./Miniconda3-py312_26.7.1-1-Linux-x86_64.sh -b
conda init bash
conda install -c conda-forge openmm matplotlib numpy scipy -y