Add steady state charge transport simulation

[ljakoster]

Would be great if the code could also do periodic boundary conditions in all three directions to calculate mobility in steady-state. That plus output of the Fermi energy and the transport energy would make it suitable for organic thermoelectrics.

[MikeHeiber]

@ljakoster Let me gather some more details about how you would want to do this simulation. Here's my initial thoughts for how it would work. Let me know if I am missing anything that you'd like to have.

It would be a simulation box without contacts, and charges would just continually flow through the volume with 3D periodic boundaries. It would be a single carrier simulation, and I would have the user define the carrier density. To make things easy, I would propose to not implement a test for steady state convergence, and instead the users could run the simulation for a different amounts of time to check for convergence, so then the user would also specify the test time. I would output the average occupation energy. The Fermi level should be calculable analytically based on whatever is the chosen DOS, but I could also output the actual value from the lattice, so that one could compare the two. Determining the transport energy would potentially be more complicated, depending on your definition. Do you have any suggestions for how to numerically calculate the transport energy? Any other capabilities you would want added?

[ljakoster]

That sounds great! That is exactly what I was trying to say. Just a simulation box without contacts but there is an electric field defined by the user. You're right about steady-state.

I think it is important to calculate the Fermi energy based on the actual realisation of whatever the energetic landscape is, especially if you implement reading the energies from a file (I think that was another feature request and it would be very useful). That would also allow the user to implement ionized dopants as fixed charges that sit off-grid. Their Coulomb wells can then be included in the energetic landscape and their effect is then automatically included in both the Fermi energy and the transport energy. That is probably much less difficult than actually implementing a set of counter ions (i.e. the ionized dopants) and update the Coulomb potentials while running the simulation. To study doped OSCs you really need to include the counter charges and it would be easiest, I think, to do it via the initial site energies.

