Finding the Perfect Diamond:
Why it’s Impossible

Gloria Chyr 
3.016 Final Project

What a perfect diamond looks like

A quick Google image search shows this as the “perfect” diamond:

Looks pretty great! But is it really that perfect?

Structure of Diamonds

Diamonds are made of carbon atoms arranged in a crystal lattice. As you may already know, crystals have long-range order which means that the carbon atoms that make up diamond can be found in specific positions. The structure of carbon atoms that makes diamond special is:

Image source: Google -> Quora

A perfect diamond, like any perfect crystal, would have no defects. In short, the structure of diamond would look exactly as it does above. But in the real world, this is impossible to find. You’ll likely find lots of defects such as vacancies pictured below.

But if crystals are these perfectly ordered objects, why do vacancies form? Why would any crystal at equilibrium have defects if it takes more energy to cause them?

Why defects form

The formation of defects in a crystal costs energy which is not favorable. But the formation of defects also increases the amount of entropy in the system which is favorable.

The entropy gain can sometimes outweigh the energy cost of forming a defect. So a perfect crystal can actually decrease its Gibb’s free energy by forming defects. This can be seen in the following equation:

ΔG = ΔH - TΔS

So formation of defects becomes a battle of enthalpy vs entropy. If the enthalpy cost of forming a defect is greater than the entropy increase of the system, then the defect does not form. If the entropy increase of the system is greater than the enthalpy of defect formation, then the defect forms.

Visualizing the formation of defects

Here, we can see two things:
1. The equilibrium number of vacancies depends on temperature.
2. The equilibrium number of vacancies is never zero.

An explanation of statement 2 can be found in the code below.

Un-hide this cell to see the code!

Creating the lattice:
A simple square lattice is used to represent the diamond lattice for ease of viewing. The concept is the same. An accurate visualization of diamond can be made simply by replacing the code that draws the lattice.

range =Table[a,{a,1,20}] ;(*this is going to be my point range*)

points =Flatten[Table[{j,range[[i]]},{j,1,20},{i,1,20}],1];
(*creates a list of points in a 20x20 lattice*)

vacancies={};

amountOfVacancies[n_]:=Module[{vacancy,disk},For[i=1,i<(n+1),i++,
vacancy = points[[RandomInteger[{1,400}]]];
disk = Graphics[{White,Disk[vacancy, 0.5]}]; 
If[MemberQ[vacancies,disk],i-=1,(*creates graphic for the vacancy*)
AppendTo[vacancies,disk];]]]

lattice =Graphics[{EdgeForm[Black],Pink,Table[Disk[{i,j},1/2],{i,20},{j,20}]},FrameLabel->"Visual example of vacancies appearing in the lattice"];(*draws the lattice*)

drawVacancies[n1_]:= Module[{n=n1}, 
vacancies ={};amountOfVacancies[n];
Show[lattice,vacancies]] (*draws the vacancies into the lattice*)

Calculating Gibb’s energy of the system:

kB = 8.6*(10^(-5));(*Boltzmann constant in eV/K*)
hV = 1.5;(*enthalpy of vacancy formation in eV*)
sVib = 0.0001; (*sVib in eV/K*)(*I'm using fake numbers (adjusted from Cu) for the purpose of teaching. The real system is much more complicated to model and finding the exact values for diamond would be nearly impossible*)

gibbsOfVacancy[n_,latticeSites_,temp_]:= (*calculates the gibbs energy*)n*hV  - n*temp*sVib-kB*temp*((latticeSites+n)*Log[latticeSites+n] -latticeSites*Log[latticeSites] - n*Log[n])

The equation that determines Gibb’s free energy depends on the number of vacancies in the lattice. Within the equation, we see that there is a Log[n] term where n is the number of vacancies. Because of this Log term, we see that there can not be a zero amount of vacancies.

latticeSites =400; (*number of sites in our lattice*)

(*plots the gibbs energy at various temperatures*)
plotManipulate[temp_]:=Plot[gibbsOfVacancy[n,latticeSites,temp],{n,1,500},Frame-> True,FrameLabel-> {"Number of Vacancies","Gibb's Free Energy (eV)"}, ImageSize-> Large,BaseStyle-> {FontSize-> 24},PlotLabel-> Style["Gibb's Energy Plot", FontSize-> 24]]

(*returns the number of vacancies present at minimum gibbs energy*)
findMinimum[temp_]:=Module[{data,position},
data = Table[{n,gibbsOfVacancy[n,400,temp]},{n,1,400,1}];
position=Position[data, Min@data];
Flatten[position][[1]]]

Manipulate[
GraphicsRow[{
plotManipulate[Temperature],
drawVacancies[findMinimum[Temperature]]
}]
,
{Temperature,1,5000,1}] (*Show gibbs energy plot and lattice next to each other. 
Temperature can be manipulated from 1 to 5000 K*)