#  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/>.


# TODO: fractional with shells

from .xyz      import _XYZ
from ..dl_py2f import utils

class Cell(_XYZ):
    '''Class of crystal lattice'''

    from numpy  import zeros, float64, int64

    _kwargs = {
                'cartesian'         : zeros(shape=(1,3), dtype=float64),
                'consts'            : zeros(shape=(6,) , dtype=float64),
                'convention'        :'conventional',
                'cutoff'            : 49.99*1.889725989, # DL_POLY segfaults given larger cutoff
                'delr'              : 0.0,
                'fractional'        : zeros(shape=(1,3), dtype=float64),
                'natoms'            : 0,
                'ndimensions'       : 3,
                'strains'           : zeros(shape=(3,3), dtype=float64),
                'strain_derivatives': zeros(shape=(3,3), dtype=float64),
                'pbc'               :'',
                'vectors'           : zeros(shape=(3,3), dtype=float64),
                'volume'            : 0.0,
              }

    _internals = { '_defer_frac2cart': False }

#    synons = { 'natoms':'npts',
#               'cartesian':'xyz' }
#    _synons = {
#                'cartesian':'xyz',
#                'natoms':'npts',
#              }

    # 'cartesian' should not appear in `_fields` as it is a reference to frag.coords
    _fields = 'fractional'

#    _priorities = 'dim', 'strains'


# CELL PROPERTIES #


    @property
    def stress_tensor(self):
        '''Stress_tensor (in GPa): currently converted from strain derivatives'''

        from ..data import unitconvert

        # 1 GPa = 1E9 J/m**3
        stress = 1.0E21*unitconvert.Ha2J*self.strain_derivatives/(self.volume*unitconvert.Bohr2Angstrom**3)

        return stress