About the transport energy: The most suited definition is: E_T = \int E \frac{\sigma(E)}{\sigma} dE, where E is the hop energy, \sigma is the total conductivity (or mobility or even current, it doesn't matter), and \sigma(E) is the conductivity flowing at a specific hop energy. The integral is then calculated over the entire volume. There exist three ways of taking E from the initial and final site energies in a jump: E=initial, E=final, and E=average of the two. It is customary to take the average. We have tried all three other simulations and the differences are, luckily, small.

What I do in my master equation code is the following (pascal!): begin trans_en:=0; {the transport energy} nrm:=0; {this is used to normalise to the total current} {loop over all grid points:} for i:=1 to nx do {this is the x-direction, the direction of the applied electric field} for j:=1 to ny do for k:=1 to nz do begin ip:=gx(i+1); {gx is a function that takes care of the PBC in x} {so away i=1 and i=nx, ip is just i+1} {Rxmin: hop rate in x-direction from i+1 -> i} {Rxpls: hop rate in x-direction from i -> i+1} {n[i,j,k]: electron density at grid point (i,j,k)} cin:=Rxmin[i,j,k]n[ip,j,k](1-n[i,j,k]); {current in} cout:=Rxpls[i,j,k]n[i,j,k](1-n[ip,j,k]); {current out} en:=0.5(eps[i,j,k] + eps[ip,j,k]); {energy of hop taken as average of both sites} trans_en:=trans_en + en(cout-cin); nrm:=nrm + cout - cin; end; if nrm <> 0 then trans_en:=trans_en/nrm else trans_en:=0; {in the end we divide by the total current in x} end; The electric field is applied in the x-direction, so the only net current flows in x. I therefore don't include the hops in y and z.

[MikeHeiber]

@ljakoster I like the idea of including fixed ionized dopants through the site energy import feature. Site energy import will definitely be a part of v1.0. I just finished adding it to the development branch this week. So yes, with this in place one would not be able to calculate the Fermi level analytically, and the actual Fermi level of the specific simulation implemented would be nice to have.

I am still not sure how I want to implement the transport energy calculation. I was also reading some ideas about how to determine it from the Baranovskii group in this paper. Maybe I can get Oelerich's opinion. @janoliver , would be willing to provide your thoughts on the transport energy calculation?

Another idea I have is to record the average of the initial and final state energies and the particle velocity in the field direction for every hop and then take a weighted average where the average energy is weighted by the particle velocity. In the end you would get an average energy that is biased towards the energies involved in the actual long range transport of the charge carriers in the field direction. Does this seem reasonable?

Anyway, if we can finalize how to best do the transport energy calculation soon, I would like implement this steady state charge transport feature in v1.0.

[ljakoster]

dear Mike,

I really like the paper you refer to but I think their aim is a bit different in that they tried to show that there is such a thing as a unique transport energy. To calculate the Seebeck coefficient, one needs to evaluate that integral. I think what you are proposing about recording the speed in the field direction is correct, but I will get back to you shortly.

Jan Anton

-- Prof. Dr. L.J.A. (Jan Anton) Koster Zernike Institute for Advanced Materials University of Groningen website http://www.koster-group.nl

On 9 November 2018 at 16:33, Michael C. Heiber notifications@github.com wrote:

I like the idea of including fixed ionized dopants through the site energy import feature. Site energy import will definitely be a part of v1.0. I just finished adding it to the development branch this week. So yes, with this in place one would not be able to calculate the Fermi level analytically, and the actual Fermi level of the specific simulation implemented would be nice to have.

I am still not sure how I want to implement the transport energy calculation. I was also reading some ideas about how to determine it from the Baranovskii group in this paper https://doi.org/10.1088/0953-8984/26/25/255801. Maybe I can get Oelerich's opinion. @janoliver https://github.com/janoliver , would be willing to provide your thoughts on the transport energy calculation?

Another idea I have is to record the average of the initial and final state energies and the particle velocity in the field direction for every hop and then take a weighted average where the average energy is weighted by the particle velocity. In the end you would get an average energy that is biased towards the energies involved in the actual long range transport of the charge carriers in the field direction. Does this seem reasonable?

Anyway, if we can finalize how to best do the transport energy calculation soon, I would like implement this steady state charge transport feature in v1.0.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/MikeHeiber/Excimontec/issues/60#issuecomment-437396073, or mute the thread https://github.com/notifications/unsubscribe-auth/AkrGVBKIt0Lm1x8cCVd1g_QVuSNYukeQks5utaBkgaJpZM4YOK8q .

[janoliver]

I am still not sure how I want to implement the transport energy calculation. I was also reading some ideas about how to determine it from the Baranovskii group in this paper. Maybe I can get Oelerich's opinion. @janoliver , would be willing to provide your thoughts on the transport energy calculation?

If you need it, I can provide a C code based on GSL to calculate the transport energy as defined in our paper. (It would be a self-consistent solution of Eq. (13) therein).

[ljakoster]

I still think it's best to directly evaluate the integral form as we want to get the Seebeck coefficient from the Peltier coefficient. This is not necessarily the same as the transport energy concept as used in analytical theories for charge transport (don't get me wrong, I really like the analytical approach!), but people call it transport energy nonetheless. Your suggestion to record the average of the initial and final state energies and the particle velocity in the field direction for every hop and then take a weighted average where the average energy is weighted by the particle velocity sounds exactly right, but I am not sure how you calculate the velocity? Could you explain how you calculate the current in the field direction in KMC? the hops are instantaneous, right? so the time that enters into the velocity is the time between the hops?

jan anton

-- Prof. Dr. L.J.A. (Jan Anton) Koster Zernike Institute for Advanced Materials University of Groningen website http://www.koster-group.nl

On Wed, 14 Nov 2018 at 09:14, Jan Oliver Oelerich notifications@github.com wrote:

