# this is always required when you launch ChemShell using the executable script "chemsh"
from chemsh import *

# use the list of charges and names (types) defined in "p450/_p450_charges.py" and "p450/_p450_names.py"
# ChemShell hint: the user can import other Python files under the same working directory just as Python modules
from _p450_charges import charges
from _p450_names   import names

# read in Tcl-ChemShell's original turorial P450 PDB and create a fragment instance
# NB: at the same time the loaded names will overwrite those in the PDB
p450 = Fragment(coords='p450.pdb', charges=charges, names=names)

# define the QM atom indices: the embeded elements will be flattened by ChemShell later
# (Python hint: https://docs.python.org/3/library/stdtypes.html#range)
indices_qm_region = [ range(4575, 4579), range(5147, 5177) ]

# the above is just equivalent to:
indices_qm_region = [4575, 4576, 4577, 4578, 5147, 5148, 5149, 5150, 5151, 5152, 5153, 5154, 5155, 5156, 5157, 5158, 5159, 5160, 5161, 5162, 5163, 5164, 5165, 5166, 5167, 5168, 5169, 5170, 5171, 5172, 5173, 5174, 5175, 5176]

# the reader can visualise the selected atoms by printing out a string of the indices:
vmd_select = 'index ' + ' or index '.join([ str(s) for s in indices_qm_region ])
print("\n To visualise with VMD:\n", vmd_select)

# write a new PDB of the QM region (just for reference)
frag_qm_region = p450.getSelected(indices_qm_region)
frag_qm_region.save('qm_region.pdb')

# get the total integer charge of the selected QM region
qm_charge = round(frag_qm_region.charges.sum())
print("\n QM charge =", qm_charge)

# create an ORCA instance
# NB: because the ORCA distro is implicitly driven by OpenMPI, it requires ChemShell to be compiled
#     as serial or we have to run ORCA in serial by setting nprocs=1
#     If running in parallel, the full path to ORCA must be set in the path option
#     (this is a requirement by ORCA rather than ChemShell)
my_qm = ORCA(method     = 'dft',          # DFT calculation
             scftype    = 'uhf',          # UHF Hamiltonian because the system is open-shell
             functional = 'b3lyp',        # DFT functional
             basis      = 'p450_tcs.bas', # use the predefined basis set file
             charge     =  qm_charge,     # the total integer charge of the QM region obtained above
             maxiter    =  100,           # maximum allowed number of SCF iterations (NB: we may change this keyword in the future)
             mult       =  2,             # spin multiplicity is 2 because there is a single unpaird electron
             path       = 'orca',         # (full) path to the ORCA executable
# YL 05/03/2021: the ORCA interface needs refactoring as the default guess is Hueckel which results in deviated energy;
#                we just switch it off for now
             guess      = '',
             nprocs     =  1,            # number of MPI processes for ORCA (NB: this option is exclusive for ORCA)
            )

# create a DL_POLY instance, using the predefined forcefield file (FIELD) generated by Tcl-ChemShell
my_mm = DL_POLY(ff='p450_tcs.ff')

# create a QM/MM instance: we only assign the fragment to the QM/MM instance, not elsewhere
my_qmmm = QMMM(frag      = p450,
               qm        = my_qm,
               mm        = my_mm,
               qm_region = indices_qm_region,
              )

# create a single-point energy evaluation instance and run
my_sp = SP(theory=my_qmmm)

# if dryrun=True the job will not be really performed but only with input files generated
my_sp.run(dryrun=False)


# energies (a.u.) calculated by Tcl-ChemShell with ORCA 4.2.1 and DL_POLY Classic:
e_tcs_orca = -1624.160155
e_tcs_dlp  = -22.822932
e_tcs_qmmm = -1646.983088