#    @stress_tensor.setter
#    def stress_tensor(self, val):
#        '''Setter of stress_tensor'''
#
#        self._stress_tensor[:] = val


    @property
    def consts(self):
        '''Lattice constants (a/b/c in a.u.)'''

        from numpy   import arccos, count_nonzero, dot, linalg, pi, isnan

        # very weird problem: self._consts and self._vectors will become unwriteable somehow!
        self._consts.setflags(write=True)

        # skip if vectors not yet set
        if count_nonzero(self._vectors) == 0:
            return self._consts

    
        # conventional definition
        self._consts[0] = linalg.norm(self._vectors[:,0])
        self._consts[1] = linalg.norm(self._vectors[:,1])
        self._consts[2] = linalg.norm(self._vectors[:,2])
        # alpha: between vec(b) and vec(c), check for nan
        if self._consts[1]*self._consts[2]:
            self._consts[3] = arccos(dot(self._vectors[:,1], self._vectors[:,2])/self._consts[1]/self._consts[2])/pi*180.0
        # beta: between vec(a) and vec(c), check for nan
        if self._consts[0]*self._consts[2]:
            self._consts[4] = arccos(dot(self._vectors[:,0], self._vectors[:,2])/self._consts[0]/self._consts[2])/pi*180.0
        # gamma: between vec(a) and vec(b), check for nan
        if self._consts[0]*self._consts[1]:
            self._consts[5] = arccos(dot(self._vectors[:,0], self._vectors[:,1])/self._consts[0]/self._consts[1])/pi*180.0

        setattr(self, 'convention', "conventional") 

        # use 0.0 instead of float('nan') for invalid constants (nan would not appear because already checked above)
        self._consts[isnan(self._consts)] = 0.0

        self._consts.setflags(write=False)
    
        return self._consts


    @consts.setter
    def consts(self, val):
        '''Setter of lattice constants (a/b/c in a.u.)'''

        from numpy import count_nonzero

        # for I/O
        try:
            self._consts.setflags(write=True)
            self._consts[:] = val
        # for initialising
        # NB: this could change the type of Cell.consts, e.g. to a list
        except:
            self._consts = val

        # if deferred, hold off frac2cart() until now
        if self._defer_frac2cart and count_nonzero(self._master['fractional']):
            self.frac2cart()

        # to protect from contamination
        self._consts.setflags(write=False)


    @property
    def constants(self):
        '''Synonym of consts'''

        return self.consts


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

        self.consts = val


    @property
    def a(self):
        '''Lattice constant a (in a.u.)'''

        return self._consts[0]


    @a.setter
    def a(self, val):
        '''Setter of lattice constant a (in a.u.)'''

        try:
            self._consts[0] = val
        except:
            pass


    @property
    def b(self):
        '''Lattice constant b (in a.u.)'''

        return self._consts[1]


    @b.setter
    def b(self, val):
        '''Setter of lattice constant b (in a.u.)'''

        try:
            self._consts[1] = val
        except:
            pass


    @property
    def c(self):
        '''Lattice constant c (in a.u.)'''

        return self._consts[2]


    @c.setter
    def c(self, val):
        '''Setter of lattice constant c (in a.u.)'''

        try:
            self._consts[2] = val
        except:
            pass


    @property
    def alpha(self):
        '''Lattice constant alpha (in degrees)'''

        return self._consts[3]


    @alpha.setter
    def alpha(self, val):
        '''Setter of lattice constant alpha (in degrees)'''

        try:
            self._consts[3] = val
        except:
            pass


    @property
    def beta(self):
        '''Lattice constant beta (in degrees)'''

        return self._consts[4]


    @beta.setter
    def beta(self, val):
        '''Setter of lattice constant beta (in degrees)'''

        try:
            self._consts[4] = val
        except:
            pass


    @property
    def gamma(self):
        '''Lattice constant gamma (in degrees)'''

        return self._consts[5]


    @gamma.setter
    def gamma(self, val):
        '''Setter of lattice constant gamma (in degrees)'''

        try:
            self._consts[5] = val
        except:
            pass


    @property
    def A(self):
        '''Vector A of cell'''

        return self.vectors[:,0]


    @property
    def B(self):
        '''Vector B of cell'''

        return self.vectors[:,1]


    @property
    def C(self):
        '''Vector C of cell'''

        return self.vectors[:,2]


    @property
    def vectors(self):
        '''Basis vectors (in a.u.)'''
 
        from numpy    import count_nonzero, pi
        from warnings import resetwarnings, simplefilter

        # YL 01/02/2021: raise error if coordinates are wrong
        simplefilter('error')

        # in Py-ChemShell (Python):       in Tcl-ChemShell (Fortran):
        # vec[x/y/z, a/b/c]              vec(a/b/c, x/y/z)
        # [ vec[0,0], vec[0,1], vec[0,2] ]     ( vec(1,1), vec(2,1), vec(3,1) )
        # [ vec[1,0], vec[1,1], vec[1,2] ]     ( vec(1,2), vec(2,2), vec(3,2) )
        # [ vec[2,0], vec[2,1], vec[2,2] ]     ( vec(1,3), vec(2,3), vec(3,3) )

        # a, b, c in a.u. and alpha, beta, gamma in radians
        # keep private as not used elsewhere
        def _fromCrystal(a, b, c, alpha, beta, gamma):
            '''Calculate cell vectors according to CRYSTAL cell parameters definition'''

            from numpy import cos, sin

            # very weird problem: self._consts and self._vectors will become unwriteable somehow!
            self._vectors.setflags(write=True)

            self._vectors[:,:] = 0.0
            self._vectors[0,0] = a*sin(gamma)
            self._vectors[0,2] = c*cos(beta)
            self._vectors[1,0] = a*cos(gamma)
            self._vectors[1,1] = b
            self._vectors[1,2] = c*(cos(alpha)-cos(beta)*cos(gamma))/sin(gamma)
            self._vectors[2,2] = c*(1+2*cos(alpha)*cos(beta)*cos(gamma)-cos(alpha)**2-cos(beta)**2-cos(gamma)**2)**0.5/sin(gamma)
            # in Tcl-ChemShell:
            #     M(3,3) = c*sqrt(1.0-((cos(beta)-cos(alpha)*cos(gamma))/sin(gamma))**2-cos(alpha)**2)
            #            = c/sin(gamma)*sqrt(sin(gamma)**2-cos(beta)**2-cos(alpha)**2*cos(gamma)**2+2*cos(alpha)*cos(beta)*cos(gamma)-cos(alpha)**2*sin(gamma)**2)
            #            = c/sin(gamma)*sqrt(1.0-cos(gamma)**2-cos(alpha)**2*(cos(gamma)**2+sin(gamma)**2)-cos(beta)**2+2*cos(alpha)*cos(beta)*cos(gamma))
            #            = c/sin(gamma)*sqrt(1.0-cos(gamma)**2-cos(alpha)**2-cos(beta)**2+2*cos(alpha)*cos(beta)*cos(gamma))
            #            = V/(a*b*sin(gamma))

        # a, b, c in a.u. and alpha, beta, gamma in radians
        # keep private as not used anywhere else
        def _fromConventional(a, b, c, alpha, beta, gamma):
            '''Calculate cell vectors according to conventional cell parameters definition'''

            from numpy import cos, sin, pi

            # very weird problem: self._consts and self._vectors will become unwriteable somehow!
            self._vectors.setflags(write=True)

            # vec(a): self._vectors[:,0]; vec(b): self._vectors[:,1]; vec(a): self._vectors[:,2]
            # vec(a) = a*{1, 0, 0}
            # vec(b) = b*{cos(gamma), sin(gamma), 0}
            # vec(c) = c*{cos(beta),
            #           [cos(alpha)-cos(beta)*cos(gamma)]/sin(gamma),
            #           [1+2*cos(alpha)*cos(beta)*cos(gamma)-cos**2(alpha)-cos**2(beta)-cos**2(gamma)]**0.5/sin(gamma)}
            self._vectors[:,:] = 0.0
            self._vectors[0,0] = a
            self._vectors[0,1] = b*cos(gamma)
            self._vectors[0,2] = c*cos(beta)
            self._vectors[1,1] = b*sin(gamma)
            self._vectors[1,2] = c*(cos(alpha)-cos(beta)*cos(gamma))/sin(gamma)
            # YL 11/02/2021: do NOT use self.dim which will trigger rewriting of self.consts
            if count_nonzero(self._consts) == 6:
                self._vectors[2,2] = c*(1+2*cos(alpha)*cos(beta)*cos(gamma)-cos(alpha)**2-cos(beta)**2-cos(gamma)**2)**0.5/sin(gamma)

        # if constants not yet defined
        if count_nonzero(self._consts) == 0:
            return self._vectors

        # don't recalculate vectors if already defined
        if count_nonzero(self._vectors) != 0:
