# This is a 1D simulation according to the model described in [1]. 
# The 'createfigures.gnuplot' file will create similar plots as the
# ones that were shown in the article
#
# [1] L. J. A. Koster, E. C. P. Smits, V. D. Mihailetchi and P. W. M. Blom,
#     "Device model for the operation of polymer/fullerene bulk heterojunction 
#     solar cells", Phys. Rev. B 72, 085205, 2005
#     http://journals.aps.org/prb/abstract/10.1103/PhysRevB.72.085205

# These come from Table I in [1]
math/def Egap 1.34
math/def eMob 2.5e-7
math/def hMob 3.0e-8
math/def Nc 2.5e25
math/def G 2.7e27
math/def epsilonTotal 3.0e-11
math/def epsRel epsilonTotal/eps0
math/def a 1.3e-9
math/def kf 1.5e6

math/def T 300 		# Not sure what this should be
math/def kT_ev kB*T/e
math/def Vbi Egap	# Just another name for the same value

# According to the article, the prefactor in the bimolecular recombination
# is the Langevin one.
math/def gamma e*0.5*(eMob+hMob)/epsilonTotal

# Einstein relation fixes the diffusion constants
math/def De eMob*kT_ev
math/def Dh hMob*kT_ev

# In the article a device of width 120 nm is used. For the simulation we'll
# use 120 grid points
math/def numPixels 120
math/def width 120e-9

# Start the simulation with the specified number of pixels and physical
# width. Also set the temperature parameter.
sim1/new numPixels width
sim1/temp/set T

# Let's define names for the left and right pixels (the contacts)
reg1/new leftContact
reg1/append/line leftContact 1 1

reg1/new rightContact
reg1/append/line rightContact NX1 NX1

# Boundary conditions as described in the article
sim1/reg/set n leftContact Nc*exp(-Egap/kT_ev)
sim1/reg/set n rightContact Nc
sim1/reg/set p leftContact Nc
sim1/reg/set p rightContact Nc*exp(-Egap/kT_ev)

# Set some other properties of the device
sim1/reg/set g ALL G
sim1/reg/set rf ALL gamma
sim1/reg/set dn ALL De
sim1/reg/set dp ALL Dh
sim1/reg/set nmob ALL eMob
sim1/reg/set pmob ALL hMob
sim1/reg/set eps ALL epsRel

# Here we specify that we're going to use the non-default recombination model
# (the one from the article) with the specified parameters.
sim1/rec/ext a kf

# Set the voltage drop over the device, this corresponds to an applied voltage
# of 0 V. We'll also initialize n and p in the interior of the device using
# interpolation.
sim1/phi/set Vbi yes
sim1/grid/init n
sim1/grid/init p

# Search for the equilibrium situation using the Newton-Raphson solver, and
# plot the data tot the file 'data_oc.dat'
sim1/rundirect
sim1/grid/fplot data_oc.dat no

# Determine the JV curve, and write it to the file 'jv.txt'
sim1/iv Vbi 0 0.9 100 jv.txt no

# Then we'll set the applied voltage to the open circuit voltage 
# (determined in a separate run by the previous command), we'll look for
# the equilibrium situation, and write the result to 'data_sc.dat'
math/def Voc 0.84318 # Obtained from a first simulation
sim1/phi/set Vbi-Voc yes
sim1/rundirect
sim1/grid/fplot data_sc.dat no