----------------------
Water and Biomolecules
----------------------

.. contents::
.. section-numbering::

.. |angst|  unicode:: U+0212B .. ANGSTROM SIGN
.. |infin|  unicode:: U+0221E .. INFINITY
.. |simeq|  unicode:: U+02243 .. ASYMPTOTICALLY EQUAL TO
.. |sigma|  unicode:: U+003C3 .. GREEK SMALL LETTER SIGMA
.. |Delta|  unicode:: U+00394 .. GREEK CAPITAL LETTER DELTA
.. |mu|     unicode:: U+003BC .. GREEK SMALL LETTER MU
.. |beta|   unicode:: U+003B2 .. GREEK SMALL LETTER BETA
.. |pi|     unicode:: U+003C0 .. GREEK SMALL LETTER PI
.. |alpha|  unicode:: U+003B1 .. GREEK SMALL LETTER ALPHA
.. |nu|     unicode:: U+003BD .. GREEK SMALL LETTER NU
.. |deg|    unicode:: U+000B0 .. DEGREE SIGN


The N-methylacetamide (NMA) dimer
=================================

In this section, we look at different models for peptide chains
and for the hydrogen bonds between peptide chains. Before polypeptides
were computationally accessible, the NMA dimer has been a popular model
for a small peptide and for the peptide hydrogen bond. As the calculations
take too long time to be completed in the exercise and as you already know
how to calculate things, you are just given the output files.

* What do amino acids look like, and how is a peptide chain composed
  of different amino acids? What does a |beta|-sheet look like?  If
  you do not remember this, you can refresh this with the
  supplementary material.

* A relaxation of a pair of NMA dimers has been made. The results are
  in:

    NMA-dimer.nc_ NMA-dimer.txt_ `NMA-dimer.traj`_

  A relaxation of a single NMA molecule has been made. The results are in:

    NMA.nc_   NMA.txt_   `NMA.traj`_

.. _NMA-dimer.nc: http://dcwww.camp.dtu.dk/campos//download/10302/NMA-dimer.nc
.. _NMA-dimer.txt: http://dcwww.camp.dtu.dk/campos/download/10302/NMA-dimer.txt
.. _NMA.nc: http://dcwww.camp.dtu.dk/campos/download/10302/NMA.nc
.. _NMA.txt: http://dcwww.camp.dtu.dk/campos/download/10302/NMA.txt

.. _NMA.traj: http://dcwww.camp.dtu.dk/campos/download/10302/NMA.traj
.. _NMA-dimer.traj : http://dcwww.camp.dtu.dk/campos/download/10302/NMA-dimer.traj

  Copy these files to your area. Look at these calculations with
  ``plottrajectory`` and understand the structure of the NMA molecule
  (see also the article, which has been handed out).

* In order to measure bond lengths and angles in a simple way (instead
  of operating directly with the positions from dacapo), one can save
  the configuration in a ``.pdb`` file.  This is done with the
  following python commands:

  >>> from Dacapo import Dacapo
  >>> atoms = Dacapo.ReadAtoms('outfile.nc')
  >>> from ASE.IO.PDB import WritePDB
  >>> WritePDB('plot.pdb', atoms)

  and the configuration is saved in the file ``plot.pdb``. This file
  can be opened with Rasmol by typing ``rasmol plot.pdb`` at the
  command line. A Rasmol window pops up and a Rasmol command line
  appears. By clicking on the atoms, their number and type are printed
  on the command line.  By typing ``pick distance`` in the command
  line and clicking onto two atoms, their distance appears. One can
  measure the angles by typing ``pick angle`` and then clicking on
  the three atoms, which span the angle of interest.

* Do you think that the NMA dimer is a good model for a peptide
  chain or a |beta|-sheet? Which features of a peptide chain does it
  include and which are not modeled well?

* Calculate the interaction energy between the two NMA molecules.
  Which hydrogen bond is mainly responsible for the attractive
  interactions? Would you expect the interactions of a real
  |beta|-sheet to be stronger or weaker than the interaction between
  two NMA molecules? Calculate the interaction energy for PW91 (the
  functional used in the calculation) and also for PBE and RPBE (just
  take the non-selfconsistent energies). Try to comment on the
  differences.

* Two static calculations have been carried out, where each of the
  two NMA molecules forming the dimer are frozen: The structures are
  not relaxed, but they are kept in the position, in which the
  molecule forms the dimer. The results are contained in the files:

    NMA_freeze1.nc_ NMA_freeze1.txt_

  for the first NMA molecule and in:

    NMA_freeze2.nc_ NMA_freeze2.txt_

  for the second NMA molecule. Copy these files to your area.