# YL 16/01/2020: returning a [2,2] matrix would be ideal but DL_PY2F can't work
#            if self.dim == '3D':
#                return self._vectors
#            elif self.dim == '2D':
#                return self._vectors[0:2,0:2]

            return self._vectors

        # conventions and corresponding methods of calculating 3x3 matrix
        selectcases = { 'conventional':_fromConventional,
                        'crystal'     :_fromCrystal,
                         None         :_fromCrystal       # ad hoc, see file.py
                      }

        # choose method to calculate vectors (3x3) according to convention
        selectcases[self.convention](self._consts[0],
                                     self._consts[1],
                                     self._consts[2],
                                     self._consts[3]/180.0*pi,
                                     self._consts[4]/180.0*pi,
                                     self._consts[5]/180.0*pi)

        # YL 01/02/2021: reset 
        simplefilter('default')
        resetwarnings()

# YL 16/01/2020: returning a [2,2] matrix would be ideal but DL_PY2F can't work
#        if self.dim == '3D':
#            return self._vectors
#        elif self.dim == '2D':
#            return self._vectors[0:2,0:2]
        return self._vectors


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

        from numpy  import append, array

        # YL 23/06/2020: allow 6 numbers for 2D and 3 numbers for 1D case
        if array(val).size == 6:
            val = append(array(val).reshape(6,), [0.0,0.0,0.0])
        if array(val).size == 3:
            val = append(array(val).reshape(3,), [0.0,0.0,0.0,0.0,0.0,0.0])

        # for I/O
        try:
            utils.nputils.setField(self._vectors, val)
        # for initialising
        except:
            self._vectors = val


    @property
    def reciprocal(self):
        '''Recipical lattice vectors'''

        from numpy import cross, dot, float64, pi, zeros

        v = self.vectors

        # YL 16/02/2020: added 2D case
        if self.ndimensions == 2:
            reci = zeros(shape=(2,2), dtype=float64)
            # increase the dimensionality to get the normal vector, because cross product only makes sense in 3D
            v0 = zeros((3,), dtype=float64)
            v1 = zeros((3,), dtype=float64)
            # YL 23/06/2020: keep the vectors array's shape as (3,3)
#            v0[0:2] = v[:,0]
#            v1[0:2] = v[:,1]
            v0[:] = v[:,0]
            v1[:] = v[:,1]
            normal = cross(v0,v1)
            reci[:,0] = (2.0*pi*cross(v1,normal)/dot(v0,cross(v1,normal)))[0:2]
            reci[:,1] = (2.0*pi*cross(normal,v0)/dot(v1,cross(normal,v0)))[0:2]
        elif self.ndimensions == 3:
            reci = zeros(shape=(3,3), dtype=float64)
            reci[:,0] = 2.0*pi*cross(v[:,1],v[:,2])/dot(v[:,0],cross(v[:,1],v[:,2]))
            reci[:,1] = 2.0*pi*cross(v[:,2],v[:,0])/dot(v[:,1],cross(v[:,2],v[:,0]))
            reci[:,2] = 2.0*pi*cross(v[:,0],v[:,1])/dot(v[:,2],cross(v[:,0],v[:,1]))

        return reci


    @property
    def fractional(self):
        '''Fractional coordinates'''

        from numpy import count_nonzero

        # YL 16/11/2019
        if self.natoms and not count_nonzero(self._master.fractional):
            self.cart2frac()

        return self._master.fractional


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

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

        self.frac2cart()


    @property
    def volume(self):
        '''Volume of the cell'''

        from numpy import absolute, cross, dot, linalg

        def __v1d(_vectors):
            return linalg.norm(_vectors[0])

        def __v2d(_vectors):
            return linalg.norm(cross(_vectors[0], _vectors[1]))

        def __v3d(_vectors):
            return absolute(dot(_vectors[0], cross(_vectors[1], _vectors[2])))

        selectcases = { 0 : __v3d,
                        1 : __v1d,
                        2 : __v2d,
                        3 : __v3d,
                      }

        return selectcases[self.ndimensions](self.vectors)


    @property
    def cartesian(self):
        '''Cartesian coordinates (computed from fractional coordinates)'''

        # self._cartesian is a reference to frag.coords
        return self._cartesian


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

        # self._cartesian is a reference to frag.coords
        utils.nputils.setField(self._cartesian, val)

	# Cell.connect() is a reference to Fragment.connect(), defined in <class Fragment>
        try:
            self.connect()
        except:
            pass

        # YL 16/11/2019
        self.cart2frac()


    @property
    def ndimensions(self):
        '''Dimensionality'''

        from numpy import count_nonzero
    
        selectcases = { 6 : 3,
                        3 : 2,
                        1 : 1,
                        0 : 0 }
    
        return selectcases[count_nonzero(self.consts)]

    @property
    def dim(self):
        '''Deprecated  dimensionality'''

        selectcases = { 3 : '3D',
                        2 : '2D',
                        1 : '1D',
                        0 : '0D' }

        print('\n >>> WARNING: String cell.dim is deprecated, please update your code to use integer cell.ndimensions')

        return selectcases[self.ndimensions]


