Author Archives: Metals, Magnets, and Miscellaneous Materials

Computationally Visualizing Crystals

Christina C. Lee, github: albi3ro

Prerequisites: none

In condensed matter, we find ourselves in the interesting middle ground of dealing with large numbers, e.g. , of extremely small particles such as atoms, or electrons.

Luckily, the particles don’t each do their own thing, but often come in nice, structured, repeated units. Lattices. As our first step into the field, we will look at the most basic type, a Bravais lattice.

In a Bravais lattice, every site looks like every other site. Mathematically, we use three vectors, $vec{a},vec{b},vec{c}$ to express how we move from one site to a neighbor.

begin{equation}
mathbf{R}_{lmn}=l vec{a} + m vec{b} + n vec{c} ;;;; text{for } l,m,n in mathbb{N}
end{equation}

For consistency, we have to put a constraint on these vectors; we cannot combine two of the vectors and obtain the third. If we could, then we couldn’t have sites in an entire three dimensional space.

Stay tuned for a later post where we explore more elaborate lattices.

# importing our packages
Pkg.add("PyPlot");
Pkg.update();
using PyPlot;

Define The Relevant Variables

Choose the lattice you want to look at, and use that string for the lattice variable.
Current options:

  • Simple Cubic = “sc”
  • Plane triangular lattice = “pt”
  • Body-Centered Cubic = “bcc”
  • Face-Centered Cubic = “fcc”

Note: Square is Simple Cubic for Nz=1

14 distinct lattice types are possible, but these common four give the important ideas.

Also, input the size of lattice you want to look at.

lattice="sc";

Nx=3;
Ny=3;
Nz=3;
# A cell to just evaluate
# This one sets the unit vectors (a,b,c) for the different unit cells
# Can you guess what a lattice will look like by looking at the vectors?
if(lattice=="sc")
    d=3;
    a=[1,0,0];
    b=[0,1,0];
    c=[0,0,1];
elseif(lattice=="pt")
    d=2;
    a=[1,0,0];
    b=[.5,sqrt(3)/2,0];
    c=[0,0,1];
elseif(lattice=="bcc")
    d=3;
    a=[.5,.5,.5];
    b=[.5,.5,-.5];
    c=[.5,-.5,.5];
elseif(lattice=="fcc")
    d=3;
    a=[.5,.5,0];
    b=[.5,0,.5];
    c=[0,.5,.5];
else
    println("Please have a correct lattice")
end
# Another cell to just evaluate
# Here we set up some numbers and matrices for our computation
N=Nx*Ny*Nz;    #The total number of sites
aM=transpose(a);
bM=transpose(repeat(b,outer=[1,Nx])); #these allow us to copy an entire row or layer at once
cM=transpose(repeat(c,outer=[1,Nx*Ny]));

X=Array{Float64}(N,3);  #where we store the positions
# Another cell to just evaluate
# Here we are actually calculating the positions for every site
for i in 1:Nx    #for the first row
    X[i,:]=(i-1)*a;
end

for j in 2:Ny    #copying the first row into the first layer
    X[Nx*(j-1)+(1:Nx),:]=X[1:Nx,:]+(j-1)*bM;
end

for j in 2:Nz    #copying the first layer into the entire cube
    X[Ny*Nx*(j-1)+(1:Nx*Ny),:]=X[1:Nx*Ny,:]+(j-1)*cM;
end

Programming Tip:

In Julia, ranges, like 1:Nx, are a special variable type that can be manipulated. We can add numbers to them:
3+(1:3)=4:6,
or add a minus sign to force it to iterate in the opposite direction, though with different start/stop:
-(1:3)=-1:-1:-3

Danger! Make sure to use the parentheses around the range if you are performing these operations.

pygui(false);  #if true, launches new window with interactive capabilities

drawcube=true;  #gives lines for a cube, helps interpret the dots
ls=2;  # how many cubes to draw
if(drawcube==true)
    v=collect(0:ls);
    zed=zeros(v);
    for i in 0:ls
        for j in 0:ls
            plot3D(zed+i,v,zed+j)
            plot3D(zed+i,zed+j,v)

            plot3D(v,zed+i,zed+j)
            plot3D(zed+j,zed+i,v)

            plot3D(v,zed+j,zed+i)
            plot3D(zed+j,v,zed+i)
        end
    end
