With the initial potential computed by assempde, HADOKEN reinserts a fraction of the total potential (which includes an additional exchange-correlation term calculated by the V_xc.m and V_total_with_xc.m m-files) into a new Schrödinger equation to initialize the self-consistent procedure. To complete one cycle of the iteration scheme shown in Fig. , HADOKEN uses a 0.01 fraction of the potential for both the charge-neutral and Fermi-pinned hexagonal geometries, and a 0.05 fraction of the potential for triangular geometries. This cyclic process of simultaneously solving the Schrödinger and Poisson equations is continued until both are self-consistent, which we define as the situation where the average energy difference over all nodes of the electrostatic potential between successive iterations is less than 0.01 eV. To maintain a stable self-consistent feedback loop, the HADOKEN code uses an under-relaxation technique set by the variable damping_factor, which is initially equal to 0.03 for hexagonal geometries and 0.05 for triangular geometries. As such, the inputted potential for the next iteration, V_poisson, is calculated as V_poisson=V_poisson_old+damping_factor*(V_poisson-V_poisson_old), where V_poisson_old is the potential just computed. When the potential is nearly converged-defined to be the average energy difference of the potential being less than 0.01C-the damping_factor variable is slowly increased to a maximum of 0.07 to further accelerate convergence. Once self-consistency is reached, the band-bending diagram, total electron density, and all occupied wavefunctions are output to the screen. Additionally, all variables used by HADOKEN are saved to a binary .mat file for further post-processing by the user.