I am still not sure how I want to implement the transport energy calculation. I was also reading some ideas about how to determine it from the Baranovskii group in this paper https://doi.org/10.1088/0953-8984/26/25/255801. Maybe I can get Oelerich's opinion. @janoliver https://github.com/janoliver , would be willing to provide your thoughts on the transport energy calculation?

If you need it, I can provide a C code based on GSL to calculate the transport energy as defined in our paper. (It would be a self-consistent solution of Eq. (13) therein).

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/MikeHeiber/Excimontec/issues/60#issuecomment-438575032, or mute the thread https://github.com/notifications/unsubscribe-auth/AkrGVEk2RQJyMml8GiV_HjEr6oWRNv1oks5uu9DkgaJpZM4YOK8q .

[MikeHeiber]

@janoliver Do you have a suggestion for calculating the transport energy numerically from the statistics of a KMC simulation. In your paper, you do this by selectively removing energies and re-running the simulations, but this method requires many transport simulations. I would prefer to determine an approximation during a single simulation. What do you think about the weighted average method that I describe above? Or maybe you have a better method in mind?

@ljakoster I was thinking that I would calculate the particle velocity using the time interval between hops. I would take the hop distance in the field direction and then divide by the residence time. Shallow traps would then have a short residence time and high velocity even if the hop distance was only 1 lattice unit. Also, with a VRH model, you could also have slower long range hops that also contribute. This is similar to what I do during the ToF simulation to numerically calculate the current over a given time interval. Does this sound like a good method?

[ljakoster]

Yes, I think that is a solid method.

Prof. Dr. L.J.A. (Jan Anton) Koster Zernike Institute for Advanced Materials University of Groningen website http://www.koster-group.nl

On Wed, 14 Nov 2018 at 16:24, Michael C. Heiber notifications@github.com wrote:

@janoliver https://github.com/janoliver Do you have a suggestion for calculating the transport energy numerically from the statistics of a KMC simulation. In your paper, you do this by selectively removing energies and re-running the simulations, but this method requires many transport simulations. I would prefer to determine an approximation during a single simulation. What do you think about the weighted average method that I describe above? Or maybe you have a better method in mind?

@ljakoster https://github.com/ljakoster I was thinking that I would calculate the particle velocity using the time interval between hops. I would take the hop distance in the field direction and then divide by the residence time. Shallow traps would then have a short residence time and high velocity even if the hop distance was only 1 lattice unit. Also, with a VRH model, you could also have slower long range hops that also contribute. This is similar to what I do during the ToF simulation to numerically calculate the current over a given time interval. Does this sound like a good method?

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/MikeHeiber/Excimontec/issues/60#issuecomment-438701050, or mute the thread https://github.com/notifications/unsubscribe-auth/AkrGVHLc8nlOynY-Za3K2PP6k81-RLC7ks5uvDXGgaJpZM4YOK8q .

[MikeHeiber]

@ljakoster The steady state charge transport feature is now implemented on the development branch and should be fully functional. I have put in a verification test for transport in a Gaussian DOS at low carrier density and low E-field, and the equilibration energy is reproduced nicely at -sigma^2/(kT) as expected analytically. Everything will soon be merged into the master branch and become part of the next major release. Let me know if you get a chance to try out the steady transport test feature and have any questions or further suggestions.

[MikeHeiber]

Resolved in #67

[ljakoster]

dear Mike, that's great, thanks! Will try it out soon. Have you also managed to implement the Fermi and transport energies?

best wishes, Jan Anton

Prof. Dr. L.J.A. (Jan Anton) Koster Zernike Institute for Advanced Materials University of Groningen website http://www.koster-group.nl

On Tue, 11 Dec 2018 at 19:13, Michael C. Heiber notifications@github.com wrote:

Closed #60 https://github.com/MikeHeiber/Excimontec/issues/60.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/MikeHeiber/Excimontec/issues/60#event-2019469973, or mute the thread https://github.com/notifications/unsubscribe-auth/AkrGVLLzZJArexxQLXkwvTHRQV5efRCEks5u3_W3gaJpZM4YOK8q .

