from scipy import *
from numpy import *
from visual import *

sumModesN = 6

def stringarraytoarray(lineread, sizeSquare):
    dummy2 = lineread.replace(',','')
    dummy3 = dummy2.replace(']','')
    dummy4 = dummy3.replace('[','')
    dummyFinal = dummy4.split()
    target = [None]*sumModesN
    m = 0
    for i in range(sumModesN):
        target[i] = [0]*sumModesN
        for j in range(sumModesN):
            if abs(float(dummyFinal[m]))>10**-15:
                target[i][j] = float(dummyFinal[m])
            else:
                target[i][j] = 0.0
            m += 1
    return array(target)


def returnSEigenValueArray(matrix,sizeSquare):       # for use in the generation of our transformation matrix
    eigenValues = linalg.eig(matrix)[0]
    newarray = zeros((sizeSquare,sizeSquare),Float)
    for p in range(sizeSquare):
        for q in range(sizeSquare):
            if p == q: newarray[p][q] = abs(eigenValues[p])**(-0.5)
            else: newarray[p][q] = 0.0
    return newarray


overlapFile = open("dumpOverlaps", mode='r')
overlapReadLine = overlapFile.readline()
overlapArray = stringarraytoarray(overlapReadLine,sumModesN)
overlapFile.close()
print "\nOverlap Array\n"
print shape(overlapArray)
print overlapArray

kineticFile = open("dumpKEoverlaps", mode='r')
kineticReadLine = kineticFile.readline()
kineticArray = stringarraytoarray(kineticReadLine,sumModesN)
kineticFile.close()
print "\nKinetic Array\n"
print shape(kineticArray)
print kineticArray

potentialFile = open("dumpNoverlaps", mode='r')
potentialReadLine = potentialFile.readline()
potentialArray = stringarraytoarray(potentialReadLine,sumModesN)
potentialFile.close()
print "\nPotential Array\n"
print shape(potentialArray)
print potentialArray

coulombicRepulsionFile = open("dumpCOoverlaps", mode='r')
coulombicRepulsionReadLine = coulombicRepulsionFile.readline()
coulombicRepulsionArray = stringarraytoarray(coulombicRepulsionReadLine,sumModesN)
coulombicRepulsionFile.close()
print "\nCoulombic Repulsion Array\n"
print shape(coulombicRepulsionArray)
print coulombicRepulsionArray

coulombicExchangeFile = open("dumpXCoverlaps", mode='r')
coulombicExchangeReadLine = coulombicExchangeFile.readline()
coulombicExchangeArray = stringarraytoarray(coulombicExchangeReadLine,sumModesN)
coulombicExchangeFile.close()
print "\nCoulombic Exchange Array\n"
print shape(coulombicExchangeArray)
print coulombicExchangeArray


Z = -1.0

U = linalg.eig(overlapArray)[1]
b = returnSEigenValueArray(overlapArray,sumModesN)
V = matrixmultiply(U,b)

realTotalArray = Z*potentialArray - 137.035999070*kineticArray + coulombicRepulsionArray + 2*coulombicExchangeArray

VmatrixDagger = transpose(conjugate(V))
transformedHamiltonian = matrixmultiply(matrixmultiply(VmatrixDagger,realTotalArray),V)

print "\nHamiltonian Array\n"
print transformedHamiltonian
output = linalg.eig(transformedHamiltonian)

print "\nEigenvalues:\n"

count = len(output[0])
eigenvalues = [None]*count
for i in range(count):
    eigenvalues[i] = float(sorted(output[0])[i])

print eigenvalues

count = len(output[0])
eigenvaluesEV = [None]*count
for i in range(count):
    eigenvaluesEV[i] = float(sorted(output[0])[i])*27.2113838685


print "\nEigenvalues in eV:\n"
print eigenvaluesEV