end

scatter3D(X[:,1],X[:,2],X[:,3],s=200*ones(X[:,1]),alpha=1)

sc

Simple Cubic: The easiest lattice out there short of the 1D chain.

pt

Point Triangular: A 2D lattice. Plotted using scatter instead of scatter3D.

bcc

Body Centered Cubic: Notice how some sites fall on the cubic lattice, but others fall in between. Generated with pygui(true) and then manipulating in 3D.

fcc

Face Centered Cubic: Here the sites either fall on the the cubic corners of in the center of the sides. Generated with pygui(true), ls=1, and then manipulating in 3D.

Go Back and Fiddle!

As you might have noticed, this isn’t just a blog where you read through the posts. Interact with it. Change some lines, and see what happens. I choose body centered cubic to display first, but what do the other lattices look like?

Chose pygui(true) to pop open a window and manipulate the plot in 3D.

Look at different lattice sizes.

Can you hand draw them on paper?

handdrawn fcc

A face centered cubic I decided to draw myself.

Let me know what you think, and enjoy the sequel as well!

Multi-site unit cells

Computationally Visualizing Crystals Pt. 2

Christina C. Lee, github: albi3ro

Prerequisites: Computationally Visualizing Crystals Pt. 1

Time to one-up the Bravais lattice from Part 1. Many beautiful lattices don’t adhere to the “every site the same” policy. They still repeat, but just take a little bit longer to get around to doing so.

Take the Kagome Lattice below,

kagome

A Kagome Lattice.

basket

A basket woven in the Japanese kagome style. Wikimedia commons

If we look at the stars at the center of triangles, we can recognize a point triangular Bravais lattice. Now each of those stars stands for a grouping of three sites in a unit cell. According to Chem Wiki, a unit cell is:

A unit cell is the most basic and least volume consuming repeating structure of any solid. It is used to visually simplify the crystalline patterns solids arrange themselves in.

I chose these triangles to be be the unit cells above and in my computational representation below, but can you think of any other ways to represent the unit cell?

Turns out, there isn’t a unique way. We can go further and define the Wigner-Seitz unit cell, which uses the Bravais translations to pick out just ONE of the various possible definitions.

In my line of work though, we often use either the easiest to write down, or the one that has the symmetries we want.

Introducing Some Lattice Options

You saw Kagome above.

The options I’ve put in now are:

  • honeycomb
  • kagome
  • shuriken (a.k.a. square-Kagome)
  • diamond
  • pyrochlore

The ones implemented here, except for diamond, are frustrated lattices that I work with in my research. Honeycomb is well known in condensed matter physics for being the structure of graphene. This is an extremely important material right now, though I work with it in terms of the Kitaev spin model (to be discussed at a later date). Kagome and pyrochlore are also popular models within my community. The shuriken lattice is more uncommon, but gaining ground in the frustration community.

Shurikens

