Assignment 3: Stat Mech and Dynamics

DUE: 11:59 PM THURSDAY, SEPTEMBER 24

For this assignment you will run molecular dynamics simulations on the same amino acid structures you used for assignment 2. Specifically, you are to run an MD simulation for a certain number of time steps, and at each timestep you need to record the same dihedral angle of the amino acids that you analyzed in assignment 2. You will plot the histogram of dihedral angles observed during the simulation.

In addition to the histogram of dihedral angles observed during MD simulation, you will also plot the Boltzmann probability density as a function of the dihedral angle on the same figure. You wrote code to compute this probability density in assignment 2. You will need to reuse this code from assignment 2 in order to generate the figures for assignment 3.

Use your script for assignment 2 as a starting point for assignment 3. Below are some directions on what you should add/remove from the assignment 2 script and some tips on how to run the MD simulation. If you did not get a working assignment 2 script, seek assistance immediately from course staff.

Replace the command line arguments of your script with the following code block:

parser = argparse.ArgumentParser(description='CompStruct Assignment 3')
parser.add_argument('pdb',help='input PDB file')
parser.add_argument('--aindex',help='index of first dihedral atom',type=int,default=1)
parser.add_argument('--bindex',help='index of second dihedral atom',type=int,default=2)
parser.add_argument('--integrator',help='integrator to use',choices=['verlet','langevin'],default='langevin')
parser.add_argument('--temp',help='temperature to simulate at',type=int,default=300)
parser.add_argument('--steps',help='number of simulation steps (1fs)',type=int,default=10000)
parser.add_argument('--output',help='output filename for graph',default='out.png')
args = parser.parse_args()

# setup openmm system using amber 14 forcefield which defines the potential energy function
pdb = PDBFile(args.pdb)
ff = ForceField('amber14-all.xml')
system = ff.createSystem(pdb.topology,ignoreExternalBonds=True)

# set integrator
if args.integrator == 'verlet':
    integrator = VerletIntegrator(1*femtosecond)
else:
    integrator = LangevinIntegrator(args.temp*kelvin, 1/picosecond, 1*femtosecond)
    integrator.setRandomNumberSeed(42)

Note the importance of setting the random seed to get reproducible results.

You will need to keep the portion of your script that computes the Boltzmann probability of each dihedral angle. However, this time, you do not need to compute the probabilities for multiple hardcoded temperatures. Rather, you only need to compute the probabilities for the temperature specified by the command line arguments (available through the variable args.temp). Once you have the code that computes the Boltzmann probability density over discrete dihedral angles, you can delete all of the subsequent code in the script. Specifically, the code that computes the probabilities of states A and B, ⟨E⟩, S, and F, and ΔG, and the code that generates plots. As in assignment 2, save the dihedral angles to a variable named angles and the corresponding probabilities in a variable named probs.

Write code in your script to run a MD simulation on the provided pdb for args.steps steps using the specified integrator and default simulation settings. If you do the dihedral angle scan before the MD simulation, make sure you set the simulation positions to the pdb positions before starting the MD simulation. Refer to the "Running Simulations" example from the OpenMM documentation to figure out how exactly to run an MD simulation.

NOTE: Do not do an energy minimization. Do not do a system equilibration. Do not solvate the system.

You need to save the dihedral angle (measured as in the last assignment) for each step of the simulation. You can do this by simulating a single step at a time and computing the dihedral angle after each step, or by implementing a custom reporter class. Save the dihedral angle values (in radians) in a list named dihedrals.

Conclude your script with the following code that will generate the necessary plots:

#convert to degrees
dihedrals_deg = (np.rad2deg(dihedrals)+360)%360
cnts,bins = np.histogram(dihedrals_deg,range(361))

#this is what the autograder will look at (with a small time step)
print(cnts)

#make plot
plt.hist(dihedrals_deg,bins=range(361),density=True)
plt.plot(angles,probs,label=r'$\frac{1}{\hat{Z}}e^{\frac{-U}{k_BT}}$')
plt.xlim(0,360)
plt.legend(fontsize=16)
plt.xlabel('Dihedral Angle (Degrees)')
plt.ylabel('Frequency/Probability')
plt.title(os.path.split(args.pdb)[-1][:-4].upper())
plt.savefig(args.output,bbox_inches='tight')

Submission

Part 1: Submit your final script to the autograder (you may submit as many times as you need).

Note that due to differences between systems (e.g. the random number generator, floating point accelerator availability), you will likely print out different numbers on your local machine than on the autograder, but you should still match output exactly when run on the autograder as long as you set the random seed as described above.

In Part 2 you will be asked to generate and upload graphs of 1 million steps of simulation. This will take much less time if run on a GPU on the cluster. Getting familiar with using SLURM now will make your life easier in the next assignment.

Part 2: Answer the questions on GradeScope (in case you are wondering, GradeScope doesn't allow combining programming assignments and question-based assignments into one assignment, hence the two parts).