#  Copyright (C) 2017 The authors of Py-ChemShell
#
#  This file is part of Py-ChemShell.
#
#  Py-ChemShell is free software: you can redistribute it and/or modify
#  it under the terms of the GNU Lesser General Public License as
#  published by the Free Software Foundation, either version 3 of the
#  License, or (at your option) any later version.
#
#  Py-ChemShell is distributed in the hope that it will be useful,
#  but WITHOUT ANY WARRANTY; without even the implied warranty of
#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#  GNU Lesser General Public License for more details.
#
#  You should have received a copy of the GNU Lesser General Public
#  License along with Py-ChemShell.  If not, see
#  <http://www.gnu.org/licenses/>.

from numpy import zeros, full, float64, int64
from .xyz  import _XYZ

class Links(_XYZ):
    '''Object of link atoms'''

# FIXME: rename to something like "_keywords"
    _kwargs = {
                'coords'   : zeros(shape=(1, 3), dtype=float64),
                'charges'  : zeros(shape=(1,)  , dtype=float64),
                'qmdiffs'  : zeros(shape=(1, 3), dtype=float64),
                'qldists'  : zeros(shape=(1,)  , dtype=float64),
                'qmdists'  : zeros(shape=(1,)  , dtype=float64),
                'gradients': zeros(shape=(1, 3), dtype=float64),
                'names'    : full((1,), 'H', dtype='S8'),
                'nlinks'   : 0,
                'znums'    : full((1,),  1, dtype=int64)
              }

    _internals = {
                 }

    synons = { 'nlinks': 'npts' }

    _synons = {
                'nlinks': 'npts',
                'coords': 'xyz'
              }

    _fields = 'names', 'znums', 'coords', 'charges', 'gradients', 'qldists', 'qmdists', 'qmdiffs'


    @property
    def names(self):
        '''BQ names'''

        return self._master.names


    @names.setter
    def names(self, val):
        '''Setter of names'''

        self._master.setField(val, field='names')


    @property
    def nlinks(self):
        '''Number of link atoms'''

        return self.natoms


    @nlinks.setter
    def nlinks(self, val):
        '''Setter of nlinks'''

        self.natoms = val


    @property
    def charges(self):
        '''Charges'''

        return self._master.charges


    @charges.setter
    def charges(self, val):
        '''Setter of charges'''

        self._master.setField(val, field='charges')


    @property
    def qldists(self):
        '''Distances between Q and L'''

        return self._master.qldists


    @qldists.setter
    def qldists(self, val):
        '''Setter of qldists'''

        self._master.setField(val, field='qldists')


    @property
    def qmdists(self):
        '''Distances between Q and M'''

        return self._master.qmdists


    @qmdists.setter
    def qmdists(self, val):
        '''Setter of qmdists'''

        self._master.setField(val, field='qmdists')


    @property
    def qmdiffs(self):
        '''Differences between Q and M'''

        return self._master.qmdiffs


    @qmdiffs.setter
    def qmdiffs(self, val):
        '''Setter of qmdiffs'''

        self._master.setField(val, field='qmdiffs')


    @property
    def gradients(self):
        '''Gradients on background point charges'''

        return self._master.gradients


    @gradients.setter
    def gradients(self, val):
        '''Setter of gradients on point charges'''

        self._master.setField(val, field='gradients')


    @property
    def znums(self):
        '''Atomic numbers: shape=(nlinks,1)'''

        self._master.setField(1, field='znums')

        return self._master.znums


    @property
    def bpairs(self):
        '''QM/MM boundaries'''

        return self._bpairs


    @bpairs.setter
    def bpairs(self, val):
        '''Setter of qatoms'''

        self._bpairs = val


    @property
    def boundariesstr(self):
        '''Boundaries represented in list of strings'''

        return [ str(bpair) for bpair in self.bpairs ]


    # "B matrix" in Tcl-ChemShell: see src/hybrid/link_atom.c
    def getMMatrix(self, bpair):
        '''Get the 3*3 matrix for Q'''

        from numpy import array, diagonal, outer

        idx = self.boundariesstr.index(str(bpair))

        qmdiffs = self.qmdiffs[idx]
        
        matrix = -self.qldists[idx]*outer(qmdiffs, qmdiffs)/self.qmdists[idx]**3

        # cannot assign function call such as: diagonal(matrix) = ...
        diag = diagonal(matrix)
        # YL: as of NumPy 1.10.1, numpy.diagonal() does not automatically return a writeable view as described in
        #     http://docs.scipy.org/doc/numpy-1.10.1/reference/generated/numpy.diagonal.html
        #     we need to change the flag
        diag.setflags(write=True)

        diag += self.qldists[idx]/self.qmdists[idx]

        return matrix


    # "C matrix" in Tcl-ChemShell: see src/hybrid/link_atom.c
    def getQMatrix(self, bpair):
        '''Get the 3*3 matrix for M'''

        from numpy import diagonal, outer

        matrix = -self.getMMatrix(bpair)

        diag = diagonal(matrix)
        diag.setflags(write=True)
        diag += 1.0

        return matrix

#    @znums.setter
#    def znums(self, val):
#        '''Setter of znums'''
#        self._znums = val


    __all__ = 'coords', 'charges', 'gradients'
