Bandstructure of κ-(BEDT-TTF)2Cu(NCS)2
Harukazu YOSHINO, April 22, 2010

<<"VectorAnalysis`"

Off[General::"spell1"]

(* Date and time for record *)
timeStart=DateList[]

Solving the Secular Equation

(* Lefthand side of the secular equation *)
	H={{h11, h12, h13, h14},
		{h21, h22, h23, h24},
		{h31, h32, h33, h34},
 	 {h41, h42, h43, h44}}

(* Diagonal elements of the 4 x 4 matrix *)
h11=-ene;h22=h11;h33=h11;h44=h11

(* Determinant of the matrix *)
J0=Simplify[Det[H]]

(* Offdiagonal elements of the matrix  Refs.1,2 *)
h12=Exp[I k.ab]Bab + Exp[I k.ad] Bad;
h21=h12 /. I->-I;
h13=Exp[I k.ac] Bac + Exp[I k.af] Baf;
h31=h13/. I->-I;
h14=Exp[I k.ae] Bae + Exp[I k.ag]Bag;
h41=h14/.I->-I;
h23=Exp[I k.bc]Bbc + Exp[I k.bi]Bbi;
h32=h23/.I->-I;
h24=Exp[I k.bg]Bbg + Exp[I k.bj] Bbj;
h42=h24/.I->-I;
h34=Exp[I k.cl]Bcl + Exp[I k.cj]Bcj;
h43=h34/.I->-I;

(* Positional vectors of each BEDT-TTF molecule   Refs.1.2 *)
(* Vector ab is defined to be (2r, 2s) *)
oa={(1/2)b+r, (1/2)c+s};
ob={(1/2)b-r, (1/2)c-s};
oc={b-r, s};
od={(3/2)b-r, (1/2)c-s};
oe={b+r, c-s};
of={b-r, c+s};
og={r, c-s};
oi={-r, s};
oj={r, -s};
ol={b+r, -s};

(* Relative vectors between neighbouring molecules *)
ab=ob-oa;
ad=od-oa;
ac=oc-oa;
af=of-oa;
ae=oe-oa;
ag=og-oa;
bc=oc-ob;
bi=oi-ob;
bg=og-ob;
bj=oj-ob;
cl=ol-oc;
cj=oj-oc;

(* Transfer integrals between neighbouring molecules  Ref.3 *)
Bab=tb1;Bcl=tb1;
Bae=tp1;Bag=tp1;
Bbc=tp2;Bbi=tp2;
Baf=tq1;Bbg=tq1;
Bac=tq2;Bbj=tq2;
Bad=tb2;Bcj=tb2;

(* k-vector *)
k={kb,kc}

(* Expand and simplify the determinant, J *)
J0=ComplexExpand[J0];

J0=FullSimplify[Collect[TrigExpand[Simplify[J0]],ene]]

(* Replace characters of transferintegrals with number in meV *)
J=J0/.{tb1->-257., tb2->-105., tp1->-114., tp2->-100., tq1->17., tq2->29.}

(* Replace b kb and c kc with theta and phi for simplicity *)
J=J/.{b kb->theta, c kc->phi}

(* Solve the equation *)
solution=Solve[J==0, ene];

(* Define the 4 bands using the above solution *)
kappa1=ene/.solution[[1]][[1]];
kappa2=ene/.solution[[2]][[1]];
kappa3=ene/.solution[[3]][[1]];
kappa4=ene/.solution[[4]][[1]];

Chop[kappa1/.{theta->-3.14159, phi->-1.15743},10^-5]

(* Band 1 *)
ParametricPlot3D[{theta, phi,kappa1}, {theta, -Pi, Pi}, {phi, -Pi, Pi}, BoxRatios->{1, 8.440/13.124, 1},
PlotRange->{{-Pi,Pi},{-Pi,Pi},{-700,700}},PlotPoints->6,
 AxesLabel->{kb,kc,E/eV}]

(* Band 2 *)
ParametricPlot3D[{theta, phi,kappa2}, {theta, -Pi, Pi}, {phi, -Pi, Pi}, BoxRatios->{1, 8.440/13.124, 1},
PlotRange->{{-Pi,Pi},{-Pi,Pi},{-700,700}},PlotPoints->6,
 AxesLabel->{kb,kc,E/eV}]

(* Band 3 *)
ParametricPlot3D[{theta, phi,kappa3}, {theta, -Pi, Pi}, {phi, -Pi, Pi}, BoxRatios->{1, 8.440/13.124, 1},
PlotRange->{{-Pi,Pi},{-Pi,Pi},{-700,700}},PlotPoints->6,
 AxesLabel->{kb,kc,E/eV}]

