Bugfix DOS integral (#88)

General: - Updated Changelog to reflect all changes - Updated parameters_default file version string to v1.0.0-rc.3 - Updated Doxygen documentation main: - Corrected DOS and DOOS data output to have the correct units - Updated version string to v1.0.0-rc.3 OSC_Sim class: - Added DOS and DOOS sampling counters to keep track of how many times the distributions have been sampled - Fixed getSteadyDOS and getSteadyDOOS functions to use new sampling counters and avoid errors averaging over multiple samplings - Changed updateSteadyData to increment the new sampling counters and removed calculation of the DOS (without Coulomb) because it does not change during the simulation tests: - Updated SteadyTransportTests by adding tests of the integral of the DOS and DOOS data
MikeHeiber · Apr 2, 2019 · dc3ba36 · dc3ba36
1 parent 886551c
commit dc3ba36
Show file tree

Hide file tree

Showing 16 changed files with 101 additions and 52 deletions.
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -6,6 +6,28 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 --------------------------------------------------------------------------------------------------------------------------------
 
+## [v1.0.0-rc.3]- 2019-04-01 - Density of States Integration Bugfix
+
+### Added
+- OSC_Sim - Steady_DOS_sampling_count and Steady_DOOS_sampling_counter to keep track of how many times the DOS and DOOS are sampled during the simulation
+- test.cpp (SteadyTransportTests) - tests to check the integral of the DOS and DOOS with and without including Coulomb interactions
+
+### Changed
+- main.cpp - Version string to v1.0.0-rc.3
+- OSC_Sim (updateSteadyData) - To increment the new DOS and DOOS sampling counters and removed calculation of the DOS because the DOS is does not change during the simulation when not including Coulomb interactions
+- OSC_Sim (getSteadyDOS) - Function is not longer const
+- parameters_default.txt - Version string to v1.0.0-rc.3
+
+### Removed
+
+### Fixed
+- main.cpp - Output of density of states and density of occupied states data to have the correct units
+- OSC_Sim (getSteadyDOOS) - To use the new sampling counter to avoid errors averaging over multiple samplings
+- OSC_Sim (getSteadyDOOS_Coulomb) - To use the new sampling counter to avoid errors averaging over multiple samplings
+- OSC_Sim (getSteadyDOS) - To calculate the DOS using only one sample at call time because the DOS (without Coulomb) does not change during the simulation
+- OSC_Sim (getSteadyDOS_Coulomb) - To use the new sampling counter to avoid errors averaging over multiple samplings
+
+
 ## [v1.0.0-rc.2]- 2019-02-09 - Steady State Charge Transport Test Update
 
 ### Added

diff --git a/docs/_o_s_c___sim_8h_source.html b/docs/_o_s_c___sim_8h_source.html
diff --git a/docs/class_excimontec_1_1_o_s_c___sim-members.html b/docs/class_excimontec_1_1_o_s_c___sim-members.html
@@ -165,7 +165,7 @@
   <tr class="even"><td class="entry"><a class="el" href="class_excimontec_1_1_o_s_c___sim.html#a1e2d85af55f9898fc4beaea3c09b2ab2">getSteadyCurrentDensity</a>() const</td><td class="entry"><a class="el" href="class_excimontec_1_1_o_s_c___sim.html">Excimontec::OSC_Sim</a></td><td class="entry"></td></tr>
   <tr><td class="entry"><a class="el" href="class_excimontec_1_1_o_s_c___sim.html#a502194201bd4e045760149bd006112c4">getSteadyDOOS</a>() const</td><td class="entry"><a class="el" href="class_excimontec_1_1_o_s_c___sim.html">Excimontec::OSC_Sim</a></td><td class="entry"></td></tr>
   <tr class="even"><td class="entry"><a class="el" href="class_excimontec_1_1_o_s_c___sim.html#a7bba196811f642f683570240ceb77e62">getSteadyDOOS_Coulomb</a>() const</td><td class="entry"><a class="el" href="class_excimontec_1_1_o_s_c___sim.html">Excimontec::OSC_Sim</a></td><td class="entry"></td></tr>
-  <tr><td class="entry"><a class="el" href="class_excimontec_1_1_o_s_c___sim.html#a6440b290be26622c04ea89cf7293575b">getSteadyDOS</a>() const</td><td class="entry"><a class="el" href="class_excimontec_1_1_o_s_c___sim.html">Excimontec::OSC_Sim</a></td><td class="entry"></td></tr>
+  <tr><td class="entry"><a class="el" href="class_excimontec_1_1_o_s_c___sim.html#a15e362e4b5b2c065e29f1cdf36235784">getSteadyDOS</a>()</td><td class="entry"><a class="el" href="class_excimontec_1_1_o_s_c___sim.html">Excimontec::OSC_Sim</a></td><td class="entry"></td></tr>
   <tr class="even"><td class="entry"><a class="el" href="class_excimontec_1_1_o_s_c___sim.html#ad819d76bf42693993bb392c48327ebf6">getSteadyDOS_Coulomb</a>() const</td><td class="entry"><a class="el" href="class_excimontec_1_1_o_s_c___sim.html">Excimontec::OSC_Sim</a></td><td class="entry"></td></tr>
   <tr><td class="entry"><a class="el" href="class_excimontec_1_1_o_s_c___sim.html#a451e3ea3667f84d6841ad91355a618df">getSteadyEquilibrationEnergy</a>() const</td><td class="entry"><a class="el" href="class_excimontec_1_1_o_s_c___sim.html">Excimontec::OSC_Sim</a></td><td class="entry"></td></tr>
   <tr class="even"><td class="entry"><a class="el" href="class_excimontec_1_1_o_s_c___sim.html#ab726ddf220bbef11061b802df50c0188">getSteadyEquilibrationEnergy_Coulomb</a>() const</td><td class="entry"><a class="el" href="class_excimontec_1_1_o_s_c___sim.html">Excimontec::OSC_Sim</a></td><td class="entry"></td></tr>

diff --git a/docs/class_excimontec_1_1_o_s_c___sim.html b/docs/class_excimontec_1_1_o_s_c___sim.html
@@ -282,9 +282,9 @@
 <tr class="memitem:a7bba196811f642f683570240ceb77e62"><td class="memItemLeft" align="right" valign="top">std::vector&lt; std::pair&lt; double, double &gt; &gt;&#160;</td><td class="memItemRight" valign="bottom"><a class="el" href="class_excimontec_1_1_o_s_c___sim.html#a7bba196811f642f683570240ceb77e62">getSteadyDOOS_Coulomb</a> () const</td></tr>
 <tr class="memdesc:a7bba196811f642f683570240ceb77e62"><td class="mdescLeft">&#160;</td><td class="mdescRight">Gets the density of occupied states calculated during the steady state charge transport test.  <a href="#a7bba196811f642f683570240ceb77e62">More...</a><br /></td></tr>
 <tr class="separator:a7bba196811f642f683570240ceb77e62"><td class="memSeparator" colspan="2">&#160;</td></tr>
-<tr class="memitem:a6440b290be26622c04ea89cf7293575b"><td class="memItemLeft" align="right" valign="top">std::vector&lt; std::pair&lt; double, double &gt; &gt;&#160;</td><td class="memItemRight" valign="bottom"><a class="el" href="class_excimontec_1_1_o_s_c___sim.html#a6440b290be26622c04ea89cf7293575b">getSteadyDOS</a> () const</td></tr>
-<tr class="memdesc:a6440b290be26622c04ea89cf7293575b"><td class="mdescLeft">&#160;</td><td class="mdescRight">Gets the density of states calculated during the steady state charge transport test.  <a href="#a6440b290be26622c04ea89cf7293575b">More...</a><br /></td></tr>
-<tr class="separator:a6440b290be26622c04ea89cf7293575b"><td class="memSeparator" colspan="2">&#160;</td></tr>
+<tr class="memitem:a15e362e4b5b2c065e29f1cdf36235784"><td class="memItemLeft" align="right" valign="top">std::vector&lt; std::pair&lt; double, double &gt; &gt;&#160;</td><td class="memItemRight" valign="bottom"><a class="el" href="class_excimontec_1_1_o_s_c___sim.html#a15e362e4b5b2c065e29f1cdf36235784">getSteadyDOS</a> ()</td></tr>
+<tr class="memdesc:a15e362e4b5b2c065e29f1cdf36235784"><td class="mdescLeft">&#160;</td><td class="mdescRight">Gets the density of states calculated during the steady state charge transport test.  <a href="#a15e362e4b5b2c065e29f1cdf36235784">More...</a><br /></td></tr>
+<tr class="separator:a15e362e4b5b2c065e29f1cdf36235784"><td class="memSeparator" colspan="2">&#160;</td></tr>
 <tr class="memitem:ad819d76bf42693993bb392c48327ebf6"><td class="memItemLeft" align="right" valign="top">std::vector&lt; std::pair&lt; double, double &gt; &gt;&#160;</td><td class="memItemRight" valign="bottom"><a class="el" href="class_excimontec_1_1_o_s_c___sim.html#ad819d76bf42693993bb392c48327ebf6">getSteadyDOS_Coulomb</a> () const</td></tr>
 <tr class="memdesc:ad819d76bf42693993bb392c48327ebf6"><td class="mdescLeft">&#160;</td><td class="mdescRight">Gets the density of states calculated during the steady state charge transport test.  <a href="#ad819d76bf42693993bb392c48327ebf6">More...</a><br /></td></tr>
 <tr class="separator:ad819d76bf42693993bb392c48327ebf6"><td class="memSeparator" colspan="2">&#160;</td></tr>
@@ -1747,8 +1747,8 @@ <h2 class="memtitle"><span class="permalink"><a href="#a7bba196811f642f683570240
 
 </div>
 </div>
-<a id="a6440b290be26622c04ea89cf7293575b"></a>
-<h2 class="memtitle"><span class="permalink"><a href="#a6440b290be26622c04ea89cf7293575b">&#9670;&nbsp;</a></span>getSteadyDOS()</h2>
+<a id="a15e362e4b5b2c065e29f1cdf36235784"></a>
+<h2 class="memtitle"><span class="permalink"><a href="#a15e362e4b5b2c065e29f1cdf36235784">&#9670;&nbsp;</a></span>getSteadyDOS()</h2>
 
 <div class="memitem">
 <div class="memproto">
@@ -1757,7 +1757,7 @@ <h2 class="memtitle"><span class="permalink"><a href="#a6440b290be26622c04ea89cf
           <td class="memname">vector&lt; std::pair&lt; double, double &gt; &gt; Excimontec::OSC_Sim::getSteadyDOS </td>
           <td>(</td>
           <td class="paramname"></td><td>)</td>
-          <td> const</td>
+          <td></td>
         </tr>
       </table>
 </div><div class="memdoc">

diff --git a/docs/class_excimontec_1_1_o_s_c___sim.js b/docs/class_excimontec_1_1_o_s_c___sim.js
@@ -58,7 +58,7 @@ var class_excimontec_1_1_o_s_c___sim =
     [ "getSteadyCurrentDensity", "class_excimontec_1_1_o_s_c___sim.html#a1e2d85af55f9898fc4beaea3c09b2ab2", null ],
     [ "getSteadyDOOS", "class_excimontec_1_1_o_s_c___sim.html#a502194201bd4e045760149bd006112c4", null ],
     [ "getSteadyDOOS_Coulomb", "class_excimontec_1_1_o_s_c___sim.html#a7bba196811f642f683570240ceb77e62", null ],
-    [ "getSteadyDOS", "class_excimontec_1_1_o_s_c___sim.html#a6440b290be26622c04ea89cf7293575b", null ],
+    [ "getSteadyDOS", "class_excimontec_1_1_o_s_c___sim.html#a15e362e4b5b2c065e29f1cdf36235784", null ],
     [ "getSteadyDOS_Coulomb", "class_excimontec_1_1_o_s_c___sim.html#ad819d76bf42693993bb392c48327ebf6", null ],
     [ "getSteadyEquilibrationEnergy", "class_excimontec_1_1_o_s_c___sim.html#a451e3ea3667f84d6841ad91355a618df", null ],
     [ "getSteadyEquilibrationEnergy_Coulomb", "class_excimontec_1_1_o_s_c___sim.html#ab726ddf220bbef11061b802df50c0188", null ],

diff --git a/docs/functions_func_g.html b/docs/functions_func_g.html
@@ -324,7 +324,7 @@ <h3><a id="index_g"></a>- g -</h3><ul>
 : <a class="el" href="class_excimontec_1_1_o_s_c___sim.html#a7bba196811f642f683570240ceb77e62">Excimontec::OSC_Sim</a>
 </li>
 <li>getSteadyDOS()
-: <a class="el" href="class_excimontec_1_1_o_s_c___sim.html#a6440b290be26622c04ea89cf7293575b">Excimontec::OSC_Sim</a>
+: <a class="el" href="class_excimontec_1_1_o_s_c___sim.html#a15e362e4b5b2c065e29f1cdf36235784">Excimontec::OSC_Sim</a>
 </li>
 <li>getSteadyDOS_Coulomb()
 : <a class="el" href="class_excimontec_1_1_o_s_c___sim.html#ad819d76bf42693993bb392c48327ebf6">Excimontec::OSC_Sim</a>

diff --git a/docs/functions_g.html b/docs/functions_g.html
@@ -327,7 +327,7 @@ <h3><a id="index_g"></a>- g -</h3><ul>
 : <a class="el" href="class_excimontec_1_1_o_s_c___sim.html#a7bba196811f642f683570240ceb77e62">Excimontec::OSC_Sim</a>
 </li>
 <li>getSteadyDOS()
-: <a class="el" href="class_excimontec_1_1_o_s_c___sim.html#a6440b290be26622c04ea89cf7293575b">Excimontec::OSC_Sim</a>
+: <a class="el" href="class_excimontec_1_1_o_s_c___sim.html#a15e362e4b5b2c065e29f1cdf36235784">Excimontec::OSC_Sim</a>
 </li>
 <li>getSteadyDOS_Coulomb()
 : <a class="el" href="class_excimontec_1_1_o_s_c___sim.html#ad819d76bf42693993bb392c48327ebf6">Excimontec::OSC_Sim</a>

diff --git a/docs/navtreeindex0.js b/docs/navtreeindex0.js
@@ -164,6 +164,7 @@ var NAVTREEINDEX0 =
 "class_excimontec_1_1_o_s_c___sim.html#a0b56db2ae5db69979f2d52ae4e6026ba":[1,0,0,1,72],
 "class_excimontec_1_1_o_s_c___sim.html#a10492db50c54e36a401beb119165243b":[1,0,0,1,6],
 "class_excimontec_1_1_o_s_c___sim.html#a1090d447104236755708cb812569cb3c":[1,0,0,1,31],
+"class_excimontec_1_1_o_s_c___sim.html#a15e362e4b5b2c065e29f1cdf36235784":[1,0,0,1,58],
 "class_excimontec_1_1_o_s_c___sim.html#a18785a34bd36bb28de60f4668a06d21f":[1,0,0,1,23],
 "class_excimontec_1_1_o_s_c___sim.html#a1b49d7df062b9184c84685fadcf944bb":[1,0,0,1,10],
 "class_excimontec_1_1_o_s_c___sim.html#a1c18548ebea94e905839f4d824d77935":[1,0,0,1,46],
@@ -191,7 +192,6 @@ var NAVTREEINDEX0 =
 "class_excimontec_1_1_o_s_c___sim.html#a53a0123e8db1299519874ca6bed8bd2b":[1,0,0,1,40],
 "class_excimontec_1_1_o_s_c___sim.html#a576de9d88dafa1a085fa9d6ae269ca87":[1,0,0,1,51],
 "class_excimontec_1_1_o_s_c___sim.html#a63bb36ed57ab8f0029ada80a20be6ed7":[1,0,0,1,70],
-"class_excimontec_1_1_o_s_c___sim.html#a6440b290be26622c04ea89cf7293575b":[1,0,0,1,58],
 "class_excimontec_1_1_o_s_c___sim.html#a6e51664e1d78a6f2afffbc9663f437a7":[1,0,0,1,22],
 "class_excimontec_1_1_o_s_c___sim.html#a72f724474d75a9864014bd0896b8fdf2":[1,0,0,1,3],
 "class_excimontec_1_1_o_s_c___sim.html#a7420ce8e9f8634ab98483fa9ea6b7972":[1,0,0,1,62],

diff --git a/docs/navtreeindex2.js b/docs/navtreeindex2.js
@@ -10,8 +10,8 @@ var NAVTREEINDEX2 =
 "dir_6f62b451fe0710dc1584a1c9c5fb6ce0.html":[2,0,0],
 "dir_d1ce2567778c05124d22b8b3609ec9ea.html":[2,0,0,0],
 "files.html":[2,0],
-"functions.html":[1,3,0,0],
 "functions.html":[1,3,0],
+"functions.html":[1,3,0,0],
 "functions_c.html":[1,3,0,1],
 "functions_d.html":[1,3,0,2],
 "functions_e.html":[1,3,0,3],

diff --git a/docs/search/all_5.js b/docs/search/all_5.js
diff --git a/docs/search/functions_5.js b/docs/search/functions_5.js
diff --git a/parameters_default.txt b/parameters_default.txt
@@ -1,4 +1,4 @@
-## OPV Parameters for Excimontec v1.0.0-rc.1
+## OPV Parameters for Excimontec v1.0.0-rc.3
 --------------------------------------------------------------
 ## Kinetic Monte Carlo Algorithm Parameters
 false //Enable_FRM

diff --git a/src/OSC_Sim.cpp b/src/OSC_Sim.cpp
@@ -2429,7 +2429,7 @@ namespace Excimontec {
 		auto hist = steady_DOOS;
 		// Normalize histogram counts to produce density
 		for (auto& item : hist) {
-			item.second /= ((params.N_tests / Steady_hops_per_DOOS_sample) + 1)*lattice.getVolume()*DOS_bin_size;
+			item.second /= Steady_DOOS_sampling_counter * lattice.getVolume()*DOS_bin_size;
 		}
 		return hist;
 	}
@@ -2438,16 +2438,38 @@ namespace Excimontec {
 		auto hist = steady_DOOS_Coulomb;
 		// Normalize histogram counts to produce density
 		for (auto& item : hist) {
-			item.second /= ((params.N_tests / Steady_hops_per_DOOS_sample) + 1)*lattice.getVolume()*DOS_bin_size;
+			item.second /= Steady_DOOS_sampling_counter * lattice.getVolume()*DOS_bin_size;
 		}
 		return hist;
 	}
 
-	vector<std::pair<double, double>> OSC_Sim::getSteadyDOS() const {
+	vector<std::pair<double, double>> OSC_Sim::getSteadyDOS() {
+		steady_DOS.clear();
+		for (long int i = 0; i < lattice.getNumSites(); i++) {
+			double energy;
+			auto site_coords = lattice.getSiteCoords(i);
+			if (lattice.isOccupied(site_coords)) {
+				if (getSiteType(site_coords) == 1) {
+					energy = params.Homo_donor + getSiteEnergy(site_coords);
+				}
+				else {
+					energy = params.Homo_acceptor + getSiteEnergy(site_coords);
+				}
+			}
+			else {
+				if (getSiteType(site_coords) == 1) {
+					energy = params.Homo_donor + getSiteEnergy(site_coords);
+				}
+				else {
+					energy = params.Homo_acceptor + getSiteEnergy(site_coords);
+				}
+			}
+			updateSteadyDOS(steady_DOS, energy);
+		}
 		auto hist = steady_DOS;
 		// Normalize histogram counts to produce density
 		for (auto& item : hist) {
-			item.second /= ((params.N_tests / Steady_hops_per_DOS_sample) + 1)*lattice.getVolume()*DOS_bin_size;
+			item.second /= lattice.getVolume()*DOS_bin_size;
 		}
 		return hist;
 	}
@@ -2456,7 +2478,7 @@ namespace Excimontec {
 		auto hist = steady_DOS_Coulomb;
 		// Normalize histogram counts to produce density
 		for (auto& item : hist) {
-			item.second /= (((params.N_tests / Steady_hops_per_DOS_sample) + 1))*lattice.getVolume()*DOS_bin_size;
+			item.second /= Steady_DOS_sampling_counter * lattice.getVolume()*DOS_bin_size;
 		}
 		return hist;
 	}
@@ -2901,39 +2923,33 @@ namespace Excimontec {
 					Steady_equilibration_energy_sum += energy;
 					Steady_equilibration_energy_sum_Coulomb += energy_C;
 				}
+				Steady_DOOS_sampling_counter++;
 			}
 			// Sample the density of states
 			if ((N_events_executed - params.N_equilibration_events) % Steady_hops_per_DOS_sample == 0) {
 				for (long int i = 0; i < lattice.getNumSites(); i++) {
-					double energy;
 					double energy_C;
 					auto site_coords = lattice.getSiteCoords(i);
 					if (lattice.isOccupied(site_coords)) {
 						auto it = getPolaronIt((*lattice.getSiteIt(site_coords))->getObjectPtr());
 						if (getSiteType(site_coords) == 1) {
-							energy = params.Homo_donor + getSiteEnergy(site_coords);
 							energy_C = params.Homo_donor + getSiteEnergy(site_coords) + calculateCoulomb(it, site_coords);
 						}
 						else {
-							energy = params.Homo_acceptor + getSiteEnergy(site_coords);
 							energy_C = params.Homo_acceptor + getSiteEnergy(site_coords) + calculateCoulomb(it, site_coords);
 						}
-						updateSteadyDOS(steady_DOS, energy);
-						updateSteadyDOS(steady_DOS_Coulomb, energy_C);
 					}
 					else {
 						if (getSiteType(site_coords) == 1) {
-							energy = params.Homo_donor + getSiteEnergy(site_coords);
 							energy_C = params.Homo_donor + getSiteEnergy(site_coords) + calculateCoulomb(true, site_coords);
 						}
 						else {
-							energy = params.Homo_acceptor + getSiteEnergy(site_coords);
 							energy_C = params.Homo_acceptor + getSiteEnergy(site_coords) + calculateCoulomb(true, site_coords);
 						}
-						updateSteadyDOS(steady_DOS, energy);
-						updateSteadyDOS(steady_DOS_Coulomb, energy_C);
 					}
+					updateSteadyDOS(steady_DOS_Coulomb, energy_C);
 				}
+				Steady_DOS_sampling_counter++;
 			}
 		}
 	}