Japanese Shurikens- a type of ninja fighting star. By kaex0r (http://www.flickr.com/photos/kaex0r414/191765028/) [CC BY 2.0 (http://creativecommons.org/licenses/by/2.0)], via Wikimedia Commons

# importing our packages
Pkg.add("PyPlot");
Pkg.update();
using PyPlot;
lattice="shuriken";

Nx=3;
Ny=3;
Nz=1;
# A cell to just evaluate
# This one sets the unit vectors (a,b,c) for the different unit cells
# Can you guess what a lattice will look like by looking at the vectors?
if(lattice=="honeycomb") #also the graphite lattice
    d=2;
    Ncell=2;
    unit=[[0 0 0]
        [sqrt(3)/2 1/2 0]];
    a=[sqrt(3),0,0];
    b=[sqrt(3)/2,3/2,0];
    c=[0,0,1];
elseif(lattice=="kagome")
    d=2;
    Ncell=3;
    unit=[[0 0 0]
          [1 0 0]
        [.5 sqrt(3)/2 0]];
    a=[2,0,0];       #Look familiar? Checkout pt from Pt. 1
    b=[1, sqrt(3), 0];
    c=[0,0,1];
elseif(lattice=="shuriken")
    d=2;
    Ncell=6;
    unit=[[0 0 0]
          [.5 0 0]
          [0 .5 0]
          [.5 .5 0]
        [.5+.25*sqrt(3) .25 0]
        [.25 .5+.25*sqrt(3) 0]];
    a=[.5+.5*sqrt(3),0,0];
    b=[0,.5+.5*sqrt(3),0];
    c=[0,0,1];
elseif(lattice=="diamond")
    d=3;
    Ncell=2;
    unit=[[0 0 0]
          [.25 .25 .25]];
    a=[.5,.5,0];    #Look familiar? Checkout fcc from Pt.1
    b=[.5,0,.5];
    c=[0,.5,.5];
elseif(lattice=="pyrochlore")
    d=3;
    Ncell=4;
    unit=[[0 0 0]
        [.25 .25 0]
        [.25 0 .25]
        [0 .25 .25]];
    a=[.5,.5,0];
    b=[.5,0,.5];
    c=[0,.5,.5];

else
    println("Please have a correct lattice")
end
"Cell finished"

Connections to Bravais Lattices

If you look at some of the comments above, and checkout the basis vectors from Crystal Shapes, like pt,
begin{equation}
a=[1,0,0];;;;;;;;; b=[.5,frac{sqrt{3}}{2},0],
end{equation}
you’ll notice they’re the same except for a scaling factor. This has to be true, since only 14 different patterns tile 3D space uniquely.

# Another cell to just evaluate
# Here we set up some numbers and matrices for our computation
N=Nx*Ny*Nz*Ncell;    #The total number of sites
aM=transpose(repeat(a,outer=[1,Ncell]));
bM=transpose(repeat(b,outer=[1,Ncell*Nx])); #these allow us to copy an entire row or layer at once
cM=transpose(repeat(c,outer=[1,Ncell*Nx*Ny]));

X=Array{Float64}(N,3);  #where we store the positions
"Cell finished"
# Another cell to just evaluate
# Here we are actually calculating the positions for every site
for i in 1:Nx    #for the first row
    X[Ncell*i-Ncell+1:Ncell*i,:]=unit+(i-1)*aM;
end

for j in 2:Ny    #copying the first row into the first layer
    X[Ncell*Nx*(j-1)+(1:Ncell*Nx),:]=X[1:Ncell*Nx,:]+(j-1)*bM;
end

for j in 2:Nz    #copying the first layer into the entire cube
    X[Ncell*Ny*Nx*(j-1)+(1:Ncell*Nx*Ny),:]=X[1:Ncell*Nx*Ny,:]+(j-1)*cM;
end
"Cell finished"
# 2D plotter
pygui(false)
w, h = plt[:figaspect](1)
figure(figsize=(w,h))
scatter(X[:,1],X[:,2])

Shuriken

3×3 Shuriken or Square-Kagome Lattice.

# 3D plotter
pygui(false)
w, h = plt[:figaspect](1)
figure(figsize=(w,h))
areas=100*ones(length(X[:,1]))
scatter3D(X[:,1],X[:,2],X[:,3])

Perdy Pictures

From these plots, some 3D structures like the pyrochlore are hard to visualize. So here’s a nice graphic I made that might help a little bit more.

Pyrochlore

Hopefully this pyrochlore is a little easier to visualize than the pyplot version. Took me long enough to make in inkscape.

honeycomb

Tikz produced Honeycomb. Coloring indicative of the lattice description of the Kitaev model.

The honeycomb, like several other lattices you see around here, is bipartite.
You can see in my image that black sites are only next to white sites, and vice versa.
This property can make the system much easier to work with. What lattices are bipartite, and which ones aren’t?

If you keep reading, these lattices will keep cropping up again and again. I’ll probably throw in some new ones as well.

Anyway, we will move onto some Quantum Mechanics to look at atomic orbitals soon!

Quantum Harmonic Oscillator

Prerequiresites: Quantum Mechanics course

Slinkies. They started out as toys. I still have one to play with on my desk.

Rubber bands What was once something useful, is now a wonderful projectile weapon.

Swings I still love them, but people seem to not make them in adult sizes for some reason.

A person’s perception of these objects starts to change as they enter their first physics class. Even in that beginning class of classical mechanics, the problems are filled with harmonic oscillators, like slinkies, rubber bands, or swings, which exert a force proportional to their displacement

and therefore a quadratic potential

This is all extremely fun and useful in the classical regime, but we add Quantum Mechanics to the mix, and LOW AND BEHOLD! we have one of the few exactly solvable models in Quantum Mechanics. Moreso, this solution demonstrates some extremely important properties of quantum mechanical systems.

The Hamiltonian

The Solution

Today, I just intend to present the form of the solution, calculate this equation numerically, and visualize the results. If you wish to know how the equation is derived, you can find a standard quantum mechanics textbook, or stay tuned till I manage to write it up.

Physicists’ Hermite Polynomials

Note: These are not the same as the “probabilists’ Hermite Polynomial”. The two functions differ by scaling factors.

Physicists’ Hermite polynomials are defined as eigenfunctions for the differential equation

I leave it as an exercise to the reader to:

  • demonstrate orthogonality with respect to the measure $e^{-x^2}$, i.e.

  • demonstrate completeness. This means we can describe every function by a linear combination of Hermite polynomials, provided it is suitably well behaved.

Though a formula exists or calculating a function at $n$ directly, the most efficient method at low $n$ for calculating polynomials relies on recurrence relationships. These recurrence relationships will also be quite handy if you ever need to show orthogonality, or expectation values.

Recurrence Relations

#Pkg.update();
#Pkg.add("PyPlot");
#Pkg.update()
#Pkg.add("Roots")
using Roots;
using PyPlot;

Programming Tip!

Since Hermite polynomials are generated recursively, I wanted to generate and save all the functions up to a designated value at once. In order to do so, I created an array, whose values are anonymous functions.

function GenerateHermite(n)
    Hermite=Function[]

    push!(Hermite,x->1);
    push!(Hermite,x->2*x);

    for ni in 3:n
        push!(Hermite,x->2.*x.*Hermite[ni-1](x).-2.*n.*Hermite[ni-2](x))
    end
    return Hermite
end

Let’s generate some Hermite polynomials and look at them.
Make sure you don’t call a Hermite you haven’t generated yet!

Hermite=GenerateHermite(5)

Programming Tip!

Since the Hermite polynomials, and the wavefunctions after them, are composed of anonymous functions, we need to use map(f,x) in order to map the function f onto the array x. Otherwise our polynomials only work on single values.

x=collect(-2:.01:2);
for j in 1:5
    plot(x,map(Hermite[j],x),label="H_$j (x)")
end
legend()
ylim(-50,50)

Hermite Polynomials

First few Hermite Polynomials

# Lets make our life easy and set all units to 1
m=1
ω=1
ħ=1

#Finally, we define Ψ
Ψ(n,x)=1/sqrt(factorial(n)*2^n)*(m*ω/(ħ*π))^(1/4)*exp(-m*ω*x^2/(2*ħ))*Hermite[n](sqrt(m*ω/ħ)*x)

Finding Zeros

The eigenvalue maps to the number of zeros in the wavefunction. Below, I use Julia’s roots package to indentify roots on the interval from -3 to 3.

zeds=Array{Array{Float64}}(1)
zeds[1]=[] #ground state has no zeros
for j in 2:4
    push!(zeds,fzeros(y->Ψ(j,y),-3,3))
end
# AHHHHH! So Much code!
# Don't worry; it's all just plotting
x=collect(-3:.01:3)  #Set some good axes

for j in 1:4    #how many do you want to view?
    plot(x,map(y->Ψ(j,y),x)+j-1,label="| $j >")
    plot(x,(j-1)*ones(x),color="black")
    scatter(zeds[j],(j-1)*ones(zeds[j]),marker="o",s=40)
end
plot(x,.5*m*ω^2*x.^2,linestyle="--",label="Potential")

scatter([],[],marker="o",s=40,label="Zeros")
xlabel("x")
ylabel("Ψ+n")
title("Eigenstates of a Harmonic Osscilator")
legend()
xlim(-3,3);
ylim(-.5,4.5);

Example Result

Eigenstates

Eigenstates of the Quantum Harmonic Osscillator

More to come

This barely scratched the surface into the richness that can be seen in the quantum harmonic oscillator. Here, we have developed a way for calculating the functions, and visualized the results. Stay tuned to hear about ground state energy, ladder operators, and atomic trapping.