#from Numeric import *
#from LinearAlgebra import *
from numpy import *

from math import *
from mymath import *
import string


class Xdatcarreader:
    def __init__(self):
        self.data = []
        self.lattmat = []

    def readFile(self, f, blocksize=1):
        if type(f) == type(""):
            f = open(f, 'r')
        else:
            return

        dummynum = 0
        self.numconfig = 0
        local_dire = None
        coordtype = 'direct'
        calctype = 'nonselective'
        calctarget = 7
        atomtag = 0
        i = 0                                 # row counter
        lattmat_tmp = zeros((3, 3), float)

        # reads scaling factor and lattice parameters
        for line in f.readlines():
            i = i+1
            line = line.split()
            if i == 2:
                # factor for scaling of the lattice vectors
                self.scale = float(line[0])
            if i > 2 and i < 6:
                # lattice-vectors matrix
                lattmat_tmp[i-3] = [float(line[0]),
                                    float(line[1]), float(line[2])]
            if i == 6:
                if self.scale < 0:          # negative scale - required volume
                    self.volume = cross(lattmat_tmp[0], lattmat_tmp[1])
                    self.volume = abs(sum(self.volume*lattmat_tmp[2]))
                    self.scale = (abs(self.scale)/self.volume)**(1.0/3.0)
                lattmat_tmp = self.scale*lattmat_tmp
                self.lattmat.append(lattmat_tmp)
                self.volume = cross(lattmat_tmp[0], lattmat_tmp[1])
                self.volume = abs(sum(self.volume*lattmat_tmp[2]))
                self.tags_atoms = line
                lattmat_tmp = zeros((3, 3), float)
            if i == 7:
                katoms = []
                for j in range(len(line)):
                    katoms.append(int(line[j]))
                self.atoms = katoms
                self.types = katoms
                self.ntypes = len(self.types)
                # self.atoms=katoms
                self.numofatoms = sum(katoms)
                lattmat_tmp = zeros((3, 3), float)
                coords = zeros((self.numofatoms, 3), float)

            # if i>8:
            if i > 7:
                indx = (i-7) % (self.numofatoms+1)
                if (indx == 0):
                    indx = self.numofatoms+1
                # print i,indx,line,self.numofatoms+8
                if indx > 1:
                    # print indx
                    # print line
                    coords[indx-2][0] = float(line[0])
                    coords[indx-2][1] = float(line[1])
                    coords[indx-2][2] = float(line[2])
                    if indx == self.numofatoms+1:
                        if (dummynum % blocksize == 0):
                            self.data.append(coords)
                            coords = zeros((self.numofatoms, 3), float)
                            self.numconfig += 1
                        dummynum += 1

        for i in range(1, len(self.data)):
            for j in range(len(self.data[0])):
                for k in range(3):
                    if self.data[i][j][k]-self.data[i-1][j][k] > 0.5:
                        self.data[i][j][k] -= 1
                    if self.data[i][j][k]-self.data[i-1][j][k] < -0.5:
                        self.data[i][j][k] += 1

        f.close()

########################################################################


def shortest_dist(cartesians, lattmat, atom1, atom2, indx=0):
    """finds the shortest non-zero distance between two atoms
    """
    cart1 = cartesians[atom1]
    cart2 = cartesians[atom2]
    dists = []
    what = []

    for i in [-1, 0, 1]:
        for j in [-1, 0, 1]:
            for k in [-1, 0, 1]:
                trans = i*lattmat[0]+j*lattmat[1]+k*lattmat[2]
                point2 = cart2+trans
                dist = sum((cart1-point2)**2)**0.5
                # if dist>0.001:
                #  dists.append(dist)
                #  what.append([i,j,k])
                dists.append(dist)
                what.append([i, j, k])

    dists_sorted = 1*dists
    dists_sorted.sort()

    target = dists_sorted[indx]
    minindex = dists.index(target)

    # minindex=0
    # for i in range(1,len(dists)):
    #  if dists[i]<dists[minindex]: minindex=i
    return [[0, 0, 0], [what[minindex][0], what[minindex][1], what[minindex][2]]]


