Band Structure Simulation using Fortran

Job ID: 38538681

Budget: €8 – €30 EUR

I need a code to perform electronic band structure calculation given a periodic one dimensional potential (Kronig-Penney).
Basically, for each value of real parameter k in a lattice, square matrix is constructed and then it's eigenvalues E_n(k) are plotted as a function of k.

Workflow of the program is roughly as follows:
1) Three real positive parameters a, b, V0 are given which determine the shape of the potential. The Kronig-Penney potential is a function of one real variable, basically a step function (see the attached picture). You can shift it if necessary. They can be loaded from stdin, file or even hardcoded.
2) The potential is expanded into a complex Fourier series. Only the coefficient V_G where G is between -G_cutoff and G_cutoff are kept and rest is discarded (the coefficients will be complex). Please let the array of coefficients actually be indexed like that for clarity. In this case this is easily done analytically (ask ChatGPT if you don't know the math), but I prefer the more general numerical approach (added challenge, not strictly necessary).
3) Lattice of values of the real variable k is constructed. It shall have k_count elements evenly spanning the interval <0 , pi/(a+b) ) .
4) For each value of k, square matrix H with 2*G_cutoff+1 rows is constructed. Rows and columns are again to be indexed from -G_cutoff to G_cutoff (i,j below runs over this range)
H_i,i = 1/2 * (k + 2*pi/(a+b) * i)^2 + V_0 on the diagonal
H_i,j = V_(i-j) off the diagonal
5) Eigenvalues of this matrix are found and saved (say as rows in a matrix BandsMatrix with k_count rows)
6) Once this has been done for all values of k in the lattice, the simulation is complete. Save the BandsMatrix and k-lattice in file(s) (can be formatted or binary).
7) Plot each column of BandMatix against the k-lattice, preferably in Octave/Matlab (if it is not your strong suit it's not strictly necessary).

90% percent of the challenge lies with steps 2 and 4, rest should be trivial.

Default parameters for testing (unitless) :
a = 3
b = 1
V0 = 1
G_cutoff = 4
k_count = 100

This problem lends itself well to parallelization, since the computation is independent for all values of k, so sneaking in some minimal MPI would be great. External libraries are OK, as long as they are reasonably low level and freely available. Result should look something like the second attached picture.Thanks!
Related categories: Fortran Physics Visualization