# END OF CELL PROPERTIES #
# CLASS METHODS#


    def strains2vectors(self, vectors=[], ndimensions=3):
        '''Apply strains to cell to generate new vectors'''

        from numpy import count_nonzero, tensordot

        # see: http://cmt.dur.ac.uk/sjc/thesis_mcg/node113.html
        if count_nonzero(self.consts) == 0:
            return

        if ndimensions == 2:
            if not len(vectors):
                self.vectors[:2] = tensordot(self.strains[:2,:2], self.vectors[:2], axes=1)
            else:
                self.vectors[:2] = tensordot(self.strains[:2,:2], vectors[:2], axes=1)

        elif ndimensions == 3:
            if not len(vectors):
                self.vectors[:] = tensordot(self.strains, self.vectors[:], axes=1)
            else:
                self.vectors[:] = tensordot(self.strains, vectors[:], axes=1)


    def consts2vectors(self):
        '''Convert lattice constants (a/b/c in a.u.) to basis vectors (in a.u.)'''

        from numpy import count_nonzero, pi

        # in Py-ChemShell (Python):       in Tcl-ChemShell (Fortran):
        # vec[x/y/z, a/b/c]                   vec(a/b/c, x/y/z)
        # [ vec[0,0], vec[0,1], vec[0,2] ]    ( vec(1,1), vec(2,1), vec(3,1) )
        # [ vec[1,0], vec[1,1], vec[1,2] ]    ( vec(1,2), vec(2,2), vec(3,2) )
        # [ vec[2,0], vec[2,1], vec[2,2] ]    ( vec(1,3), vec(2,3), vec(3,3) )

        # a, b, c in a.u. and alpha, beta, gamma in radians
        # keep private as not used anywhere else
        def _fromCrystal(a, b, c, alpha, beta, gamma):
            '''Calculate cell vectors according to CRYSTAL cell parameters definition'''

            from numpy import cos, sin

            self.vectors[:,:] = 0.0
            self.vectors[0,0] = a*sin(gamma)
            self.vectors[0,2] = c*cos(beta)
            self.vectors[1,0] = a*cos(gamma)
            self.vectors[1,1] = b
            self.vectors[1,2] = c*(cos(alpha)-cos(beta)*cos(gamma))/sin(gamma)
            self.vectors[2,2] = c*(1+2*cos(alpha)*cos(beta)*cos(gamma)-cos(alpha)**2-cos(beta)**2-cos(gamma)**2)**0.5/sin(gamma)
            # in Tcl-ChemShell:
            #     M(3,3) = c*sqrt(1.0-((cos(beta)-cos(alpha)*cos(gamma))/sin(gamma))**2-cos(alpha)**2)
            #            = c/sin(gamma)*sqrt(sin(gamma)**2-cos(beta)**2-cos(alpha)**2*cos(gamma)**2+2*cos(alpha)*cos(beta)*cos(gamma)-cos(alpha)**2*sin(gamma)**2)
            #            = c/sin(gamma)*sqrt(1.0-cos(gamma)**2-cos(alpha)**2*(cos(gamma)**2+sin(gamma)**2)-cos(beta)**2+2*cos(alpha)*cos(beta)*cos(gamma))
            #            = c/sin(gamma)*sqrt(1.0-cos(gamma)**2-cos(alpha)**2-cos(beta)**2+2*cos(alpha)*cos(beta)*cos(gamma))
            #            = V/(a*b*sin(gamma))

        # a, b, c in a.u. and alpha, beta, gamma in radians
        # keep private as not used anywhere else
        def _fromConventional(a, b, c, alpha, beta, gamma):
            '''Calculate cell vectors according to conventional cell parameters definition'''

            from numpy import cos, sin

            # vec(a): self.vectors[:,0]; vec(b): self.vectors[:,1]; vec(a): self.vectors[:,2]
            # vec(a) = a*{1, 0, 0}
            # vec(b) = b*{cos(gamma), sin(gamma), 0}
            # vec(c) = c*{cos(beta),
            #           [cos(alpha)-cos(beta)*cos(gamma)]/sin(gamma),
            #           [1+2*cos(alpha)*cos(beta)*cos(gamma)-cos**2(alpha)-cos**2(beta)-cos**2(gamma)]**0.5/sin(gamma)}
            self.vectors[:,:] = 0.0
            self.vectors[0,0] = a
            self.vectors[0,1] = b*cos(gamma)
            self.vectors[0,2] = c*cos(beta)
            self.vectors[1,1] = b*sin(gamma)
            self.vectors[1,2] = c*(cos(alpha)-cos(beta)*cos(gamma))/sin(gamma)
            self.vectors[2,2] = c*(1+2*cos(alpha)*cos(beta)*cos(gamma)-cos(alpha)**2-cos(beta)**2-cos(gamma)**2)**0.5/sin(gamma)

        if count_nonzero(self.consts) == 0:
            return

        # conventions and corresponding methods of calculating 3x3 matrix
        selectcases = { "conventional":_fromConventional,
                        "crystal"     :_fromCrystal,
                         None         :_fromCrystal       # ad hoc, see file.py
                      }

        # choose method to calculate vectors (3x3) according to convention
        selectcases[self.convention](self.consts[0],
                                     self.consts[1],
                                     self.consts[2],
                                     self.consts[3]/180.0*pi,
                                     self.consts[4]/180.0*pi,
                                     self.consts[5]/180.0*pi)

        self.calcProperties()


