from scipy import *
from numpy import *

sumModesN = 16

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**-3:
                target[i][j] = float(dummyFinal[m])
            else:
                target[i][j] = 0.0
            m += 1
    return array(target)


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

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

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

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


realTotalArray = potentialArray + kineticArray + coulombicRepulsionArray + coulombicExchangeArray
print "\nHamiltonian Array\n"
print realTotalArray
output = linalg.eig(realTotalArray)

print "\nEigenvalues:\n"

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

print eigenvalues
##print "\nEigenvalues in eV:\n"
##print output[0]*27.2113838685
