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

sumModesN = 6
Z = -1.0
alpha = 1/137.035999070
c = 300000000
hbar = 1.05457148*10**-34
e = 1.60217646*10**-19
bohrradius = 5.2917720859*10**-11
scale = hbar*c/(e*bohrradius)

print scale


def replacestructure(thefile, sizeSquare):
    dummyint = ''
    for i in xrange(sizeSquare):
        dummyint += thefile.readline()
    dummy2 = dummyint.replace(',','')
    dummy3 = dummy2.replace(']','')
    dummy4 = dummy3.replace('[','')
    return dummy4.split()

def stringarraytoarray(thefile, sizeSquare):
    dummy = replacestructure(thefile, sizeSquare)
    target = [None]*sizeSquare
    m = 0
    for i in xrange(sizeSquare):
        target[i] = [0]*sizeSquare
        for j in xrange(sizeSquare):
            if abs(float(dummy[m]))>10**-15:
                target[i][j] = float(dummy[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')
overlapArray = stringarraytoarray(overlapFile,sumModesN)
overlapFile.close()
print "\nOverlap Array\n"
print shape(overlapArray)
print overlapArray

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

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

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

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


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

realTotalArray = potentialArray + 2*kineticArray + 2*coulombicRepulsionArray #+ 4*coulombicExchangeArray

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

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

print "\nEigenvalues in eV:\n"

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

print eigenvalues1

realTotalArray = potentialArray + 2*kineticArray + 2*coulombicRepulsionArray - 4*coulombicExchangeArray

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

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

print "\nEigenvalues in eV:\n"

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

print eigenvalues2

print "\nDifference in Spectra:\n"

print array(eigenvalues2)-array(eigenvalues1)