[MikeHeiber]

@ljakoster Yes everything is implemented. I also allow the user to specify the number of equilibration iterations before the simulation starts collecting data to calculate the steady state results. The simulation creates the appropriate number of polarons in the lattice based on the specified charge carrier density by filling up the DOS, so basically it starts at a 0K configuration. It then starts running the simulation for the specified number of equilibration iterations. Finally it continues running for another set of iterations specified by the N_tests parameter and during those iterations it calculates the mobility, the transport energy, and the equilibration energy (avg DOOS). The final results file contains all of this information and also the numerically determined Fermi energy.

[ljakoster]

wonderful! I have managed to get everything up and running (took a while to figure out how to set up googletest and KMC_Lattice, but it's easy really) and I am looking forward to doing some tests!

Prof. Dr. L.J.A. (Jan Anton) Koster Zernike Institute for Advanced Materials University of Groningen website http://www.koster-group.nl

On Wed, 12 Dec 2018 at 17:00, Michael C. Heiber notifications@github.com wrote:

@ljakoster https://github.com/ljakoster Yes everything is implemented. I also allow the user to specify the number of equilibration iterations before the simulation starts collecting data to calculate the steady state results. The simulation creates the appropriate number of polarons in the lattice based on the specified charge carrier density by filling up the DOS, so basically it starts at a 0K configuration. It then starts running the simulation for the specified number of equilibration iterations. Finally it continues running for another set of iterations specified by the N_tests parameter and during those iterations it calculates the mobility, the transport energy, and the equilibration energy (avg DOOS). The final results file contains all of this information and also the numerically determined Fermi energy.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/MikeHeiber/Excimontec/issues/60#issuecomment-446639846, or mute the thread https://github.com/notifications/unsubscribe-auth/AkrGVLta8d_eibnCMH0eGXOvgVngT8LKks5u4SgAgaJpZM4YOK8q .

[ljakoster]

hi Mike,

thanks again. I have the code up and running. I have a few questions though: 1) the mobility, EF, Etransport; at which time are they displayed? or is this a time-average? the system needs to relax/thermalise first, right? How can I check if the system is sufficiently thermalised? I would have thought that I needed to track mobility and the energies over time? 2) the file 'analysis_summary.txt' lists the mobility and the Fermi energy, but not the transport and equilibration energies (nan +- nan)?

jan anton

-- Prof. Dr. L.J.A. (Jan Anton) Koster Zernike Institute for Advanced Materials University of Groningen website http://www.koster-group.nl

On Wed, 12 Dec 2018 at 17:00, Michael C. Heiber notifications@github.com wrote:

@ljakoster https://github.com/ljakoster Yes everything is implemented. I also allow the user to specify the number of equilibration iterations before the simulation starts collecting data to calculate the steady state results. The simulation creates the appropriate number of polarons in the lattice based on the specified charge carrier density by filling up the DOS, so basically it starts at a 0K configuration. It then starts running the simulation for the specified number of equilibration iterations. Finally it continues running for another set of iterations specified by the N_tests parameter and during those iterations it calculates the mobility, the transport energy, and the equilibration energy (avg DOOS). The final results file contains all of this information and also the numerically determined Fermi energy.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/MikeHeiber/Excimontec/issues/60#issuecomment-446639846, or mute the thread https://github.com/notifications/unsubscribe-auth/AkrGVLta8d_eibnCMH0eGXOvgVngT8LKks5u4SgAgaJpZM4YOK8q .

[MikeHeiber]

@ljakoster Once all of the polarons are created the simulation runs for a number of equilibration iterations. You can modify how many iterations are used for the equilibration stage by changing the N_equilibration_events parameter. After the equilibration stage is done, the simulation enters the test stage and will run for another number iterations set by the N_tests parameter. The results are then a time average over the test stage. You can increase the N_tests number to get better statistics for the final averages. If you want to check whether the system is thermalized, I would keep N_tests constant and change N_equilibration_events and then see if the results change.