## FIXME: obsolete
#    def calcProperties(self, maxcoords=None):
#
#        from chemsh.base import boundary
#
#        if maxcoords:
#            boundary.boundaries[self.dim](self, maxcoords)


    @property
    def _f2c_matrix(self):
        '''Matrix for converting fractional to cartesian'''

        from numpy import array, cos, dot, sin

        # constants in a.u. and Rad
        a     = self.a
        b     = self.b
        c     = self.c
        alpha = self.alpha/180.0
        beta  = self.beta/180.0
        gamma = self.gamma/180.0
        v     = self.volume

        m = [[a  , b*cos(gamma), c*cos(beta)],
             [0.0, b*sin(gamma), c*(cos(alpha)-cos(beta)*cos(gamma))/sin(gamma)],
             [0.0, 0.0         , v/a/b/sin(gamma)]]

        return array(m)


    def frac2cart(self):
        '''Fractional to cartesian'''

        from numpy import array, cos, count_nonzero, dot, linalg, sin

        # do nothing if fractional is not ready
        if not count_nonzero(self._master['fractional']):
            return

        # defer calculation of Cartesian if constants/vectors not yet ready
        if not count_nonzero(self.consts):
            self._defer_frac2cart = True
 
        # YL 16/11/2019
        self.cartesian = dot(self.fractional, self.vectors.T)
        #self.cartesian = dot(self._master['fractional'], self._f2c_matrix)

        # YL 05/05/2018: removed use of findCloserSites() in Fragment
        self.findCloserSites()


    # YL 16/11/2019
    def cart2frac(self):
        '''Cartesian to fractional'''

        from numpy import array, cos, count_nonzero, dot, float64, pi, sin, hstack, zeros

        # do nothing if cartesian is not ready
        if not count_nonzero(self.cartesian):
            return

        if not count_nonzero(self.consts):
            self.vectors2consts()

        # resize if necessary
        if self._master.shape[0] != self.cartesian.shape[0]:
            self.resize(self.cartesian.shape[0])

        if self.ndimensions == 3:
            fractional = dot(self.cartesian, self.reciprocal)/pi/2.0
            self._master.setField(fractional, field='fractional')
        elif self.ndimensions == 2:
            fractional = dot(self.cartesian[:,0:2], self.reciprocal)/pi/2.0
#            z = zeros(shape=(fractional.shape[0],1),dtype=float64)
            # YL 20/01/2020: the z-axis is not periodic so 2D fractional coordinates are actually cartesian in the z-direction
            z = self.cartesian[:,2].reshape((fractional.shape[0],1))
            self._master.setField(hstack((fractional,z)), field='fractional')


    @property
    def _c2f_matrix(self):
        '''Matrix for converting from cartesian to fractional'''

        from numpy import array, cos, dot, sin

        # constants in a.u. and Rad
        a     = self.a
        b     = self.b
        c     = self.c
        alpha = self.alpha/180.0
        beta  = self.beta/180.0
        gamma = self.gamma/180.0
        v     = self.volume

        m = [[1/a, -cos(gamma)/a/sin(gamma), b*c*(cos(alpha)*cos(gamma)-cos(beta))/v/sin(gamma)],
             [0.0, 1.0/b/sin(gamma)        , a*c*(cos(beta)*cos(gamma)-cos(alpha))/v/sin(gamma)],
             [0.0, 0.0                     , a*b*sin(gamma)/v]]

        return array(m)



    def vectors2consts(self):
        '''Convert basis vectors (in a.u.) to lattice constants (a/b/c in a.u.)'''

        from numpy import arccos, count_nonzero, dot, linalg, pi
    
        if count_nonzero(self.vectors) == 0:
            return
    
        # conventional definition
        self.consts[0] = linalg.norm(self.vectors[:,0])
        self.consts[1] = linalg.norm(self.vectors[:,1])
        self.consts[2] = linalg.norm(self.vectors[:,2])
        # alpha: between vec(b) and vec(c)
        self.consts[3] = arccos(dot(self.vectors[:,1], self.vectors[:,2])/self.consts[1]/self.consts[2])/pi*180.0
        # beta: between vec(a) and vec(c)
        self.consts[4] = arccos(dot(self.vectors[:,0], self.vectors[:,2])/self.consts[0]/self.consts[2])/pi*180.0
        # gamma: between vec(a) and vec(b)
        self.consts[5] = arccos(dot(self.vectors[:,0], self.vectors[:,1])/self.consts[0]/self.consts[1])/pi*180.0
    
        setattr(self, 'convention', "conventional") 
    
        self.calcProperties()