(* Band 4 *)
ParametricPlot3D[{theta, phi,kappa4}, {theta, -Pi, Pi}, {phi, -Pi, Pi}, BoxRatios->{1, 8.440/13.124, 1},
PlotRange->{{-Pi,Pi},{-Pi,Pi},{-700,700}},PlotPoints->6,
 AxesLabel->{kb,kc,E/eV}]

ParametricPlot3D[{{theta, phi,kappa1},{theta, phi,kappa2},{theta, phi,kappa3},{theta, phi,kappa4}}, {theta, -Pi, Pi}, {phi, -Pi, Pi}, BoxRatios->{1, 8.440/13.124, 1},PlotRange->{{-Pi,Pi},{-Pi,Pi},{-700,700}},
 AxesLabel->{kb,kc,E/eV},PlotPoints->6]

Calculate Fermi Energy

band={kappa1,kappa2};

Counting version

(* Division *)
nk=7;  (* increase this to get a finer value of Ef *)
div=2^nk;
(* Cutoff *)
dSkCut = 10^-5;
(* Area to be calculated and 1st Burillouin zone *)
Sk=0;
(* Change this if you want another filling, for example, 1/4-filled -> 0.5, 1/2-filled->1, 3/4-filled-> 1.5 *)
SkF=1;(* 1/2-filled of 2 bands = 3/4-filled of 4 bands *)

e1=Chop[band[[1]]/.{theta->0,phi->0},10^-5]; (* Chop is needed to truncate small imaginary numbers. *)
e2=Chop[band[[1]]/.{theta->Pi,phi->Pi},10^-5];
e3=Chop[band[[2]]/.{theta->0,phi->0},10^-5];
e4=Chop[band[[2]]/.{theta->Pi,phi->Pi},10^-5];
eneMin=Min[e1,e2,e3,e4];
eneMax=Max[e1,e2,e3,e4];
Ef1={eneMin,eneMax,(eneMin+eneMax)/2}; (* Just to define the list *)
Sk1={0,0,0};(* Just to define the list *)

(* Make table of energy *)
eneB1=Table[Chop[band[[1]]/.{theta->Pi/div*(i+.5),phi->Pi/div*(j+.5)},10^-5],{i,0,div-1},{j,0,div-1}];
eneB2=Table[Chop[band[[2]]/.{theta->Pi/div*(i+.5),phi->Pi/div*(j+.5)},10^-5],{i,0,div-1},{j,0,2*div-1}];

(* Find Ef of half filling for total of bands 1 and 2 *)
dSk=1;
count2=0;
While[dSk>dSkCut && count2<30,
count2++;
(* Calculate three kinds of Sk *)
Ef1[[3]]=(Ef1[[1]]+Ef1[[2]])/2;
For[j=1,j≤3,j++,
Sk1[[j]]=0;
Ef=Ef1[[j]];
Sk=0;
For[p=1,p<=div,p++,
For[q=1,q<=div,q++,
(* Sk *)
If[eneB1[[p]][[q]]≤Ef,Sk++];
If[eneB2[[p]][[q]]≤Ef,Sk++];
];(* For(phi1) *)
];(* For(theta1) *)
Sk=Sk/(div^2);
Sk1[[j]]=Sk-SkF; (* difference between calculated S and SF *)
];(* For(j) *)
Print["EF ",N[Ef1[[1]],8],"  ",N[Ef1[[3]],8],"  ",N[Ef1[[2]],8],"  ",
"Sk ",N[Sk1[[1]]+SkF,8],"  ",N[Sk1[[3]]+SkF,8],"  ",N[Sk1[[2]]+SkF,8]];
If[Sk1[[1]]*Sk1[[3]]<0,Ef1[[2]]=Ef1[[3]],Ef1[[1]]=Ef1[[3]]];
dSk=Abs[Sk1[[3]]];
];
(* While *)
Print["Ef=",Ef1[[3]]," meV"]

ParametricPlot3D[{{theta, phi,kappa1},{theta, phi,kappa2}}, {theta, -Pi, Pi}, {phi, -Pi, Pi}, BoxRatios->{1, 8.440/13.124, 1},
 AxesLabel->{kb,kc,E/eV},PlotPoints->6,
PlotRange->{{-Pi,Pi},{-Pi,Pi},{0,Ef}}]

(* Date and time for record *)
timeEnd=DateList[];
timeEnd-timeStart

References
[1]	H. Urayama, H. Yamochi, G. Saito, S. Sato,  A. Kawamoto, J. Tanaka, T. Mori, Y. Maruyama and H. Inokuchi,
	Chem. Lett. (1988) 463.
[2]	D. Jung, M. Evain, J. J. Nova, M.-H. Whangbo,  M. A. Beno, A. M. Kini, A. J. Schultz, J. M. Williams and P. J. Nigrey,
	 Inorg. Chem. 28 (1989) 4516.
[3]	K. Oshima, T. Mori, H. Inokuchi, H. Urayama, H. Yamochi and G. Saito, Phys. Rev. B38 (1988) 938.