Phononen Dispersion von Hexagonal Close Packed

Einführung

Hexagonal Close Packed unterscheidet sich von manchen anderen Krystallstrukturen dadurch, dass es 2 Atome in der primitiven Einheitszelle hat. Das macht die Sache also komplizierter. Wer sich sowas das erste mal ansieht sollte mit einem einfacheren Beispiel wie fcc beginnen.

Parameter

(*Federkonstante zum nächsten Nachbarn*)
springconst=1;
(*Atommasse*)
m=1;
(*Gitterkonstanten*)
a=1;
c=2*Sqrt[6]/3*a;
(*Definition der primitiven Basisvektoren*)
a1={a/2,-a*Sqrt[3]/2,0};
a2={a/2,a*Sqrt[3]/2,0};
a3={0,0,c};
(*Sitz der Atome, siehe Abbildung unten*)
atom1={0,0,0};
atom2={1/3,2/3,1/2};
(*Sitz der nächsten Nachbarn von Atom1 und Atom2, die ersten 6 sind jeweils u Terme (gehören also in die gleiche Gruppe wie Atom1), die letzen 6 sind v Terme (also die Gruppe von Atom2)*)
pos1={{0,1,0},{1,1,0},{1,0,0},{0,-1,0},{-1,-1,0},{-1,0,0},{(1/3),-(1/3),(1/2)},{(1/3),(2/3),(1/2)},{-(2/3),-(1/3),(1/2)},{(1/3),(2/3),-(1/2)},{-(2/3),-(1/3),-(1/2)},{(1/3),-(1/3),-(1/2)}};
pos2={{0,1,0},{0,0,0},{1,1,0},{0,0,1},{1,1,1},{0,1,1},{(1/3),-(1/3),(1/2)},{-(2/3),-(1/3),(1/2)},{-(2/3),(2/3),(1/2)},{(1/3),(5/3),(1/2)},{(4/3),(5/3),(1/2)},{(4/3),(2/3),(1/2)}};
(*Symmetriepunkte, entlang dieser Punkte wird die Dispersionsrelation geplottet*)
symPoints={{0,0,0},1/3 *b1+1/3 *b2,1/2*b2,{0,0,0},1/2*b3,1/2 *b2+1/2* b3,1/2*b2};
(*Der Name der Symmetriepunkte für den Plot*)
label={"Γ","K","M","Γ","A","L","M"};

Berechnen der reziproken Basisvektoren:

VE=a1.(a2\[Cross]a3);
b1=2*Pi/VE*(a2\[Cross]a3);
b2=2*Pi/VE*(a3\[Cross]a1);
b3=2*Pi/VE*(a1\[Cross]a2);

Hier eine Abbildung der 1. Brillouin Zone:

http://cst-www.nrl.navy.mil/users/mehl/phonons/hcp/

Graphische Darstellung der räumlichen Verteilung der nächsten Nachbarn

In rot sind die primitiven Basisvektoren eingezeichnet. Ein Eckpunkt der primitiven Einheitszelle ist frei weil dieses Atom kein nächster Nachbar ist von Atom1 oder Atom2.
Das rote Atom ist Atom1 (1. Atom in der primitiven Einheitszelle), das blaue ist Atom2. Die restlichen rötlichen Atome sind "u" Atome, spielen also die gleiche Rolle wie Atom1 in ihrer Einheitszelle, die bläulichen Atome sind "v" Atome und sind mit Atom2 in ihrer Einheitszelle äquivalent.

Graphics3D[{{Magenta,PointSize[.1],Point[Table[pos1[[i]][[1]]*a1+pos1[[i]][[2]]*a2+pos1[[i]][[3]]*a3,{i,1,6}]]},{Magenta,PointSize[.1],Point[Table[pos2[[i]][[1]]*a1+pos2[[i]][[2]]*a2+pos2[[i]][[3]]*a3,{i,1,6}]]},{Cyan,PointSize[.1],Point[Table[pos1[[i]][[1]]*a1+pos1[[i]][[2]]*a2+pos1[[i]][[3]]*a3,{i,7,12}]]},{Cyan,PointSize[.1],Point[Table[pos2[[i]][[1]]*a1+pos2[[i]][[2]]*a2+pos2[[i]][[3]]*a3,{i,7,12}]]},{Red,PointSize[.1],Point[atom1[[1]]*a1+atom1[[2]]*a2+atom1[[3]]*a3]},{Blue,PointSize[.1],Point[atom2[[1]]*a1+atom2[[2]]*a2+atom2[[3]]*a3]},{Red,Arrow[{atom1,a1}],Arrow[{atom1,a2}],Arrow[{atom1,a3}]},{Black,Line[{a1+a2,a2,a2+a3,a3,a1+a3,a1,a1+a2,a1+a2+a3,a1+a3,a1+a2+a3,a2+a3}]},{Black,Text["a1",a1/2],Text["a2",a2/2],Text["a3",a3/2]}},Axes->True]