# END OF CLASS METHODS#


    # original DL_Poly method to calculate cosines: bbb[4] = cos(gamma), bbb[5] = cos(beta), and bbb[6] = cos(alpha)
    # having height vector h parallel to vectors cross product a×b: __h_ab = V/|a×b|
    @property
    def __h_bc(self):
        '''Perpendicular heights against surface (area) |b×c|'''

        from numpy import cross, linalg

        def __h1d(_vectors): 
            return 0.0

        def __h2d(_vectors): 
            return self.volume/linalg.norm(_vectors[1,:])

        def __h3d(_vectors): 
            return self.volume/linalg.norm(cross(_vectors[1,:], _vectors[2,:]))


        selectcases = { 0 : __h3d,
                        1 : __h1d,
                        2 : __h2d,
                        3 : __h3d }

        # y and z
        return selectcases[self.ndimensions](self.vectors)
    

    @property
    def __h_ca(self):
        '''Perpendicular heights against surface (area) |a×c|'''

        from numpy import cross, linalg

        def __h1d(_vectors): 
            return 0.0

        def __h2d(_vectors): 
            return self.volume/linalg.norm(_vectors[0,:])

        def __h3d(_vectors): 
            return self.volume/linalg.norm(cross(_vectors[2,:], _vectors[0,:]))

        selectcases = { 0 : __h3d,
                        1 : __h1d,
                        2 : __h2d,
                        3 : __h3d }

        # x and z
        return selectcases[self.ndimensions](self.vectors)
    

    @property
    def __h_ab(self):
        '''Perpendicular heights against surface (area) |a×b|'''

        from numpy import cross, linalg

        def __h2d(_vectors): 
            return 0.0 

        def __h3d(_vectors): 
            return self.volume/linalg.norm(cross(_vectors[0,:], _vectors[1,:]))


        selectcases = { 0 : __h3d,
                        1 : __h2d,
                        2 : __h2d,
                        3 : __h3d }

        # x and y
        return selectcases[self.ndimensions](self.vectors)


    @property
    def conncutoff(self):
        '''Cutoff for connectivity'''

        return self._conncutoff


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

        self._conncutoff = val
    

    # YL 04/02/2021: only used by DL_POLY so far
    @property
    def cutoff(self):
        '''Non-bonded cutoff (in a.u.) according to number of dimensions (different from Fragment.cutoff which is for connectivity!)'''

        from ..data import unitconvert

        # in an older version we used a dict to predefine the formulas, but it is not suitable now
        # since we have adapted @property and __h_bc/ca/ab are also calculated for even an empty (0D) cell (devision by zero)
        if self.ndimensions == 0:
            # YL 03/09/2020: used to be hardcoded 50.0
            # __radius is defined by Fragbase.cell
            cutoff = (self.__radius+1.0)*unitconvert.Angstrom2Bohr
        elif self.ndimensions == 3:
            cutoff = min(self.__h_bc, self.__h_ca, self.__h_ab)/2.0
        elif self.ndimensions == 2:
            # see tcl-chemsh/src/dl_poly/param.c
            cutoff = min(self.__h_bc, self.__h_ca)/2.0
        # TODO
        elif self.ndimensions == 1:
            cutoff = (self.__h_bc, self.__h_ca, self.__h_ab)/2.0
             
        # NB: type(self.cutoff) is float64 (rather than Python float!)
        cutoff -= 0.01*unitconvert.Angstrom2Bohr
        cutoff += self.delr
    
        return cutoff


    # YL 03/02/2021: this cutoff is only used for switching on/off the cell lists/verlet list algorithm of calculating connectivity
    #                not to be confused with Cell.cutoff!
    @property
    def cell_lists_cutoff(self):
        '''Cutoff for computing connectivity (normally is a reference to Fragment.cell_lists_cutoff)'''

        return self._cell_lists_cutoff


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

        self._cell_lists_cutoff = val
    

# YL: TODO: height_cutoff_ratio is fixed to 1 in Tcl-ChemShell as well
# in Tcl-ChemShell the way to obtain it is commented out:
#     irat = nint((rcut+delr)/rlink)
#     irat = min(max(irat,1),5)
# where rlink is not define anywhere

    @property
    def _nlinkcells(self):
        ''''''

        return self._nlinkcells_x*self._nlinkcells_y*self._nlinkcells_z


    @property
    def _use_cell_lists(self):
        '''Whether to use the cell lists algorithm (O(N)) for calculating the connectivity (rather than the O(N**2) verlet list method)'''

        from numpy import any, array

        # YL 01/02/2021: in Tcl-ChemShell it is switched off when self._nlinkcells == 27
        #                and it was commented by PS that there was some problems (we don't know what)
        # but here we just use it until we find real problems
        # there is a redundant check of enough linkcells in parlink() of nlist_builders_module.f in DL_POLY Classic but we don't need to do it
        if any(array([self._nlinkcells_x,self._nlinkcells_y,self._nlinkcells_z])<3) \
            or self.pbc in [ 'truncated_octahedron', 'rhombic_dodecahedron', 'hexagonal_prism' ]:
            return False
        else:
            return True


    @property
    def _height_cutoff_ratio(self):

        return 1


    @property
    def _nlinkcells_x(self):
        ''''''

        return max((2*self._height_cutoff_ratio+1), int(self.__h_bc*float(self._height_cutoff_ratio)/self.cell_lists_cutoff))


    @property
    def _nlinkcells_y(self):
        ''''''

        return max((2*self._height_cutoff_ratio+1), int(self.__h_ca*float(self._height_cutoff_ratio)/self.cell_lists_cutoff))


    @property
    def _nlinkcells_z(self):
        ''''''

        return max((2*self._height_cutoff_ratio+1), int(self.__h_ab*float(self._height_cutoff_ratio)/self.cell_lists_cutoff))


    @property
    def _nlinkcells_xyz(self):
        '''Array of number of linkcells along X, Y, and Z'''

        from numpy import array

        # rule of Python int() is the same as that of Fortran
        # in chemsh3/src/dl_poly/parset.f it is:
        #     ilx = max(3,int(celprp(7)/cut)), while
        # in chemsh3/src/dl_poly/parlink.f it is done in two steps:
        #     ilx = int(max(1.d0,celprp(7)*dble(irat)/(rcut+delr)))
        #     ilx = max(ilx,2*irat+1)
        # it can be proven that they are mathematically identical when irat=1:
        #     http://www.wolframalpha.com/input/?i=y%3Dmax%283%2Cfloor%28x%29%29
        #     http://www.wolframalpha.com/input/?i=y%3Dmax%283%2Cfloor%28max%281.0%2Cx%29%29%29
        # so we can reuse this method calcNumCells()
        # ncells along three axises
        try:
            return self.__nlinkcells_xyz
        except:
            self.__nlinkcells_xyz = array([ self._nlinkcells_x, self._nlinkcells_y, self._nlinkcells_z ], dtype=int)
            return self.__nlinkcells_xyz

    # mxcell in DL_POLY Classic
    @property
    def _max_nlinkcells(self):
        ''''''

        if self._use_cell_lists:
            return (self._nlinkcells_x+2)*(self._nlinkcells_y+2)*(self._nlinkcells_z+2)
        else:
            return 1

    # see: parset() of setup_module.f of DL_POLY Classic
    @property
    def _maxlist(self):

        from numpy import pi
        from ..data import periodicdata

        if self._use_cell_lists:
            # periodicdata.nsubcells[1] = 14
            return periodicdata.nsubcells[1]*round(1.5*float(self.natoms)/float(self._nlinkcells))
        else:
            ratio=1.5*self.natoms/self.volume*(4.0/3.0*pi)*self.cell_lists_cutoff**3
            return min(int(round(ratio)), int((self.natoms+1)/2))
           

    # YL 01/02/2021: same as Fragbase._max_abs_coords
    @property
    def _max_abs_coords(self):
        '''Greatest absolute values of Cartesian coords'''

        from numpy import array, fabs

        _max_abs_coords = []

        for i in range(3):
            _max_abs_coords.append(fabs(self.cartesian)[:,i].max())

        return array(_max_abs_coords)

