# Cut solvent shell around protein system from MD snapshot
#
# Created by YL Jul 2021
# Revised by YL, KS, TWK Oct-Dec 2021

# Initialise ChemShell
from chemsh import *

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

# If you have run the full (long) solvation workflow example, you can instead
# use your own snapshot by commenting the line above and uncommenting:
#p450 = Fragment(coords='../1_solvation/_snapshots/_snapshot_19.pqr')

# Select a solvent shell around the protein to reduce the system's size

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

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

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

# 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')