Berechnung der Matrix

(*Funktion zur Berechnung eines Einheitsvektors zum nächsten Atom.*)
elmn[pos_,atom_]:=Normalize[(pos[[1]]-atom[[1]])*a1+(pos[[2]]-atom[[2]])*a2+(pos[[3]]-atom[[3]])*a3];
(*Funktionen zur Berechnung der Eigenfunktionen des Translationsoperators der u und v Atome*)
ulmn[pos_,k_]:=Exp[I*(pos[[1]]*a1.k+pos[[2]]*a2.k+pos[[3]]*a3.k)];
vlmn[pos_,k_]:=ulmn[pos-atom2,k];

Hier gehts ans Eingemachte, die Erzeugung der Matrixelemente.
Die Eigenfunktionen des Translationsoperators haben die Form:
(Subscript[u, lmn]^x)=(Subscript[u, (OverVector[k])]^x)*(e^(i(l*(Subscript[(OverVector[a]), 1])*(OverVector[k])+m*(Subscript[(OverVector[a]), 2])*(OverVector[k])+n*(Subscript[(OverVector[a]), 3])*(OverVector[k])-ω*t))) 
Die Kraft von einem u Atom in x Richtung hat folgende Form, denn sie ist proportional zur Auslenkung:
(Subscript[F, x])=(∑_lmn)C*(((Subscript[(OverVector[u]), lmn])-(Subscript[(OverVector[u]), 000]))*(Subscript[(OverVector[e]), lmn]))*((OverVector[(Subscript[e, x])])*(Subscript[(OverVector[e]), lmn]))+(∑_lmn)C*(((Subscript[(OverVector[v]), lmn])-(Subscript[(OverVector[u]), 000]))*(Subscript[(OverVector[e]), lmn]))*((OverVector[(Subscript[e, x])])*(Subscript[(OverVector[e]), lmn]))
Wobei hier über die nächsten u und v Atome summiert wird.  Es wird von kleinen Auslenkungen um die Ruhelage ausgegangen, deswegen zeigt der Kraftvektor in Richtung (Subscript[(OverVector[e]), lmn]). Dies wird in ein Gleichungssystem folgender Form gebracht:
-ω2*m*(({{(Subscript[(OverVector[u]), 000])}, {(Subscript[(OverVector[v]), 000])}}))=C*M*(({{(Subscript[(OverVector[u]), 000])}, {(Subscript[(OverVector[v]), 000])}})) 
Diese Matrix M wird hier unten berechnet. Da (Subscript[(OverVector[u]), 000]) und (Subscript[(OverVector[v]), 000]) jeweils 3 Komponenten haben ist es eine 6 * 6 Matrix. Obige Gleichung entsteht einfach durch Einsetzen des Translationsoperators in das Newton'sche Kraftgesetz:
m*((∂^2)/(∂(x^2)))(OverVector[x])=(OverVector[F])
Aus den Eigenwerten der Matrix M folgen dann die Frequenzen:
ω=(Sqrt[-((λ*C)/m)])
Die Matrix wird in Form von vier 3*3 Matrizen zusammengesetzt.

