'''
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 *

# construct a ChemShell fragment from a file describing the MgO unit cell
mgo_cell = Fragment(coords   = 'mgo_shells_unitcell.cjson',
                    connect_mode = 'ionic')                    # connectivity mode

# check the lattice parameters
print("MgO lattice constants:", mgo_cell.cell.consts)
print("MgO lattice vectors:", mgo_cell.cell.vectors)

# save to an XYZ file for convenience of visualisation
# NB: the saved XYZ file doesn't contain any cell information
mgo_cell.save('mgo_cell_nonperiodic.xyz')

# radii parameters for cluster cutting and QM/MM partitioning:
#
# * origin    |<- cutoff_boundary ->|                         |<- bq_margin ->|
#
#  QM         | Interface           | Active MM  | Frozen MM  | Point charges |
#  (region 1) | (region 2)          | (region 3) | (region 4) | (region 5)    |
#
# |<-------------- radius_active --------------->|
#
# |<---------------------- radius_cluster ------------------->|
#
# (qm_radius, and cutoff_boundary will be defined in a later stage)
#
# the returned value is a fragment instance of the production cluster
# NB: the internal units of ChemShell are always a.u. (i.e., bohr)
mgo_cluster = mgo_cell.construct_cluster(origin=1,            # can be a set of fractional coodinates, e.g., [0.0,0.0,0.0]
                                         radius_active=15.0,  # 5.0 a.u., only used by BQs fitting in this stage
                                         radius_cluster=20.0, # 20.0 Angstrom
                                         bq_margin=5.0,       # BQs thickness
                                         bq_density=3,        # BQs density
                                         adjust_charge='coordination_scaled')

# save the production cluster to a new file
mgo_cluster.save('mgo_shells_cluster.pun', 'pun')

# XYZ for visualisation
mgo_cluster.save('mgo_shells_cluster.xyz', 'xyz')