def read_irccar(f):
    f = open(f, 'r')
    i = 0
    ircdata = []
    dimX = 0
    for line in f.readlines():
        line = line.split()
        if i == 0:
            npoints = int(line[0])
        else:
            dummy = []
            for i in range(len(line)):
                dummy.append(float(line[i]))
            if dimX == 0:
                dimX = len(dummy)
                if dimX == 0:
                    print(koko)
            else:
                if len(dummy) != dimX:
                    print(koko)
            ircdata.append(dummy)
        i += 1
    ircdata = array(ircdata)
    f.close()
    return ircdata


def read_definition(f, crt, lattmat):
    """here we read allowed intervals 
    for reactant and product states
    """
    if type(f) == type(""):
        f = open(f, "r")
    inttags = []
    intwhat = []
    intwhere = []
    intcoefs = []
    intcnref = []
    intcncoefs = []
    complextype = []
    ctag = 0
    for line in f.readlines():
        if ctag == 0:
            line = string.split(line, chr(58))
            if ((line[0][0] == 'C') or (line[0][0] == 'c')):
                ctag = 1
                continue
            tagline = line[0]
            inttags.append(tagline)
            whatline = line[1].split()
            for i in range(len(whatline)):
                whatline[i] = int(whatline[i])
            intwhat.append(whatline)

            line2_dummy = line[2].split()
            # print line2_dummy[0]

            if len(line) == 4:
                line3_dummy = line[3].split()
                intcnref.append(float(line3_dummy[0]))
                intcncoefs.append(
                    [float(line3_dummy[1]), float(line3_dummy[2])])
            else:
                intcnref.append(None)
                intcncoefs.append([None, None])

            delement = array([0, 0, 0])

            #! minimal image convention
            if tagline == 'X' or tagline == 'Y' or tagline == 'Z':
                if line2_dummy[0] == 'S' or line2_dummy[0] == 'S1' or line2_dummy[0] == 'S2':
                    arr1 = [0, 0, 0]
                    delement = array([arr1])
                else:
                    arr1 = [int(line2_dummy[0]), int(
                        line2_dummy[1]), int(line2_dummy[2])]
                    delement = array([arr1])
            if tagline == 'R' or tagline == 'IR1' or tagline == 'IR6' or tagline == 'RX' or tagline == 'RY' or tagline == 'RZ':
                if line2_dummy[0] == 'S' or line2_dummy[0] == 'S1':
                    delement = array(shortest_dist(
                        crt, lattmat, whatline[0], whatline[1]))
                elif line2_dummy[0] == 'S2':
                    delement = array(shortest_dist(
                        crt, lattmat, whatline[0], whatline[1], 1))
                else:
                    arr1 = [int(line2_dummy[0]), int(
                        line2_dummy[1]), int(line2_dummy[2])]
                    arr2 = [int(line2_dummy[3]), int(
                        line2_dummy[4]), int(line2_dummy[5])]
                    delement = array([arr1, arr2])
                # print delement
            if tagline == 'W' or tagline == 'WP':
                if line2_dummy[0] == 'S' or line2_dummy[0] == 'S1':
                    delement = array(shortest_dist(
                        crt, lattmat, whatline[0], whatline[1]))
                elif line2_dummy[0] == 'S2':
                    delement = array(shortest_dist(
                        crt, lattmat, whatline[0], whatline[1], 1))
                else:
                    arr1 = [int(line2_dummy[0]), int(
                        line2_dummy[1]), int(line2_dummy[2])]
                    arr2 = [int(line2_dummy[3]), int(
                        line2_dummy[4]), int(line2_dummy[5])]
                    delement = array([arr1, arr2])

            if tagline == 'acosZR':
                if line2_dummy[0] == 'S' or line2_dummy[0] == 'S1':
                    delement = array(shortest_dist(
                        crt, lattmat, whatline[0], whatline[1]))
                elif line2_dummy[0] == 'S2':
                    delement = array(shortest_dist(
                        crt, lattmat, whatline[0], whatline[1], 1))
                else:
                    arr1 = [int(line2_dummy[0]), int(
                        line2_dummy[1]), int(line2_dummy[2])]
                    arr2 = [int(line2_dummy[3]), int(
                        line2_dummy[4]), int(line2_dummy[5])]
                    delement = array([arr1, arr2])
            if tagline == 'Xr' or tagline == 'Yr' or tagline == 'Zr':
                if line2_dummy[0] == 'S' or line2_dummy[0] == 'S1':
                    delement = array(shortest_dist(
                        crt, lattmat, whatline[0], whatline[1]))
                elif line2_dummy[0] == 'S2':
                    delement = array(shortest_dist(
                        crt, lattmat, whatline[0], whatline[1], 1))
                else:
                    arr1 = [int(line2_dummy[0]), int(
                        line2_dummy[1]), int(line2_dummy[2])]
                    arr2 = [int(line2_dummy[3]), int(
                        line2_dummy[4]), int(line2_dummy[5])]
                    delement = array([arr1, arr2])
            if tagline == 'M':
                if line2_dummy[0] == 'S' or line2_dummy[0] == 'S1':
                    delement1 = array(shortest_dist(
                        crt, lattmat, whatline[0], whatline[1]))
                    delement2 = shortest_dist(
                        crt, lattmat, whatline[1], whatline[2])
                    delement = [delement1[0], delement1[1],
                                delement1[1]+delement2[1]]
                elif line2_dummy[0] == 'S2':
                    delement1 = array(shortest_dist(
                        crt, lattmat, whatline[0], whatline[1]))
                    delement2 = shortest_dist(
                        crt, lattmat, whatline[1], whatline[2], 1)
                    delement = [delement1[0], delement1[1],
                                delement1[1]+delement2[1]]
                else:
                    arr1 = [int(line2_dummy[0]), int(
                        line2_dummy[1]), int(line2_dummy[2])]
                    arr2 = [int(line2_dummy[3]), int(
                        line2_dummy[4]), int(line2_dummy[5])]
                    arr3 = [int(line2_dummy[6]), int(
                        line2_dummy[7]), int(line2_dummy[8])]
                    delement = array([arr1, arr2, arr3])
            if tagline == 'A':
                if line2_dummy[0] == 'S' or line2_dummy[0] == 'S1':
                    delement1 = shortest_dist(
                        crt, lattmat, whatline[1], whatline[0])
                    delement2 = shortest_dist(
                        crt, lattmat, whatline[1], whatline[2])
                    delement = [delement1[1], delement1[0], delement2[1]]
                elif line2_dummy[0] == 'S2':
                    delement1 = shortest_dist(
                        crt, lattmat, whatline[1], whatline[0], 1)
                    delement2 = shortest_dist(
                        crt, lattmat, whatline[1], whatline[2], 1)
                    delement = [delement1[1], delement1[0], delement2[1]]
                else:
                    arr1 = [int(line2_dummy[0]), int(
                        line2_dummy[1]), int(line2_dummy[2])]
                    arr2 = [int(line2_dummy[3]), int(
                        line2_dummy[4]), int(line2_dummy[5])]
                    arr3 = [int(line2_dummy[6]), int(
                        line2_dummy[7]), int(line2_dummy[8])]
                    delement = array([arr1, arr2, arr3])
            if tagline == 'AL':
                arr4 = [0, 0, 0]
                if line2_dummy[0] == 'S' or line2_dummy[0] == 'S1':
                    delement1 = shortest_dist(
                        crt, lattmat, whatline[1], whatline[0])
                    delement2 = shortest_dist(
                        crt, lattmat, whatline[1], whatline[2])
                    if len(line2_dummy) == 4:
                        arr4 = [int(line2_dummy[1]), int(
                            line2_dummy[2]), int(line2_dummy[3])]
                    delement = [delement1[1], delement1[0], delement2[1], arr4]
                elif line2_dummy[0] == 'S2':
                    delement1 = shortest_dist(
                        crt, lattmat, whatline[1], whatline[0], 1)
                    delement2 = shortest_dist(
                        crt, lattmat, whatline[1], whatline[2], 1)
                    if len(line2_dummy) == 4:
                        arr4 = [int(line2_dummy[1]), int(
                            line2_dummy[2]), int(line2_dummy[3])]
                    delement = [delement1[1], delement1[0], delement2[1], arr4]
                else:
                    arr1 = [int(line2_dummy[0]), int(
                        line2_dummy[1]), int(line2_dummy[2])]
                    arr2 = [int(line2_dummy[3]), int(
                        line2_dummy[4]), int(line2_dummy[5])]
                    arr3 = [int(line2_dummy[6]), int(
                        line2_dummy[7]), int(line2_dummy[8])]
                    arr4 = [int(line2_dummy[9]), int(
                        line2_dummy[10]), int(line2_dummy[11])]
                    delement = array([arr1, arr2, arr3, arr4])

            if tagline == 'T' or tagline == 'Amid':
                if line2_dummy[0] == 'S' or line2_dummy[0] == 'S1':
                    delement1 = shortest_dist(
                        crt, lattmat, whatline[1], whatline[0])
                    delement2 = shortest_dist(
                        crt, lattmat, whatline[1], whatline[2])
                    delement3 = shortest_dist(
                        crt, lattmat, whatline[2], whatline[3])
                    delement = [delement1[1], delement1[0], delement2[1], array(
                        delement3[1])+array(delement2[1])]
                elif line2_dummy[0] == 'S2':
                    delement1 = shortest_dist(
                        crt, lattmat, whatline[1], whatline[0], 1)
                    delement2 = shortest_dist(
                        crt, lattmat, whatline[1], whatline[2], 1)
                    delement3 = shortest_dist(
                        crt, lattmat, whatline[2], whatline[3], 1)
                    delement = [delement1[1], delement1[0], delement2[1], array(
                        delement3[1])+array(delement2[1])]
                else:
                    arr1 = [int(line2_dummy[0]), int(
                        line2_dummy[1]), int(line2_dummy[2])]
                    arr2 = [int(line2_dummy[3]), int(
                        line2_dummy[4]), int(line2_dummy[5])]
                    arr3 = [int(line2_dummy[6]), int(
                        line2_dummy[7]), int(line2_dummy[8])]
                    arr4 = [int(line2_dummy[9]), int(
                        line2_dummy[10]), int(line2_dummy[11])]
                    delement = array([arr1, arr2, arr3, arr4])
            if tagline == 'RatioR':
                if line2_dummy[0] == 'S' or line2_dummy[0] == 'S1':
                    delement1 = array(shortest_dist(
                        crt, lattmat, whatline[0], whatline[1]))
                    delement2 = array(shortest_dist(
                        crt, lattmat, whatline[2], whatline[3]))
                    delement = [delement1[0], delement1[1],
                                delement2[0], delement2[1]]
                elif line2_dummy[0] == 'S2':
                    delement1 = array(shortest_dist(
                        crt, lattmat, whatline[0], whatline[1], 1))
                    delement2 = array(shortest_dist(
                        crt, lattmat, whatline[2], whatline[3], 1))
                    delement = [delement1[0], delement1[1],
                                delement2[0], delement2[1]]
                else:
                    arr1 = [int(line2_dummy[0]), int(
                        line2_dummy[1]), int(line2_dummy[2])]
                    arr2 = [int(line2_dummy[3]), int(
                        line2_dummy[4]), int(line2_dummy[5])]
                    arr3 = [int(line2_dummy[6]), int(
                        line2_dummy[7]), int(line2_dummy[8])]
                    arr4 = [int(line2_dummy[9]), int(
                        line2_dummy[10]), int(line2_dummy[11])]
                    delement = array([arr1, arr2, arr3, arr4])
            intwhere.append(delement)

        if ctag == 1:
            line = string.split(line, chr(58))
            comptype = string.split(line[1])
            complextype.append(comptype[0])
            coefs = string.split(line[0])
            for i in range(len(coefs)):
                coefs[i] = float(coefs[i])
            intcoefs.append(coefs)
    f.close()
    return inttags, intwhat, intwhere, intcoefs, intcnref, intcncoefs, complextype