@@ -2976,7 +2992,7 @@ namespace Excimontec {
 		}
 		// Calculate start of data range
 		double min_val = smallest_bin_int * DOS_bin_size - 0.5*DOS_bin_size;
-		// Add to state_energy to histogram
+		// Add the state_energy to histogram
 		int index = (int)floor((state_energy - min_val) / DOS_bin_size);
 		density_of_states[index].second += 1.0;
 		return;

diff --git a/src/OSC_Sim.h b/src/OSC_Sim.h
@@ -299,7 +299,7 @@ namespace Excimontec {
 
 		//! \brief Gets the density of states calculated during the steady state charge transport test.
 		//! \return A pair vector where the first value is the state energy and the second in the density of states.
-		std::vector<std::pair<double, double>> getSteadyDOS() const;
+		std::vector<std::pair<double, double>> getSteadyDOS();
 
 		//! \brief Gets the density of states calculated during the steady state charge transport test.
 		// This function calculates the state energies including the Coulomb potential due to interactions between the polarons.
@@ -557,6 +557,8 @@ namespace Excimontec {
 		std::vector<int> transient_hole_counts;
 		int Steady_hops_per_DOS_sample = 1000000;
 		int Steady_hops_per_DOOS_sample = 1000;
+		int Steady_DOS_sampling_counter = 0;
+		int Steady_DOOS_sampling_counter = 0;
 		double DOS_bin_size = 1e-2;
 		double Steady_equilibration_time = 0.0;
 		double Steady_equilibration_energy_sum = 0.0;