.. _p450_shell:

********************************
Preparing the non-periodic model
********************************

.. note:: The scripts for this step of the tutorial can be found in 
  ``p450_enzyme/2_qmmm``.


In this stage we prepare a model suitable for QM/MM calculations based on 
an MD snapshot generated from the solvation workflow in the previous stage.


Cutting a solvent shell
=======================

Unlike the initial MD runs in the solvation workflow, the QM/MM
calculation is non-periodic and does not require a cubic water box.
We can reduce the system size and therefore also the computational
expense of the calculation by cutting out a solvent shell around the 
protein.

The script to cut the solvent shell is ``cut_solvent_shell.py``.
Firstly, we read in a pre-prepared snapshot in PQR format (to include 
charges)::

 # Read in pre-prepared snapshot PQR file (including charges)
 p450 = Fragment(coords='p450_full_snapshot_sample.pqr')

Alternatively, if you have successfully run the long solvation workflow
example in the previous stage, you can read in a snapshot from that
instead, e.g.::

 p450 = Fragment(coords='../1_solvation/_snapshots/_snapshot_19.pqr')

.. note :: A PDB file may also be used instead of a PQR file, if you
 also provide a list of charges, i.e.::

   p450 = Fragment(coords='snapshot.pdb', charges=[1.00,...])

Once the Fragment has been read in, we identify a list of non-solvent
atoms using an atom selection command::

  # Identify non-solvent atoms to cut shell around
  indices_non_solvent = p450.selectByCommonNames('non-solvent')

Next, we can cut the shell itself, using a special selection
command defined for this purpose::

  # Cut the shell
  indices_p450_w_water_shell = p450.selectByShell(around=indices_non_solvent,
        padding=7.0, unit='angstrom', boundary='incl', convex=True)

The ``selectByShell`` method is controlled by the following options:

* ``around`` specifies the system to cut around, defined by the previous 
  atom selection.
* ``padding`` defines the size of the shell to cut, in units defined by 
  ``unit``.
* ``boundary`` - if ``'incl'`` (default) will keep whole residues in
  the shell, if any part of them is within the padding distance.
  Alternatively, ``excl`` will remove whole residues unless all of the 
  atoms are within the padding distance.
* ``convex`` defines the selection algorithm. If ``True`` (default), the
  convex hull algorithm will be used which is the most efficient option, 
  but results in a slightly different (smoother) shell than an explicit 
  distance criterion (``convex=False``).

Once the shell selection has been made, a new Fragment object is created
using the selection::

 # Create a new Fragment object from the shell
 p450_w_water_shell = p450.getSelected(indices_p450_w_water_shell)

Finally, the new fragment needs to be saved to disk for use in the
following QM/MM calculation::

 # Save the new Fragment for use in the QM/MM calculation
 p450_w_water_shell.save('p450_w_water_shell.pqr')
 p450_w_water_shell.save('p450_w_water_shell.pdb')

Before moving on to the next step, open ``p450_w_water_shell.pdb``
in your visualiser and see how it differs from the system used 
in the solvation workflow.