# Cut a cluster from bulk shell model MgO, with coordination-scaled charges

from numpy import array, arange
from chemsh import *
from chemsh.utils import testutils

# Note as we're reading coordinates directly in rather than using 
# a c_create-like routine there is no testing for equivalent sites
# so haven't shifted any atoms or shells around as the Tcl test did
mncpco3 = Fragment(coords='bulk_molecular.c', connect_mode='covalent')

cluster = mncpco3.construct_cluster(origin_atom=0,
                                    radius_cluster=40.0,
                                    crystal_type='molecular')

cluster.save('bulk_molecular_saved.pun', 'pun')

# Tcl-ChemShell reference
ref = Fragment(coords='tcl-chemshell/ref/ref_bulk_molecular.c', connect_mode='covalent')

# Test atoms
testutils.validate(cluster.natoms , ref.natoms)
testutils.validate(cluster.znums  , ref.znums, summarise='inf')
testutils.validate(cluster.coords , ref.coords , tol=1.0e-10)

# relaxed tolerance due to numerical differences in bq least squares fit
#testutils.validate(cluster.charges, ref.charges, tol=1.0e-6)


# Test shells
testutils.validate(cluster.shells.coreatoms, ref.shells.coreatoms)
testutils.validate(cluster.shells.coords   , ref.shells.coords , tol=1.0e-10)
testutils.validate(cluster.shells.charges  , ref.shells.charges, tol=1.0e-10)

# Test connectivity
testutils.validate(cluster.conn, ref.conn, summarise=100)
