In this tutorial we will simulate the operation of an organic donor/acceptor bulk heterojunction under irradation of light. Output will be the current voltage characteristics of the bulk heterojunction. To create the model of the bulk heterojunction we will use a morphology file and energy level file as can be e.g. created using cgkmc. Transfer integrals will be taken from a Quantum Patch calculation but could also be setup up parametrically.

Tutorial files:

The settings file can be created manually with the instructions given below or be downloaded here. The Quantum Patch input file can be downloaded here, the pdb files to identify donor and acceptor here. The morphology can be downloaded here, the energy level distribution here.


Make sure the environment varibales neccesary to run lightforge are set as shown here.

System setup

Mixed donor acceptor bulk hetero junction. Electrode workfunctions chosen so that built in voltage leads to charge separation.


Create a file called "settings" with the following content:

pbc: [False, True, True]
connect_electrodes: True
coulomb_mesh: True
mesh_scale: [1.0, 1.0, 1.0]
excitonics: "use presets"
electrode_stack_distance: 0.8

Periodic boundary conditions are switched off in transport direction as electrodes are attached. Exciton transport is treated by chosing from typical predefined types of molecules. Coulomb_mesh is required to treat long range periodic coulomb interaction, if two metallic contacts are present. The distance of electrodes and bhj is set to 0.8nm.

irradiance: 100000
light_on_time: 1.0e-2
light_off_time: 1.0e-2
rates: Miller

We choose an irradiance of 100000 w/m^2. This will result in the generation of excitons. Illumination will be switched on and off periodically every 0.01 s ( far longer than the simulation time). We use Miller-Abrahams type resonance condition for rate calculations to avoid the inverted Marcus regime for charge separation.

 holes: True
 electrons: True
 excitons: True

morphology_width: 50

We require all types of particles. The system has a width of 50 nm.

- name: molAmolB
  QP_output.zip: ICT_QP_output_0.zip

The Quantum Patch output with the transfer integrals.

- name: acceptor
  input_mode_transport: "QP: sig PAR: eaip,l"
    molecule_pdb: ICT_molecule_0.pdb
    QP_output_sigma: molAmolB
    exciton preset: fluorescent
    - [6.0,2.5]
    - [0.2,0.2]

- name: donor
  input_mode_transport: "QP: sig PAR: eaip,l"
    molecule_pdb: ICT_molecule_1.pdb
    QP_output_sigma: molAmolB
    exciton preset: fluorescent
    - [5.0,1.5]
    - [0.2,0.2]

We define two materials. A donor and an acceptor. Energy level distributions will be read in, but the mean will be shifted to a parametric EA and IP values. Reorganization energies are also supplied parametrically (0.2 eV in all cases). The molecule has to be identified with a pdb. QP_output_sigma will be ignored as energy level distributions will be read in from seperate file (see below). Excitonic properties of the material will be taken as typical properties for an fluorescent emitter. See here how to change exciton presets.

- thickness: 40
  morphology_input_mode: "custom Files: morphology, homolumo"
  morphology file: BHJ_COM.dat
  homolumo file: BJH_homolumo.dat
  - material: acceptor
    concentration: 1.0
  - material: donor
    concentration: 1.0

We define one mixed layer. Morphology and energy level distributions are supplied by file. The concentration field is overwritten by the custom morphology which contains molecule type information.

neighbours: 26
transfer_integral_source: QP_output
expansion_scheme: "no expansion (!dimensions fixed!)"
td_scaling: 0.4

Charge transport is possible to the 26 nearest neighbours. Transfer integrals are extrapolated using Quantum Patch output, the morphology is read in 1 to 1 and not increased in size. Make sure that morphology size of the read in file is equal to the specified width and thickness of the layer (here 50nm and 40 nm). Transition dipoles are estimated from exciton lifetimes as speciefied in the chosen exciton preset and the S1 energy. This estimatin is scaled by a factor 0.4. The transition dipoles are later used to calculate FRET rates. The lower transition dipole results in lower Förster radii, helping the effiency of the simulation.

- electrode_workfunction: -4.4
  coupling_model: QP_output

- electrode_workfunction: -3.1
  coupling_model: QP_output
std_scale: 1.0

Two electrodes are defined. The first electrode is attached before the layer that was defined first, the second one after the last layer. Transfer integrals to the electrode are estimated as being on average 0.8 standard deviations bigger than the distribution of organic organic transfer integrals. This leads to efficient charge ejection for the simulation.

- molecule 1: acceptor
  molecule 2: acceptor
  QP_output: molAmolB

- molecule 1: donor
  molecule 2: acceptor
  QP_output: molAmolB

- molecule 1: donor
  molecule 2: donor
  QP_output: molAmolB

We need to define the parameters for the transfer integrals between all kind of neighbouring pairs in the OLED. The QP output is defined at the top of the settings file. For each pair we need to supply a the name of QP_output file in which the transfer integrals can be found. In our case this is the output of a mixed QP calculation.

- simulations: 3
  measurement: pv_transient
  Temperature: 300
  field_direction: [1, 0, 0]
  field_strength: [-0.02,0.01,0.03,0.05]
  initial_holes: 0
  initial_electrons: 0

We do a simulation at constant Voltages. The field strength has to compensate for the in-built potential which is created due to the different workfunctions. 4.4V-3.1V = 1.3V. The length of the device is 40 nm, and electrodes are attached 0.8 nm from the layer. 1.3V/41.66nm = 0.03 V/nm. This means 0.03 V roughly compensates the inbuilt voltage and the current will change direction and the device will start to operate as a diode at larger fields (here 0.05 eV/nm). At the smallest field of -0.02 eV/nm energy is invested to pull charges out of the device. At 0.01 eV/nm and 0.03 eV/nm energy is extracted from the device. Thus we simulate all three quadrants of the operation.

iv_fluctuation: 0.001
new_wano: True
max_iterations: 1000000

We run the simulation for 1000000 Monte-Carlo steps

Simulation using lightforge

Run lightforge by typing:

# For V4
$OPENMPI_PATH/bin/mpirun -x OMP_NUM_THREADS --bind-to none -n 13 --mca btl self,vader,tcp python -m mpi4py $LFPATH/lightforge.py -s settings
# For V3 and below
$MPI_PATH/bin/mpirun -np 13 python -m mpi4py $LFPATH/lightforge.py -s settings

to simulate everything simultaneously:

Monitoring the calculation

In case of a parallel calculation the status of the simulation of each data point is written to the files output_job_i.

Inspecting the results

In the folder materials we get a overview of the energy levels in the device in the files energy_crossection_x.png The energy levels for the negative field of -0.02 eV to pull out charge: The energy levels at a field of 0.03 eV/nm close to the open circuit voltage:

The resulting IV looks like this:

We can see that the IV has the usual shape of an OPV device but a quite bad fill-factor.

The results of the search are