This is an old revision of the document!
Potential energy surface of glyala dipeptide
Ramachandran plot
Glyala is one of the simplest molecules that exhibits some important features common to larger biomolecules. In particular, it has more than one long-lived conformation, which we will identify in this exercise by mapping out its potential energy surface.
The conformations of alanine dipeptide are characterized by the dihedral angles of the backbone. Below, we color carbons in green, hydrogens in white, oxygen in red and nitrogen in blue, i.e. the torsional angle $\phi$ is N-C-C-N , while $\psi$ is C-N-C-C along the backbone.
tar -xvf glyala-epot.tar.gz
glyala.pdb with VMD and determine the atomic indices of the atoms defining the dihedral angles.
With this knowledge at hand, we will fix the dihedral angles and perform geometry optimization for all remaining degrees of freedom.
- The atomic indices defining the dihedral indices in the input file
geo.inare missing. ReplaceI1toI4by the atomic indices determined previously. Note: While VMD starts counting atoms from 0, CP2K starts counting from 1, i.e. the VMD indices need to be increased by 1. - Use
perform-gopt.shto perform the grid of geometry optimizations. - Use gnuplot to plot the potential energy surface (we have provided a script
epot.gp). Which are the two most favoured conformations?$ gnuplot
gnuplot >> load "epot.gp"
Water
We have prepared a CP2K input file md.in for running a MD simulation of liquid water using the force field from the first exercise (parametrized by Praprotnik et al.).
- Check that the MD is energy conserving and well-behaved.
Repeat the MD using initial temperatures 200 and 400 K. In order not to overwrite any of your previous files, it is advisable to run the new simulations in different folders.
- What are the final average temperatures in each of the simulations?
- Why are they different from the initial ones?
- The initial atomic configuration stems from an equilibration run. At which temperature was the system (approximately) equilibrated?
Next we are going to analyze the trajectories in order to calculate the radial distribution function (rdf, $g(r)$) as a function of temperature.
VMD comes with an extension for exactly this purpose: In the VMD Main window open “Extensions → Analysis” click on “Radial Pair Distribution function $g(r)$”. In the appearing window use “Utilities → Set unit cell dimensions” to let VMD know the simulation box you used. After that use Selection 1 and 2 to define the atomic types that you want to calculate the rdf for, for example “element H”.
- Plot $g_{O-O}(r)$ at 200, 300 and 400 K into the same graph.
- What are the differences in the height of the first peak?
- What does this say about the structure of the liquid and is this expected? (2P)
- Compare to experimental data
goo.ALStaken at 300 K.
Then we will calculate diffusion coefficient. The diffusion coefficient is a proportionality constant between the molar flux due to molecular diffusion and the gradient in the concentration of the species (or the driving force for diffusion), which is defined by: $6D=\lim_{t\to\infty} \ \frac{\delta <r^2(t)>}{\delta t}$
Glyala in water
Now, we will move to a more realistic system - Glyala in water. We will preformed a MD of glyala in water and save the trajectory.
The initial geometry provided in the PDB file is a glyala molecule solvated by 73 water molecules. The geometry is not equilibrated. You need first to equilibrate the system at 300K. When the system is equilibrated, you need to analysis the result.
- Perform the molecular dynamics simulation using NVT ensemble at 300K.
- Re-run the calculation using NVT ensemble with different TIMECON (500, 2000 fs) in the &THERMOSTAT section, and plot the total energy, temperature against time. Explain what you observe.
- Determine from which step the system is equilibrated, plot the calculated properties and explain why.
- Compute the O-O radial distribution function for water with acceptable statistics using 20 ps (after equilibration) of simulated time.
In last exercise, one already knew how to calculate the RDF for the Argon system. In TASK3, you need to calculate the RDF only for water instead of whole system. Since the glyala contain two oxygen atoms, it is not reasonable to include the oxygen atoms in glyala molecule if we are only interested in O-O RDF for water. Using VMD, the O-O RDF for the water can be easily calculated. In the
Selection 1, Selection 2
, one need to specify
element O and not same residue as element C
The frames should start from the beginning of production run.