Anders this is the code i have used. check it and suggest some solution please .
from readTotalEnergy import readTotalEnergy
import pylab
# define the work function of gold
w = 5.28
#gate bias interval
v_g_interval = numpy.linspace(-10,10,201)
#source-drain bias interval
v_sd_interval = numpy.linspace(-15,15,301)
#----------finished define parameters ----------#
#read in the total energy data from the files: naphathalene_set[-2,-1,0,1,2].nc
voltage_list= []
energy_list = []
#loop over the charge states
charge_states = [-2,-1,0,1,2]
for q in charge_states:
voltage,energy = readTotalEnergy('naphthalene_set'+str(q)+'.nc')
voltage_list = voltage_list + [voltage]
# add energy of additional/missing electrons on benzene
energy = energy -q*w
energy_list = energy_list + [energy]
# make numpy arrays instead of lists
voltage_list = numpy.array(voltage_list)
energy_list = numpy.array(energy_list)
#make function that can return the charging energies
charging_energy = []
for i in range(len(charge_states)-1):
f = SplineInterpolation1D(voltage_list,energy_list[i+1]-energy_list)
charging_energy = charging_energy + [f]
#calculate number of charge states in the bias window
def conductionChannels(v_g,v_sd):
channels = 0
for f in charging_energy:
channels = channels + (abs(f(v_g)) <= abs(v_sd/2) )
return channels
#Generate the mesh points of the contour plot
X, Y = numpy.meshgrid(v_g_interval,v_sd_interval)
#evaluate number of charge states for each mesh point
Z = [ conductionChannels(X[i,j],Y[i,j])
for i in range(numpy.shape(X)[0])
for j in range(numpy.shape(X)[1])]
Z = numpy.array(Z).reshape(numpy.shape(X))
#make the plot
pylab.contour(X,Y,Z)
pylab.contourf(X,Y,Z)
pylab.xlabel("Gate voltage (Volt)")
pylab.ylabel("Source-Drain bias (Volt)")
pylab.show()
