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

from numpy import array
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
mgo = Fragment(coords='bulk_valence_matrix.c', connect_mode='ionic')

# has to be a NumPy array rather than Python list
# it must be the transpose of Tcl-ChemShell input matrix because NumPy is row-major!
valence_matrix = array([[0.0,-1/3],[1/3,0.0]])

cluster = mgo.construct_cluster(origin_atom=1,
                                radius_cluster=20.0,
                                radius_active=10.0,
                                bq_margin=7.0,
                                bq_density=3,
                                coordinated='cations',
                                adjust_charge='valence_model',
                                valence_matrix=valence_matrix)

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

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

# Test atoms
testutils.validate(cluster.natoms , ref.natoms)
testutils.validate(cluster.znums  , ref.znums)
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-5)

# 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)

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