Script used for calculation is mentioned below
read Total Energy script.
---
from NL.IO.NLSaveUtilities import nlinspect
from NanoLanguage import *
#define function for reading in total energy and gate voltage data from a file
def readTotalEnergy(filename):
"""
function for reading a list of total energies and gate voltages from a file
@param filename : filename of the netcdf file
@return voltages, energies : list of gate voltages and energies
found in filename
"""
if filename == None:
return
#get list of data in file
file_data = nlinspect(filename)
if file_data == None:
raise ValueError, "Wrong file format"
total_energy_list = numpy.array([])
gate_voltage_list = numpy.array([])
for id in file_data:
if (id[0] == 'TotalEnergy'):
total_energy=nlread(filename, object_id = id[1])[0]
total_energy_list = numpy.append(total_energy_list,
total_energy.evaluate().inUnitsOf(eV))
#extract the gate voltage from the id of the data
if id[1][0:3] == 'gID':
gate_voltage = 0.0
else:
gate_voltage = float(id[1].strip(
"abcdefghijklmnopqrstuvxyzABCDEFGHIJKLMNOPQRSTUVXYZ"))
gate_voltage_list = numpy.append(gate_voltage_list,gate_voltage)
if len(total_energy_list) ==0:
raise ValueError, "No total energy in NC file"
#sort the data
index_list = numpy.argsort(gate_voltage_list)
return gate_voltage_list[index_list], total_energy_list[index_list]
------------
Total Energy VS gate voltage script
-------
from readTotalEnergy import readTotalEnergy
# define the work function of gold
w = 5.28
#read in the total energy data from the files: pentacene_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('pentacene_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]
#plot the total energies
import pylab
pylab.figure()
for i in range(len(voltage_list)):
pylab.plot(voltage_list,energy_list)
pylab.xlabel("Gate voltage (Volt)")
pylab.ylabel("Total Energy (eV)")
pylab.show()
---
I was facing problem in attaching script file that's why copying script in reply tab.
Thanking you!