I am not sure why the analysis_summary file is showing nan for the transport and equilibration energies. What are you using for the DOS and for the Steady_carrier_density, N_equilibration_events and N_tests parameters?

[ljakoster]

ah, thanks, that makes sense. I thought that N_equilibration_events set the overall simulation time and I didn't understand N_tests. I have increased the latter (from 10) and decreased the former. I will play with these.

About the analysis file: it's a bit weird as the resultsXXX.txt files do show the EF and E_transport. Looks like these numbers are not read properly?

-- Prof. Dr. L.J.A. (Jan Anton) Koster Zernike Institute for Advanced Materials University of Groningen website http://www.koster-group.nl

On Thu, 13 Dec 2018 at 16:01, Michael C. Heiber notifications@github.com wrote:

@ljakoster https://github.com/ljakoster Once all of the polarons are created the simulation runs for a number of equilibration iterations. You can modify how many iterations are used for the equilibration stage by changing the N_equilibration_events parameter. After the equilibration stage is done, the simulation enters the test stage and will run for another number iterations set by the N_tests parameter. The results are then a time average over the test stage. You can increase the N_tests number to get better statistics for the final averages. If you want to check whether the system is thermalized, I would keep N_tests constant and change N_equilibration_events and then see if the results change.

I am not sure why the analysis_summary file is showing nan for the transport and equilibration energies. What are you using for the DOS and for the Steady_carrier_density, N_equilibration_events and N_tests parameters?

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/MikeHeiber/Excimontec/issues/60#issuecomment-446999144, or mute the thread https://github.com/notifications/unsubscribe-auth/AkrGVK12dl0WdUAC09NWnHej04QblIUFks5u4mvHgaJpZM4YOK8q .

[ljakoster]

with the more sensible settings (N_tests), the analysis file does display EF and Etr. Wonderful! Will run more tests to make sure I am using the code correctly.

Prof. Dr. L.J.A. (Jan Anton) Koster Zernike Institute for Advanced Materials University of Groningen website http://www.koster-group.nl

On Thu, 13 Dec 2018 at 16:39, Koster, L.J.A. l.j.a.koster@rug.nl wrote:

ah, thanks, that makes sense. I thought that N_equilibration_events set the overall simulation time and I didn't understand N_tests. I have increased the latter (from 10) and decreased the former. I will play with these.

About the analysis file: it's a bit weird as the resultsXXX.txt files do show the EF and E_transport. Looks like these numbers are not read properly?

-- Prof. Dr. L.J.A. (Jan Anton) Koster Zernike Institute for Advanced Materials University of Groningen website http://www.koster-group.nl

On Thu, 13 Dec 2018 at 16:01, Michael C. Heiber notifications@github.com wrote:

@ljakoster https://github.com/ljakoster Once all of the polarons are created the simulation runs for a number of equilibration iterations. You can modify how many iterations are used for the equilibration stage by changing the N_equilibration_events parameter. After the equilibration stage is done, the simulation enters the test stage and will run for another number iterations set by the N_tests parameter. The results are then a time average over the test stage. You can increase the N_tests number to get better statistics for the final averages. If you want to check whether the system is thermalized, I would keep N_tests constant and change N_equilibration_events and then see if the results change.

I am not sure why the analysis_summary file is showing nan for the transport and equilibration energies. What are you using for the DOS and for the Steady_carrier_density, N_equilibration_events and N_tests parameters?

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/MikeHeiber/Excimontec/issues/60#issuecomment-446999144, or mute the thread https://github.com/notifications/unsubscribe-auth/AkrGVK12dl0WdUAC09NWnHej04QblIUFks5u4mvHgaJpZM4YOK8q .