import pylab as pl
import numpy as np

# Number of runs
num_run = 10

# Energy range for transmission
energies = numpy.linspace(-0.2, 0.2, 51)*eV
evec = energies.inUnitsOf(eV)

glist = []
for run in range(num_run):
    filename = 'nA_3-nB_3-nMD_3-T_300-run_%d.hdf5' % (run) 
    transmission = nlread(filename, TransmissionSpectrum)[0]
    t = transmission.evaluate()
    if  run == 0: t_av=t
    t_av=t_av+t
    pl.plot(t,evec,linewidth=1, color='k')
    g = transmission.conductance().inUnitsOf(Siemens)
    glist.append(g)

t_av = t_av/(num_run+1)
g_av = np.mean(glist)

device_config = nlread(filename, DeviceConfiguration)[0]
lele = device_config.electrodes()[0]
rele = device_config.electrodes()[1]
creg = device_config.centralRegion()
lattice = creg.bravaisLattice()
A = lattice.unitCellVolume()/lattice.primitiveVectors()[2][2]
L_lele = lele.bravaisLattice().primitiveVectors()[2,2].inUnitsOf(Ang)
L_rele = rele.bravaisLattice().primitiveVectors()[2,2].inUnitsOf(Ang)
L_creg = creg.bravaisLattice().primitiveVectors()[2,2].inUnitsOf(Ang)
L = L_creg - (L_lele + L_rele)

print()
print("+-------------------------------------------------------------------------+")
print("| MD-Landauer Reprot                                                      |")
print("+-------------------------------------------------------------------------+")
print("| Lenght of the MD region =", L,      "Ang")
print("| Averaged conductance    =", g_av,   "Siemens")
print("| Averaged resistance     =", 1/g_av, "Ohm")
print("+-------------------------------------------------------------------------+")

pl.plot(t_av,evec,linewidth=3, color='r', label='Average')
pl.tick_params(labelsize=14)
pl.legend(loc='lower right')
pl.ylim(-0.2,0.2)
pl.ylabel('Energy (eV)',fontsize=20)
pl.xlabel('Transmission coefficient',fontsize=20)
pl.savefig('transmission.png')
pl.show()