#    @property
#    def __radius(self):
#        ''''''
#    
#        return self.___radius
#
#
#    @__radius.setter
#    def __radius(self, val):
#        ''''''
#
#        self.___radius = val


    @property
    def _inverse_vecs(self):
        '''Inverse matrix of the unit cell vectors (3*3)'''

        from numpy  import amax, append, asarray, fabs, linalg
        from ..base import errors

        # YL 24/06/2020: Cell.vectors changed to shape (3,3) in all dimensionalities
        # 1D not supported in Tcl-ChemShell
        if self.ndimensions == 1:
#            vecs = append(asarray(self._vectors), [ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 ]).reshape((3,3))
             pass
        # 2D: see celldim() of tcl-chemsh/src/fragment/connectf.f
        # in TCS there are 3 vectors (the 3rd being dummy); whereas in PCS there are only 2
        elif self.ndimensions == 2:
            maxz = fabs(self.cartesian)[:,2].max()
#            vecs = append(asarray(self._vectors), [ 0.0, 0.0, 1.0 ]).reshape((3,3))
            vecs = self._vectors
            vecs[2,2] = 1.0
            # NB: this is the cutoff for connectivity rather than Cell.cutoff
            cutoff = self.conncutoff
            vecs[2,2] = amax([ maxz+2.0*cutoff, 3.0*cutoff, vecs[2,2]])
        else:
            vecs = self._vectors

        try:
            return linalg.inv(vecs)
        except:
            raise errors.ChemShCellError('failed to calculate the inverse matrix of unit cell vectors:\n', self._vectors)


    # a Python rewrite of find_equiv_coords (see: chemsh/src/cluster/coord_transform.f90 in Tcl-ChemShell)
    # in Tcl-ChemShell tol**2 = 1.D-4 (in a.u.)
    def findCloserSites(self, origins=0.0, tolerance=0.1):
        '''Find equivalent sites that are closer to the origin'''

        from itertools import product
        from numpy     import argwhere, dot, zeros
        from numpy     import sum    as np_sum
        from numpy     import square as np_square

        dists = self.cartesian - origins

        sum_of_squares = np_sum(np_square(dists), axis=1)

        selectcases = { 3 : 3,
                        2 : 2,
                        1 : 1,
                        0 : 3,
                      }

        # translation in all directions
        for translation in product(range(-1,2), repeat=selectcases[self.ndimensions]):

            # paste 1D/2D/3D dot product into an array for 3D (shape=(3,))
            displacement = zeros((3,))
            displacement[:selectcases[self.ndimensions]] = dot(self.vectors[:,:selectcases[self.ndimensions]], translation)

# YL
#            displacement = dot(self.vectors, translation)

            # coodinates after translated
            coords = self.cartesian + displacement

            # indices of closer sites
            closers = argwhere(sum_of_squares - np_sum(np_square(coords-origins), axis=1) > tolerance).ravel()

            # replace
            if closers.size:
                self.cartesian[closers] = coords[closers]
                # recalculate the distance wrt the updated coords
                sum_of_squares = np_sum(np_square(self.cartesian-origins), axis=1)


    def applyConnectPBC(self, version='DL_POLY'):
