'''
Disclaimer: This sample is only for demonstration purpose
and does not have to make chemical/physical sense.

Nov 2019
you.lu@ukri.stfc.org
'''

from chemsh import *

# read in the cluster from the previous stage
mgo_cluster = Fragment(coords='../01_cluster_cutting/mgo_shells_cluster.pun', connect_mode='ionic')

# QM atoms can be selected by complicated methods, for example by radius of 4.0 a.u.
qm_list = mgo_cluster.select(radius=4.0, centre=1)
# here we keep it simple as an explicit list for demonstration
qm_list = [0,1]

# NB: if "implicit" interface (default) is chosen here, interface atoms (region "2" by default,
#     "Mg2" in this sample's output fragment) have to be contained in the QM region (see keyword
#     "qm_region" in the later stage of calculations) and defined in ECPs
#     "explicit" does not require QM region inclusion but have to be defined as BQ EPCs
# returned value is a ChemShell fragment with partitioned labels ("1","2","3","4",&"5" by default)
partitioned = mgo_cluster.partition(# use the QM region selection
                                    qm_region=qm_list,
                                    # interface layer for ECPs: 7.5 a.u. thickness
                                    cutoff_boundary=7.5,
                                    # 15.0 a.u. radius from origin
                                    radius_active=15.0,
                                    # index of atom or X|Y|Z coordinates such as [0.0,0.0,0.0]
                                    origin=1,
                                    # optional: use only magnesium as QM termination (outermost atoms of region "1")
                                    #qm_termination='Mg',
                                    # can be overloaded by any ChemShell-supporting file type
                                    qm_save_to='mgo_partitioned_qm_region.xyz',
                                    # implicit interface (default)
                                    qmmm_interface='implicit',
                                    # oxygen (or normally other anions) must be excluded since there are no proper ECPs
                                    interface_exclude=["O"])

# a string "OxxMgyy" of QM region (region "1", excluding interface)
qm_formula = partitioned.getSelected(partitioned.select(regions=1)).formula

# this step is not necessary because ChemShell writes out a new punch file anyway
# save to a new file for later use (better clarify the "implicit" status!)
partitioned.save('mgo_shells_qmmm_implicit.pun')