.. _NMA_freeze1.nc: http://dcwww.camp.dtu.dk/campos/download/10302/NMA_freeze1.nc
.. _NMA_freeze1.txt: http://dcwww.camp.dtu.dk/campos/download/10302/NMA_freeze1.txt
.. _NMA_freeze2.nc: http://dcwww.camp.dtu.dk/campos/download/10302/NMA_freeze2.nc
.. _NMA_freeze2.txt: http://dcwww.camp.dtu.dk/campos/download/10302/NMA_freeze2.txt

* Now we are interested in, how the electronic density is changed,
  when two NMA molecules interact, i.e. how the two molecules are
  polarized. For this we plot the difference between the electron
  density of the NMA dimer and the sum of the densities of the two
  single NMA molecules alone.  This is done by the script:

    `diff_density_plot.py`_

  Try to understand, what this script does.  Then run it interactively
  in a python prompt.

* In the plot, blue and yellow areas show up. Can you understand from
  the script, which color stands for electron density depletion and
  which for electron density addition?

* From the plot one can see which hydrogen bond is mainly responsible
  for the interaction.  Which one is it?

..
  * The unit of the electron density output from dacapo is a bit
    awkward. Summing the electron density up over the whole grid gives
    $N_{el}\cdot 80 \cdot 110 \cdot 80$, where N$_{el}$ are the number
    of electrons in the system 80x110x80 are the dimensions of the hard
    grid used for the density (can be found in the ``*.txt`` output
    file.  Thus, to get the value in SI units from the dacapo output,
    one has to divide by $\cdot 80 \cdot 110 \cdot 80$.

.. _diff_density_plot.py : attachment : diff_density_plot.py

Parallel and antiparallel two-stranded beta-sheets
====================================================

A realistic extended model of |beta|-sheets is only marginally
accessible to DFT or other quantum-mechanical calculations. Therefore,
calculations are often done with force fields, which are optimized to
reproduce experimental |beta|-sheet geometries.  Two |beta|-sheet
structures, parallel and antiparallel have been calculated and are in
the files:

  `gly-gly10_par.pdb`_ `gly-gly10_anti.pdb`_

for chains consisting of glycine and:

  `ala-ala10_par.pdb`_ `ala-ala10_anti.pdb`_

.. _gly-gly10_par.pdb : attachment : gly-gly10_par.pdb
.. _gly-gly10_anti.pdb : attachment : gly-gly10_anti.pdb 
.. _ala-ala10_par.pdb : attachment : ala-ala10_par.pdb 
.. _ala-ala10_anti.pdb : attachment : ala-ala10_anti.pdb 


for chains consisting of alanine. Look at the structures and try to
understand them. How have the peptide chains been terminated?

* The most important hydrogen bond in |beta|-sheets is the N-H...O=C
  hydrogen bond. How well do these bonds fit together in the parallel
  and antiparallel |beta|-sheets, where do they fit better?  Which
  type of |beta|-sheet does the NMA dimer resemble?

* Measure the bond lengths of the bonds in the peptide chain (C-N
  and C-C bonds) and the C=O bond.  From other organic compounds the
  following approximate bond lengths are known

  - C=O double bond: 1.20 Å
  - C-C single bond: 1.54 Å
  - C-N single bond: 1.47 Å
  - C=N double bond: 1.27 Å

Which bonds in the peptide chain have single bond character and which
have double bond character?  What does that mean for the flexibility
of the peptide chain?

* Measure the bond geometries (lengths and angles) of the
  N-H...O=C hydrogen bond both in the NMA dimer and in
  the alanine |beta|-sheets. How well do they agree?

* What effect does water as a solvent have on the stability of
  |beta|-sheets. Which type of amino acids would you expect to find
  in a |beta|-sheet?


The water dimer
===============

The file 2H2O.nc_ is the output from a Dacapo calculation of a
water dimer.  Study, and understand the script `2H2O.py`_ that
was used for the calculation (don't repeat the calculation).

.. _2H2O.nc: http://dcwww.camp.dtu.dk/campos/download/10302/2H2O.nc

* The dimer structure has already been relaxed.  What is the length
  of the hydrogen bond?

* Make a density-difference plot like you did for the NMA
  molecule.  Now you will have to do the two frozen water molecule
  calculations yourself.  The first one can be done like this::

    atoms = Dacapo.ReadAtoms('2H2O.nc') # Restart from the old calculation
    del atoms[0:3]                      # Delete the first three atoms
    calc = atoms.GetCalculator()        # Get a hold on the calculator
    calc.SetNumberOfBands(6)            # Reduce the number of bands
    calc.SetNetCDFFile('H2Oa.nc')       # Write output to 'H2Oa.nc'
    print atoms.GetPotentialEnergy()    # Go!

* What is the binding energy of the dimer?  Compare with the
  binding energy of the NMA dimer.


.. _2H2O.py : attachment : 2H2O.py