from chemsh import *

# this example shows how to port the P450 system prepared by Py-ChemShell
# to Tcl-ChemShell via PSF
# for Tcl-ChemShell calculation, see:
# appendix_1_porting_py-_to_tcl-chemshell/

# you.lu@stfc.ukri.org
# Oct 2021

# read in the cut-out shell-shape structure
p450_w_water_shell = Fragment(coords='../2_qmmm/p450_w_water_shell.pdb', title='p450')

# create a NAMD instance, optionally with user files (but let's here use files
# same as for the solvation and QM/MM setup)
# (parameters file is not needed by PSFGEN)
my_namd = NAMD(top = [ '../1_solvation/camphor.rtf',
#                       '../top_all22_heme.inp'
                     ],
               frag = p450_w_water_shell,
              )

# call PSFGEN to generate a PSF file
my_namd.psfgen(return_frag=False, guess_coord=True)


# * we need to concatenate the forcefield files since Tcl-ChemShell takes
#   only one parameter file
# * par_all36m_prot.prm, toppar_water_ions.str, and haems-ligands.prm are
#   Py-ChemShell's built-in data files; camphor.prm is user provided
NAMD.mergeParameters([ 'par_all36m_prot.prm',
                       'toppar_water_ions.str',
                       'haems-ligands.prm',
                       '../1_solvation/camphor.prm' ],
                     save_to='all36m_prot_water_ions_haems_ligands_camphor.prm', # save to filename (default: _namd.prm)
                     ff='charmm',                                                # forcefield type (default: 'charmm')
                     reverse=False,                                              # reverse the files' order to make former ones overriding the latter (default: False)
                     fill_with_zeros=False,                                      # write dummy zeros for missing coefficients (default: False)
                                                                                 # do not use it for Tcl-ChemShell!
                     )
