# Single-point QM/MM calculation on the P450 system

# Created by YL Jun 2022

# Initialise ChemShell
from chemsh import *

# Read in protein system with solvent shell in PDB format (without charges)
p450_w_water_shell = Fragment(coords='p450_w_water_shell.pdb')

# Define the QM region by building up atom lists of the haem, oxygen and 
# cysteinate ligand
# NB: Remember that Py-ChemShell atom numbering starts from zero!

# Haem (without side chains) 
list_indices_haem_wo_sidechain = list(range(6384, 6413))
# Oxygen: atoms 6355, 6356
list_indices_o2 = p450_w_water_shell.select(residues='O2').tolist()
# CYS347 (truncated as SCH3): atoms 5445, 5446, 5447, 5448
list_indices_cys347 = p450_w_water_shell.select(residues='CYS', resIDs=347, 
        truncate_by='alpha-beta', visualise='vmd').tolist()

# Concatenate the lists to form the QM region
indices_qm_region = list_indices_haem_wo_sidechain + list_indices_o2 + list_indices_cys347

# Save for viewing
frag_qm_region = p450_w_water_shell.getSelected(indices_qm_region)
frag_qm_region.save('qm_region.pdb')

# Define QM method
# The QM total charge will be deferred by setting charge='automatic' and determined later
# by DL_FIELD as the sum of the formal charges
my_qm = GAMESS_UK(method='dft', functional='b3lyp', atom_names='symbol',
                  basis='hcno_6-31gp_fe_lanl2dz_gamess-uk.bas', ecp='fe_lanl2dz_gamess-uk.ecp',
                  charge='automatic', mult=3, scftype='uhf', maxiter=2000, restart=False)

# Define MM forcefield using DL_FIELD

# First, convert user-defined CHARMM files for camphor to DL_FIELD format
udff = DL_FIELD.charmm2udff(rtf='../1_solvation/camphor.rtf', 
        prm='../1_solvation/camphor.prm', version='CHARMM36_prot')

# Now use DL_FIELD to define the forcefield using built-in CHARMM36_prot
# NB: cutoff distance defined in a.u.
my_ff = DL_FIELD(ff=['CHARMM36_prot', udff],
                 override_builtin=False, original_tip3p=False)

# Define the MM method (DL_POLY) using the forcefield defined by DL_FIELD
# Default: dl_field_charges=True
my_mm = DL_POLY(ff=my_ff, dl_field_charges=True)

# Define the QM/MM calculation
my_qmmm = QMMM(frag=p450_w_water_shell, qm=my_qm, mm=my_mm, 
               qm_region=indices_qm_region)

# Run a single-point QM/MM energy evaluation
my_sp = SP(theory=my_qmmm)
my_sp.run(dryrun=False)