Fuu[k_,i_,j_]:=Simplify[ExpToTrig[Sum[(ulmn[pos1[[n]],k]-ulmn[atom1,k])*elmn[pos1[[n]],atom1][[i]]*elmn[pos1[[n]],atom1][[j]],{n,1,6}]]]-Simplify[ExpToTrig[Sum[(ulmn[atom1,k])*elmn[pos1[[n]],atom1][[i]]*elmn[pos1[[n]],atom1][[j]],{n,7,12}]]];
Fuv[k_,i_,j_]:=Simplify[ExpToTrig[Sum[(vlmn[pos1[[n]],k])*elmn[pos1[[n]],atom1][[i]]*elmn[pos1[[n]],atom1][[j]],{n,7,12}]]];
Fvu[k_,i_,j_]:=Simplify[ExpToTrig[Sum[(ulmn[pos2[[n]],k])*elmn[pos2[[n]],atom2][[i]]*elmn[pos2[[n]],atom2][[j]],{n,1,6}]]];
Fvv[k_,i_,j_]:=Simplify[ExpToTrig[Sum[(vlmn[pos2[[n]],k]-vlmn[atom2,k])*elmn[pos2[[n]],atom2][[i]]*elmn[pos2[[n]],atom2][[j]],{n,7,12}]]]-Simplify[ExpToTrig[Sum[(vlmn[atom2,k])*elmn[pos2[[n]],atom2][[i]]*elmn[pos2[[n]],atom2][[j]],{n,1,6}]]];

M[k_]:=Join[Join[Table[Fuu[k,i,j],{i,1,3},{j,1,3}],Table[Fuv[k,i,j],{i,1,3},{j,1,3}],2],Join[Table[Fvu[k,i,j],{i,1,3},{j,1,3}],Table[Fvv[k,i,j],{i,1,3},{j,1,3}],2]];
M[{kx,ky,kz}]//MatrixForm
Mnum[k1_]=N[M[k1]];

Lösen des Eigenwertproblems

Zunächst wird lediglich der Weg im reziproken Raum erzeugt, auf dem die ω Werte errechnet werden sollen. Hierzu werden einfach Stützstellen zwischen den Symmetriepunkten erzeugt und in k abgespeichert.

leng=Length[symPoints];
Clear[kvec]
kvec={};
Do[kvec=Append[kvec,symPoints[[i+1]]-symPoints[[i]]],{i,1,leng-1}];
Clear[k]
k={{0,0,0}};
Do[Do[k=Append[k,N[symPoints[[j]]+kvec[[j]]*i/100]],{i,1,100}],{j,1,leng-1}];
Clear[kleng]
kleng={0};
Do[kleng=Append[kleng,kleng[[i-1]]+Norm[k[[i-1]]-k[[i]]]//N],{i,2,(leng-1)*100+1}];

Hier wird das eigentliche Eigenwertproblem entlang von k gelöst. Zur graphischen Ausgabe werden die verschiedenen Lösungen in mehrere Listen abgespeichert. Da es eine 6*6 Matrix ist erhält man 6 Eigenwerte.

list=Table[Re[Sqrt[-springconst*Eigenvalues[Mnum[k[[i]]]]/m]],{i,1,(leng-1)*100+1}];
Clear[list1]
list1={};
Do[list1=Append[list1,{kleng[[i]],list[[i]][[1]]}],{i,1,(leng-1)*100+1}]
Clear[list2]
list2={};
Do[list2=Append[list2,{kleng[[i]],list[[i]][[2]]}],{i,1,(leng-1)*100+1}]
Clear[list3]
list3={};
Do[list3=Append[list3,{kleng[[i]],list[[i]][[3]]}],{i,1,(leng-1)*100+1}]
Clear[list4]
list4={};
Do[list4=Append[list4,{kleng[[i]],list[[i]][[4]]}],{i,1,(leng-1)*100+1}]
Clear[list5]
list5={};
Do[list5=Append[list5,{kleng[[i]],list[[i]][[5]]}],{i,1,(leng-1)*100+1}]
Clear[list6]
list6={};
Do[list6=Append[list6,{kleng[[i]],list[[i]][[6]]}],{i,1,(leng-1)*100+1}]

Plotten der Dispersionsrelation

Graphics[{Line[list1],Line[list2],Line[list3],Line[list4],Line[list5],Line[list6],Red,Table[Line[{{kleng[[100*i+1]],0},{kleng[[100*i+1]],3}}],{i,1,leng-1}],Black,Text[label[[1]],{-.1,-.1}],Table[Text[label[[i+1]],{kleng[[100*i+1]]-.1,-.1}],{i,1,leng-1}]},Axes->True,AspectRatio->.5,AxesLabel->{"","ω/(Sqrt[C/m])"}]