# Load Dynamical matrix
dynamical_matrix = nlread('../dyn-mat/ff-device-si.hdf5',DynamicalMatrix)[-1]

# Load Hamiltonian derivatives
hamiltonian_derivatives = nlread('ham-der.hdf5',HamiltonianDerivatives)[-1]

# Load DeviceConfiguration list
device_configuration_list = []
import glob
for filename in glob.glob("../0K-iv/ivcurve_selfconsistent_configurations_*.hdf5"):
    print('file read: ',filename)
    device_configuration_list.append(nlread(filename, DeviceConfiguration)[-1])

# Define energies; take only Fermi energy.
energies = [0.0]*eV

# Define q-points; take only points in first quadrant in (qx,qy) space
numQx = 1
numQy = 1

dQx = 0.5
dQy = 0.5

qx = numpy.linspace(0,dQx,numQx)
qy = numpy.linspace(0,dQy,numQy)

q_points = []
for x in qx:
    for y in qy:
        q_points.append([float(x),float(y),0.0])

# Define k-points; take only k=(0,0)
k_points = [[0.0,0.0,0.0]]

for i in range(0,len(device_configuration_list)):#range(1,len(device_configuration_list)):
    V_L = device_configuration_list[i].calculator().electrodeVoltages()[0].inUnitsOf(Volt)
    V_R = device_configuration_list[i].calculator().electrodeVoltages()[1].inUnitsOf(Volt)
    bias = V_L - V_R

    # -------------------------------------------------------------
    # Inelastic transmission spectrum
    # -------------------------------------------------------------
    inelastic_transmission_spectrum = InelasticTransmissionSpectrum(
        configuration=device_configuration_list[i],
        dynamical_matrix=dynamical_matrix,
        hamiltonian_derivatives=hamiltonian_derivatives,
        energies=energies,
        kpoints=k_points,
        qpoints=q_points,
        self_energy_calculator=RecursionSelfEnergy(),
        energy_zero_parameter=AverageFermiLevel,
        infinitesimal=1e-06*eV,
        phonon_modes=None,
        method=LOE,
    )

    filename = 'loe_bias_%4.4f_kpts_1x1_qpts_'%bias + str(numQx)+'x'+str(numQy)+'.hdf5'
    nlsave(filename,inelastic_transmission_spectrum)