|
点電荷モデルと結晶場パラメータ
結晶場を実際の点電荷による静電ポテンシャルで考える．周囲の点電荷の配置が一見複雑で，どんな結晶場パラメータが現れるのかすぐには判断できないが，対称性はきちんと考えて計算したいときに有用．実際の対称性のとおりに点電荷を置き，点電荷の位置と有効電荷を入力すると，結晶場の展開係数q_lmを計算する，これが，球テンソル型の結晶場演算子に対するパラメータとなる．あらゆる対称性に適用可能．最後に，出てきたq_lmを立方調和型に変換し，結晶場を多極子表記すると，対称性が低い場合でも見通しがよくなる．また，従来型結晶場演算子O_lmの係数B_lmに変換したり，結晶場を多極子場表記で表示したりするルーチンも用意した．ただし，従来型結晶場演算子はl=even, m=evenの項のみに対応． Point charges and the CEF parameters
-
(* point charges *)
pcsites={{2,0,0}, {0,2,0}, {0,0,2}, {-2,0,0}, {0,-2,0}, {0,0,-2}}*0.1;
charges={-1,-1,-1,-1,-1,-1}; (* Radial averages of <r^2>, <r^4>, <r^6> *)
rave[2]=0.2188; rave[4]=0.118; rave[6]=0.1328; (* Expansion coefficients q_lm for spherical tensor CEF operators *)
pctp=Module[{xd,yd,zd,r,phi,theta},
Table[{xd=pcsites[[i,1]];yd=pcsites[[i,2]];zd=pcsites[[i,3]];
r=Sqrt[xd^2+yd^2+zd^2];
Which[xd==0 && yd>0,phi=Pi/2,xd==0 && yd<0,phi=3*Pi/2,yd==0 && xd>0,phi=0,
yd==0 && xd<0,phi=Pi,xd==0 && yd==0,phi=0,True,phi=ArcTan[xd,yd]];
Which[zd>0 && xd==0 && yd==0,theta=0,zd<0 && xd==0 && yd==0,theta=Pi,True,
theta=ArcCos[zd/r]]; r,theta,phi},{i,1,Length[pcsites]}]];
Print["Ion Positions (r,theta,phi)="]; Print[pctp];
Do[ Do[ qfac[k,m] = Chop[N[
Sqrt[4*Pi/(2*k+1)]*Sum[ -charges[[i]]*(-1)^m
*SphericalHarmonicY[k,-m,pctp[[i,2]],pctp[[i,3]]]/pctp[[i,1]]^(k+1)
,{i,1,Length[pctp]}] ] ], {m,k,-k,-1}], {k,0,6}];
Print["Expansion coefficient of the CEF potential. q(odd,m)=zero when there is inversion symmetry."]
Print["q(0,m)=",Table[qfac[0,m],{m,0,0}]];
Print["q(1,m)=",Table[qfac[1,m],{m,1,-1,-1}]];
Print["q(2,m)=",Table[qfac[2,m],{m,2,-2,-1}]];
Print["q(3,m)=",Table[qfac[3,m],{m,3,-3,-1}]];
Print["q(4,m)=",Table[qfac[4,m],{m,4,-4,-1}]];
Print["q(5,m)=",Table[qfac[5,m],{m,5,-5,-1}]];
Print["q(6,m)=",Table[qfac[6,m],{m,6,-6,-1}]]; (* factors to connect conventional and spherical tensor CEF operators *)
facB20=2; facB22=2/Sqrt[6];
facB40=8; facB42=8/Sqrt[5*8]; facB44=8/Sqrt[2*5*7];
facB60=16; facB62=16/Sqrt[3*5*7]; facB64=16/Sqrt[2*7*9]; facB66=16/Sqrt[3*7*11];
-
(* B_lm parameters for the conventional CEF operators: even terms only *) b20=Chop[qfac[2,0]*rave[2]*stvfac[n4f,2]/facB20/redfac3[2]];
b22=Chop[qfac[2,2]*rave[2]*stvfac[n4f,2]/facB22/redfac3[2]];
b40=Chop[qfac[4,0]*rave[4]*stvfac[n4f,4]/facB40/redfac3[4]];
b42=Chop[qfac[4,2]*rave[4]*stvfac[n4f,4]/facB42/redfac3[4]];
b44=Chop[qfac[4,4]*rave[4]*stvfac[n4f,4]/facB44/redfac3[4]];
b60=Chop[qfac[6,0]*rave[6]*stvfac[n4f,6]/facB60/redfac3[6]];
b62=Chop[qfac[6,2]*rave[6]*stvfac[n4f,6]/facB62/redfac3[6]];
b64=Chop[qfac[6,4]*rave[6]*stvfac[n4f,6]/facB64/redfac3[6]];
b66=Chop[qfac[6,6]*rave[6]*stvfac[n4f,6]/facB66/redfac3[6]]; Print["CEF coefficients B(l,m) (K) (for cubic, tetragonal, orthorhombic)"]
Print["B(2,0)=",b20,Tab,"B(2,2)=",b22]
Print["B(4,0)=",b40,Tab,"B(4,2)=",b42,Tab,"B(4,4)=",b44]
Print["B(6,0)=",b60,Tab,"B(6,2)=",b62,Tab,"B(6,4)=",b64,Tab,"B(6,6)=",b66]
-
(* transformation from B_lm to qfac[l, m]: even terms only *) q20 = b20/rave[2]/stvfac[n4f, 2]*facB20*redfac3[2];
q22 = b22/rave[2]/stvfac[n4f, 2]*facB22*redfac3[2];
q40 = b40/rave[4]/stvfac[n4f, 4]*facB40*redfac3[4];
q42 = b42/rave[4]/stvfac[n4f, 4]*facB42*redfac3[4];
q44 = b44/rave[4]/stvfac[n4f, 4]*facB44*redfac3[4];
q60 = b60/rave[6]/stvfac[n4f, 6]*facB60*redfac3[6];
q62 = b62/rave[6]/stvfac[n4f, 6]*facB62*redfac3[6];
q64 = b64/rave[6]/stvfac[n4f, 6]*facB64*redfac3[6];
q66 = b66/rave[6]/stvfac[n4f, 6]*facB66*redfac3[6]; {qfac[2,2],qfac[2,1],qfac[2,0],qfac[2,-1],qfac[2,-2]}={q22,0,q20,0,q22};
{qfac[4,4],qfac[4,3],qfac[4,2],qfac[4,1],qfac[4,0],
qfac[4,-1],qfac[4,-2],qfac[4,-3],qfac[4,-4]} = {q44,0,q42,0,q40,0,q42,0,q44};
{qfac[6,6],qfac[6,5],qfac[6,4],qfac[6,3],qfac[6,2],
qfac[6,1],qfac[6,0],qfac[6,-1],qfac[6,-2],qfac[6,-3],
qfac[6,-4],qfac[6,-5],qfac[6,-6]}={q66,0,q64,0,q62,0,q60,0,q62,0,q64,0,q66};
Print["Expansion coefficient of the CEF potential. Even terms only."]
Print["q(2,m)=",Table[qfac[2,m],{m,2,-2,-1}]];
Print["q(4,m)=",Table[qfac[4,m],{m,4,-4,-1}]];
Print["q(6,m)=",Table[qfac[6,m],{m,6,-6,-1}]];
(* Multipolar expansion of the CEF using cubic harmonic representation *)
-
(* rank 0 : dipole field *)
a1g[0, 1] = qfac[0, 0]; (* rank 1 *)
ceft1u[1,1]=(qfac[1,-1]-qfac[1,1])/Sqrt[2];
ceft1u[1,2]=-(I*(qfac[1,-1]+qfac[1,1]))/Sqrt[2];
ceft1u[1,3]=qfac[1,0]; (* rank 2 : quadrupole field *)
cefeg[2,1]=qfac[2,0];
cefeg[2,2]=(qfac[2,-2]+qfac[2,2])/Sqrt[2];
ceft2g[2,1]=-(I*(qfac[2,-1]+qfac[2,1]))/Sqrt[2];
ceft2g[2,2]=(qfac[2,-1]-qfac[2,1])/Sqrt[2];
ceft2g[2,3]=-(I*(qfac[2,-2]-qfac[2,2]))/Sqrt[2]; (* rank 3 : octupole field *)
cefa2u[3,1]=-(I*(qfac[3,-2]-qfac[3,2]))/Sqrt[2];
ceft1u[3,1]=(Sqrt[5]*qfac[3,-3]-Sqrt[3]*qfac[3,-1]+Sqrt[3]*qfac[3,1]-Sqrt[5]*qfac[3,3])/4;
ceft1u[3,2]=I/4*(Sqrt[5]*qfac[3,-3]+Sqrt[3]*qfac[3,-1]+Sqrt[3]*qfac[3,1]+Sqrt[5]*qfac[3,3]);
ceft1u[3,3]=qfac[3,0];
ceft2u[3,1]=(-(Sqrt[3]*qfac[3,-3])-Sqrt[5]*qfac[3,-1]+Sqrt[5]*qfac[3,1]+Sqrt[3]*qfac[3,3])/4;
ceft2u[3,2]=I/4*(Sqrt[3]*qfac[3,-3]-Sqrt[5]*qfac[3,-1]-Sqrt[5]*qfac[3,1]+Sqrt[3]*qfac[3,3]);
ceft2u[3,3]=(qfac[3,-2]+qfac[3,2])/Sqrt[2]; (* rank 4 : hexadecapole field *)
cefa1g[4,1]=(Sqrt[30]*qfac[4,-4]+2*Sqrt[21]*qfac[4,0]+Sqrt[30]*qfac[4,4])/12;
cefeg[4,1]=(-(Sqrt[42]*qfac[4,-4])+2*Sqrt[15]*qfac[4,0]-Sqrt[42]*qfac[4,4])/12;
cefeg[4,2]=-((qfac[4,-2]+qfac[4,2])/Sqrt[2]);
ceft1g[4,1]=I/4*(qfac[4,-3]+Sqrt[7]*qfac[4,-1]+Sqrt[7]*qfac[4,1]+qfac[4,3]);
ceft1g[4,2]=(-qfac[4,-3]+Sqrt[7]*qfac[4,-1]-Sqrt[7]*qfac[4,1]+qfac[4,3])/4;
ceft1g[4,3]=-(I*(qfac[4,-4]-qfac[4,4]))/Sqrt[2];
ceft2g[4,1]=-I/4*(Sqrt[7]*qfac[4,-3]-qfac[4,-1]-qfac[4,1]+Sqrt[7]*qfac[4,3]);
ceft2g[4,2]=(-(Sqrt[7]*qfac[4,-3])-qfac[4,-1]+qfac[4,1]+Sqrt[7]*qfac[4,3])/4;
ceft2g[4,3]=-(I*(qfac[4,-2]-qfac[4,2]))/Sqrt[2];
dipCEF=Chop[{ceft1u[1,1],ceft1u[1,2],ceft1u[1,3]}];
quadCEF=Chop[{cefeg[2,2],cefeg[2,1],ceft2g[2,1],ceft2g[2,2],ceft2g[2,3]}];
octCEF=Chop[{cefa2u[3,1],ceft1u[3,1],ceft1u[3,2],ceft1u[3,3],
ceft2u[3,1],ceft2u[3,2],ceft2u[3,3]}];
hexCEF=Chop[{cefa1g[4,1],cefeg[4,1],cefeg[4,2],ceft1g[4,1],ceft1g[4,2],ceft1g[4,3],
ceft2g[4,1],ceft2g[4,2],ceft2g[4,3]}];
Print["dip:", dipCEF]; Print["quad:", quadCEF];
Print["oct:", octCEF]; Print["hex:", hexCEF];
|