# YL: check the literature if this should be max(2.0*cellObj.greatest_abs_coords[2]+cellObj.cutoff, 3.0*cellObj.cutoff) or max(2.0*cellObj.greatest_abs_coords[2]+cellObj.cutoff, 3.0*cellObj.cutoff, cellObj.vectors[2,2])! though normally this does not matter since cellObj.vectors[2,2] can hardly be greater than 3.0*cellObj.cutoff
        # imcon = 0 (0D/no)

        # throws read-only error
        self._vectors.setflags(write=True)

        # imcon = 0
        def _NoBound():
            '''Apply no boundary condition'''
       
            # will be 566.86110492033 if applied to an empty cell. check your result!
       
            # here was a mistake (applying @property decorators in <cell.Cell>)!
            # after self.vectors[0,0] is reset, self.vectors[1,1] and [2,2] cannot be defined!
            # because self.consts has new values which affect self.ndims and self.cell_lists_cutoff!
            # so modify no self.vectors but only change self._vectors!
            lbuff = [ max(2.0*self._max_abs_coords[i]+self.cell_lists_cutoff, 3.0*self.cell_lists_cutoff) for i in range(3) ]
       
            self._vectors[0,0] = lbuff[0]
            self._vectors[1,1] = lbuff[1]
            self._vectors[2,2] = lbuff[2]

        # imcon = 1
        def _Cubic():
            pass
       
        # imcon = 2
        def _Orthorhombic():
            pass
       
        # imcon = 3 (3D/parallelepiped)
        def _Parallelepiped():
            pass
    
        # imcon = 4
        def _TruncatedOctahedron():
            self._loglink = False
       
        # imcon = 5
        def _RhombicDodecahedron():
            self._loglink = False
       
        # imcon = 6 (2D/x-y_parallelogram)
        def _XYParallelogram():

            from numpy        import asarray, dot, round
            from numpy.linalg import det

            vecs = self._vectors
            v = asarray([ [ vecs[0,0], vecs[0,1] ], [ vecs[1,2], vecs[1,1] ] ])
            d = det(v)
            r0 =  1.0/d*vecs[1,2]
            r1 = -1.0/d*vecs[0,1]
            r3 = -1.0/d*vecs[1,1]
            r4 =  1.0/d*vecs[0,0]
            v2 = asarray([ [ r0, r1 ], [ r3, r4 ] ])
            prod = dot(self.cartesian[:,:2], v2)
            prod = prod - round(prod)
            self.cartesian[:,:2] = dot(prod, v2)

            # YL 01/02/2021: parlink() in nlist_builders_module.f of DL_POLY Classic and Tcl-ChemShell's connect();
            #                the latter looks like a duplicate
            # it's self._max_abs_coords[2]+2.0*self.cell_lists_cutoff in Tcl-ChemShell's celldim()
            if version == 'DL_POLY':
                self._vectors[2,2] = max(2.0*self._max_abs_coords[2]+self.cell_lists_cutoff, 3.0*self.cell_lists_cutoff, self._vectors[2,2])
            else:
                self._vectors[2,2] = max(self._max_abs_coords[2]+2.0*self.cell_lists_cutoff, 3.0*self.cell_lists_cutoff, self._vectors[2,2])

        # imcon = 7
        def _HexagonalPrism():
        # this appears in PF's code. is it true in DLP?
            self._loglink = False    
       

        # boundary conditions
        selectcases = { "no"                  : _NoBound,
                        "0D"                  : _NoBound,
                        "cubic"               : _Cubic,
                        "orthorhombic"        : _Orthorhombic,
                        "parallelepiped"      : _Parallelepiped,
                        "3D"                  : _Parallelepiped,
                        "truncated_octahedron": _TruncatedOctahedron,
                        "rhombic_dodecahedron": _RhombicDodecahedron,
                        "2D"                  : _XYParallelogram,     # no periodicity in z
                        "x-y_parallelogram"   : _XYParallelogram,     # no periodicity in z
                        "xy_parallelogram"    : _XYParallelogram,     # no periodicity in z
                        "hexagonal_prism"     : _HexagonalPrism }

        selectcases[self.pbc]()
       

    @property
    def pbc(self):
        '''Determine the periodic boundary condition'''

        from numpy import allclose
        from math  import isclose

        # user-defined
        if self._pbc:
            return self._pbc

        # automatic detection
        else:
            tol = 1.0e-3
            orthogonal  = allclose(self.consts[3:6], 90.0, rtol=0.0, atol=tol)
            equilateral = allclose(self.consts[0:3], self.consts[0:3].mean(), rtol=0.0, atol=tol)
            hexagonal   = allclose(self.consts[0:2], self.consts[0:2].mean(), rtol=0.0, atol=tol) and allclose(self.consts[5], 120.0, rtol=0.0, atol=tol)

            # YL FIXME: this is incomplete!
            if self.ndimensions == 3:
                if orthogonal and equilateral:
                    return 'cubic'
                # YL 10/05/2018: this is definitely wrong! but since we have done nothing for the PBC it does not affect the result for now
                elif orthogonal and not equilateral:
                    return 'orthorhombic'
                #
                #    return 'tetragonal'
                elif self.a != self.b != self.c != self.a and self.alpha != self.beta != self.gamma != self.alpha != 90.0:
                    return 'parallelepiped'
                # YL 31/01/2021: added hexagonal
                elif hexagonal:
                    return 'hexagonal_prism'
                # YL TODO: more PBCs
                else:
                    return 'parallelepiped'
    
            elif self.ndimensions == 2:
                return 'x-y_parallelogram'
    
            else:
                return 'no'


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

        from ..base import periodic

        self._pbc = val.strip().lower()

        self._pdc = periodic.verifyRequestedPBC(self._pbc, self.ndimensions)


    @property
    def imcon(self):
        '''DL_POLY's notation code for PBC

           imcon=0 no boundary conditions
                 1 standard cubic
                 2 orthorhombic
                 3 parallelepiped
                 4 truncated octahedron
                 5 rhombic dodecahedron
                 6 x-y parallelogram: no periodicity in z
                 7 hexagonal prism'''

        selectcases = { 'no'                  :0,
                        'cubic'               :1,
                        'orthorhombic'        :2,
                        'parallelepiped'      :3,
                        'truncated_octahedron':4,
                        'rhombic_dodecahedron':5,
                        'x-y_parallelogram'   :6,
                        'hexagonal_prism'     :7,
                      }

        return selectcases[self.pbc]


#    # YL 04/11/2020: we forgot to define `Cell.natoms`. `_XYZ.natoms` returns the shape of self._master so it was always 3
#    @property
#    def natoms(self):
#        '''Number of atoms'''
#
#        return self.cartesian.shape[0]
#
#
#    @natoms.setter
#    def natoms(self, val):
#        '''Dummy setter of natoms'''
#








