[ANN] SimpleQuantum: A package for calculating properties of solids + a tutorial blog on how to do these calculations in Julia

SimpleQuantum is now available for easy-to-setup calculations of (for now just electronic) properties of solids. It provides ways to define and visualize crystal structures and a simple way to do common tasks, like calculate band structures. For now, it can do calculations based on tight binding and nearly free electron models. Associated with it is the Computational Physics for the Masses blog series that documents, line-by-line, how the algorithms are implemented and explains the physics to people with minimal background.

A few quick snippets showing the basic use:

Band diagram of toy model graphene

using SimpleQuantum

# Define the graphene crystal.
graphene = Crystal(
    Lattice(2.468Å, 2.468Å, 120),
    UnitCell(:C, [2/3, 1/3], [1/3, 2/3])
)

# Create an empty hopping list based on graphene.
grhops = Hoppings(graphene)

# Iterate through unique nearest neighbor pairs and add them to the hopping list.
for hop ∈ SimpleQuantum.unique_neighbors(graphene)
    addhop!(grhops, -1.0Ha, hop.i, hop.j, hop.δ)
end

# Assemble the tight binding problem.
grprob = TightBindingProblem(grhops, [
    :K => [1/3,1/3],
    :Γ => [0,0],
    :M => [1/2,0],
    :K => [1/3,1/3]
], 0.005)

# Solve the problem.
sol = solve(grprob)

# Plot the band diagram.
plotSolution(sol)

Which outputs:

Visualization of the diamond structure:

using SimpleQuantum
using GLMakie

diamond = Crystal(
    Lattice(2.527Å, 2.527Å, 2.527Å, 60, 60, 60),
    UnitCell(:C, [0.0, 0.0, 0.0], [1/4, 1/4, 1/4])
)

fig3d = Figure()
axds = Axis3(fig3d, xlabel="x/a₀", ylabel="y/a₀", zlabel="z/a₀", xgridvisible=false, ygridvisible=false, zgridvisible=false, aspect = :data)
axds_single = Axis3(fig3d, xlabel="x/a₀", ylabel="y/a₀", zlabel="z/a₀", xgridvisible=false, ygridvisible=false, zgridvisible=false, aspect = :data)

plotcrystal!(axds, diamond; ncells=2, showcell=false, showbonds=true)
plotcrystal!(axds_single, diamond; ncells=1, showcell=true, showbonds=true)
fig3d[1,1] = axds
fig3d[1,2] = axds_single

current_figure()

Outputs:

15 Likes

Oooh, very nice package and especially blog. I wish I had this when starting out! It also makes for a nice background before tackling more quantitative methods such as DFT, I might give this as a resource to students. I’m looking forward to the next posts :slight_smile: May I suggest something more descriptive as a package name? This is a bit generic… Also note that the geometry stuff you can also offload to Brillouin.jl (but I understand the desire to be self contained!)

3 Likes

The package right now is not meant to be pushed to registry and is mostly there to accompany the tutorials. Once it gets all the functionality I want out of it (or close enough), I’ll probably make a fork that isn’t burdened with the pedagogical baggage of the current iteration and do a lot of optimization. At that point, I’ll pick a less generic sounding name as that will be more “production-ready”.

I haven’t even thought of looking if there is something for the geometry but Brillouin.jl does look interesting. It is a tad bit too specialized to be a drop-in replacement, but that’s something I might be able to work around.

Is it possible to calculate transition dipole moments with this package?

It is very beautiful!!!

Thanks for this, I will read and try your codes.

This is really great!

No. That would require it to have information about the shape of the orbitals, which tight binding doesn’t need or have at all, and the nearly free electron is just select plane waves.

If you know your orbitals, you can represent your dipole operator in that basis and use the tight binding solver to find the projection of the bands onto that basis (see the blog post, I do that there for the full graphene structure). The moment will then be just <band_a| dipole_op |band_b>.